Risoluzione di un'equazione lineare
https://stackoverflow.com/questions/769
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianDomanda
Devo risolvere a livello di codice un sistema di equazioni lineari in C, Objective C o (se necessario) C++.
Ecco un esempio delle equazioni:
-44.3940 = a * 50.0 + b * 37.0 + tx
-45.3049 = a * 43.0 + b * 39.0 + tx
-44.9594 = a * 52.0 + b * 41.0 + tx
Da questo, mi piacerebbe ottenere la migliore approssimazione per a, b, E tx.
Soluzione
Regola di CramerEEliminazione gaussianasono due buoni algoritmi di uso generale (vedi anche Equazioni lineari simultanee).Se stai cercando il codice, dai un'occhiata GiNaC, Massimi, E C++ simbolico (a seconda dei requisiti di licenza, ovviamente).
MODIFICARE:So che lavori in C, ma devo anche mettere una buona parola per questo SymPy (un sistema di algebra informatica in Python).Puoi imparare molto dai suoi algoritmi (se sai leggere un po' di Python).Inoltre, è sotto la nuova licenza BSD, mentre la maggior parte dei pacchetti matematici gratuiti sono GPL.
Altri suggerimenti
Puoi risolverlo con un programma esattamente nello stesso modo in cui lo risolvi manualmente (con moltiplicazione e sottrazione, quindi reimmettendo i risultati nelle equazioni).Questa è matematica piuttosto standard a livello di scuola secondaria.
-44.3940 = 50a + 37b + c (A)
-45.3049 = 43a + 39b + c (B)
-44.9594 = 52a + 41b + c (C)
(A-B): 0.9109 = 7a - 2b (D)
(B-C): 0.3455 = -9a - 2b (E)
(D-E): 1.2564 = 16a (F)
(F/16): a = 0.078525 (G)
Feed G into D:
0.9109 = 7a - 2b
=> 0.9109 = 0.549675 - 2b (substitute a)
=> 0.361225 = -2b (subtract 0.549675 from both sides)
=> -0.1806125 = b (divide both sides by -2) (H)
Feed H/G into A:
-44.3940 = 50a + 37b + c
=> -44.3940 = 3.92625 - 6.6826625 + c (substitute a/b)
=> -41.6375875 = c (subtract 3.92625 - 6.6826625 from both sides)
Quindi ti ritroverai con:
a = 0.0785250
b = -0.1806125
c = -41.6375875
Se inserisci nuovamente questi valori in A, B e C, scoprirai che sono corretti.
Il trucco sta nell'usare una semplice matrice 4x3 che si riduce a sua volta in una matrice 3x2, quindi una matrice 2x1 che è "a = n", dove n è un numero reale.Una volta ottenuto questo, inseriscilo nella matrice successiva per ottenere un altro valore, quindi questi due valori nella matrice successiva finché non hai risolto tutte le variabili.
A condizione che tu abbia N equazioni distinte, puoi sempre risolvere per N variabili.Dico distinti perché questi due non lo sono:
7a + 2b = 50
14a + 4b = 100
Loro sono il Stesso equazione moltiplicata per due in modo da non poterne ottenere una soluzione - moltiplicando la prima per due e poi sottraendo si ottiene l'affermazione vera ma inutile:
0 = 0 + 0
A titolo di esempio, ecco del codice C che risolve le equazioni simultanee che ti vengono inserite nella tua domanda.Innanzitutto alcuni tipi necessari, variabili, una funzione di supporto per stampare un'equazione e l'inizio di main:
#include <stdio.h>
typedef struct { double r, a, b, c; } tEquation;
tEquation equ1[] = {
{ -44.3940, 50, 37, 1 }, // -44.3940 = 50a + 37b + c (A)
{ -45.3049, 43, 39, 1 }, // -45.3049 = 43a + 39b + c (B)
{ -44.9594, 52, 41, 1 }, // -44.9594 = 52a + 41b + c (C)
};
tEquation equ2[2], equ3[1];
static void dumpEqu (char *desc, tEquation *e, char *post) {
printf ("%10s: %12.8lf = %12.8lfa + %12.8lfb + %12.8lfc (%s)\n",
desc, e->r, e->a, e->b, e->c, post);
}
int main (void) {
double a, b, c;
Successivamente, la riduzione delle tre equazioni in tre incognite a due equazioni in due incognite:
// First step, populate equ2 based on removing c from equ.
dumpEqu (">", &(equ1[0]), "A");
dumpEqu (">", &(equ1[1]), "B");
dumpEqu (">", &(equ1[2]), "C");
puts ("");
// A - B
equ2[0].r = equ1[0].r * equ1[1].c - equ1[1].r * equ1[0].c;
equ2[0].a = equ1[0].a * equ1[1].c - equ1[1].a * equ1[0].c;
equ2[0].b = equ1[0].b * equ1[1].c - equ1[1].b * equ1[0].c;
equ2[0].c = 0;
// B - C
equ2[1].r = equ1[1].r * equ1[2].c - equ1[2].r * equ1[1].c;
equ2[1].a = equ1[1].a * equ1[2].c - equ1[2].a * equ1[1].c;
equ2[1].b = equ1[1].b * equ1[2].c - equ1[2].b * equ1[1].c;
equ2[1].c = 0;
dumpEqu ("A-B", &(equ2[0]), "D");
dumpEqu ("B-C", &(equ2[1]), "E");
puts ("");
Successivamente, la riduzione delle due equazioni con due incognite ad un'equazione con un'incognita:
// Next step, populate equ3 based on removing b from equ2.
// D - E
equ3[0].r = equ2[0].r * equ2[1].b - equ2[1].r * equ2[0].b;
equ3[0].a = equ2[0].a * equ2[1].b - equ2[1].a * equ2[0].b;
equ3[0].b = 0;
equ3[0].c = 0;
dumpEqu ("D-E", &(equ3[0]), "F");
puts ("");
Ora che abbiamo una formula del tipo number1 = unknown * number2, possiamo semplicemente calcolare il valore sconosciuto con unknown <- number1 / number2.Poi, una volta calcolato il valore, sostituiscilo in una delle equazioni a due incognite e calcola il secondo valore.Quindi sostituisci entrambe le incognite (ora note) in una delle equazioni originali e ora avrai i valori per tutte e tre le incognite:
// Finally, substitute values back into equations.
a = equ3[0].r / equ3[0].a;
printf ("From (F ), a = %12.8lf (G)\n", a);
b = (equ2[0].r - equ2[0].a * a) / equ2[0].b;
printf ("From (D,G ), b = %12.8lf (H)\n", b);
c = (equ1[0].r - equ1[0].a * a - equ1[0].b * b) / equ1[0].c;
printf ("From (A,G,H), c = %12.8lf (I)\n", c);
return 0;
}
L'output di quel codice corrisponde ai calcoli precedenti in questa risposta:
>: -44.39400000 = 50.00000000a + 37.00000000b + 1.00000000c (A)
>: -45.30490000 = 43.00000000a + 39.00000000b + 1.00000000c (B)
>: -44.95940000 = 52.00000000a + 41.00000000b + 1.00000000c (C)
A-B: 0.91090000 = 7.00000000a + -2.00000000b + 0.00000000c (D)
B-C: -0.34550000 = -9.00000000a + -2.00000000b + 0.00000000c (E)
D-E: -2.51280000 = -32.00000000a + 0.00000000b + 0.00000000c (F)
From (F ), a = 0.07852500 (G)
From (D,G ), b = -0.18061250 (H)
From (A,G,H), c = -41.63758750 (I)
Per un sistema 3x3 di equazioni lineari immagino che andrebbe bene implementare i propri algoritmi.
Tuttavia, potresti doverti preoccupare della precisione, della divisione per zero o dei numeri molto piccoli e di cosa fare con infinite soluzioni.Il mio suggerimento è di utilizzare un pacchetto di algebra lineare numerica standard come LAPACK.
Dai un'occhiata a Fondazione risolutore Microsoft.
Con esso potresti scrivere codice come questo:
SolverContext context = SolverContext.GetContext();
Model model = context.CreateModel();
Decision a = new Decision(Domain.Real, "a");
Decision b = new Decision(Domain.Real, "b");
Decision c = new Decision(Domain.Real, "c");
model.AddDecisions(a,b,c);
model.AddConstraint("eqA", -44.3940 == 50*a + 37*b + c);
model.AddConstraint("eqB", -45.3049 == 43*a + 39*b + c);
model.AddConstraint("eqC", -44.9594 == 52*a + 41*b + c);
Solution solution = context.Solve();
string results = solution.GetReport().ToString();
Console.WriteLine(results);
Ecco l'output:
===Rapporto sul servizio Solver Foundation===
Appuntamento:20/04/2009 23:29:55
Nome del modello:Predefinito
Capacità richieste:LP
Tempo di risoluzione (ms):1027
Tempo totale (ms):1414
Stato di completamento risoluzione:Ottimale
Risolutore selezionato:Microsoft.SolverFoundation.Solvers.SimplexSolver
Direttive:
Microsoft.SolverFoundation.Services.Directive
Algoritmo:Primordiale
Aritmetica:Ibrido
Prezzo (esatto):Predefinito
Prezzo (doppio):Bordo più ripido
Base:Lento
Conteggio pivot:3
===Dettagli della soluzione===
Obiettivi:
Decisioni:
UN:0,0785250000000004
B:-0,180612500000001
C:-41.6375875
Stai cercando un pacchetto software che svolga il lavoro o esegua effettivamente le operazioni sulla matrice e simili ed esegua ogni passaggio?
Il primo, l'ha appena usato un mio collega OcamlGLPK.È solo un involucro per il GLPK, ma rimuove molti passaggi di configurazione.Sembra che dovrai restare con il GLPK, in C, però.Per quest'ultimo, grazie a Delicious per aver salvato un vecchio articolo con cui ho imparato l'LP qualche tempo fa, PDF.Se hai bisogno di aiuto specifico per impostare ulteriormente, faccelo sapere e sono sicuro che io o qualcuno torneremo ad aiutarti, ma penso che sia abbastanza semplice da qui.Buona fortuna!
Kit di strumenti numerici modello del NIST dispone degli strumenti per farlo.
Uno dei modi più affidabili è utilizzare a Decomposizione QR.
Ecco un esempio di wrapper in modo che io possa chiamare "GetInverse(A, InvA)" nel mio codice e inserirà l'inverso in InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D è definito nella libreria.
Dalla formulazione della tua domanda, sembra che tu abbia più equazioni che incognite e desideri ridurre al minimo le incongruenze.Questo viene tipicamente fatto con la regressione lineare, che minimizza la somma dei quadrati delle incoerenze.A seconda della dimensione dei dati, puoi farlo in un foglio di calcolo o in un pacchetto statistico.R è un pacchetto gratuito di alta qualità che esegue la regressione lineare, tra molte altre cose.C'è molto nella regressione lineare (e molti trucchi), ma poiché è semplice da fare per casi semplici.Ecco un esempio R che utilizza i tuoi dati.Nota che "tx" è l'intercettazione del tuo modello.
> y <- c(-44.394, -45.3049, -44.9594)
> a <- c(50.0, 43.0, 52.0)
> b <- c(37.0, 39.0, 41.0)
> regression = lm(y ~ a + b)
> regression
Call:
lm(formula = y ~ a + b)
Coefficients:
(Intercept) a b
-41.63759 0.07852 -0.18061
In termini di efficienza di runtime, altri hanno risposto meglio di me.Se avrai sempre lo stesso numero di equazioni come variabili, mi piace Regola di Cramer poiché è facile da implementare.Basta scrivere una funzione per calcolare il determinante di una matrice (o usarne una già scritta, sono sicuro che puoi trovarne una là fuori) e dividere i determinanti di due matrici.
Personalmente, ho un debole per gli algoritmi di Ricette numeriche.(Sono affezionato all'edizione C++.)
Questo libro ti insegnerà perché gli algoritmi funzionano, oltre a mostrarti alcune implementazioni di tali algoritmi con un buon debug.
Certo, potresti semplicemente usarlo alla cieca CLAPACK (L'ho usato con grande successo), ma vorrei prima digitare un algoritmo di eliminazione gaussiana per avere almeno una vaga idea del tipo di lavoro svolto per rendere stabili questi algoritmi.
Più tardi, se stai facendo un'algebra lineare più interessante, guardando il codice sorgente di Ottava risponderà a molte domande.
function x = LinSolve(A,y)
%
% Recursive Solution of Linear System Ax=y
% matlab equivalent: x = A\y
% x = n x 1
% A = n x n
% y = n x 1
% Uses stack space extensively. Not efficient.
% C allows recursion, so convert it into C.
% ----------------------------------------------
n=length(y);
x=zeros(n,1);
if(n>1)
x(1:n-1,1) = LinSolve( A(1:n-1,1:n-1) - (A(1:n-1,n)*A(n,1:n-1))./A(n,n) , ...
y(1:n-1,1) - A(1:n-1,n).*(y(n,1)/A(n,n)));
x(n,1) = (y(n,1) - A(n,1:n-1)*x(1:n-1,1))./A(n,n);
else
x = y(1,1) / A(1,1);
end
