mnf5'06.wpd

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

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:

którego rozwi zaniem b dzie oczywi cie:

Je li bł dy wyst puj ce w wektorze prawych stron oznaczymy jako b dziemy mogli zapisa :

b ;

A x

b ; x

b ;

Wynika st d (zakładaj c, e norm liczymy jako pierwiastek z sumy kwardatów):

A x

Mo na wi c zapisa nast puj ce nierówno ci:

;

Wprowadzimy teraz oznaczenie:

które umo liwi nam poł czenie zapisanych powy ej nierówno ci:

Nierówno ta okre la „granice” bł du (wynikaj cego z bł dów ) jakim obarczone mo e by rozwi zanie . Wida ,

e „im wi ksze K tym gorzej”. Macierze, dla których K przybiera du e warto ci nazywane s „ 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”.

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

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:

to nale y obliczy jego wyznacznik główny :

oraz wszystkie pozostałe wyznaczniki postaci (np.):

;

Po czym wyznaczy niewiadome z wzorów:

; x

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):

załó my, e element a

0, je li podzielimy pierwsze równanie układu przez a

otrzymamy:

gdzie:

j 2,3,4,5

Wykorzystuj c powy sze równanie mo na łatwo wyeliminowa zmienn x

z pozostałych równa . W tym celu

mno ymy je odpowiednio przez a

, a

i a

i odejmujemy od odpowiednich równa . Otrzymujemy układ trzech

równa :

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

(1)

gdzie współczynniki:

(1)

i 2,3,4

j 2,3,4,5

Znowu pierwsze równanie układu dzielimy przez jego pierwszy współczynnik (a

(1)

) i otrzymujemy:

(1)

gdzie:

(1)

j 3,4,5

W taki sam jak poprzednio sposób eliminujemy zmienn x

i otrzymujemy układ:

(2)

gdzie współczynniki:

(2)

(1)

i 3,4

j 3,4,5

Pierwsze z równa znowu dzielimy przez pierwszy współczynnik:

(2)

gdzie:

(2)

j 4,5

i za pomoc tego równania eliminujemy zmienn x

(3)

gdzie:

(3)

(2)

W ten sposób przekształcili my wyj ciowy układ równa do postaci:

(1)

(2)

(3)

który mo emy ju łatwo rozwi za obliczaj c kolejno:

(3)

(2)

(1)

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.

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Jak łatwo zauwa y opisana powy ej metoda ma pewne mankamenty:

— je li który z kolejnych elementów macierzy wykorzystywanych do normalizowania wierszy (a

(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)

i,j 1,2, ,n

w rezultacie element a

(0)

ma najwi ksz warto bezwzgl dn .

Nast pnie, tak samo jak poprzednio, dzielimy pierwsze równanie przez a

(0)

i eliminujemy x

z pozostałych równa .

W nowootrzymanym układzie równa znowu wyszukujemy najwi kszy element:

(1)

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 :

j 1

j i

i 1,2, ,n

Szacuje si , e w metodzie Gaussa ilo ci koniecznych do wykonania mno e i dziele oraz dodawa i odejmowa

s rz du n

Przykład: Rozwa my układ równa :

kolejne przekształcenia daj nast puj ce układy:

7
2

9
2

7
2

5
2

11x

7
2

9
2

16
25

Rozwi zaniem s : x

= 4, x

= 1, x

= 2.

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

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:

j n m, n m 1, ,k

j n m,n m , ,k

i 1,2, ,n i !n

k 1,2, ,n

Przykład. Dla ilustracji rozpatrzymy ten sam układ równa :

Macierz współczynników równania, wraz z wyrazami wolnymi ma posta :

7 4 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 a

wynosz zero.

7
2

9
2

7
2

5
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 a

88
25

0 1

16
25

0 0

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:

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 ½n

, 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)

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

Rzeczywista symetryczna macierz A (n×n) jest dodatnio okre lona je li istnieje rzeczywista,

nieosobliwa macierz M taka, e

. Albo: je li dla ka dego wektora x zachodzi

. Albo:

wszystkie wyznaczniki postaci

, dla i od 1 do n, s dodatnie (jest to tzw. kryterium Sylwestra).

Wszystkie warto ci własne macierzy dodatnio okre lonej s dodatnie.

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

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 wyra 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 lona

Je li metod t mo na zastosowa to jest około 2 razy szybsza od opisanych wcze niej.

1,n 1

2,n 1

n,n 1

i przekształcimy macierz

A do postaci BC, gdzie:

1 c

0 1

0 0

Dzi ki czemu mamy:

BCx

Elementy macierzy

B i C obliczane s nast puj co:

i 1

k 1

i j>1

oraz:

i 1

k 1

1<i<j

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 znale 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:

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

1,n 1

i,n 1

i 1

k 1

i >j

a potem:

k i 1

i<n

Dla przykładu metod t zastosujemy do nast puj cego układu równa :

macierze

B i C s postaci:

1 2

5 1

5 3

0 0

5 2

2
3

0 0

2
3

2 0

1
3

6 2

1
2

1
3

2
3

0 1

1
2

1
4

0 0 1

1
4

0 0 0

ostatecznie obliczamy wektory

y i x (strzałki odpowiadaj kolejno ci oblicze ):

3
4

;

2
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

i 1

n 1

n 2

n 1

Rozwi zanie takiego układu równa przeprowadza si w dwóch etapach. W etapie pierwszym oblicza si kolejno

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

współczynniki i :

i 1

i 2,3, ,n

i i 1

i 2,3, ,n

a w drugim na ich podstawie znajduje si rozwi zania:

i 1

i n 1,n 2, ,1

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

s obliczane jako sumy (ró nice), a nast pnie pojawiaj si w mianowniku wzoru na

kolejny współczynnik

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

mo e zosta przekształcony do postaci:

C x

gdzie

C jest pewn , nieznan jeszcze, macierz , a jest przekształconym wektorem prawych stron prawych (w

szczególno ci mo e to by niezmieniony wektor — zale nie od przekształce które zastosowano aby przej od

A do

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

C).

Wychodz c od dowolnie wybranego wektora:

(0)

mo na okre li proces iteracyjny:

(i 1)

(i )

i 0,1,2,

czyli:

(i 1)

(i )

(i 1)

(i )

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:

j 1

< 1

i 1,2, ,n

lub:

i 1

< 1

j 1,2, ,n

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:

(k)

max

j 1,2, ,n

(k)

(k 1)

lub:

(k)

j 1

(k)

(k 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:

n,n 1

n 1

w tym przypadku elementy macierzy

C okre la si nast puj co:

i j

a wektora :

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Warunki zbie no ci b d wówczas spełnione je li zachodzi:

i j

i 1,2, ,n

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 x

na lew

stron , a po prawej pozostawiaj c wszystkie (wł cznie z wyci gni t na lewo zmienn x

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 x

wykorzystuje si wszystkie nowoobliczone warto ci zmiennych: x

(k+1)

, x

(k+1)

, ..., x

i–1

(k+1)

oraz „starsze” warto ci

nast pnych zmiennych: x

(k)

, ..., x

(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:

A x

j i

(k 1)

j i 1

(k)

i 1, ,N

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:

(k 1)

i 1

j 1

(k 1)

j i 1

(k)

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

—

(k)

dostaniemy wówczas nieco inny wzór iteracyjny na

(k 1)

(k)

i 1

j 1

(k 1)

j i

(k)

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)

–1

U oraz c = (L+D)

–1

to przedstawiony algorytm mo na zapisa w takiej postaci:

(i+1)

(i)

Metoda „nadrelaskacji”

Z przedstawionym powy ej przekształceniem wi e si kolejna metoda, zwana „nadrelaksacj ” (over-relaxation) —

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

wyznaczon poprzednio poprawk powi kszamy o jaki czynnik :

(k 1)

(k)

i 1

j 1

(k 1)

j i

(k)

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:

, które mo na oczywi cie bez problemu znale innymi metodami, ale tu posłu y nam

1 ; x

2 ; x

3 ; x

jako dobra ilustracja.

Pierwszy układ ma posta :

— jak wida jest to układ z dominuj c główn

9 1

1 2

1 8 3 2
2 1

8 3

3 2

1 7

16
32

przek tn .

W drugim:

— zmieniony został ostatni wiersz — nie mamy ju dominacji

9 1

1 2

1 8

1 1

8 3

3 2

11 7

16
32

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

pierwszego wiersza:

9 1

1 12

1 8 3

2 1

8 3

3 2

1 7

56
32

Układ czwarty stanowi poł czenie zmian wprowadzonych w układach drugim i trzecim:

9 1

1 12

1 8

2 1

8 3

3 2

11 7

56
32

W układzie pi tym dodatkowo zmieniono wiersz drugi:

9 1

1 12

1 8 13

2 1

8 3

3 2

11 7

56
62

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):

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

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:

x,y

musimy przekształci go do postaci:

x,y

wybrawszy jakie warto ci x

i y

mo na posłu y si wzorami iteracyjnymi:

i 1

i 1,2,3,

Je li spełnione s warunki:

a) funkcje F

i F

maj ci głe pochodne w pewnym obszarze R,

b) przybli enia pocz tkowe (x

i y

) oraz wszystkie kolejne le w tym e obszarze R,

c) w całym obszarze R spełnione s nierówno ci:

< 1

to kolejne przybli enia zbiegaj si do rozwi za układu.

Osobnym problemem jest sposób znalezienia funkcji F

i F

. Jeden z mo liwych sposobów jest nast puj cy

(niestety — dla układu dwóch równa ):

Funkcje F

i F

maj posta :

x,y

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

gdzie współczynniki , , i stanowi przybli one rozwi zania nast puj cego układu równa :

Je li tylko pochodne cz stkowe f

i f

nie zmieniaj si zbyt szybko w otoczeniu punktu (x

) to warunki zbie no ci

b d spełnione.

Przykład. Rozwa my układ równa :

x,y

1
2

sin xy

x
2

x,y

2ex

Przepiszmy je w postaci wła ciwej dla oblicze iteracyjnych:

x,y

sin xy

x,y

2 x

1
4

2x 1

Je li wybierzemy nast puj ce warto ci pocz tkowe: x

= 0.4 i y

= 3.0 to otrzymamy nast puj cy ci g rozwi za :

0.455 0.498 0.505 0.499 0.499 0.500

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):

0.455 0.493 0.500 0.500

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:

0.454 0.493 0.500 0.500 0.500

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:

2 sin xy

1
2

2x 1

0.394 0.444 0.525 0.579 0.519 0.438 0.407 0.439 0.509

3.343 3.607 3.489 2.767 2.640 2.895 3.244 3.530 3.526

0.187 0.241 0.107 0.034 –0.181 –0.302 –0.097 –0.348 –0.159

3.343 2.497 2.044 0.688 –0.067 1.217 –0.359 0.831 0.396

0.394 –0.488 –0.790 0.062 –0.042 –0.827 –0.386 0.445

–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

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

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 m

/s) pomi dzy wszystkimi w złami (i ich kierunki).

Wiadomo, e spadek ci nienia w poziomej rurze jest dany wzorem:

1
2

gdzie f

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 :

przy czym wszystkie stałe (niezale ne od konkretnej rury) mo na zebra razem (ozn. C ):

wówczas dla ka dej pary w złów otrzymamy równanie:

gdzie:

ostatnie równanie mo na przekształci tak aby automatycznie miało prawidłowy znak:

W ka dym w 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:

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 .

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

Wzór iteracyjny b dzie mie posta :

i j

gdzie (dla skrócenia zapisu):

1
2

Dane do oblicze :

= 0.056

= 820 kg/m

= 3.8 atm

= 0

= L

= 45 m

= L

= 30 m

= D

= 3"

= D

= 4"

Metoda Newtona-Raphsona

Załó my, e dany jest nast puj cy układ równa nieliniowych:

, ,x

Aby zapisa wzory w postaci macierzowej najpierw zdefiniujemy nast puj ce macierze:

, ,x

, ,f

Zdefiniujemy:

i okre limy macierz pochodnych cz stkowych:

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 (

= (x

,...,x

)) mo emy zapisa wzór iteracyjny:

k 1

gdzie

jest rozwi zaniem nast puj cego układu równa algebraicznych:

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

było „dostatecznie” blisko pierwiastka .

Przykład. Rozpatrzymy znowu ten sam układ równa :

Metody Numeryczne w Fizyce 1

Rozwi zywanie układów równa

x,y

1
2

sin xy

x
2

x,y

2ex

Obliczymy odpowiednie pochodne:

1
2

y cos xy

x cos xy

Kolejne poprawki zmiennych x i y b dziemy oblicza rozwi zuj c układ równa :

(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 :

–0.431 –0.243 –0.258 –0.261 –0.261

1.693 1.181 0.538 0.623 0.623

czyli inn par rozwi za . Jednak przyj cie warto ci pocz tkowych x

= 0.6 i y

= 3.0 prowadzi do rozwi za

uzyskanych poprzednio:

0.502 0.500 0.500 0.500 0.500 0.500 0.500 0.500

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):