Resolvendo uma equação linear
https://stackoverflow.com/questions/769
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPergunta
Preciso resolver programaticamente um sistema de equações lineares em C, Objective C ou (se necessário) C++.
Aqui está um exemplo das equações:
-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 disso, gostaria de obter a melhor aproximação para a, b, e tx.
Solução
Regra de CramereEliminação gaussianasão dois bons algoritmos de uso geral (veja também Equações Lineares Simultâneas).Se você está procurando por código, confira GiNaC, Máxima, e SimbólicoC++ (dependendo dos seus requisitos de licenciamento, é claro).
EDITAR:Eu sei que você está trabalhando no país C, mas também tenho que dar uma boa palavra para SymPy (um sistema de álgebra computacional em Python).Você pode aprender muito com seus algoritmos (se souber ler um pouco de python).Além disso, está sob a nova licença BSD, enquanto a maioria dos pacotes matemáticos gratuitos são GPL.
Outras dicas
Você pode resolver isso com um programa exatamente da mesma maneira que resolve manualmente (com multiplicação e subtração e, em seguida, inserindo os resultados nas equações).Esta é uma matemática bastante comum no nível do ensino médio.
-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)
Então você acaba com:
a = 0.0785250
b = -0.1806125
c = -41.6375875
Se você inserir esses valores novamente em A, B e C, verá que eles estão corretos.
O truque é usar uma matriz 4x3 simples que, por sua vez, se reduz a uma matriz 3x2, depois a 2x1 que é "a = n", sendo n um número real.Depois de ter isso, você o alimenta na próxima matriz para obter outro valor e, em seguida, esses dois valores na próxima matriz até resolver todas as variáveis.
Desde que você tenha N equações distintas, você sempre poderá resolver N variáveis.Digo distinto porque esses dois não são:
7a + 2b = 50
14a + 4b = 100
Eles são o mesmo equação multiplicada por dois para que você não possa obter uma solução delas - multiplicar a primeira por dois e depois subtrair deixa você com a afirmação verdadeira, mas inútil:
0 = 0 + 0
A título de exemplo, aqui está um código C que resolve as equações simultâneas que você colocou em sua pergunta.Primeiro, alguns tipos necessários, variáveis, uma função de suporte para imprimir uma equação e o início da 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 seguir, a redução das três equações com três incógnitas para duas equações com duas 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 seguir, a redução das duas equações com duas incógnitas para uma equação com uma 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 ("");
Agora que temos uma fórmula do tipo number1 = unknown * number2, podemos simplesmente calcular o valor desconhecido com unknown <- number1 / number2.Então, depois de descobrir esse valor, substitua-o em uma das equações com duas incógnitas e calcule o segundo valor.Em seguida, substitua ambas as incógnitas (agora conhecidas) em uma das equações originais e agora você terá os valores para todas as três 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;
}
A saída desse código corresponde aos cálculos anteriores nesta resposta:
>: -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 um sistema 3x3 de equações lineares, acho que não há problema em implementar seus próprios algoritmos.
No entanto, você pode ter que se preocupar com a precisão, a divisão por zero ou números realmente pequenos e o que fazer com um número infinito de soluções.Minha sugestão é usar um pacote padrão de álgebra linear numérica, como LAPACK.
Dê uma olhada no Fundação Microsoft Solver.
Com ele você poderia escrever um 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);
Aqui está a saída:
===Relatório de serviço do Solver Foundation===
Data hora:20/04/2009 23:29:55
Nome do modelo:Padrão
Capacidades solicitadas:LP
Tempo de resolução (ms):1027
Tempo total (ms):1414
Resolva o status de conclusão:Ótimo
Solucionador selecionado:Microsoft.SolverFoundation.Solvers.SimplexSolver
Diretivas:
Microsoft.SolverFoundation.Services.Directive
Algoritmo:Primitivo
Aritmética:Híbrido
Preço (exato):Padrão
Preço (duplo):Borda mais íngreme
Base:Folga
Contagem dinâmica:3
===Detalhes da solução===
Metas:
Decisões:
a:0,0785250000000004
b:-0,180612500000001
c:-41.6375875
Você está procurando um pacote de software que faça o trabalho ou realmente faça as operações da matriz e outras coisas e execute cada etapa?
O primeiro, um colega meu acabou de usar Ocaml GLPK.É apenas um invólucro para o GLPK, mas remove muitas etapas de configuração.Parece que você terá que ficar com o GLPK, em C.Por este último, obrigado ao delicioso por salvar um artigo antigo que aprendi LP há algum tempo, PDF.Se você precisar de ajuda específica para configurar ainda mais, avise-nos e tenho certeza de que eu ou alguém voltaremos e ajudaremos, mas acho que é bastante simples a partir daqui.Boa sorte!
Kit de ferramentas numéricas de modelo do NIST tem ferramentas para fazer isso.
Uma das maneiras mais confiáveis é usar um Decomposição QR.
Aqui está um exemplo de wrapper para que eu possa chamar "GetInverse(A, InvA)" em meu código e colocar o inverso em InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D é definido na biblioteca.
Pela redação da sua pergunta, parece que você tem mais equações do que incógnitas e deseja minimizar as inconsistências.Isso normalmente é feito com regressão linear, que minimiza a soma dos quadrados das inconsistências.Dependendo do tamanho dos dados, você pode fazer isso em uma planilha ou em um pacote estatístico.R é um pacote gratuito de alta qualidade que faz regressão linear, entre muitas outras coisas.Há muito na regressão linear (e muitas pegadinhas), mas é simples de fazer em casos simples.Aqui está um exemplo R usando seus dados.Observe que o "tx" é a interceptação do seu 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
Em termos de eficiência de tempo de execução, outros responderam melhor do que eu.Se você sempre terá o mesmo número de equações que variáveis, eu gosto Regra de Cramer pois é fácil de implementar.Basta escrever uma função para calcular o determinante de uma matriz (ou usar uma que já esteja escrita, tenho certeza que você encontrará uma por aí) e dividir os determinantes de duas matrizes.
Pessoalmente, gosto dos algoritmos de Receitas Numéricas.(Gosto da edição C++.)
Este livro ensinará por que os algoritmos funcionam, além de mostrar algumas implementações bem depuradas desses algoritmos.
Claro, você poderia usar cegamente CLAPACK (Eu o usei com grande sucesso), mas eu digitaria primeiro um algoritmo de eliminação gaussiana para ter pelo menos uma vaga ideia do tipo de trabalho realizado para tornar esses algoritmos estáveis.
Mais tarde, se você estiver fazendo álgebra linear mais interessante, olhando o código-fonte do Oitava responderá muitas perguntas.
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
