Kurs:Gewöhnliche Differentialgleichungen/Eingebettete Runge-Kutta Paare

4.2 Eingebetette Runge-Kutta Paare

Auch mit expliziten Runge-Kutta Verfahren ist die Schrittweitensteuerung möglich. Im Unterschied zu dem in Absatz 3.2 beschriebenem Ansatz für die Schrittweitensteuerung des expliziten Eulerverfahrens, wird hier für die Abschätzung des lokalen Fehlers die Schrittweite nicht halbiert, sondern es werden zwei eRKV mit unterschiedlichen Konsistenzordnungen verwendet. Um die Anzahl der Funktionsauswertungen in der Berechnung der Zwischenstufen $k_j$ so gering wie möglich zu halten, verwenden diese zwei ’eingebetteten’ eRKV gemeinsame Stützstellen. Hierbei ist zu erwähnen, dass nach den Erkentnissen über die Stufenzahl und Ordnung eines eRKV das Verfahren höherer Ordnung weitere zusätzliche Zwischenstellen benötigt, um diese höhere Ordnung zu erreichen.

Als eingebetettes Runge-Kutta Paar verstehen wir zwei eRKV mit $s$ gemeinsamen Stützstellen (und Stufen). Das erste der eRKV hat dabei die Konsistenzordnung $p$ und das zweite Verfahren erreicht mit den $s$ gemeinsamen Stützstellen und den zusätzlichen $s+1,s+2,\dots ,m$ Stützsstellen, die Konsistenzordnung $q>p$. Die Gewichtsvektoren dieser beider eRKV sind im Allgemeinen unterschiedlich und die Iterationsvorschrifte (die Runge-Kutta-Updates) der beiden Verfahren sind: $$\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}$$ Wegen den s gemeinsamen Stufen werden in beiden Verfahren die Funktionsauswertungen $k_1,\dots k_s$ verwendet. Ein Beispiel eines eingebetetten Runge-Kutta Paares ist das
Fehlenberg RK-Paar der Ordnung 2 und 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 Schrittweitensteuerung bei eingebetetten Runge-Kutta Verfahren

Im Folgenden beschränken wir uns auf den Fall $q=p+1$ und leiten die Formel der Schrittweitensteuerung her. Seien ${\eta }_{i+1},{\tilde{\eta }}_{i+1}$ die numerischen Lösungen in $t_{i+1}$, die aus der exakten Lösung $y(t_i)$ statt ${\eta }_i,{\tilde{\eta }}_i$ ausgehen und die Verfahren haben die Konsistenzordnung $p$ und $p+1$. Dann gilt für die lokalen Fehler nach Anwendung der Taylorentwicklung (mit Taylorrest) $$\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}$ nennen wir die führenden Fehlerglieder der beiden Verfahren. Diese sind die Glieder mit der kleinsten Potenz von $h$ und enthalten bis zu $(p+1)$-te und $(p+2)$-te partielle Ableitungen der Funktion $f$ in $t_i$. Nach dem Subtrahieren der beiden Gleichungen erhalten wir ${\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)$

d.h. der lokale Fehler des ersten Verfahrens (eRKV1) lässt sich als Unterschied der beiden numerischen Lösungen abschätzen. Vernachlässigen wir in (4.9) den Rest $\widetilde{𝒪}(h^{p+2})$ mit höheren Potenzen von $h$, erhalten wir die Abschätzung des lokalen Fehlers mithilfe der numerischen Lösungen der zwei eingebetteten eRKV, $$|\epsilon (t_{i+1},h)|\approx |{\tilde{\eta }}_{i+1}-{\eta }_{i+1}|\text{für}h\to 0.$$

Gleichzeitig gilt für ein Verfahren der Konsistenzordnung $p$, dass $$|\epsilon (t_{i+1},h)|\le K{h}^{p+1}\text{mit Fehlerkonstante}K=\underset{t\in I}{\max }a(t).$$

Danach lässt sich die Fehlerkonstante $K$ bei einer festen Schrittweite $h$ berechnen als $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)$

Diese Abschätzung der Fehlerkonstanten ist die Grundidee der Schrittweitensteuerung bei eingebetteten Runge-Kutta Verfahren.

ALGORITHMUS:

1. $t=t_i$, wähle eine Probeschrittweite $H$ und die Toleranz $tol$ für den lokalen Fehler.
   Berechne die numerische Lösung im nächsten Schritt $t_i+H$ mit Probeschrittweite $H$:
   • mithilfe eines Schrittes des eRKV1 mit Schrittweite $H$ und $s$ Stufen: ${\eta }_{i+1}^H$,
   • mithilfe eines Schrittes des eRKV2 mit Schrittweite $H$ und $m>s$ Stufen: ${\tilde{\eta }}_{i+1}^H$.
2. Schätze die Fehlerkonstante $K$ anhand (4.10) ab, als $K\gtrsim \frac{|{\tilde{\eta }}_{i+1}^H-{\eta }_{i+1}^H|}{H^{p+1}}.$
3. Berechne die optimale Schrittweite $h_{opt}$, sodass der lokale Fehler $\epsilon (t_i,h_{opt})$ etwa der vorgegebenen Toleranz $tol$ entspricht, wobei die Fehlerkonstante wie oben geschätzt wird: $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)$
4. Bestimme (und speichere) die numerische Lösung im nächsten Schritt $t_i+h_{opt}$ mithilfe der Schrittweite $h_{opt}$ und des Verfahrens eRKV2 der Ordnung $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.$$

Weitere Beispiele der eingebetteten Runge-Kutta Verfahren sind die Formelpaare: Dorman $\mathrm{\&}$ Prince (4, 5 Ordnung), oder Bogacki $\mathrm{\&}$ Shampine (2, 3 Ordnung), die als Matlab-Funktionen ’ODE45’ und ’ODE23’ bekannt sind. Oft werden die Runge-Kutta Paare als FSAL-Typ (First Same As Last) Verfahren implementiert, wobei in jedem Zeitschritt die erste Funktionsauswertung $k_1$ von $f$ gespart wird, da diese mit der letzten Auswertung $k_s$ des vorherigen Schrittes übereinstimmt, wie es in folgendem Verfahren der Fall ist:

Hier gilt $$\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{(nächstes Schritt)}\\ 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}$$