« 3. Metody dla równań różniczkowych zwyczajnych  Spis treści  5. Metody dla równań różniczkowych zwyczajnych - teoria zbieżności »  

Zagadnienia

4. Metody dla równań różniczkowych zwyczajnych - rząd schematów

W tym rozdziale zajmiemy się pewnymi własnościami schematów dla równań różniczkowych zwyczajnych. W szczególności przedstawimy pojęcie rzędu schematu oraz zdefiniujemy, co oznacza zbieżność schematu z odpowiednim rzędem.

4.1. Kilka kolejnych schematów

Można postawić pytanie, czy istnieją schematy o wyższej dokładności niż schematy Eulera. Okazuje się, że tak jest i w tym rozdziale przedstawimy kolejne schematy, które dokładniej przybliżają rozwiązanie wyjściowego problemu różniczkowego.

Dość niska dokładność schematów Eulera, którą zaobserwowaliśmy w eksperymentach z rozdziału 3 wynika z tego, że pochodną rozwiązania przybliżyliśmy najprostszym ilorazem różnicowym. W schematach Eulera przybliżamy pochodną poprzez iloraz różnicowy dla parametru h>0 i otrzymujemy:

x⁢t+h-x⁢th-d⁢xd⁢t⁢t=O⁢h (4.1)

o ile x ma ciągłą drugą pochodną w otoczeniu t.

Jeśli x jest bardziej regularna, to pochodną można przybliżyć dokładniej, np. poprzez iloraz różnicowy centralny (pochodna różnicowa centralna)

x⁢t+h-x⁢t-h2⁢h-d⁢xd⁢t⁢t=O⁢h2. (4.2)

Dowód pozostawiamy jako zadanie.

Otrzymujemy w ten sposób:

xn+1=xn-1+2⁢⁢h⁢⁢fn (4.3)

czyli schemat kroku środkowego (midpoint) dla (3.1) .

Schemat midpoint, czyli kroku środkowego, jest dwu-krokowy, tzn. że aby obliczyć xn+1 musimy znać xn i xn-1, czyli trzeba znać x0 i x1.

\
Rys. 4.1. Schematy midpoint i Eulera otwarty na odcinku [0,1].
\
Rys. 4.2. Schematy midpoint i Eulera otwarty na odcinku [0,10].

Policzmy przy pomocy tego schematu rozwiązanie zagadnienia początkowego:

d⁢xd⁢t=a⁢⁢x⁢⁢x⁢0=1

Na początek weźmy a=1 i porównajmy z rozwiązaniem; za x1 do naszych testów schematu midpoint weźmiemy dokładną wartość rozwiązania: exp⁡h, por. rysunek 4.1. Wyraźnie dokładniejszym okazuje się schemat midpoint.

Można się zastanowić, co się stanie na dłuższym odcinku czasu, por. rysunek 4.2. Okazuje się, że schemat midpoint dokładniej działa także w tym przypadku.

Schemat ten nie jest jednak w ogóle używany. W kolejnym rozdziale wyjaśnimy dlaczego.

Inną drogą wprowadzenia nowych schematów jest skorzystanie z rozwinięcia rozwiązania w szereg Taylora: (3.6), tak jak dla schematu Eulera, ale z większą ilością członów. Otrzymujemy w ten sposób np. schemat Taylora:

x⁢t+h≈x⁢t+d⁢xd⁢t⁢t⁢⁢h+12⁢d2⁢xd⁢t2⁢t⁢⁢h2
=x⁢t+f⁢t,x⁢t⁢⁢h+h22⁢⁢∂⁡f∂⁡x⁢t,x⁢t⁢f⁢t,x⁢t+∂⁡f∂⁡t⁢t,x⁢t.

Skorzystaliśmy tu z tego, że d2⁢xd⁢t2=dd⁢t⁢f⁢t,x⁢t=∂⁡f∂⁡x⁢t,x⁢t⁢d⁢xd⁢t⁢t+∂⁡f∂⁡t⁢t,x⁢t.

Schemat Taylora, a dokładniej schemat Taylora rzędu dwa, wygląda następująco:

xn+1=xn+h⁢⁢fn+h22⁢⁢∂x⁡fn⁢⁢fn+∂t⁡fn, (4.4)

gdzie ∂x⁡fn=∂⁡f∂⁡x⁢tn,xn i ∂t⁡fn=∂⁡f∂⁡t⁢tn,xn. W przypadku równania autonomicznego (f⁢t,x=f⁢x) schemat się upraszcza i otrzymujemy:

xn+1=xn+h⁢⁢fn+h22⁢⁢∂x⁡fn⁢⁢fn.

Proszę zauważyć, że ogólnie ∂x⁡f⁢t,x jest macierzą m×m, a ∂t⁡f⁢t,x jest wektorem wymiaru m, czyli koszt schematu Taylora w przypadku wielowymiarowym dla m>1 jest dość duży. Musimy obliczyć w każdym kroku dwa wektory tzn. fn i ∂t⁡fn oraz macierz ∂x⁡fn, wymnożyć tę macierz przez fn i przemnożyć odpowiednie wektory przez h,h2 i dodać je do siebie. Możemy w ten sposób tworzyć kolejne schematy Taylora o coraz większej dokładności - jeśli f jest funkcją dostatecznie gładką. Będą to schematy coraz droższe, szczególnie w przypadku dużego wymiaru m.

Na rysunkach 4.3 i 4.4 widać, że podobnie jak dla schematu midpoint, schemat Taylora jest dokładniejszy niż schemat Eulera otwarty.

\
Rys. 4.3. Rozwiązanie dokładne y′=y z y⁢0=1, rozwiązanie schematem Taylora i otwartym schematem Eulera na [0,1] z h=0.1.
\
Rys. 4.4. Rozwiązanie dokładne y′=y z y⁢0=1, rozwiązanie schematem Taylora i otwartym schematem Eulera na [0,10] z h=0.1.

4.1.1. Zbieżność metod - idea

Błąd schematu (ang. global error) np. Eulera otwartego, czy zamkniętego, czy schematu midpoint zastosowanych do rozwiązywania przybliżonego (3.1) możemy zdefiniować dla ustalonego t∈t0,T jako

Eh⁢t=xnh-x⁢t

dla h=t-t0/n, x rozwiązania (3.1) na t0,T. Zbadajmy, jak zachowuje się błąd Eh⁢t wraz ze zmniejszaniem h w ustalonym t. W szczególności, czy maleje do zera.

Popatrzmy, co pokazują eksperymenty - zastosowaliśmy otwarty schemat Eulera z różnymi krokami do policzenia przybliżenia rozwiązania równania d⁢x/d⁢t=x z x⁢0=1 dla czasu t=1, czyli znamy dokładną wartość rozwiązania x⁢1=e. Ustaliliśmy h0=0.5, a następnie kolejno je połowiliśmy tzn. 2-1⁢h0,…,2-5⁢h0. Wyniki są w tabeli 4.1.

ho⁢t⁢w⁢a⁢r⁢t⁢y⁢⁢E⁢u⁢l⁢e⁢rm⁢i⁢d⁢p⁢o⁢i⁢n⁢tT⁢a⁢y⁢l⁢o⁢r5.0⁢e-01-4.7⁢e-01-7.8⁢e-02-7.0⁢e-022.5⁢e-01-2.8⁢e-01-2.3⁢e-02-2.4⁢e-021.2⁢e-01-1.5⁢e-01-6.4⁢e-03-6.6⁢e-036.2⁢e-02-8.0⁢e-02-1.7⁢e-03-1.7⁢e-033.1⁢e-02-4.1⁢e-02-4.3⁢e-04-4.4⁢e-041.6⁢e-02-2.1⁢e-02-1.1⁢e-04-1.1⁢e-047.8⁢e-03-1.1⁢e-02-2.7⁢e-05-2.8⁢e-05
Tabela 4.1. Błąd dla schematów: otwartego, Eulera, schematu midpoint i schematu Taylora, przybliżających rozwiązanie d⁢x/d⁢t=x z x⁢0=1 dla t=1 czyli exp⁡1.

Widać, że dla schematu Eulera błąd dla zmniejszonego dwukrotnie h maleje dwukrotnie co sugeruje, że błąd zachowuje się jak O⁢h, gdy dla schematów midpoint i schematu Taylora błąd maleje czterokrotnie, czyli zachowuje się jak O⁢h2.

W schemacie midpoint przybliżamy pochodną różnicą centralną, dla której zachodzi:

x⁢t+h-x⁢t-h/2⁢⁢h=d⁢xd⁢t⁢t+O⁢h2

dla dostatecznie gładkiej funkcji, a w przypadku otwartego schematu Eulera - zwykłym ilorazem różnicowym

x⁢t+h-x⁢t/h=d⁢xd⁢t⁢t+O⁢h.

Przy konstrukcji schematu Taylora wykorzystujemy więcej członów z rozwinięcia rozwiązania w szereg Taylora (3.6). Każdy dodatkowy człon z szeregu Taylora powinien podwyższyć dokładność danego schematu.

Dlatego też wprowadza się pojęcie rzędu lokalnego błędu schematu (ang. local truncation error), czyli rzędu schematu. Badamy lokalny błąd schematu względem parametru h, jeśli wstawimy za xn dokładną wartość rozwiązania x⁢tn. Najpierw zdefiniujmy samo pojęcie schematu rozwiązywania (3.1), potem zbieżności schematu i rzędu schematu.

Definicja 4.1

Schematem k krokowym rozwiązywania zadania początkowego (3.1) ze stałym krokiem h>0 na odcinku [t0,T]) nazywamy równanie różnicowe:

xn=Φ⁢h,tn,xn-k⁢…,xn-1,xn⁢⁢n≥k (4.5)

z warunkami startowymi x0,…,xk-1 dla tn=t0+n⁢⁢h. Jeśli Φ nie zależy od xn, to schemat nazywamy otwartym (ang. explicit). W przeciwnym razie - schemat nazywamy zamkniętym (ang. implicit).

Schematy konstruujemy tak, aby dla ustalonego h zachodziło xn≈x⁢tn.

Definicja 4.2

Niech x∈C1⁢t0,T rozwiązania zagadnienia początkowego (3.1). Błąd schematu k krokowego postaci (4.5) dla t=t0+n⁢⁢h∈t0,T definiujemy jako

Eh⁢t=xnh-x⁢t⁢t=t0+n⁢⁢h,

a błąd globalny (ang. global error) na t0,T jako

Eh=maxn=0,…,N⁡Eh⁢tnh

dla N=T-t0/h. Schemat jest zbieżny na t0,T, jeśli

Eh→0⁢⁢h→0,

a jest zbieżny z rzędem p (ang. convergent with order p) (rząd błędu globalnego wynosi p), jeśli dodatkowo

Eh≤C⁢⁢hp

dla pewnej stałej C>0 niezależnej od h (zazwyczaj zależnej od rozwiązania x (3.1) i T-t0).

Definicja 4.3

Niech x∈C1⁢t0,T będzie rozwiązaniem zagadnienia początkowego (3.1). Dla parametru h>0 i schematu k krokowego postaci (4.5) błąd lokalny (ang. local truncation error) definiujemy jako

eh=maxt∈t0,T-k⁢⁢h⁡x⁢t+k⁢⁢h-Φ⁢t,…,t+k⁢⁢h,x⁢t,…,x⁢t+k⁢⁢h.
Definicja 4.4

Schemat (4.5) jest rzędu p (ang. local truncation error is of order p), jeśli dla x∈Cp+1⁢t0,T rozwiązania zagadnienia początkowego (3.1) zachodzi

eh≤C⁢⁢hp+1

dla pewnej dodatniej stałej C niezależnej od h.

Dla otwartego schematu Eulera lokalny błąd schematu jest równy:

eh=maxt∈t0,T-h⁡x⁢t+h-x⁢t-h⁢f⁢t,x⁢t.

Z rozwinięcia w szereg Taylora widzimy, że:

eh=h⁢⁢maxt∈t0,T-h⁡x⁢t+h-x⁢th-d⁢xd⁢t⁢t=O⁢h2

o ile x rozwiązanie (3.1) jest klasy C2, czyli schemat ma rząd jeden. Analogicznie można pokazać, że rząd zamkniętego schematu Eulera jest też jeden, a rząd schematów midpoint i Taylora wynosi dwa. Wykazanie tego, pozostawimy jako zadanie.

4.1.2. Schematy Adamsa

Możemy też wyprowadzić schematy korzystając z równoważnej całkowej wersji zagadnienia początkowego (3.1):

x⁢t+h=x⁢t+∫tt+hd⁢x/d⁢t⁢s⁢d⁢s=x⁢t+∫tt+hf⁢s,x⁢s⁢d⁢s. (4.6)

To prowadzi do konstrukcji całej rodziny schematów (tzw. schematów Adamsa). Jeśli wprowadzimy siatkę równomierną z krokiem h>0, tzn. wprowadzamy tkk=1N dla tkh≡tk=t0+k⁢⁢h, to możemy przybliżyć wartość rozwiązania x⁢tk zastępując w (4.6) całką z jakiejś aproksymacji funkcji f, którą daje się wyliczyć znając wartości fk=f⁢tk,xk dla ustalonej ilości k=n+1,n,n-1,…, np. k=n,n-1,n-2. Wtedy

xn+1=xn+∫tntn+1P⁢s⁢d⁢s,

gdzie P⁢s jest jakimś wielomianem przybliżającym f⁢s,x⁢s zdefiniowanym poprzez wartości odpowiednie fk dla k≤n+1.

W przypadku schematów Adamsa, P⁢t definiujemy jako odpowiedni wielomian interpolacyjny Lagrange'a dla funkcji f z węzłami w punktach tn+j dla j=1,0,-1,.. dla schematu zamkniętego (lub j=0,-1,.. dla schematu otwartego), spełniający odpowiednie warunki interpolacyjne:

P⁢tn+j=fn+j=f⁢tn+j,xn+j

dla p+1 kolejnych indeksów j≤1 dla schematów Adamsa zamkniętych i j<1 dla schematów Adamsa otwartych. Wtedy otrzymujemy klasę zamkniętych schematów Adamsa-Moultona postaci:

xn+1=xn+∑j=-p+11β⌃j⁢fn+j

lub otwartych schematów Adamsa-Bashfortha:

xn+1=xn+∑j=-p0β⌃j⁢fn+j

Przenumerowując indeksy uzyskujemy schemat zamknięty Adamsa-Moultona p krokowy:

xn+p=xn+p-1+h⁢⁢∑j=0pβj⁢fn+j

lub otwarty p+1 krokowy Adamsa-Bashfortha:

xn+p+1=xn+p+h⁢⁢∑j=0pβj⁢fn+j.

Oczywiście w obu przypadkach βj nie zależą od rozwiązania x, ani od f.

W szczególności dla p=1, P⁢t jest wielomianem interpolacyjnym stałym, zdefiniowanym przez wartość w jednym punkcie odpowiednio tn czy tn-1. Dla

p=1⁢⁢P⁢s=fn,

otrzymujemy schemat otwarty Eulera:

xn+1=xn+∫tntn+hfn⁢d⁢s=xn+h⁢⁢fn

Biorąc wartość w punkcie tn+1 uzyskujemy schemat zamknięty Eulera. A dla p=2, P⁢s jest wielomianem liniowym interpolującym f w punktach tn i tn+1. Wtedy otrzymujemy schemat trapezów (ang. trapezoidal scheme):

P⁢s=h-1⁢⁢tn+1-s⁢⁢fn+s-tn⁢⁢fn+1,

czyli

xn+1=xn+∫tntn+hP⁢s⁢d⁢s=xn+0.5⁢⁢h⁢⁢fn+fn+1. (4.7)

Można pokazać, że schemat trapezów jest rzędu dwa.

W przypadku, gdy punkt tn+1 czyli fn+1 nie jest uwzględniony w definicji P⁢s tzn. βp+1=0 rozpatrujemy otwarte schematy Adamsa, które też nazywamy schematami Adamsa-Bashforda, np. schemat otwarty Eulera. W przeciwnym przypadku otrzymujemy schematy zamknięte, które nazywamy schematami Adamsa-Moultona: np. schemat zamknięty Eulera lub schemat trapezów.

4.2. Schematy liniowe wielokrokowe

Definicja 4.5

Dla zadania początkowego (3.1) schematem liniowym wielokrokowym (ang. linear multistep) - dokładniej k krokowym dla stałego kroku dla h=T-t0n nazywamy równanie różnicowe:

∑j=0kαj⁢xn+jh=h⁢⁢∑j=0kβj⁢fj+nh⁢⁢n≥0 (4.8)

z αk≠0 i fjh=f⁢tj,xjh dla tj=t0+j⁢⁢h.

Jeśli βk≠0, to schemat nazywamy zamkniętym, a w przeciwnym wypadku mówimy o schemacie otwartym.

Jeśli znamy x0h,x1h,…,xk-1h to możemy wyliczyć rozwiązanie schematu xjh≈x⁢tj dla tj=t0+j⁢⁢h i j≥k (o ile ono istnieje, co w przypadku schematów zamkniętych nie jest oczywiste).

Zgodnie z Definicją 4.4 schemat liniowy k-krokowy ma rząd p≥1 jeśli dla x∈Cp+1⁢t0,T rozwiązania zagadnienia (3.1) dla t∈[t0,T]) takich, że t+k⁢⁢h≤T lokalny błąd schematu spełnia

eh⁢t:=∑j=0kαj⁢x⁢t+j⁢⁢h-h⁢⁢∑j=0kd⁢xd⁢t⁢t+j⁢⁢h≤C⁢⁢hp+1 (4.9)

ze stałą niezależną od C, czyli eh=O⁢hp+1.

Jeśli za xn weźmiemy wartości rozwiązania w punktach czasu tn, to błąd schematu wynosi O⁢hp+1 (dla gładkiego rozwiązania).

Oczywiście schematy Adamsa opisane w rozdziale 4.1.2 są szczególnym przypadkiem schematów liniowych wielokrokowych. Tak więc schematy: otwarty i zamknięty Eulera, schemat midpoint, lub schemat trapezów są schematami wielokrokowymi liniowymi - w myśl naszej definicji.

4.3. Schematy jednokrokowe

W tym podrozdziale wprowadzimy pojęcie schematu jednokrokowego:

Definicja 4.6

Dla zadania początkowego (3.1) schematem jednokrokowym dla stałego kroku h=T-t0N nazywamy równanie różnicowe:

xn+1=xn+h⁢⁢ϕ⁢h,tn,xn,xn+1⁢⁢n=0,…,N (4.10)

gdzie tj=t0+j⁢⁢h a ϕ jest funkcją ciągłą określoną na 0,H×t0,T×Ux0×Ux0 dla Ux0 otoczenia x0. Dodatkowo, jeśli ϕ nie zależy od xn+1, to schemat jednokrokowy nazywamy otwartym, a w przeciwnym wypadku mówimy o schemacie zamkniętym.

W przypadku schematów otwartych możemy wyliczyć xn+1 znając xn, natomiast w przypadku schematów zamkniętych musimy rozwiązać liniowy, bądź nieliniowy układ równań. Do tej pory poznaliśmy dwa schematy jednokrokowe (które zarazem są schematami liniowymi wielokrokowymi) - czyli oba schematy Eulera i schemat trapezów.

Analogicznie do przypadku schematów liniowych wielokrokowych, zgodnie z Definicją 4.4, schemat jednokrokowy ma rząd p≥1, jeśli dla x∈Cp+1⁢t0,T rozwiązania zagadnienia (3.1) dla 1≤p, h=T-t0N i t∈[t0,T]) lokalny błąd schematu spełnia:

eh⁢t:=x⁢t+h-x⁢t-h⁢⁢ϕ⁢h,t,d⁢xd⁢t⁢t,d⁢xd⁢t⁢t+h≤C⁢⁢hp+1,

ze stałą C niezależną od t, czyli eh=O⁢hp+1 dla dostatecznie gładkiego rozwiązania.

4.3.1. Schematy Rungego-Kutty

Podstawową klasą schematów jednokrokowych są tzw. schematy Rungego-Kutty lub - mówiąc krótko - schematy Rungego. Idea ich jest prosta.

Załóżmy, że znamy xn, i chcemy wyliczyć wartość xn+1 ze wzoru uwzględniającego wartość pola wektorowego nie tylko w xn, ale również w dodatkowym punkcie x~. Wtedy

xn+1=F⁢h,tn,xn,x~.

Biorąc schemat otwarty Eulera z krokiem h~ otrzymujemy punkt

x~=xn+h~⁢⁢f⁢t,x⁢t,

który, jak wiemy, przybliża x⁢t+h~, ale niedokładnie. Możemy policzyć wartość pola wektorowego f w tym punkcie i następnie, wykorzystując wartość yn=f⁢tn,xn i y~=f⁢tn+h~,x~, znaleźć lepsze przybliżenie x⁢tn+1 - czyli np. za przybliżenie pola wektorowego wziąć ważoną średnią obu wartości b⁢⁢yn+c⁢y~ dla pewnych ustalonych wag b,c. Możliwości jest wiele. Pojawia się pytanie: jak oceniać różne konstrukcje F? Można tak dobierać F, aby rząd schematu był możliwie duży.

Załóżmy, że h~=a⁢⁢h. Wtedy szukamy schematu postaci:

xn+1=xn+b⁢⁢h⁢⁢fn+c⁢⁢h⁢⁢f⁢tn+a⁢⁢h,xn+a⁢⁢h⁢⁢f⁢tn,xn (4.11)

tak, aby schemat miał maksymalny rząd.

Rozwijamy rozwiązanie x w szereg Taylora:

x⁢t+h-x⁢t=h⁢⁢d⁢xd⁢t⁢t+0.5⁢⁢h2⁢⁢d2⁢xd⁢t2⁢t+O⁢h3

i rozwijając ostatni z członów (4.11) w punkcie t,x otrzymujemy:

f⁢t+a⁢⁢h,x+a⁢⁢h⁢⁢f⁢t,x=f+a⁢⁢h⁢⁢fx⁢⁢f+a⁢⁢h⁢⁢ft
=f+a⁢⁢h⁢fx⁢⁢f+ft
=d⁢xd⁢t+a⁢⁢h⁢⁢d2⁢xd⁢t2.

Skorzystaliśmy z tego, że d2⁢xd⁢t2=dd⁢t⁢f⁢t,x⁢t=fx⁢⁢f+ft. Zatem, wstawiając dwa ostatnie równania do (4.11) otrzymujemy warunki na to, aby schemat był rzędu dwa:

b+c=1
c⁢⁢a=0.5
\
Rys. 4.5. Schemat Heuna w porównaniu ze schematem otwartym Eulera na [0,1] z h=0.1.
\
Rys. 4.6. Zmodyfikowany schemat Euler w porównaniu ze schematem otwartym Eulera na [0,1] z h=0.1.

Tak więc otrzymaliśmy całą rodzinę schematów Rungego-Kutty rzędu dwa, np.:

  1. Zmodyfikowany schemat Eulera

    xn+1=xn+h⁢⁢f⁢tn+h2,xn+h2⁢⁢fn (4.12)

    dla c=1,b=0,a=0.5,

  2. Schemat Heuna

    xn+1=xn+h2⁢⁢fn+f⁢tn+1,xn+h⁢⁢fn, (4.13)

    dla b=c=0.5;a=1.

Warto zauważyć, że w niektórych publikacjach wszystkie schematy otwarte Rungego-Kutty rzędu dwa nazywane są zmodyfikowanym schematem Eulera.

\
Rys. 4.7. Schemat Heuna na [0,100] i exp⁡t. Wykres dla t bliskich 100.

Na rysunkach 4.6 i 4.5 pokazano rozwiązania uzyskane tymi dwoma schematami dla zadania d⁢x/d⁢t=x z x⁢0=1 na 0,1. Widać, że wykresy się pokrywają z wykresem rozwiązania. W porównaniu do otwartego schematu Eulera widzimy znaczącą poprawę. Zobaczmy, co się dzieje na dłuższym odcinku czasu w przypadku schematu Heuna, por. rysunek 4.7.

4.3.1.1. Graficzne wytłumaczenie zmodyfikowanego schematu Eulera

Na rysunku 4.8 zawarliśmy graficzne wytłumaczenie jednego kroku zmodyfikowanego schematu Eulera.

\
Rys. 4.8. Graficzne wytłumaczenie zmodyfikowanego schematu Eulera.

Punkt A - oznacza t,x, czyli wcześniej obliczone przybliżenie dla czasu t. Chcemy wyznaczyć przybliżenie dla czasu t+h. Otwarty schemat Eulera w jednym kroku przyjmuje za różnicę między kolejnymi punktami h⁢⁢f⁢t,x, czyli idzie w kierunku pola wektorowego w punkcie t, czyli na naszym rysunku daje to punkt C. Z kolei w zmodyfikowanym schemacie Eulera przyjmujemy za różnicę h pomnożone przez kierunek pola wyznaczonego w dodatkowym pomocniczym punkcie x⌃=x+0.5⁢⁢fn oznaczonym jako B, tzn. idziemy w kierunku f⁢t+0.5⁢⁢h,x⌃ i otrzymujemy w efekcie punkt D.

Na rysunku widać, że jeśli nachylenie pola wektorowego mocno się zmienia, to pole w punkcie B powinno mieć lepszy kierunek niż w A, czy w E. Jest to oczywiście argument heurystyczny.

\
Rys. 4.9. Schematy Rungego rzędu cztery i schemat Heuna rzędu dwa na 0,1.
\
Rys. 4.10. Schematy Rungego rzędu cztery i schemat Heuna rzędu dwa na 0,50 w skali pół-logarytmicznej.

Analogicznie konstruuje się schematy Rungego wyższych rzędów poprzez wprowadzenie większej ilości kroków pośrednich, jak również schematy zamknięte Rungego - dopuszczając wartość fn+1 w schemacie.

Podamy kilka wzorów na powszechnie używany otwarty schemat Rungego-Kutty czwartego rzędu.

Najpierw definiujemy cztery wartości:

K1=f⁢tn,xnK2=f⁢tn+h2,xn+h2⁢⁢K1K3=f⁢tn+h2,xn+h2⁢⁢K2K4=f⁢tn+h,xn+h⁢⁢K3 (4.14)

i otrzymujemy ostateczny wzór:

xn+1=xn+h6⁢K1+2⁢⁢K2+2⁢⁢K3+K4 (4.15)

schematu rzędu cztery (co oczywiście wymaga dowodu).

Istnieje oczywiście cała rodzina otwartych schematów Rungego-Kutty rzędu cztery. Tu podaliśmy tylko przykładowy schemat z tej rodziny.

Popatrzmy jak działa ten schemat w porównaniu ze schematem Heuna dla naszego modelowego zagadnienia d⁢x/d⁢t=x z x⁢0=1.

Na odcinku 0,1 dla h=0.1 oba schematy praktycznie pokrywają się z rozwiązaniem, por. rysunek 4.9.

\
Rys. 4.11. Schematy Rungego rzędu cztery i schemat Heuna rzędu dwa na 44,50.

Na rysunku 4.10 widać wykresy rozwiązań na 0,50 (w skali pół-logarytmicznej), a na rysunku 4.11 wykres na odcinku 44,50 - tu już widać, że schemat rzędu cztery jest jednak znacznie lepszy.

W rozdziale 5.4 przedstawimy metodę eksperymentalną badania rzędu schematu. Przy okazji zobaczymy ogromną różnicę w dokładności obu schematów, por. Tabela 5.2.

4.4. Zadania

Ćwiczenie 4.1

Udowodnij, że wzory (4.1) i (4.2) są prawdziwe dla funkcji dostatecznie gładkich.

Wskazówka: 

Rozwiń funkcje w szereg Taylora.

Ćwiczenie 4.2

Pokaż, że rząd schematów Eulera wynosi jeden, o ile rozwiązanie zadania Cauchy'ego jest klasy C2.

Ćwiczenie 4.3

Pokaż, że rząd schematu kroku środkowego wynosi k, o ile rozwiązanie zadania Cauchy'ego jest klasy Ck+1 dla k=1,2.

Ćwiczenie 4.4

Znajdź rząd schematu Taylora (4.4) dla rozwiązania dostatecznie gładkiego.

Ćwiczenie 4.5

Znajdź wzór na schemat Taylora rzędu trzy.

Ćwiczenie 4.6

Rozpatrzmy rodzinę schematów:

xn+1=xn+h⁢⁢a⁢⁢fn+b⁢⁢fn+1

Określ rząd schematu w zależności od wartości parametrów a,b. Dla jakich a,b rząd jest największy? Dla jakich wartości parametrów schemat będzie zamknięty, a dla jakich otwarty?

Ćwiczenie 4.7

Wyprowadź otwarty dwukrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia jeden (żeby policzyć xn+2 potrzebujemy h,tn,xn+1, fn i fn+1, tak jak opisano w rozdziale 4.1.2). Zbadaj rząd schematu.

Rozwiązanie: 

Zgodnie z zasadą konstrukcji schematów Adamsa musimy scałkować na odcinku tn,tn+h wielomian liniowy interpolacyjny P1⁢s=fn+1-h⁢⁢fn-fn-1⁢s-tn. Otrzymujemy schemat xn+1=xn+0.5⁢⁢h⁢⁢3⁢⁢fn-fn-1.

Ćwiczenie 4.8

Wyprowadź otwarty trzykrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia dwa (żeby policzyć xn+3 potrzebujemy h,tn,xn+2 i fn,fn+1,fn+2, tak jak opisano w rozdziale 4.1.2). Zbadaj rząd schematu.

Ćwiczenie 4.9

Wyprowadź zamknięty dwukrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia dwa (żeby policzyć xn+2 potrzebujemy h,tn,xn+1 i fn,fn+1,fn+2, tak jak opisano w rozdziale 4.1.2). Zbadaj rząd schematu.

Ćwiczenie 4.10 (średnio trudne)

Rozpatrzmy rodzinę schematów:

xn+2=xn+1+h⁢⁢a⁢⁢fn+b⁢⁢fn+1+c⁢⁢fn+2.

Określ rząd schematu w zależności od wartości parametrów a,b,c. Dla jakich wartości rząd jest największy? Dla jakich wartości parametrów schemat będzie zamknięty, a dla jakich otwarty?

Ćwiczenie 4.11 (trudne)

Udowodnij, że schemat (4.15) ma rząd cztery.

Ćwiczenie 4.12 (laboratoryjne)

Zbadaj eksperymentalnie metodą połowienia kroków rząd lokalnego błędu schematu dla schematów:

  • otwartego schematu Eulera (3.4),

  • zamkniętego schematu Eulera (3.5),

  • schematu midpoint (4.3),

  • schematu Taylora (4.4),

  • schematu Heuna (4.13),

  • schematu trapezów (4.7),

  • zmodyfikowanego schematu Eulera (4.12),

  • schematu Rungego rzędu cztery (4.15).

zastosowanych do modelowego zadania d⁢xd⁢t=1+x2 z x⁢0=0 z rozwiązaniem x⁢t=tan⁡t. Tzn. dla ustalonego t, np. t=1 i kolejnych połowionych kroków hk=hk-12=2-k⁢h0 z h0=1⁢e-1 liczymy lokalny błąd schematu ehk⁢t, por. (4.9), i następnie stosunek ehk⁢tehk+1⁢t. Jeśliby ten stosunek wynosił w przybliżeniu 2p+1, to lokalny błąd schematu zachowuje się jak O⁢hp+1, oznacza to, że schemat posiada rząd p przynajmniej dla tego zadania początkowego.

« 3. Metody dla równań różniczkowych zwyczajnych  Spis treści  5. Metody dla równań różniczkowych zwyczajnych - teoria zbieżności »  

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

strona główna | webmaster | o portalu | pomoc

© Wydział Matematyki, Informatyki i Mechaniki UW, 2009-2010. Niniejsze materiały są udostępnione bezpłatnie na licencji Creative Commons Uznanie autorstwa-Użycie niekomercyjne-Bez utworów zależnych 3.0 Polska.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.