Rozwiązywanie równania liniowego

głosy
35

Muszę programowo rozwiązać układ równań liniowych C, cel C i (w razie potrzeby), 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ć najlepszą aproksymacji a, boraz tx.

Utwórz 03/08/2008 o 19:14
źródło użytkownik
W innych językach...                            


10 odpowiedzi

głosy
17

Wzory Cramera i eliminacji Gaussa są dwa algorytmy dobra, ogólnego przeznaczenia (patrz również Jednoczesna równań liniowych ). Jeśli szukasz kodu, sprawdź GiNaC , Maxima i SymbolicC ++ (w zależności od wymogów licencyjnych, oczywiście).

EDIT: Wiem, że pracujesz w C ziemi, ale mam również umieścić w dobrym słowem dla SymPy (system komputerowy algebra w Pythonie). Można się wiele nauczyć od swoich algorytmów (jeśli można przeczytać trochę python). Również, to w ramach nowej licencji BSD, podczas gdy większość z darmowych pakietów matematycznych są GPL.

Odpowiedział 03/08/2008 o 19:37
źródło użytkownik

głosy
15

Można rozwiązać ten problem z programem dokładnie w ten sam sposób można rozwiązać go ręcznie (z mnożenia i odejmowania, a następnie z powrotem do karmienia wyników równań). Jest to dość standardowe matematyki szkół średnich poziomu.

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

Trick jest użycie prostego macierzy 4x3, który z kolei redukuje się do macierzy 3x2, 2x1, który następnie jest „= n”, n oznacza wartość liczby. Kiedy już, że karmisz go do następnej matrycy je, aby uzyskać inną wartość, to te dwie wartości do następnej matrycy do góry, dopóki nie rozwiązuje wszystkich zmiennych.

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

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

Są to samo równanie 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 stwierdzeniem:

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 i początku 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, gdy mamy formułę typu number1 = unknown * number2, możemy po prostu wypracować z nieznaną wartość unknown <- number1 / number2. Potem, kiedy już się zorientowali, że wartość, zastąpił go w jednym z równań z dwiema niewiadomymi i wypracować 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)
Odpowiedział 26/02/2009 o 11:59
źródło użytkownik

głosy
7

W przypadku układu równań liniowych 3x3 Chyba byłoby dobrze, aby rozciągnąć swoje własne algorytmy.

Jednakże, może trzeba się martwić o dokładności, dzielenie przez zero lub bardzo małe liczby i co zrobić nieskończenie wiele rozwiązań. Moja propozycja jest taka, aby przejść ze standardowym numeryczne algebry liniowej opakowania takie jak LAPACK .

Odpowiedział 08/08/2008 o 18:59
źródło użytkownik

głosy
6

Spójrz na Solver Fundacji Microsoft .

Dzięki niemu można napisać kod jak poniżej:

  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 wniosek: LP
Rozwiąż czasu (ms): 1027
Całkowity czas (ms): 1414
Rozwiąż Zakończenie status: Optymalne
Solver Wybrane: Microsoft.SolverFoundation.Solvers.SimplexSolver
dyrektyw:
Microsoft.SolverFoundation.Services.Directive
Algorytm: Primal
arytmetyczne: Hybrid
cennik (dokładna): Domyślny
tydzień (podwójne): SteepestEdge
Podstawa: Slack
Pivot Count: 3
== = Szczegóły Rozwiązanie ===
cele:

Decyzje:
a: +0,0785250000000004
b: -0,180612500000001
c: -41.6375875

Odpowiedział 21/04/2009 o 04:33
źródło użytkownik

głosy
3

Template Numerical Toolkit z NIST ma narzędzia do robić.

Jednym z bardziej niezawodnych sposobów jest wykorzystanie rozkładu QR .

Oto przykład z owijki tak, że mogę nazwać „GetInverse (A, INVA)” w moim kodu i będzie można umieścić odwrotność do INVA.

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

Array2D jest zdefiniowana w bibliotece.

Odpowiedział 25/08/2008 o 19:53
źródło użytkownik

głosy
3

Szukasz pakietu oprogramowania, który będzie pracować czy rzeczywiście robi operacje macierzowe i takie i zrobić każdy krok?

The pierwszy współpracownik mój tylko używane SML glpk . To jest po prostu otoki dla glpk , ale usuwa wiele etapów ustawiania rzeczy. Wygląda na to, że będziemy musieli się trzymać z glpk, w C, choć. Dla tych ostatnich, dzięki smaczne do zapisywania starą artykuł kiedyś nauczyć LP chwilę z powrotem, PDF . Jeśli potrzebujesz konfigurowania dalej szczególnej pomocy, daj nam znać i jestem pewien, ja czy ktoś będzie wędrować z powrotem i pomóc, ale myślę, że to dość proste stąd. Powodzenia!

Odpowiedział 03/08/2008 o 19:27
źródło użytkownik

głosy
2

Pod względem wydajności run-time, inni odpowiedział lepiej ode Jeśli zawsze będą miały taką samą liczbę równań jako zmienne, lubię Wzory Cramera jak to jest łatwe do wdrożenia. Wystarczy napisać funkcję, aby obliczyć wyznacznik macierzy (lub użyć jednego, że już napisane, jestem pewien, że można znaleźć tam) i podzielić determinanty dwóch macierzy.

Odpowiedział 19/09/2008 o 04:36
źródło użytkownik

głosy
2

Z treści pytania, wydaje się, że masz więcej równań niż niewiadomych i chcesz zminimalizować niespójności. Zazwyczaj wykonuje się z regresji liniowej, która minimalizuje sumę kwadratów niezgodności. W zależności od rozmiaru danych, można to zrobić w arkuszu kalkulacyjnym lub w pakiecie statystycznym. R jest wysokiej jakości, wolny pakiet, który robi regresji liniowej, wśród wielu innych rzeczy. Jest dużo do regresji liniowej (i dużo Gotcha'S), ale jak to jest proste do zrobienia dla prostych przypadków. Oto przykład R używając swoich danych. Zauważ, że „TX” jest osią do modelu.

> 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  
Odpowiedział 17/09/2008 o 15:22
źródło użytkownik

głosy
1
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
Odpowiedział 30/12/2014 o 15:24
źródło użytkownik

głosy
1

Osobiście jestem częściowe do algorytmów Numerical Recipes . (Jestem lubiący C ++ edycji).

Ta książka nauczy Cię, dlaczego algorytmy pracy, a także pokazać kilka całkiem-dołkowych debugowany implementacje tych algorytmów.

Oczywiście, można po prostu ślepo użyć CLAPACK (Używałem go z wielkim sukcesem), ale najpierw ręcznie wpisać algorytm eliminacji Gaussa przynajmniej mają słabe pojęcie o rodzaju pracy, która poszła do dokonywania tych algorytmów stabilny.

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

Odpowiedział 25/08/2008 o 19:22
źródło użytkownik

Cookies help us deliver our services. By using our services, you agree to our use of cookies. Learn more