zatem mamy
co potwierdza tym samym dużą wrażliwość układów równań z taką macierzą współ-czynników na zaburzenie danych.
Zbadanie, czy układ równań jest źle uwarunkowany może być bardzo trudne. Jeśli istnieje przypuszczenie, że układ równań jest źle uwarunkowany to można zastosować np. jego skalowanie [6], polegające na znalezieniu takiej macierzy diagonalnej D lub też dwóch macierzy diagonalnych
oraz
żeby wskaźniki uwarunkowania albo
były możliwie małe.
2.3. Metoda eliminacji Gaussa
Układ równań (2.35) zapiszemy w postaci:
(2.49)
w której macierz układu A została oznaczona jako a wyrazy wolne zostały umieszczone w dodatkowej (n +1) - szej kolumnie macierzy
Metoda eliminacji Gaussa polega na przekształceniu rozwiązywanego układu równań (2.49) do równoważnego układu o macierzy trójkątnej górnej. Zakładając, że i odejmując pierwsze równanie pomnożone przez od pozostałych równań (i = 2, 3, ..., n) otrzymujemy pierwszy układ zredukowany:
(2.50)
gdzie
Następnie, jeśli
, to odejmujemy drugie równanie układu (2.50) pomnożone przez od dalszych równań tego układu. Uzyskujemy drugi układ zredukowany:
(2.51)
gdzie
Kontynuując ten sposób obliczeń po wykonaniu kroków eliminacyjnych otrzymujemy układ równań o górnej macierzy trójkątnej
(2.52)
(2.52cd.)
którego współczynniki dla każdego
są określone zależnościami:
(2.53)
Rozwiązanie układu (2.52) można już łatwo wyznaczyć ze wzorów:
(2.54)
Etap pierwszy obliczeń prowadzący do układu (2.52) - (2.53) nazywa się eliminacją zmiennych, etap drugi (2.54) - postępowaniem odwrotnym.
Metoda eliminacji Gaussa w swej najprostszej postaci jest niezawodna, jeśli wszystkie współczynniki
(k = 1, 2, ..., n) są różne od zera oraz jeśli moduły elementów macierzy i moduły elementów kolejnych macierzy zredukowanych są tego samego rzędu. W przypadku bowiem gdy moduły współczynników układu zredukowanego są dużo większe od modułów współczynników układu przed redukcją może nastąpić znaczna utrata dokładności obliczeń.
Dla zapewnienia pełnej niezawodności metody eliminacji Gaussa stosowane są dwie jej modyfikacje:
1) eliminacja z częściowym wyborem elementu podstawowego,
2) eliminacja z pełnym wyborem elementu podstawowego.
Częściowy wybór elementu podstawowego polega na wyszukiwaniu przed każdym krokiem eliminacji zmiennych największego co do modułu elementu
w ciągu elementów Po ustaleniu numeru r równania, w którym ten element występuje, wiersz o numerze r zamieniany jest miejscem z wierszem o numerze k. Gdy istnieje jednoznaczne rozwiązanie układu równań żaden element podstawowy nie będzie równy zeru, gdyż oznaczałoby to liniową zależność wierszy, z których wybierany był ten element. Czynności przesta-wiania wierszy musi towarzyszyć zmiana znaków w dowolnym spośród przestawianych równań, jeśli wyznacznik macierzy współczynników układu powstałego w wyniku przestawiania wierszy ma być równy wyznacznikowi macierzy współczynników układu wyjściowego.
Przy pełnym wyborze elementu podstawowego wyszukiwany jest największy co do modułu element
w całej macierzy
Po ustaleniu numeru r wiersza i numeru s kolumny następuje przestawienie wiersza o numerze r z wierszem o numerze k, a następnie kolumny o numerze s z kolumną o numerze k. Liczba przestawień wierszy i kolumn jest zapamiętywana i uwzględniana przy obliczaniu wyznacznika macierzy współczynników.
W celu wyznaczenia rozwiązania układu równań liniowych metodą eliminacji Gaussa należy wykonać [9]
mnożeń i dzieleń oraz
dodawań. Zatem w przypadku układu 15 równań: M = 1345, D = 1235, a więc komputer wykonujący milion operacji na sekundę rozwiąże ten układ w czasie ok. 0.01 sekundy.
Metoda eliminacji Gaussa jest metodą bardzo efektywną i zapewniającą najwyższą jakość ze względu na błędy obliczeń przy wykonywaniu obliczeń z wyborem elementu podstawowego.
Obliczony wektor rozwiązań dla układu równań liniowych może być poprawiony przy wykorzystaniu następującego procesu iteracyjnego:
- obliczamy wektor reszt:
- rozwiązujemy układ równań:
- wyznaczamy nowe przybliżenie:
dla i = 0, 1, .... Obliczenia powtarzamy, jeśli maleje miara błędu:
Metodę eliminacji Gaussa można stosować do równoczesnego rozwiązywania wielu układów równań, różniących się jedynie wektorem prawych stron. W szczególności można wyznaczyć macierz odwrotną.
Oznaczając elementy nieosobliwej macierzy kwadratowej A przez , a elementy macierzy odwrotnej przez - z definicji macierzy odwrotnej (2.29) otrzymujemy następujące równanie macierzowe
(2.55)
Mnożąc macierze A i uzyskujemy n układów równań względem
niewiadomych
(2.56)
przy czym wszystkie elementy wektorów prawych stron są zerami - oprócz jednego, który jest jedynką przypadającą w wierszu o numerze j.
Zastosowanie metody eliminacji Gaussa do rozwiązywania układu równań liniowych (2.36) pozwala niejako „przy okazji” na obliczenie wyznacznika macierzy A tego układu. Łatwo bowiem można sprawdzić, że wyznacznik macierzy trójkątnej
górnej w przekształconym układzie równań (2.52) jest równy iloczynowi jej elementów diagonalnych. Uzyskujemy więc zależność
(2.57)
przy czym przy obliczaniu wyznacznika może wystąpić nadmiar bądź niedomiar zmiennopozycyjny, jeśli rząd elementów diagonalnych jest bardzo duży lub bardzo mały.
*
Program 2.1 jest przeznaczony do rozwiązywania układu równań liniowych (2.36) metodą eliminacji Gaussa z częściowym wyborem elementu podstawowego. W programie tym wykonywane są następujące czynności:
1) wczytywanie danych zapisanych w pliku Pr_2_1.dan mających następujące znaczenie:
n, m - liczby równań i prawych stron,
A[1..n,1..n+m] - tablica zawierająca zarówno współczynniki układu równań jak i wektory prawych stron,
2) rozwiązywanie układu równań za pomocą procedury
ElimGaussa(n,m,A,det), (2.58)
w której obliczane wektory rozwiązań są wpisywane na miejsce odpowiadających im wektorów prawych stron, a wartość wyznacznika podstawiana jest pod zmienną rzeczywistą det,
3) obliczanie macierzy odwrotnej z układów równań (2.55) - (2.56),
4) drukowanie danych i wyników obliczeń.
Tabulogram modułu Obliczenia programu 2.1 jest następujący:
{Program 2.1}
unit Obliczenia;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls,
Forms, Dialogs, StdCtrls, Buttons;
const nmax = 20;
type
Tabl = array[1..nmax,1..2*nmax] of Real;
. . . . . . . . . . . . . . . . . . . . .
var
Form3: TForm3;
i,j,k,n,m: Integer;
nazwa: String;
plik: Text;
det: Real;
A,B: Tabl;
implementation
uses Ustawienia, Informacje, Grafika, Podglad;
{$R *.DFM}
procedure ElimGaussa(n,m: Integer; var A: Tabl; var det: Real);
var
i,j,k,r: Integer;
max,zp: Real;
label omin;
begin
for k:=1 to n-1 do begin
max:=0; r:=k;
for i:=k to n do
if Abs(A[i,k])>max then begin
max:=Abs(A[i,k]);; r:=i;
end;
if max<1e-10 then goto omin;
for j:=k to n+m do begin
zp:=A[k,j];
A[k,j]:=-A[r,j];
A[r,j]:=zp;
end;
for i:=k+1 to n do begin
zp:=A[i,k]/A[k,k];
for j:=k to n+m do
A[i,j]:=A[i,j]-zp*A[k,j];
end;
end;
det:=1;
for i:=1 to n do
det:=det*A[i,i];
if (Abs(det)<1e-10) or (m<1) then goto omin;
for k:=n+1 to n+m do begin
A[n,k]:=A[n,k]/A[n,n];
for i:=n-1 downto 1 do begin
zp:=0;
for j:=i+1 to n do
zp:=zp+A[i,j]*A[j,k];
A[i,k]:=(A[i,k]-zp)/A[i,i];
end;
end;
omin:
end;
. . . . . . . . . . . . . . . . . . . . . . .
procedure TForm3.BitBtn1Click(Sender: TObject);
begin
Form2.Show;
AssignFile(plik,nazwa);
Reset(plik);
Readln(plik,n,m);
for i:=1 to n do
for j:=1 to n+m do begin
Read(plik,A[i,j]);
B[i,j]:=A[i,j];
end;
CloseFile(plik);
AssignFile(plik,Edit1.Text);
Rewrite(plik); Writeln(plik,'PROGRAM 2.1.');
Writeln(plik,'Rozwiązywanie układu równań liniowych.');
Writeln(plik,'Metoda eliminacji Gaussa z częściowym wyborem');
Writeln(plik,'elementu podstawowego.');
Writeln(plik,'');
Writeln(plik,'Liczba równań układu - n =',n:3);
Writeln(plik,'Liczba prawych stron - m =',m:3);
Writeln(plik,'');
Writeln(plik,'Macierz współczynników:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:2);
k:=0; Write(plik,' ');
for j:=1 to n do begin
k:=k+1;
if k=5 then begin
k:=0; Writeln(plik);
Write(plik,' ');
end;
Write(plik,' ',A[i,j]:16);
end;
Writeln(plik);
end;
Writeln(plik);
Writeln(plik,'Wektory prawych stron:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
Write(plik,' ');
for j:=n+1 to n+m do
Write(plik,' ',A[i,j]:16);
Writeln(plik);
end;
Writeln(plik);
ElimGaussa(n,m,A,det);
Writeln(plik,'Wyznacznik - det = ',det:16);
Writeln(plik);
Writeln(plik,'Rozwiązania układów równań:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
Write(plik,' ');
for j:=n+1 to n+m do
Write(plik,' ',A[i,j]:16);
Writeln(plik);
end;
Writeln(plik);
for i:=1 to n do
for j:=1 to n do
if i=j then B[i,n+j]:=1
else B[i,n+j]:=0;
ElimGaussa(n,n,B,det);
Writeln(plik,'Macierz odwrotna:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
k:=0; Write(plik,' ');
for j:=1 to n do begin
k:=k+1;
if k=5 then begin
k:=0; Writeln(plik);
Write(plik,' ');
end;
Write(plik,' ',B[i,n+j]:16);
end;
Writeln(plik);
end;
Writeln(plik); CloseFile(plik);
Form2.Wyniki.Lines.LoadFromFile(Edit1.Text);
end;
. . . . . . . . . . . . . . . . . . . . . . . . .
procedure TForm3.SpeedButton1Click(Sender: TObject);
begin
if SaveDialog1.Execute then begin
Dane.Lines.SaveToFile(SaveDialog1.FileName);
nazwa:=SaveDialog1.FileName;
end;
end;
end.
Przy wykorzystaniu programu 2.1 rozwiązano układ równań:
(2.59)
i wyznaczono macierz odwrotną do macierzy współczynników. Rozwiązaniem układu (2.59) są pierwiastki:
macierz odwrotna do macierzy współczynników jest następująca
Po wczytaniu odpowiednich danych z pliku Pr_2_1.dan otrzymujemy plik Pr_2_1.wyn zawierający rozwiązanie układu równań (2.59) oraz elementy macierzy odwrotnej który ma następującą postać:
68 2. Układy równań liniowych
2.3. Metoda eliminacji Gaussa 65