Resolviendo una ecuación lineal

Question

Necesito resolver programáticamente un sistema de ecuaciones lineales en C, Objetivo C, o (si es necesario) C++.

He aquí 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 

partir de esto, me gustaría obtener la mejor aproximación para a, b, y tx.

Answer 1

Cramer's Rule y Gaussian Elimination son dos algoritmos bueno, de propósito general (también ver Simultaneous Linear Equations). Si está buscando un código, consulte GiNaC, Maxima y SymbolicC++ (según sus requisitos de licencia, por supuesto).

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

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

El primero, un compañero de trabajo que acabo de utilizar Ocaml GLPK. Es solo un contenedor para el GLPK, pero elimina muchos de los pasos de configuración. Sin embargo, parece que tendrás que quedarte con el GLPK, en C. Para este último, gracias a delicioso para guardar un artículo antiguo solía aprender LP hace un tiempo atrás, PDF. Si necesita ayuda específica para la configuración, infórmenos y estoy seguro de que yo o alguien volveremos y lo ayudaremos, pero creo que es bastante directo desde aquí. ¡Buena suerte!

Answer 3

Para un sistema 3x3 de ecuaciones lineales, supongo que sería correcto implementar sus propios algoritmos.

Sin embargo, puede que tenga que preocuparse por la precisión, división por cero o números realmente pequeños y qué hacer con infinitas soluciones. Mi sugerencia es ir con un paquete de álgebra lineal numérico estándar como LAPACK.

Answer 4

Personalmente, soy parcial a los algoritmos de Numerical Recipes. (Soy aficionado a la edición en C++).

Este libro le enseñará por qué los algoritmos funcionan, además de mostrarle algunas implementaciones bastante depuradas de esos algoritmos.

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

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

Answer 5

Template Numerical Toolkit de NIST tiene herramientas para hacerlo.

Una de las maneras más confiables es usar un QR Decomposition.

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 se define en la biblioteca.

Answer 6

Según la redacción de su pregunta, parece que tiene más ecuaciones que incógnitas y desea minimizar las incoherencias. 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 hacer esto en una hoja de cálculo o en un paquete estadístico. R es un paquete gratuito de alta calidad que realiza regresiones lineales, entre muchas otras cosas. Hay mucho en la regresión lineal (y muchos errores), pero es fácil de hacer en casos simples. Aquí hay un ejemplo R usando sus datos. Tenga en cuenta que el "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 

Answer 7

En términos de eficiencia en tiempo de ejecución, otros han respondido mejor que yo Si usted siempre tendrá el mismo número de ecuaciones como variables, me gusta Cramer's rule 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 encontrar una) y divida los determinantes de dos matrices.

Answer 8

Puede resolver esto con un programa exactamente de la misma manera que lo resuelve a mano (con multiplicación y sustracción, y luego devolviendo los resultados a las ecuaciones). Esto es bastante estándar en matemáticas de nivel secundario.

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

por lo que terminan con:

a = 0.0785250 
b = -0.1806125 
c = -41.6375875 

Si conecta estos valores de nuevo en A, B y C, se encontrará que son correctos.

El truco es usar una matriz 4x3 simple que a su vez se reduce a una matriz de 3x2, luego un 2x1 que es "a = n", siendo n un número real. Una vez que tienes eso, lo pasas a la siguiente matriz para obtener otro valor, luego esos dos valores en la siguiente matriz hasta que hayas resuelto todas las variables.

Siempre que tenga N ecuaciones distintas, siempre puede resolver para N variables. Digo distinta porque estos dos no son:

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

Ellos son la ecuación misma multiplicado por dos, así que no se puede obtener una solución de ellos - la multiplicación de la primera por dos y luego restando te deja con la afirmación verdadera, pero inútil :

0 = 0 + 0 

---

a modo de ejemplo, aquí hay algo de código C que se resuelve el sistema de ecuaciones que está colocado en su pregunta.Primeros algunos tipos necesarios, variables, una función de soporte para la impresión de 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 de el 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 resuelve el segundo valor. A continuación, sustituir esas dos incógnitas (ahora conocidos) 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 9

Eche un vistazo al Microsoft Solver Foundation.

Con él podría 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:
=== Solver Foundation Service Informe ===
de fecha y hora: 04/20/2009 23: Nombre 29:55
Modelo: Por defecto
capacidades solicitadas: LP
Resolver Tiempo (ms): 1027
tiempo total (ms): 1414
Solución de estado de finalización: Óptimo
Solver seleccionada: Microsoft.SolverFoundation.Solvers.SimplexSolver
Directivas:
Microsoft.SolverFoundation.Services.Directive
algoritmo: Primal
aritmético: híbrido
Precio: (exacta): Por defecto
precios (doble): SteepestEdge
Base: Slack
pivote Count: 3
=== === Detalles Solución
Objetivos:

Decisiones:
a: ,0785250000000004
b: -0.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