Eine lineare Gleichung lösen
https://stackoverflow.com/questions/769
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich muss ein System linearer Gleichungen in C, Objective C oder (falls erforderlich) C++ programmgesteuert lösen.
Hier ist ein Beispiel für die 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 erhalten a, b, Und tx.
Lösung
Cramers RegelUndGaußsche Eliminierungsind zwei gute Allzweckalgorithmen (siehe auch Simultane lineare Gleichungen).Wenn Sie nach Code suchen, schauen Sie hier vorbei GiNaC, Maxima, Und SymbolicC++ (abhängig natürlich von Ihren Lizenzanforderungen).
BEARBEITEN:Ich weiß, dass Sie in C-Land arbeiten, aber ich muss auch ein gutes Wort dafür einlegen SymPy (ein Computeralgebrasystem in Python).Sie können viel von seinen Algorithmen lernen (wenn Sie sich ein wenig mit Python auskennen).Außerdem steht es unter der neuen BSD-Lizenz, während die meisten kostenlosen Mathematikpakete unter der GPL stehen.
Andere Tipps
Sie können dies mit einem Programm genauso lösen, wie Sie es von Hand lösen (mit Multiplikation und Subtraktion und anschließender Rückführung der Ergebnisse in die Gleichungen).Das ist 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)
Am Ende haben Sie also:
a = 0.0785250
b = -0.1806125
c = -41.6375875
Wenn Sie diese Werte wieder in A, B und C einsetzen, werden Sie feststellen, dass sie korrekt sind.
Der Trick besteht darin, eine einfache 4x3-Matrix zu verwenden, die sich wiederum auf eine 3x2-Matrix und dann auf eine 2x1-Matrix reduziert, was „a = n“ ist, wobei n eine tatsächliche Zahl ist.Sobald Sie das haben, geben Sie es in die nächste Matrix ein, um einen weiteren Wert zu erhalten, und dann fügen Sie diese beiden Werte in die nächste Matrix ein, bis Sie alle Variablen gelöst haben.
Vorausgesetzt, Sie haben N verschiedene Gleichungen, können Sie immer nach N Variablen auflösen.Ich sage unterschiedlich, weil diese beiden nicht sind:
7a + 2b = 50
14a + 4b = 100
Sie sind die Dasselbe Gleichung mit zwei multipliziert, sodass Sie daraus keine Lösung erhalten können. Wenn Sie die erste Gleichung mit zwei multiplizieren und dann subtrahieren, erhalten Sie die wahre, aber nutzlose Aussage:
0 = 0 + 0
Als Beispiel sehen Sie hier einen C-Code, der die simultanen Gleichungen berechnet, die Sie in Ihre Frage einfügen.Zuerst einige notwendige Typen, Variablen, eine Unterstützungsfunktion zum Ausdrucken einer Gleichung und den Beginn von 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 erfolgt die Reduktion der drei Gleichungen mit drei Unbekannten auf 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 erfolgt die Reduktion der beiden Gleichungen mit zwei Unbekannten auf eine 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 haben wir eine Formel dieser Art number1 = unknown * number2, können wir den unbekannten Wert einfach mit berechnen unknown <- number1 / number2.Sobald Sie diesen Wert herausgefunden haben, setzen Sie ihn in eine der Gleichungen mit zwei Unbekannten ein und ermitteln Sie den zweiten Wert.Setzen Sie dann diese beiden (jetzt bekannten) Unbekannten in eine der ursprünglichen Gleichungen ein und Sie haben nun die Werte für alle drei Unbekannten:
// 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;
}
Die Ausgabe 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)
Für ein 3x3-System linearer Gleichungen wäre es meiner Meinung nach in Ordnung, eigene Algorithmen einzuführen.
Möglicherweise müssen Sie sich jedoch Gedanken über die Genauigkeit, die Division durch Null oder sehr kleine Zahlen machen und darüber, was mit unendlich vielen Lösungen zu tun ist.Mein Vorschlag ist, ein Standardpaket für numerische lineare Algebra zu verwenden, z LAPACK.
Werfen Sie einen Blick auf die Microsoft Solver Foundation.
Damit könnte man Code wie diesen 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-Bericht===
Terminzeit:20.04.2009 23:29:55
Modellname:Standard
Gewünschte Fähigkeiten:LP
Lösungszeit (ms):1027
Gesamtzeit (ms):1414
Abschlussstatus lösen:Optimal
Ausgewählter Löser:Microsoft.SolverFoundation.Solvers.SimplexSolver
Richtlinien:
Microsoft.SolverFoundation.Services.Directive
Algorithmus:Ursprünglich
Arithmetik:Hybrid
Preise (genau):Standard
Preis (doppelt):Steilste Kante
Basis:Locker
Pivot-Anzahl:3
===Lösungsdetails===
Ziele:
Entscheidungen:
A:0,0785250000000004
B:-0,180612500000001
C:-41.6375875
Suchen Sie nach einem Softwarepaket, das die Arbeit oder tatsächlich die Matrixoperationen usw. erledigt und jeden Schritt erledigt?
Das erste, ein Kollege von mir hat es gerade benutzt Ocaml GLPK.Es ist nur eine Hülle für die GLPK, aber es entfällt viele Schritte zum Einrichten.Es sieht jedoch so aus, als müssten Sie beim GLPK bleiben, allerdings in C.Für Letzteres danke ich Delicious für die Speicherung eines alten Artikels, den ich vor einiger Zeit über LP gelernt habe. PDF.Wenn Sie bei der weiteren Einrichtung spezielle Hilfe benötigen, lassen Sie es uns wissen. Ich bin sicher, dass ich oder jemand zurückkommen und Ihnen helfen wird, aber ich denke, von hier aus ist es ziemlich einfach.Viel Glück!
Vorlage für numerisches Toolkit von NIST verfügt über Tools dafür.
Eine der zuverlässigeren Methoden ist die Verwendung von a QR-Zerlegung.
Hier ist ein Beispiel für einen Wrapper, mit dem ich „GetInverse(A, InvA)“ in meinem Code aufrufen kann und der die Umkehrung in InvA einfügt.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D ist in der Bibliothek definiert.
Aus dem Wortlaut Ihrer Frage geht hervor, dass Sie mehr Gleichungen als Unbekannte haben und die Inkonsistenzen minimieren möchten.Dies geschieht typischerweise mit einer linearen Regression, die die Summe der Quadrate der Inkonsistenzen minimiert.Abhängig von der Größe der Daten können Sie dies in einer Tabellenkalkulation oder in einem Statistikpaket tun.R ist ein hochwertiges, kostenloses Paket, das unter anderem lineare Regression durchführt.Die lineare Regression hat viel zu bieten (und viele Fallstricke), ist aber für einfache Fälle einfach durchzuführen.Hier ist ein R-Beispiel mit Ihren Daten.Beachten Sie, dass „tx“ der Schnittpunkt 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
Bezüglich der Laufzeiteffizienz haben andere besser geantwortet als ich.Wenn Sie immer die gleiche Anzahl von Gleichungen wie Variablen haben, gefällt mir das Cramers Regel da es einfach umzusetzen ist.Schreiben Sie einfach eine Funktion, um die Determinante einer Matrix zu berechnen (oder verwenden Sie eine bereits geschriebene Funktion, ich bin mir sicher, dass Sie eine finden können) und dividieren Sie die Determinanten zweier Matrizen.
Persönlich bin ich ein Fan der Algorithmen von Numerische Rezepte.(Ich mag die C++-Edition.)
In diesem Buch erfahren Sie, warum die Algorithmen funktionieren, und es zeigt Ihnen einige ziemlich gut debuggte Implementierungen dieser Algorithmen.
Natürlich könnte man es einfach blind verwenden CLAPACK (Ich habe es mit großem Erfolg verwendet), aber ich würde zunächst einen Gaußschen Eliminierungsalgorithmus manuell eingeben, um zumindest eine ungefähre Vorstellung davon zu bekommen, welche Art von Arbeit in die Stabilität dieser Algorithmen gesteckt wurde.
Wenn Sie später interessantere lineare Algebra betreiben, schauen Sie sich den Quellcode von an Oktave wird viele Fragen beantworten.
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
