Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

1

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Rozwi zywanie równa ró niczkowych zwyczajnych

Równaniem ró niczkowym zwyczajnym n-tego rz du b dziemy nazywa równanie postaci:

F x,y, dy

dt

, d

2

y

dx

2

, d

3

y

dx

3

, , d

n

y

dx

n

0

Ogólnie równanie takie mo e by spełnione przez wiele funkcji y (x ) — aby mo liwe było jednoznaczne wyznaczenie

rozwi zania konieczne jest podanie dodatkowych informacji — warto ci funkcji y (x) i ewentualnie jej pochodnych dla

jakiej warto ci x. Koniecznych b dzie n takich warto ci. Je li podano je wszystkie dla tego samego punktu x

0

to mamy

zagadnienie pocz tkowe, je li wchodzi w gr wi cej ni jedna warto x to zagadnienie takie nazywamy brzegowym.

Zazwyczaj w przypadku równa ró niczkowych zwyczajnych wy szych rz dów mo liwe jest dokonanie podstawienia

n –1 nowych zmiennych (za wy sze pochodne) i wówczas równanie wy szego rz du zostaje zamienione na układ n równa

rz du pierwszego. Dlatego te omawia si zwykle tylko metody rozwi zywania równa rz du pierwszego.

Rozwi zania numeryczne, które b dziemy otrzymywa nie b d oczywi cie w postaci funkcji analitycznych — b d

to warto ci zbli one do warto ci rozwi zania obliczone dla pewnej sko czonej liczby punktów x

i

. Odległo ci pomi dzy

dwoma s siednimi punktami b dziemy oznacza przez h (na razie zakładamy, e s one stałe, mo liwo zmiany kroku

b dzie omawiana oddzielnie).

Zagadnienia pocz tkowe

Rozwini cie w szereg Taylora

Jedn z metod rozwi zywania równa ró niczkowych zwyczajnych postaci

jest wyra enie

dy dx

f x,y

rozwi zania y (x) w postaci szeregu Taylora w otoczeniu punktu x

0

:

y x

0

h

y x

0

h f x

0

,y x

0

h

2

2!

f x

0

,y x

0

h

3

3!

f x

0

,y x

0

Je li jako warunek pocz tkowy podano y (x

0

) to f (x

0

,y (x

0

)) otrzymamy bezpo rednio z równania ró niczkowego, a

pozostałe pochodne mo na wyliczy . Po dan ilo wyrazów nale y wstawi do badanego równania i scałkowa . Niestety

w ogólnym przypadku obliczenia kolejnych pochodnych staj si zbyt skomplikowane.

Metoda Eulera

Jest to najprostsza (w swej klasycznej postaci) z metod numerycznych. Odpowiada ona rozwi zaniu za pomoc

rozwini cia w szereg Taylora, w którym zachowano tylko pierwsze pochodne. Rozwi zywane równanie ró niczkowe

przekształcimy do postaci:

dy

dx

f x,y

i zało ymy, e podano nam pocz tkow warto funkcji y

0

(x

0

).

Działanie metod Eulera (podstawowej i jej obu modyfikacji) b dzie ilustrowane na dwóch przykładowych równaniach

ró niczkowych:

dy
dx

y

2x

y

i

dy
dx

x

y

Równania te b d rozwi ywane wszystkimi trzema metodami w przedziale <0,1>, dla warunku pocz tkowego y(0) = 1

i z trzena ró nymi krokami h wynosz cymi 0.2, 0.1 i 0.05.

Klasyczna metoda Eulera

Mo emy zastosowa nast puj cy algorytm:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

2

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

1

y x

0

h f x

0

,y x

0

y

i 1

y

i

h f x

i

,y

i

i >0

za pomoc którego wychodz c od znanych warto ci pocz tkowych x

0

,y

0

mo na obliczy kolejno warto ci szukanej funkcji

we wszystkich dalszych punktach x

i

. Mo na udowodni , e bł d (zdefiniowany jako ró nica rozwi zania analitycznego

i obliczonych warto ci przybli onych) d y do zera dla h d

cego do zera. Metoda Eulera jest metod zbie n o

dokładno ci rz du O(h) (chodzi tu o bł d globalny — tzn. po przej ciu wszystkich kroków — bł d lokalny —

odpowiadaj cy przej ciu od jednej warto ci do s siedniej — jest rz du O(h

2

)).

Zmodyfikowane metody Eulera

Pierwsza modyfikacja metody Eulera polega na obliczeniu najpierw warto ci po rednich:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

3

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

x

i

1
2

x

i

h
2

y

i

1
2

y

i

h
2

f x

i

,y

i

f

i

1
2

f x

i

1
2

,y

i

1
2

a nast pnie wykorzystaniu ich do obliczenia kolejnej warto ci y

i+1

:

y

i 1

y

i

h f

i

1
2

Druga modyfikacja metody Eulera (zwana te metod Heuna) polega na obliczeniu najpierw warto ci „wst pnie

przybli onych”:

y

i 1

y

i

h f x

i

,y

i

f

i 1

f x

i 1

,y

i 1

po czym oblicza si warto ci y

i+1

:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

4

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

i 1

y

i

h

f

i

f

i 1

2

W obu modyfikacjach odrzucona reszta rozwini cia w szereg Taylora jest rz du O(h

3

) (jest to bł d lokalny — tzn.

odpowiadaj cy przej ciu pomi dzy dwiema kolejnymi warto ciami przybli enia).

W praktyce przyjmuje si , e bł d uzyskanego rozwi zania jest równy ró nicy pomi dzy rozwa anym rozwi zaniem,

a innym rozwi zaniem tego samego zagadnienia otrzymanym dla dwa razy wi kszego kroku h.

Bł d rozwi zania mo na te szacowa na podstawie tzw.

kryterium Rungego — wynosi on ró nicy rozwa anego

rozwi zania i rozwi zania otrzymanego dla kroku dwukrotnie wi kszego.

Porównanie metod Eulera

Poni ej znajduj si wykresy przedstawiaj ce rozwi zania oraz warto ci bł dów rozwi za otrzymanych

poszczególnymi metodami dla ró nych długo ci kroków.

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

5

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Jak wida nie mo na jednoznacznie stwierdzi , która z modyfikacji metody Eulera jest lepsza — zale y to od

konkretnego rozwi zywanego zagadnienia.

Metoda Eulera z iteracj

Metod Eulera mo na udokładni stosuj c dla ka dej obliczonej warto ci y

i

dodatkowo proces iteracyjny. Traktuj c

warto obliczon klasyczn metod Eulera jako zerowe przybli enie:

y

0

i 1

y

i

h f x

i

,y

i

zastosujemy nast puj cy proces iteracyjny:

y

k

i 1

y

i

h
2

f x

i

,y

i

f x

i 1

,y

k 1

i 1

który b dzie powtarzany, a dwa kolejne wyniki b d zgodne z zadan ilo ci cyfr znacz cych, wtedy ostatni z nich

przyjmiemy za y

i+1

.

Zazwyczaj proces taki jest szybkozbie ny (tzn. zadowalaj cy wynik powini my otrzyma ju po kilku iteracjach) —

je li po 3-4 iteracjach nie uzyskujemy zbie no ci nale y zmniejszy krok h.

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

6

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Rozwi zywanie równa ruchu

Bardzo cz stym przypadkiem jest rozwi zywanie równa ruchu (dla jednego lub bardzo wielu ciał) — s to oczywi cie

równania drugiego rz du. Ka de takie równanie (drugiego rz du) mo na zapisa jako układ dwóch równa pierwszego

rz du. Poni sze równania zostały zapisane wektorowo — z punktu widzenia rozwi zania numerycznego liczba

wyst puj cych w problemie wymiarów sprowadza si do wymiaru tablic wykorzystywanych do przechowywania

rozwi zania.

Podstawow posta równania ruchu:

d

2

R

dt

2

F

m

a

Przepiszemy w nieco innej postaci:

dv

dt

a r,v

dr

dt

v

Typowy algorytm Eulera b dzie mie wówczas tak posta :

v

i 1

v

i

a

i

r

i 1

r

i

v

i

Algorytm Eulera-Cromera

Pewn modyfikacj algorytmu Eulera jest algorytm nazywany algorytmem Eulera-Cromera. Modyfikacja polega na tym,

e do obliczenia kolejnej warto ci poło enia wykorzystuje si obliczon przed chwil warto pr dko ci:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

7

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

v

i 1

v

i

a

i

r

i 1

r

i

v

i 1

Jest to bardzo niewielka zmiana (co do zasady identyczna z modyfikacj Seidla, któr poznali my w przypadku układów

równa ), jednak w wi kszo ci przypadków algorytm Eulera-Cromera zdecydowanie góruje nad typowym algorytmem

Eulera.

Algorytm Verlet’a

Załó my, e mamy do rozwi zania równanie ruchu (albo równania — dla wielu ciał). Poniewa jest to równanie 2-go

stopnia rozwi zujemy je jako układ dwóch równa rz du pierwszego. Je li interesuje nas przede wszystkim trajektoria

ciała, to mo na byłoby twierdzi , e rozwi zanie pierwszego z równa — daj ce pr dko — jest dla nas strat czasu.

W sytuacji takiej zastosowanie znajduje algorytm Verlet’a.

Jeszcze raz równanie ruchu:

d 2

R

dt

2

F

m

a

Drug pochodn zapiszemy w postaci ró nicowej:

d

2

R

dt

2

1

2

R

n 1

2R

n

R

n 1

gdzie n oznacza numer kolejny kroku czasowego, a jest jego długo ci .

Proste przekształcenie doprowadza nas do:

R

n 1

2R

n

R

n 1

2

a

n

Jak wida z wzoru na R

n+1

do „uruchomienia” tego algorytmu potrzebne s dwie kolejne warto ci pocz tkowe poło e .

W problemach fizycznych z jakimi mo emy si spotka , zazwyczaj b dziemy zna tylko jedn warto poło enia — w

chwili pocz tkowej i odpowiadaj c jej warto pr dko ci. Dla „uruchomienia algorytmu konieczne b dzie

„wygenerowanie” drugiej potrzebnej warto ci — najpro ciej mo na to zrobi zakładaj c stało przyspieszenia w

pierwszym odcinku czasu . W takim przypadku warto R

1

b dzie taka jak dla ruchu jednostajnie przyspieszonego:

R

1

R

0

V

0

2

2

a

0

Gdyby potrzebna była równie pr dko cz stki to j równie mo emy zapisa w postaci ró nicowej (skorzystamy

z wzoru centralnego):

V

dR

dt

1

2

R

n 1

R

n 1

sk d otrzymamy:

V

n

R

n 1

R

n 1

2

Jak wida wzór ten mo e w pewnych sytuacjach by „niepor czny” — obliczana pr dko jest „opó niona” w

stosunku do poło enia o jeden krok czasowy. Mo na temu zaradzi obliczaj c j wprost z działaj cej siły:

V

n 1

V

n

a

n

znacznie dokładniejsze wyniki uzyskamy wykorzystuj c do oblicze zamiast a

n

redni wyznaczon z a

n

i a

n+1

(a

n+1

b dziemy zna bo „przed chwil ” obliczyli my poło enie R

n+1

).

Poni ej przedstawiono rozwi zania prostego zagadnienia: harmonicznego ruchu wahadła bez oporów, tłumienia i

wymuszenia. Dragj ca masa jest jednostkowa, współczynnik spr ysto ci te jest jednostkowy.

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

8

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Pierwsze trzy rysunki przedstawiaj wykres poło enia drgaj cej masy w funkcji czasu. W przypadku podstawowej

metody Eulera ju „na oko” wida , e wyniki s złe — wiemy, e energia, a co za tym idzie i amplituda, opisywanego

układu powinny by stałe. Jednak wybór której z pozostałych dwóch metod mo e by trudny. Aby lepiej je oceni

wyznaczymy całkowit energi układu i zbadamy jej zmiany z czasem. Przy rzyj tych parametrach (jednostkowa masa

i współczynik spr ysto ci) całkowita energia wynikaj ca z warunków pocz tkowych równie jest jednostkowa. Fakt, e

energia ta ewidentnie nie jest stała podczas trwania ruchu, nie musi od razu dyskwalifikowa metody — ze wzgl du na

ró nego rodzaju bł dy nie mo emy oczekiwa rzeczywistej stało ci całkowitej energii — decyduj cym kryterium b dzie

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

9

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

to czy obliczona energia całkowita waha si w ograniczonym zakresie czy te amplituda jej waha wzrasta. Oczywi cie

celem wyboru metody nale ałoby sporz dzi te wykresy dla du o dłu szego czasu — poni ej czas przedłu ono

dziesi ciokrotnie:

Metody Rungego-Kutty

Metod Eulera otrzymujemy dzi ki rozwini ciu szukanej funkcji w szereg Taylora i odrzuceniu wyrazów

zawieraj cych pochodne rz du wy szego ni pierwszy — bł d lokalny jest wówczas rz du O(h

2

). Zachowanie wy szych

pochodnych jest mo liwe tylko dla bardzo szczególnych postaci równa , w ogólnym przypadku b d one zbyt

skomplikowane. Mo liwe jest jednak stworzenie metody jednokrokowej (tzn. wymagaj cej znajomo ci tylko jednego

wcze niejszego punktu) o dokładno ci lokalnej wy szego rz du. S to metody Rungego-Kutty.

Ogólnie algorytm metod Rungego-Kutty jest nast puj cy:

y

i 1

y

i

h

x

i

,y

i

,h

funkcja (x,y,h) jest odpowiednio dobran aproksymacj funkcji f (x,y ) na przedziale (x

i

,x

i+1

).

Metoda Rungego-Kutty rz du drugiego

Załó my, e funkcj (x,y,h ) zdefiniujemy jako redni z dwóch przybli e pochodnej f (x,y ), oznaczymy je przez

k

1

i k

2

: (x,y,h) = ak

1

+ bk

2

. Wówczas algorytm Rungego-Kutty (drugiego rz du) b dzie mie posta :

y

i 1

y

i

h ak

1

bk

2

Załó my dalej, e:

k

1

f x

i

,y

i

k

2

f x

i

ph,y

i

qhf x

i

,y

i

f x

i

ph,y

i

qhk

1

gdzie p i q s stałymi, które wyznaczymy pó niej.

Rozwi my k

2

w szereg Taylora opuszczaj c wszystkie wyrazy zawieraj ce pot gi h wy sze ni pierwsza:

k

2

f x

i

ph,y

i

qhf x

i

,y

i

f x

i

,y

i

phf

x

x

i

,y

i

qhf x

i

,y

i

f

y

x

i

,y

i

O h

2

Szereg Taylora dla funkcji dwóch zmiennych:

f (x r,y s)

f (x,y)

rf

x

sf

y

r

2

2

f

xx

rsf

xy

s

2

2

f

yy

Podstawiaj c definicj k

1

i rozwini cie k

2

do algorytmu otrzymamy:

y

i 1

y

i

h af x

i

,y

i

bf x

i

,y

i

h

2

bpf

x

x

i

,y

i

bqf x

i

,y

i

f

y

x

i

,y

i

O h

3

Nast pnie rozwiniemy w szereg Taylora funkcj b d c rozwi zaniem dokładnym:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

10

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y x

i

h

y x

i 1

y x

i

hf x

i

,y x

i

h

2

2

f x

i

,y x

i

h

3

3!

f

,y

gdzie jest liczb z przedziału (x

i

,x

i+1

). Pochodna f (x

i

,y(x

i

)):

f x

i

,y x

i

f

x

x

i

,y x

i

f

y

x

i

,y x

i

f x

i

,y x

i

Je li teraz za damy, aby współczynniki przy odpowiednich pot gach h były równe to otrzymamy:

y x

i

y

i

f x

i

,y x

i

a b f x

i

,y

i

1
2

f

x

x

i

,y x

i

f

y

x

i

,y x

i

f x

i

,y x

i

bpf

x

x

i

,y

i

bqf

y

x

i

,y

i

f x

i

,y

i

Je li zało ymy y

i

= y (x

i

) i za damy, aby dla wszystkich funkcji współczynniki przy h

2

były równe otrzymamy

a + b = 1, bp = ½ oraz bq = ½, a zatem:

a

1

b

p

q

1

2b

Poniewa mamy tylko trzy równania, a cztery niewiadome, mo emy jedn z nich wybra dowolnie, najcz ciej

przyjmuje si b = ½ albo b = 1.

W przypadku b = ½ otrzymamy a = ½, p = 1 i q = 1. Algorytm b dzie mie wówczas posta :

y

i 1

y

i

h
2

f x

i

,y

i

f x

i

h,y

i

hf x

i

,y

i

co mo na zapisa równie w postaci:

y

i 1

y

i

hf x

i

,y

i

y

i 1

y

i

h
2

f x

i

,y

i

f x

i 1

,y

i 1

Jest to metoda identyczna z podan wcze niej drug modyfikacj metody Eulera. Metody tego typu nazywane s

metodami predyktor-korektor, poniewa najpierw za pomoc metody daj cej gorsze przybli enie przewidujemy warto

y

i+1

(jest to warto jak y

i+1

miałoby gdyby nachylenie (pochodna) wynosiło k

1

), a nast pnie uzyskujemy skorygowan

warto y

i+1

(jako redni wa on z warto ci f na obu ko cach przedziału (x

i

,x

i+1

). Metoda ta jest znana równie pod nazw

metody Heuna. Ilustruje j lewy rysunek powy ej.

Dla b = 1 otrzymamy a = 0, p = ½ i q = ½ i algorytm przyjmie posta :

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

11

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

i

1
2

y

i

h
2

f x

i

,y

i

y

i 1

y

i

hf x

i

h
2

,y

i

1
2

Ten algorytm z kolei odpowiada przedstawionej ju pierwszej modyfikacji metody Eulera. Najpierw za pomoc

normalnej metody Eulera otrzymujemy warto y w połowie przedziału (x

i

,x

i+1

), a nast pnie znajdujemy warto f w tym

punkcie i wykorzystujemy j do obliczenia y na ko cu przedziału. Ilustruje to prawy rysunek powy ej.

Metoda Rungego-Kutty rz du trzeciego

Metody Rungego-Kutty wy szych rz dów wyprowadza si analogicznie. Dla metody rz du trzeciego funkcja ma

posta : = ak

1

+ bk

2

+ ck

3

, gdzie k

1

, k

2

i k

3

s przybli eniami warto ci pochodnej w ró nych punktach wewn trz przedziału

(x

i

,x

i+1

):

k

1

f x

i

,y

i

k

2

f x

i

ph,y

i

phk

1

k

3

f x

i

rh,y

i

shk

2

r s hk

1

Algorytm jest wi c dany wzorem:

y

i 1

y

i

h ak

1

bk

2

ck

3

Aby wyznaczy stałe a, b, c, p, r i s rozwijamy k

2

i k

3

oraz funkcj y w szeregi Taylora wokół punktu (x

i

,y

i

) odrzucaj c

wyrazy zawieraj ce pot gi h wy sze ni trzecia i przyrównujemy odpowiednie wyrazy. Otrzymamy układ równa :

a

b

c

1

bp

cr

1
2

bp

2

cr

2

1
3

cps

1
6

Dwie z tych stałych mo na przyj dowolnie. Kutta dokonał wyboru, przy którym algorytm dla metody trzeciego rz du

jest nast puj cy:

y

i 1

y

i

h
6

k

1

4k

2

k

3

k

1

f x

i

,y

i

k

2

f x

i

1
2

h,y

i

1
2

hk

1

k

3

f x

i

h,y

i

2hk

2

hk

1

W powy szych wzorach przyj to a = c = 1/6, b = 4/6, p = ½, r = 1 i s = 2.

Poni ej przedstawiam bł dy metody Runge-Kutty trzeciego rz du zastosowanej do identycznych jak poprzednio dwóch

równa przykładowych:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

12

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Metoda Rungego-Kutty rz du czwartego

Algorytm metod rz du czwartego jest oczywi cie nast puj cy:

y

i 1

y

i

h ak

1

bk

2

ck

3

dk

4

U ywanych jest kilka algorytmów, autorem dwu poni szych jest Kutta:

y

i 1

y

i

h
6

k

1

2k

2

2k

3

k

4

k

1

f x

i

,y

i

k

2

f x

i

1
2

h,y

i

1
2

hk

1

k

3

f x

i

1
2

h,y

i

1
2

hk

2

k

4

f x

i

h,y

i

hk

3

albo

y

i 1

y

i

h
8

k

1

3k

2

3k

3

k

4

k

1

f x

i

,y

i

k

2

f x

i

1
3

h,y

i

1
3

hk

1

k

3

f x

i

2
3

h,y

i

1
3

hk

1

hk

2

k

4

f x

i

h,y

i

hk

1

hk

2

hk

3

Jednak najcz ciej wykorzystywan metod czwartego rz du jest metoda Gilla:

y

i 1

y

i

h
6

k

1

2 1

1

2

k

2

2 1

1

2

k

3

k

4

k

1

f x

i

,y

i

k

2

f x

i

1
2

h,y

i

1
2

hk

1

k

3

f x

i

1
2

h,y

i

1
2

1

2

hk

1

1

1

2

hk

2

k

4

f x

i

, h,y

i

1

2

hk

2

1

1

2

hk

3

Poni ej przedstawiam bł dy metody Runge-Kutty czwartego rz du zastosowanej do identycznych jak poprzednio

dwóch równa :

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

13

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

k

1

f x,y

k

2

f x c

2

h, y ha

2,1

k

1

k

s

f x c

s

h, y h

s 1

j 1

a

s,j

k

j

x,y,h

s

i 1

b

i

k

i

y

k 1

y

k

h

x

k

,y

k

,h

Ogólny schemat metod Runge-Kutta

Ogólny wzór podany przez Rungego:

Odnosi si on w gruncie rzeczy do wszystkich metod Runge-Kutta. Indeks s odpowiada tu rz dowi metody. W

literaturze mo na spotka taki zapis zmiennych wyst puj cych w tym schmacie:

0

c

2

a

2,1

c

3

a

3,1

a

3,2

c

s

a

s,1

a

s,2

a

s,s 1

b

1

b

2

b

s 1

b

s

Ogólnie rzecz bior c podczas „własnor cznego” tworzenia jakiego schematu tego typu potrzebne b d rozwini cia

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

14

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

f x

1

,x

2

, ,x

n

f a

1

,a

2

, ,a

n

n

i 1

f

x

i a

1

,a

2

, ,a

n

x

i

a

i

1

2!

n

i 1

n

j 1

2

f

x

i

x

j a

1

,a

2

, ,a

n

x

i

a

i

x

j

a

j

R

m

x

1

,x

2

, ,x

n

Taylora funkcji kilku zmiennych. Ogólny wzór na takie rozwini cie funkcji f w pobli u punktu a jest nast puj cy:

Bł dy metod Rungego-Kutty

Ogólnie bł d w metodach Rungego-Kutty mo na okre li nast puj co. Załó my, e bł d zdefiniujemy jako ró nic

pomi dzy warto ci rozwi zania analitycznego y (gdyby było znane), a warto ci otrzyman przy kroku h

1

lub h

2

w

wyniku scałkowania równania na przedziale (x

i

,x

i+1

):

1

y

y

1

m

h

m 1

1

x

i 1

x

i

h

1

2

y

y

2

m

h

m 1

2

x

i 1

x

i

h

2

Dziel c oba powy sze równania stronami otrzymamy:

y

y

1

y

2

h

1

h

2

m

1

h

1

h

2

m

co dla h

2

= h

1

/2 przechodzi w:

y

y

1

2

m

y

2

1 2

m

Podstawiaj c t warto do pierwszego lub drugiego wyra enia na bł d otrzymamy:

1

2

m

y

2

y

1

2

m

1

2

y

2

y

1

2

m

1

gdzie m jest rz dem metody, a y

1

i y

2

s warto ciami otrzymanymi odpowiednio dla kroków h

1

i h

2

= h

1

/2. Oczywi cie taki

sposób obliczania bł du wymaga praktycznie potrojenia ilo ci oblicze , nie nadaje si wi c do praktycznego zastosowania.

W przypadku metod nieparzystego rz du istnieje mo liwo przeprowadzenia po ka dym cyklu cyklu odwrotnego (tzn.

ponownego obliczenia warto ci y

i

na podstawie wła nie obliczonej warto ci y

i+1

). Po takim obliczeniu jako bł d mo na

przyj połow ró nicy mi dzy obliczonymi w obu kierunkach warto ciami y

i

.

Wybór długo ci kroku

Collatz podał praktyczne jako ciowe kryterium doboru kroku całkowania dla pierwszego z podanych algorytmów

Kutty czwartego rz du: po ka dym kroku nale y obliczy warto wyra enia: k

3

– k

2

/ k

2

– k

1

. Je eli przekracza ona

kilka setnych to długo kroku nale y zmniejszy .

Cz sto stosowana jest metoda polegaj ca na obliczaniu (co jak liczb kroków) warto ci przybli onej bł du:

= y

1

– y

2

/ 15. Je li warto jest mniejsza od po danej dokładno ci obliczenia prowadzone s dalej. Je li jest

mniejsza ni 1/50 zało onej dokładno ci krok powi ksza si dwukrotnie. Je li jest wi ksza — obliczenia powtarza si (od

miejsca gdzie ostatnio sprawdzano ) z krokiem dwukrotnie mniejszym.

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

15

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Najcz ciej spotykan metod oblicze ze zmiennym krokiem jest metoda Runge-Kutta-Fehlberg (o której dalej).

W przypadku programu maj cego mo no zmiany długo ci kroku musimy oczywi cie poda dopuszczalny zakres

tych zmian, a wi c poda jego maksymaln i minimaln

dopuszczaln długo . Oczywi cie wektor wynikowy takiego

programu (zawieraj cy obliczone warto ci y

i

) ma nieprzewidywaln długo . Co wi cej, konieczne b dzie zapisanie w

pami ci równie wektora „chwil czasu”, którym odpowiadaj wyznaczone warto ci funkcji — zatem cały program mo e

okaza si bardziej „pami cio erny”. Z drugiej jednak strony, je li charakter funkcji nie wymusi skracania kroku, a

umo liwi jego wydłu anie, mo emy zaoszcz dzi i na pami ci i na czasie (a przynajmniej na cyklach oblicze ). Je li oka e

si , e uzyskane rozwi zania nie wypadaj w punktach, na których nam zale ało b dziemy musieli przeprowadzi

interpolacj .

Stabilno metod jednokrokowych

Przykład: Rozwa my znowu równanie:

dy
dx

f x,y

x

y

którego rozwi zanie analityczne ma posta :

y x

x 1 1 x

0

y x

0

e

x x

0

Przy warunku pocz tkowym y (0) = –1 jego rozwi zaniem analitycznym jest: y (x ) = –x – 1. Jak wida przy takim

warunku pocz tkowym składnik eksponencjalny w rozwi zaniu analitycznym znika, ale nawet minimalna zmiana warunku

pocz tkowego, np. na –0,99999 spowoduje w ko cu znaczn zmian warto ci y dla du ych x. Tymczasem rozwi zywanie

równa ró niczkowych metodami jednokrokowymi sprowadza si w ka dym kroku do rozwi zywania kolejnego

zagadnienia pocz tkowego — wynik uzyskany w kroku i staje si warunkiem pocz tkowym dla obliczenia wyniku i +1.

Tego rodzaju niestabilno jest cech konkretnego równania i jego warunku pocz tkowego, a nie zale y od tego jaki

konkretny algorytm jednokrokowy zostanie wykorzystany.

Metoda Runge-Kuta-Fehlberg

Metoda ta jest znana równie jako „RKF45" (funkcja „ODE45" Matlaba wykorzystuje t wła nie metod ). Jest to

metoda o zmiennym kroku, jednak w celu doboru długo ci kroku nie stosuje si w niej porównywania wyników

otrzymanych przy dwóch ró nych długo ciach kroków — zamiast tego porównuje wyniki trzymane dla dwóch ró nych

rz dów metody obliczeniowej (czwartego i pi tego).

Algorytm.

Kolejne współczynniki:

k

1

hf t

i

,y

i

k

2

hf t

i

1
4

h,y

i

1
4

k

1

k

3

hf t

i

3
8

h,y

i

3

32

k

1

9

32

k

2

k

4

hf t

i

12
13

h,y

i

1932
2197

k

1

7200
2197

k

2

7296
2197

k

3

k

5

hf t

i

h,y

i

439
216

k

1

8k

2

3680

513

k

3

845

4104

k

4

k

6

hf t

i

1
2

h,y

i

8

27

k

1

2k

2

3544
2565

k

3

1859
4104

k

4

11
40

k

5

Kolejne przybli enie (obliczenia rz du 4):

y

i 1

y

i

25

216

k

1

1408
2565

k

3

2197
4104

k

4

1
5

k

5

Lepsze rozwi zanie (obliczenia rz du 5):

z

i 1

y

i

16

135

k

1

6656

12825

k

3

28561
56430

k

4

9

50

k

5

2

55

k

6

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

16

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Nast pny krok wyznacza si badaj c nast puj c wieko :

s

Tol h

2 z

i 1

y

i 1

1
4

0.84

Tol h

z

i 1

y

i 1

1
4

Długo kroku mo emy regulowa na kilka sposobów:

• metod „konserwatywn ”: Je li s < 3/4 — krok zmniejszamy o połow , je li s > 3/2 — krok zwi kszamy dwukrotnie

— jak łatwo zauwa y przy zastosowaniu tej metody przez wi kszo czasu krok pozostanie stały

• metod „dynamiczn ”: Sprowadzamy warto s do wybranego zakresu (np. <0.1,4.0> — tzn. je li s < 0.1 to

podstawiamy s = 0.1 i je li s > 4 to podstawiamy s = 4), a nast pnie mno ymy dotychczasowy krok przez s (post pujemy

tak w ka dym kroku).

• metod „mieszan ”: Post pujemy tak jak w metodzie „konserwatywnej”, ale zmian kroku uzyskujemy mno c jego

aktualn długo przez obliczon warto s.

Oczywi cie niezale nie od zastosowanej metody w programie powinny by narzucone granice dopuszczalnych zmian

długo ci kroku — tzn. maksymalna i minimalna nieprzekraczalna długo kroku.

Podane warto ci granic przedziałów warto ci współczynnika s wynikaj z praktyki, mo emy je dobra do konkretnego

rozwi zywanego zagadnienia

Istnieje pogl d, e jako wynik procedury nale y, zamiast wektora y, przedstawia obliczany równolegle wektor z.

Poni sze wykresy przedstawiaj porównanie metody Runge-Kutta rz du 4 i metody Runge-Kutta-Fehlebrg

(wykorzystuj c przedstawion wy ej „konserwatywn ” metod zmiany kroku) zastosowanych do ju wcze niej podanych

dwóch równa przykładowych. Pocz tkowy krok zastosowany w metodzie RKF był taki sam jak krok u ywany przez cały

czas oblicze w metodzie RK4. Warto zauwa y , e zastosowanie w metodzie RKF zbyt du ej tolerancji mo e prowadzi

do znacznych bł dów.

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

17

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Metoda Runge-Kutta-Nyström

Jak wcze niej powiedziano liniowe równania wy szych rz dów mo emy zast powa układami równa ni szych

rz dów i rozwi zywa je kolejno. Jednak w przypadku równa rz du drugiego opracowano metod czwartego rz du,

b d c adaptacj metody Runge-Kutta. Metoda ta stosuje si do równa postaci:

y

f x,y,y

Zakładamy, e znane s warto ci funkcji i jej pochodnej w punkcie pocz tkowym. Dla równoodległych punktów x

i

(krok

równy h) algorytm jest nast puj cy:

y

n 1

y

n

h y

n

1
3

k

1

k

2

k

3

k

1

1
2

hf x

n

,y

n

,y

n

k

2

1
2

hf x

n

1
2

h , y

n

1
2

h y

n

1
2

k

1

, y

n

k

1

k

3

1
2

hf x

n

1
2

h , y

n

1
2

h y

n

1
2

k

1

, y

n

k

2

k

4

1
2

hf x

n

h , y

n

h y

n

k

3

, y

n

2k

3

y

n 1

y

n

1
3

k

1

2k

2

2k

3

k

4

W metodzie tej mo na stosowa opisany ju wcze niej dla metody Runge-Kutta sposób zmieniania długo ci kroku

h w trakcie oblicze , z tym, e oprócz ró nicy warto ci y

1

i y

2

nale y kontrolowa równie ró nice pomi dzy

odpowiednimi pochodnymi (i tak przecie obliczanymi); wydłu enie kroku jest mo liwe gdy zarówno warto ci funkcji

jak i pochodnych spełniaj warunek.

Zastosowanie metod jednokrokowych do układów równa

Jak ju wcze niej wspominałem, podr czniki nie omawiaj sposobów rozwi zywania równa ró niczkowych rz dów

wy szych ni pierwszy, poniewa zawsze mo na je sprowadzi do układu odpowiedniej liczby równa rz du pierwszego.

Rozwa my równanie ruchu (czyli równanie rz du drugiego):

2

r

t

2

a

podstawiaj c now zmienn (pr dko ) otrzymamy układ dwóch równa (wektorowych) rz du pierwszego (albo sze ciu

równa skalarnych):

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

18

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

r

t

v

v

t

a

x

t

v

x

y

t

v

y

z

t

v

z

v

x

t

a

x

v

y

t

a

y

v

z

t

a

z

W dalszym ci gu b d omawia skalarne równanie drugiego rz du, które po przekształceniu do układu dwóch równa

rz du pierwszego zapiszemy wektorowo:

d

2

r

dt

2

a

dr

dt

v

dv

dt

a

d

dt

r
v

v

a

Gdyby my rozwa ali ruch w trzech wymiarach we wzorach powy ej mieliby my wektory nie dwu- ale sze cioelementowe.

Aby zastosowa metod Runge-Kutty trzeciego rz du wła ciwe dla niej wzory zapiszemy w postaci wektorowej:

k

1

f t

i

, y

i

f t

i

,

r

i

v

i

v

i

a(r

i

)

k

2

f t

i

1
2

h, y

i

1
2

h k

1

f t

i

1
2

h,

r

i

v

i

1
2

h

v

i

a(r

i

)

f t

i

1
2

h,

r

i

1
2

hv

i

v

i

1
2

ha(r

i

)

v

i

1
2

ha(r

i

)

a(r

i

1
2

hv

i

)

k

3

f t

i

h, y

i

2h k

2

h k

1

f t

i

h,

r

i

v

i

2h

v

i

1
2

ha(r

i

)

a(r

i

1
2

hv

i

h

v

i

a(r

i

)

f t

i

h,

r

i

2hv

i

h

2

a(r

i

) hv

i

v

i

2ha(r

i

1
2

hv

i

) ha(r

i

)

v

i

2ha(r

i

1
2

hv

i

) ha(r

i

)

a(r

i

2hv

i

h

2

a(r

i

) hv

i

I ostateczne rozwi zanie:

y

i 1

y

i

h
6

k

1

4 k

2

k

3

Warto zwróci uwag , e do zapisania tych wzorów w Scilabie wyniki przekształce widoczne po prawej stronie

przedstawionych powy ej wzorów nie b d potrzebne — funkcj f zdefiniujemy jako przyjmuj c dwa argumenty, z

których pierwszym b dzie skalar (czas), a drugim wektor zawieraj cy wszystkie niewiadome.

Analogicznie b dziemy post powa w przypadku ka dej innej metody.

Metody wielopunktowe

Wszystkie metody przedstawione powy ej s nazywane jednopunktowymi — a to dlatego, e do prowadzenia oblicze

wystarczy znajomo tylko jednej warto ci szukanej funkcji. Podanie warunku pocz tkowego jest wystarczaj ce. W

przeciwie stwie do nich metody wielopunktowe wymagaj znajomo ci warto ci funkcji w kilku punktach pocz tkowych.

Zazwyczaj jednak warunek pocz tkowy zadany jest w jednym punkcie, wówczas mo liwe jest rozpocz cie oblicze któr

z przedstawionych powy ej metod jednopunktowych, obliczenie niezb dnej liczby warto ci pocz tkowych, a nast pnie

kontynuowanie oblicze metod wielopunktow .

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

19

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Metoda Adamsa

Załó my, e dla zagadnienia y = f(x,y) z warunkiem pocz tkowym y(x

0

) = y

0

znamy trzy kolejne warto ci

y

1

= y(x

1

) = y(x

0

+h), y

2

= y(x

2

) = y(x

0

+2h), y

3

= y(x

3

) = y(x

0

+3h) — obliczone w inny sposób.

Rozwi zanie znajdujemy w dwóch etapach: w pierwszym dokonujemy ekstrapolacji znanych punktów i znajdujemy

przybli on warto punktu kolejnego:

y

i 1

y

i

h

24

55f

i

59f

i 1

37f

i 2

9f

i 3

A nast pnie, wykorzystuj c obliczon warto y

i+

*

1

do obliczenia f

i+

*

1

w etapie drugim dokonujemy jej udokładnienia:

y

i 1

y

i

h

24

9f

i 1

19f

i

5f

i 1

f

i 2

Poni ej przedstawiam bł dy metody Adamsa zastosowanej do identycznych jak poprzednio dwóch równa , do

obliczenia punktów pocz tkowych wykorzystano metod Runge-Kutty drugiego, trzeciego i czwartego rz du:

Metoda Milne'a

Metoda ta działa analogicznie jak przedstawiona powy ej metoda Adamsa — obliczenia przeprowadzamy równie

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

20

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

w dwóch etapach:

y

i 1

y

i 3

4h

3

2f

i 2

f

i 1

2f

i

y

i 1

y

i 1

h
3

f

i 1

4f

i

f

i 1

Jest to metoda 4-go rz du. Przyjmuje si , e bł d bezwzgl dny warto ci udokładnionej (obliczonej w drugim etapie)

jest równy w przybli eniu 1/29 ró nicy pomi dzy warto ci uzyskan w pierwszym etapie a warto ci z etapu drugiego

— mo liwe jest wi c wygodne kontrolowanie dokładno ci podczas oblicze .

Poni ej przedstawiam bł dy metody Milne’a czwartego rz du zastosowanej do identycznych jak poprzednio dwóch

równa , do obliczenia punktów pocz tkowych wykorzystano metod Runge-Kutty drugiego, trzeciego i czwartego rz du:

Istnieje tak e metoda Milne'a rz du 6-go:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

21

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

i 1

y

i 5

3h

10

11f

i

14f

i 1

26f

i 2

14f

i 3

11f

i 4

y

i 1

y

i 3

2h

45

7f

i 1

32f

i

12f

i 1

32f

i 2

7f

i 3

Poni ej przedstawiam bł dy metody Milne’a szóstego rz du zastosowanej do identycznych jak poprzednio dwóch

równa , do obliczenia punktów pocz tkowych wykorzystano metod Runge-Kutty drugiego, trzeciego i czwartego rz du:

Dokładno rozwi zania a dokładno punktów pocz tkowych

Jak wida na powy szych wykresach, w przypadku metody Milne’a (i ogólnie w przypadku wszystkich metod

wielokrokowych) obliczenie warto ci pocz tkowych ze zbyt mał dokładno ci mo e uczyni całe rozwi zanie

bezu ytecznym.

Metody "półanalityczne"

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

22

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Rozwa my zagadnienie Cauchy-ego dla równania ró niczkowego z warunkiem pocz tkowym:

y

f x,y

;

y x

0

y

0

Metoda kolejnych przybli e polega na znalezieniu y(x) jako granicy ci gu funkcji y

n

(x) znajdowanych za pomoc

wzoru rekurencyjnego:

y

n

x

y

0

x

x

0

f x ,y

n 1

x dx

Udowodniono, e je li w pewnym obszarze R (okre lonym jako: x – x

0

a oraz y – y

0

b) spełniony jest warunek:

f x,y

1

f x,y

2

N y

1

y

2

gdzie N jest jak stał , to niezale nie od wyboru funkcji pocz tkowej y

0

(x) kolejne przybli enia y

n

(x) w pewnym

przedziale (x,x+h) d

do poszukiwanego rozwi zania zagadnienia.

Bł d n-tego przybli enia w tym przedziale mo na oszacowa nast puj co:

n

y x

y

n

x

MN

n

x x

0

n 1

n 1 !

gdzie M jest okre lone jako:

M

max

x,y R

f x,y

a h jest okre lone jako:

h

min a, b

M

Metod t mo na równie łatwo zastosowa do układów równa ró niczkowych.

Przykład: Znale kilka kolejnych przybli e rozwi zania równania: y = x

2

+ y

2

dla warunku pocz tkowego y(0) = 0.

Zapiszmy równanie w postaci całkowej:

y x

x

0

x

2

y

2

dx

Jako przybli enie wyj ciowe przyjmiemy y

0

(x) 0. Kilka kolejnych przybli e b dzie mie nast puj c posta :

y

0

0

y

1

x

3

3

y

2

x

3

3

x

7

63

y

3

x

3

3

x

7

63

2x

11

2079

x

15

59535

y

4

x

3

3

x

7

63

2x

11

2079

13x

15

218295

82x

19

37328445

662x

23

10438212015

4x

27

3341878155

x

31

109876902975

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

23

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

5

x

3

3

x

7

63

2x

11

2079

13x

15

218295

46x

19

12442815

7382x

23

39665205657

428x

27

50998121559

17853193x

31

52289648987481675

738067298x

35

60394544580541334625

10307579354x

39

26851414520508677374275

6813116x

43

667791107219128471425

89797289962x

47

407191917287292044684671125

1159084x

51

308193360994502430591375

721012x

55

15328564533673936679413125

8x

59

21664518085681213656375

x

63

760594829864786522589375

Okre limy bł d przybli enia numer 3. Przyjmijmy a = 1 i b = 0,5 (mog to by dowolne liczby jako e funkcja x

2

+ y

2

jest okre lona i ci gła na całej płaszczy nie xy). Otrzymamy wówczas:

M

max f x,y

max x

2

y

2

1,25

N

max f x,y

max 2y

1

sk d wynika, e h = 0,4. Dla przedziału (0;0,4) otrzymamy:

y x y

3

x

1,25 1

3

x

4

4!

5

96

x

4

i dla całego przedziału:

m a x

(0;0,4)

y x y

3

x

5 0,4

3

96

0,00133

Obliczone w powy szy sposób bł dy poszczególnych przybli e s nast puj ce: n=2 — 0.1, n=3 — 0.01(3), n=4

— 0.001(3), n=5 — 0.00010(6). Jak wida dla ka dego kolejnego przybli enia bł d maleje o rz d wielko ci.

Zastosowany wzór na bł d cz sto daje wynik zbyt ostro ny (tzn. obliczony bł d jest zbyt du y). W praktyce obliczenia

prowadzi si do chwili, gdy kolejne y

n

i y

n–1

zgadzaj si z wymagan dokładno ci .

Na rysunku przedstawiaj cym ró nice kolejnych przybli e znacznie powi kszono skal osi y. O x obejmuje nieco

zbyt du y zakres — do oblicze przyj to zakres x-ów od –1 do +1. Je li chcemy aby rozwi zanie było dokładne dla

wi kszego zakresu x musimy posłu y si przybli eniem wy szego rz du. W sytuacji dost pno ci programów całkuj cych

metoda ta mo e by bardzo u yteczna.

Przykład: Zamiast przyj jako funkcj wyj ciow funkcj to samo ciowo równ zeru (jak w poprzednim przykładzie)

mo na znale funkcj bardziej zbli on do poszukiwanego rozwi zania. Na przykład dla równania y = x + 0,1y

2

z

warunkiem pocz tkowym y(0) = 1 funkcj wyj ciow mo e by :

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

24

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

0

x

y

0

y 0

1!

x

y 0

2!

x

2

Powy szy wzór uzyskali my rozwijaj c y(x) wokół punktu x

0

=0 (wówczas x jest odległo ci od x

0

) w szereg Taylora

ograniczony do drugich pochodnych. Z warunków pocz tkowych i z równania otrzymamy wówczas:

y

0

1

y 0

0,1y

2

0

0,1

y 0

1

0,2y

0

y

0

1,02

zatem ostatecznie b dzie to funkcja:

y

0

x

1

0,1x

0,51x

2

Zagadnienia brzegowe

W przypadku zagadnie pocz tkowych dysponowali my pełn informacj o szukanej funkcji w jednym wybranym

punkcie pocz tkowym (tzn. znali my warto funkcji i warto ci odpowiedniej liczby jej pochodnych — zale nie od

stopnia rozwi zywanego równania). W przypadku zagadnie brzegowych tak informacj nie dysponujemy, jest ona

„rozło ona” pomi dzy dwa punkty ograniczaj ce interesuj cy nas obszar, w którym chcemy znale rozwi zanie.

Wprawdzie rozwi zywanie zagadnie pocz tkowych jest w pewnym sensie łatwiejsze (mo liwe jest zastosowanie du o

prostszych metod) to jednak kolejne obliczane warto ci funkcji s obarczone coraz to wi kszymi bł dami, w kra cowych

przypadkach mo e to powodowa konieczno porzucenia prostych metod rozwi zywania. W przypadku zagadnie

brzegowych ograniczenia narzucone na oba kra ce przedziału powoduj bardziej równomierne rozło enie bł dów

znalezionych rozwi za .

Metoda strzałów

Metoda ta jest logiczn konsekwencj ch ci zastosowania do zagadnie brzegowych metod opracowanych dla

zagadnie pocz tkowych (a wi c w zasadzie nie nadaj cych si do takiego zastosowania). Czasami nazywana jest metod

prób i bł dów.

Załó my, e interesuje nas zagadnienie brzegowe dla równania drugiego rz du. W przypadku zagadnienia

pocz tkowego oczekiwaliby my warto ci funkcji i jej pochodnej w punkcie pocz tkowym, natomiast w przypadku

zagadnienia brzegowego b dziemy dysponowa warto ciami funkcji na obu kra cach przedziału (oznaczmy je

odpowiednio przez a i b oraz przez f(a) i f(b)) natomiast nie b dzie dana adna warto pochodnej szukanej funkcji..

Je li odgadniemy (oboj tne czy bior c warto z sufitu czy

te szacuj c na podstawie jakich sensownych przesłanek)

warto pochodnej na jednym z kra ców (oznaczmy j przez

c

1

), to b dziemy mogli zastosowa któr z wcze niej

opisanych metod i za jej pomoc obliczy wszystkie warto ci

szukanej funkcji a do drugiego brzegu wł cznie (warto t

oznaczymy przez f

1

). Nie oczekujemy oczywi cie, e obliczona

w ten sposób warto funkcji w punkcie b b dzie zgadza si

z wynikaj c z warunków zadania warto ci brzegow . Jednak

uzyskamy w ten sposób par warto ci: zało ona warto

pochodnej i otrzymana warto funkcji na drugim brzegu.

Teraz mo emy przyj inn warto pochodnej na brzegu (c

2

)

i znów przeprowadzi obliczenia (otrzymuj c f

2

).

Je li b dziemy mie szcz cie jedna z otrzymanych

warto ci b dzie za du a a druga za mała. Wtedy warto c

3

wyznaczymy cho by za pomoc interpolacji liniowej:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

25

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

c

3

c

1

c

2

c

1

f

1

f b

f

1

f

2

(przy zało eniu, e f

1

było za du e a f

2

za małe).

Pierwszy przykład nie jest wła ciwie zagadnieniem brzegowym — rozwi zywane jest równanie pierwszego rz du,

jednak podano warunek (warto funkcji szukanej) nie dla pocz tku lecz dla ko ca przedziału. Mo na byłoby oczywi cie

rozwi za to zagadnienie któr kolwiek z poznanych ju metod stosuj c krok o ujemnej warto ci jednak i w takim

przypadku mo na zastosowa metod strzałów.

Kolejne przykłady dotycz równania drugiego rz du z podanymi warto ciami funkcji szukanej na obu kra cach

przedziału (jest to ju wi c zagadnienie brzegowe). Jak wida całkiem niezł dokładno rozwi zania mo na uzyska

nawet bez nadmiernej dbało ci o dokładno poszukiwanej warto ci pocz tkowej pierwszej pochodnej.

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

26

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Metody ró nicowe

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

27

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

Metody ró nicowe w odniesieniu do zagadnie brzegowych dla równa ró niczkowych zwyczajnych s takie same

jak w przypadku równa ró niczkowych cz stkowych, rozwi zania s jednak prostsze, bo mamy do czynienia tylko z

jedn zmienn niezale n . Zastosowane mog by praktycznie wszystkie wzory podane dla równa ró niczkowych

cz stkowych.

Metoda passage

Dla równania liniowego drugiego rz du mo liwe jest wyprowadzenie wzorów z uwzgl dnieniem dowolnych warunków

brzegowych (trzeciego rodzaju — zawieraj cych i warto funkcji i warto pochodnej na brzegu przedziału). Załó my,

e rozwi zywany problem ma posta :

y

p(x)y

q(x)y

f(x)

0

y(a)

1

y (a)

A

0

y(b)

1

y (b)

B

gdzie p(x), q(x) i f(x) s znanymi funkcjami, a i b s warto ciami granicznymi zmiennej x za

0

,

1

,

0

,

1

, A i B s

znanymi stałymi. Zakładamy dalej, e odcinek a–b został podzielony na N kroków o długo ci x, przy czym w zeł nr 0

odpowiada lewej granicy, a w zeł nr N — granicy prawej.

Wówczas oznaczaj c:

m

i

2q

i

x

2

4

2 p

i

x

;

k

i

2 p

i

x

2 p

i

x

mo emy przeprowadzi obliczenia jak nast puje. Najpierw obliczamy stałe c

1

i d

1

:

c

1

1

0

x

m

1

1

0

x

k

1 1

d

1

2f

1

x

2

2 p

1

x

k

1

A x

1

0

x

nast pnie wyznaczamy pozostałe stałe a do c

N

i d

N

:

c

i

1

m

i

k

i

c

i 1

d

i

2f

i

x

2

2 p

i

x

k

i

c

i 1

d

i 1

Na podstawie c

N

i d

N

obliczamy y

N

:

y

N

2B x

1

d

N

c

N 1

d

N 1

2

0

x

1

c

N 1

1

c

N

po czym kolejne warto ci y

i

(to, e wzory przewiduj mo liwo obliczenia warto ci y

o

wynika z faktu, e warunki

brzegowe mog nie zawiera warto ci funkcji na brzegu — gdy nie s to warunki pierwszego rodzaju):

y

i

c

i

d

i

y

i 1

y

0

1

y

1

A x

1

0

x

Wyprowadzenie: Rozwi zywane równanie ma posta :

y

p x y

q x y

f x

a warunki brzegowe:

0

y a

1

y a

A

0

y b

1

y b

B

W postaci ró nicowej:

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

28

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

i 1

2y

i

y

i 1

h

2

p

i

y

i 1

y

i 1

2h

q

i

y

i

f

i

0

y

0

1

y

1

y

0

h

A

0

y

N

1

y

N 1

y

N 1

2h

B

Warto zauwa y , e do zapisania prawego warunku brzegowego (dla x=b) wykorzystano ró nic centraln — zapis

ró nicowy nie jest wi c symetryczny wzgl dem x i nale y oczekiwa , e wyniki mog by ró ne w zale no ci od tego od

którego ko ca przedziału zaczniemy rozwi zywa .

Po przekształceniu otrzymamy:

y

i 1

m

i

y

i

k

i

y

i 1

i

gdzie:

m

i

2q

i

h

2

4

2 hp

i

;

k

i

2 hp

i

2 hp

i

;

i

2h

2

f

i

2 hp

i

Jest to oczywi cie układ równa opisany macierz trójdiagonaln . Załó my, e dla i = 1,2,...,n–1 mo liwe jest zapisanie

rozwi za w postaci:

y

i

c

i

d

i

y

i 1

Je li jest to prawd , to mo emy podstawi taki wzór na y

i–1

do ró nicowej postaci rozwi zywanego równania:

y

i 1

c

i 1

d

i 1

y

i

y

i 1

m

i

y

i

k

i

c

i 1

d

i 1

y

i

i

co po przekształceniach daje:

y

i

1

m

i

k

i

c

i 1

i

k

i

c

i 1

d

i 1

y

i 1

czyli:

c

i

1

m

i

k

i

c

i 1

;

d

i

i

k

i

c

i 1

d

i 1

Dla i = 1 otrzymamy:

y

2

m

1

y

1

k

1

y

0

1

Z warunku brzegowego dla x = a mamy:

0

y

0

1

y

1

y

0

h

A

y

0

Ah

1

y

1

0

h

1

po podstawieniu:

y

2

m

1

y

1

k

1

Ah

1

y

1

0

h

1

1

i po przekształceniu:

m

1

1

k

1

o

h

1

y

1

y

2

1

k

1

Ah

0

h

1

y

1

0

h

1

m

1

0

h

1

1

k

1

1

k

1

Ah

0

h

1

y

2

co daje nam warto ci c

1

i d

1

. Na podstawie wzorów na c

i

i d

i

mo emy obliczy wszystkie kolejne warto ci współczynników

c i d. Nast pnie rozwa ymy układ równa (który umo liwi nam wyeliminowanie z oblicze warto ci y

N+1

):

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

29

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

y

N

c

N

d

N

y

N 1

y

N 1

c

N 1

d

N 1

y

N

0

y

N

1

y

N 1

y

N 1

2h

B

po przekształceniach otrzymamy:

y

N 1

d

N

y

N

c

N

0

y

N

1

2h

d

N

y

N

c

N

c

N 1

d

N 1

c

N 1

y

N

B

i dalej:

2h

0

1

c

N 1

1

c

N

y

N

2hB

1

c

N 1

d

N 1

d

N

Ostatecznie warto y

N

wyniesie:

y

N

2hB

1

c

N 1

d

N 1

d

N

2h

0

1

c

N 1

1

c

N

Po obliczeniu y

N

mo emy „cofaj c” si obliczy wszystkie pozostałe warto ci y, a do y

1

.

Lewy warunek brzegowy zapisany w postaci ró nicowej mo emy przekształci do postaci umo liwiaj cej obliczenie

y

0

:

y

0

Ah

1

y

1

1

0

h

Metoda ró nicowa dla równa nieliniowych drugiego rz du

Rozwa my teraz zagadnienie brzegowe dla nieliniowego równania drugiego rz du:

y

f x,y,y

z warunkami brzegowymi:

0

y a

1

y a

A

0

y b

1

y b

B

Równanie ró niczkowe zast pimy nast puj cym równaniem ró nicowym:

y

i 1

2y

i

y

i 1

h

2

f x

i

,y

i

,

y

i 1

y

i 1

2h

0

y

0

1

y

1

y

0

h

A

;

0

y

N

1

y

N

y

N 1

h

B

Jest to równanie nieliniowe, rozwi zujemy je metod iteracyjn . Oznaczmy:

0

y

0

y

0

1

y

1

y

0

h

N

y

0

y

N

1

y

N

y

N 1

h

iteracje b d mie posta (indeks górny (r) oznacza numer iteracji):

y

(r 1)

i 1

2y

(r 1)

i

y

(r 1)

i 1

h

2

f x

i

,y

(r )

i

,

y

(r )

i 1

y

(r )

i 1

2h

;

i 1, ,N 1

0

y

(r 1)

i

A

(i 0)

;

N

y

(r 1)

i

B

(i N )

Metody Numeryczne w Fizyce 1

Rozwi zywanie ODE

30

© Andrzej Brozi, Instytut Fizyki Politechniki Łódzkiej

W ramach ka dej z iteracji musimy rozwi za układ równa algebraicznych. Rozwi zanie to mo na przedstawi w

postaci:

y

(r 1)

i

h A

0

b a

A

1

1

B

i

0

B A

0

h

2

N 1

j 1

g

ji

f

(r)

j

gdzie g

ji

i obliczane s jako:

1
h

0 0

b a

0 1

1 0

g

ji

1 j

0

1

h

i

0

0

N

1

h

j i

1 i

0

1

h

j

0

0

N

1

h

j>i

Warto zauwa y , e we wzorze iteracyjnym tylko warto f

j

zale y od numeru iteracji.