Zagadnienia

6.1 Sztywne równania różniczkowe zwyczajne
- 6.1.1 Przypadek skalarny
- 6.1.2 Przypadek wielowymiarowy
6.2 Przykłady schematów sztywnych
6.3 Schematy zamknięte. Predyktor-korektor
6.4 Adaptacyjny dobór kroku całkowania
6.5 Metoda strzałów
6.6 Zadania

6. Sztywność, zmienny krok całkowania i metoda strzałów

W tym rozdziale zajmiemy się ważnymi schematami rozwiązywania tzw. sztywnych układów równań różniczkowych zwyczajnych. Omówimy schematy ze zmiennym krokiem całkowania i metodę strzałów rozwiązywania zadań brzegowych.

6.1. Sztywne równania różniczkowe zwyczajne

Dość trudno jest podać precyzyjnie poprawną matematycznie definicję sztywności dla dowolnego zadania różniczkowego zwyczajnego. My przyjmiemy definicję pragmatyczną za [13]:

Definicja 6.1

Równanie różniczkowe zwyczajne nazywamy sztywnym (ang. stiff), jeśli numeryczne schematy zamknięte, w szczególności metody zamknięte Adamsa, działają zdecydowanie lepiej niż schematy otwarte przybliżonego rozwiązywania zagadnień początkowych.

Oczywiście definicja nie jest do końca precyzyjna. Podamy też inne definicje sztywności (ang. stiffness) np. dla zagadnień liniowych, jakkolwiek powyższa definicja jest dla nas wygodna, ponieważ podkreśla to, że równania sztywne rozwiązujemy przy pomocy schematów zamkniętych. Za chwilę podamy kilka przykładów sztywnych równań różniczkowych, aby przekonać się, że pojawiają się one dość często w realistycznych modelach nauk przyrodniczych, por. rozdział 6.2.

6.1.1. Przypadek skalarny

W rozdziale 3.3.1 już zobaczyliśmy, że dla modelowego zadania skalarnego:

$\frac{dx}{dt}=a*x\quad x(0)=1\qquad a<0.$

rozwiązanie uzyskane przy pomocy otwartego schematu Eulera zachowuje własności rozwiązania $x(t)=\exp(a*t)$ , tzn. jest dodatnie i malejące do zera dla $t\rightarrow+\infty$ tylko wtedy, gdy jest spełniony warunek: $h<-1/a.$ W przypadku gdy $|a|$ jest bardzo duże, warunek ten wymusza, że musimy stosować bardzo małe $h$ . Natomiast rozwiązanie uzyskane przy pomocy zamkniętego schematu Eulera zachowuje własności powyższe rozwiązania dla dowolnego $h>0$ .

Schemat zamknięty Eulera można zatem uznać za lepszy od schematu otwartego dla tego zagadnienia dla ujemnego $a$ o dużym module.

6.1.2. Przypadek wielowymiarowy

Załóżmy, że rozpatrujemy jednorodne liniowe równanie różniczkowe zwyczajne ze stałymi współczynnikami:

$\frac{dx}{dt}=A*x\qquad x(t_{0})=\vec{\beta},$

gdzie $A$ - to macierz o stałych współczynnikach taka, że w pewnej bazie jest ona diagonalizowalna, tzn. istnieje macierz nieosobliwa $C$ taka, że

$A=C\Lambda C^{{-1}}$

$\Lambda=\begin{pmatrix}\lambda _{1}&&&\\ &\lambda _{2}&&\\ &&\ddots&\\ &&&\lambda _{N}\end{pmatrix}$

Załóżmy, że wszystkie $\lambda _{k}<0$ , wtedy oczywiście rozwiązaniem jest

$x(t)=Ce^{{\Lambda(t-t_{0})}}C^{{-1}}\vec{\beta}.$

Oznaczmy $k$ -tą kolumnę macierzy $C$ przez $\vec{c}_{k}$ , tzn. $C=(\vec{c}_{1},\ldots,\vec{c}_{N})$ , wtedy otrzymujemy:

$x(t)=\sum _{k}\alpha _{k}\vec{c}_{k}e^{{\lambda _{k}(t-t_{0})}}.$

dla $\vec{\alpha}=(\alpha _{1},\ldots,\alpha _{N})^{T}=C^{{-1}}\vec{\beta}$ . Jeśli zastosujemy otwarty schemat Eulera do tego zagadnienia, to analogicznie jak w poprzednim rozdziale otrzymujemy, że

$\displaystyle x_{n}$

$\displaystyle=$

$\displaystyle(I+h*A)x_{{n-1}}=C(I+h*\Lambda)C^{{-1}}x_{{n-1}}=C(I+h*\Lambda)^{n}\vec{\alpha}=\sum _{k}\alpha _{k}\vec{c}_{k}(1+h\lambda _{k})^{n}.$

Czyli: o ile $|\lambda _{k}|\gg 1$ ( $\lambda _{k}<0$ ) tym odpowiednia składowa rozwiązania $\vec{c}_{k}e^{{\lambda _{k}(t-t_{0})}}$ szybciej dąży do zera. Z drugiej strony warunek na to, aby odpowiadająca składowa rozwiązania dyskretnego otrzymanego otwartym schematem Eulera nie zmieniała znaku i zbiegała do zera wynosi:

$h<1/|\lambda _{k}|,$

czyli jest to warunek ograniczający dopuszczalny zakres wartości $h$ .

Widzimy zatem, że otwarty schemat Eulera dla takiego równania jest zupełnie niepraktyczny na dłuższych odcinkach czasu, ponieważ ograniczenie na $h$ związane jest ze składowymi rozwiązań, które najszybciej zanikają, czyli na dłuższym odcinku czasu nie mają większego wpływu na rozwiązanie. Z kolei dla schematu zamkniętego Eulera nie otrzymujemy żadnych ograniczeń na $h$ , ponieważ:

$x_{n}=C(I-h*\Lambda)^{{-n}}\vec{\alpha}=\sum _{k}\alpha _{k}(1-h\lambda _{k})^{{-n}}.$

Widzimy, że dla tego typu równań schemat otwarty zachowuje się gorzej niż odpowiedni schemat zamknięty, co jest zgodne z naszą oryginalną definicją sztywności. W ogólnym przypadku, gdy wszystkie części rzeczywiste wartości własnych macierzy $A$ są ujemne sytuacja jest analogiczna.

Zatem definiujemy zadanie liniowe jednorodne jako sztywne, jeśli:

$Re\,\sigma(A)\subset\{ x:x<0\}$
$\frac{\max _{{\lambda _{k}\in\sigma(A)}}|Re\,\lambda _{k}|}{\min _{{\lambda _{k}\in\sigma(A)}}|Re\,\lambda _{k}|}$ jest duże.

Tutaj $\sigma(A)$ oznacza zbiór wartości własnych macierzy $A$ . Oczywiście w tej definicji nie jest doprecyzowane co oznacza <<duże>>, ale łatwo określić, że jeśli stosunek w drugim punkcie jest równy dziesięć, to układ nie jest sztywny, a jeśli $10^{{20}}$ to układ jest sztywny. Rozszerza się powyższą definicję sztywności na układy równań nieliniowych przyjmując, że układ:

$\frac{dx}{dt}=F(t,x)$

jest sztywny w obszarze $G$ i dla $t$ z odcinka $(a,b)$ jeśli Jakobian $F$ , czyli $D_{x}F(t,x)$ spełnia powyższą definicję dla każdego $x\in G$ i $t\in(a,b)$ .

Wadą powyższej definicji sztywności jest to, że nie obejmuje np. układu skalarnego $\stackrel{.}{x}=a*x$ dla $a<0$ .

6.2. Przykłady schematów sztywnych

W tym rozdziale podamy kilka przykładów równań sztywnych, por. [13].

6.2.1. Oscylator Van der Pola

Równanie Van der Pola opisujące oscylator z nieliniowym tłumieniem:

$\frac{d^{2}y}{dx^{2}}-a*(1-y*y)*\frac{dy}{dx}+y=0\qquad a>0,$

gdzie $a>0$ jest parametrem, dla dużego $a>1$ np. $a=1000$ powyższe równanie jest sztywne.

6.2.2. Reakcje chemiczne

Rozważmy następujące reakcje chemiczne, które symbolicznie opiszemy następująco:

$\displaystyle A$	$\displaystyle\stackrel{0.04}{\rightarrow}$	$\displaystyle B\qquad\mathrm{(wolna)}$
$\displaystyle B+B$	$\displaystyle\stackrel{3*10^{7}}{\rightarrow}$	$\displaystyle C+B\qquad\mathrm{(bardzo\; szybka)}$
$\displaystyle B+C$	$\displaystyle\stackrel{10^{4}}{\rightarrow}$	$\displaystyle A+C\qquad\mathrm{(szybka)}$

co prowadzi do następującego układu równań różniczkowych zwyczajnych:

$\begin{array}[]{lccrrrrrr}A:&x_{1}^{{\prime}}&=&-0.04*x_{1}&+&10^{4}x_{2}*x_{3}&&\qquad\ x_{1}(0)=1\\ B:&x_{2}^{{\prime}}&=&0.04*x_{1}&-&10^{4}x_{2}*x_{3}&-3*10^{7}x_{2}^{2}&x_{2}(0)=1\\ C:&x_{3}^{{\prime}}&=&&&&3*10^{7}x_{2}^{2}&x_{3}(0)=1\end{array}$

Czytelnikowi pozostawiamy sprawdzenie z pomocą octave'a, że np. schematy otwarte nie działają najlepiej dla tego problemu.

6.2.3. Równania paraboliczne

Ten przykład jest szczególnym przypadkiem dyskretyzacji równań, których metody dyskretyzacji omawiane są później dokładniej w rozdziale 14.

Rozpatrzmy równanie paraboliczne:

$\frac{\partial u}{\partial t}(t,x)=\frac{\partial^{2}u}{\partial x^{2}}(t,x)+f(t,x)\qquad x\in(0,1)\quad t\in(0,T]$

z warunkami brzegowymi $u(t,0)=u(t,1)=0$ i początkowym $u(0,x)=u_{0}(x)$ .

Dyskretyzując je względem zmiennej przestrzennej $x$ metodą różnic skończonych, tzn. wprowadzając siatkę $x_{k}=k*h$ dla $k=0,\ldots,N$ z $h=1/N$ i przybliżając drugą pochodną przez

$\frac{\partial^{2}u}{\partial x^{2}}(t,x)\approx\frac{u(t,x-h)-2*u(t,x)+u(t,x+h)}{h^{2}},$

otrzymujemy następujący układ równań różniczkowych zwyczajnych:

$\frac{du_{k}}{dt}(t)=\frac{u_{{k-1}}(t)-2*u_{k}(t)+u_{{k+1}}(t)}{h^{2}}+f_{k}(t)\qquad k=1,\ldots,N-1$

z $u_{0}(t)=u_{N}(t)=0$ (warunki brzegowe) i warunkiem początkowym $u_{k}(0)=u_{0}(k*h)$ i $f_{k}(t)=f(k*h,t)$ , por. rozdział 14.

Oczekujemy, że $u_{k}(t)\approx u(t,x_{k})$ gdzie $u(t,x)$ jest rozwiązaniem wyjściowego problemu. Nietrudno zauważyć, że powyższy układ równań zwyczajnych można zapisać jako

$\vec{u}^{{\prime}}=-h^{{-2}}A_{N}\vec{u}+\vec{f},\qquad\vec{u}(0)=(u_{0}(x_{1}),\ldots,u_{0}(x_{{N-1}}))^{T}$

dla $\vec{u}=(u_{1}(t),\ldots,u_{{N-1}}(t))^{T}$ , $\vec{f}(t)=(f_{1}(t),\ldots,f_{{N-1}}(t))^{T}$ i

$A_{N}=\begin{pmatrix}2&-1&0&&\\ -1&2&-1&\\ &\ddots&\ddots&\ddots&&\\ &&-1&2&-1\\ &&0&-1&2\end{pmatrix}$

(6.1)

Można pokazać, że wartości własne $A_{N}$ są dodatnie i

$\frac{\max _{{\lambda\in\sigma(A_{N})}}\lambda}{\min _{{\lambda\in\sigma(A_{N})}}\lambda}=O(N^{2}),$

czyli układ jest sztywny dla dużych $N$ , czyli małych $h$ .

Dla przykładowych dużych $N$ możemy sprawdzić, czy rzeczywiście tak jest w octave, że stosunek największej do najmniejszej wartości własnej wynosi ok. $4000$ dla $N=100$ , a dla $N=10^{3}$ ten stosunek wynosi ok. $4*10^{5}$ . Tu podajemy kod octave obliczający ten stosunek:

6.3. Schematy zamknięte. Predyktor-korektor

Schematy zamknięte stosujemy dla zadań sztywnych.

Aby obliczyć kolejne przybliżenie $x_{n}$ w schematach zamkniętych, musimy rozwiązać liniowy lub nieliniowy układ równań np. dla zamkniętego schematu Eulera:

$x_{n}=x_{{n-1}}+h*f_{n},$

czyli w każdym kroku musimy rozwiązać układ równań:

$g(x_{n})=0$

dla funkcji $g(x)=\textcolor{red}{x-h*f(t_{n},x)}-x_{{n-1}}$ \textcolor.

W ogólności dla zamkniętego schematu $k$ -krokowego liniowego musimy w każdym kroku rozwiązać układ równań względem $x_{n}$ :

$0=g(x_{n})=\textcolor{red}{\alpha _{k}x_{n}-h*\beta _{k}f(t_{n},x_{n})}+\sum _{{j=0}}^{{k-1}}\alpha _{j}x_{{n+j-k}}-h\sum _{{j=0}}^{{k-1}}\beta _{j}f_{{n+j-k}}.$

Analogiczna sytuację widzimy dla zamkniętych schematów jednokrokowych. Powyższe równanie (układ równań) możemy rozwiązać przy pomocy różnych metod np. metody Newtona, czy jakiejś wersji metody iteracji prostej, zob. np. [18], [17]. Zauważmy, że w przypadku otwartego schematu Eulera $x_{n}$ jest punktem stałym dla funkcji $G_{n}(x)=h*f(t_{n},x)-x_{{n-1}}$ , czyli naturalne jest zastosowanie następującej metody iteracyjnej: dla danego $x^{0}$ liczymy

$x^{k}=h*f(t_{n},x^{k})+x_{{n-1}}=G_{n}(x^{k})\qquad k=1,\ldots.$

Z odpowiedniej gładkości pola wektorowego $f$ wynika, że $x\mapsto f(t_{n},x)$ jest funkcją lipschizowską względem $x$ (lokalnie na $\overline{K}=\overline{K}(x_{{n-1}},\epsilon)$ ) ze stałą Lipschitza $L_{f}$ i $f$ jest ograniczona na kuli $\overline{K}$ , tzn. dla $x\in\overline{K}$ zachodzi $\| f(t_{n},x)\|\leq M$ . Wtedy

$\forall x\in\overline{K}\qquad\| G_{n}(x)-x_{{n-1}}\|\leq h*\| f(t_{n},x)\|\leq h*M,$

i stała Lipschitza $G_{n}$ wynosi $h*L_{f}$ , ponieważ

$\forall x,y\in\overline{K}\qquad\| G_{n}(x)-G_{n}(y)\|=\| h*f(t_{n},x)-h*f(t_{n},y)\|\leq h*L_{f}\| x-y\|.$

Zatem jeśli dla odpowiednio małego $h$ zachodzi $h*M<\epsilon$ i $h*L_{f}\leq\alpha<1$ , to $G_{n}$ jest kontrakcją na $\overline{K}$ . Z tego wynika istnienie $x_{n}$ i zbieżność metody iteracji prostej. W przypadku innych schematów zamkniętych możemy skonstruować analogiczne wersje metody iteracji prostych.

Postawmy kwestię - jak dobierać startowe przybliżenie $x^{0}$ .

Pierwsza opcja to: wziąć $x^{0}=x_{{n-1}}$ . Z ciągłości rozwiązania wynika, że jeśli $h$ jest dostatecznie małe, to $x_{n}\approx x(t_{n})$ i $x_{{n-1}}\approx x(t_{n}-h)$ są sobie bliskie, a dokładnie zachodzi $\| x(t_{n}-h)-x_{n}\|\approx O(h)$ .

Zastanówmy się, czy można dobrać $x^{0}$ lepiej?

Istnieje możliwość, żeby za $x^{0}$ brać przybliżenie $x(t_{n})$ obliczone jednym krokiem otwartego schematu tego samego rzędu $p$ co schemat zamknięty (oczywiście zbieżnym z tym samym rzędem) tzn.

$x^{0}=\hat{x}_{n}=\Phi(h,t_{n},x_{{n-l}}\ldots,x_{{n-1}})$

gdzie $x_{{n-1}},\ldots,x_{{n-l}}$ są obliczone wczesniej naszym schematem zamkniętym, a $\hat{x}_{n}=\Phi(h,t_{n},\hat{x}_{{n-l}}\ldots,\hat{x}_{{n-1}})$ jest dowolnym $l$ - krokowym otwartym schematem zbieżnym z rzędem $p$ , por. (4.5).

Wtedy widzimy, że $\| x^{0}-x(t_{n})\|\approx O(h^{p})$ więc i o ile $h$ dostatecznie małe $\| x^{0}-x_{n}\|\approx O(h^{p})$ .

W takim przypadku schemat otwarty nazywamy predyktorem, a schemat zamknięty, który naprawdę stosujemy do rozwiązania zadania początkowego - korektorem. Podsumowując; nazwy schemat predyktor-korektor używa się względem schematu zamkniętego rzędu $p$ , zaimplementowanego w ten sposób, że kolejne $x_{n}^{h}$ przybliżenie $x(t_{n}^{h})$ obliczone jest poprzez zastosowanie w każdym kroku czasowym jakiejś metody iteracyjnej rozwiązywania nieliniowego równania (układu równań) z przybliżeniem startowym obliczonym odpowiednim pojedyńczym krokiem schematu otwartego tego samego rzędu (predyktorem). Metoda iteracyjna niekoniecznie musi być taka, jak opisana powyżej. Do rozwiązywania nieliniowego układu równań można stosować też np. metodę Newtona, czy jeszcze inną metodę iteracyjną, por. np. [18] lub [25].

W praktyce bierze się odpowiednie pary schematów tego samego rzędu: np. otwarty schemat Eulera za predyktor i zamknięty schemat Eulera za korektor, czy ogólniej - schemat otwarty Adamsa-Bashfordsa rzędu $k$ za predyktor ze schematem zamkniętym Adamsa-Moultona rzędu $k$ jako korektorem. Popatrzmy, jak wygląda przykładowa implementacja schematu predyktor-korektor w przypadku schematów Eulera: otwartego schematu Eulera wziętego jako predyktor i zamkniętego schematu Eulera, który tu pełni rolę korektora dla równania $\stackrel{.}{x}=1+x*x$ z $x(0)=0$ z rozwiązaniem $x(t)=\tan(t)$ . Zaimplementowaliśmy powyższy schemat w octave biorąc jako metodę iteracyjnego rozwiązywania równania nieliniowego w każdym kroku funkcję octave fsolve():

 function [X,t]=predkoreuler(f,t0=0,x0=1,N=100,h=1.0/N)
 # Parametry funkcji:
 #f - wskaznik do pola wektorowego - funkcji dwóch argumentów f(x,t)
 #    przy czym x0 - wektor pionowy dlugosc M;
 # przyklad definicji wskaznika do prostego pola wekt.: f=@(x,t) -x;
 #t0 - czas początkowy
 #h - stały krok dla schematu Eulera
 #N - ilość kroku schematu
 #Funkcja zwraca macierz X wymiaru (N+1)xM długości N+1 taka ze
 #X(k,:) jest przybliżeniem rozwiazania w punkcie czasu t0+(k-1)*h
 #oraz wektor t dlugosci N+1 z dyskretnymi punktami czasowymi
 global xx hh tt
 hh=h;
 M=length(x0);
 X=zeros(N+1,M);
 t=zeros(N+1,1);
 xx=X(1,:)=x0;
 tt=t(1)=t0;
 for k=2:N+1,
   xp=xx+h*f(xx,tt); #predyktor
   g=@(x) x - hh*f(x,tt) - xx; #funkcja pomocnicza dla zamkniętego schematu Eulera
   X(k,:)=xx=fsolve(g,xp); #rozwiązujemy równanie - korektor
   tt+=hh;
   t(k)=tt;
 endfor
 endfunction

6.4. Adaptacyjny dobór kroku całkowania

Stałe w twierdzeniach o zbieżności schematów są znacznie zawyżone i dobór kroku całkowania w oparciu o szacowania z tych twierdzeń jest niepraktyczny. Czy można jakoś oszacować błąd na bieżąco i zmieniać krok całkowania w zależności od tych oszacowań?

Załóżmy, że chcemy użyć konkretnego schematu jednokrokowego rzędu $k$ przybliżonego rozwiązywania zadania początkowego (3.1) takiego, że przy założeniu odpowiedniej gładkości pola wektorowego $f$ otrzymujemy, że błąd metody spełnia dla $0<h<1$ :

$e_{n}^{h}=x_{n}^{h}-x(t_{n}^{h})=e(t_{n}^{h})h^{k}+O(h^{{k+1}}).$

dla $t_{n}^{h}=t_{0}+n*h\in[a,b]$ , $x$ rozwiązania (3.1), $x_{n}^{h}$ rozwiązania przybliżonego obliczonego naszym schematem i pewnej funkcji $e(t)$ . Można pokazać, że tak rzeczywiście jest, i że funkcja $e(t)$ jest rozwiązaniem odpowiedniego równania różniczkowego. Najważniejsze zaś jest to, że $e(t)$ nie zależy od $h$ . Oznaczmy przez $x(t;h):=x_{n}^{h}$ rozwiązanie otrzymane przy pomocy schematu dla ustalonego $h$ i dla $t=t_{n}^{h}$ , a przez $e(t;h):=e_{n}^{h}$ oznaczmy błąd metody w punkcie $t=t_{n}^{h}=t_{0}+n*h$ .

Wtedy

	$\displaystyle e(t;h)=e(t)h^{k}+O(h^{{k+1}}),$
	$\displaystyle e(t;h/2)=e(t)\frac{h^{k}}{2^{k}}+O(h^{{k+1}}).$

Postępując podobnie jak w ekstrapolacji Richardsona, jeśli odejmiemy stronami te równości to otrzymamy:

$\displaystyle e(t;h)-e(t;h/2)=x(t;h)-x(t;h/2)=e(t)\frac{h^{k}}{2^{k}}(2^{k}-1)+O(h^{{k+1}}),$

czyli otrzymujemy:

$E(x;h/2):=\frac{x(t;h)-x(t;h/2)}{2^{k}-1}=e(t)\frac{h^{k}}{2^{k}}+O(h^{{k+1}})\approx e(t;h/2).$

dla $h\ll 1$ , a zatem:

$\| e(t)\|\approx\| E(x;h/2)\|*\frac{2^{k}}{h^{k}}.$

Otrzymaliśmy w ten sposób estymator błędu.

Jeśli chcemy otrzymać błąd na poziomie $\epsilon$ i dla pewnego $t$ obliczyliśmy:

$\| E(x;h/2)\|\approx\| e(t)\|\frac{h^{k}}{2^{k}}\approx\| e(t;h/2)\|,$

to możemy wyliczyć $h_{1}$ , dla którego błąd będzie na poziomie $\epsilon$ , tzn. przyjmując

$\| e(t;h_{1})\|\approx\| e(t)\| h_{1}^{k}=\epsilon$

otrzymujemy

$h_{1}=\left(\frac{\epsilon}{\| e(t)\|}\right)^{{1/k}}\approx\frac{h}{2}\left(\frac{\epsilon}{\| E(x;h/2)\|}\right)^{{1/k}}.$

(6.2)

Następnie możemy zastosować ten schemat z krokiem $h_{1}$ .

Oczywiście adapatacyjną zmianę kroku całkowania (ang. adaptive step control) można stosować i do zwiększania kroku w celu obniżania kosztu obliczeń.

To znaczy, że jeśli $\| E(x;h/2)\|>\epsilon$ , to możemy zmniejszyć krok zgodnie z powyższym wzorem i wtedy powtarzamy obliczenia z mniejszym krokiem $h_{1}$ . Jeśli $\| E(x;h/2)\|\leq\epsilon$ to za przybliżenie $x(t)$ weźmiemy $x(t;h/2)$ , a do następnego kroku możemy przyjąć nowy większy krok $h_{1}$ z (6.2).

Oczywiście zamiast połowienia kroku możemy obliczać $x(t;h/q)$ dla $q=3$ lub $4$ i wtedy otrzymujemy analogiczne wzory.

Można też, zamiast stosowania tego samego schematu dwa razy z krokiem $h$ i potem $h/2$ , obliczać przybliżenie $x(t)$ , schematem rzędu $k$ , a potem większego rzędu np. $k+1$ , jak to się dzieje np. w metodzie Rungego-Fehlberga, gdzie stosuje się schematy Rungego-Kutty czwartego rzędu i Rungego-Kutty piątego rzędu, por. rozdział 17.2 w [24].

6.5. Metoda strzałów

Metoda strzałów (ang. shooting method) służy rozwiązywaniu zadań brzegowych. Rozpatrujemy w tym przypadku równanie różniczkowe zwyczajne, w którym część warunków początkowych zastępujemy liniowymi lub nieliniowymi warunkami brzegowymi, tzn. szukamy funkcji klasy $C^{1}$ na odcinku $[a,b]$ spełniającej:

	$\displaystyle\frac{dx}{dt}=f(t,x),\qquad t\in(a,b)$				(6.3)
	$\displaystyle g(x(a),x(b))=0$

dla $g:U\subset\mathbb{R}^{m}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{{2m}}$ danej funkcji co najmniej ciągłej.

Proszę zauważyć, że ogólnie takie zadanie nie musi mieć rozwiązania nawet w prostym przypadku np.

$\frac{d^{2}x}{dt^{2}}=-x\qquad x(0)=0\quad x(\pi)=1.$

Rozwiązanie ogólne tego równania to $c_{1}\sin(t)+c_{2}\cos(t)$ i z powyższych warunków brzegowych otrzymujemy sprzeczne warunki na $c_{2}$ : $c_{1}\sin(0)+c_{2}\cos(0)=c_{2}=0$ i $c_{1}\sin(\pi)+c_{2}\cos(\pi)=-c_{2}=1$ .

Jeśli istnieje rozwiązanie zadania brzegowego (6.3), to oczywiście jest to szczególny przypadek rozwiązania zadania początkowego (dla pewnej wartości $s$ ):

	$\displaystyle\frac{d^{2}x}{dt^{2}}$	$\displaystyle=$	$\displaystyle f(t,x)\qquad t\in(a,b),$		(6.4)
	$\displaystyle x(a)$	$\displaystyle=$	$\displaystyle s.$

Dodatkowo wiemy, że jeśli $f$ jest funkcją ciągłą, to wartość rozwiązania powyższego zadania początkowego dla $t=b$ , tzn. $x(b;s)$ jest funkcją ciągłą względem $s$ . A jeśli $f$ jest klasy $C^{k}$ , to $x(b;s)$ ma taką samą gładkość jak $f$ , por. np. [23].

Jeśli istnieje rozwiązanie (6.3), to dla pewnego $s^{*}$ zachodzi $g(s,x(b;s^{*}))=0$ . Sprowadziliśmy zadanie brzegowe do zadania nieliniowego znalezienia pierwiastka układu:

$F(s):=g(s;x(b;s))=0.$

$\$

Rys. 6.1. Metoda strzałów - przybliżone rozwiązanie zadania brzegowego: $\frac{d^{2}y}{dx^{2}}=\sin(y)$ z $y(0)=1$ i $y(1)=2$ .

$\$

Rys. 6.2. Kilka strzałów tzn. wykresów przybliżonych rozwiązań zadania początkowego: $\frac{d^{2}y_{k}}{dx^{2}}=\sin(y_{k})$ z $y_{k}(0)=1$ i $\frac{dy_{k}}{dx}(0)=s_{k}$ dla kolejnych iteracji $s_{k}$ .

Do rozwiązania tego układu możemy zastosować jakąś metodę rozwiązywania układów równań nieliniowych, np. metodę bisekcji (o ile zadanie jest skalarne), czy metodę Newtona lub iteracji prostych, por. np. [18].

Można się spytać: jak obliczyć wartość $F(s)$ dla danego $s$ . Trzeba obliczyć $x(b;s)$ , które jest wartością rozwiązania zadania początkowego (6.4) dla $t=b$ z warunkiem początkowym $x(a)=s$ .

Zwykle nie znamy rozwiązań ogólnych tego równania, więc musimy zastosować jakiś schemat rozwiązywania zadania początkowego.

Dla przykładu zastosowaliśmy metodę strzałów do rozwiązania zadania:

$\frac{d^{2}y}{dx^{2}}=\sin(y)\qquad y(0)=1\quad y(1)=2.$

Wykorzystaliśmy metodę rozwiązywania równań zwyczajnych w octave lsode(), w połączeniu z funkcją octave'a rozwiązywania równań nieliniowych fsolve().

Otrzymaliśmy, że dla $s=0.53595$ błąd wynosi w przybliżeniu $10^{{-8}}$ . Na rysunku 6.1 widzimy wykres rozwiązania, a na rysunku 6.2 widać wykresy przybliżeń rozwiązania, tzn. rozwiązania zadania początkowego z $\stackrel{.}{x}=s_{k}$ dla $s_{k}$ wartości kolejnych iteracji metody Newtona. Na rysunku 6.3 widzimy wykres funkcji $F(s)$ .

$\$

Rys. 6.3. Wykres funkcji $s\mapsto y(1;s)$ dla $y$ rozwiązania zadania początkowego: $\frac{d^{2}y}{dx^{2}}=\sin(y)$ z $y(0)=1$ i $\frac{dy}{dx}(0)=s$ .

Niestety metoda strzałów w wielu przypadkach może być bardzo niestabilna. Rozpatrzmy bardzo proste liniowe zadanie: $\frac{d^{2}y}{dx^{2}}+y(x)=0$ z warunkami brzegowymi $y(0)=y(20)=1$ , dla którego znamy rozwiązanie $y(t)=(e^{{t-20}}+e^{{-t}})/(1+e^{{-20}})$ . Zastosowanie metody strzałów z wykorzystaniem standardowej metody rozwiązywania równań zwyczajnych octave'a, czyli funkcją lsode(), daje rozwiązanie przybliżone, dla którego błąd w $t=20$ wynosi ok. $190$ . Wynika to z tego, że wartość rozwiązania zadania początkowego $\frac{d^{2}y}{dx^{2}}+y(x)=0$ $y(0)=1$ $y^{{\prime}}(0)=s$ dla $t=20$ , tzn. $y(20;s)$ , jest bardzo niestabilna. Małe zaburzenie $s$ powoduje ogromną zmianę wyniku. W tym przypadku możemy funkcję $y(t;s)$ , wyliczyć analitycznie, co pozostawiamy jako zadanie. Natomiast zastosowanie metody różnic skończonych daje dobre wyniki, por. rozdział 7.

6.6. Zadania

Ćwiczenie 6.1

Rozpatrzmy następujące zadanie brzegowe:

$\frac{d^{2}x}{dt^{2}}(t)-a(t)x(t)=f(t),\quad x(0)=\alpha,\quad x(b)=\beta.$

dla $f$ funkcji $C^{\infty}$ .

Pokaż, że to zadanie ma jednoznaczne rozwiązanie dla stałego współczynnika $a=Const\geq 0$ .
Przy założeniu, że współczynnik $a$ jest stały, wyznacz wszystkie wartości $a$ , dla których powyższe zadanie może nie mieć rozwiązania.
Pokaż, że jeśli znamy rozwiązania zadania początkowego $x_{1}$ i $x_{2}$ :

$\frac{d^{2}x}{dt^{2}}(t)-a(t)x(t)=f(t)\quad x(0)=\alpha\quad x^{{\prime}}(0)=s_{k}$

dla różnych $s_{1}\not=s_{2}$ takie, że $x_{1}(b)\not=x_{2}(b)$ , to możemy wyznaczyć wzór na $s$ od $s_{1},s_{2},x_{1}(b),x_{2}(b)$ takie, że rozwiązanie zadania początkowego dla tego równania z warunkiem początkowym $x(0)=\alpha\quad x^{{\prime}}(b)=s$ będzie rozwiązaniem wyjściowego zadania brzegowego.

Ćwiczenie 6.2 (laboratoryjne)

Rozpatrzmy następujące zadanie brzegowe:

$\frac{d^{2}x}{dt^{2}}(t)-a(t)x(t)=f(t),\quad x(0)=\alpha,\quad x(b)=\beta$

dla $a(t)\geq 0$ .

Zaimplementuj w octave metodę rozwiązywania zadania brzegowego metodą strzałów korzystając z rozwiązania poprzedniego zadania. W szczególności przetestuj dla $a(t)=1$ i $f(t)=0$ z warunkami brzegowymi $x(0)=x(b)=1$ dla $b=1,10,20$ . Porównaj wynik z rozwiązaniem dokładnym $x(t)=(e^{{t-b}}+e^{{-t}})/(1+e^{{-b}})$ .

Wskazówka:

Rozwiąż na odcinku $[0,b]$ korzystając z funkcji octave lsode() zadanie początkowe dla tego równania z dwoma różnymi warunkami początkowym $x_{k}(0)=\alpha$ i $x_{k}^{{\prime}}(0)=s_{k}$ dla $k=1,2$ i $s_{1}\not=s_{2}$ . Następnie oblicz $s$ takie, że rozwiązanie zadania początkowego z $x(0)=\alpha$ i $x^{{\prime}}(0)=s$ będzie rozwiązaniem wyjściowego zadania brzegowego.

Ćwiczenie 6.3 (laboratoryjne)

Bazując na otwartym schemacie Eulera zaimplementuj w octave schemat z adaptacyjnym krokiem całkowania korzystający ze wzoru (6.2) w rozdziale 6.4. Następnie dla równania $y^{{\prime}}=-y$ z $y(0)=1$ sprawdź błąd tego schematu dla $t=1$ i $t=100$ .

Ćwiczenie 6.4 (częściowo laboratoryjne)

Wyprowadź wzór analogiczny do (6.2) w rozdziale 6.4 dla schematu drugiego rzędu obliczając $x(t;h/q)$ dla $q=3$ , tzn. wzór na oszacowanie błędu bazujący na przybliżeniach rozwiązania otrzymanych danym schematem dla $h$ i $h/3$ . Zastosuj otrzymane wzory dla zadania początkowego z poprzedniego zadania i schematu Heuna.

Ćwiczenie 6.5

Udowodnij, że macierz $A_{N}$ dana wzorem (6.1) jest symetryczna, nieosobliwa i nieujemnie określona, czyli dodatnio określona.

Wskazówka:

Nieujemną określoność najprościej udowodnić z twierdzenia Gerszgorina, por. [17]. A nieosobliwość macierzy - wprost zakładając, że istnieje niezerowy wektor w jądrze macierzy i dochodząc do sprzeczności.

Ćwiczenie 6.6 (laboratoryjne)

Zaimplementuj w octave schemat predyktor-korektor biorąc za korektor schemat trapezów rzędu dwa, a za predyktor schemat Heuna. Do rozwiązywania nieliniowego układu równań zastosuj funkcję octave'a fsolve(). Przetestuj rząd takiego schematu metodą połowionego kroku, jak opisano w rozdziale 5.4, dla równania $\stackrel{.}{x}=1+x*x$ z $x(0)=0$ dla $t=1$ z rozwiązaniem $x(t)=\tan(t)$ i dla równania wahadła porównując z rozwiązaniem otrzymanym dla równania wahadła przy pomocy funkcji octave'a lsode().

Ćwiczenie 6.7 (laboratoryjne)

Zaimplementuj w octave schemat predyktor-korektor biorąc za korektor zamknięty schemat Eulera, a za predyktor otwarty schemat Eulera. Do rozwiązywania nieliniowego układu równań zastosuj swoje metody rozwiązywania równań nieliniowych tzn.:

metodę iteracji prostych, jak opisano w rozdziale 6.3,
wielowymiarową metodę Newtona.

Przetestuj rząd takiego schematu metodą połowionego kroku, jak opisano w rozdziale 5.4, dla równania $\stackrel{.}{x}=1+x*x$ z $x(0)=0$ z rozwiązaniem $x(t)=\tan(t)$ dla $t=1$ i dla równania wahadła porównując z rozwiązaniem otrzymanym dla równania wahadła przy pomocy funkcji octave'a lsode(). Porównaj czas i ilość iteracji potrzebne do wyliczenia $x_{n}$ każdą z tych metod przy tym samym warunku stopu metody.

Wskazówka:

Metoda Newtona rozwiązywania $G(x)=0$ jest zdefiniowana następująco: $x^{{k+1}}=x^{k}+h^{k}$ dla $DG(x^{k})h^{k}=-G(x^{k})$ . Układ równań liniowych $Ay=b$ możemy rozwiązać w octave przy użyciu operatora backslash tzn. y=A \ b.

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

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

Zagadnienia

6. Sztywność, zmienny krok całkowania i metoda strzałów

6.1. Sztywne równania różniczkowe zwyczajne

Definicja 6.1

6.1.1. Przypadek skalarny

6.1.2. Przypadek wielowymiarowy

6.2. Przykłady schematów sztywnych

6.2.1. Oscylator Van der Pola

6.2.2. Reakcje chemiczne

6.2.3. Równania paraboliczne

6.3. Schematy zamknięte. Predyktor-korektor

6.4. Adaptacyjny dobór kroku całkowania

6.5. Metoda strzałów

6.6. Zadania

Ćwiczenie 6.1

Ćwiczenie 6.2 (laboratoryjne)

Ćwiczenie 6.3 (laboratoryjne)

Ćwiczenie 6.4 (częściowo laboratoryjne)

Ćwiczenie 6.5

Ćwiczenie 6.6 (laboratoryjne)

Ćwiczenie 6.7 (laboratoryjne)