Kurs:Gewöhnliche Differentialgleichungen/4.2 Eingebetette 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 ${\textstyle 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 ${\textstyle s}$ gemeinsamen Stützstellen (und Stufen). Das erste der eRKV hat dabei die Konsistenzordnung ${\textstyle p}$ und das zweite Verfahren erreicht mit den ${\textstyle s}$ gemeinsamen Stützstellen und den zusätzlichen ${\textstyle s+1,s+2,\ldots ,m}$ Stützsstellen, die Konsistenzordnung ${\textstyle q>p}$. Die Gewichtsvektoren dieser beider eRKV sind im Allgemeinen unterschiedlich und die Iterationsvorschrifte (die Runge-Kutta-Updates) der beiden Verfahren sind:

$${\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}}$$

Wegen den s gemeinsamen Stufen werden in beiden Verfahren die Funktionsauswertungen ${\textstyle k_{1},\ldots k_{s}}$ verwendet. Ein Beispiel eines eingebetetten Runge-Kutta Paares ist das
Fehlenberg RK-Paar der Ordnung 2 und 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 Schrittweitensteuerung bei eingebetetten Runge-Kutta Verfahren

Im Folgenden beschränken wir uns auf den Fall ${\textstyle q=p+1}$ und leiten die Formel der Schrittweitensteuerung her. Seien ${\textstyle {\eta }_{i+1},{\tilde {\eta }}_{i+1}}$ die numerischen Lösungen in ${\textstyle t_{i+1}}$, die aus der exakten Lösung ${\textstyle y(t_{i})}$ statt ${\textstyle \eta _{i},{\tilde {\eta }}_{i}}$ ausgehen und die Verfahren haben die Konsistenzordnung ${\textstyle p}$ und ${\textstyle p+1}$. Dann gilt für die lokalen Fehler nach Anwendung der Taylorentwicklung (mit Taylorrest)

$${\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 ${\textstyle \theta _{1},\theta _{2}\in [t_{i},t_{i}+h]}$. ${\textstyle 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 ${\textstyle h}$ und enthalten bis zu ${\textstyle (p+1)}$-te und ${\textstyle (p+2)}$-te partielle Ableitungen der Funktion ${\textstyle f}$ in ${\textstyle t_{i}}$. Nach dem Subtrahieren der beiden Gleichungen erhalten wir ${\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)$

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 ${\textstyle {\tilde {\cal {O}}}(h^{p+2})}$ mit höheren Potenzen von ${\textstyle h}$, erhalten wir die Abschätzung des lokalen Fehlers mithilfe der numerischen Lösungen der zwei eingebetteten eRKV,

$$|\varepsilon (t_{i+1},h)|\approx |{\tilde {\eta }}_{i+1}-{\eta }_{i+1}|\ {\mbox{für}}\ h\to 0.$$

Gleichzeitig gilt für ein Verfahren der Konsistenzordnung ${\textstyle p}$, dass

$$|\varepsilon (t_{i+1},h)|\leq Kh^{p+1}\ {\mbox{mit Fehlerkonstante}}\ K=\max _{t\in I}a(t).$$

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

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

ALGORITHMUS:

1. ${\textstyle t=t_{i}}$, wähle eine Probeschrittweite ${\textstyle H}$ und die Toleranz ${\textstyle tol}$ für den lokalen Fehler.
   Berechne die numerische Lösung im nächsten Schritt ${\textstyle t_{i}+H}$ mit Probeschrittweite ${\textstyle H}$:
   • mithilfe eines Schrittes des eRKV1 mit Schrittweite ${\textstyle H}$ und ${\textstyle s}$ Stufen: ${\textstyle \eta _{i+1}^{H}}$,
   • mithilfe eines Schrittes des eRKV2 mit Schrittweite ${\textstyle H}$ und ${\textstyle m>s}$ Stufen: ${\textstyle {\tilde {\eta }}_{i+1}^{H}}$.
2. Schätze die Fehlerkonstante ${\textstyle K}$ anhand (4.10) ab, als ${\textstyle K\gtrsim {\frac {|{\tilde {\eta }}_{i+1}^{H}-\eta _{i+1}^{H}|}{H^{p+1}}}.}$
3. Berechne die optimale Schrittweite ${\textstyle h_{opt}}$, sodass der lokale Fehler ${\textstyle \varepsilon (t_{i},h_{opt})}$ etwa der vorgegebenen Toleranz ${\textstyle tol}$ entspricht, wobei die Fehlerkonstante wie oben geschätzt wird: $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)$
4. Bestimme (und speichere) die numerische Lösung im nächsten Schritt ${\textstyle t_{i}+h_{opt}}$ mithilfe der Schrittweite ${\textstyle h_{opt}}$ und des Verfahrens eRKV2 der Ordnung ${\textstyle 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 ${\textstyle \&}$ Prince (4, 5 Ordnung), oder Bogacki ${\textstyle \&}$ 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 ${\textstyle k_{1}}$ von ${\textstyle f}$ gespart wird, da diese mit der letzten Auswertung ${\textstyle k_{s}}$ des vorherigen Schrittes übereinstimmt, wie es in folgendem Verfahren der Fall ist:

Hier gilt

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