Obyčajné differenciálne rovnice/Runge-Kutta páry

4.2 Vložené Runge-Kutta páry

Riadenie veľkosti kroku je možné aj pri explicitných Runge-Kutta metódach. Na rozdiel od prístupu opísaného v časti 3.2 pre riadenie veľkosti kroku explicitnej Eulerovej metódy sa tu veľkosť kroku pre odhad lokálnej chyby neznižuje na polovicu; namiesto toho sa používajú dve eRKV s rôznymi poradiami konzistencie. S cieľom minimalizovať počet vyhodnotení funkcie pri výpočte medzikrokov $k_{j}$ tieto dve „vložené“ eRKV používajú spoločné body vzorkovania. Na tomto mieste je potrebné uviesť, že podľa zistení o počte krokov a poradí eRKV si metóda vyššieho rádu vyžaduje ďalšie dodatočné medziľahlé body, aby sa dosiahol tento vyšší rád.

Ako vloženú Runge-Kutta dvojicu chápeme dve eRKV so $s$ spoločnými interpolačnými bodmi (a krokmi). Prvá z eRKV má rád konzistencie $p$ a druhá metóda dosahuje rád konzistencie $s$ s $s+1,s+2,\ldots ,m$ opornými bodmi a ďalšími $s+1,s+2,\ldots ,m$ opornými bodmi. Váhové vektory týchto dvoch eRKV sa vo všeobecnosti líšia a iteračné špecifikácie (aktualizácie Runge-Kutta) týchto dvoch metód sú ${\begin{aligned}&&{\mbox{eRKV1:}}\ \eta _{i+1}&=\eta _{i}+h\sum _{j=1}^{s}b_{j}k_{j}\\&&{\mbox{eRKV2:}}\ {\tilde {\eta }}_{i+1}&={\tilde {\eta }}_{i}+h\sum _{j=1}^{s}{\tilde {b}}_{j}k_{j}+h\sum _{j=s+1}^{m}{\tilde {b}}_{j}k_{j}.\end{aligned}}$ Kvôli spoločným krokom s sa v oboch metódach používajú vyhodnotenia funkcií $k_{1},\ldots k_{s}$. Príkladom vloženej dvojice Runge-Kutta je
Fehlenbergova RK dvojica rádov 2 a 3:

${\begin{aligned}k_{1}=f(t_{i},\eta _{i})\qquad \qquad \qquad \ \ \ \qquad \qquad \\k_{2}=f{\big (}t_{i}+h,\eta _{i}+hk_{1}{\big )}\qquad \qquad \qquad \\k_{3}=f{\big (}t_{i}+{\frac {h}{2}},\eta _{i}+{\frac {h}{4}}k_{1}+{\frac {h}{4}}k_{2}{\big )}\qquad \\{\eta }_{i+1}=\eta _{i}+h{\Big (}{\frac {1}{2}}k_{1}+{\frac {1}{2}}k_{2}{\Big )}\qquad \qquad \\{\tilde {\eta }}_{i+1}={\tilde {\eta }}_{i}+h{\Big (}{\frac {1}{6}}k_{1}+{\frac {1}{6}}k_{2}+{\frac {2}{3}}k_{3}{\Big )}.\end{aligned}}$

4.2.1 Regulácia kroku pre Runge-Kutta páry

V nasledujúcom texte sa obmedzíme na prípad $q=p+1$ a odvodíme vzorec pre riadenie šírky kroku. Nech ${\eta }_{i+1},{\tilde {\eta }}_{i+1}$ sú numerické riešenia v $t_{i+1}$, ktoré sú odvodené z presného riešenia $y(t_{i})$ namiesto $\eta _{i},{\tilde {\eta }}_{i}$ a metódy majú poradie konzistencie $p$ a $p+1$. Potom pre lokálne chyby po použití Taylorovho rozšírenia (s Taylorovým zvyškom) platí ${\begin{aligned}&&\varepsilon (t_{i+1},h)=h\tau (t_{i+1},h)=y(t_{i+1})-\eta _{i+1}=a(t_{i})h^{p+1}+{\cal {O}}(h^{p+2})=a(\theta _{1})h^{p+1},\\&&{\tilde {\varepsilon }}(t_{i+1},h)=h{\tilde {\tau }}(t_{i+1},h)=y(t_{i+1})-{\tilde {\eta }}_{i+1}=b(t_{i})h^{p+2}+{\cal {O}}(h^{p+3})=b(\theta _{2})h^{p+2},\end{aligned}}$

wobei $\theta _{1},\theta _{2}\in [t_{i},t_{i}+h]$. $a(t_{i})h^{p+1},b(t_{i})h^{p+1}$ nazývame vedúce členy chyby oboch metód. Sú to členy s najmenšou mocninou $h$ a obsahujú až $(p+1)$-tú a $(p+2)$-tú parciálnu deriváciu funkcie $f$ v $t_{i}$. Po odčítaní týchto dvoch rovníc dostaneme ${\tilde {\eta }}_{i+1}-{\eta }_{i+1}=a(t_{i})h^{p+1}+{\tilde {\cal {O}}}(h^{p+2})\approx h\tau (t_{i+1},h),$ $\qquad \qquad (4.9)$

t. j. lokálnu chybu prvej metódy (eRKV1) možno odhadnúť ako rozdiel medzi dvoma numerickými riešeniami. Vo (4.9) zanedbáme zvyšok ${\tilde {\cal {O}}}.(h^{p+2})$ s vyššími mocninami $h$, získame odhad lokálnej chyby pomocou numerických riešení dvoch vložených eRKV, $|\varepsilon (t_{i+1},h)|\approx |{\tilde {\eta }}_{i+1}-{\eta }_{i+1}|\ {\mbox{für}}\ h\to 0.$

Zároveň pre metódu rádu konzistencie $p$ platí $|\varepsilon (t_{i+1},h)|\leq Kh^{p+1}\ {\mbox{s konštantou chyby}}\ K=\max _{t\in I}a(t).$

Chybovú konštantu $K$ potom možno vypočítať pre pevnú veľkosť kroku $h$ ako $K\geq {\frac {|\varepsilon (t_{i+1},h)|}{h^{p+1}}}\gtrsim {\frac {|{\tilde {\eta }}_{i+1}-{\eta }_{i+1}|}{h^{p+1}}}.\qquad \qquad (4.10)$

Tento odhad konštánt chyby je základnou myšlienkou riadenia šírky kroku vo vložených Runge-Kutta metódach.
ALGORITHMUS:

1. $t=t_{i}$, zvolíme skúšobnú veľkosť kroku $H$ a toleranciu $tol$ pre lokálnu chybu.
   Vypočítame numerické riešenie v ďalšom kroku $t_{i}+H$ so skúšobnou veľkosťou kroku $H$:
   • pomocou kroku eRKV1 s veľkosťou kroku $H$ a $s$ krokov: $\eta _{i+1}^{H}$,
   • pomocou kroku eRKV2 s veľkosťou kroku $H$ a $m>s$ krokov: ${\tilde {\eta }}_{i+1}^{H}$.
2. Odhadnite konštantu chyby $K$ pomocou (4.10) ako $K\gtrsim {\frac {|{\tilde {\eta }}_{i+1}^{H}-\eta _{i+1}^{H}|}{H^{p+1}}}.$

1. Vypočítajte optimálnu veľkosť kroku $h_{opt}$ tak, aby lokálna chyba $\varepsilon (t_{i},h_{opt})$ približne zodpovedala zadanej tolerancii $tol$, pričom konštanta chyby sa odhaduje podľa vyššie uvedeného: $tol\approx |\varepsilon (t_{i},h_{opt})|\approx Kh_{opt}^{p+1}\gtrsim {\frac {|{\tilde {\eta }}_{i+1}^{H}-\eta _{i+1}^{H}|}{H^{p+1}}}h_{opt}^{p+1}\implies$ $h_{opt}\leq \left({\frac {tol}{K}}\right)^{\frac {1}{p+1}}<\left({\frac {tol}{|{\tilde {\eta }}_{i+1}^{H}-\eta _{i+1}^{H}|}}\right)^{\frac {1}{p+1}}H.$ $\qquad \qquad \qquad (4.11)$
2. Určte (a uložte) numerické riešenie v ďalšom kroku $t_{i}+h_{opt}$ pomocou veľkosti kroku $h_{opt}$ a metódy eRKV2 rádu $p+1$: $\eta _{i+1}:={\tilde {\eta }}_{i+1}^{h_{opt}}={\tilde {\eta }}_{i}+h_{opt}\sum _{j=1}^{m}{\tilde {b}}_{j}k_{j}.$

Ďalšími príkladmi vložených Runge-Kutta metód sú dvojice vzorcov: Dorman $\&$ Prince (4., 5. rádu) alebo Bogacki $\&$ Shampine (2., 3. rádu), ktoré sú známe ako funkcie Matlabu 'ODE45' a 'ODE23'. Runge-Kutta dvojice sa často implementujú ako metóda typu FSAL (First Same As Last), pričom v každom časovom kroku sa uloží prvé vyhodnotenie funkcie $k_{1}$ $f$, pretože to zodpovedá poslednému vyhodnoteniu $k_{s}$ predchádzajúceho kroku, ako je to v nasledujúcej metóde:

Tu platí nasledujúce ${\begin{aligned}k_{4}&=f{\big (}t_{i}+{h},\eta _{i}+h{\big (}{\frac {2}{9}}k_{1}+{\frac {1}{3}}k_{2}+{\frac {4}{9}}k_{3}{\big )}{\big )}\\{\tilde {\eta }}_{i+1}&={\tilde {\eta }}_{i}+h{\Big (}{\frac {2}{9}}k_{1}+{\frac {1}{3}}k_{2}+{\frac {4}{9}}k_{3}{\Big )}\\&\vdots &\ \ t=t_{i+1}\ {\mbox{(Ďalší krok)}}\\k_{1}&=f(t_{i+1},{\eta }_{i+1})=f{\Big (}t_{i}+{h},\eta _{i}+h{\Big (}{\frac {2}{9}}k_{1}+{\frac {1}{3}}k_{2}+{\frac {4}{9}}k_{3}{\Big )}{\Big )}\equiv k_{4}.\end{aligned}}$