Lösen einer linearen Gleichung

Question

Ich muss ein System linearer Gleichungen in C, Objective C oder (falls erforderlich) C++ programmatisch lösen.

Hier ist ein Beispiel der Gleichungen:

-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 

Daraus möchte ich die beste Annäherung bekommen für a, b und tx.

Answer 1

Cramer's Rule und Gaussian Elimination sind zwei gute, universelle Algorithmen (auch Simultaneous Linear Equations sehen). Wenn Sie nach Code suchen, sehen Sie sich GiNaC, Maxima und SymbolicC++ an (natürlich abhängig von Ihren Lizenzanforderungen).

EDIT: Ich weiß, dass Sie in C Land arbeiten, aber ich muss auch ein gutes Wort für SymPy (ein Computer-Algebra-System in Python). Sie können eine Menge von seinen Algorithmen lernen (wenn Sie ein wenig Python lesen können). Es ist auch unter der neuen BSD-Lizenz, während die meisten kostenlosen Mathe-Pakete GPL sind.

Answer 2

Sind Sie auf der Suche nach einem Softwarepaket, das die Arbeit oder tatsächlich die Matrix-Operationen und dergleichen tun und jeden Schritt tun?

Das erste, ein Mitarbeiter von mir gerade verwendet Ocaml GLPK. Es ist nur ein Wrapper für die GLPK, aber es entfernt eine Menge der Schritte der Einrichtung. Es sieht so aus, als ob du in C bei der GLPK bleiben musst. Für letzteres, dank lecker zum Speichern eines alten Artikels, lernte ich LP eine Weile zurück, . Wenn Sie spezielle Hilfe benötigen, lassen Sie es uns wissen, und ich bin mir sicher, dass ich oder jemand wieder reingehen und helfen wird, aber ich denke, dass es von hier aus ziemlich einfach ist. Viel Glück!

Answer 3

Für ein 3x3-System linearer Gleichungen wäre es wohl in Ordnung, eigene Algorithmen einzuführen.

Allerdings müssen Sie sich Sorgen um Genauigkeit, Division durch Null oder wirklich kleine Zahlen und was Sie tun, um unendlich viele Lösungen zu machen. Mein Vorschlag ist, mit einem numerischen linearen Standardalgebra-Paket wie LAPACK zu gehen.

Answer 4

Persönlich bin ich teilweise zu den Algorithmen von Numerical Recipes. (Ich mag die C++ Edition.)

Dieses Buch wird Ihnen zeigen, warum die Algorithmen funktionieren, und Ihnen einige ziemlich gut getestete Implementierungen dieser Algorithmen zeigen.

Natürlich könnten Sie blind CLAPACK verwenden (ich habe es mit großem Erfolg verwendet), aber ich würde zuerst einen Gaussian Elimination-Algorithmus eingeben, um zumindest eine schwache Vorstellung von der Art der Arbeit zu haben, die gegangen ist um diese Algorithmen stabil zu machen.

Später, wenn Sie mehr interessante lineare Algebra tun, wird die Suche nach dem Quellcode von Octave eine Menge Fragen beantworten.

Answer 5

Template Numerical Toolkit von NIST hat Tools dafür.

Eine der zuverlässigeren Methoden ist die Verwendung eines QR Decomposition.

Hier ist ein Beispiel für einen Wrapper, so dass ich "GetInverse (A, InvA)" in meinem Code aufrufen kann, und es wird die Umkehrung in InvA setzen.

void GetInverse(const Array2D<double>& A, Array2D<double>& invA) 
    { 
    QR<double> qr(A); 
    invA = qr.solve(I); 
    } 

Array2D ist in der Bibliothek definiert.

Answer 6

Aus dem Wortlaut Ihrer Frage scheint es, als hätten Sie mehr Gleichungen als Unbekannte und Sie möchten die Inkonsistenzen minimieren. Dies geschieht typischerweise mit linearer Regression, die die Summe der Quadrate der Inkonsistenzen minimiert. Abhängig von der Größe der Daten können Sie dies in einer Tabelle oder in einem statistischen Paket tun. R ist ein hochwertiges, kostenloses Paket, das neben vielen anderen Dingen lineare Regressionen durchführt. Es gibt eine Menge an linearer Regression (und eine Menge Fehler), aber für einfache Fälle ist es einfach. Hier ist ein R-Beispiel, das Ihre Daten verwendet. Beachten Sie, dass der "tx" der Abschnitt zu Ihrem Modell ist.

> 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 

Answer 7

In Bezug auf die Laufzeiteffizienz, andere haben beantwortet besser als ich Wenn Sie immer die gleiche Anzahl von Gleichungen als Variablen haben, Ich mag Cramer's rule wie es einfach ist, zu implementieren. Schreiben Sie einfach eine Funktion, um die Determinante einer Matrix zu berechnen (oder verwenden Sie eine, die bereits geschrieben wurde, ich bin mir sicher, dass Sie eine finden können), und teilen Sie die Determinanten von zwei Matrizen.

Answer 8

Sie können dies mit einem Programm genau so lösen, wie Sie es von Hand lösen (mit Multiplikation und Subtraktion, dann führen Sie die Ergebnisse wieder in die Gleichungen ein). Dies ist eine ziemlich normale Mathematik auf Sekundarschulniveau.

-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) 

So Sie am Ende mit:

a = 0.0785250 
b = -0.1806125 
c = -41.6375875 

Wenn Sie diese Werte wieder in A, B und C stecken, finden Sie sie richtig sind.

Der Trick besteht darin, eine einfache 4x3-Matrix zu verwenden, die wiederum zu einer 3x2-Matrix reduziert wird, dann zu einer 2x1, die "a = n" ist, wobei n eine tatsächliche Zahl ist. Sobald Sie das haben, füttern Sie es in die nächste Matrix, um einen anderen Wert zu erhalten, dann diese zwei Werte in die nächste Matrix, bis Sie alle Variablen gelöst haben.

Vorausgesetzt, Sie haben N verschiedene Gleichungen, können Sie immer für N Variablen lösen. Ich sage anders, weil diese beide sind nicht:

7a + 2b = 50 
14a + 4b = 100 

Sie das mit zwei multipliziert gleichen Gleichung sind so können Sie keine Lösung von ihnen erhalten - Multiplikation der ersten von zwei dann Subtrahieren Blätter Sie mit der wahren, aber nutzlos Aussage :

0 = 0 + 0 

---

Als Beispiel, hier einig C-Code, der die simultanen Gleichungen ausarbeitet, die Sie in Ihrer Frage gestellt sind.Zunächst einige notwendige Typen, Variablen, eine Stützfunktion für eine Gleichung auszudrucken, und der Beginn der 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; 

Als nächstes wird die Reduktion der drei Gleichungen mit drei Unbekannten zu zwei Gleichungen mit zwei Unbekannten:

// 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 (""); 

Als nächstes wird die Reduktion der beiden Gleichungen mit zwei unbekannten zu einer Gleichung mit einer unbekannten:

// 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 (""); 

Jetzt, da wir eine Formel vom Typ number1 = unknown * number2 können wir den unbekannten Wert einfach mit unknown <- number1/number2 ausrechnen. Sobald Sie diesen Wert herausgefunden haben, ersetzen Sie ihn in eine der Gleichungen mit zwei Unbekannten und berechnen Sie den zweiten Wert. Dann ersetzen diese beiden (jetzt bekannt) Unbekannten in einer der ursprünglichen Gleichungen und Sie nun die Werte für alle drei Unbekannten haben:

// 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; 
} 

Der Ausgang dieses Codes entspricht den früheren Berechnungen in dieser Antwort:

  >: -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) 

Answer 9

Werfen Sie einen Blick auf die Microsoft Solver Foundation.

Mit ihm können Sie Code wie folgt schreiben:

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); 

Hier ist die Ausgabe:
=== Solver Foundation Service Report ===
für Datum und Uhrzeit: 2009.04.20 23: 29:55
Modellname: Standard
Capabilities angefordert: LP
Lösen Zeit (ms): 1027
Gesamtspielzeit (ms): 1414
LösenCompletion Status: Optimal
Solver ausgewählt: Microsoft.SolverFoundation.Solvers.SimplexSolver
Richtlinien:
Microsoft.SolverFoundation.Services.Directive
Algorithmus: Primal
Arithmetik: Hybrid
Pricing (exakt):
Standard Pricing (double): SteepestEdge
Basis: Slack
Pivot Count: 3
=== Lösungsdetails ===
Ziele:

Entscheidungen:
a: ,0785250000000004
b: -,180612500000001
c: -41,6375875

Answer 10

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