1. Metody rozwiązywania układów równań liniowych

1.1 Wstęp teoretyczny

Układ równań może być zapisany w postaci:

lub w postaci sumy:

lub w postaci wektorowej:

AX=B, gdzie:

, ,

Rozwiązanie takiego układu równań polega na znalezieniu wartości zmiennych x_i (i = 1,2,...n) pod warunkiem, że znane są wartości współczynników a_ij i b_i (i = 1,2,...,n j = 1,2,...,n) oraz, że a_ij i b_i należą do zbioru liczb Rzeczywistych.

Metody rozwiązywania dzieli się na:

a) dokładne (tj. po wykonaniu skończonej liczby operacji, otrzymuje się rozwiązanie dokładne - przy zaniedbaniu błędów zaokrągleń), np. metoda Gaussa, Gaussa-Jordana,

b) iteracyjne (tj. wychodząc z pewnego przybliżonego rozwiązania, oblicza się kolejne, lepsze przybliżenie), np. metoda Jacobiego.

1.2 Rozkład LU

Rozkład LU macierzy, to inaczej przedstawienie jej w postaci iloczynu macierzy trójkątnej dolnej, mającej na diagonali głównej jedynki i trójkątnej górnej.

A = L * U

Schemat obliczania:

1. należy sprawdzić, czy macierz jest nieosobliwa (posiada rozwiązanie) oraz czy na głównej diagonali nie ma 0, jeśli są należy tak poprzestawiać wiersze macierzy, aby je wyeliminować,

Obliczanie macierzy L

1. w pracy wszystkie obliczenia są wykonywane dla macierzy kwadratowych, wiec można założyć, że macierz L ma taki sam rozmiar co macierz A, czyli n,

2. zakładamy, że elementy dla i = 1, 2, ..., n, j = 1, 2, ..., n, i = j,

3. dla i = 1, 2, ..., n, j = i+1, i+2, ..., n, wykonujemy następujące obliczenia:

a dla p = 1, 2, ..., n:

Obliczanie macierzy U

1. do obliczeń macierzy U korzystamy z oryginalnej macierzy A, a nie tej zniekształconej podczas obliczeń macierzy L,

2. obliczenia wykonujemy dla k wierszy dla k = 2, 3, ..., n wg poniższych wzorów:

3. otrzymana macierz jest macierzą U.

1.3 Metoda eliminacji Gaussa

Metoda ta polega na sprowadzeniu układu równań AX=B do A⁽ⁿ⁾X=B⁽ⁿ⁾, gdzie macierz A⁽ⁿ⁾ jest macierzą trójkątną górną a następnie za pomocą postępowania odwrotnego (zwanego również podstawieniem wstecznym) oblicza się wartości x_i.

Mając następujący układ równań:

należy go przekształcić do takiej postaci:

Macierz trójkątna górna, jest to macierz, która poniżej głównej diagonali ma same 0.

Aby macierz A przekształcić w macierz A⁽ⁿ⁾ należy wyeliminować zmienne spod głównej diagonali, stąd pochodzi nazwa metody – metoda eliminacji.

Schemat postępowania jest następujący:

1. należy sprawdzić, czy macierz jest nieosobliwa, czyli czy posiada rozwiązanie oraz, w przypadku, gdy na głównej diagonali występują 0, poprzestawiać wiersze tak, aby tych 0 tam nie było;

2. działania są wykonywane od drugiego wiersza układu (pierwszy element pierwszego wiersza a₁₁ leży na głównej diagonali, więc nie można go eliminować), czyli dla
   i = 2,3,...,n;

3. należy policzyć iloraz , dlatego trzeba było tak poprzestawiać wiersze, żeby nie było 0 na diagonali, dzielenie przez zero jest zabronione,

4. otrzymaną wartość odejmujemy od j (j – 1,2,..., n) elementu w wierszu i, czyli kolejne elementy w i-tym wierszu obliczamy za pomocą wzoru:

5. wykonując n-1 działań otrzymujemy układ równań, w którym ostatni wiersz ma postać:

jak widać, obliczenie x_n nie stanowi problemu, wystarczy podzielić przez ;

1. mając obliczoną wartość x_n można obliczyć x_n-1, czyli niewiadomą z wyższego wiersza:

2. obliczamy zmienne ‘od końca’, stad nazwa tego schematu – postępowanie odwrotne. Do obliczania wartości x_i służy wzór:

1.4 Metoda eliminacji zupełnej Gaussa-Jordana

Jest to metoda, która przekształca macierz A z układ równań AX=B do macierzy diagonalnej, czyli takiej, która na głównej diagonali ma 1, a w pozostałych miejscach 0. Docelowo układ równań:

ma mieć postać:

Jak widać końcowa postać układu równań jest jednocześnie jego rozwiązaniem, ponieważ wszystkie wiersze mają taka sama postać: , dla i = 1, 2, ..., n.

Schemat rozwiązania w części pokrywa się z metodą eliminacji Gaussa:

1. sprawdzamy, czy macierz jest nieosobliwa oraz czy nie ma zer na głównej diagonali, jeśli ma to je eliminujemy zmieniając kolejność wierszy,

2. poniższe instrukcje wykonujemy dla wszystkich wierszy, czyli dla k = 1, 2,..., n,

3. element musi mieć wartość 1, więc wszystkie elementy wiersza k dzielimy przez (o ile ),

, dla j = 1, 2, ..., n

1. teraz, gdy trzeba wyzerować współczynniki przy pozostałych , obliczeń dokonuje się wg wzoru:

,

dla i = 1, 2, ..., n, , j = 1, 2, ..., n

1.5 Metoda iteracji prostej Jacobiego

Metoda ta należy do grupy metod iteracyjnych, czyli obliczających względnie dokładny wynik. Względnie, ponieważ podaje się dokładność z jaka ma być obliczony wynik. Posługiwanie się tą metodą wymaga od liczącego dobrej znajomości działań na macierzach, zwłaszcza mnożenia macierzy.

Schemat obliczeń:

1. należy sprawdzić, czy macierz jest nieosobliwa oraz, czy na głównej diagonali nie ma 0, jeśli występują, dokonuje się przestawienia wierszy w układzie równań tak, aby je wyeliminować,

2. zakładamy dokładność z jaką chcemy obliczyć wynik,

3. następnie macierz główną układu – A rozkładamy na sumę macierzy L + D + U, gdzie L jest macierzą poddiagonalną, D – macierzą diagonalną a U – macierzą naddiagonalną:

= + +

1. do dalszych obliczeń potrzebna jest macierz , czyli należy odwrócić macierz D, jednym ze sposobów jest przekształcenie macierzy:

w macierz:

w taki sam sposób, jak w metodzie Gaussa-Jordana przekształca się macierz A w macierz diagonalną. Otrzymana macierz:

jest szukana macierzą ,

1. teraz należy macierz D pomnożyć przez , jeśli jako wynik powstanie macierz diagonalna, to macierz została odwrócona prawidłowo:

* =

1. teraz na podstawie wzoru:

obliczamy wektor X. Obliczeń dokonujemy iteracyjnie, warunkiem stopu jest spełnienie warunku:

założona dokładność obliczeń

lub warunkiem stopu jest założona maksymalna ilość iteracji.

Przykład programu w Delphi

{$APPTYPE CONSOLE}

uses
  SysUtils;

const
  nmax=10;

var
  x,b:array[1..nmax] of extended;
  A:array[1..nmax,1..nmax] of extended;
  n,z,w,k,gw:integer;
  m,d,s:extended;           // n - rozmiar macierzy
                            // w - liczba wierszy
                            // k - liczba kolumn
                            // gw - liczba kolejejnych krokow wyzerowan kolumn
                            // m - mnoznik
                            // d - mnozenie wspol*x
                            // s -suma wspol.*x

begin
  { TODO -oUser -cConsole Main : Insert code here }

  writeln('Program do rozwiazywania rownania macierzowego A*x=b metoda Gaussa !!!');

  // 1.Podac wymiar macierzy do n = 10 -> (n)
  write('Podaj wymiar macierzy A: ');
  readln(n);
  writeln;

  // 2.Podac dokladnosc zaokraglnia
  writeln('Podaj dokladnosc zaokraglania: ');
  readln(z);
  writeln;

  // 3.Podac elementy macierzy A -> akw k=1..n , w=1..n
   writeln('Podaj elementy macierzy A: ');
   for w:=1 to n do
     begin
       for k:=1 to n do
         begin
            write('A[',w,',',k,']: ');
            readln(A[w,k]);
         end;
     end;
     writeln;

  // 3.Podac elementy wyrazow wolnych b -> bw w=1..n
  writeln('Podaj elementy macierzy b: ');
  for w:=1 to n do
    begin
      write('b[',w,']: ');
      readln(b[w]);
    end;
    writeln;

  // 5.Doprowadzic maciez [A|b] do postaci gornotrojkatnej metoda eliminacji Gausa (postepowanie proste)
   for gw:=1 to n-1 do
     begin
       for w:=1+gw to n do
         begin
            m:=A[w,gw]/A[gw,gw];
            for k:=gw to n do
                begin
                  A[w,k]:=A[w,k]-m*a[gw,k];
                end;
            b[w]:=b[w]-m*b[gw];
         end;
     end;

{  // 6. Wyswietlenie macierzy gornotrojkatnej
  writeln('Macierz gornotrojkatna [A|b] (eliminacja Gaussa) wyglada tak: ');
  for w:=1 to n do
   begin
     for k:=1 to n do
       begin
       writeln('A[',w,',',k,'] = ',a[w,k]:0:z);
       end;
     writeln('b[',w,'] = ',b[w]:0:z);
     writeln;
   end;  }

  // 7. Wyznaczyc rozwiazania x -> xw w=n..1
   writeln('Rozwiazaniami sa liczby: ');
   for w:=n downto 1 do
     begin
       s:=0;
       for k:=w+1 to n do
         begin
           d:=a[w,k]*x[k];
           s:=s+d;
         end;
         x[w]:=(b[w]-s)/a[w,w];
         writeln('x[',w,'] = ',x[w]:0:z);
     end;
  readln;
end.