Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 1
Rozwiązywanie układów równań
Metody dokładne
Metody, w których rozwiązanie otrzymuje się po skończonej liczbie kroków. Gdyby wszystkie działania
arytmetyczne były przeprowadzane dokładnie otrzymane rozwiązanie również byłoby dokładne  oczywiście z
powodu błędów obcięcia dokładne nie będzie (należy jednak pamiętać, że w przypadku większości zastosowań błędy
rozwiązań wynikające z błędów obcięcia są niezauważalnie małe).
Przykłady i wzory będą zapisywane dla układu czterech równań z czterema niewiadomymi, ale oczywiście równań i
niewiadomych może być n.
Warunki istnienia rozwiązań
Aby rozwiązanie układu równań było możliwe macierz współczynników musi być nieosobliwa, ale w przypadku
obliczeń numerycznych jeśli któreś z równań będzie bliskie liniowej kombinacji pozostałych to już może to
uniemożliwić rozwiązanie. Poza tym w trakcie obliczeń kumulować się będą błędy obcięcia (zwłaszcza gdy macierz
jest bliska osobliwej) i mogą bardzo zniekształcić wyniki.
Sposób w jaki błędy (czy niedokładności) występujące w wektorze prawych stron wpływają na wynik można
określić następująco. Zakładamy, że dane jest równanie postaci:
Ax b
którego rozwiązaniem będzie oczywiście:
x A 1b
JeÅ›li bÅ‚Ä™dy wystÄ™pujÄ…ce w wektorze prawych stron oznaczymy jako ´b bÄ™dziemy mogli zapisać:
Ax b ; A´x ´b ; x A 1b ; ´x A 1´b
Wynika stąd (zakładając, że normę liczymy jako pierwiastek z sumy kwardatów):
A x b
A ´x ´b
x A 1 b
´x A 1 ´b
Można więc zapisać następujące nierówności:
1 ´b ´x 1 ´x
A ; A 1 ´b
A
A 1 b x x b
Wprowadzimy teraz oznaczenie:
K A A 1
które umożliwi nam połączenie zapisanych powyżej nierówności:
1 ´b ´x ´b
K
K
b x b
Nierówność ta okreÅ›la  granice bÅ‚Ä™du (wynikajÄ…cego z bÅ‚Ä™dów ´b) jakim obarczone może być rozwiÄ…zanie x. Widać,
że  im większe K tym gorzej . Macierze, dla których K przybiera duże wartości nazywane są  z
le uwarunkowanymi ,
w szczególności dla macierzy osobliwej K ma wartość nieskończoną.
Warto tu jeszcze raz powtórzyć, że o ile rozwiązanie  na karteczce (analityczne) nie powiedzie się (tzn. nie uda się
przedstawić wzorów określających rozwiązanie, kwestię obliczenia wartości liczbowych tego rozwiązania pozostawiam
 na boku ) tylko gdy macierz jest osobliwa, o tyle w przypadku rozwiązań numerycznych do braku rozwiązania
wystarczy, że macierz będzie  bliska osobliwej .
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 2
Wzory Cramera
Oczywiście jeszcze ze szkoły wszyscy pamiętają wzory Cramera (metodę wyznaczników). Jeśli dany jest układ
równań postaci:
a11x1 a12x2 a13x3 a14x4 b1
a21x1 a22x2 a23x3 a24x4 b2
a31x1 a32x2 a33x3 a34x4 b3
a41x1 a42x2 a43x3 a44x4 b4
to należy obliczyć jego wyznacznik główny ":
a11 a12 a13 a14
a21 a22 a23 a24
"
a31 a32 a33 a34
a41 a42 a43 a44
oraz wszystkie pozostałe wyznaczniki postaci (np.):
a11 b1 a13 a14 a11 a12 b1 a14
a21 b2 a23 a24 a21 a22 b2 a24
"2 ; "3
a31 b3 a33 a34 a31 a32 b3 a34
a41 b4 a43 a44 a41 a42 b4 a44
Po czym wyznaczyć niewiadome z wzorów:
"1 "2 "3 "4
x1 ; x2 ; x3 ; x4
" " " "
Jak łatwo zauważyć metoda ta będzie bardzo pracochłonna w przypadku dużych układów równań (dla układu n
równań wymaga obliczenia n+1 wyznaczników stopnia n).
Metoda eliminacji Gaussa
Załóżmy, że mamy układ równań (w przykładowych zapisach poniżej wektor wolnych wyrazów będę oznaczać
przez a z indeksami odpowiadającymi n+1 kolumnie (ma to podkreślać fakt, że podczas rozwiązywania układu równań
będziemy operować na całej tablicy, włącznie z wektorem wyrazów wolnych):
a11x1 a12x2 a13x3 a14x4 a15
a21x1 a22x2 a23x3 a24x4 a25
a31x1 a32x2 a33x3 a34x4 a35
a41x1 a42x2 a43x3 a44x4 a45
załóżmy, że element a11 0, jeśli podzielimy pierwsze równanie układu przez a11 otrzymamy:
x1 b12x2 b13x3 b14x4 b15
gdzie:
a1j
b1j j 2,3,4,5
a11
Wykorzystując powyższe równanie można łatwo wyeliminować zmienną x1 z pozostałych równań. W tym celu
mnożymy je odpowiednio przez a21, a31 i a41 i odejmujemy od odpowiednich równań. Otrzymujemy układ trzech
równań:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 3
(1) (1) (1) (1)
a22 x2 a23 x3 a24 x4 a25
(1) (1) (1) (1)
a32 x2 a33 x3 a34 x4 a35
(1) (1) (1) (1)
a42 x2 a43 x3 a44 x4 a45
gdzie współczynniki:
(1)
aij aij ai1b1i i 2,3,4 j 2,3,4,5
Znowu pierwsze równanie układu dzielimy przez jego pierwszy współczynnik (a22(1)) i otrzymujemy:
(1) (1) (1)
x2 b23 x3 b24 x4 b25
gdzie:
(1)
a2j
(1)
b2j j 3,4,5
(1)
a22
W taki sam jak poprzednio sposób eliminujemy zmienną x2 i otrzymujemy układ:
(2) (2) (2)
a33 x3 a34 x4 a35
(2) (2) (2)
a43 x3 a44 x4 a45
gdzie współczynniki:
(2) (1) (1)
aij aij ai(1)b2j i 3,4 j 3,4,5
2
Pierwsze z równań znowu dzielimy przez pierwszy współczynnik:
(2) (2)
x3 b34 x4 b35
gdzie:
(2)
a3j
(2)
b3j j 4,5
(2)
a33
i za pomocą tego równania eliminujemy zmienną x3:
(3) (3)
a44 x4 a45
gdzie:
(3) (2) (2) (2)
a4j a4j a43 b3j
W ten sposób przekształciliśmy wyjściowy układ równań do postaci:
x1 b12x2 b13x3 b14x4 b15
(1) (1) (1)
x2 b23 x3 b24 x4 b25
(2) (2)
x3 b34 x4 b35
(3) (3)
a44 x4 a45
który możemy już łatwo rozwiązać obliczając kolejno:
(3) (3)
x4 a45 a44
(2) (2)
x3 b35 b34 x4
(1) (1) (1)
x2 b25 b24 x4 b23 x3
x1 b15 b14x4 b13x3 b12x2
Jak widać metoda Gaussa składa się z dwóch części: w pierwszej z nich następuje eliminacja kolejnych
niewiadomych  w efekcie następuje przekształcenie macierzy współczynników w macierz trójkątną; w drugiej części
następuje obliczenie kolejnych niewiadomych.
Wszystkie układy, które drogą zamiany kolejności równań mogą zostać sprowadzone do przypadku macierzy
trójkątnej można rozwiązać bardzo łatwo  stosując tylko drugą część metody Gaussa.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 4
Jak łatwo zauważyć opisana powyżej metoda ma pewne mankamenty:
 jeśli któryś z kolejnych elementów macierzy wykorzystywanych do normalizowania wierszy (aii(i)) będzie
wynosić zero to nie uda się przeprowadzić obliczeń,
 jeśli któryś z powyższych elementów będzie wystarczająco mały to spowoduje to wystąpienie dużych
błędów,  dowiedziono, że przedstawiona powyżej metoda jest numerycznie niestabilna  zawsze należy
stosować przedstawioną poniżej modyfikację  wybór elementu głównego.
Z powyższych powodów stosuje się raczej pewną modyfikację metody Gaussa, polegającą na zmianie kolejności
równań i zmiennych w taki sposób, aby za każdym razem do normalizacji kolejnych wierszy wykorzystywać
największy co do wartości element przekształcanej macierzy.
Najpierw zamieniając kolejność kolumn macierzy wyjściowej doprowadzamy do spełnienia warunku:
(0) (0)
a11 aij i,j 1,2, ,n
w rezultacie element a11(0) ma największą wartość bezwzględną.
Następnie, tak samo jak poprzednio, dzielimy pierwsze równanie przez a11(0) i eliminujemy x1 z pozostałych równań.
W nowootrzymanym układzie równań znowu wyszukujemy największy element:
(1) (1)
a22 aij i,j 2,3, n
i tak przestawiamy wiersze i kolumny aby znalazł się on na pierwszej pozycji. Koniecznie należy przy tym pamiętać o
odpowiednim przestawieniu równań otrzymanych wcześniej (oczywiście dotyczy to przypadku gdy przestawiamy
kolumny macierzy) oraz (w przypadku przestawiania wierszy) o odpowiednim przestawieniu wierszy wektora
niewiadomych x.
Dalsze obliczenia sÄ… analogiczne jak poprzednio, przy czym zawsze przy okazji wyboru kolejnego maksymalnego
elementu pamiętać należy o odpowiednim przestawieniu poprzednich (już gotowych) równań (jeśli dokonywano
zamiany kolumn) i wierszy (również wektora x).
W praktyce najczęściej stosuje się tzw. częściowy wybór elementu głównego  tzn. bada się tylko aktualną
pierwszą kolumnę i dokonuje ewentualnej zamiany tylko dla wierszy. Jest to znacznie prostsze od wyboru pełnego, a
najczęściej wystarcza.
Dokonywanie wyboru jest zbędne w przypadku macierzy z dominującą przekątną:
n
aii aij i 1,2, ,n
j 1
j i
Szacuje się, że w metodzie Gaussa ilości koniecznych do wykonania mnożeń i dzieleń oraz dodawań i odejmowań
są rzędu n3.
Przykład: Rozważmy układ równań:
2x1 7x2 4x3 9
x1 9x2 6x3 1
3x1 8x2 5x3 6
kolejne przekształcenia dają następujące układy:
7 9
x1 x2 2x3
2 2
25 7
x2 8x3
2 2
5 39
x2 11x3
2 2
7 9
x1 x2 2x3
2 2
16 7
x2 x3
25 25
47 94
x3
5 5
RozwiÄ…zaniem sÄ…: x1 = 4, x2 = 1, x3 = 2.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 5
Metoda eliminacji Jordana (Gaussa-Jordana, eliminacji zupełnej)
W metodzie tej, będącej modyfikacją zwykłej metody Gaussa, łączy się przekształcenie macierzy opisującej układ
równań z obliczaniem wyników.
Algorytm jest następujący:
akj
akj j n m, n m 1, ,k
akk
k 1,2, ,n
m,n
aij aij aikakj j i n 1,2, ,nm , ,k
i !n
Przykład. Dla ilustracji rozpatrzymy ten sam układ równań:
2x1 7x2 4x3 9
x1 9x2 6x3 1
3x1 8x2 5x3 6
Macierz współczynników równania, wraz z wyrazami wolnymi ma postać:
2 7 4 9
1 9 6 1
3 8 5 6
Tak jak poprzednio pierwszy wiersz macierzy normalizujemy za pomocą pierwszego jej elementu, a następnie
od wszystkich następnych wierszy odejmujemy otrzymany właśnie pierwszy rząd pomnożony przez pierwszy
współczynnik w danym wierszu  w rezultacie w pozostałych wierszach współczynniki ai1 wynoszą zero.
7 9
1 2
2 2
25 7
0 8
2 2
5 39
0 11
2 2
Następnie normalizujemy wiersz drugi, dzieląc wszystkie jego elementy przez 25/2 i redukujemy pozostałe
równania (także już przekształcone równanie pierwsze) odejmując od nich wiersz drugi pomnożony przez
odpowiedni współczynnik ai2:
6 88
1 0
25 25
16 7
0 1
25 25
47 94
0 0
5 5
Na koniec (w tym przykładzie) normalizujemy wiersz trzeci dzieląc jego współczynniki (teraz już tylko wyraz
wolny) przez 47/5 i redukujemy pozostałe równania tak jak poprzednio. Otrzymujemy:
1 0 0 4
0 1 0 1
0 0 1 2
Ostatnia kolumna (w której na początku umieściliśmy wyrazy wolne) zawiera rozwiązanie.
Szacuje siÄ™, że w metodzie Jordana ilość mnożeÅ„ i dzieleÅ„ oraz dodawaÅ„ i odejmowaÅ„ sÄ… rzÄ™du ½n3, a wiÄ™c nieco
więcej niż w metodzie Gaussa, przewaga metody Jordana polega na nieco mniejszej zajętości pamięci.
Obie powyżej opisane metody można zastosować do rozwiązywania układów równań o większej liczbie  prawych
stron . W wielu przypadkach mamy do czynienia z zagadnieniami, w których macierz współczynników wynika z
charakteru równań opisujących modelowane zjawisko, natomiast postać prawej strony układu równań wynika np. z
warunków brzegowych. Często możemy spotkać się z sytuacją gdy będziemy chcieli rozwiązać kilka (lub nawet wiele)
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 6
układów równań różniących się między sobą tylko prawymi stronami. W takim przypadku nasza macierz będzie mieć
odpowiednio większą liczbę kolumn (po prawej stronie dopiszemy wszystkie  dodatkowe prawe strony)  natomiast
sam algorytm nie zmieni się. Oczywiście równoczesne rozwiązywanie kilku układów równań daje wyraz
nÄ…
oszczędność czasu w stosunku do rozwiązywania każdego zagadnienia osobno (jest jednak możliwe tylko w przypadku
układów równań różniących się tylko prawymi stronami).
Metoda Cholesky'ego (Banachiewicza)
Rozważmy znowu układ równań (zakładamy, że macierz współczynników jest symetryczna i dodatnio określona1).
Jeśli metodę tę można zastosować to jest około 2 razy szybsza od opisanych wcześniej.
a11 a12 a1n x1 a1,n 1
a21 a22 a2n x2 a2,n 1
Ax w
a a a x a
n1 n2 nn n n,n 1
i przekształcimy macierz A do postaci BC, gdzie:
b11 0 0
1 c12 c1n
b21 b22 0
0 1 c2n
B i C
b b b
0 0 1
n1 n2 nn
Dzięki czemu mamy:
Ax w Cx y
BCx w
A BC By w
Elementy macierzy B i C obliczane są następująco:
bi1 ai1
i 1
bji aji bjkcki i j>1
k 1
oraz:
a1j
c1j
b11
i 1
1
cij aij bikckj 1bii k 1
Obliczenia prowadzimy naprzemiennie wierszami i kolumnami: najpierw kolumna 1-sza macierzy B, po niej 1-szy
wiersz macierzy C; następnie 2-ga kolumna macierzy B i po niej 2-gi wiersz macierzy C itd.
Mając dane macierze B i C możemy znalez
ć poszukiwany wektor rozwiązań x  rozwiązując kolejno dwa
równania:
B y = w i C x = y
Ponieważ obie macierze (B i C ) są trójkątne rozwiązanie znajduje się następująco:
1
Rzeczywista symetryczna macierz A (n×n) jest dodatnio okreÅ›lona jeÅ›li istnieje rzeczywista,
T T
nieosobliwa macierz M taka, że A MM . Albo: jeśli dla każdego wektora x zachodzi x Ax 0. Albo:
a11
wszystkie wyznaczniki postaci , dla i od 1 do n, sÄ… dodatnie (jest to tzw. kryterium Sylwestra).
aii
Wszystkie wartości własne macierzy dodatnio określonej są dodatnie.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 7
a1,n 1
y1
b11
i 1
1
yi ai,n 1 bikyk i >j
bii
k 1
a potem:
x yn
n
n
xi yi cikxk ik i 1
Dla przykładu metodę tę zastosujemy do następującego układu równań:
3x1 x2 x3 2x4 6
5x1 x2 3x3 4x4 12
2x1 x3 x4 1
x1 5x2 3x3 3x4 3
macierze B i C sÄ… postaci:
1 1 2
3 0 0 0
1
3 3 3
3 1 1 2
2
5 2 0 0
1 1
3
0 1
5 1 3 4
2 4
2
2 2 0
2 0 1 1
1
3
0 0 1 1
4
1 5 3 3
1 1
1 5 6 2
0 0 0 1
3 2
ostatecznie obliczamy wektory y i x (strzałki odpowiadają kolejności obliczeń):
2
1
3
1
4
y ; x
3
2
1
4
3
3
Metoda przemiatania
Metoda ta ma zastosowanie do rozwiązywania układów równań, których współczynniki tworzą macierz
trójdiagonalną (tzn. mającą elementy niezerowe tylko na trzech środkowych przekątnych). Znana jest również pod
nazwÄ… the sweep method albo <5B>4 ?@>3>=:8.
b1x1 c1x2 d1
a2x1 b2x2 c2x3 d2
a3x2 b3x3 c3x4 d3
aixi 1 bixi cixi 1 di
a x b x c x d
n 1 n 2 n 1 n 1 n 1 n n 1
a x b x d
n n 1 n n n
Rozwiązanie takiego układu równań przeprowadza się w dwóch etapach. W etapie pierwszym oblicza się kolejno
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 8
współczynniki ² i Å‚:
²1 b1
d1
Å‚1
²1
aici 1
²i bi i 2,3, ,n
²i 1
di aiłi 1
Å‚i i 2,3, ,n
²i
a w drugim na ich podstawie znajduje siÄ™ rozwiÄ…zania:
x Å‚n
n
cixi 1
xi Å‚i i n 1,n 2, ,1
²i
Jest to metoda bardzo szybka, zaś macierze trójdiagonalne na szczęście dość często występują w zagadnieniach
obliczeń numerycznych (zwłaszcza podczas rozwiązywania układów równań różniczkowych cząstkowych metodą
różnic skończonych).
Ponieważ współczynniki ²i sÄ… obliczane jako sumy (różnice), a nastÄ™pnie pojawiajÄ… siÄ™ w mianowniku wzoru na
kolejny współczynnik łi należy liczyć się z możliwością wystąpienia dzielenia przez wartość bliską zeru (albo nawet
przez zero)  mogłoby to uniemożliwić obliczenia. Sytuacja taka na pewno nie wystąpi jeśli macierz współczynników
równania jest diagonalnie dominująca (tzn. moduł wartości diagonalnej jest większy lub równy od sumy modułów obu
pozostałych współczynników).
Metody przybliżone  iteracyjne
Metody iteracyjne określane są jako przybliżone z tego względu, że dają rozwiązanie dokładne dopiero po
nieskończonej liczbie kroków  kolejne wyniki są coraz dokładniejszymi przybliżeniami  prawdziwego rozwiązania.
Są one właściwie jedynymi metodami skutecznymi w przypadku układów równań nieliniowych (podobnie jak w
przypadku nieliniowych funkcji jednej zmiennej). Można jednak stosować je także do układów równań liniowych.
Metody iteracyjne mają pewną przewagę nad metodami dokładnymi  polega ona na tym, że nie trzeba
przejmować się przenoszeniem błędów pomiędzy kolejnymi iteracjami. W przypadku metod dokładnych błędy obcięcia
popełnione w dowolnym momencie kumulują się i mogą nawet bardzo znacznie zniekształcić ostateczne rozwiązanie.
Tymczasem w metodach iteracyjnych, nawet gdy kolejne przybliżenia są obarczone dodatkowo błędami wynikającym z
błędów obcięcia itp., to są one traktowane jako kolejne podstawy do wyznaczenia następnego, dokładniejszego
przybliżenia. Z tego powodu mówimy, że metody przybliżone (iteracyjne) można stosować do udokładniania wyników
otrzymanych metodami dokładnymi.
Układy równań liniowych
Metoda iteracyjna (prosta)
Można powiedzieć, że metoda ta jest adaptacją metody iteracyjnej przedstawionej już wcześniej dla równań jednej
zmiennej. Załóżmy, że układ równań:
A x b
może zostać przekształcony do postaci:
x C x f
gdzie C jest pewną, nieznaną jeszcze, macierzą, a f jest przekształconym wektorem prawych stron prawych b (w
szczególności może to być niezmieniony wektor b  zależnie od przekształceń które zastosowano aby przejść od A do
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 9
C).
WychodzÄ…c od dowolnie wybranego wektora:
(0)
x1
(0)
x2
(0)
x
(0)
xn
można określić proces iteracyjny:
(i 1) (i )
x Cx f i 0,1,2,
czyli:
(i (i (i (i
x1 1) c11x1 ) c12x2 ) c1nxn ) f1
(i (i (i (i
xn 1) c x1 ) c x2 ) c xn ) fn
n1 n2 nn
Przeprowadzając kolejne iteracje otrzymamy ciąg wektorów x(i), który powinien być zbieżny do rozwiązania
dokładnego.
Warunkiem zbieżności tego ciągu jest spełnienie jednego z następujących warunków:
n
cij Ä… < 1 i 1,2, ,n
j 1
n
lub: cij ² < 1 j 1,2, ,n
i 1
Należy zwrócić uwagę, że spełnienie któregoś z wymienionych warunków gwarantuje zbieżność, natomiast jego
niespełnienie nie gwarantuje braku zbieżności.
W przypadku spełnienia któregoś z tych warunków możliwe jest też wyznaczenie błędu jakim obarczone jest k-te
przybliżenie rozwiązania:
Ä…
xi xi(k) max xj(k) xj(k 1)
1 Ä… j 1,2, ,n
lub:
n
²
xi xi(k) xj(k) xj(k 1)
1 ²
j 1
Problem może stanowić znalezienie odpowiedniej macierzy C (spełniającej warunki zbieżności metody).
Najczęściej stosuje się jeden z dwóch sposobów:
1) Jeśli diagonalne elementy macierzy są niezerowe to układ równań można przepisać w postaci:
1
x1 b1 a12x2 a13x3 a1nx
n
a11
1
x2 b2 a21x1 a23x3 a2nx
n
a22
1
x bn a x1 a x2 a x
n n1 n2 n,n 1 n 1
a
nn
w tym przypadku elementy macierzy C określa się następująco:
aij
cij i j cii 0
aii
a wektora f:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 10
bi
fi
aii
Warunki zbieżności będą wówczas spełnione jeśli zachodzi:
aii > aij i 1,2, ,n
i j
Jest to tzw. macierz z dominujÄ…cÄ… przekÄ…tnÄ….
2) Druga metoda ma zastosowanie gdy elementy na przekątnej macierzy współczynników równania są bliskie jedności,
a wszystkie pozostałe są znacznie od jedności mniejsze. Wówczas równania przekształcamy przenosząc xi na lewą
stronę, a po prawej pozostawiając wszystkie (włącznie z wyciągniętą na lewo zmienną xi ze współczynnikiem
zmniejszonym o 1) zmienne i wyrazy wolne. Taki układ już bez dalszych przekształceń jest wykorzystywany do
iteracji. Warunki zbieżności są w takim przypadku z reguły spełnione.
Metoda Seidel'a (Gaussa-Seidel'a)
Jest to pewna modyfikacja prostej metody iteracyjnej. Polega ona na tym, że do obliczania kolejnego przybliżenia
wykorzystuje się  najnowsze dostępne wartości zmiennych. Np. do obliczenia (k+1)-go przybliżenia zmiennej xi
wykorzystuje się wszystkie nowoobliczone wartości zmiennych: x1(k+1), x2(k+1), ..., xi 1(k+1) oraz  starsze wartości
następnych zmiennych: xi(k), ..., xn(k). Inaczej mówiąc zakłada się, że układ równań jest spełniony przez taki właśnie
wektor zawierający mieszaninę dwóch kolejnych przybliżeń rozwiązania:
i N
A x b aijxj(k 1) aijxj(k) bi i 1, ,N
j i j i 1
Oczywiście dopóki nie osiągniemy  docelowego wektora będącego rozwiązaniem dokładnym równość po prawej
stronie zapisanego powyżej równania nie jest spełniona. Powyższe równanie można przekształcić do postaci:
i 1 N
1
xi(k 1) bi aijxj(k 1) aijxj(k)
aii j 1
j i 1
Stosując tę metodę uzyskuje się szybszą zbieżność niż w metodzie prostej (chociaż nie zawsze), a dodatkową
korzyścią jest możliwość posługiwania się tylko jedną macierzą, w której nowoobliczone elementy natychmiast
zastępują stare  daje to oszczędność pamięci.
Zapisany powyżej wzór można przepisać nieco inaczej  w nawiasie możemy dodać i odjąć składnik aiixi(k) 
dostaniemy wówczas nieco inny wzór iteracyjny na xi(k 1):
i 1 N
1
xi(k 1) xi(k) bi aijxj(k 1) aijxj(k)
aii j 1
j i
którego drugi składnik można traktować jak poprawkę dla kolejnej obliczanej wartości.
Metodę tę można znacznie krócej zapisać w notacji macierzowej. Zakładamy, że macierz A można zapisać w
postaci sumy trzech macierzy:
A = L + D + U
gdzie D jest macierzą diagonalną, której jedynymi elementami niezerowymi są elementy macierzy A leżące na jej
głównej przekątnej, L jest dolną macierzą trójkątną zawierającą elementy macierzy A leżące wyłącznie pod jej główną
przekątną, a U jest górną macierzą trójkątną zawierającą elementy macierzy A leżące wyłącznie nad jej główną
przekątną. Jeżeli teraz następująco zdefiniujemy macierz B oraz wektor c :
B =  (L+D) 1U oraz c = (L+D) 1b
to przedstawiony algorytm można zapisać w takiej postaci:
x (i+1) = Bx (i) + c
Metoda  nadrelaskacji
Z przedstawionym powyżej przekształceniem wiąże się kolejna metoda, zwana  nadrelaksacją (over-relaxation) 
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 11
wyznaczonÄ… poprzednio poprawkÄ™ powiÄ™kszamy o jakiÅ› czynnik É:
i 1 N
É
xi(k 1) xi(k) bi aijxj(k 1) aijxj(k)
aii j 1
j i
Zmiana ta, wydawałoby się niczym nie uzasadniona, powoduje (może powodować) kolejne przyspieszenie
zbieżności oraz większą  odporność na złe uwarunkowanie macierzy współczynników układu równań.
Przykłady
Dla przykładu rozwiążemy kilka prostych układów równań. Wszystkie one mają rozwiązania dokładne postaci:
x1 1 ; x2 2 ; x3 3 ; x4 4, które można oczywiście bez problemu znalez
ć innymi metodami, ale tu posłuży nam
jako dobra ilustracja.
x1 16
9 1 1 2
x2 32
1 8 3 2
Pierwszy układ ma postać:  jak widać jest to układ z dominującą główną
2 1 8 3 x3 12
3 2 1 7
x4 26
przekÄ…tnÄ….
x1 16
9 1 1 2
x2 32
1 8 3 2
W drugim:  zmieniony został ostatni wiersz  nie mamy już dominacji
1 1 8 3 x3 11
3 2 11 7
x4 4
głównej przekątnej.
Układ trzeci jest analogiczny jak drugi, z tą tylko różnicą, że brak dominacji głównej przekątnej wynika z postaci
x1 56
9 1 1 12
x2 32
1 8 3 2
pierwszego wiersza: .
2 1 8 3 x3 12
3 2 1 7
x4 26
Układ czwarty stanowi połączenie zmian wprowadzonych w układach drugim i trzecim:
x1 56
9 1 1 12
x2 32
1 8 3 2
.
2 1 8 3 x3 12
3 2 11 7
x4 4
x1 56
9 1 1 12
x2 62
1 8 13 2
W układzie piątym dodatkowo zmieniono wiersz drugi: .
2 1 8 3 x3 12
3 2 11 7
x4 4
I w końcu w układzie szóstym żadna z wartości z głównej przekątnej nie jest dominująca w swoim wierszu (jest to
więc układ najgorszy z punktu widzenia możliwości rozwiązania iteracyjnego):
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 12
x1 56
9 1 1 12
x2 62
1 8 13 2
.
2 11 8 3 x3 8
3 2 11 7
x4 4
Błędy rozwiązania (wektor czteroelementowy) obliczone jako normę różnicy wektorów  aktualnego przybliżenia i
rozwiązania dokładnego. Dla każdego układu przedstawione będą dwa wykresy  początkowe 10 oraz 300 ietracji.
Układ pierwszy:
Widać, że poczynajÄ…c od wartoÅ›ci É ok. 1.5 metoda nadrelaksacji nie dziaÅ‚a (wyniki  rozbiegajÄ… siÄ™ ), najszybsza
zbieżność wsytÄ™puje jednak dla É=1 (czyli bez nadrelaksacji).
Układ drugi:
Tym razem widzimy, że najszybszÄ… zbieżność otrzymujemy dla É=1.25 (chociaż wartość 1 jest niewiele gorsza).
Metoda jest rozbieżna tylko dla jednej (spoÅ›ród badanych) wartoÅ›ci É.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 13
Układ trzeci:
W przypadku ukÅ‚adu trzeciego É=1 nadal jest  na drugim miejscu , jednak najlepszÄ… zbieżność uzyskujemy
dla wartoÅ›ci É mniejszej od 1  É=0.75. Metoda jest zbieżna tylko dla É 1 (spoÅ›ród badanych wartoÅ›ci).
Układ czwarty:
W układzie czwartym już w dwóch wierszach maksymalne wartości współczynników występują poza główną
przekÄ…tnÄ…, jednak metoda nadal dziaÅ‚a (jest rozbieżna tylko dla dwóch maksymalnych badanych wartoÅ›ci É).
Układ piąty:
W układzie piątym maksymalne wartości występują poza główną przekątną w trzech wierszach, nadal jednak
uzyskujemy zbieżność dla wszystkich (poza jednÄ…) badanych wartoÅ›ci É.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 14
Układ szósty:
W układzie szóstym w ogóle nie można już mówić o dominacji głównej przekątnej. Metoda jest rozbieżna dla
wiÄ™kszoÅ›ci badanych wartoÅ›ci É, ale dla wartoÅ›ci mniejszych od 1 nadal jest zbieżna.
Układy równań nieliniowych
Metoda iteracyjna
Załóżmy, że mamy układ dwóch równań nieliniowych:
f1 x,y 0
f2 x,y 0
musimy przekształcić go do postaci:
x F1 x,y
y F2 x,y
wybrawszy jakieś wartości x0 i y0 można posłużyć się wzorami iteracyjnymi:
i 1 i i
x F1 x ,y
i 1,2,3,
i 1 i i
y F2 x ,y
Jeśli spełnione są warunki:
a) funkcje F1 i F2 mają ciągłe pochodne w pewnym obszarze R,
b) przybliżenia początkowe (x0 i y0) oraz wszystkie kolejne leżą w tymże obszarze R,
c) w całym obszarze R spełnione są nierówności:
F1 F2
q1 < 1
x x
F1 F2
q2 < 1
y y
to kolejne przybliżenia zbiegają się do rozwiązań układu.
Osobnym problemem jest sposób znalezienia funkcji F1 i F2. Jeden z możliwych sposobów jest następujący
(niestety  dla układu dwóch równań):
Funkcje F1 i F2 mają postać:
F1 x,y x Ä…f1 x,y ²f2 x,y
F2 x,y y Å‚f1 x,y ´f2 x,y
Ä…´ ²Å‚
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 15
gdzie współczynniki Ä…, ², Å‚ i ´ stanowiÄ… przybliżone rozwiÄ…zania nastÄ™pujÄ…cego ukÅ‚adu równaÅ„:
f1 x0,y0 f2 x0,y0
1 Ä… ² 0
x x
f1 x0,y0 f2 x0,y0
Ä… ² 0
y y
f1 x0,y0 f2 x0,y0
Å‚ ´ 0
x x
f1 x0,y0 f2 x0,y0
1 Å‚ ´ 0
y y
Jeśli tylko pochodne cząstkowe f1 i f2 nie zmieniają się zbyt szybko w otoczeniu punktu (x0,y0) to warunki zbieżności
będą spełnione.
Przykład. Rozważmy układ równań:
y x
1
f1 x,y sin xy 0
2
4Ä„ 2
1 ey
2x
f2 x,y 1 e e 2ex 0
4Ä„ Ä„
Przepiszmy je w postaci właściwej dla obliczeń iteracyjnych:
y
x F1 x,y sin xy
2Ä„
1
2x 1
y F2 x,y 2Ä„x Ä„ e 1
4
Jeśli wybierzemy następujące wartości początkowe: x0 = 0.4 i y0 = 3.0 to otrzymamy następujący ciąg rozwiązań:
x i : 0.455 0.498 0.505 0.499 0.499 0.500
y i : 3.037 3.107 3.141 3.144 3.142 3.142
Jeśli posłużymy się metodą Gaussa-Seidel'a (tzn. w każdym obliczeniu będziemy wykorzystywać najnowsze
dostępne wartości) to przy tych samych wartościach wyjściowych otrzymamy taki ciąg rozwiązań (gdy najpierw
obliczamy x, a potem obliczając y wykorzystujemy właśnie obliczoną wartość x):
x i : 0.455 0.493 0.500 0.500
y i : 3.107 3.138 3.142 3.142
Jak widać zbieżność będzie nieco szybsza. Możemy też najpierw obliczać y i wykorzystywać jego wartość do
obliczenia x:
x i : 0.454 0.493 0.500 0.500 0.500
y i : 3.037 3.107 3.138 3.142 3.142
Możliwy jest też odwrotny wybór równań  tzn. z pierwszego z nich obliczamy y, a z drugiego x:
y 2Ä„ sin xy x
y
1 1
2x 1
x e 1
2 8Ä„
2Ä„
x i : 0.394 0.444 0.525 0.579 0.519 0.438 0.407 0.439 0.509
y i : 3.343 3.607 3.489 2.767 2.640 2.895 3.244 3.530 3.526
x i : 0.187 0.241 0.107 0.034  0.181  0.302  0.097  0.348  0.159
y i : 3.343 2.497 2.044 0.688  0.067 1.217  0.359 0.831 0.396
x i : 0.394  0.488  0.790 0.062  0.042  0.827  0.386 0.445
y i :  2.513  2.476 3.066 4.967  3.280 0.262 5.193 2.428
Przebieg obliczeń dla obu przypadków ilustrują poniższe wykresy (w przypadku dolnego zamieniono kolejność
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 16
równań):
Przebieg zbiegania się rozwiązania można również analizować obserwując zmiany normy kolejnych przybliżeń
rozwiÄ…zania:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 17
Linia pozioma na takim wykresie oznacza osiągnięcie wartości stabilnej  czyli znalezienie rozwiązania (jakiegoś
 niekoniecznie akurat tego, na którym na zależało).
Przykład zastosowania prostej metody iteracyjnej
Dany jest układ poziomych rur, o znanych długościach i średnicach, połączonych w n węzłach. Dla niektórych
węzłów znane jest ciśnienie. Problem polega na znalezieniu ciśnień we wszystkich nieznanych węzłach oraz prędkości
przepływów (w m3/s) pomiędzy wszystkimi węzłami (i ich kierunki).
Wiadomo, że spadek ciśnienia w poziomej rurze jest dany wzorem:
1 L
2
"p fM Áu
2 D
gdzie fM jest tzw. współczynnikiem Moody'ego, Á jest gÄ™stoÅ›ciÄ… cieczy, u jest prÄ™dkoÅ›ciÄ… Å›redniÄ…, a D i L sÄ…
odpowiednio średnicą i długością rury. Wykorzystując objętościową prędkość przepływu Q można zapisać:
2
8fM ÁQ L
"p
5
Ä„2D
przy czym wszystkie stałe (niezależne od konkretnej rury) można zebrać razem (ozn. C ):
2
LijQij
pi pj C
5
Dij
wówczas dla każdej pary węzłów otrzymamy równanie:
2
pi pj cijQij
gdzie:
CLij
cij
5
Dij
ostatnie równanie można przekształcić tak aby automatycznie miało prawidłowy znak:
Qij pi pj 1
cij pi pj
W każdym węz
le, w którym nie określono ciśnienia (np. o numerze j), suma przepływów pochodzących z
wszystkich sąsiednich węzłów musi wynosić zero:
Qij pi pj 1 0
cij pi pj
i j i j
(j  węzeł badany, i  węzły sąsiednie) i te właśnie równania, zapisane dla wszystkich węzłów o nieznanym ciśnieniu
stanowią układ, który będziemy rozwiązywać.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 18
Wzór iteracyjny będzie mieć postać:
aijpi
i j
pj
aij
i j
gdzie (dla skrócenia zapisu):
1
aij cij pi pj 2
Dane do obliczeń:
1 2 3
fM = 0.056 L12 = L23 = 45 m
Á = 820 kg/m3 L14 = L45 = L25 = 30 m
p1 = 3.8 atm D12 = D23 = 3"
4 5
p3 = 0 D14 = D45 = D25 = 4"
Metoda Newtona-Raphsona
Załóżmy, że dany jest następujący układ równań nieliniowych:
f1 x1,x2, ,x 0
n
f2 x1,x2, ,x 0
n
fn x1,x2, ,x 0
n
Aby zapisać wzory w postaci macierzowej najpierw zdefiniujemy następujące macierze:
Ç x1,x2, ,x
n
f Ç f1 Ç ,f2 Ç , ,fn Ç
Zdefiniujemy:
fi Ç
fij Ç
xj
i określimy macierz pochodnych cząstkowych:
Åš Ç fij Ç
dla: 1 i n, 1 j n.
Przy takim zapisie, zakładając że wybraliśmy wektor zawierający początkowe przybliżenia wszystkich
niewiadomych (Ç0 = (x10,x20,...,xn0)) możemy zapisać wzór iteracyjny:
Çk 1 Çk ´k
gdzie ´k jest rozwiÄ…zaniem nastÄ™pujÄ…cego ukÅ‚adu równaÅ„ algebraicznych:
Åš Çk ´k f Çk
gdzie k jest oczywiście numerem iteracji.
Metoda ta ma pewną przewagę nad poprzednią  warunek zbieżności wymaga tylko aby składowe wektora
pochodnych Ś były ciągłe w otoczeniu pierwiastka ą (ą jest oczywiście wektorem), wyznacznik det(Ś(ą)) 0 oraz by
przybliżenie poczÄ…tkowe Ç0 byÅ‚o  dostatecznie blisko pierwiastka Ä….
Przykład. Rozpatrzymy znowu ten sam układ równań:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 19
y x
1
f1 x,y sin xy 0
2
4Ä„ 2
1 ey
2x
f2 x,y 1 e e 2ex 0
4Ä„ Ä„
Obliczymy odpowiednie pochodne:
f1 1 y cos xy
x 2 2
f1 1 x cos xy
y 4Ä„ 2
f2
1
2x
2e 2 e
x 2Ä„
f2 e
y Ä„
Kolejne poprawki zmiennych x i y będziemy obliczać rozwiązując układ równań:
f1 f1
k k
´x ´y f1 x ,y
x k k y k k
x ,y x ,y
f2 f2
k k
´x ´y f2 x ,y
x k k y k k
x ,y x ,y
(oczywiście w każdym przypadku chodzi o wartość funkcji lub wartość jej odpowiedniej pochodnej w punkcie
będącym aktualnym przybliżeniem rozwiązania).
Przy takich samych jak poprzednio (w metodzie iteracyjnej) danych poczÄ…tkowych (tzn. 0.4 i 3) otrzymamy
następujący ciąg rozwiązań:
x i :  0.431  0.243  0.258  0.261  0.261
y i : 1.693 1.181 0.538 0.623 0.623
czyli inną parę rozwiązań. Jednak przyjęcie wartości początkowych x0 = 0.6 i y0 = 3.0 prowadzi do rozwiązań
uzyskanych poprzednio:
x i : 0.502 0.500 0.500 0.500 0.500 0.500 0.500 0.500
y i : 3.082 3.094 3.105 3.113 3.124 3.128 3.131 3.133
Poniżej znajdują się dwie ilustracje metody Newtona-Raphsona dla punktów początkowych (0.4,3.0) i (2.0,1.0):
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Rozwiązywanie układów równań 20
Także i w tym przypadku spojrzymy jeszcze na przebieg zmian norm wektorów będących kolejnymi
przybliżeniami rozwiązania:
Dla drugiego punktu początkowego na wykresie norm pokazałem 50 iteracji (podczas gdy na wykresie rozwiązań
było ich 20)  dzięki temu wyraz
niej widać, że osiągamy zbieżność.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej