Resolver una ecuación lineal
https://stackoverflow.com/questions/769
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
Necesito resolver mediante programación un sistema de ecuaciones lineales en C, Objective C o (si es necesario) C++.
A continuación se muestra un ejemplo de las ecuaciones:
-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
A partir de esto, me gustaría obtener la mejor aproximación para a, b, y tx.
Solución
La regla de CrameryEliminación gaussianaHay dos buenos algoritmos de propósito general (ver también Ecuaciones lineales simultáneas).Si estás buscando un código, consulta GiNaC, máxima, y SimbólicoC++ (dependiendo de los requisitos de su licencia, por supuesto).
EDITAR:Sé que estás trabajando en tierra C, pero también tengo que hablar bien de ti. SymPy (un sistema de álgebra informática en Python).Puedes aprender mucho de sus algoritmos (si puedes leer un poco de Python).Además, está bajo la nueva licencia BSD, mientras que la mayoría de los paquetes matemáticos gratuitos son GPL.
Otros consejos
Puedes resolver esto con un programa exactamente de la misma manera que lo resuelves a mano (con multiplicación y resta, y luego reintroduciendo los resultados en las ecuaciones).Estas son matemáticas bastante estándar de nivel de escuela secundaria.
-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)
Entonces terminas con:
a = 0.0785250
b = -0.1806125
c = -41.6375875
Si vuelve a introducir estos valores en A, B y C, encontrará que son correctos.
El truco consiste en utilizar una matriz simple de 4x3 que a su vez se reduce a una matriz de 3x2 y luego a 2x1 que es "a = n", siendo n un número real.Una vez que tenga eso, lo introduce en la siguiente matriz para obtener otro valor, luego esos dos valores en la siguiente matriz hasta que haya resuelto todas las variables.
Siempre que tenga N ecuaciones distintas, siempre podrá resolver N variables.Digo distintos porque estos dos no lo son:
7a + 2b = 50
14a + 4b = 100
Ellos son las mismo ecuación multiplicada por dos, por lo que no puedes obtener una solución de ellas; multiplicar la primera por dos y luego restar te deja con la afirmación verdadera pero inútil:
0 = 0 + 0
A modo de ejemplo, aquí hay un código C que resuelve las ecuaciones simultáneas que colocaste en tu pregunta.Primero, algunos tipos necesarios, variables, una función de soporte para imprimir una ecuación y el inicio 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;
A continuación, la reducción de las tres ecuaciones con tres incógnitas a dos ecuaciones con dos incógnitas:
// 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 ("");
A continuación, la reducción de las dos ecuaciones con dos incógnitas a una ecuación con una incógnita:
// 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 ("");
Ahora que tenemos una fórmula del tipo number1 = unknown * number2, simplemente podemos calcular el valor desconocido con unknown <- number1 / number2.Luego, una vez que hayas calculado ese valor, sustitúyelo en una de las ecuaciones con dos incógnitas y calcula el segundo valor.Luego sustituya ambas incógnitas (ahora conocidas) en una de las ecuaciones originales y ahora tendrá los valores de las tres incógnitas:
// 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;
}
El resultado de ese código coincide con los cálculos anteriores en esta respuesta:
>: -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)
Para un sistema de ecuaciones lineales de 3x3, supongo que estaría bien implementar sus propios algoritmos.
Sin embargo, es posible que deba preocuparse por la precisión, la división por cero o números realmente pequeños y qué hacer con infinitas soluciones.Mi sugerencia es optar por un paquete estándar de álgebra lineal numérica como LAPACK.
Échale un vistazo al Fundación Microsoft Solver.
Con él podrías escribir código como este:
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);
Aquí está el resultado:
===Informe de servicio de Solver Foundation===
Fecha y hora:20/04/2009 23:29:55
Nombre del modelo:Por defecto
Capacidades solicitadas:LP
Tiempo de resolución (ms):1027
Tiempo total (ms):1414
Estado de finalización de la resolución:Óptimo
Solucionador seleccionado:Microsoft.SolverFoundation.Solvers.SimplexSolver
Directivas:
Directiva.de.servicios.de.Microsoft.SolverFoundation.
Algoritmo:Primitivo
Aritmética:Híbrido
Precio (exacto):Por defecto
Precio (doble):Borde más empinado
Base:Flojo
Conteo de pivotes:3
===Detalles de la solución===
Objetivos:
Decisiones:
a:0.0785250000000004
b:-0,180612500000001
C:-41.6375875
¿Está buscando un paquete de software que haga el trabajo o que realmente haga las operaciones matriciales y demás y realice cada paso?
El primero, lo acaba de usar un compañero de trabajo. Ocaml GLPK.Es sólo un envoltorio para el GLPK, pero elimina muchos de los pasos necesarios para configurar las cosas.Sin embargo, parece que tendrás que seguir con el GLPK en C.Para esto último, gracias a Delicious por guardar un artículo antiguo que solía aprender LP hace un tiempo. PDF.Si necesita ayuda específica para seguir configurando, háganoslo saber y estoy seguro de que alguien o yo volveremos a entrar y le ayudaremos, pero creo que es bastante sencillo a partir de aquí.¡Buena suerte!
Kit de herramientas numéricas de plantilla del NIST tiene herramientas para hacerlo.
Una de las formas más confiables es utilizar un Descomposición QR.
Aquí hay un ejemplo de un contenedor para que pueda llamar a "GetInverse(A, InvA)" en mi código y colocará el inverso en InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D está definido en la biblioteca.
Por la redacción de su pregunta, parece que tiene más ecuaciones que incógnitas y desea minimizar las inconsistencias.Esto normalmente se hace con regresión lineal, que minimiza la suma de los cuadrados de las inconsistencias.Dependiendo del tamaño de los datos, puedes hacerlo en una hoja de cálculo o en un paquete estadístico.R es un paquete gratuito de alta calidad que realiza regresión lineal, entre muchas otras cosas.Hay mucho en la regresión lineal (y muchas trampas), pero es sencillo de hacer en casos simples.Aquí hay un ejemplo de R usando sus datos.Tenga en cuenta que "tx" es la intersección de su modelo.
> 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 términos de eficiencia en tiempo de ejecución, otros han respondido mejor que yo.Si siempre tendrás el mismo número de ecuaciones que variables, me gusta regla de cramer ya que es fácil de implementar.Simplemente escribe una función para calcular el determinante de una matriz (o usa una que ya esté escrita, estoy seguro de que puedes encontrar una) y divide los determinantes de dos matrices.
Personalmente, soy partidario de los algoritmos de Recetas Numéricas.(Me gusta la edición C++).
Este libro le enseñará por qué funcionan los algoritmos y además le mostrará algunas implementaciones bastante bien depuradas de esos algoritmos.
Por supuesto, puedes usar a ciegas CLAPACK (Lo he usado con gran éxito), pero primero escribiría un algoritmo de eliminación gaussiana para al menos tener una idea vaga del tipo de trabajo que se ha realizado para estabilizar estos algoritmos.
Más adelante, si estás haciendo álgebra lineal más interesante, revisa el código fuente de Octava responderá muchas preguntas.
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
