Rozwiązywanie równania liniowego

Question

Potrzebuję programowo rozwiązać układ równań liniowych w C, Cel C lub (jeśli potrzeba) C++.

Oto przykład z równań:

-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 

Z tego, chciałbym, aby uzyskać najlepsze przybliżenie dla a, b i tx.

Answer 1

Cramer's Rule i Gaussian Elimination dwa algorytmy dobrze ogólnego zastosowania (patrz również Simultaneous Linear Equations). Jeśli szukasz kodu, sprawdź: GiNaC, Maxima i SymbolicC++ (w zależności od twoich wymagań licencyjnych, oczywiście).

EDYCJA: Wiem, że pracujesz w kraju C, ale muszę też napisać dobre słowo dla SymPy (system algebry komputerowej w Pythonie). Możesz wiele nauczyć się z jego algorytmów (jeśli potrafisz odczytać trochę pythona). Ponadto jest objęty nową licencją BSD, podczas gdy większość darmowych pakietów matematycznych to GPL.

Answer 2

Poszukujesz pakietu oprogramowania, który wykona pracę lub faktycznie wykona operacje macierzy i takie i zrobi każdy krok?

Pierwszy, mój współpracownik właśnie użył Ocaml GLPK. To jest tylko opakowanie dla GLPK, ale usuwa wiele czynności związanych z konfiguracją. Wygląda jednak na to, że będziesz musiał trzymać się GLPK, ale w C. Dla tego ostatniego, dzięki pysznym do zapisania starego artykułu, nauczyłem się LP z powrotem, PDF. Jeśli potrzebujesz konkretnej pomocy w konfiguracji, daj nam znać, a jestem pewien, że ja lub ktoś wróci z powrotem i pomożemy, ale myślę, że jest to dość proste. Powodzenia!

Answer 3

Dla układu równań liniowych 3x3, wydaje mi się, że byłoby dobrze wprowadzić własne algorytmy.

Możesz jednak martwić się o dokładność, dzielenie przez zero lub bardzo małe liczby i co zrobić z nieskończenie wieloma rozwiązaniami. Moja sugestia dotyczy standardowego pakietu liczbowej algebry liniowej, na przykład LAPACK.

Answer 4

Osobiście jestem stronniczy od algorytmów Numerical Recipes. (Lubię wersję C++).

Ta książka nauczy Cię, dlaczego działają algorytmy, oraz pokazuje całkiem dobrze zdebugowane implementacje tych algorytmów.

Oczywiście można po prostu użyć na ślepo CLAPACK (użyłem go z wielkim sukcesem), ale najpierw ręcznie opracowałem algorytm gaussowskich eliminacji, aby przynajmniej mieć słabe pojęcie o rodzaju pracy, która zniknęła do stabilizacji tych algorytmów.

Później, jeśli robisz bardziej interesującą algebrę liniową, rozglądanie się po kodzie źródłowym Octave odpowie na wiele pytań.

Answer 5

Template Numerical Toolkit z NIST ma narzędzia do robienia tego.

Jednym z bardziej niezawodnych sposobów jest użycie QR Decomposition.

Oto przykład opakowania, dzięki któremu mogę wywołać w moim kodzie "GetInverse (A, InvA)", co spowoduje odwrotność do InvA.

void GetInverse(const Array2D<double>& A, Array2D<double>& invA) 
    { 
    QR<double> qr(A); 
    invA = qr.solve(I); 
    } 

Array2D jest zdefiniowany w bibliotece.

Answer 6

Z treści Twojego pytania wynika, że ​​masz więcej równań niż niewiadomych i chcesz zminimalizować niespójności. Zazwyczaj dokonuje się tego za pomocą regresji liniowej, która minimalizuje sumę kwadratów niespójności. W zależności od rozmiaru danych można to zrobić w arkuszu kalkulacyjnym lub w pakiecie statystycznym. R to wysokiej jakości, bezpłatny pakiet, który dokonuje regresji liniowej, między wieloma innymi rzeczami. Istnieje wiele regresji liniowej (i dużo błędów), ale jak to jest proste w prostych przypadkach. Oto przykład R z wykorzystaniem twoich danych. Zwróć uwagę, że "tx" jest punktem przecięcia z twoim modelem.

> 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

Pod względem wydajności run-time, inni odpowiedział lepiej niż I. Jeśli zawsze będą miały taką samą liczbę równań jak zmienne, lubię Cramer's rule jak to jest łatwe do wdrożenia. Wystarczy napisać funkcję, aby obliczyć wyznacznik macierzy (lub użyć takiego, który już został napisany, na pewno można go znaleźć) i podzielić determinanty dwóch macierzy.

Answer 8

Możesz rozwiązać ten problem za pomocą programu dokładnie w ten sam sposób, w jaki rozwiązałeś go ręcznie (mnożenie i odejmowanie, a następnie podawanie wyników do równania). Jest to dość standardowa matematyka na poziomie szkoły średniej.

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

więc skończyć z:

a = 0.0785250 
b = -0.1806125 
c = -41.6375875 

Jeśli podłączysz te wartości z powrotem do A, B i C, a przekonasz się, że są prawidłowe.

Sztuką jest użycie prostej macierzy 4x3, która z kolei zmniejsza z kolei do matrycy 3x2, a następnie 2x1, która jest "a = n", gdzie n jest liczbą rzeczywistą. Kiedy już to zrobisz, wprowadzasz go do następnej macierzy, aby uzyskać kolejną wartość, następnie te dwie wartości do następnej macierzy, aż rozwiążesz wszystkie zmienne.

Pod warunkiem, że masz N różnych równań, zawsze możesz rozwiązać N zmiennych. Mówię wyraźny, ponieważ te dwa nie są:

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

Są równanie samo pomnożyć przez dwa, więc nie można uzyskać rozwiązanie z nimi - mnożąc pierwsze dwa potem odejmowanie pozostawia nam z prawdziwego ale bezużytecznego rachunku :

0 = 0 + 0 

---

tytułem przykładu, oto niektóre kod C, która działa na zewnątrz równań że jesteś umieszczone w swoim pytaniu.Pierwsze pewne niezbędne typy, zmienne, funkcji wsparcia dla drukowania równanie, a początkiem 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; 

Następnie redukcja trzech równań z trzema niewiadomymi do dwóch równań z dwiema niewiadomymi:

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

Następnie redukcja dwóch równań z dwiema niewiadomymi do jednego równania z jedną niewiadomą:

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

teraz, że mamy formułę typ number1 = unknown * number2, możemy po prostu wypracować nieznaną wartość za pomocą unknown <- number1/number2. Następnie, gdy już zorientujesz się, że wartość się zmieniła, podstaw ją na jedno z równań z dwoma niewiadomymi i wypracuj drugą wartość. Następnie zastąpić oba te (obecnie znane) niewiadomych w jednej z oryginalnych wzorów i teraz masz wartości dla wszystkich trzech niewiadomych:

// 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; 
} 

Wyjście z tego kodu pasuje do wcześniejszych obliczeń w tej odpowiedzi:

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

Spójrz na Microsoft Solver Foundation.

Dzięki niemu można napisać kod tak:

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

Oto wynik:
=== Solver Foundation Obsługa Zgłoś ===
Datetime: 20.04.2009 23: 29:55
Nazwa modelu: domyślne
Możliwości wnioskowane: LP
Rozwiąż czasu (ms): 1027
Całkowity czas (ms): 1414
Rozwiąż Zakończenie Status: Optimal
Solver Wybrano: Microsoft.SolverFoundation.Solvers.SimplexSolver
dyrektywach:
Microsoft.SolverFoundation.Services.Directive
algorytmu: Primal
arytmetyczne: Hybrid
cenowe (dokładna): Domyślnie
tydzień (podwójne): SteepestEdge
Podstawa: Slack
Pivot Count: 3
=== Szczegóły Rozwiązanie ===
cele:

Decyzje:
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