Zagadnienia

5.1 Teoria zbieżności schematów jednokrokowych
5.2 Teoria zbieżności schematów liniowych wielokrokowych
- 5.2.1 Stabilność, zgodność
  - 5.2.1.1 Zgodność schematu
  - 5.2.1.2 Zbieżność
- 5.2.2 Stabilność, a silna stabilność
5.3 Wartości startowe schematów wielokrokowych
5.4 Eksperymentalne badanie rzędu zbieżności schematów
5.5 Zadania

5. Metody dla równań różniczkowych zwyczajnych - teoria zbieżności

W tym rozdziale przedstawimy teorię zbieżności schematów jednokrokowych i wielokrokowych. Rozpatrzymy tylko przypadek skalarny, ale teoria dla układów równań jest analogiczna. Wystarczy tylko moduł zastąpić przez jakąś normę w $\mathbb{R}^{m}$ np. normę euklidesową.

5.1. Teoria zbieżności schematów jednokrokowych

Teoria zbieżności schematów jednokrokowych jest teorią odrębną od teorii dla schematów wielokrokowych liniowych.

Okaże się, że kluczowym pojęciem jest tu zgodność schematu jednokrokowego - inaczej konsystentność, którą definiujemy następująco:

Definicja 5.1

Schemat jednokrokowy (4.10) jest zgodny (konsystentny), (ang. consistent) jeśli:

$\phi$ jest ciągłą ze względu na wszystkie zmienne
$\phi(0,t,x,x)=f(t,x)$ dla wszystkich $(t,x)$ .
$\phi$ jest lipschitzowska ze względu na zmienne $x_{n}$ i $x_{{n+1}}$ tzn. istnieje $L>0$ takie, że dla wszystkich $x_{1},x_{2},y_{1},y_{2}\in U_{{x_{0}}}$

$|\phi(h,t,x_{1},x_{2})-\phi(h,t,y_{1},y_{2})|\leq L\sum _{{k=1}}^{2}|x_{k}-y_{k}|.$

Twierdzenie 5.1 (o zbieżności schematów jednokrokowych)

Jeśli rozwiązanie zagadnienia początkowego (3.1) $x\in C^{{p+1}}([t_{0},T])$ , schemat jednokrokowy jest zgodny i jest rzędu $p\geq 1$ , to ten schemat jest zbieżny z rzędem $p$ .

Dowód zostanie tutaj przedstawiony dla prostoty tylko dla schematów otwartych z rzędem $p$ tj. $\phi(h,t,x,y)=\phi(h,t,x)$ . Dowód w całej ogólności można znaleźć np. w [23].

Oznaczmy przez $E_{n}=x_{n}-x(t_{n})$ , czyli błąd pomiędzy obliczonym schematem przybliżeniem rozwiązania dla czasu $t_{n}$ , a dokładną wartością rozwiązania $x(t_{n})$ . Niech $\tau _{n}=e_{h}(t_{n})=x(t_{{n+1}})-x(t_{n})-h\phi(h,t_{n},x(t_{n}))$ , czyli $\tau _{n}$ to lokalny błąd schematu dla czasu $t_{n}$ . Wtedy otrzymujemy, że

$E_{n}=E_{{n-1}}+h(\phi(h,t_{{n-1}},x_{{n-1}})-\phi(h,t_{{n-1}},x(t_{{n-1}})))-\tau _{{n-1}},$

a stąd, korzystając ze zgodności schematu, a dokładniej lipschitzowskości funkcji $\phi$ , por. Definicja 5.1, otrzymujemy

$|E_{n}|\leq(1+h*L)|E_{{n-1}}|+|\tau _{{n-1}}|.$

Dalej, poprzez indukcję matematyczną otrzymujemy:

$|E_{n}|\leq(1+h*L)^{n}|E_{0}|+\sum _{{k=0}}^{{n-1}}(1+h*L)^{{n-k-1}}|\tau _{k}|$

Korzystając z tego, że ( $|1+x|\leq e^{{|x|}}$ )

$(1+h*L)^{n}\leq e^{{n*h*L}}\leq e^{{L*(T-t_{0})}}$

dla $n$ takich, że $h*n\leq T-t_{0}$ widzimy, że

$|E_{n}|\leq e^{{L*(T-t_{0})}}(|E_{0}|+\sum _{{k=0}}^{{n-1}}|\tau _{k}|)$

Zauważmy, że $E_{0}=0$ . Widzimy też, że

$|\tau _{n}|\leq e_{h}.$

Ponieważ schemat ma rząd $p$ i $x\in C^{{p+1}}$ , to $e_{h}=O(h^{{p+1}})$ zatem

$|E_{n}|\leq e^{{L*(T-t_{0})}}n*e_{h}\leq e^{{L*(T-t_{0})}}\frac{T-t_{0}}{h}O(h^{{p+1}})=O(h^{p}).$

∎

W szczególności zbieżność z rzędem $p$ oznacza dla ustalonego $t\in[t_{0},T]$ i $n*h=t$ , że

$|x_{n}^{h}-x(t)|=O(h^{p})\rightarrow 0\qquad h\rightarrow 0\quad(n\rightarrow\infty).$

Przykład 5.1

Zgodność otwartego schematu Eulera. W zasadzie sprawdzenie konsystentności (zgodności) schematu jest w tym przypadku oczywiste, ponieważ funkcja $\Phi(h,t,x,y)=f(t,x)$ , czyli spełnia założenia zgodności. Dodatkowo wiemy, że schemat ma rząd jeden, więc jeśli $x$ rozwiązanie zadania początkowego (3.1) należy do $C^{2}$ , to $|x_{n}-x(t)|=O(h)$ dla $n*h=t-t_{0}$ , co potwierdzają wyniki eksperymentów z tabeli 4.1 w przypadku skalarnym dla równania $\frac{dx}{dt}=x$ z $x(0)=1$ .

5.2. Teoria zbieżności schematów liniowych wielokrokowych

Teoria zbieżności dla schematów wielokrokowych liniowych różni się od teorii zbieżności dla schematów jednokrokowych. Precyzyjniej pisząc - używamy tak samo jak w przypadku schematów jednokrokowych pojęcie rzędu schematu, ale równie ważne jest pojęcie stabilności schematu. Używając nieformalnego języka - stabilność schematu oznacza to, że błędy z poprzednich kroków się nie kumulują, a wręcz zanikają.

5.2.1. Stabilność, zgodność

Na początku rozważmy dowolny schemat liniowy wielokrokowy (4.8) zastosowany do równania z polem wektorowym równym zero tzn. do zagadnienia początkowego:

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

którego jedynym rozwiązaniem jest rozwiązania stałe $x(t)=1$ . Nasz schemat staje się wtedy równaniem różnicowym liniowym jednorodnym o stałych współczynnikach:

$\alpha _{{k}}x_{{n+k}}+\ldots\alpha _{0}x_{n}=0.$

(5.1)

Wprowadzając wielomian

$\rho(\lambda)=\sum _{{j=0}}^{k}\alpha _{j}\lambda^{j}$

(5.2)

otrzymujemy, że jeśli $\xi\not=0$ jest zerem (pierwiastkiem) wielomianu $\rho(\lambda)$ , to ciąg

$x_{n}=\xi^{n}\qquad n=0,1,2,\ldots,\infty$

jest rozwiązaniem równania różnicowego. Jeśli zero jest dwukrotne , to otrzymujemy nowe rozwiązanie związane z tym pierwiastkiem tzn.:

$y_{n}=n*\xi^{{n-1}},\qquad n=0,1,2,\ldots,\infty.$

i tak dalej - jeśli zero jest trzykrotne, to kolejne rozwiązanie:

$z_{n}=n*(n-1)*\xi^{{n-2}}.$

W każdym razie - co ważne - jeśli jakiś pierwiastek $\rho(\lambda)$ ma moduł $|\xi|>1$ lub $|\xi|=1$ , ale krotność pierwiastka jest większa od jeden, to istnieją rozwiązania nieograniczone, a dokładnie, zachodzi dla pewnych rozwiązań:

$|x_{n}|\rightarrow+\infty\qquad n\rightarrow\infty.$

Jeśli pierwiastek $|\xi|<1$ to zawsze wszystkie rozwiązania z nim związane spełniają

$|x_{n}|\rightarrow 0\qquad n\rightarrow\infty.$

Podsumowując: jeśli jakiś pierwiastek ma moduł mniejszy od jeden, albo ma moduł jeden i krotność jeden, to rozwiązania z nim związane są ograniczone. Wydaje się, że wymaganie żeby wszystkie rozwiązania schematu zastosowanego do równania z prawą stroną równą zero były ograniczone jest jak najbardziej uzasadnione.

Dlatego w ten sposób definiujemy stabilność (ang. stability) schematu liniowego wielokrokowego (4.8):

Definicja 5.2

Schemat (4.8) jest stabilny, jeśli każdy pierwiastek $\xi$ wielomianu $\rho(\lambda)=\sum _{{j=0}}^{k}\alpha _{j}\lambda^{j}$ spełnia

$|\xi|\leq 1,$

a w przypadku jeśli $|\xi|=1$ , to krotność pierwiastka $\xi$ wynosi jeden.

Dodatkowo wprowadzamy pojęcie silnej stabilności schematu (ang. strong stability):

Definicja 5.3

Schemat (4.8) jest silnie stabilny, jeśli spełnia Definicję 5.2 stabilności, i jeśli $\xi$ jest zerem wielomianu $\rho(\lambda)$ takim, że $|\xi|=1,$ to $\xi=1$ .

Pojęcie silnej stabilności nie wpływa na samą teorię zbieżności schematów, ale można się spodziewać, że schematy silnie stabilne zachowują się lepiej, por. rozdział 5.2.2.

W praktycznych obliczeniach wszystkie używane schematy liniowe wielokrokowe są silnie stabilne.

5.2.1.1. Zgodność schematu

Możemy oczekiwać, że jednym z rozwiązań naszego równania różnicowego (5.1) będzie rozwiązanie stałe: $x_{n}=1$ , czyli oczekujemy żeby $\xi=1$ było pierwiastkiem wielomianu $\rho(\lambda)$ , tzn. $\rho(1)=0$ . Ten warunek nazwiemy prezgodnością (ang. preconsistency) schematu.

Dodatkowo rozpatrzmy drugie proste równanie z warunkiem początkowym:

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

którego rozwiązaniem jest $x(t)=t$ . Jeśli zastosujemy schemat do tego równania (zawsze $f_{n}=1$ ), to otrzymujemy następujące liniowe równanie różnicowe o stałych współczynnikach:

$\sum _{{j=0}}^{k}\alpha _{j}x_{{n+j}}^{h}-h*\sum _{{j=0}}^{k}\beta _{j}=0\qquad n\geq 0.$

Dodatkowo wprowadzimy wielomian:

$\sigma(\lambda)=\sum _{{j=0}}^{k}\beta _{j}\lambda^{j}.$

Zakładając, że schemat jest dokładny w $t_{n}=n*h$ dla tego zagadnienia początkowego tzn. wstawiając $x_{n}=x(n*h)=n*h$ otrzymujemy:

$\sum _{{j=0}}^{k}\alpha _{j}(n+j)-\sum _{{j=0}}^{k}\beta _{j}=n*\rho(1)+\rho^{{\prime}}(1)-\sigma(1)=0.$

Jeśli założymy, że schemat jest prezgodny, to powyższe równanie daje nam:

$\rho(1)=0\qquad\rho^{{\prime}}(1)-\sigma(1)=0.$

(5.3)

Wprowadzamy pojęcie zgodności (konsystentności) schematu liniowego wielokrokowego:

Definicja 5.4

Schemat liniowy wielokrokowy (4.8) jest zgodny (konsystentny), (ang.consistent) jeśli zachodzi (5.3).

Nie jest trudno pokazać, że ten warunek jest równoważny temu, że rząd schematu jest co najmniej jeden, tzn. prawdziwe jest następujące stwierdzenie:

Stwierdzenie 5.1

Schemat liniowy wielokrokowy (4.8) jest zgodny wtedy i tylko wtedy, gdy rząd schematu wynosi co najmniej jeden.

Dowód pozostawiamy jako zadanie, por. ćwiczenie 5.14.

5.2.1.2. Zbieżność

Twierdzenie 5.2 (o zbieżności schematów wielokrokowych)

Jeśli rozwiązanie zagadnienia początkowego $x\in C^{{p+1}}([t_{0},T])$ , schemat liniowy $k$ -krokowy jest rzędu $p$ dla $p\geq 1$ i jest stabilny, oraz wartości startowe $x_{j}^{h}$ dla $j=0,\ldots,k-1$ spełniają nierówność:

$\max _{{j=0,\ldots,k-1}}|x(t_{j}^{h})-x_{j}^{h}|\leq Ch^{p}$

dla pewnej stałej nieujemnej $C$ niezależnej od $h$ , to ten schemat jest zbieżny z rzędem $p$ .

Dowód twierdzenia można znaleźć w [5] lub [23].

Przykład 5.2

Zbieżność otwartego schematu Eulera jako schematu liniowego wielokrokowego. Wiemy już, że rząd otwartego schematu Eulera wynosi jeden. Wystarczy więc zbadać stabilność schematu. Wielomian $\rho(\lambda)=\lambda-1$ zatem ma jeden pierwiastek $\lambda _{1}=1$ więc schemat jest stabilny, a nawet silnie stabilny. Jeśli rozwiązanie jest klasy $C^{2}$ , to zbieżność zachodzi z rzędem jeden, co potwierdzają wyniki eksperymentów z tabeli 4.1 w przypadku skalarnym dla równania $\frac{dx}{dt}=x$ z $x(0)=1$ .

5.2.2. Stabilność, a silna stabilność

Rozpatrzmy schemat kroku środkowego (4.3), który jest dwukrokowym schematem liniowym rzędu dwa (zadanie). Nietrudno sprawdzić, że jest to schemat stabilny, ale nie silnie stabilny.

$\$

Rys. 5.1. Schematy midpoint i Eulera otwarty na odcinku [0,1] dla $\frac{dy}{dx}=-y;y(0)=1$ z $h=1e-1$ .

Rozpatrzmy nasze modelowe równanie liniowe:

$\frac{dx}{dt}=-x\quad t>0,\quad x(0)=1,$

którego rozwiązanie $x(t)$ jest równe $\exp(-t)$ . Zastosujmy schemat midpoint biorąc za $x_{1}$ dokładną wartość rozwiązania $x_{1}^{h}=\exp(-h)$ . Narysujmy wykres rozwiązania przybliżonego otrzymanego tym schematem i schematem Eulera otwartym dla $h=0.1$ na odcinku $[0,1]$ , por. rysunek 5.1, i potem na $[0,10]$ , por. rysunek 5.2. W tym drugim przypadku rozwiązanie przybliżone otrzymane schematem midpoint od pewnego momentu przestaje zachowywać się jak rozwiązanie dokładne $\exp(-t)$ , tzn. widzimy coraz większe oscylacje.

$\$

Rys. 5.2. Schematy midpoint i Eulera otwarty na odcinku [0,10] dla $\frac{dy}{dx}=-y;y(0)=1$ z $h=1e-1$ .

$\$

Rys. 5.3. Schematy midpoint i Eulera otwarty na odcinku [0,10] dla $\frac{dy}{dx}=-y;y(0)=1$ z $h=1e-2$ .

$\$

Rys. 5.4. Schematy midpoint i Eulera otwarty na odcinku [0,20] dla $\frac{dy}{dx}=-y;y(0)=1$ z $h=1e-2$ .

Weźmy mniejsze $h=0.01$ . Wtedy na $[0,10]$ wszystko jest w porządku, ale na odcinku $[0,20]$ znów rozwiązanie otrzymane schematem midpoint zaczyna od pewnego momentu oscylować z coraz większą amplitudą, zamiast maleć do zera. Im $h$ jest mniejsze, tym odcinek na którym rozwiązanie otrzymane schematem midpoint dobrze aproksymuje rozwiązanie równania wyjściowego jest większy, ale zawsze w pewnym momencie pojawi się zaburzenie numeryczne, por. rysunki 5.3 i 5.4.

Numeryczne zaburzenie wynika właśnie z tego, że schemat nie jest silnie stabilny. Można wypisać wzór na rozwiązania równania różnicowego zadanego przez schemat midpoint dla tego równania i wykazać, że po jakimś czasie zawsze pojawią się oscylacje.

Wyniki tu otrzymane nie stoją w jakiekolwiek sprzeczności z teorią zbieżności schematów liniowych wielokrokowych. Można sprawdzić, że schemat midpoint spełnia założenia tej teorii dla tego równania, i że na ustalonym odcinku $[0,T]$ zachodzi zbieżność z rzędem dwa. W praktyce nie należy stosować schematów, które nie są silnie stabilne, ponieważ zawsze mogą pojawić się oscylacje po a priori nie znanym czasie, szczególnie kiedy chcemy rozwiązać równanie na długim odcinku czasu.

5.3. Wartości startowe schematów wielokrokowych

W tym miejscu warto zwrócić uwagę, jak w praktycznych obliczeniach implementować schematy liniowe wielokrokowe. Aby przy pomocy danego schematu $k$ -krokowego rzędu $p$ obliczyć $x_{{n+k}}$ , należy znać poprzednie $k$ przybliżenia $x_{{n+k-1}},\ldots,x_{{n}}$ . Czyli aby wystartować schemat musimy znaleźć odpowiednio dobre (tzn. takie, że $E_{j}=O(h^{p})$ , por. Twierdzenie 5.2) startowe przybliżenia $x_{0},x_{1},\ldots,x_{{k-1}}$ . Wartość $x_{0}$ jest zadana jako warunek początkowy, tzn. możemy przyjąć dany warunek początkowy $x_{0}=x(t_{0})$ . Natomiast musimy wyliczyć $x_{1},\ldots,x_{{k-1}}$ . W praktyce do wyliczenia tych wartości możemy np. zastosować $k-1$ razy schemat jednokrokowy tego rzędu co najmniej $p$ .

Dla schematu Adamsa-Bashfortha dwukrokowego rzędu dwa za $x_{1}$ możemy przyjąć $x_{1}$ obliczone jednym krokiem schematu Heuna. Tu warto zauważyć że musimy zastosować tylko jeden krok schematu Heuna. Uwzględniając to widzimy, że $x_{1}$ można by też obliczyć stosując dowolny schemat jednokrokowy rzędu jeden np. otwarty schemat Eulera. Ogólniej do obliczenia startowych wartości $x_{1},\ldots,x_{{k-1}}$ dla schematu $k$ -krokowego rzędu $p$ można zastosować schemat jednokrokowy rzędu $p-1$ lub większego.

5.4. Eksperymentalne badanie rzędu zbieżności schematów

$\begin{array}[]{|c|c|c|c|c|c|c|}\hline h&otwarty\; Euler&|e_{n}/e_{{n+1}}|&Heun&|e_{n}/e_{{n+1}}|&zmod.\; Euler&|e_{n}/e_{{n+1}}|\\ \hline 2.5e-01&4.61e-02&0.00&1.90e-03&0.00&6.57e-03&0.00\\ 1.1e-01&1.98e-02&2.33&3.51e-04&5.40&1.19e-03&5.53\\ 5.3e-02&9.23e-03&2.14&7.63e-05&4.60&2.56e-04&4.64\\ 2.6e-02&4.47e-03&2.07&1.78e-05&4.28&5.97e-05&4.29\\ 1.3e-02&2.20e-03&2.03&4.32e-06&4.13&1.44e-05&4.14\\ 6.3e-03&1.09e-03&2.02&1.06e-06&4.07&3.55e-06&4.07\\ 3.1e-03&5.44e-04&2.01&2.63e-07&4.03&8.79e-07&4.03\\ 1.6e-03&2.71e-04&2.00&6.55e-08&4.02&2.19e-07&4.02\\ \hline\end{array}$

Tabela 5.1. Eksperymentalne badanie rzędu zbieżności schematów: otwartego Eulera, Heuna i zmodyfikowanego schematu Eulera poprzez połowienie kroków dla równania $\frac{dx}{dt}=\cos(x)^{2}$ z $x(0)=0$ dla $t=1$ , którego rozwiązaniem jest $\mathrm{arctan(x)}$ . W kolumnach o indeksach parzystych jest podany błąd odpowiedniego schematu dla $t=1$ , a w kolumnach 3, 5 i 7 są stosunki błędów dla kroku $2*h$ podzielone przez błędy dla kroku $h$ dla odpowiednich schematów.

$\begin{array}[]{|c|c|c|c|c|}\hline h&Heun&|e_{n}/e_{{n+1}}|&rk4&|e_{n}/e_{{n+1}}|\\ \hline 5.3e-01&5.96e+03&&9.09e+01&\\ 2.6e-01&1.91e+03&3.12&6.41e+00&14.18\\ 1.3e-01&5.29e+02&3.61&4.24e-01&15.12\\ 6.3e-02&1.38e+02&3.83&2.73e-02&15.56\\ 3.1e-02&3.52e+01&3.92&1.73e-03&15.78\\ 1.6e-02&8.88e+00&3.96&1.09e-04&15.89\\ 7.8e-03&2.23e+00&3.98&6.81e-06&15.95\\ 3.9e-03&5.59e-01&3.99&4.27e-07&15.97\\ 2.0e-03&1.40e-01&4.00&2.66e-08&16.02\\ \hline\end{array}$

Tabela 5.2. Eksperymentalne badanie rzędu zbieżności schematów: Heuna i Rungego Kutty rzędu cztery poprzez połowienie kroków dla równania $\frac{dx}{dt}=x$ z $x(0)=1$ dla $t=10$ . W kolumnach indeksach parzystych jest podany błąd odpowiedniego schematu dla $t=10$ , a w kolumnach 3, i 5 są stosunki błędów dla kroku $2*h$ podzielone przez błędy dla kroku $h$ dla odpowiednich schematów.

Rząd zbieżności schematu dla czasu $t$ można badać eksperymentalnie dla ustalonego zagadnienia początkowego $\frac{dx}{dt}=f(t,x)$ z warunkiem początkowym $x(0)=x_{0}$ i ze znanym rozwiązaniem $x(t)$ określonym na odcinku $[t_{0},T]$ . Możemy wtedy przy pomocy naszego schematu, np. schematu Heuna, zmodyfikowanego schematu Eulera i otwartego schematu Eulera obliczać przybliżone wartości $x_{n}^{h}$ dla ustalonego $t=t_{n}^{h}=t_{0}+n*h$ , z połowionym krokiem $h$ : $h_{0},h_{0}/2,\ldots,2^{{-k}}h_{0}$ dla ustalonego $h_{0}$ . Będziemy badać jak zmienia się błąd w kolejnych krokach, a dokładniej, jak zmienia się stosunek błędów dla kolejnych połowionych $h$ .

Jeśli błąd dla ustalonego $t$ zachowywałby się jak $O(h^{p})$ , a dokładniej, gdyby $e_{n}=|x_{n}^{h}-x(t)|=c*h^{p}+O(h^{p})$ , to stosunek błędu powinien się zachowywać jak:

$\frac{\| x_{n}^{{2h}}-x(t)\|}{\| x_{{2n}}^{h}-x(t)\|}=\frac{c2^{p}h^{p}+O(h^{{p+1}})}{ch^{p}+O(h^{{p+1}})}\approx 2^{p}\qquad h<1$

dla dostatecznie małych $h$ , czyli dla schematów Heuna i zmodyfikowanego schematu Eulera jak cztery, a dla otwartego schematu Eulera jak dwa, a dla schematu rzędu cztery jak $2^{4}$ - czyli szesnaście.

W tabeli 5.1 widzimy wyniki eksperymentu dla $t=1$ dla schematów: otwartego schematu Eulera, Heuna i zmodyfikowanego schematu Eulera zastosowanych do zagadnienia początkowego $\frac{dx}{dt}=\cos(x)^{2}$ z $x(0)=0$ , którego rozwiązaniem jest $\mathrm{arctan}(x)$ , czyli $x(1)=\mathrm{arctan}(1)$ .

W tabeli 5.2 przedstawiliśmy wyniki tego samego eksperymentu, ale dla schematu Rungego-Kutty rzędu cztery dla zadania początkowego $\frac{dx}{dt}=x$ z $x(0)=1$ i dla $t=10$ , czyli z rozwiązaniem $x(10)=\exp(10)$ . Skoro schemat jest rzędu cztery, możemy oczekiwać, że błąd będzie jak $O(h^{4})$ dla ustalonego $t$ .

Wyniki uzyskane w tabeli 5.2 potwierdzają nasze przypuszczenie. Dla kroku o połowę mniejszego błąd maleje około $2^{4}=16$ razy dla dostatecznie małych $h$ (im $h$ mniejsze, tym ten stosunek bliższy jest szesnastu). Przy okazji zauważmy, jak ogromna jest różnica w błędzie dla schematu rzędu cztery (Rungego-Kutty), a dla schematu Heuna rzędu dwa. W przypadku tego ostatniego - błąd bezwzględny dla $h=2*10^{{-3}}$ jest rzędu $1/10$ , a dla tego pierwszego błąd jest rzędu $10^{{-8}}$ . Z kolei patrząc na błędy względne, w przypadku schematu Heuna błąd jest na poziomie $10^{{-5}}$ , a dla schematu Rungego $10^{{-12}}$ , czyli poziom błędu w arytmetyce podwójnej precyzji praktycznie jest wystarczająco dokładny.

Błędem bezwzględnym nazywamy $e_{n}$ , a względnym $e_{n}/|x(t_{n})|$ . W przypadku gdy $|x(t_{n})|$ jest bardzo duże lub bardzo małe należy rozpatrywać błąd względny choćby z powodu własności arytmetyki zmiennopozycyjnej.

5.5. Zadania

Ćwiczenie 5.1

Zbadaj rząd zbieżności zamkniętego schematu Eulera korzystając z teorii zbieżności dla schematów jednokrokowych.

Ćwiczenie 5.2

Zbadaj rząd zbieżności zamkniętego schematu Eulera korzystając z teorii zbieżności dla schematów wielokrokowych liniowych. Czy schemat jest silnie stabilny?

Ćwiczenie 5.3

Zbadaj rząd zbieżności schematu trapezów dany wzorem (4.7).

korzystając z teorii zbieżności dla schematów jednokrokowych,
korzystając z teorii zbieżności dla schematów wielokrokowych liniowych.

Czy schemat jest silnie stabilny?

Ćwiczenie 5.4

Zbadaj rząd zbieżności zmodyfikowanego Eulera i schematu Heuna.

Ćwiczenie 5.5

Zbadaj rząd zbieżności schematu punktu środkowego (ang. midpoint). Czy schemat jest silnie stabilny?

Ćwiczenie 5.6 (laboratoryjne)

Zaimplementuj schematy: otwarty i zamknięty Eulera w octave i przetestuj rzędy zbieżności eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4 dla równania skalarnego $\frac{dy}{dx}=y^{2};y(0)=1$ oraz dla równania drugiego rzędu $\frac{d^{2}y}{dx^{2}}=-y$ z $y(0)=0,\frac{dy}{dx}(0)=1$ .

Wskazówka:

Do rozwiązywania nieliniowego równania przy implementacji schematu zamkniętego w każdym kroku możesz użyć funkcji fsolve().

Ćwiczenie 5.7 (laboratoryjne)

Zaimplementuj schemat trapezów, por. (4.7), w octave i przetestuj rząd zbieżności eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4 dla równania skalarnego $\frac{dy}{dx}=-y;y(0)=1$ oraz dla równania drugiego rzędu $\frac{d^{2}y}{dx^{2}}=-y$ z $y(0)=0,\frac{dy}{dx}(0)=1$ .

Wskazówka:

Do rozwiązywania nieliniowego równania w każdym kroku możesz użyć funkcji fsolve().

Ćwiczenie 5.8 (laboratoryjne)

Zaimplementuj zmodyfikowany schemat Eulera oraz schemat Heuna w octave i przetestuj rzędy zbieżności eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4.

Ćwiczenie 5.9 (laboratoryjne)

Zaimplementuj schemat Rungego rzędu cztery, dany wzorem (4.15) w octave i przetestuj rząd zbieżności schematu eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4.

Ćwiczenie 5.10 (laboratoryjne)

Zaimplementuj schemat (ang. midpoint) w octave i przetestuj rząd zbieżności schematu eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4, przyjmując, że $x_{1}=y(t_{0}+h)$ , czyli jest równe dokładnemu rozwiązaniu. Dla równania $\frac{dy}{dx}=-y$ z $y(0)=1$ zastosuj schemat na długim odcinku czasu dla ustalonego $h$ . Czy rozwiązanie zachowuje się tak jak tego oczekujemy? Zmniejsz $h$ i zobacz czy sytuacja się poprawia?

Ćwiczenie 5.11

Zbadaj rząd, stabilność i rząd zbieżności otwartego dwukrokowego schematu Adamsa, por. Ćwiczenie 4.7.

Ćwiczenie 5.12 (laboratoryjne)

Zaimplementuj w octave otwarty dwukrokowy schemat Adamsa, por. Ćwiczenie 4.7. Następnie przetestuj eksperymentalnie rząd zbieżności tego schematu metodą połowionego kroku, jak opisano w rozdziale 5.4 dla $\frac{dy}{dt}=-y$ z $y(0)=-1$ . Biorąc za $x_{1}$ :

$x_{1}=y(h)=e^{{-h}}$ czyli rozwiązanie dokładne dla $t=h$ .
$x_{1}$ obliczone schematem Heuna czyli schematem Rungego-Kutty rzędu dwa.
$x_{1}$ obliczone schematem Eulera otwartym czyli schematem rzędu jeden.

Ćwiczenie 5.13

Znajdź dokładne rozwiązanie $(y_{k})_{{k=1}}^{\infty}$ schematu różnicowego, który jest schematem punktu środkowego zastosowanym dla równania $\frac{dy}{dx}=-y$ . Pokaż, że jeśli $y_{0}=1$ , a $y_{1}$ różne od jednej konkretnej wartości (co np. odpowiada warunkom startowym $y_{0}=1,y_{1}=\exp(-h)$ ), to $\lim _{{n\rightarrow\infty}}|y_{n}|=+\infty$ , czyli pojawią się niepotrzebne kumulujące się zaburzenia numeryczne.

Ćwiczenie 5.14 (trudne)

Udowodnij stwierdzenie 5.1.

Wskazówka:

Z tego, że rząd schematu jest co najmniej jeden, wynika od razu warunek zgodności schematu (wystarczy rozważyć oba równania $\frac{dy}{dx}=0$ z $y(0)=1$ i $\frac{dy}{dx}=0$ z $y(0)=0$ ). Natomiast aby pokazać, że ze zgodności schematu wynika to, że rząd schematu jest większy od jeden, należy zastosować wzór Taylora.

Ćwiczenie 5.15

Czy schemat Adamsa (por. rozdział 4.1.2) może nie być stabilny, ewentualnie nie być silnie stabilny? Uzasadnij podając ewentualnie kontrprzykład niestabilnego (czy nie będącego silnie stabilnym) schematu Adamsa.

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

Numeryczne rozwiązywanie równań różniczkowych

Zagadnienia

5. Metody dla równań różniczkowych zwyczajnych - teoria zbieżności

5.1. Teoria zbieżności schematów jednokrokowych

Definicja 5.1

Twierdzenie 5.1 (o zbieżności schematów jednokrokowych)

Przykład 5.1

5.2. Teoria zbieżności schematów liniowych wielokrokowych

5.2.1. Stabilność, zgodność

Definicja 5.2

Definicja 5.3

5.2.1.1. Zgodność schematu

Definicja 5.4

Stwierdzenie 5.1

5.2.1.2. Zbieżność

Twierdzenie 5.2 (o zbieżności schematów wielokrokowych)

Przykład 5.2

5.2.2. Stabilność, a silna stabilność

5.3. Wartości startowe schematów wielokrokowych

5.4. Eksperymentalne badanie rzędu zbieżności schematów

5.5. Zadania

Ćwiczenie 5.1

Ćwiczenie 5.2

Ćwiczenie 5.3

Ćwiczenie 5.4

Ćwiczenie 5.5

Ćwiczenie 5.6 (laboratoryjne)

Ćwiczenie 5.7 (laboratoryjne)

Ćwiczenie 5.8 (laboratoryjne)

Ćwiczenie 5.9 (laboratoryjne)

Ćwiczenie 5.10 (laboratoryjne)

Ćwiczenie 5.11

Ćwiczenie 5.12 (laboratoryjne)

Ćwiczenie 5.13

Ćwiczenie 5.14 (trudne)

Ćwiczenie 5.15