Die Mathe-Redaktion - 22.08.2017 14:59 - Registrieren/Login
Auswahl
Schwarzes Brett
Aktion im Forum
Suche
Stichwortsuche in Artikeln und Links von Matheplanet
Suchen im Forum
Suchtipps

Bücher
Englische Bücher
Software
Suchbegriffe:
Mathematik bei amazon
Naturwissenschaft & Technik
In Partnerschaft mit Amazon.de
Kontakt
Mail an Matroid
[Keine Übungsaufgaben!]
Impressum

Bitte beachten Sie unsere Nutzungsbedingungen, die Distanzierung, unsere Datenschutzerklärung und
die Forumregeln.

Sie können Mitglied werden oder den Newsletter bestellen.

Der Newsletter Apr. 2017

Für Mitglieder
Mathematisch für Anfänger
Wer ist Online
Aktuell sind 758 Gäste und 18 Mitglieder online.

Sie können Mitglied werden:
Klick hier.

Über Matheplanet
 
Mathematik: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
Freigegeben von matroid am So. 03. Juli 2016 20:56:22
Verfasst von Higlav -   950 x gelesen [Gliederung] [Statistik] Druckbare Version Druckerfreundliche Version
Mathematik



Vorwort

Im Rahmen einer kleineren Projektarbeit zur Verbesserung meines Notenschnittes in meinem Numerik-Modul entwickelte ich eher am Rande und per Zufall eine Schrittweitensteuerung für eingebettete Runge-Kutta-Verfahren, welche eine gewisse Verbesserung zum klassischen Algorithmus bietet. In diesem Artikel werde ich diese Optimierung vorstellen.
Vielleicht ist "Verbesserung" etwas unglücklich gewählt. Bei Bedarf ändere ich es auch auf "Modifikation".



Navigation

  1. Einführung
  2. Grundlagen
    1. Schrittweitensteuerung
    2. Eingebettete Verfahren
  3. Die verbesserte Schrittweitenabschätzung
  4. Programmierung
  5. Vergleich zum herkömmlichen Algorithmus
  6. Fazit




1: Einführung

So manch physikalischer Zusammenhang lässt sich mathematisch am einfachsten anhand einer Differentialgleichung beschreiben. Beispielsweise lässt sich ein freier Fall leichter mit einer gewöhnlichen Differentialgleichung(nachfolgend "ODE"="ordinary differential equation") beschreiben als gleich mit einer kompletten Bahnfunktion:

Hier gilt Folgendes: <math>m_2</math>-Attraktor ist eine Kugel mit Radius <math>r</math>, auf deren Oberfläche der Nullpunkt des Koordinatensystems gewählt wurde. Der Einfluss der <math>x</math>- und <math>y</math>-Koordinaten auf den Abstand der beiden Schwerpunkte sowie die Massenveränderung durch die Relativitätstheorie werden nicht berücksichtigt.
Beobachtung: Beschleunigung in <math>z</math>-Richtung <math>\approx-G\frac{m_1+m_2}{m_1\cdot(r+z)^2}</math> - bzgl. der anderen Koordinaten bleibt der Impuls gleich dem Startimpuls.
Differentialgleichung: <math>z''=-G\frac{m_1+m_2}{m_1\cdot(r+z)^2},\quad x'=v_x(0),\quad y'=v_y(0)</math>
Lösungsfunktion: <math>z(t)=-G\frac{m_1+m_2}{2m_1\cdot(r+z(t))^2}t^2+v_z(0)\cdot t+z(0),\quad x(t)=v_x(0)\cdot t,\quad y(t)=v_y(0)\cdot t</math>

Die Lösungsfunktion wird nach dem Notieren der Differentialgleichung anhand schwarzer Magie der Analysis klassifiziert und gelöst. Jedoch gestalten sich viele reale Probleme als solche, die zu lösen kaum möglich sind. In diesem Falle ist eine Approximation die schärfste Klinge der Alternativen. Hierbei wird anhand eines geeigneten Verfahrens, welches auf die ODE angewandt wird, die Lösung approximiert. Bei dieser Approximation stellen sich jedoch einige Fragen:

  • Wie gross ist der Fehler?
  • Wie genau kann approximiert werden?
  • Welches Verfahren ist geeignet?
  • Welche Parameter müssen wie geregelt werden?

2: Grundlagen

Es gibt verschiedene Runge-Kutta-Verfahren, die sich in Fehlerordnung und ihrem jeweils charakteristischen Butcher-Tableau unterscheiden. Letzteres sagt alles über die jeweilige Methode aus und ist wie folgt definiert:

Je nach Quelle sind <math>c</math> und <math>b</math> miteinander vertauscht. Die Bezeichnungen hier entsprechen der Darstellung in Karl Strehmels "Numerik gewöhnlicher Differentialgleichungen : Nichtsteife, steife und differential-algebraische Gleichungen"(2. Auflage, 2012, 978-3-8348-1847-8) sowie derjenigen auf Wikipedia.
<math>     \begin{equation}\nonumber     \begin{array}{c|c}     c&A\\\hline &b^{{}_T}     \end{array} = \begin{array}{c|cccc}     c_1&a_{11}&a_{12}&\hdots&a_{1s}\\     c_2&a_{21}&a_{22}&\hdots&a_{2s}\\     \vdots&\vdots&\vdots&\ddots&\vdots\\     c_s&a_{s1}&a_{s2}&\hdots&a_{ss}\\     \hline     &b_1&b_2&\hdots&b_s     \end{array}     \end{equation}     </math>

Die Matrix <math>A</math> ist bei expliziten Verfahren eine untere Dreiecksmatrix mit der Hauptdiagonale gleich 0. <math>s</math> heisst die Stufe des Verfahrens. Es gibt zu jeder Fehlerordnung eine entsprechende Stufe <math>s</math>:


<math>     \begin{tabular}{r||cccccccc|l}     Ordnung p&1&2&3&4&5&6&7&8&p\geq9\\\hline     Stufe s&1&2&3&4&6&7&9&11&s\geq p+3     \end{tabular}     </math>

Als Beispiel sei hier das Butcher-Tableau der Heun-Methode aufgeführt:

< <math>     \begin{equation}\nonumber     B_{Heun} = \begin{array}{c|cc}     0&0&0\\     1&1&0\\     \hline     &\frac12&\frac12     \end{array}     \end{equation}     </math>

Das Vorgehen bei der Berechnung eines Integrationsschrittes ist wie folgt: Zuerst werden die Vektoren <math>b</math> und <math>c</math> und die Matrix <math>A</math> benötigt. Dies sei hier am vorherigen Beispiel des Heun-Verfahrens demonstriert:

<math>     \begin{equation}\nonumber     c=\begin{pmatrix}0\\1\end{pmatrix},\qquad b=\begin{pmatrix}\frac12\\\frac12\end{pmatrix},\qquad A=\begin{pmatrix}0&0\\1&0\end{pmatrix}     \end{equation}     </math>

Nun errechnet sich der nächste Wert mit der Formel

<math>     \begin{equation}\nonumber     x_{k+1}=x_k+h\sum_{j=1}^sb_jk_j,\qquad k_j=f\left(t_k+h\cdot c_j,\ x_k+h\sum_{l=1}^ja_{jl}k_l\right)     \end{equation}     </math>

Daraus ergibt sich in Fall des Heun-Verfahrens mit <math>s=|b|=|c|=2</math> dann


<math>     \begin{array}{rl}\nonumber     x_{k+1}&=x_k+h(b_1k_1+b_2k_2)\\\nonumber     k_1&=f(t_k+h\cdot c_1,\ x_k+h\cdot c_1(a_{11}k_1))\\\nonumber     k_2&=f(t_k+h\cdot c_2,\ x_k+h\cdot c_2(a_{21}k_1+a_{22}k_2))\\\nonumber\Leftrightarrow     k_1&=f(t_k+h\cdot0,\ x_k+h\cdot0(0\cdot k_1))\\\nonumber     k_2&=f(t_k+h\cdot1,\ x_k+h\cdot1(1\cdot k_1+0\cdot k_2))     \end{array}     </math>

Hier wird auch gleich ersichtlich, wieso die Hauptdiagonale gleich 0 sein muss. Denn ist sie es nicht, ergäbe dies eine Selbstreferenz der Zwischenwerte <math>k_j</math>. Bei den impliziten Runge-Kutta-Verfahren ist dies wiederum anders, aber das ist nicht Gegenstand dieses Artikels.
Da es sich hierbei jedoch immer um Näherungen handelt, sind die errechneten Werte fehlerbehaftet. Dieser Fehler wird wie folgt dargestellt:

<math>     \begin{equation}\nonumber     \mathcal O\left(h^{p+1}\right)=C\cdot h^{p+1}=\alpha_1h^{p+1}+\alpha_2h^{p+2}+\alpha_3h^{p+3}+\cdots     \end{equation}     </math>

Dabei ist der Koeffizientenvektor <math>\vec\alpha</math> unbekannt. Charakteristisch für dieses Fehlerpolynom ist <math>\alpha_1h^p</math>, wenn <math>h<1</math>. Für die Approximationen wird dieses Glied der niedersten Ordnung betrachtet. Da aber <math>\vec\alpha</math> sich bei veränderter Fehlerordnung <math>p</math> ebenfalls verändert, heisst das, dass ein Verfahren höherer Ordnung nicht zwangsläufig den kleineren Fehler hat.


2.1: Schrittweitensteuerung

Nun ergibt sich bei ODE's, die stark schwanken sowie bei ODE's, die eine sehr flache Topologie haben ein Problem: Entweder ist die gewählte, feste Schrittweite zu gross oder zu klein. Welcher Fall vorliegt kann ohne Analyse der ODE nicht festgestellt werden. Aus diesem Grund bietet es sich aus Performance- und Genauigkeitsgründen an, die Schrittweite dynamisch zu verändern. Im optimalsten Fall lässt sich der LDE und/oder der GDE unter einem bestimmten Wert halten und so die nächstbeste Schrittweite schätzen.
Nun ist diese Vorgehensweise keine richtige Lösung, sondern lediglich eine Approximation. Diese Approximation kommt mit einem unbestimmbaren(sofern die exakte Lösung nicht gegeben ist) Fehler <math>\varepsilon</math>:

<math>     \begin{tikzpicture}[scale=10]         \newcommand\Xe{1}         \newcommand\Ye{1}         \newcommand\dx{0.02}         \clip (-3*\dx, -5*\dx) rectangle ({\Xe+7*\dx}, {\Ye+3*\dx});         \draw[color=black!20] (-\dx, -\dx) grid ({\Xe+\dx}, {\Ye+\dx});         \draw[->, thick] (-\dx, 0) -- ({\Xe+\dx}, 0) node[anchor=north west]{$t$};         \draw[->, thick] (0, -\dx) -- (0, {\Ye+\dx}) node[anchor=south east]{$x$};          \newcommand\ex{2/3}         \newcommand\X{0.5}         \newcommand\h{0.1}         \newcommand\F[1]{-0.0217374*(#1)^8+0.465713*(#1)^7-3.99*(#1)^6+17.594*(#1)^5-42.11*(#1)^4+53.48*(#1)^3-32.4*(#1)^2+8*(#1)}         \newcommand\f[1]{-0.1738992*(#1)^7+3.259991*(#1)^6-23.94*(#1)^5+87.97*(#1)^4-168.44*(#1)^3+160.44*(#1)^2-64.8*(#1)+8}          \fill[fill=red!40,draw=none,variable=\x,samples=50,fill opacity=0.4]         plot[domain=\X:\Xe] ({\x},        {\F{\x}       +exp(\ex*(\x-\X))-1}) --         plot[domain=\X:\Xe] ({\X+\Xe-\x}, {\F{\X+\Xe-\x}-exp(\ex*(\Xe-\x))+1});         \draw[draw=green,variable=\x,samples=50]         plot[domain=\X:\Xe] ({\x},        {\F{\x}       +exp(\ex*(\x-\X))-1});         \draw[draw=green,variable=\x,samples=50]         plot[domain=\X:\Xe] ({\X+\Xe-\x}, {\F{\X+\Xe-\x}-exp(\ex*(\Xe-\x))+1});          \draw[color=red!70!yellow,variable=\x,samples=50] plot[domain=0:\Xe] ({\x},{\F{\x}});         \draw (0.23, {\F{0.23}}) node[anchor=south]{$x(t)$};          \draw[color=blue] (\X, {\F{\X}}) circle ({\dx/2.5});         \draw[color=blue,->] (\X, {\F{\X}}) -- ({\X+\h},{\F{\X}+\h*(\f{\X})});          \draw ({\X+\h}, {\F{\X+\h}})         --++ ({\Xe-\X-\h+2*\dx}, 0)         node[pos=1,anchor=south west]{$x_{k+1}$}         --++ (0, {-(\F{\X+\h}-(\F{\X}+\h*(\f{\X})))})         node[pos=0.5,anchor=west]{$\textcolor{red}\varepsilon$}         --+ ({-(\Xe-\X-\h+2*\dx)}, 0)         node[pos=0,anchor=north west]{$\tilde x_{k+1}$}         -- ({\X+\h}, -2*\dx)         node[pos=1,anchor=north west]{$t_{k+1}$}         -- (\X, -2*\dx)         node[pos=0.5,anchor=north]{$h$}         -- (\X, {\F{\X}})         node[pos=0,anchor=north east]{$t_k$};     \end{tikzpicture}     </math>

Dieser lokale Fehler(LDE) summiert sich vorzeichenbehaftet über das abgeschlossene Integrationsintervall zu eine globalen Fehler(GDE) auf. Beim Vergleich zur Taylorreihe ergeben sich für die Runge-Kutta-Methoden folgendes Restglied, welches dem Fehler entspricht(<math>p</math> sei die Fehlerordnung des jeweiligen Verfahrens):

<math>     \begin{equation}\nonumber     \varepsilon=\mathcal O(h^p)     \end{equation}     </math>

Hierbei ist <math>\mathcal O(h^p)</math> ein Polynom von <math>h^{p+n}</math>. Die Koeffizienten dieses Polynoms ist dabei abhängig von der ODE selbst an der jeweiligen Stellen sowie den Werten. Das heisst, dieser Fehler variiert je nach Topologie und eingesetzten Werten. Demnach kann sich der Fehler noch weiter vergrössern, wenn "falsche" Werte - Approximationen anstatt den exakten Werten - eingesetzt werden.


2.2: Eingebettete Verfahren

Die eingebetteten Verfahren berechnen jeden Schritt einmal mit einem Verfahren der Ordnung <math>p</math> und einmal mit einem der Ordnung <math>p+1</math>.
Eingebettete RK-Verfahren werden mit dem erweiterten Butcher-Tableau beschrieben und sind in der Form:

<math>     \begin{equation}\nonumber         \begin{array}{c|c}             \vec c & A\\\hline             & \vec b^{\ T}\\             & \vec b^{\ *T}         \end{array}     \end{equation}     </math>

Dabei steht der Vektor <math>\vec b^{\ T}</math> für das Verfahren der Fehlerordnung <math>p</math> und der Vektor <math>\vec b^{\ *T}</math> für dasjenige der Fehlerordnung <math>q=p+1</math>.
Für die Berechnung der beiden Approximationswerte bleiben die <math>k</math>-Werte gleich und werden dann entweder mit dem <math>\vec b^{\ T}</math>- oder dem <math>\vec b^{\ *T}</math>-Vektor linear kombiniert. Betrachtet man den Fehler der beiden Werte zeigt sich folgendes(<math>\textcolor{green}{\hat\circ}</math> = Verfahren höherer Ordnung, <math>\textcolor{blue}{\tilde\circ}</math> = Verfahren niederer Ordnung):

<math>     \newcommand\ape[1]{\textcolor{green}{#1}}     \newcommand\apz[1]{\textcolor{blue}{#1}}     \begin{array}{cc}     x_{k+1}&=\apz{\hat x_{k+1}}+\mathcal O(h^{p+2})\\     x_{k+1}&=\ape{\tilde x_{k+1}}+\mathcal O(h^{p+1})     \end{array}     </math>

Daraus ergibt sich bei Betrachtung des führenden Fehlertermes(demjenigen niederster Ordnung bei <math>h<1</math>) folgendes Bild:

<math>     \begin{tikzpicture}[scale=5]         \newcommand\ape[1]{\textcolor{green}{#1}}         \newcommand\apz[1]{\textcolor{blue}{#1}}         \newcommand\Hk{0.55}         \newcommand\Hs{0.45}         \draw[color=black!20] (0, 0) grid (1, 1);         \draw[->,thick] (0, 0) node[anchor=north east]{$0$} -- (1, 0) node[anchor=north west]{$h$};         \draw[->,thick] (0, 0) -- (0, 1) node[anchor=south east]{$\varepsilon$};         \fill[fill=red!5,draw=orange!50!red,domain=0:1,variable=\x,samples=20] plot (\x, \x^2-\x^3) -- cycle;         \draw[->] (1.1, 0.3) node[anchor=west]{$e(h)=\left|\ape{bh^{p+1}}-\apz{ah^{p+2}}\right|$} -- (0.7, 0.1);         \draw[color=blue,domain=0:1,variable=\x,samples=20] plot (\x, \x^3);         \draw (0.8, 0.7) node[anchor=south east]{$\ape{bh^{p+1}}$};         \draw[color=green!70!blue,domain=0:1,variable=\x,samples=20] plot (\x, \x^2);         \draw (0.75, 0.5) node[anchor=north west]{$\apz{ah^{p+2}}$};         \draw[dotted] (\Hk, 0) node[anchor=north]{$h_k$} --++ (0, \Hk^2) --+ (-\Hk, 0) node[anchor=east]{$\ape{\tilde\varepsilon}$};         \draw[dotted] (\Hk, \Hk^3) --+ (-\Hk, 0) node[anchor=south east]{$\apz{\hat\varepsilon}$};         \draw[dotted] (\Hs, 0) node[anchor=north]{$h_k^*$} --++ (0, \Hs^3) --+ (-\Hs, 0) node[anchor=east]{$\varepsilon_0^*$};     \end{tikzpicture}     </math>

Und um den Fehler abschätzen zu können, macht man eine Abschätzung über die Differenz der beiden Fehlerglieder niederster Ordnung:

<math>     \newcommand\ape[1]{\textcolor{green}{#1}}     \newcommand\apz[1]{\textcolor{blue}{#1}}     \begin{equation}\nonumber         \varepsilon=\left|\apz{\hat x_{k+1}}-\ape{\tilde x_{k+1}}\right|     \end{equation}     </math>

Das heisst, dass der Fehler um eine weitere ganze Ecke geschätzt wird und deshalb ein schlechterer Schätzer ist.
Damit lässt sich nun die neue Schrittweite abschätzen:

<math>     \newcommand\ape[1]{\textcolor{green}{#1}}     \newcommand\apz[1]{\textcolor{blue}{#1}}     \begin{equation}\nonumber         h_k^*=h_k\beta\cdot\left\{\begin{matrix}         \sqrt[\underline{p+}1]{\frac{\varepsilon_0}\varepsilon},\qquad&\varepsilon\ge\varepsilon_0=\varepsilon_\mathrm{abs}+\max{\left(\apz{\hat x_{k+1}},\ape{\tilde x_{k+1}}\right)}\cdot\varepsilon_\mathrm{rel}\\         \sqrt[\underline{p+}2]{\frac{\varepsilon_0}\varepsilon},\qquad&\varepsilon<\varepsilon_0=\varepsilon_\mathrm{abs}+\max{\left(\apz{\hat x_{k+1}},\ape{\tilde x_{k+1}}\right)}\cdot\varepsilon_\mathrm{rel}         \end{matrix}\right.     \end{equation}     </math>

Der Verfahrweg sieht vereinfacht so aus:

Der Ablauf entspricht nicht ganz dem später implementierten oder der in Matlab eingebauten Lösung.
<math>     \newcommand\ape[1]{\textcolor{green!70!blue}{#1}}     \newcommand\apz[1]{\textcolor{blue}{#1}}     \begin{tikzpicture}[scale=1, text height=0.4cm,         every node/.style={column sep=0.5cm, row sep=1cm, thick, align=center},         decision/.style = {diamond, text centered,         inner sep=1pt, rounded corners,         top color=white,         bottom color=blue!80!cyan!30, draw=blue!80!cyan},         block/.style    = {rectangle, text centered,         minimum height=2em, rounded corners,         top color=white,         bottom color=orange!80!yellow!30, draw=orange!80!yellow},         line/.style     = {draw, ->, shorten >=2pt}]          \node (Start) at (0, 0) {$t_0,t_e,x_0,h_0,\varepsilon_0,ode(t,x,h)^p$};         \path (Start) ++ (0, -2) node(Initial)[block, text height=1.2cm] {$x_k:=x_0$\\$t_k:=t_0$\\$h_k:=h_0$};         \path (Initial) ++ (0, -3) node(FrEnde)[decision] {$t_k\leq t_e$};         \path (FrEnde) ++ (0, -3) node(Iter)[block, text height=1cm] {$\apz{\hat x_{k+1}}=ode\left(t_k,x_k,h_k\right)^{p+1}$\\         $\ape{\tilde x_{k+1}}=ode\left(t_k,x_k,h_k\right)^p$};         \path (Iter) ++ (0, -3) node(Fehler)[block] {$\varepsilon=\left|\apz{\hat x_{k+1}}-\ape{\tilde x_{k+1}}\right|$};         \path (Fehler) ++ (0, -2.5) node(Entsch)[decision] {$\varepsilon<\varepsilon_0$};         \path (Entsch) ++ (0, -2) node(DecJa)[block] {$h_{k+1}=h_k\beta\cdot\mathrm{min}\left\lbrace\alpha_\mathrm{max},\sqrt[p+1]{\frac{\varepsilon_0}\varepsilon}\right\rbrace$};         \path (Entsch) ++ (5, 0) node(DecNein)[block] {$h_k:=h_k\beta\cdot\mathrm{max}\left\lbrace\alpha_\mathrm{min},\sqrt[p+2]{\frac{\varepsilon_0}\varepsilon}\right\rbrace$};         \path (DecJa) ++ (0, -2) node(SchrNext)[block, text height=0.8cm]{$x_{k+1}=\apz{\hat x_{k+1}}$\\$t_{k+1}=t_k+h_k$};         \path (FrEnde) ++ (-3.5, -8) node(Reset)[block, text height=1.2cm] {$x_k:=x_{k+1}$\\$t_k:=t_{k+1}$\\$h_k:=h_{k+1}$};         \path (DecNein) ++ (0, -2) coordinate(knot);         \path (FrEnde) ++ (2.5, 0) node(Ende){Ende};         \path (Iter) --++ (5, 0) coordinate(Reentry);          \draw[->,font=\scriptsize,thick]         (Start) edge (Initial)         (Initial) edge node[anchor=west]{Iterationsstart} (FrEnde)         (FrEnde) edge node[near start, anchor=west]{Ja} (Iter)         (FrEnde) edge node[near start, anchor=south]{Nein} (Ende)         (Iter) edge node[anchor=west]{Fehlerschätzung} (Fehler)         (Fehler) edge (Entsch)         (Entsch) edge node[anchor=south east]{Ja} (DecJa)         (Entsch) edge node[anchor=south west]{Nein} (DecNein)         (DecJa) edge (SchrNext)         (DecNein) -- (Reentry) edge (Iter)         (SchrNext) -| (Reset) --++ (0, 8) -- (FrEnde)         ;          \node(leyend) at (8, -5.5){         \begin{tabular}{>{\sffamily}r@{ : }p{7cm}}             \multicolumn{2}{c}{\textbf{Variablen}} \\             $ode(t,x,h)^p$ & Iteration an $t$ und $x$ mit bestimmter Runge-Kutta-Methode der Ordnung $p$\\             $\apz{\hat\circ}$ & Verfahren höherer Ordnung\\             $\ape{\tilde\circ}$ & Verfahren niederer Ordnung\\             $\alpha_\mathrm{min}$ & $\in[0.2,0.5]$, untere Sicherheitsschranke\\             $\alpha_\mathrm{max}$ & $\in[1.5,2]$, obere Sicherheitsschranke\\             $\beta$ & $\in[0.8,0.95]$, Sicherheitsfaktor         \end{tabular}         };     \end{tikzpicture}     </math>

3: Die verbesserte Schrittweitenabschätzung

Im Laufe der Arbeit wurde eine andere Herangehensweise an das Problem der Fehlerschätzung erarbeitet. Dabei wird der Fehler nicht über einen Umweg geschätzt sondern lediglich mit der Annahme, dass <math>a\approx b\approx c</math> gilt.

<math>     \begin{tikzpicture}[scale=5]         \newcommand\ape[1]{\textcolor{green}{#1}}         \newcommand\apz[1]{\textcolor{blue}{#1}}         \newcommand\Hk{0.55}         \newcommand\Hs{0.45}         \draw[color=black!20] (0, 0) grid (1, 1);         \draw[->,thick] (0, 0) node[anchor=north east]{$0$} -- (1, 0) node[anchor=north west]{$h$};         \draw[->,thick] (0, 0) -- (0, 1) node[anchor=south east]{$\varepsilon$};         \fill[fill=red!5,draw=orange!50!red,domain=0:1,variable=\x,samples=20] plot (\x, \x^2-\x^3) -- cycle;         \draw[->] (1.1, 0.3) node[anchor=west]{$e(h)=\left|\ape{bh^{p+1}}-\apz{ah^{p+2}}\right|$} -- (0.7, 0.1);         \draw[color=blue,domain=0:1,variable=\x,samples=20] plot (\x, \x^3);         \draw (0.8, 0.7) node[anchor=south east]{$\ape{bh^{p+1}}$};         \draw[color=green!70!blue,domain=0:1,variable=\x,samples=20] plot (\x, \x^2);         \draw (0.75, 0.5) node[anchor=north west]{$\apz{ah^{p+2}}$};         \draw[dotted] (\Hk, 0) node[anchor=north]{$h_k$} --++ (0, \Hk^2) --+ (-\Hk, 0) node[anchor=east]{$\ape{\tilde\varepsilon}$};         \draw[dotted] (\Hk, \Hk^3) --+ (-\Hk, 0) node[anchor=south east]{$\apz{\hat\varepsilon}$};         \draw[dotted] (\Hs, 0) node[anchor=north]{$h_k^*$} --++ (0, \Hs^3) --+ (-\Hs, 0) node[anchor=east]{$\varepsilon_0^*$};     \end{tikzpicture}     </math>

Damit ergeben sich folgende Gleichungen(<math>\textcolor{red}{Unbekannte}</math>):

<math>     \newcommand\ape[1]{\textcolor{green}{#1}}     \newcommand\apz[1]{\textcolor{blue}{#1}}     \newcommand\unbk[1]{\textcolor{red}{#1}}     \begin{array}{rl}         \unbk{a}h_k^{p+2}&=\unbk{(}\apz{\tilde\varepsilon}\unbk{)}=\left|\apz{\hat x_{k+1}}-\unbk{x_{k+1}}\right|\\         \unbk{b}h_k^{p+1}&=\unbk{(}\ape{\hat\varepsilon}\unbk{)}=\left|\ape{\tilde x_{k+1}}-\unbk{x_{k+1}}\right|\\         \unbk{ch_k^*}^{p+2}&=\varepsilon_0^*=\left|\unbk{\overset{opt}{x_{k+1}}}-\unbk{x_{k+1}}\right|\\         e(h_k)=\left|\unbk{b}h_k^{p+1}-\unbk{a}h_k^{p+2}\right|&=\left|\unbk{(}\ape{\tilde\varepsilon}\unbk{)}-\unbk{(}\apz{\hat\varepsilon}\unbk{)}\right|=\left|\ape{\tilde x_{k+1}}-\apz{\hat x_{k+1}}\right|     \end{array}     </math>

Das liegt der Annahme zugrunde, dass alle Koeffizienten der Fehlergliedern höherer Ordnung <math>\alpha_{n+1}=0</math> sind. Nun sind immernoch sechs Unbekannte in drei Gleichungen. Nun wird angenommen, dass <math>a=b=c</math> gilt und erhält damit zwei Gleichungen mit zwei Unbekannten, welche gelöst werden kann:

<math>     \newcommand\ape[1]{\textcolor{green}{#1}}     \newcommand\apz[1]{\textcolor{blue}{#1}}     \newcommand\unbk[1]{\textcolor{red}{#1}}     \begin{array}{rl}\nonumber         \left|\ape{\tilde x_{k+1}}-\apz{\hat x_{k+1}}\right|&=\left|\unbk{a}h_k^{p+1}-\unbk{a}h_k^{p+2}\right|\\         \unbk{ah_k^*}^{p+2}&=\varepsilon_0^*\\\nonumber         \Leftrightarrow\left|\ape{\tilde x_{k+1}}-\apz{\hat x_{k+1}}\right|&=\left|\unbk{a}h_k^{p+1}-\unbk{a}h_k^{p+2}\right|\\         \unbk{a}&=\frac{\varepsilon_0^*}{\unbk{h_k^*}^{p+2}}\\\nonumber         \Leftrightarrow\left|\ape{\tilde x_{k+1}}-\apz{\hat x_{k+1}}\right|&=\frac{\varepsilon_0^*}{\unbk{h_k^*}^{p+2}}\cdot\left|h_k^{p+1}-h_k^{p+2}\right|\\         \Leftrightarrow\unbk{h_k^*}^{p+2}\cdot&=\frac{\varepsilon_0^*}{\left|\frac{\ape{\tilde x_{k+1}}-\apz{\hat x_{k+1}}=\varepsilon_k}{h_k^{p+1}-h_k^{p+2}}\right|}\\         \Leftrightarrow\unbk{h_k^*}&=\sqrt[p+2]{\frac{\varepsilon_0^*\left|h_k^{p+1}-h_k^{p+2}\right|}{\varepsilon_k}}     \end{array}     </math>

Nun werden wieder die Sicherheitsfaktoren eingesetzt:

<math>     \newcommand\unbk[1]{\textcolor{red}{#1}}     \begin{equation}\nonumber         \unbk{h_k^*}=\beta\cdot\sqrt[p+2]{\frac{\varepsilon_0^*\left|h_k^{p+1}-h_k^{p+2}\right|}{\varepsilon_k}},\quad|\varepsilon_k<\varepsilon_0^*     \end{equation}     </math>

Damit kann nun die neue Schrittweite besser geschätzt werden.


4: Programmierung

Der Algorithmus wurde nach Vorlage der Matlab-Lösung aufgebaut("...\MATLAB\R20XXx\toolbox\matlab\funfun\ode45.m"):


matlab
  1. properties (Access = public)
  2. RKMethod % The Runge-Kutta-method
  3. AbsTol = 1e3; % Absolute tolerance
  4. RelTol = 1e3; % Relative tolerance
  5. AlphaMin = 0.1; % Smallest possible factor for new stepsize [0.1, 0.5]
  6. AlphaMax = 2; % Biggest possible factor for new stepsize [1.5, 2]
  7. Beta = 0.8; % Factor of safety for new stepsize [0.8, 0.95]
  8. MaxStep = 0.3; % Maximal stepsize
  9. MinStep = 16*eps; % Minimal stepsize
  10. DispStats = false;
  11. end
  12. %...
  13. function [X, T] = Embedded(self, t_0, t_e, x_0, h_0)
  14. nSuccess = 0;
  15. nFailure = 0;
  16. nFuncEval = 1;
  17. t_k = t_0;
  18. T = t_k;
  19. x_k = x_0;
  20. X = x_k';
  21. p = self.RKMethod.order_of_error;
  22. threshold = 1000*self.AbsTol;
  23. B1 = self.RKMethod.BtchrTableau.b(:, 1);
  24. B2 = self.RKMethod.BtchrTableau.b(:, 2);
  25. B = B1 - B2;
  26. k0 = self.RKMethod.ODS.f(t_0, x_0)';
  27. if (h_0 == 0)
  28. h_0 = min(self.MaxStep, abs(t_e - t_0));
  29. rh = norm(k0' ./ max(abs(X),threshold),inf) / (self.Beta * nthroot(self.RelTol, p));
  30. if h_0 * rh > 1
  31. h_0 = 1 / rh;
  32. end
  33. h_0 = max(h_0, self.MinStep);
  34. end
  35. h_k = max(min(h_0, self.MaxStep), self.MinStep);
  36. l = size(x_k, 1);
  37. while (t_k < t_e - eps)
  38. accepted = false;
  39. h_k = min(self.MaxStep, max(self.MinStep, h_k));
  40. % Stretch the step if within 10% of tfinal-t.
  41. if 1.1*h_k >= abs(t_e - t_k)
  42. h_k = abs(t_e - t_k);
  43. end
  44. noFailBefore = true;
  45.  
  46. while (~accepted)
  47. k = self.RKMethod.CalculateStep(t_k, x_k(1:l, 1:end), h_k, k0);
  48. x_kp1 = x_k + h_k*k*B1;
  49. nFuncEval = nFuncEval + length(k) - 1;
  50.  
  51. % error estimation:
  52. eps_k = h_k * norm((k*B) ./ max(max(abs(x_k),abs(x_kp1)),threshold),inf);
  53. accepted = (eps_k <= self.RelTol);
  54. if (~accepted)
  55. h_k = max(self.MinStep, h_k * max(self.AlphaMin, self.Beta*nthroot(self.RelTol/eps_k, p)));
  56. if h_k < self.MinStep
  57. error('stepsize below limit');
  58. end
  59. h_kp1 = h_k;
  60. noFailBefore = false;
  61. nFailure = nFailure + 1;
  62. elseif (noFailBefore)
  63. temp = 1.25*nthroot(eps_k/self.RelTol, p);
  64. if temp > 0.2
  65. h_kp1 = h_k/temp;
  66. else
  67. h_kp1 = 5*h_k;
  68. end
  69. h_kp1 = 0.9*self.Beta*nthroot(self.RelTol/eps_k*(h_k^(p+1)-h_k^(p+2))*h_k, p+2);
  70. nSuccess = nSuccess + 1;
  71. else
  72. noFailBefore = true;
  73. h_kp1 = h_k;
  74. nSuccess = nSuccess + 1;
  75. end
  76. end
  77.  
  78. x_k = x_kp1;
  79. t_k = t_k + h_k;
  80. h_k = h_kp1;
  81. T = [T; t_k];
  82. X = [X; x_kp1(1:l, 1:end)'];
  83. k0 = k(:, end);
  84. end
  85. if (self.DispStats)
  86. strsetting1 = ['%-24s:\n'];
  87. strsetting2 = ['%-24s:% 6d\n'];
  88. fprintf(strsetting1, 'Embedded-Stats');
  89. fprintf(strsetting2, 'Function-evaluations', nFuncEval);
  90. fprintf(strsetting2, 'Calculated steps', nSuccess + nFailure);
  91. fprintf(strsetting2, 'Successful steps', nSuccess);
  92. fprintf(strsetting2, 'Failed steps', nFailure);
  93. end
  94. end
DiffSolver.m

Die Schrittweitenabschätzung geschieht in Zeile 69, wobei der Koeffizient noch etwas stärker geschwächt werden muss, um gleich viele Funktionsauswertungen zu errechnen.


5: Vergleich zum herkömmlichen Algorithmus

Gegeben ist die folgende ODE:

<math>     x'(t)=-\left(\sin(t^3)+3t^3\cos(t^3)\right)\cdot x,\qquad x(0)=1     </math>

Mit der exakten Lösung:

<math>     x(t)=e^{-t\sin(t^3)}     </math>

Die Approximation sowie die Lösungskurve werden im Intervall <math>I=[0,3]</math> betrachtet, worauf die Lösungskurve wie folgt aussieht:

<math>     \begin{tikzpicture}     \begin{axis}[         xmin=0, xmax=3,         ymin=0, ymax=18,         xtick={0,1,2,3},         ytick={0,1,3,6,9,12,15,18},         ymajorgrids=true,         xmajorgrids=true,         x post scale=2.5,         grid style=dashed,     ]     \addplot[color=red,domain=0:3,samples=1000] {e^(-x*sin(x^3*180/pi))};     \end{axis}     \end{tikzpicture}     </math>

Um nun einen Vergleich der beiden Algorithmen gewährleisten zu können müssen ihre Parameter so gleich wie möglich sein. Die Matlab-Lösungsroutine ode45() interpoliert beispielsweise noch weitere Punkte für jeden Schritt(Option "Refine"), was im erarbeiteten Algorithmus nicht implementiert wurde.
Der absolut zugelassene Fehler pro Schritt wird auf <math>\varepsilon_{abs}=1e^{-2}</math> und der relativ zugelassene Fehler pro Schritt auf <math>\varepsilon_{rel}=1e^{-3}</math> gesetzt. Die maximale Schrittweite wird wie bei ode45() auf einem Zehntel der Intervalllänge gesetzt: <math>h_{max}=0.3</math> und <math>\beta=0.8</math> sowie <math>\alpha=0.9</math>.
Damit ergibt sich folgendes Bild:

<math>     \begin{tikzpicture}[scale=0.438]     \begin{axis}[     width=15.5in,     height=9.533802in,     at={(2.6in,1.286771in)},     scale only axis,     unbounded coords=jump,     xmin=0,     xmax=3,     xlabel={$t$},     xmajorgrids,     ymin=0,     ymax=18,     ylabel={$x$},     ymajorgrids,     zmin=0,     zmax=0.5,     zlabel={$h$},     zmajorgrids,     view={0}{90},     title style={font=\bfseries},     title={\huge Approximation der ODE},     legend style={at={(0.417535,0.684476)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}     ]     \addplot [color=green,solid,domain=0:3,samples=1000] {e^(-x*sin(x^3*180/pi))};     \addlegendentry{\huge Lösungskurve};      \addplot[color=red,solid,mark=+,mark options={solid}]      table[row sep=crcr] {     0    1\\     0.3    0.991932924766409\\     0.6    0.87933185610863\\     0.9    0.549081765262921\\     1.2    0.305681366853469\\     1.5    1.42351902547222\\     1.59112023546321    3.45296258345408\\     1.68224047092642    5.4009710822195\\     1.77861581831144    2.98241187017865\\     1.87499116569646    0.571468444052614\\     1.95380909747185    0.164763291550242\\     2.03262702924724    0.175886415612087\\     2.12778414952011    1.55336341165315\\     2.1777521597779    5.52978846648203\\     2.22772017003569    9.23940844001603\\     2.27417869297541    5.14666113557036\\     2.32063721591513    1.17490340132937\\     2.35670450986101    0.309968210586736\\     2.39277180380688    0.114935934087523\\     2.43735073432531    0.10119541955808\\     2.49448823995335    0.634499873215585\\     2.53722079661798    4.43714306652752\\     2.56914833441594    11.5041032213757\\     2.60107587221391    11.8808972616692\\     2.62445604581586    6.27996464568087\\     2.64447035978651    2.52756471085025\\     2.66451843688945    0.841005166516692\\     2.68129364621518    0.331700570762945\\     2.69955083413251    0.139024799049777\\     2.7244809924526    0.0696176487709387\\     2.75751777716675    0.095672320599614\\     2.80686551262708    1.40249837912733\\     2.83398542489711    7.14752549880928\\     2.86110533716714    16.892422412315\\     2.88313958616546    14.0871947666446\\     2.90099531625651    6.71238309342761\\     2.91591238096774    2.63030414153191\\     2.92916005019101    0.997056246858232\\     2.94116129266428    0.406429387349556\\     2.9541823857941    0.167038674428572\\     2.97049570038561    0.0726493052605177\\     3    0.0567595448147536\\     };      \addlegendentry{\huge Verbesserter Algorithmus};      \addplot [color=blue,solid,mark=o,mark options={solid}]      table[row sep=crcr] {     0    1\\     0.3    0.991932924766409\\     0.6    0.87933185610863\\     0.9    0.549081765262921\\     1.2    0.305681366853469\\     1.5    1.42351902547222\\     1.5769312492176    3.04954984174219\\     1.65386249843519    5.10772340192303\\     1.72490187369762    4.86225832514755\\     1.79594124896004    2.34406286869554\\     1.87112568208866    0.614739465605984\\     1.94682346430978    0.177279664791787\\     2.04638575186337    0.215160432093097\\     2.1413834738008    2.28375419495105\\     2.18558234832801    6.42098303580685\\     2.22978122285522    9.23850598171357\\     2.28193548460479    4.2345023630199\\     2.31678392689203    1.35728195818149\\     2.35163236917926    0.370922000616778\\     2.3847010160304    0.136118883462149\\     2.46689030997514    0.209068362319873\\     2.52550956585548    2.69056376187815\\     2.55655325828498    8.66410150528848\\     2.58759695071449    13.3377693102215\\     2.62536958771215    6.0621338280467\\     2.64700636270726    2.21425908537444\\     2.66864313770236    0.666939704541487\\     2.6885937567247    0.228704687052871\\     2.7147440643487    0.0838338995667424\\     2.76065890886312    0.107108382808794\\     2.80610321600487    1.34866363536941\\     2.83385307407508    7.1793236261912\\     2.8616029321453    17.1690326348336\\     2.88194655789135    14.6971333178647\\     2.9022901836374    6.31398708542791\\     2.9195270099869    2.05863694659545\\     2.93471219996752    0.664574170184959\\     2.94957894850768    0.228058645717422\\     2.96880107122159    0.0786675019359497\\     3    0.0576746694542346\\     };      \addlegendentry{\huge Herkömmlicher Algorithmus};      \end{axis}     \end{tikzpicture}     </math>

<math>     \begin{tikzpicture}[scale=0.438]     \begin{axis}[     width=15.5in,     height=9.533802in,     at={(2.6in,1.286771in)},     scale only axis,     unbounded coords=jump,     xmin=0,     xmax=3,     tick align=outside,     xlabel={\huge $t$},     xmajorgrids,     ymin=0,     ymax=18,     ylabel={\huge $x$},     ymajorgrids,     zmin=0,     zmax=0.5,     zlabel={\huge $h$},     zmajorgrids,     view={15}{40},     title style={font=\bfseries},     title={\huge Schrittweite beim Lösen der ODE},     legend style={at={(0.536806,0.617691)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}     ]     \addplot [color=green,solid,domain=0:3,samples=1000] {e^(-x*sin(x^3*180/pi))};     \addlegendentry{\huge Lösungskurve};      \addplot3 [color=red,solid,mark=+,mark options={solid}]      table[row sep=crcr] {     0    1    0.3\\     0.3    0.991932924766409    0.3\\     0.6    0.87933185610863    0.3\\     0.9    0.549081765262921    0.3\\     1.2    0.305681366853469    0.3\\     1.5    1.42351902547222    0.0911202354632106\\     1.59112023546321    3.45296258345408    0.0911202354632108\\     1.68224047092642    5.4009710822195    0.0963753473850195\\     1.77861581831144    2.98241187017865    0.0963753473850195\\     1.87499116569646    0.571468444052614    0.0788179317753912\\     1.95380909747185    0.164763291550242    0.0788179317753912\\     2.03262702924724    0.175886415612087    0.0951571202728658\\     2.12778414952011    1.55336341165315    0.0499680102577882\\     2.1777521597779    5.52978846648203    0.0499680102577882\\     2.22772017003569    9.23940844001603    0.0464585229397243\\     2.27417869297541    5.14666113557036    0.0464585229397243\\     2.32063721591513    1.17490340132937    0.0360672939458753\\     2.35670450986101    0.309968210586736    0.0360672939458753\\     2.39277180380688    0.114935934087523    0.0445789305184259\\     2.43735073432531    0.10119541955808    0.0571375056280372\\     2.49448823995335    0.634499873215585    0.0427325566646304\\     2.53722079661798    4.43714306652752    0.0319275377979671\\     2.56914833441594    11.5041032213757    0.0319275377979671\\     2.60107587221391    11.8808972616692    0.0233801736019523\\     2.62445604581586    6.27996464568087    0.0200143139706479\\     2.64447035978651    2.52756471085025    0.0200480771029365\\     2.66451843688945    0.841005166516692    0.0167752093257278\\     2.68129364621518    0.331700570762945    0.0182571879173361\\     2.69955083413251    0.139024799049777    0.0249301583200854\\     2.7244809924526    0.0696176487709387    0.0330367847141564\\     2.75751777716675    0.095672320599614    0.0493477354603238\\     2.80686551262708    1.40249837912733    0.0271199122700314\\     2.83398542489711    7.14752549880928    0.0271199122700314\\     2.86110533716714    16.892422412315    0.0220342489983185\\     2.88313958616546    14.0871947666446    0.0178557300910551\\     2.90099531625651    6.71238309342761    0.0149170647112205\\     2.91591238096774    2.63030414153191    0.0132476692232752\\     2.92916005019101    0.997056246858232    0.0120012424732652\\     2.94116129266428    0.406429387349556    0.01302109312982\\     2.9541823857941    0.167038674428572    0.0163133145915095\\     2.97049570038561    0.0726493052605177    0.0295042996143948\\     3    0.0567595448147536    0\\     };      \addlegendentry{\huge Verbesserter Algorithmus};      \addplot3 [color=red,solid]      table[row sep=crcr] {     0    18    0\\     0.3    18    7.68014643215409e-07\\     0.6    18    1.46356185304342e-06\\     0.9    18    2.08700339310308e-06\\     1.2    18    2.83973062631393e-06\\     1.5    18    0.00878616992773185\\     1.59112023546321    18    0.0214851266855014\\     1.68224047092642    18    0.0338933565228299\\     1.77861581831144    18    0.020936589752214\\     1.87499116569646    18    0.00556739471477963\\     1.95380909747185    18    5.18680170079611e-05\\     2.03262702924724    18    0.000233443436315484\\     2.12778414952011    18    0.000638491108475359\\     2.1777521597779    18    0.00214514374834085\\     2.22772017003569    18    0.00213973165588577\\     2.27417869297541    18    0.00140976023462258\\     2.32063721591513    18    0.00163915546825777\\     2.35670450986101    18    0.00174246966404568\\     2.39277180380688    18    0.000439128327804059\\     2.43735073432531    18    0.000525371136222433\\     2.49448823995335    18    0.00416429994357193\\     2.53722079661798    18    0.020944018200967\\     2.56914833441594    18    0.0552723633990055\\     2.60107587221391    18    0.0608768018862111\\     2.62445604581586    18    0.0308066383614412\\     2.64447035978651    18    0.0129317604930326\\     2.66451843688945    18    0.00575525057378168\\     2.68129364621518    18    0.00241720723981259\\     2.69955083413251    18    0.00101406475135191\\     2.7244809924526    18    0.000485808373559646\\     2.75751777716675    18    0.000719688201776481\\     2.80686551262708    18    0.00752634098617788\\     2.83398542489711    18    0.042499091823835\\     2.86110533716714    18    0.097992285104219\\     2.88313958616546    18    0.0793316399369175\\     2.90099531625651    18    0.0397736293541842\\     2.91591238096774    18    0.0145018069067935\\     2.92916005019101    18    0.00473414262304805\\     2.94116129266428    18    0.00176314820309537\\     2.9541823857941    18    0.000685920951040253\\     2.97049570038561    18    0.00032260119277612\\     3    18    1.11430193948728e-05\\     };      \addlegendentry{\huge Fehler des verbesserten Algorithmus'};      \addplot3 [color=blue,solid,mark=o,mark options={solid}]      table[row sep=crcr] {     0    1    0.3\\     0.3    0.991932924766409    0.3\\     0.6    0.87933185610863    0.3\\     0.9    0.549081765262921    0.3\\     1.2    0.305681366853469    0.3\\     1.5    1.42351902547222    0.076931249217596\\     1.5769312492176    3.04954984174219    0.076931249217596\\     1.65386249843519    5.10772340192303    0.0710393752624232\\     1.72490187369762    4.86225832514755    0.0710393752624232\\     1.79594124896004    2.34406286869554    0.0751844331286231\\     1.87112568208866    0.614739465605984    0.0756977822211189\\     1.94682346430978    0.177279664791787    0.0995622875535851\\     2.04638575186337    0.215160432093097    0.0949977219374341\\     2.1413834738008    2.28375419495105    0.0441988745272108\\     2.18558234832801    6.42098303580685    0.0441988745272108\\     2.22978122285522    9.23850598171357    0.0521542617495649\\     2.28193548460479    4.2345023630199    0.034848442287239\\     2.31678392689203    1.35728195818149    0.034848442287239\\     2.35163236917926    0.370922000616778    0.0330686468511385\\     2.3847010160304    0.136118883462149    0.0821892939447406\\     2.46689030997514    0.209068362319873    0.0586192558803345\\     2.52550956585548    2.69056376187815    0.0310436924295039\\     2.55655325828498    8.66410150528848    0.0310436924295039\\     2.58759695071449    13.3377693102215    0.0377726369976656\\     2.62536958771215    6.0621338280467    0.0216367749951059\\     2.64700636270726    2.21425908537444    0.0216367749951059\\     2.66864313770236    0.666939704541487    0.0199506190223384\\     2.6885937567247    0.228704687052871    0.0261503076239973\\     2.7147440643487    0.0838338995667424    0.0459148445144164\\     2.76065890886312    0.107108382808794    0.0454443071417523\\     2.80610321600487    1.34866363536941    0.0277498580702167\\     2.83385307407508    7.1793236261912    0.0277498580702167\\     2.8616029321453    17.1690326348336    0.0203436257460474\\     2.88194655789135    14.6971333178647    0.0203436257460474\\     2.9022901836374    6.31398708542791    0.0172368263495053\\     2.9195270099869    2.05863694659545    0.0151851899806217\\     2.93471219996752    0.664574170184959    0.0148667485401548\\     2.94957894850768    0.228058645717422    0.0192221227139098\\     2.96880107122159    0.0786675019359497    0.0311989287784127\\     3    0.0576746694542346    0\\     };      \addlegendentry{\huge Herkömmlicher Algorithmus};      \addplot3 [color=blue,solid]      table[row sep=crcr] {     0    18    0\\     0.3    18    7.68014643215409e-07\\     0.6    18    1.46356185304342e-06\\     0.9    18    2.08700339310308e-06\\     1.2    18    2.83973062631393e-06\\     1.5    18    0.00878616992773185\\     1.5769312492176    18    0.0188963278929903\\     1.65386249843519    18    0.0317297461519352\\     1.72490187369762    18    0.0305116748382028\\     1.79594124896004    18    0.0135604811306149\\     1.87112568208866    18    0.00527878306961527\\     1.94682346430978    18    0.000180072115700536\\     2.04638575186337    18    0.00170564231453019\\     2.1413834738008    18    0.00597081723820825\\     2.18558234832801    18    0.0171668246299932\\     2.22978122285522    18    0.02555710042264\\     2.28193548460479    18    0.00699407763693749\\     2.31678392689203    18    0.00367610422135733\\     2.35163236917926    18    0.00252866154746584\\     2.3847010160304    18    0.000799487911716384\\     2.46689030997514    18    0.00327188355843921\\     2.52550956585548    18    0.0153377400744725\\     2.55655325828498    18    0.0496367764092316\\     2.58759695071449    18    0.0777815043020968\\     2.62536958771215    18    0.0289832078561361\\     2.64700636270726    18    0.0116443543315126\\     2.66864313770236    18    0.00565794011313936\\     2.6885937567247    18    0.00214201910171316\\     2.7147440643487    18    0.000620181411517234\\     2.76065890886312    18    0.00149674536931059\\     2.80610321600487    18    0.00782329857711317\\     2.83385307407508    18    0.0365254049835375\\     2.8616029321453    18    0.0900146980785195\\     2.88194655789135    18    0.07831152924226\\     2.9022901836374    18    0.0284281680131864\\     2.9195270099869    18    0.0125361638244166\\     2.93471219996752    18    0.00554848720549117\\     2.94957894850768    18    0.00215687331124195\\     2.96880107122159    18    0.000637566660123781\\     3    18    0.000926267658875879\\     };      \addlegendentry{\huge Fehler des herkömmlichen Algorithmus'};      \end{axis}     \end{tikzpicture}     </math>

<math>     \begin{tikzpicture}[scale=0.438]     \pgfplotsset{     width=15.5in,     height=9.533802in,     at={(2.6in,1.286771in)},     scale only axis,     unbounded coords=jump,     xmin=0,     xmax=3,     tick align=outside,     title style={font=\bfseries},     legend style={at={(0.536806,0.617691)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}     }      \begin{axis}[     ymin=0,     ymax=18,     axis y line*=right,     axis x line=none,     ]     \addplot [color=green,solid,domain=0:3,samples=1000] {e^(-x*sin(x^3*180/pi))};     \label{Lkurve}     \end{axis}      \begin{axis}[     ymin=0,     ymax=1,     axis y line*=left,     xmajorgrids,     ymajorgrids,     xlabel={\huge $t$},     ylabel={\huge $h$}     ]     \addlegendimage{/pgfplots/refstyle=Lkurve}\addlegendentry{\huge Lösungskurve}     \addplot [color=red,solid,mark=+,mark options={solid}]      table[row sep=crcr] {     0    0.3\\     0.3    0.3\\     0.6    0.3\\     0.9    0.3\\     1.2    0.3\\     1.5    0.0911202354632106\\     1.59112023546321    0.0911202354632108\\     1.68224047092642    0.0963753473850195\\     1.77861581831144    0.0963753473850195\\     1.87499116569646    0.0788179317753912\\     1.95380909747185    0.0788179317753912\\     2.03262702924724    0.0951571202728658\\     2.12778414952011    0.0499680102577882\\     2.1777521597779    0.0499680102577882\\     2.22772017003569    0.0464585229397243\\     2.27417869297541    0.0464585229397243\\     2.32063721591513    0.0360672939458753\\     2.35670450986101    0.0360672939458753\\     2.39277180380688    0.0445789305184259\\     2.43735073432531    0.0571375056280372\\     2.49448823995335    0.0427325566646304\\     2.53722079661798    0.0319275377979671\\     2.56914833441594    0.0319275377979671\\     2.60107587221391    0.0233801736019523\\     2.62445604581586    0.0200143139706479\\     2.64447035978651    0.0200480771029365\\     2.66451843688945    0.0167752093257278\\     2.68129364621518    0.0182571879173361\\     2.69955083413251    0.0249301583200854\\     2.7244809924526    0.0330367847141564\\     2.75751777716675    0.0493477354603238\\     2.80686551262708    0.0271199122700314\\     2.83398542489711    0.0271199122700314\\     2.86110533716714    0.0220342489983185\\     2.88313958616546    0.0178557300910551\\     2.90099531625651    0.0149170647112205\\     2.91591238096774    0.0132476692232752\\     2.92916005019101    0.0120012424732652\\     2.94116129266428    0.01302109312982\\     2.9541823857941    0.0163133145915095\\     2.97049570038561    0.0295042996143948\\     3    0\\     };      \addlegendentry{\huge Verbesserter Algorithmus};      \addplot [color=red,solid]      table[row sep=crcr] {     0    0 0\\     0.3    7.68014643215409e-07\\     0.6    1.46356185304342e-06\\     0.9    2.08700339310308e-06\\     1.2    2.83973062631393e-06\\     1.5    0.00878616992773185\\     1.59112023546321    0.0214851266855014\\     1.68224047092642    0.0338933565228299\\     1.77861581831144    0.020936589752214\\     1.87499116569646    0.00556739471477963\\     1.95380909747185    5.18680170079611e-05\\     2.03262702924724    0.000233443436315484\\     2.12778414952011    0.000638491108475359\\     2.1777521597779    0.00214514374834085\\     2.22772017003569    0.00213973165588577\\     2.27417869297541    0.00140976023462258\\     2.32063721591513    0.00163915546825777\\     2.35670450986101    0.00174246966404568\\     2.39277180380688    0.000439128327804059\\     2.43735073432531    0.000525371136222433\\     2.49448823995335    0.00416429994357193\\     2.53722079661798    0.020944018200967\\     2.56914833441594    0.0552723633990055\\     2.60107587221391    0.0608768018862111\\     2.62445604581586    0.0308066383614412\\     2.64447035978651    0.0129317604930326\\     2.66451843688945    0.00575525057378168\\     2.68129364621518    0.00241720723981259\\     2.69955083413251    0.00101406475135191\\     2.7244809924526    0.000485808373559646\\     2.75751777716675    0.000719688201776481\\     2.80686551262708    0.00752634098617788\\     2.83398542489711    0.042499091823835\\     2.86110533716714    0.097992285104219\\     2.88313958616546    0.0793316399369175\\     2.90099531625651    0.0397736293541842\\     2.91591238096774    0.0145018069067935\\     2.92916005019101    0.00473414262304805\\     2.94116129266428    0.00176314820309537\\     2.9541823857941    0.000685920951040253\\     2.97049570038561    0.00032260119277612\\     3    1.11430193948728e-05\\     };      \addlegendentry{\huge Fehler des verbesserten Algorithmus'};      \addplot [color=blue,solid,mark=o,mark options={solid}]      table[row sep=crcr] {     0    0.3\\     0.3    0.3\\     0.6    0.3\\     0.9    0.3\\     1.2    0.3\\     1.5    0.076931249217596\\     1.5769312492176    0.076931249217596\\     1.65386249843519    0.0710393752624232\\     1.72490187369762    0.0710393752624232\\     1.79594124896004    0.0751844331286231\\     1.87112568208866    0.0756977822211189\\     1.94682346430978    0.0995622875535851\\     2.04638575186337    0.0949977219374341\\     2.1413834738008    0.0441988745272108\\     2.18558234832801    0.0441988745272108\\     2.22978122285522    0.0521542617495649\\     2.28193548460479    0.034848442287239\\     2.31678392689203    0.034848442287239\\     2.35163236917926    0.0330686468511385\\     2.3847010160304    0.0821892939447406\\     2.46689030997514    0.0586192558803345\\     2.52550956585548    0.0310436924295039\\     2.55655325828498    0.0310436924295039\\     2.58759695071449    0.0377726369976656\\     2.62536958771215    0.0216367749951059\\     2.64700636270726    0.0216367749951059\\     2.66864313770236    0.0199506190223384\\     2.6885937567247    0.0261503076239973\\     2.7147440643487    0.0459148445144164\\     2.76065890886312    0.0454443071417523\\     2.80610321600487    0.0277498580702167\\     2.83385307407508    0.0277498580702167\\     2.8616029321453    0.0203436257460474\\     2.88194655789135    0.0203436257460474\\     2.9022901836374    0.0172368263495053\\     2.9195270099869    0.0151851899806217\\     2.93471219996752    0.0148667485401548\\     2.94957894850768    0.0192221227139098\\     2.96880107122159    0.0311989287784127\\     3    0\\     };      \addlegendentry{\huge Herkömmlicher Algorithmus};      \addplot [color=blue,solid]      table[row sep=crcr] {     0    0\\     0.3    7.68014643215409e-07\\     0.6    1.46356185304342e-06\\     0.9    2.08700339310308e-06\\     1.2    2.83973062631393e-06\\     1.5    0.00878616992773185\\     1.5769312492176    0.0188963278929903\\     1.65386249843519    0.0317297461519352\\     1.72490187369762    0.0305116748382028\\     1.79594124896004    0.0135604811306149\\     1.87112568208866    0.00527878306961527\\     1.94682346430978    0.000180072115700536\\     2.04638575186337    0.00170564231453019\\     2.1413834738008    0.00597081723820825\\     2.18558234832801    0.0171668246299932\\     2.22978122285522    0.02555710042264\\     2.28193548460479    0.00699407763693749\\     2.31678392689203    0.00367610422135733\\     2.35163236917926    0.00252866154746584\\     2.3847010160304    0.000799487911716384\\     2.46689030997514    0.00327188355843921\\     2.52550956585548    0.0153377400744725\\     2.55655325828498    0.0496367764092316\\     2.58759695071449    0.0777815043020968\\     2.62536958771215    0.0289832078561361\\     2.64700636270726    0.0116443543315126\\     2.66864313770236    0.00565794011313936\\     2.6885937567247    0.00214201910171316\\     2.7147440643487    0.000620181411517234\\     2.76065890886312    0.00149674536931059\\     2.80610321600487    0.00782329857711317\\     2.83385307407508    0.0365254049835375\\     2.8616029321453    0.0900146980785195\\     2.88194655789135    0.07831152924226\\     2.9022901836374    0.0284281680131864\\     2.9195270099869    0.0125361638244166\\     2.93471219996752    0.00554848720549117\\     2.94957894850768    0.00215687331124195\\     2.96880107122159    0.000637566660123781\\     3    0.000926267658875879\\     };      \addlegendentry{\huge Fehler des herkömmlichen Algorithmus'};      \end{axis}     \end{tikzpicture}     </math>

Es ergeben sich folgende Daten:

Algorithmus Berechnete Schritte Funktionsauswertungen Erfolgreiche Schritte Fehlerhafte Schritte
Herkömmlich 50 301 39 11
Verbessert 50 301 41 9

Der verbesserte Algorithmus erbringt demnach bessere Resultate als der herkömmliche Algorithmus, wobei der Fehler etwas kleiner wird.
Die Resultate variieren, je nach dem wie stark der Koeffizient für die nächste Schrittweite geändert wird. Werte wie <math>0.8</math> sind klein genug, um gleich viele Funktionsauswertungen zu bekommen und nurnoch 4 fehlerhafte Iterationen zu haben zum Preis eines leicht erhöhten Fehlers, wohingegen z.B. <math>0.7</math> einige mehr Auswertungen benötigt und der Fehler trotzdem grösser wird.
Als beide Algorithmen die Van-der-Pol-ODE lösen sollten(mit gleichen Parametern) explodierte der Fehler der herkömmlichen ode45-Methode, während die Verbesserte zuverlässig eine gute Approximation ablieferte.


6: Fazit



Die selbst-entwickelte Schrittweitenabschätzung funktioniert besser als die herkömmliche, wenn auch nur gering. Trotzdem könnte eine Berechnungszeiteinsparung gerade bei grossen Kalkulationen auch mit ein paar wenigen Prozent erstrebenswert sein.
<math>\square</math>



Link auf diesen Artikel Link auf diesen Artikel  Druckbare Version Druckerfreundliche Version  Einen Freund auf diesen Artikel aufmerksam machen Weitersagen Kommentare zeigen Kommentare  
pdfFür diesen Artikel gibt es keine pdf-Datei


Arbeitsgruppe Alexandria Dieser Artikel ist im Verzeichnis der Arbeitsgruppe Alexandria eingetragen:
: Mathematik :: automatisch eingefügt und unbearbeitet :
Verbesserung von Eingebetteten Runge-Kutta-Verfahren [von Higlav]  
div.article { position: relative; font-size: 10pt; text-align: left; } div.article h2 { padding-left: 10px; font-size: 16pt; color: #2e2e30; background-color: #dcdcdc; } div.article h3 { padding-left: 3
[Die Arbeitsgruppe Alexandria katalogisiert die Artikel auf dem Matheplaneten]

 
Verwandte Links
 
Besucherzähler 950
 
Aufrufstatistik des Artikels
Insgesamt 2 externe Besuche in 2017 [Anzeigen]
DomainAnzahlProz
http://google.de150%50 %
http://google.at150%50 %

[Seitenanfang]

" Mathematik: Verbesserung von Eingebetteten Runge-Kutta-Verfahren" | 6 Kommentare
 
Für den Inhalt der Kommentare sind die Verfasser verantwortlich.

Re: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
von cis am Mo. 04. Juli 2016 20:59:15


Tolle Graphiken  wink

 [Bearbeiten]

Re: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
von Slash am Di. 05. Juli 2016 00:46:38


Hi, zwei kleine Bugs: 1) Der Vorwort-Balken geht über das Kugelsymbol - sieht nicht schön aus. 2) Der Navi-Link zu Punkt 5 funktioniert nicht.

Sonst macht der Artikel aber optisch einen sehr guten Eindruck. smile

Gruß, Slash

EDIT: Bugs sind jetzt behoben. Danke.

 [Bearbeiten]

Re: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
von Delastelle am Fr. 08. Juli 2016 19:49:30


Hallo Higlav!

Ein Beispiel mit einer Differentialgleichung ist vielleicht etwas wenig
um ein Verfahren zu beurteilen.

Viele Grüße
Ronald

 [Bearbeiten]

Re: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
von Higlav am Fr. 15. Juli 2016 13:20:12


Hallo Delastelle,

Du hast schon Recht, dass eine einzige Differentialgleichung als Messlatte etwas ungeschickt ist. Jedoch habe ich das Ganze auch an verschiedenen Differentialgleichungen getestet:
VanDerPol:
<math>\vec r=\begin{pmatrix}r_1=r_2'\\r_2=8r_2(1-r_1^2)-r_1\end{pmatrix},\qquad\vec r^{(0)}=\begin{pmatrix}2\\0\end{pmatrix}</math>
Satellit:
<math>\vec r=\begin{pmatrix}u_1=u_2'\\u_2=a_1\\v_1=v_2'\\v_2=a_2\\w_1=w_2'\\w_2=a_3\end{pmatrix},\qquad
\vec x=\begin{pmatrix}u_1\\v_1\\w_1\end{pmatrix},\
\vec a=\begin{cases}
-\vec x\cdot G\frac M{|\vec x|^3}&,|\vec x|\geq r_E\\
\vec x\cdot G\frac M{r_E^4}&
\end{cases},\
\vec r^{(0)}=\begin{pmatrix}10e6\\2000\\-3e6\\-1000\\15e6\\2000\end{pmatrix}</math>
Rigid Body Problem:
<math>\vec r=\begin{pmatrix}r_1=r_2r_3\\r_2=-r_1r_3\\r_3=-0.51r_1r_2\end{pmatrix},\qquad\vec r^{(0)}=\begin{pmatrix}0\\1\\1\end{pmatrix}</math>
Und die modifizierte Methode erzielte immer bessere Resultate: Weniger Schritte und/oder kleinerer Fehler.
Klar ist auch mit diesen DGL nicht die ganze Palette abgedeckt, aber mir fiel da nicht mehr viel ein. Oder kennst du noch eine DGL, an der man die beiden Algorithmen "messen" könnte?

Mit freundlichen Grüssen,

Higlav


PS: @cis: Danke. Ich mag es, solche Grafiken zu erstellen. Ich visualisiere gerne. (Auch wenn die Plots beim Vergleich leider etwas klein und somit schwer lesbar geworden sind)

 [Bearbeiten]

Re: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
von Delastelle am So. 17. Juli 2016 21:13:05


Hallo Higlav!

Mir würde noch einfallen:
Räuber-Beute-Differentialgleichungen.

Neben der Schrittzahl ist auch die Rechenzeit wichtig.
In Matlab kann man mittels
TIC
Programm
TOC
die Rechenzeit ermitteln.

Auch ist der Vergleich mit anderen Lösern interessant.
Manche Probleme löst man besser mit einem ODE45. Es gibt aber auch
ODE23 oder ODE78 (nicht MATLAB).

Viele Grüße
Ronald

 [Bearbeiten]

Re: Verbesserung von Eingebetteten Runge-Kutta-Verfahren
von Higlav am Mo. 18. Juli 2016 11:22:56


Hallo Ronald,

Mit Räuber-Beute-Differentialgleichung meinst du höchstwahrscheinlich die von Lotka und Volterra, nicht wahr?
<math>
\vec r=\begin{pmatrix}r_1=r_1(\varepsilon_1-\gamma_1r_2)\\r_2=-r_2(\varepsilon_2-\gamma_2N_1)\end{pmatrix},\qquad\vec r^{(0)}=\begin{pmatrix}300\\300\end{pmatrix}
</math>

Die Matlab-Befehle TIC und TOC kannte ich noch nicht. Ich benutze immer die Run-and-Time-Funktion, um zu sehen welche Teile meines Programmes wie lange brauchen.
Allerdings verstehe ich deine Argumentation nicht: Es sind doch die Anweisungen an die CPU, die für die Rechenzeit massgebend sind. Die Zeit zu messen bei so ähnlichen Programmen finde ich unangemessen, da die Kausalkette bis zur entgültigen Ausführzeitspanne heutzutage mit Mehrkern-Prozessoren und Multi-Threading(bzw. allgemeiner: Parallelisierung) so eine Messung relativ unbrauchbar macht. Ich mag mich nicht mehr ganz erinnern, aber irgendwo hatten wir mal zwei Programme per Zeitmessung verglichen und die Fehler zur zweiten Messung waren enorm.
Was sich aber sagen lässt ist, dass meine Routine fünf Anweisungen mehr kompiliert: Zweimal die Potenzierung der Schrittweite, die Subtraktion der beiden Resultate, die Vorzeichen-Löschung und die Multiplikation mit dem Fehler-Quotient. Das sind ein paar Full-Subtractors, Shifts, Signs und derlei. Es hat nicht mal eine Wurzel.
Ausserdem ist der wesentliche Teil ja die Differentialgleichung selbst. Sie ist die grosse Unbekannte, wenn es zur Laufzeitberechnung kommt. Je weniger Funktionsauswertungen(ohne grossen Fehlerzuwachs) desto besser. Deshalb ist das klassische Runge-Kutta-Verfahren ja auch so beliebt, obwohl es 4 Stufen(also 4 Funktionsauswertungen pro Schritt) beinhaltet: Für eine höhere Genauigkeit erfordert es bereits ein 6-stufiges Verfahren. Das Euler-Verfahren ist zu ungenau und nurnoch zu Schulungszwecken zu verwenden. Zweistufige Verfahren wie das von Heun sind zwar stabiler und energieerhaltend, aber immernoch mit einem relativ grossen Fehler verbunden. Verfahren der 3-ten Ordnung sind zwar besser, aber nicht energieerhaltend. Nur die 4-Stufigen vereinigen Energieerhaltung und Präzision so gekonnt.
Meine Schrittweitenschätzung ist von der Stufe ja unabhängig, solange das Verfahren eingebettet ist. Es stimmt allerdings, dass ich es nur mit ode45 getestet habe - ich habe mir keine Mühe gemacht, noch das Butcher-Tableau für ode23 oder ode78 abzutippeln. Das deshalb, weil mein Auftrag, in dessem Rahmen ich mich mit dem Thema überhaupt befasst habe, sich nur auf den Vergleich zur Matlab-ode45-Routine bezog.
Ich werde allerdings mal diese Butcher-tableaus mal heraussuchen und die Verfahren so miteinander vergleichen - auch wenn ich einen vergleichbaren Effekt erwarte.

Mit freundlichen Grüssen,

Raphael

 [Bearbeiten]

 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2017 by Matroids Matheplanet
This web site was made with PHP-Nuke, a web portal system written in PHP. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]