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,\dots ,m$ opornými bodmi a ďalšími $s+1,s+2,\dots ,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{array}{cccc}& & \text{eRKV1:}{\eta }_{i+1}& ={\eta }_i+h\sum _{j=1}^s{b}_j{k}_j\\ & & \text{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{array}$ Kvôli spoločným krokom s sa v oboch metódach používajú vyhodnotenia funkcií $k_1,\dots k_s$. Príkladom vloženej dvojice Runge-Kutta je
Fehlenbergova RK dvojica rádov 2 a 3:

$\begin{array}{c}{k}_1=f(t_i,{\eta }_i)\\ k_2=f(t_i+h,{\eta }_i+h{k}_1)\\ k_3=f(t_i+\frac{h}{2},{\eta }_i+\frac{h}{4}{k}_1+\frac{h}{4}{k}_2)\\ {\eta }_{i+1}={\eta }_i+h(\frac{1}{2}{k}_1+\frac{1}{2}{k}_2)\\ {\tilde{\eta }}_{i+1}={\tilde{\eta }}_i+h(\frac{1}{6}{k}_1+\frac{1}{6}{k}_2+\frac{2}{3}{k}_3).\end{array}$

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{array}{ccc}& & \epsilon (t_{i+1},h)=h\tau (t_{i+1},h)=y(t_{i+1})-{\eta }_{i+1}=a(t_i)h^{p+1}+𝒪(h^{p+2})=a({\theta }_1)h^{p+1},\\ & & \tilde{\epsilon }(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}+𝒪(h^{p+3})=b({\theta }_2)h^{p+2},\end{array}$

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}+\widetilde{𝒪}(h^{p+2})\approx h\tau (t_{i+1},h),$ $(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 $\widetilde{𝒪}.(h^{p+2})$ s vyššími mocninami $h$, získame odhad lokálnej chyby pomocou numerických riešení dvoch vložených eRKV, $|\epsilon (t_{i+1},h)|\approx |{\tilde{\eta }}_{i+1}-{\eta }_{i+1}|\text{für}h\to 0.$

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

Chybovú konštantu $K$ potom možno vypočítať pre pevnú veľkosť kroku $h$ ako $K\ge \frac{|\epsilon (t_{i+1},h)|}{h^{p+1}}\gtrsim \frac{|{\tilde{\eta }}_{i+1}-{\eta }_{i+1}|}{h^{p+1}}.(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 $\epsilon (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 |\epsilon (t_i,h_{opt})|\approx K{h}_{opt}^{p+1}\gtrsim \frac{|{\tilde{\eta }}_{i+1}^H-{\eta }_{i+1}^H|}{H^{p+1}}{h}_{opt}^{p+1}⟹$ $h_{opt}\le {\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.$ $(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 $\mathrm{\&}$ Prince (4., 5. rádu) alebo Bogacki $\mathrm{\&}$ 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{array}{cc}{k}_4& =f(t_i+h,{\eta }_i+h(\frac{2}{9}{k}_1+\frac{1}{3}{k}_2+\frac{4}{9}{k}_3\left)\right)\\ {\tilde{\eta }}_{i+1}& ={\tilde{\eta }}_i+h(\frac{2}{9}{k}_1+\frac{1}{3}{k}_2+\frac{4}{9}{k}_3)\\ & ⋮& t=t_{i+1}\text{(Ďalší krok)}\\ k_1& =f(t_{i+1},{\eta }_{i+1})=f(t_i+h,{\eta }_i+h(\frac{2}{9}{k}_1+\frac{1}{3}{k}_2+\frac{4}{9}{k}_3\left)\right)\equiv k_4.\end{array}$