Numerické riešenie diferenciálnych rovníc

Ciele

1. Naučiť sa určovať približné riešenie diferenciálnej rovnice Eulerovou a Heunovou metódou.
2. Naučiť sa určovať približné riešenie diferenciálnej rovnice Rungeho-Kuttovou metódou štvrtého rádu.

Úvod

Riešenie problémov fyzikálneho, technického a popri tom aj matematického charakteru nás často privedie k úlohám, ktoré sa dajú opísať len pomocou obyčajných alebo parciálnych diferenciálnych rovníc. Riešenie daných fyzikálnych a technických problémov je potom v podstate závislé na vyriešení týchto diferenciálnych rovníc. Vyriešiť diferenciálnu rovnicu, respektíve nájsť jej presné riešenie, je však možné len vo veľmi málo prípadoch. Z tohto dôvodu sa na riešenie diferenciálnych rovníc používa celá škála približných metód. V tejto kapitole vynecháme analytické metódy a budeme sa zaoberať len numerickými metódami riešenia diferenciálnych rovníc a ich systémov. Tieto metódy sú založené na hľadaní aproximácie presného riešenia \(y(x) \) v izolovaných bodoch \[x_{0},x_{1},\ldots,x_{n}, \] kde \[x_{i}=x_{0}+ih,\ i=0,1,\ldots,n. \] Číslo \(h \) sa nazýva krok metódy. Pre presné riešenie diferenciálnej rovnice v bode \(x_{i}\) budeme používať označenie \(y\left(x_{i}\right) \) a jeho približnú hodnotu budeme označovať \(y_{i}. \) V rámci štúdia jednotlivých numerických metód sa budeme zaoberať obyčajne jednou diferenciálnou rovnicou \[y'=f\left(x,y\left(x\right)\right) \] so začiatočnou podmienkou \[y\left(x_{0}\right)=y_{0}. \] Úlohu nájsť riešenie diferenciálnej rovnice, ktoré vyhovuje začiatočnej podmienke, budeme ďalej nazývať Cauchyho úlohou pre danú diferenciálnu rovnicu. Presným riešením Cauchyho úlohy pre jednu diferenciálnu rovnicu \[y'=f\left(x,y\left(x\right)\right) \] je funkcia \(y\left(x\right)\) , pre ktorú platí: \[ y(x)=y(x_{0}) +\int^{x}_{x_{0}} f\left(t,y\left(t\right)\right) dt . \] Ak bod \(x\) nahradíme bodom \(x_{i+1}\) a za začiatočný bod zvolíme \(x_{i}\) , tak pomocou tejto integrálnej rovnice vypočítame presnú hodnotu riešenia \( y(x)\) v uzlovom bode \(x_{i+1}\) , teda \[ y(x_{i+1})=y(x_{i}) +\int^{x_{i+1}}_{x_{i}} f\left(t,y\left(t\right)\right) dt . \] Označme \(F\left(x\right)=f\left(x,y\left(x\right)\right)\). Na intervale \(\left\langle x_{i},x_{i+1}\right\rangle\) ľahko vypočítame presnú hodnotu funkcie \(F\) len v uzlovom bode \(x_{i}\), pretože \(F\left(x_{i}\right)=f\left(x_{i},y\left(x_{i}\right)\right)\) a \(y\left(x_{i}\right)=y_{i}\) je známa začiatočná hodnota riešenia \(y\left(x\right)\) v začiatočnom bode \(x_{i}\). V ostatných bodoch intervalu nepoznáme presnú hodnotu funkcie \(F\), lebo nepoznáme presnú hodnotu riešenia \(y\left(x\right)\). Vypočítame približnú hodnotu riešenia \(y\left(x\right)\) v uzlovom bode \(x_{i+1}\). Určitý integrál \(\int^{x_{i+1}}_{x_{i} }f\left(t,y\left(t\right)\right)dt\) nahradíme súčinom \(hF\left(x_{i}\right)\), ktorý je vlastne plošným obsahom obdĺžnika s výškou \(F\left(x_{i}\right)\) a šírkou \(h\). Na základe tejto aproximácie integrálu vypočítame približnú hodnotu \(y_{i+1}\) riešenia \(y\left(x\right)\) v uzlovom bode \(x_{i+1}\) pomocou vzorca \[y_{i+1}=y_{i}+hf\left(x_{i},y_{i}\right)\] alebo \[y_{i+1}=y_{i}+k_{1},\] kde \[k_{1}=hf\left(x_{i},y_{i}\right).\]Tento vzorec je známy ako rekurentný vzorec Eulerovej metódy. Ak určitý integrál \(\int^{x_{i+1}}_{x_{i} }f\left(t,y\left(t\right)\right)dt\) nahradíme určitým integrálom \(\int^{x_{i+1}}_{x_{i} }L_{1}\left(t,x_{i},x_{i+1}\right)dt\) teda na aproximáciu funkcie \(f\left(t,y\left(t\right)\right)\) použijeme Lagrangeov interpolačný polynóm prvého stupňa, tak pre výpočet \(y_{i+1}\) dostaneme vzorec \[y_{i+1}=y_{i}+\frac{h}{2}\left[f\left(x_{i},y_{i}\right)+f\left(x_{i+1},y\left(x_{i+1}\right)\right)\right].\] Hodnotu \(y\left(x_{i+1}\right)\) nahradíme hodnotou \(y_{i+1}\), ktorú vypočítame pomocou Eulerovho rekurentného vzorca \[y_{i+1}=y_{i}+k_{1},\] kde \[k_{1}=hf\left(x_{i},y_{i}\right).\] Vo vzorci funkciu \(f\left(x_{i+1},y\left(x_{i+1}\right)\right)\) zameníme funkciou \(f\left(x_{i}+h,y_{i}+k_{1}\right)\) a keď označíme \(k_{2}=f\left(x_{i}+h,y_{i}+k_{1}\right)\), tak dostaneme rekurentný vzorec Heunovej metódy \[y_{i+1}=y_{i}+\frac{k_{1}+k_{2}}{2}.\] Ak určitý integrál \(\int^{x_{i+1}}_{x_{i} }f\left(t,y\left(t\right)\right)dt\) nahradíme určitým integrálom \(\int^{x_{i+1}}_{x_{i} }L_{2}\left(t,x_{i},x_{i+1}\right)dt\), teda na aproximáciu funkcie \(f\left(t,y\left(t\right)\right)\) použijeme Lagrangeov interpolačný polynóm druhého stupňa, tak pre výpočet \(y_{i+1}\) dostaneme vzorec \[y_{i+1}=y_{i}+\frac{h}{6}\left[f\left(x_{i},y_{i}\right)+4f\left(x_{i}+\frac{h}{2},y \left(x_{i}+\frac{h}{2}\right)\right)+f\left(x_{i}+h,y \left(x_{i}+h\right)\right)\right].\] Presné hodnoty \( y\left(x_{i}+\frac{h}{2}\right) \) a \( y \left(x_{i}+h\right)\) nahradíme približnými. Stredný člen \[ 4f\left(x_{i}+\frac{h}{2},y \left(x_{i}+\frac{h}{2}\right)\right)\] rozpíšeme ako súčet dvoch členov s násobkami \(2.\) Keď označíme \[k_{1}=hf\left(x_{i},y_{i}\right),\ k_{2}=hf\left(x_{i}+\frac{h}{2},y_{i}+\frac{k_{1}}{2}\right), \] \[k_{3}=hf\left(x_{i}+\frac{h}{2},y_{i}+\frac{k_{2}}{2}\right),\ k_{4}=hf\left(x_{i}+h,y_{i}+k_{3}\right), \] tak dostaneme rekurentný vzorec Rungeho-Kuttovej metódy 4. rádu \[y_{i+1}=y_{i}+\frac{1}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right).\]
Poznámka: Medzi Rungeho-Kuttove metódy zaraďujeme aj Eulerovu a Heunovu metódu.

Postup

1. Eulerova a Heunova metóda.
   
   Určíme približné hodnoty riešenia Cauchyho úlohy pomocou Eulerovej metódy.
   Na konečnom intervale \( \left\langle x_{0},x_{n}\right\rangle\) v uzlových bodoch \( x_{i}\) vypočítame pomocou rekurentného vzorca Eulerovej metódy \[y_{i+1}=y_{i}+k_{1},\ i=0,1,\ldots,n-1,\] približné hodnoty \( y_{i}\) riešenia Cauchyho úlohy.
   Vypočítané približné hodnoty riešenia Cauchyho úlohy zapíšeme do tabuľky.
   Nájdeme riešenie Cauchyho úlohy pomocou Heunovej metódy.
   Na konečnom intervale \( \left\langle x_{0},x_{n}\right\rangle\) v uzlových bodoch \( x_{i}\) vypočítame pomocou rekurentného vzorca Heunovej metódy \[y_{i+1}=y_{i}+\frac{k_{1}+k_{2}}{2},\ \ i=0,1,\ldots,n-1,\] približné hodnoty \( y_{i}\) riešenia Cauchyho úlohy.
   Vypočítané približné hodnoty riešenia Cauchyho úlohy zapíšeme do tabuľky.

   Príklad: Cauchyho úlohu \[y'=1+2\left(y-x\right),\ y\left(0\right)=1 \] riešte približne na intervale \[\left\langle 0;0,2\right\rangle \] najprv Eulerovou a potom Heunovou metódou s krokom \(h=0,1. \) Výsledky porovnajte s hodnotami presného riešenia danej úlohy.

   Poznámka: V tomto prípade vieme vypočítať hodnoty presného riešenia Cauchyho úlohy, takže bez väčších problémov sa nám podarí porovnať v tabuľke približné a presné riešenie. V bežnej praxi však vo väčšine prípadov nevieme nájsť presné riešenie Cauchyho úlohy.

   
   

   Riešenie:
   Použijeme rekurentný vzorec Eulerovej metódy \[y_{i+1}=y_{i}+k_{1}, \] kde \[k_{1}=hf\left(x_{i},y_{i}\right) \] a \[f\left(x,y\right)=1+2\left(y-x\right). \] Teda \[y_{1}=y_{0}+hf\left(x_{0},y_{0}\right)=1+0,1\left(1+2\right)=1+0,3=1,3, \] \[y_{2}=y_{1}+hf\left(x_{1},y_{1}\right)=1,3+0,1\left(1+2\left(1,3-0,1\right) \right)=1,3+0,34=1,64. \] Porovnanie približne vypočítaných hodnôt riešenia s hodnotami presného riešenia je uvedené v tabuľke \[\begin{array}{|c|c|c|c|c|} \def\a{&} \hline i\a x_i \a y_i \a y(x_i)\a e_i=\left|y(x_i)-y_i\right| \\ \hline 0 \a 0 \a 1 \a 1 \a 0 \\ \hline 1 \a 0{,}1 \a 1{,}3 \a 1{,}3214 \a 0{,}0214 \\ \hline 2 \a 0{,}2 \a 1{,}64 \a 1{,}6918 \a 0{,}0518 \\ \hline \end{array} \]
   Ukážeme, že Heunova metóda je oveľa presnejšia ako predchádzajúca Eulerova metóda. Jedinou nevýhodou použitia Heunovej metódy v porovnaní s Eulerovou metódou je zložitejší rekurentný vzorec, ktorý obsahuje o jednu smernicu viac ako rekurentný vzorec Eulerovej metódy. Vypočítame približné hodnoty \(y_{i}\) riešenia Cauchyho úlohy pomocou rekurentného vzorca Heunovej metódy \[y_{i+1}=y_{i}+\frac{k_{1}+k_{2}}{2},\] kde \[k_{1}=hf\left(x_{i},y_{i}\right),\ k_{2}=hf\left(x_{i}+h,y_{i}+k_{1}\right).\] Porovnanie približne vypočítaných hodnôt \(y_{i}\) riešenia s hodnotami \(y\left(x_{i}\right)\) presného riešenia Cauchyho úlohy je uvedené v tabuľke \[\begin{array}{|c|c|c|c|c|} \def\a{&} \hline i\a x_i \a y_i \a y(x_i)\a e_i=\left|y(x_i)-y_i\right| \\ \hline 0 \a 0 \a 1 \a 1 \a 0 \\ \hline 1 \a 0{,}1 \a 1{,}32 \a 1{,}3214 \a 0{,}0014 \\ \hline 2 \a 0{,}2 \a 1{,}6884 \a 1{,}6918 \a 0{,}0034 \\ \hline \end{array} \]

   Zobraziť riešenieSkryť riešenie
2. Rungeho-Kuttova metóda 4. rádu.
   
   V prvom kroku pomocou rekurentného vzorca Rungeho-Kuttovej metódy \[y_{i+1}=y_{i}+\frac{1}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\] na konečnom intervale \(\left\langle x_{0},x_{n}\right\rangle \) v uzlových bodoch \( x_{i}\) vypočítame približné hodnoty \( y_{i}\) riešenia Cauchyho úlohy.
   V druhom kroku z vypočítaných približných hodnôt riešenia Cauchyho úlohy Rungeho-Kuttovou metódou 4. rádu vytvoríme prehľadnú tabuľku.

   Príklad: Cauchyho úlohu \[y'=\left(xy-5\right)e^{2x+3}, \qquad y\left(0\right)=1\] riešte približne na intervale\( \left\langle 0;0,2\right\rangle \) s krokom \( h=0,1 \) Rungeho-Kuttovou metódou 4. rádu.
   

   Riešenie:
   Použijeme rekurentný vzorec Rungeho-Kuttovej metódy \[y_{i+1}=y_{i}+\frac{1}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\] kde \[k_{1}=hf\left(x_{i},y_{i}\right),\ k_{2}=hf\left(x_{i}+\frac{h}{2},y_{i}+\frac{k_{1}}{2}\right), \] \[k_{3}=hf\left(x_{i}+\frac{h}{2},y_{i}+\frac{k_{2}}{2}\right),\ k_{4}=hf\left(x_{i}+h,y_{i}+k_{3}\right). \] Vypočítané jednotlivé smernice \[ k_{i},\ i=1,2,3,4 \] a približné hodnoty \(y_{i}\) riešenia Cauchyho úlohy sú uvedené v tabuľke: \[\begin{array}{|c|c|c|c|c|c|c|} \def\a{&} \hline i\a x_i \a y_i \a k_1\a k_2\a k_3\a k_4 \\ \hline 0 \a 0 \a 1 \a -10{,}042768 \a -11{,}545308\a -11{,}628692\a -14{,}873752 \\ \hline 1 \a 0{,}1 \a -10{,}3877420 \a -14{,}934772 \a -21{,}016962 \a -22{,}253744\a -34{,}836960\\ \hline 2 \a 0{,}2 \a -33{,}596279 \\ \hline \end{array} \]

   Zobraziť riešenieSkryť riešenie

Zdroje

1. Buša, Pirč, Schrötter: Numerické metódy, pravdepodobnosť a matematická štatistika, 2006, ISBN 80-8073-632-4. Stiahnuť obrazovkovú alebo tlačovú verziu.
2. Daňo, Ostertagová: Vybrané kapitoly z numerických metód, pravdepodobnosti a matematickej štatistiky, Equilibria s.r.o., Košice, 2012, ISBN 978-80-8143-012-1.

Doplňujúce úlohy

Úloha:
Cauchyho úlohu \[ y'=xy,\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,5\right\rangle \) Eulerovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3\a 0,4\a 0,5 \\ \hline y_i \a 1 \a 1 \a 1,01\a 1,0302\a 1,061106\a 1,103550 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=\frac{x}{y},\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,5\right\rangle \) Eulerovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3\a 0,4\a 0,5 \\ \hline y_i \a 1 \a 1 \a 1,01\a 1,029802\a 1,058934\a 1,096708 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=\frac{x}{y^{2}}-1,\ y(0)=1 \] riešte na intervale \(\left\langle 1;1,4\right\rangle \) Eulerovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|c|} \hline x_i \a 1 \a 1,1 \a 1,2 \a 1,3\a 1,4 \\ \hline y_i \a 1 \a 1 \a 1,01\a 1,027636\a 1,050738 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=\frac{\cos y}{1,5+x}+0,1y^{2},\ y(0)=0 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) Eulerovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 0 \a 0,066667 \a 0,129072\a 0,187573 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=xy,\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) Heunovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 1 \a 1,005 \a 1,020175\a 1,045986 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=y^{2}e^{x}-2y,\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) Heunovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 1 \a 0,904759 \a 0,818582\a 0,740606 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=x^{2}(1+y),\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) Heunovou metódou s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 1 \a 1,001 \a 1,006007\a 1,019082 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=xy,\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 1 \a 1,005013 \a 1,020201\a 1,046028 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=x^{2}(1+y),\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 1 \a 1,000667 \a 1,005340\a 1,018081 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=y^{2}-3xy,\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,3\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \a 0,3 \\ \hline y_i \a 1 \a 1,093964 \a 1,171454\a 1,225463 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=2xy-1,\ y(1)=0 \] riešte na intervale \(\left\langle 1;1,3\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|c|} \hline x_i \a 1 \a 1,1 \a 1,2 \a 1,3 \\ \hline y_i \a 0 \a -0,111460 \a -0,252903\a -0,438528 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=-y+\ln x,\ y(1)=1 \] riešte na intervale \(\left\langle 1;1,4\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,2\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|} \hline x_i \a 1 \a 1,2 \a 1,4 \\ \hline y_i \a 1 \a 0,836311 \a 0,732562 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=-(y+1)\cos x,\ y(0)=1 \] riešte na intervale \(\left\langle 0;0,2\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,1\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|} \hline x_i \a 0 \a 0,1 \a 0,2 \\ \hline y_i \a 1 \a 0,809976 \a 0,639642 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=2x^{2}-y,\ y(0,5)=1 \] riešte na intervale \(\left\langle 0,5;0,6\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,05\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|} \hline x_i \a 0,5 \a 0,55 \a 0,6 \\ \hline y_i \a 1 \a 0,978156 \a 0,962744 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok
Úloha:
Cauchyho úlohu \[ y'=x^{2}-y,\ y(0,5)=1 \] riešte na intervale \(\left\langle 0,5;0,6\right\rangle \) štandardnou metódou Rungeho-Kutta 4.rádu s krokom \(h=0,05\).

Výsledok:

Výsledok: \[ \begin{array}{|c|c|c|c|} \hline x_i \a 0,5 \a 0,55 \a 0,6 \\ \hline y_i \a 1 \a 0,964693 \a 0,933791 \\ \hline \end{array} \]
Zobraziť výsledokSkryť výsledok