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:

\left|\frac{x(t+h)-x(t)}{h}-\frac{dx}{dt}(t)\right|=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)

\left|\frac{x(t+h)-x(t-h)}{2h}-\frac{dx}{dt}(t)\right|=O(h^{2}). (4.2)

Dowód pozostawiamy jako zadanie.

Otrzymujemy w ten sposób:

x_{{n+1}}=x_{{n-1}}+2\, h\, f_{n} (4.3)

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

Schemat midpoint, czyli kroku środkowego, jest dwu-krokowy, tzn. że aby obliczyć x_{{n+1}} musimy znać x_{n} i x_{{n-1}}, czyli trzeba znać x_{0} i x_{1}.

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

\frac{dx}{dt}=a\, x\qquad x(0)=1

Na początek weźmy a=1 i porównajmy z rozwiązaniem; za x_{1} 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:

\displaystyle x(t+h) \displaystyle\approx \displaystyle x(t)+\frac{dx}{dt}(t)\, h+\frac{1}{2}\frac{d^{2}x}{dt^{2}}(t)\, h^{2}
\displaystyle= \displaystyle x(t)+f(t,x(t))\, h+\frac{h^{2}}{2}\,(\frac{\partial f}{\partial x}(t,x(t))f(t,x(t))+\frac{\partial f}{\partial t}(t,x(t))).

Skorzystaliśmy tu z tego, że \frac{d^{2}x}{dt^{2}}=\frac{d}{dt}f(t,x(t))=\frac{\partial f}{\partial x}(t,x(t))\frac{dx}{dt}(t)+\frac{\partial f}{\partial t}(t,x(t)).

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

x_{{n+1}}=x_{n}+h\, f_{n}+\frac{h^{2}}{2}\,(\partial _{x}f_{n}\, f_{n}+\partial _{t}f_{n}), (4.4)

gdzie \partial _{x}f_{n}=\frac{\partial f}{\partial x}(t_{n},x_{n}) i \partial _{t}f_{n}=\frac{\partial f}{\partial t}(t_{n},x_{n}). W przypadku równania autonomicznego (f(t,x)=f(x)) schemat się upraszcza i otrzymujemy:

x_{{n+1}}=x_{n}+h\, f_{n}+\frac{h^{2}}{2}\,(\partial _{x}f_{n}\, f_{n}).

Proszę zauważyć, że ogólnie \partial _{x}f(t,x) jest macierzą m\times m, a \partial _{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. f_{n} i \partial _{t}f_{n} oraz macierz \partial _{x}f_{n}, wymnożyć tę macierz przez f_{n} i przemnożyć odpowiednie wektory przez h,h^{2} 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^{{\prime}}=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^{{\prime}}=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\in[t_{0},T] jako

E_{h}(t)=|x_{n}^{h}-x(t)|

dla h=(t-t_{0})/n, x rozwiązania (3.1) na [t_{0},T]. Zbadajmy, jak zachowuje się błąd E_{h}(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 dx/dt=x z x(0)=1 dla czasu t=1, czyli znamy dokładną wartość rozwiązania x(1)=e. Ustaliliśmy h_{0}=0.5, a następnie kolejno je połowiliśmy tzn. 2^{{-1}}h_{0},\ldots,2^{{-5}}h_{0}. Wyniki są w tabeli 4.1.

\begin{array}[]{|c|c|c|c|}\hline h&otwarty\; Euler&midpoint&Taylor\\
\hline 5.0e-01&-4.7e-01&-7.8e-02&-7.0e-02\\
2.5e-01&-2.8e-01&-2.3e-02&-2.4e-02\\
1.2e-01&-1.5e-01&-6.4e-03&-6.6e-03\\
6.2e-02&-8.0e-02&-1.7e-03&-1.7e-03\\
3.1e-02&-4.1e-02&-4.3e-04&-4.4e-04\\
1.6e-02&-2.1e-02&-1.1e-04&-1.1e-04\\
7.8e-03&-1.1e-02&-2.7e-05&-2.8e-05\\
\hline\end{array}
Tabela 4.1. Błąd dla schematów: otwartego, Eulera, schematu midpoint i schematu Taylora, przybliżających rozwiązanie dx/dt=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(h^{2}).

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

(x(t+h)-x(t-h))/(2\, h)=\frac{dx}{dt}(t)+O(h^{2})

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

(x(t+h)-x(t))/(h)=\frac{dx}{dt}(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 x_{n} dokładną wartość rozwiązania x(t_{n}). 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 [t_{0},T]) nazywamy równanie różnicowe:

x_{n}=\Phi(h,t_{n},x_{{n-k}}\ldots,x_{{n-1}},x_{n})\qquad n\geq k (4.5)

z warunkami startowymi x_{0},\ldots,x_{{k-1}} dla t_{n}=t_{0}+n\, h. Jeśli \Phi nie zależy od x_{n}, to schemat nazywamy otwartym (ang. explicit). W przeciwnym razie - schemat nazywamy zamkniętym (ang. implicit).

Schematy konstruujemy tak, aby dla ustalonego h zachodziło x_{n}\approx x(t_{n}).

Definicja 4.2

Niech x\in C^{1}([t_{0},T)) rozwiązania zagadnienia początkowego (3.1). Błąd schematu k krokowego postaci (4.5) dla t=t_{0}+n\, h\in[t_{0},T] definiujemy jako

E_{h}(t)=|x_{n}^{h}-x(t)|\qquad(t=t_{0}+n\, h),

a błąd globalny (ang. global error) na [t_{0},T] jako

E_{h}=\max _{{n=0,\ldots,N}}E_{h}(t_{n}^{h})

dla N=(T-t_{0})/h. Schemat jest zbieżny na [t_{0},T], jeśli

E_{h}\rightarrow 0\qquad h\rightarrow 0,

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

E_{h}\leq C\, h^{p}

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

Definicja 4.3

Niech x\in C^{1}([t_{0},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

e_{h}=\max _{{t\in[t_{0},T-k\, h]}}|x(t+k\, h)-\Phi(t,\ldots,t+k\, h,x(t),\ldots,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\in C^{{p+1}}([t_{0},T)) rozwiązania zagadnienia początkowego (3.1) zachodzi

e_{h}\leq C\, h^{{p+1}}

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

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

e_{h}=\max _{{t\in[t_{0},T-h]}}|x(t+h)-x(t)-hf(t,x(t))|.

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

e_{h}=h\,\max _{{t\in[t_{0},T-h]}}\left|\frac{x(t+h)-x(t)}{h}-\frac{dx}{dt}(t)\right|=O(h^{2})

o ile x rozwiązanie (3.1) jest klasy C^{2}, 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)+\int _{t}^{{t+h}}dx/dt(s)ds=x(t)+\int _{t}^{{t+h}}f(s,x(s))ds. (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 \{ t_{k}\} _{{k=1}}^{N} dla t_{k}^{h}\equiv t_{k}=t_{0}+k\, h, to możemy przybliżyć wartość rozwiązania x(t_{k}) zastępując w (4.6) całką z jakiejś aproksymacji funkcji f, którą daje się wyliczyć znając wartości f_{k}=f(t_{k},x_{k}) dla ustalonej ilości k=n+1,n,n-1,\ldots, np. k=n,n-1,n-2. Wtedy

x_{{n+1}}=x_{n}+\int _{{t_{n}}}^{{t_{{n+1}}}}P(s)ds,

gdzie P(s) jest jakimś wielomianem przybliżającym f(s,x(s)) zdefiniowanym poprzez wartości odpowiednie f_{k} dla k\leq 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 t_{{n+j}} dla j=1,0,-1,.. dla schematu zamkniętego (lub j=0,-1,.. dla schematu otwartego), spełniający odpowiednie warunki interpolacyjne:

P(t_{{n+j}})=f_{{n+j}}=f(t_{{n+j}},x_{{n+j}})

dla p+1 kolejnych indeksów j\leq 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:

x_{{n+1}}=x_{n}+\sum _{{j=-p+1}}^{1}\hat{\beta}_{j}f_{{n+j}}

lub otwartych schematów Adamsa-Bashfortha:

x_{{n+1}}=x_{n}+\sum _{{j=-p}}^{0}\hat{\beta}_{j}f_{{n+j}}

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

x_{{n+p}}=x_{{n+p-1}}+h\,\sum _{{j=0}}^{p}\beta _{j}f_{{n+j}}

lub otwarty p+1 krokowy Adamsa-Bashfortha:

x_{{n+p+1}}=x_{{n+p}}+h\,\sum _{{j=0}}^{p}\beta _{j}f_{{n+j}}.

Oczywiście w obu przypadkach \beta _{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 t_{n} czy t_{{n-1}}. Dla

p=1\qquad P(s)=f_{n},

otrzymujemy schemat otwarty Eulera:

x_{{n+1}}=x_{n}+\int _{{t_{n}}}^{{t_{n}+h}}f_{n}ds=x_{n}+h\, f_{n}

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

P(s)=h^{{-1}}\,((t_{{n+1}}-s)\, f_{n}+(s-t_{n})\, f_{{n+1}}),

czyli

x_{{n+1}}=x_{n}+\int _{{t_{n}}}^{{t_{n}+h}}P(s)ds=x_{n}+0.5\, h\,(f_{n}+f_{{n+1}}). (4.7)

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

W przypadku, gdy punkt t_{{n+1}} czyli f_{{n+1}} nie jest uwzględniony w definicji P(s) tzn. \beta _{{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=\frac{T-t_{0}}{n} nazywamy równanie różnicowe:

\sum _{{j=0}}^{k}\alpha _{j}x_{{n+j}}^{h}=h\,\sum _{{j=0}}^{k}\beta _{j}f^{h}_{{j+n}}\qquad n\geq 0 (4.8)

z \alpha _{k}\not=0 i f^{h}_{j}=f(t_{j},x_{j}^{h}) dla t_{j}=t_{0}+j\, h.

Jeśli \beta _{k}\not=0, to schemat nazywamy zamkniętym, a w przeciwnym wypadku mówimy o schemacie otwartym.

Jeśli znamy x^{h}_{0},x^{h}_{1},\ldots,x^{h}_{{k-1}} to możemy wyliczyć rozwiązanie schematu x^{h}_{j}\approx x(t_{j}) dla t_{j}=t_{0}+j\, h i j\geq 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\geq 1 jeśli dla x\in C^{{p+1}}([t_{0},T]) rozwiązania zagadnienia (3.1) dla t\in[t_{0},T]) takich, że t+k\, h\leq T lokalny błąd schematu spełnia

e_{h}(t):=\left|\sum _{{j=0}}^{k}\alpha _{j}x(t+j\, h)-h\,\sum _{{j=0}}^{k}\frac{dx}{dt}(t+j\, h)\right|\leq C\, h^{{p+1}} (4.9)

ze stałą niezależną od C, czyli e_{h}=O(h^{{p+1}}).

Jeśli za x_{n} weźmiemy wartości rozwiązania w punktach czasu t_{n}, to błąd schematu wynosi O(h^{{p+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=\frac{T-t_{0}}{N} nazywamy równanie różnicowe:

x_{{n+1}}=x_{n}+h\,\phi(h,t_{n},x_{n},x_{{n+1}})\qquad n=0,\ldots,N (4.10)

gdzie t_{j}=t_{0}+j\, h a \phi jest funkcją ciągłą określoną na [0,H)\times[t_{0},T)\times U_{{x_{0}}}\times U_{{x_{0}}} dla U_{{x_{0}}} otoczenia x_{0}. Dodatkowo, jeśli \phi nie zależy od x_{{n+1}}, to schemat jednokrokowy nazywamy otwartym, a w przeciwnym wypadku mówimy o schemacie zamkniętym.

W przypadku schematów otwartych możemy wyliczyć x_{{n+1}} znając x_{n}, 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\geq 1, jeśli dla x\in C^{{p+1}}([t_{0},T)) rozwiązania zagadnienia (3.1) dla 1\leq p, h=\frac{T-t_{0}}{N} i t\in[t_{0},T]) lokalny błąd schematu spełnia:

e_{h}(t):=|x(t+h)-x(t)-h\,\phi(h,t,\frac{dx}{dt}(t),\frac{dx}{dt}(t+h))|\leq C\, h^{{p+1}},

ze stałą C niezależną od t, czyli e_{h}=O(h^{{p+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 x_{n}, i chcemy wyliczyć wartość x_{{n+1}} ze wzoru uwzględniającego wartość pola wektorowego nie tylko w x_{n}, ale również w dodatkowym punkcie \tilde{x}. Wtedy

x_{{n+1}}=F(h,t_{n},x_{n},\tilde{x}).

Biorąc schemat otwarty Eulera z krokiem \tilde{h} otrzymujemy punkt

\tilde{x}=x_{n}+\tilde{h}\, f(t,x(t)),

który, jak wiemy, przybliża x(t+\tilde{h}), ale niedokładnie. Możemy policzyć wartość pola wektorowego f w tym punkcie i następnie, wykorzystując wartość y_{n}=f(t_{n},x_{n}) i \tilde{y}=f(t_{n}+\tilde{h},\tilde{x}), znaleźć lepsze przybliżenie x(t_{{n+1}}) - czyli np. za przybliżenie pola wektorowego wziąć ważoną średnią obu wartości b\, y_{n}+c\tilde{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 \tilde{h}=a\, h. Wtedy szukamy schematu postaci:

x_{{n+1}}=x_{n}+b\, h\, f_{n}+c\, h\, f(t_{n}+a\, h,x_{n}+a\, h\, f(t_{n},x_{n})) (4.11)

tak, aby schemat miał maksymalny rząd.

Rozwijamy rozwiązanie x w szereg Taylora:

x(t+h)-x(t)=h\,\frac{dx}{dt}(t)+0.5\, h^{2}\,\frac{d^{2}x}{dt^{2}}(t)+O(h^{3})

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

\displaystyle f(t+a\, h,x+a\, h\, f(t,x)) \displaystyle= \displaystyle f+a\, h\, f_{x}\, f+a\, h\, f_{t}
\displaystyle= \displaystyle f+a\, h(f_{x}\, f+f_{t})
\displaystyle= \displaystyle\frac{dx}{dt}+a\, h\,\frac{d^{2}x}{dt^{2}}.

Skorzystaliśmy z tego, że \frac{d^{2}x}{dt^{2}}=\frac{d}{dt}f(t,x(t))=f_{x}\, f+f_{t}. Zatem, wstawiając dwa ostatnie równania do (4.11) otrzymujemy warunki na to, aby schemat był rzędu dwa:

\displaystyle b+c \displaystyle= \displaystyle 1
\displaystyle c\, a \displaystyle= \displaystyle 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

    x_{{n+1}}=x_{n}+h\, f\left(t_{n}+\frac{h}{2},x_{n}+\frac{h}{2}\, f_{n}\right) (4.12)

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

  2. Schemat Heuna

    x_{{n+1}}=x_{n}+\frac{h}{2}\,\left(f_{n}+f(t_{{n+1}},x_{n}+h\, f_{n})\right), (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 dx/dt=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 \hat{x}=x+0.5\, f_{n} oznaczonym jako B, tzn. idziemy w kierunku f(t+0.5\, h,\hat{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ść f_{{n+1}} w schemacie.

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

Najpierw definiujemy cztery wartości:

\begin{array}[]{l}K_{1}=f(t_{n},x_{n})\\
K_{2}=f(t_{n}+\frac{h}{2},x_{n}+\frac{h}{2}\, K_{1})\\
K_{3}=f(t_{n}+\frac{h}{2},x_{n}+\frac{h}{2}\, K_{2})\\
K_{4}=f(t_{n}+h,x_{n}+h\, K_{3})\end{array} (4.14)

i otrzymujemy ostateczny wzór:

x_{{n+1}}=x_{n}+\frac{h}{6}(K_{1}+2\, K_{2}+2\, K_{3}+K_{4}) (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 dx/dt=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 C^{2}.

Ćwiczenie 4.3

Pokaż, że rząd schematu kroku środkowego wynosi k, o ile rozwiązanie zadania Cauchy'ego jest klasy C^{{k+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:

x_{{n+1}}=x_{n}+h\,(a\, f_{n}+b\, f_{{n+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ć x_{{n+2}} potrzebujemy h,t_{n},x_{{n+1}}, f_{n} i f_{{n+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 [t_{n},t_{n}+h] wielomian liniowy interpolacyjny P_{1}(s)=f_{n}+\frac{1}{-h}\,(f_{n}-f_{{n-1}})(s-t_{n}). Otrzymujemy schemat x_{{n+1}}=x_{n}+0.5\, h\,(3\, f_{n}-f_{{n-1}}).

Ćwiczenie 4.8

Wyprowadź otwarty trzykrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia dwa (żeby policzyć x_{{n+3}} potrzebujemy h,t_{n},x_{{n+2}} i f_{n},f_{{n+1}},f_{{n+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ć x_{{n+2}} potrzebujemy h,t_{n},x_{{n+1}} i f_{n},f_{{n+1}},f_{{n+2}}, tak jak opisano w rozdziale 4.1.2). Zbadaj rząd schematu.

Ćwiczenie 4.10 (średnio trudne)

Rozpatrzmy rodzinę schematów:

x_{{n+2}}=x_{{n+1}}+h\,(a\, f_{n}+b\, f_{{n+1}}+c\, f_{{n+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 \frac{dx}{dt}=1+x^{2} z x(0)=0 z rozwiązaniem x(t)=\tan(t). Tzn. dla ustalonego t, np. t=1 i kolejnych połowionych kroków h_{k}=\frac{h_{{k-1}}}{2}=2^{{-k}}h_{0} z h_{0}=1e-1 liczymy lokalny błąd schematu e_{{h_{k}}}(t), por. (4.9), i następnie stosunek \frac{e_{{h_{k}}}(t)}{e_{{h_{{k+1}}}}(t)}. Jeśliby ten stosunek wynosił w przybliżeniu 2^{{p+1}}, to lokalny błąd schematu zachowuje się jak O(h^{{p+1}}), oznacza to, że schemat posiada rząd p przynajmniej dla tego zadania początkowego.

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

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.