math - sistema - Resolviendo una ecuación lineal

resolver ecuaciones lineales

Question

Necesito resolver mediante programación un sistema de ecuaciones lineales en C, Objetivo C o (si es necesario) C ++.

Aquí hay 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 tx .

Answer 1

¿Está buscando un paquete de software que haga el trabajo o realmente haga las operaciones de la matriz y tal y haga cada paso?

El primero, un compañero de trabajo mío acaba de usar Ocaml GLPK . Es solo una envoltura para el GLPK , pero elimina muchos de los pasos para configurar las cosas. Parece que vas a tener que seguir con el GLPK, en C, sin embargo. Para este ú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 la configuración adicional, avísenos y estoy seguro de que alguien o yo volveremos y ayudaremos, pero creo que es bastante sencillo desde aquí. ¡Buena suerte!

Answer 2

Eche un vistazo a la Microsoft Solver Foundation .

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á la salida:
=== Informe del servicio de la Fundación Solver ===
Fecha y hora: 20/04/2009 23:29:55
Nombre del modelo: Predeterminado
Capacidades solicitadas: LP
Tiempo de resolución (ms): 1027
Tiempo total (ms): 1414
Resolver estado de finalización: óptimo
Solucionador seleccionado: Microsoft.SolverFoundation.Solvers.SimplexSolver
Directivas:
Microsoft.SolverFoundation.Services.Directive
Algoritmo: Primal
Aritmética: Híbrida
Precio (exacto): Predeterminado
Precio (doble): SteepestEdge
Base: holgura
Conteo de pivotes: 3
=== Detalles de la solución ===
Metas:

Decisiones:
a: 0.0785250000000004
b: -0.180612500000001
c: -41.6375875

Answer 3

En términos de eficiencia en tiempo de ejecución, otros han respondido mejor que yo. Si siempre tendrá la misma cantidad de ecuaciones que variables, me gusta la regla de Cramer, ya que es fácil de implementar. Simplemente escriba una función para calcular el determinante de una matriz (o use una que ya esté escrita, estoy seguro de que puede encontrarla), y divida los determinantes de dos matrices.

Answer 4

Para un sistema 3x3 de ecuaciones lineales, 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 ir con un paquete estándar de álgebra lineal numérica como LAPACK .

Answer 5

Personalmente, soy parcial a los algoritmos de recetas numéricas . (Soy aficionado a la edición C ++).

Este libro le enseñará por qué funcionan los algoritmos, y le mostrará algunas implementaciones depuradas bastante bien de esos algoritmos.

Por supuesto, podría usar CLAPACK ciegas (lo he usado con gran éxito), pero primero escribiría a mano un algoritmo de eliminación gaussiana para al menos tener una idea tenue del tipo de trabajo que se ha dedicado a hacer estos algoritmos estable.

Más tarde, si está haciendo un álgebra lineal más interesante, mirar alrededor del código fuente de Octave responderá muchas preguntas.

Answer 6

Puede resolver esto con un programa exactamente de la misma manera que lo resuelve a mano (con multiplicación y sustracción, luego ingresando los resultados en las ecuaciones). Esta es una matemática bastante estándar de 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 conectar estos valores a A, B y C, encontrará que son correctos.

El truco es usar una matriz simple de 4x3 que se reduce a su vez a una matriz de 3x2, luego un 2x1 que es "a = n", siendo n un número real. Una vez que tiene 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 puede resolver N variables. Digo distinto porque estos dos no son:

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

Son la misma ecuación multiplicada por dos, por lo que no puede obtener una solución: multiplicar la primera por dos y luego restar le 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 se colocan en su 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;

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

Luego, 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 unknown <- number1 / number2 con unknown <- number1 / number2 . Luego, una vez que haya calculado ese valor, sustitúyalo en una de las ecuaciones con dos incógnitas y calcule el segundo valor. Luego sustituya ambas incógnitas (ahora conocidas) en una de las ecuaciones originales y ahora tiene los valores para 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;
}

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

Answer 7

Según la redacción de su pregunta, parece que tiene más ecuaciones que incógnitas y desea minimizar las inconsistencias. Esto se hace típicamente con regresión lineal, que minimiza la suma de los cuadrados de las inconsistencias. Dependiendo del tamaño de los datos, puede hacerlo en una hoja de cálculo o en un paquete estadístico. R es un paquete gratuito de alta calidad que hace regresión lineal, entre muchas otras cosas. Hay mucho en la regresión lineal (y muchos problemas), pero como es fácil de hacer para casos simples. Aquí hay un ejemplo de R con sus datos. Tenga en cuenta que el "tx" es la intercepció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  

Answer 8

La regla de Cramer y la eliminación gaussiana son dos buenos algoritmos de propósito general (también vea Ecuaciones lineales simultáneas ). Si está buscando código, consulte GiNaC , Maxima y SymbolicC++ (dependiendo de sus requisitos de licencia, por supuesto).

EDITAR: Sé que estás trabajando en C land, pero también tengo que poner una buena palabra para SymPy (un sistema de álgebra de computadora 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.

Answer 9

Template Numerical Toolkit de NIST tiene herramientas para hacerlo.

Una de las formas más confiables es usar una descomposición QR .

Aquí hay un ejemplo de un contenedor para que pueda llamar "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 se define en la biblioteca.

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