Résoudre une équation linéaire
https://stackoverflow.com/questions/769
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Je dois résoudre par programme un système d'équations linéaires en C, Objective C ou (si nécessaire) C++.
Voici un exemple des équations :
-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
À partir de là, j'aimerais obtenir la meilleure approximation pour a, b, et tx.
La solution
La règle de CrameretÉlimination gaussiennesont deux bons algorithmes à usage général (voir aussi Équations linéaires simultanées).Si vous cherchez du code, consultez GiNaC, Maxime, et SymboliqueC++ (en fonction de vos exigences en matière de licence, bien sûr).
MODIFIER:Je sais que vous travaillez en Terre C, mais je dois aussi dire un bon mot pour SymPy (un système de calcul formel en Python).Vous pouvez apprendre beaucoup de ses algorithmes (si vous savez lire un peu de python).De plus, il est sous la nouvelle licence BSD, alors que la plupart des packages mathématiques gratuits sont sous GPL.
Autres conseils
Vous pouvez résoudre ce problème avec un programme exactement de la même manière que vous le résolvez à la main (avec multiplication et soustraction, puis réintégration des résultats dans les équations).Il s’agit de mathématiques assez classiques au niveau secondaire.
-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)
Vous vous retrouvez donc avec :
a = 0.0785250
b = -0.1806125
c = -41.6375875
Si vous rebranchez ces valeurs dans A, B et C, vous constaterez qu'elles sont correctes.
L'astuce consiste à utiliser une simple matrice 4x3 qui se réduit à son tour à une matrice 3x2, puis à un 2x1 qui vaut "a = n", n étant un nombre réel.Une fois que vous avez cela, vous l'introduisez dans la matrice suivante pour obtenir une autre valeur, puis ces deux valeurs dans la matrice suivante jusqu'à ce que vous ayez résolu toutes les variables.
À condition que vous ayez N équations distinctes, vous pouvez toujours résoudre N variables.Je dis distinct car ces deux-là ne le sont pas :
7a + 2b = 50
14a + 4b = 100
Ils sont les même équation multipliée par deux afin que vous ne puissiez pas en obtenir une solution - multiplier le premier par deux puis soustraire vous laisse avec l'énoncé vrai mais inutile :
0 = 0 + 0
À titre d'exemple, voici du code C qui résout les équations simultanées que vous avez placées dans votre question.D'abord quelques types nécessaires, des variables, une fonction de support pour imprimer une équation, et le début de 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;
Ensuite, la réduction des trois équations à trois inconnues à deux équations à deux inconnues :
// 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 ("");
Ensuite, la réduction des deux équations à deux inconnues à une équation à une inconnue :
// 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 ("");
Maintenant que nous avons une formule du type number1 = unknown * number2, nous pouvons simplement calculer la valeur inconnue avec unknown <- number1 / number2.Ensuite, une fois que vous avez déterminé cette valeur, remplacez-la dans l'une des équations à deux inconnues et calculez la deuxième valeur.Remplacez ensuite ces deux inconnues (maintenant connues) dans l’une des équations d’origine et vous avez maintenant les valeurs des trois inconnues :
// 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;
}
La sortie de ce code correspond aux calculs précédents dans cette réponse :
>: -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)
Pour un système 3x3 d'équations linéaires, je suppose que ce serait bien de déployer vos propres algorithmes.
Cependant, vous devrez peut-être vous soucier de l'exactitude, de la division par zéro ou de très petits nombres et de ce qu'il faut faire avec une infinité de solutions.Ma suggestion est d'opter pour un package d'algèbre linéaire numérique standard tel que LAPACK.
Jetez un oeil à Fondation Microsoft Solver.
Avec cela, vous pourriez écrire du code comme celui-ci :
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);
Voici le résultat :
===Rapport du service Solver Foundation===
Dateheure :20/04/2009 23:29:55
Nom du modèle:Défaut
Capacités demandées :LP
Temps de résolution (ms) :1027
Temps total (ms) :1414
Résoudre l'état d'achèvement :Optimal
Solveur sélectionné :Microsoft.SolverFoundation.Solvers.SimplexSolver
Directives :
Microsoft.SolverFoundation.Services.Directive
Algorithme:Primitif
Arithmétique:Hybride
Prix (exact) :Défaut
Tarif (double) :Bord le plus raide
Base:Mou
Nombre de pivots :3
===Détails de la solution===
Objectifs:
Les décisions:
un:0,0785250000000004
B :-0,180612500000001
c :-41.6375875
Êtes-vous à la recherche d'un progiciel qui fera le travail ou effectuera réellement les opérations matricielles et autres et effectuera chaque étape ?
Le premier, un de mes collègues vient d'utiliser Ocaml GLPK.C'est juste un emballage pour le GLPK, mais cela supprime de nombreuses étapes de configuration.Il semble que vous devrez vous en tenir au GLPK, en C.Pour ce dernier, merci à délicieux d'avoir sauvegardé un vieil article que j'avais l'habitude d'apprendre le LP il y a quelque temps, PDF.Si vous avez besoin d'une aide spécifique pour vous installer davantage, faites-le-nous savoir et je suis sûr que moi ou quelqu'un reviendra vous aider, mais je pense que c'est assez simple à partir de là.Bonne chance!
Boîte à outils numérique modèle du NIST dispose d'outils pour ce faire.
L'un des moyens les plus fiables consiste à utiliser un Décomposition QR.
Voici un exemple de wrapper pour que je puisse appeler "GetInverse(A, InvA)" dans mon code et qu'il mettra l'inverse dans InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D est défini dans la bibliothèque.
D'après la formulation de votre question, il semble que vous ayez plus d'équations que d'inconnues et que vous souhaitiez minimiser les incohérences.Cela se fait généralement par régression linéaire, qui minimise la somme des carrés des incohérences.Selon la taille des données, vous pouvez le faire dans une feuille de calcul ou dans un logiciel statistique.R est un package gratuit de haute qualité qui effectue, entre autres, une régression linéaire.Il y a beaucoup de choses à faire dans la régression linéaire (et beaucoup de pièges), mais c'est simple à faire pour des cas simples.Voici un exemple R utilisant vos données.Notez que le "tx" est l'interception de votre modèle.
> 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
En termes d'efficacité d'exécution, d'autres ont mieux répondu que moi.Si vous avez toujours le même nombre d'équations que de variables, j'aime bien La règle de Cramer car il est facile à mettre en œuvre.Écrivez simplement une fonction pour calculer le déterminant d'une matrice (ou utilisez-en une déjà écrite, je suis sûr que vous pouvez en trouver une) et divisez les déterminants de deux matrices.
Personnellement, je suis partisan des algorithmes de Recettes numériques.(J'aime l'édition C++.)
Ce livre vous apprendra pourquoi les algorithmes fonctionnent et vous montrera quelques implémentations assez bien déboguées de ces algorithmes.
Bien sûr, vous pouvez simplement utiliser aveuglément CLAPACK (Je l'ai utilisé avec beaucoup de succès), mais je taperais d'abord un algorithme d'élimination gaussienne pour avoir au moins une vague idée du type de travail qui a été consacré à rendre ces algorithmes stables.
Plus tard, si vous faites de l'algèbre linéaire plus intéressante, en parcourant le code source de Octave répondra à beaucoup de questions.
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
