![]() |
Wyjście Spis treści Poprzedni Następny
Autor artykułu: mgr Jerzy Wałaszek, wersja 2.0 |
©2014 mgr Jerzy Wałaszek |
Znaleźć rozwiązanie układu równań liniowych postaci:
a1,1x1 + a1,2x2
+ ... + a1,nxn = b1
a2,1x1 + a2,2x2
+ ... + a2,nxn = b2
...
an,1x1 + an,2x2
+ ... + an,nxn = bn
Wykorzystując metodę eliminacji Gaussa-Crouta.
Opisany w poprzednim rozdziale algorytm eliminacji Gaussa posiada dosyć istotną wadę. Bazuje on na dzieleniu przez elementy przekątnej głównej macierzy współczynników. Otóż może się zdarzyć, iż element przekątnej głównej jest równy zero lub zostanie wyzerowany w wyniku obliczeń. W takim przypadku algorytm zostanie przerwany z błędem dzielenia przez zero. Jako przykład niech posłuży poniższy układ równań:
0x1 + 2x2 + 3x3 + 4x4
= 49
1x1 + 0x2 + 3x3 + 4x4 = 45
1x1 + 2x2 + 0x3 + 4x4 = 36
1x1 + 2x2 + 3x3 + 0x4 = 23
Dane dla programów dla tego układu są następujące:
4
0 2 3 4 49
1 0 3 4 45
1 2 0 4 36
1 2 3 0 23
Po wprowadzeniu tych danych do przykładowego programu umieszczonego w poprzednim rozdziale otrzymujemy komunikat "DZIELNIK ZERO". Tymczasem rozwiązanie równania istnieje i jest następujące:
x1 = 2
x2 = 3
x3 = 5
x4 = 7
Ulepszeniem algorytmu eliminacji Gaussa jest metoda Crouta, która polega na tym, iż na początku eliminacji wyszukujemy w wierszu macierzy AB element o największym module, po czym zamieniamy miejscami kolumnę ze znalezionym elementem z kolumną zawierającą element głównej przekątnej. W ten sposób dzielnik będzie posiadał największą na moduł wartość i pozbędziemy się sytuacji, gdy może on posiadać wartość 0 (chyba że cały wiersz macierzy AB jest zerowy – ale wtedy równania i tak się nie da rozwiązać).
Tak zmodyfikowany algorytm eliminacji Gaussa wymaga zapamiętania, które kolumny zostały zamienione, ponieważ na etapie wyznaczania niewiadomych musimy znać numery wyliczanych niewiadomych. Osiągniemy to wprowadzając do algorytmu dodatkowy wektor W przechowujący informację o numerach kolumn – zamiana miejscami kolumn nie wymaga faktycznego ich przemieszczania w pamięci – wystarczy odwoływać się do nich poprzez wektor kolumn. Jeśli w wektorze W zamienimy ze sobą numery dwóch kolumn (są to tylko dwa elementy, zatem wymiana jest w czasie stałym), to efekt będzie taki, jakby te kolumny zostały zamienione miejscami w macierzy AB:
Na powyższym rysunku widzimy odwołanie do kolumny 2 macierzy poprzez wektor numerów kolumn, oznaczony tutaj kolorem żółtym. Gdy kolumny nie były zamieniane, to odwołanie prowadzi do tej samej kolumny – na pozycji 2 w wektorze kolumn jest numer 2. Jednakże po zamianie (tutaj kolumny nr 2 z nr 5) odwołanie do kolumny nr 2 przenosi nas do kolumny nr 5, chociaż w macierzy docelowej kolumny te znajdują się dalej na swoich miejscach. Rozwiązanie to zaoszczędzi nam czasu przy przenoszeniu elementów. Dodatkowo zawartość wektora numerów kolumn może być wykorzystana w końcowej fazie algorytmu przy wyznaczaniu wartości zmiennych.
n | – | liczba niewiadomych |
AB | – | macierz n × (n + 1) zawierająca współczynniki ai,j, i = 1...n; j = 1...n oraz w kolumnie n + 1 współczynniki bi układu równań |
X | – | wektor n elementowy, w którym zostaną umieszczone niewiadome xi. |
Jeśli zmienna r ma wartość true, to wektor X zawiera rozwiązanie układu równań. W przeciwnym razie nastąpiło dzielenie przez zero i algorytm nie może znaleźć poprawnych rozwiązań układu równań. Pierwotna macierz współczynników AB zostaje zniszczona – zastąpiona macierzą trójkątną, powstałą w czasie redukcji współczynników.
W | – | wektor n + 1 elementowy zawierający numery kolumn macierzy AB.
W
![]() |
i,j,k | – | zmienne sterujące pętli. i,j,k
![]() |
m | – | mnożnik, przez który będą mnożone elementy macierzy AB.
m
![]() |
s | – | zlicza sumę iloczynów. s
![]() |
ε | – | określa dokładność porównania z zerem; ε = 10-12, ε
![]() |
K01: | r ← false | ; zakładamy porażkę | |||
K02: | ε ← 10-12 | ; określamy dokładność porównania z zerem | |||
K03: | Dla i = 1,2,...,n+1: wykonuj W[i] ← i | ; do wektora wpisujemy kolejne numery kolumn | |||
K04: | Dla i = 1,2,...,n-1 wykonuj K05...K11 | ; dokonujemy eliminacji współczynników | |||
K05: | k ← i | ||||
K06: |
Dla j = i+1, i+2, ...,n
wykonuj: jeśli | AB[i,W[k]] | < | AB[i,W[j]] |, to k ← j |
; wyszukujemy element o największym module | |||
K07: | W[k] ↔ W[i] | ; w wektorze zamieniamy numery kolumn wg ; znalezionego elementu |
|||
K08: | Jeśli | AB[i,W[i]] | < ε, to zakończ | ; jeśli znaleziony element jest zerem, kończymy | |||
K09: | Dla j = i+1, i+2, ..., n wykonuj K10...K11 | ||||
K10: |
|
; wyznaczamy mnożnik | |||
K11: | Dla
k = i+1, i+2, ..., n+1 wykonuj: AB[j,W[k]] ← AB[j,W[k]] + m × AB[i,W[k]] |
; sumujemy wiersz j-ty z wierszem i-tym ; przemnożonym przez m |
|||
K12: | Dla i = n, n-1, ..., 1 wykonuj K13...K16 | ; wyliczamy kolejne niewiadome | |||
K13: | Jeśli | AB[i,W[i]] | < ε, to zakończ | ||||
K14: | s ← AB[i,n+1] | ; do s trafia współczynnik b | |||
K15: | Dla j = n, n-1, ..., i+1 wykonuj: s ← s - AB[i,W[j]] × X[W[j]] | ; zliczamy sumę iloczynów współczynników przez niewiadome | |||
K16: |
|
; obliczamy kolejną niewiadomą | |||
K17: | r ← true | ; zaznaczamy sukces | |||
K18: | Zakończ |
Ważne: Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich. |
Program wczytuje dane definiujące macierz współczynników AB i wylicza niewiadome X. Jeśli obliczenie nie może zostać wykonane, program wypisuje komunikat "DZIELNIK ZERO".
Dane dla macierzy AB są zdefiniowane następująco:
Pierwsza liczba określa liczbę niewiadomych n.
Następne n × (n + 1) liczb określa zawartość macierzy AB wierszami – pierwsze n liczb każdego wiersza odnosi się do współczynników a, a (n + 1)-sza liczba odnosi się do wyrazu wolnego b. Poniżej mamy przykładowe dane wejściowe dla programu:
Równanie:
2x2 + 3x3 + 4x4
= 49
x1 + 3x3 + 4x4 = 45
x1 + 2x2 + 4x4 = 36
x1 + 2x2 + 3x3 = 23
Dane dla programu:
4
0 2 3 4 49
1 0 3 4 45
1 2 0 4 36
1 2 3 0 23
Lazarus |
// Eliminacja Gaussa-Crouta // Data: 3.03.2010 // (C)2012 mgr Jerzy Wałaszek //----------------------------- program gausscrout; type NInt = array of integer; // typ wektora kolumn NReal = array of double; // typ tablic wierszy MReal = array of NReal; // typ tablicy wskaźników const eps = 1e-12; // Funkcja realizuje algorytm eliminacji Gaussa-Crouta //---------------------------------------------------- function Gauss(n : integer; AB : MReal; X : NReal; W : NInt) : boolean; var i,j,k : integer; m,s : double; begin for i := 0 to n do W[i] := i; // eliminacja współczynników for i := 0 to n - 2 do begin k := i; for j := i + 1 to n - 1 do if abs(AB[i][W[k]]) < abs(AB[i][W[j]]) then k := j; j := W[k]; W[k] := W[i]; W[i] := j; for j := i + 1 to n - 1 do begin if abs(AB[i][W[i]]) < eps then exit(false); m := -AB[j][W[i]]/ AB[i][W[i]]; for k := i + 1 to n do AB[j][W[k]] := AB[j][W[k]] + m * AB[i][W[k]]; end; end; // wyliczanie niewiadomych for i := n - 1 downto 0 do begin if abs(AB[i][W[i]]) < eps then exit(false); s := AB[i][n]; for j := n - 1 downto i + 1 do s := s - AB[i][W[j]] * X[W[j]]; X[W[i]] := s / AB[i][W[i]]; end; Gauss := true; end; // Program główny var AB : MReal; X : NReal; W : NInt; n,i,j : integer; begin // odczytujemy liczbę niewiadomych read(n); // tworzymy macierze AB, X i W o odpowiednich rozmiarach SetLength(AB,n); SetLength(X,n); SetLength(W,n+1); for i := 0 to n - 1 do SetLength(AB[i],n + 1); // odczytujemy dane dla macierzy AB for i := 0 to n - 1 do for j := 0 to n do read(AB[i][j]); if Gauss(n,AB,X,W) then begin for i := 0 to n - 1 do writeln('x',i+1,' = ',X[i]:9:4); end else writeln('DZIELNIK ZERO'); // usuwamy macierze z pamięci for i := 0 to n - 1 do SetLength(AB[i],0); SetLength(AB,0); SetLength(X,0); SetLength(W,0); end. |
Code::Blocks |
// Eliminacja Gaussa-Crouta // Data: 3.03.2010 // (C)2012 mgr Jerzy Wałaszek //----------------------------- #include <iostream> #include <iomanip> #include <cmath> using namespace std; const double eps = 1e-12; // Funkcja realizuje algorytm eliminacji Gaussa-Crouta //---------------------------------------------------- bool gauss(int n, double ** AB, double * X, int * W) { int i,j,k; double m,s; for(i = 0; i <= n; i++) W[i] = i; // eliminacja współczynników for(i = 0; i < n - 1; i++) { k = i; for(j = i + 1; j < n; j++) if(fabs(AB[i][W[k]]) < fabs(AB[i][W[j]])) k = j; swap(W[k],W[i]); for(j = i + 1; j < n; j++) { if(fabs(AB[i][W[i]]) < eps) return false; m = -AB[j][W[i]] / AB[i][W[i]]; for(k = i + 1; k <= n; k++) AB[j][W[k]] += m * AB[i][W[k]]; } } // wyliczanie niewiadomych for(i = n - 1; i >= 0; i--) { if(fabs(AB[i][W[i]]) < eps) return false; s = AB[i][n]; for(j = n - 1; j >= i + 1; j--) s -= AB[i][W[j]] * X[W[j]]; X[W[i]] = s / AB[i][W[i]]; } return true; } // Program główny int main() { double **AB, *X; int n,i,j,*W; cout << setprecision(4) << fixed; // odczytujemy liczbę niewiadomych cin >> n; // tworzymy macierze AB, X i W AB = new double * [n]; X = new double [n]; W = new int [n + 1]; for(i = 0; i < n; i++) AB[i] = new double[n + 1]; // odczytujemy dane dla macierzy AB for(i = 0; i < n; i++) for(j = 0; j <= n; j++) cin >> AB[i][j]; if(gauss(n,AB,X,W)) { for(i = 0; i < n; i++) cout << "x" << i + 1 << " = " << setw(9) << X[i] << endl; } else cout << "DZIELNIK ZERO\n"; // usuwamy macierze z pamięci for(i = 0; i < n; i++) delete [] AB[i]; delete [] AB; delete [] X; delete [] W; return 0; } |
Free Basic |
' Eliminacja Gaussa-Crouta ' Data: 3.03.2010 ' (C)2012 mgr Jerzy Wałaszek '----------------------------- Const eps As Double = 1e-12 ' Funkcja realizuje algorytm eliminacji Gaussa-Crouta '---------------------------------------------------- Function Gauss(n As Integer, AB As Double Ptr Ptr, X As Double Ptr, W As Integer Ptr) As Integer Dim As Integer i,j,k Dim As Double m,s For i = 0 To n W[i] = i Next ' eliminacja współczynników For i = 0 To n - 2 k = i For j = i + 1 To n - 1 If Abs(AB[i][W[k]]) < Abs(AB[i][W[j]]) Then k = j Next Swap W[k], W[i] For j = i + 1 To n - 1 If Abs(AB[i][W[i]]) < eps Then Return 0 m = -AB[j][W[i]]/ AB[i][W[i]] For k = i + 1 To n AB[j][W[k]] += m * AB[i][W[k]] Next Next Next ' wyliczanie niewiadomych For i = n - 1 To 0 Step -1 If Abs(AB[i][W[i]]) < eps Then Return 0 s = AB[i][n] For j = n - 1 To i + 1 Step -1 s -= AB[i][W[j]] * X[W[j]] Next X[W[i]] = s / AB[i][W[i]] Next Gauss = 1 End Function ' Program główny '--------------- Dim AB As Double Ptr Ptr Dim X As Double Ptr Dim W As Integer Ptr Dim As Integer n,i,j Open Cons For Input As #1 ' odczytujemy liczbę niewiadomych Input #1,n ' tworzymy macierze AB, X i W o odpowiednich rozmiarach AB = New Double Ptr [n] X = New Double [n] W = New Integer [n + 1] For i = 0 To n - 1 AB[i] = New Double [n + 1] Next ' odczytujemy dane dla macierzy AB For i = 0 To n - 1 For j = 0 To n Input #1,AB[i][j] Next Next Close #1 If Gauss(n,AB,X,W) = 1 Then For i = 0 To n - 1 Print Using "x# = ####.####"; i + 1; X[i] Next Else Print "DZIELNIK ZERO" End If ' usuwamy macierze z pamięci For i = 0 To n - 1 Delete [] AB[i] Next Delete [] AB Delete [] X Delete [] W End |
Wynik |
4 0 2 3 4 49 1 0 3 4 45 1 2 0 4 36 1 2 3 0 23 x1 = 2.0000 x2 = 3.0000 x3 = 5.0000 x4 = 7.0000 |
W związku z dużą liczbą listów do naszego serwisu edukacyjnego nie będziemy udzielać odpowiedzi na prośby rozwiązywania zadań, pisania programów zaliczeniowych, przesyłania materiałów czy też tłumaczenia zagadnień szeroko opisywanych w podręcznikach.
![]() | I Liceum Ogólnokształcące |