Rozwiązywanie układów równań liniowych

Zagadnienie

Niech zadane będą A = (a_ij) - macierz n razy n, i b = (b_j) - wektor o długości n. Szukamy rozwiązanie x = (x_j) układu równań Ax = b, czyli

Wiemy, że jeśli , to istnieje jednoznaczne rozwiązanie powyższego układu.

Wzór Kramera

Rozwiązanie można przedstawić wzórem Kramera:

gdzie

czyli i-ta kolumna A zastąpiona została przez wektor b.

Wzór Leibnitza i wzór Laplace'a

Wyznacznik macierzy można obliczyć w nastepujący sposób:

gdzie S_n jest zbiorem wszystkich permutacji zbioru . Koszt takiej operacji wynosi aż O(n n!) operacji i dlatego on się nie nadaje do obliczeń numerycznych.

Podobnie jest ze wzorem Laplace'a

gdzie A_1i jest macierzą o rozmiarach (n-1)x(n-1) powstałą z A poprzez wykreślenie pierwszej kolumny i i-tego wiersza. Koszt obliczenia wyznacznika jest rzędu O(2ⁿ).

Rozwiązywanie układu równań liniowych z macierzą górnotrójkątną

Układ

rozwiązuje się w sposób oczywisty poprzez

Koszt wyliczenia rozwiązania wynosi O(n²).

Eliminacja Gaussa

Poprzez eliminację Gaussa układ równań Ax = b przekształcimy w łatwy do rozwiązywania układ Rx = b', przy czym R jest macierzą górnotrójkątną. Znalezienie rozkładu A = L R, gdzie L z kolei jest macierzą dolnotrójkątną, nazywa się również faktoryzacja.

Załóżmy, że mamy już

Stosując przekształcenia

otrzymujemy macierz A^{(k + 1)}, przy czym i

Koszt wyliczenia rozwiązania układu eliminacją Gaussa szacuje sie na O(n³) operacji.

Eliminacja Gaussa dla macierzy trójprzekątniowej

Niech Ax = f, gdzie

Wtedy

gdzie

Faktoryzacja Choleskiego

Niech macierz A będzie symetryczna i dodatnio określona. Wtedy istnieje macierz dolnotrójkątna L taka, że A = LL^T, czyli

Z tego wynika następujący algorytm:

Algorytm: Faktoryzacja Choleskiego
1. Dla k = 1 do n: 2. 3. Dla i = k+1 do n: 4. 5. Koniec pętli 6. Koniec pętli

Koszt rozwiązywania układu równań metodą Choleskiego to połowa kosztu rozwiązywania eliminacją Gaussa.

Przykładowy program w Pascalu:

program równania_liniowe;

uses crt;
const max=20;
type   matrix=array[1..max+1,1..max+1]of extended;
       vector=array[1..max+1]of extended;
procedure gauss(n: integer; a: matrix;var x:vector);
{ Optymalna szybkosc dzialania }
{ Eliminacja Gaussa wierszami }
   var
      i, j, k, wm: integer;
      p, w, max  : extended;
   begin
      w:= 1;
      for i:= 1 to n-1 do
         begin
            max:= abs (a[i,i]);
            wm:= i;
            for j:= i+1 to n do         { wybor maxymalnego elementu }
               if abs (a[j,i]) > max then
                  begin
                     max:= abs (a[j,i]);
                     wm:= j
                  end;
            if wm > i then              { zamiana wierszy }
               begin
                  for k:= i to n+1 do
                     begin
                        p:= a[i,k];
                        a[i,k]:= a[wm,k];
                        a[wm,k]:= p
                     end;
                  w:= -w                { zamiana wierszy zmienia znak }
               end;
            if a[i,i] = 0 then          { macierz osobliwa }
               w:= 0;
            if w <> 0 then              { eliminacja Gaussa }
               for j:= i+1 to n do
                  begin
                     p:= a[j,i] / a[i,i];
                     for k:= n+1 downto i+1 do
                        a[j,k]:= a[j,k] - p * a[i,k]
                  end;
            w:= w * a[i,i]              { iloczyn elementow na przekatnej }
         end;
      w:= w * a[n,n];
x[n+1]:=w;
if w<>0 then
for i := n downto 1 do
    begin
      p := A [i,n+1];
      for j := i + 1 to n do
        p := p - A [i,j] * x [j];
      x [i] := p/A[i,i];
    end;
   end;
 var i,j,n:integer;
 a:matrix;
 x:vector;
 ch:char;
 begin
 clrscr;
 repeat
 writeln('Obliczanie ukladow rownan.');
 writeln;
 write('Podaj ilosc rownan n=');
 readln(n);
 for i:=1 to n do
 begin
 writeln('Wprowadz ',i,'-te rownanie');
 for j:=1  to n+1 do
 begin
 write('a[',i,',',j,']=');
 readln(a[i,j]);
 end;
 writeln;
 end;
 writeln;
 gauss(n,a,x);
 for i:=1 to n+1 do
 if i=n+1 then writeln('det(A)=',x[i]:10:6)
 else writeln('x[',i,']=',x[i]:10:6);
 ch:=readkey;
 until ch=#27;
 end.