Zagadnienia

4.1 Kilka kolejnych schematów
- 4.1.1 Zbieżność metod - idea
- 4.1.2 Schematy Adamsa
4.2 Schematy liniowe wielokrokowe
4.3 Schematy jednokrokowe
- 4.3.1 Schematy Rungego-Kutty
  - 4.3.1.1 Graficzne wytłumaczenie zmodyfikowanego schematu Eulera
4.4 Zadania

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

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$ ,
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.