Zagadnienia

3. Metody dla równań różniczkowych zwyczajnych

W tym i kilku następnych rozdziałach zajmiemy się schematami rozwiązywania równań różniczkowych zwyczajnych. Ten rozdział jest poświęcony wprowadzeniu najprostszych metod (schematów) rozwiązywania równań różniczkowych zwyczajnych.

3.1. Wprowadzenie

Załóżmy, że rozpatrujemy zagadnienie początkowe pierwszego rzędu (zagadnienie Cauchy'ego) :

\left\{\begin{array}[]{lcl}\frac{dx}{dt}&=&f(t,x)\\
x(t_{0})&=&x_{0}\end{array}\right. (3.1)

gdzie G\subset\mathbb{R}^{m}, f:(a,b)\times G\rightarrow\mathbb{R}^{m} jest funkcją ciągłą, a t_{0}\in(a,b),x_{0}\in G jest ustalone.

Z ogólnej teorii równań różniczkowych, por. [23] wiadomo, że

Twierdzenie 3.1 (Peano)

Jeśli f jest funkcją ciągłą na otoczeniu (t_{0},x_{0}), to istnieje rozwiązanie (3.1) określone na pewnym otoczeniu t_{0}.

Jeśli dodatkowo założymy, że f jest funkcją lipschitzowska na otoczeniu (t_{0},x_{0}) względem zmiennej x, to możemy pokazać jednoznaczności rozwiązania, tzn. zachodzi twierdzenie:

Twierdzenie 3.2 (Picarda-Lindelöfa)

Jeśli f jest funkcją ciągłą na otoczeniu (t_{0},x_{0}) oraz f jest funkcją lipschizowską względem x w pewnej kuli B((t_{0},x_{0}),\delta), tzn.

\exists L\geq 0\quad\forall(t,x),(t,y)\in B((t_{0},x_{0}),\delta)\qquad\| f(t,x)-f(t,y)\|\leq L\| x-y\|,

to istnieje c>0 i x\in C^{1}((t_{0}-c,t_{0}+c),R^{n}) takie, że x jest jednoznacznym rozwiązaniem (3.1).

Od tej pory będziemy przyjmować, że funkcja f, zwana też polem wektorowym, spełnia założenia twierdzenia Picarda-Lindelöfa, tzn. Twierdzenia 3.2, czyli że istnieje jednoznaczne rozwiązanie zadania Cauchy'ego na odcinku [t_{0},T].

Zauważmy, że każde rozwiązanie jest krzywą styczną do pola wektorowego.

3.2. Równania liniowe ze stałymi współczynnikami

W tym podrozdziale krótko przypomnimy teorie dla ważnej klasy równań różniczkowych zwyczajnych, tzn. jednorodnych równań liniowych ze stałymi współczynnikami, czyli równań postaci:

\frac{dx}{dt}=Ax

gdzie A - to stała macierz n×n, dla których znamy rozwiązania zadania Cauchy'ego:

\frac{dx}{dt}=Ax\qquad x(t_{0})=x_{0} (3.2)

Znamy wzór na rozwiązanie tego zadania:

x(t)=e^{{A\,(t-t_{0})}}\, x_{0},

gdzie eksponent od macierzy zdefiniowany jest wzorem

\exp(B)=e^{B}=\sum _{{k=0}}^{\infty}\frac{B^{k}}{k!}.

Skorzystamy ze znajomości postaci rozwiązania tej klasy równań w rozdziale 6.

W zależności od postaci Jordana macierzy A można wypisać postać \exp(A\, t), w szczególności jeśli macierz A jest diagonalizowalna, tzn. istnieje baza wektorów własnych, które zapisane jako kolumny macierzy V dają:

V\Lambda V^{{-1}}=A

gdzie \Lambda - to macierz diagonalna z wartościami własnymi macierzy A na diagonali:

\Lambda=\begin{pmatrix}\lambda _{1}&&\\
&\ddots&\\
&&\lambda _{n}\end{pmatrix}.

Wtedy wiadomo, że

e^{{A\, t}}=V\,\begin{pmatrix}e^{{\lambda _{1}t}}&&\\
&\ddots&\\
&&e^{{\lambda _{n}t}}\end{pmatrix}V^{{-1}}.

W szczególności jeśli Re\,\lambda _{k}<0 dla wszystkich k=1,\ldots,n to każde rozwiązanie spełnia

\lim _{{t\rightarrow+\infty}}\| x(t)\|=0,

a z kolei jeśli istnieje \lambda _{k} takie, że Re\lambda _{k}>0, to istnieje rozwiązanie zagadnienia Cauchy'go z niezerowym warunkiem brzegowym, dla którego

\lim _{{t\rightarrow+\infty}}\| x(t)\|=+\infty.

3.3. Kilka prostych schematów

Załóżmy, że rozpatrujemy zadanie skalarne tzn. m=1. Chcemy w jakiś sposób przybliżyć rozwiązanie równania (3.1). Przybliżamy pochodną poprzez iloraz różnicowy dla pewnego parametru h>0:

\frac{x(t+h)-x(t)}{h}\approx\frac{dx}{dt}

i otrzymujemy otwarty schemat Eulera:

\frac{x_{h}(t+h)-x_{h}(t)}{h}=f(t,x_{h}(t))

czy inaczej:

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

znając rozwiązanie w punkcie t_{0} x_{h}(t_{0})=x_{0} możemy wyznaczyć przybliżone rozwiązanie x_{h} w kolejnych punktach t_{n}=t_{0}+n\, h z powyższego wzoru. Ale można też przybliżyć pochodną biorąc parametr -h w tył:

\frac{x_{h}(t)-x_{h}(t-h)}{h}\approx\frac{dx}{dt}

i wtedy zastępując pochodną przez taki iloraz otrzymujemy zamknięty schemat Eulera:

\frac{x_{h}(t)-x_{h}(t-h)}{h}=f(t,x_{h}(t))

czy inaczej:

x_{h}(t+h)=x_{h}(t)+h\, f(t,x_{h}(t+h)).

Proszę zauważyć, że jeśli znamy rozwiązanie w punkcie t_{0}, tzn. x_{h}(t_{0})=x_{0}, to aby wyznaczyć kolejne przybliżenia rozwiązania w punktach t=t_{n}=t_{0}+nh należy rozwiązać równania postaci:

g(y):=y-h\, f(t,y)-x_{h}(t)=0 (3.3)

względem y, co sprawia, że zamknięty schemat Eulera może wydać się mało praktyczny w porównaniu z otwartym schematem Eulera. Dla niektórych równań jest to pozorne. Zauważmy tylko, że im h mniejsze, tym potencjalnie równanie (3.3) jest łatwiejsze do rozwiązania (dlaczego?). Przyjrzymy się temu problemowi dokładniej w kolejnych rozdziałach.

W dalszej części wykładu założymy, że chcemy przybliżyć rozwiązanie x(t) na odcinku [t_{0},T], na którym x\in C^{k}([t_{0},T]), w dyskretnych punktach czasu:

t_{n}\equiv t_{n}^{h}=t_{0}+n\, h\qquad h>0.

Często będziemy opuszczali indeks h, o ile to nie będzie powodowało niejasności. Wartość rozwiązania w punkcie t_{n}^{h}, czyli x(t_{n}) będzie przybliżana przez x_{n}^{h}, spełniające odpowiedni schemat. Wygodnie jest też oznaczać f(t_{n}^{h},x_{n}^{h})=f_{n}^{h}. Górny indeks h będziemy często opuszczali, jeśli h będzie ustalone.

Tak więc otwarty schemat Eulera możemy zapisać jako:

x_{{n+1}}=x_{n}+h\, f_{n}\qquad n>0\qquad x_{0}=x(t_{0}), (3.4)

a zamknięty schemat Eulera możemy zapisać jako:

x_{{n+1}}=x_{n}+h\, f_{{n+1}}\qquad n>0\qquad x_{0}=x(t_{0}). (3.5)

Kolejnym wyprowadzeniem otwartego schematu Eulera (3.4) jest obcięcie rozwinięcia szeregu Taylora rozwiązania:

x(t+h)=x(t)+\frac{dx}{dt}(t)\, h+\frac{1}{2}\frac{d^{2}x}{dt^{2}}(t)\, h^{2}+\ldots (3.6)

Zostawiamy tylko pierwsze dwa człony i otrzymujemy

x(t+h)\approx x(t)+\frac{dx}{dt}(t)\, h=x(t)+f(x(t),t)\, h

czyli wstawiając x_{n} za przybliżenie x(t_{n}), a x_{{n+1}} za przybliżenie x(t_{n}+h)=x(t_{{n+1}}) otrzymujemy znów otwarty schemat Eulera (3.4). Analogicznie możemy wyprowadzić zamknięty schemat Eulera (3.5) rozwijając rozwiązanie w t dla h<0.

Jeszcze innym intuicyjnym wyprowadzeniem schematu otwartego Eulera jest podążanie za polem wektorowym. Jak wiemy, wykresem rozwiązania równania różniczkowego jest krzywa styczna do zadanego pola wektorowego f(t,x) spełniająca odpowiedni warunek początkowy. Zatem znając rozwiązanie przybliżone dla t_{n} tzn. mając x_{n}, możemy wyznaczyć x_{{n+1}}, przybliżenie rozwiązania x(t_{{n+1}}), biorąc poprawkę w kierunku pola wektorowego tzn.:

x_{{n+1}}=x_{{n}}+h\, f(x_{n},t_{n}).

Czyli znów otrzymujemy otwarty schemat Eulera.

Zadajmy pytanie, czy takie schematy są wystarczająco dokładne. Czy one działają stabilnie na dłuższych odcinkach czasu, na których istnieje rozwiązanie?

Sprawdźmy, co się dzieje dla modelowego zadania:

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

którego rozwiązaniem jest x(t)=e^{{at}}.

Otwarty schemat Eulera daje nam ciąg:

x_{n}^{h}=x_{{n-1}}+h\, a\, x_{{n-1}}=(1+h\, a)x_{{n-1}}=(1+h\, a)^{n}.

Ustalmy t=h\, n, czyli h=t/n. Wtedy

x_{n}^{h}=(1+h\, a)^{n}=(1+t\, a/n)^{n}\rightarrow e^{{at}}\qquad n\rightarrow\infty.

Dla zamkniętego schematu Eulera otrzymujemy analogicznie:

x_{n}^{h}=x_{{n-1}}^{h}+h\, a\, x_{n}^{h},

czyli

x_{n}=(1-h\, a)^{{-1}}\, x_{{n-1}}=(1-a\, h)^{{-n}}=\frac{1}{(1-a\, t/n)^{n}}\rightarrow\frac{1}{e^{{-a\, t}}}=e^{{a\, t}}.

Popatrzmy jak te dwa schematy działają (w praktyce) na wykresach dla a=1 i x(0)=x_{0}=1, por. rysunki 3.1 - 3.4 dla otwartego schematu Eulera i rysunki 3.5 - 3.8 dla zamkniętego schematu Eulera.

\
Rys. 3.1. Otwarty Euler - cztery punkty.
\
Rys. 3.2. Otwarty Euler - 10 punktów.
\
Rys. 3.3. Otwarty Euler - 100 punktów.
\
Rys. 3.4. Otwarty Euler - na [0,100], 100 punktów.
\
Rys. 3.5. Zamknięty Euler - cztery punkty.
\
Rys. 3.6. Zamknięty Euler - 10 punktów.
\
Rys. 3.7. Zamknięty Euler - 100 punktów.
\
Rys. 3.8. Zamknięty Euler - na [0,100], 100 punktów.

Zauważmy, że wykres rozwiązania ze schematu Eulera otwartego jest poniżej wykresu dokładnego rozwiązania, a dla schematu zamkniętego - powyżej, co widać lepiej na rysunku 3.9.

\
Rys. 3.9. Schematy Euler - na [0,1], 20 punktów.

Popatrzmy na przypadek dwuwymiarowy. Weźmy modelowe zadanie wahadła. Dla małych prędkości możemy przyjąć, że \sin(x)\approx x, stąd otrzymujemy równanie liniowe ze stałymi współczynnikami (zlinearyzowane równanie wahadła):

\frac{d^{2}x}{dt^{2}}=-a\, x,

gdzie x to prędkość kątowa, a a=g/l>0 dla g wartości przyspieszenia ziemskiego i l długości wahadła.

Zapisując to równanie jako układ dwóch równań pierwszego rzędu otrzymujemy:

\begin{array}[]{l}\frac{dx}{dt}=y\\
\frac{dy}{dt}=-a\, x.\\
\end{array}

Przyjmijmy, że a=1.

Znamy rozwiązania:

x(t)=c_{1}\,\sin(t)+c_{2}\,\cos(t),

czyli trajektorie rozwiązania zawarte są w okręgach.

A teraz zastosujmy otwarty schemat Eulera do tego równania z warunkiem początkowym (x(0),y(0))^{T}=(0,1)^{T}, którego rozwiązaniem jest x(t)=\sin(t) z y(t)=\cos(t):

\left\{\begin{array}[]{l}x_{{n+1}}=x_{n}+h\, y_{n}\\
y_{{n+1}}=y_{n}-h\, x_{n}\end{array}\right.\quad n=0,1,\ldots,N

dla ustalonego h>0 i T=N\, h.

Zatem: x_{n}\approx x(t_{n})=\sin(n\, h), a y_{n}\approx y(t_{n})=\cos(n\, h) z x_{0}=0 i y_{0}=1

Dla zamkniętego schematu Eulera jest analogicznie:

\left\{\begin{array}[]{l}x_{{n+1}}=x_{n}+h\, y_{{n+1}}\\
y_{{n+1}}=y_{n}-h\, x_{{n+1}}\end{array}\right.\quad n=0,1,\ldots,N,

czy równoważnie

\left\{\begin{array}[]{l}x_{{n+1}}-h\, y_{{n+1}}=x_{n}\\
y_{{n+1}}+h\, x_{{n+1}}=y_{n}\end{array}\right.\quad n=0,1,\ldots,N

czyli w każdym kroku dla ustalonego n musimy rozwiązać układ dwóch równań liniowych.

Popatrzmy teraz na rysunki 3.10 - 3.12.

\
Rys. 3.10. Schematy Eulera dla równania 2-wymiarowego - rozwiązanie.
\
Rys. 3.11. Schematy Eulera dla równania 2-wymiarowego - pochodna rozwiązania.
\
Rys. 3.12. Schematy Eulera dla równania 2-wymiarowego - trajektoria.

Widać, że mimo małego kroku rzędu (h=1e-2) wyniki są wyraźnie gorsze niż w przypadku skalarnym, mimo że wyjściowe równanie różniczkowe jest liniowe.

Rozważmy wyjściowe równanie wahadła, por. Przykład 2.5. Znów przyjmijmy, że g=l i warunek początkowy x(0)=0 i y(0)=1. Wtedy schematy Eulera przybierają odpowiednio formę:

schemat otwarty Eulera:

\left\{\begin{array}[]{l}x_{{n+1}}=x_{n}+h\, y_{n}\\
y_{{n+1}}=y_{n}+h\,\sin(x_{n})\end{array}\right.\quad n=0,1,\ldots,N

z x_{0}=x(0)=0 i y_{0}=y(0)=1. Znając x_{n},y_{n} otrzymujemy natychmiast wzór na x_{{n+1}},y_{{n+1}}.

W przypadku schematu zamkniętego Eulera:

\left\{\begin{array}[]{l}x_{{n+1}}=x_{n}+h\, y_{{n+1}}\\
y_{{n+1}}=y_{n}-h\,\sin(x_{{n+1}})\end{array}\right.\quad n=0,1,\ldots,N

z x_{0}=x(0)=0 i y_{0}=y(0)=1, musimy w każdym kroku rozwiązać układ równań nieliniowych:

\left\{\begin{array}[]{l}x_{{n+1}}-h\, y_{{n+1}}=x_{n}\\
y_{{n+1}}+h\,\sin(x_{{n+1}})=y_{n}\end{array}\right.\quad n=0,1,\ldots,N

Im h jest bliższe zera, tym układ jest łatwiejszy do rozwiązania.

Można pokazać, że rozwiązanie wyjściowego równania ma trajektorie okresowe, co potwierdza wykres na rysunku 3.15 (tu wyliczony przy pomocy dużo dokładniejszego schematu niż schematy Eulera). W kolejnych rysunkach 3.10- 3.13 - prezentujemy przybliżone rozwiązania dla nieliniowego równania wahadła, otrzymane przy pomocy obu schematów Eulera.

\
Rys. 3.13. Schematy Eulera dla równania 2-wymiarowego - rozwiązanie.
\
Rys. 3.14. Schematy Eulera dla równania 2-wymiarowego - pochodna rozwiązania.
\
Rys. 3.15. Schematy Eulera dla równania 2-wymiarowego - trajektoria.

3.3.1. Absolutna stabilność schematów Eulera

Rozpatrzmy ponownie modelowe zadanie skalarne, ale na długich odcinkach czasu:

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

Rozwiązaniem jest x(t)=\exp(a\, t) i wtedy \lim _{{t\rightarrow+\infty}}x(t)=0. Im |a| większe, tym rozwiązanie szybciej zbiega do zera.

Rozpatrzmy teraz zastosowanie otwartego i zamkniętego schematu Eulera do rozwiązania tego zagadnienia. Dla otwartego schematu Eulera wiemy już, że:

x_{n}=(1+a\, h)^{n}.

Zauważmy, że przy ustalonym h ciąg przybliżeń x_{n} jest dodatni i zbiega do zera dla n\rightarrow+\infty o ile zachodzi warunek:

h<-1/a.

W przypadku gdy parametr a jest ujemny i o dużym module, warunek ten wymusza to, że musimy wziąć bardzo małe h, aby otrzymać schematem otwartym Eulera rozwiązanie przybliżone, które jest dodatnie i malejące do zera, czyli zachowujące się jak rozwiązanie zagadnienia początkowego: \exp(a\, t).

Natomiast dla zamkniętego schematu Eulera widzimy, że:

x_{n}=(1-a\, h)^{{-n}}.

Otrzymujemy wtedy, że dla dowolnego a<0 zachodzi x_{n}>0 i x_{n}\rightarrow 0 dla n\rightarrow+\infty, czyli nie otrzymujemy żadnego ograniczenia na krok h, co jest istotne, jeśli chcemy rozwiązywać równanie na długim odcinku czasu.

Schemat zamknięty Eulera można uznać za lepszy od schematu otwartego dla tego zagadnienia dla ujemnego a o bardzo dużym module, szczególnie na długim odcinku czasu, ponieważ nie wymusza żadnych ograniczeń na krok h. Wrócimy do tego problemu w następnych rozdziałach.

3.4. Zadania

Ćwiczenie 3.1

Czy rozwiązanie y(x) zagadnienia początkowego:

\frac{dy}{dx}=x^{{2/3}}\qquad y(0)=0.

jest wyznaczone jednoznacznie? Znajdź wszystkie rozwiązania y(x) tego zagadnienia początkowego.

Wskazówka: 

Jest to równanie o zmiennych rozdzielonych (autonomiczne) \frac{dy}{dx}=f(y)g(x) z warunkiem początkowym y(x_{0})=y_{0}, więc w postaci uwikłanej rozwiązanie ma postać \int _{{y_{0}}}^{{y}}1/f(y)dy=\int _{{x_{0}}}^{{x}}g(x)dx.

Ćwiczenie 3.2 (laboratoryjne)

Zaimplementuj w octavie otwarty schemat Eulera i zastosuj go do równania:

\frac{dy}{dx}=x^{{2/3}}\qquad y(0)=0.

na różnych odcinkach czasu np. [0,1] lub [0,10] i różnych wartości h np. h=1e-1,1e-2,1e-4. Zmniejszając h sprawdź, czy ten schemat znajdzie rozwiązanie różne od zera. Następnie weź przybliżenie startowe na poziomie błędu zaokrągleń np. x_{0}=10^{{-16}} i sprawdź, jakie schemat znajduje rozwiązania; w szczególności, czy są one różne od zera.

Ćwiczenie 3.3

Rozpatrzmy równanie różniczkowe zwyczajne liniowe jednorodne rzędu n o stałych współczynnikach:

\frac{d^{n}x}{dt^{n}}(t)+a_{{n-1}}\frac{d^{{n-1}}x}{dt^{{n-1}}}(t)+\ldots+a_{0}x(t)=0.
  1. Poprzez podstawienie x_{k}(t)=\frac{d^{k}x}{dt^{k}} k=0,\ldots,n-1 sprowadź to równanie do równania liniowego jednorodnego ze stałą macierzą:

    \frac{dx}{dt}=A\,\vec{x},
  2. Znajdź wielomian charakterystyczny A oraz dla n=2 postać Jordana tej macierzy w zależności od tego, jakie wartości własne ma A,

  3. Przy założeniu, że A ma n jednokrotnych wartości własnych rzeczywistych, znajdź e^{{A\, t}} i dla n=2 znajdź rozwiązanie zadania początkowego dla tego równania z warunkami początkowymi: \frac{d^{k}x}{dt^{k}}(0)=b_{k}\qquad k=0,\ldots,n-1.

Ćwiczenie 3.4 (częściowo laboratoryjne)

Dla n=2 i macierzy A kolejno [a,1;0,a], [a,0;0,b], [a,-b;b,a] dla rożnych wartości parametrów a,b, np. a=1,b=10, naszkicuj na kartce portrety fazowe (wykresy trajektorii) równania jednorodnego:

\frac{dx}{dt}=A\,\vec{x}

w otoczeniu zera. Naszkicuj pole wektorowe na ekranie korzystając z funkcji octave'a quiver() i portrety fazowe z pomocą funkcji lsode().

Ćwiczenie 3.5 (laboratoryjne)

Zaimplementuj w octavie otwarty schemat Eulera i zastosuj go do równania \frac{dy}{dx}=a\, y z y(0)=1 dla różnych wartości parametru a np. a=-1e-3,-100,-1,1,10. Narysuj na monitorze wykresy przybliżonych rozwiązań razem z wykresem rozwiązania dokładnego y(t)=\exp(a\, t).

Ćwiczenie 3.6 (częściowo laboratoryjne)

Rozpatrzmy równanie \frac{dy}{dx}=\begin{pmatrix}-20&1\\
-21&1\end{pmatrix}y. Policz wartości własne macierzy A=\begin{pmatrix}-20&1\\
-21&1\end{pmatrix} i porównaj z wynikiem obliczonym w octavie z użyciem odpowiedniej funkcji np. eig(). Znajdź rozwiązanie ogólne tego równania. Przy pomocy otwartego schematu Eulera i funkcji octave'a lsode() rozwiąż to równanie z y(0)=(1,1)^{T} na odcinku [0,100] z h=0.1. Porównaj wyniki rysując wykresy na ekranie obu rozwiązań i rozwiązania dokładnego, które należy też wyznaczyć.

Wskazówka: 

Rozwiązanie ogólne - to \exp(A\, t)\, c, gdzie c wektor stałych a funkcja expm() octave'a pozwala obliczyć eksponent macierzy.

Ćwiczenie 3.7 (częściowo laboratoryjne)

Udowodnij, że przybliżenia rozwiązania układu równań \frac{dx}{dt}=y;\frac{dy}{dt}=-x z x(0)=1,y(0)=0, otrzymane za pomocą otwartego (lub zamkniętego) schematu Eulera, mają normę drugą zbieżną do jeden, tzn. \sqrt{(x_{n})^{2}+(y_{n})^{2}} zbiegają do jeden, dla ustalonego czasu t=n\, h z h dążącym do zera. Zaimplementuj oba schematy Eulera dla tego równania w octave (w przypadku zamkniętego schematu Eulera użyj operatora backslash: \ w każdym kroku czasowym do rozwiązania odpowiedniego układu dwóch równań liniowych). Naszkicuj na ekranie monitora portret fazowy przy pomocy plot(), lsode() i obu schematów dla różnych wartości h. Policz wartości normy drugiej rozwiązań otrzymanych przy pomocy tych schematów i lsode() dla ustalonego czasu np. t=1 czy t=1000 i różnych wartości h.

Ćwiczenie 3.8 (laboratoryjne)

Naszkicuj na ekranie monitora portrety fazowe równań liniowych \frac{dy}{dx}=Ay=[a,b;c,d]y dla macierzy A o różnych postaciach Jordana przy pomocy plot(), lsode().

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.