v1.2 – Explizite Einschrittverfahren

 

1.2.1 Eulersches Polygonzugverfahren

Das einfachste explizite Einschrittverfahren ist das Eulersche Polygonzugverfahren. Man approximiert damit die Differentialgleichung

\dot y\left( t \right) = f\left( {t,y\left( t \right)} \right)

im Punkt {t_k} durch die Differenzengleichung

\frac{{{u_{k+1}}-{u_k}}}{{{\tau _k}}} = f\left( {{t_k},{u_k}} \right)

Äquivalent umgeformt bedeutet das:

{u_0} = {y_0}

{u_{k+1}} = {u_k}+{\tau _k}f\left( {{t_k},{u_k}} \right),\quad k = 0,1,2, \ldots

Beispiel:

\dot y\left( t \right) = 2y\left( t \right)+{t^2},\quad y\left( 0 \right) = 1,\quad \tau = 0,1

{u_0} = 1,\quad {u_{k+1}} = {u_k}+\tau \left( {2{u_k}+t_k^2} \right)

{u_1} = 1+0,1\left( {2 \cdot 1+{0^2}} \right) = 1,2

{u_2} = 1,2+0,1\left( {2 \cdot 1,2+{{0,1}^2}} \right) = 1,441

{u_3} = 1,441+0,1\left( {2 \cdot 1,441+{{0,2}^2}} \right) = 1,7332

{u_4} = 2,084\quad {u_5} = 2,523\quad {u_6} = 3,053\quad {u_7} = 3,699

1.2.2 Allgemeine Definition explizite Einschrittverfahren

Explizite Einschrittverfahren sind Verfahren der Form

{u_0} = {y_0}

{u_{k+1}} = {u_k}+{\tau _k}\Phi \left( {{t_k},{u_k},{\tau _k},f} \right),\quad k = 0,1,2, \ldots

zur Lösung von Anfangswertaufgaben. Die Funktion \Phi heißt Inkrement- bzw. Zuwachsfuntion.

Die Inkrementfunktion des expliziten Euler-Verfahrens lautet \Phi = f\left( {{t_k},{u_k}} \right).

1.2.3 Beispiel: Feder-Masse-Schwinger

num-102-feder-masse-dampfer-schwinger

Die Bewegung der Pendelmasse kann durch eine Differentialgleichung zweiter Ordnung beschrieben werden:

m\ddot x+c\dot x+kx = 0

Dabei ist m die Masse der Kugel, \ddot x die Beschleunigung, c die Dämpfungskonstante, \dot x die Geschwindigkeit, k die Federkonstante und x die Auslenkung. Die DGL entspricht dem zweiten Newtonschen Gesetz.

Die Differentialgleichung zweiter Ordnung kann mit Hilfe der neuen Variablen v = \dot x in ein System erster Ordnung umgeformt werden:

\begin{array}{*{20}{c}}{\dot x = v,} & {x\left( 0 \right) = {x_0}} \\{\dot v = -\frac{c}{m}v-\frac{k}{m}x,} & {v\left( 0 \right) = {v_0}} \\   \end{array}

In Matrixschreibweise ergibt sich:

\dot y = \frac{1}{m}\left( {\begin{array}{*{20}{c}}  0 & m \\{-k} & {-c} \\   \end{array} } \right)y,\quad \quad y = \left( {\begin{array}{*{20}{c}}  x \\  v \\   \end{array} } \right),\quad {y_0} = \left( {\begin{array}{*{20}{c}}{{x_0}} \\{{v_0}} \\   \end{array} } \right)

Für das explizite Eulerverfahren ergibt sich:

{u_0} = \left( {\begin{array}{*{20}{c}}{{x_0}} \\{{v_0}} \\   \end{array} } \right)

{u_{k+1}} = {u_k}+{\tau _k}f\left( {t,y} \right) = {u_k}+\frac{\tau }{m}\left( {\begin{array}{*{20}{c}}  0 & m \\{-k} & {-c} \\   \end{array} } \right){u_k} = \left( {\begin{array}{*{20}{c}}  1 & \tau \\{-\frac{{\tau k}}{m}} & {1-\frac{{\tau c}}{m}} \\   \end{array} } \right){u_k}

Dieses einfache Verfahren liefert erstaunlich gute Ergebnisse. Für kleinere Schrittweite \tau erhöht sich die Genauigkeit.

1.2.4 Lokaler Diskretisierungsfehler, Konsistenzfehler

Man bezeichnet die Funktion

{\ell _\tau }\left( t \right) = \frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau }-\Phi \left( {t,y\left( t \right),\tau ,f} \right)

als lokalen Diskretisierungsfehler oder Konsistenzfehler. \frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau } wäre genau der Wert der Ableitung, den man benutzen müsste, um bei gegebener Schrittweite von einem exakten Wert zum nächsten zu kommen. \Phi ist der Wert der Ableitung, den man stattdessen benutzt. Der lokale Diskretisierungsfehler ist also der Fehler, der bei genau einem Schritt des Einschrittverfahrens mit exaktem Startwert entsteht.

Gilt \left| {{\ell _\tau }\left( t \right)} \right| \leq M{\tau ^p},\quad p \geq 1 für \tau \leq {\tau _0} und \forall t \in \left[ {{t_0},T} \right], so heißt das entsprechende Diskretisierungsverfahren konsistent mit Ordnung p. Statt \left| {{\ell _\tau }\left( t \right)} \right| \leq M{\tau ^p} ist auch die Schreibweise {\ell _\tau }\left( t \right) = O\left( {{\tau ^p}} \right) gebräuchlich.

Konsistenz heißt insbesondere, dass

0 = \mathop {\lim }\limits_{\tau \to 0} {\ell _\tau }\left( t \right) = \underbrace {\dot y\left( t \right)}_{ = f\left( {t,y} \right)}-\mathop {\lim }\limits_{\tau \to 0} \Phi \left( {t,y\left( t \right),\tau ,f} \right)

Es wurde hier \mathop {\lim }\limits_{\tau \to 0} \frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau } = \dot y\left( t \right) gesetzt, da dieser Grenzwert ja genau die Definition des Differentialquotienten ist. Aus obiger Gleichung folgt:

\mathop {\lim }\limits_{\tau \to 0} \Phi \left( {t,y\left( t \right),\tau ,f} \right) = f\left( {y,t} \right)

1.2.5 Beispiel: Eulerverfahren

{\ell _\tau }\left( t \right) = \frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau }-\Phi \left( {t,y\left( t \right),\tau ,f} \right)

\quad = \frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau }-f\left( {t,y\left( t \right)} \right)

\quad = \frac{1}{\tau }\left[ {y\left( {t+\tau } \right)-y\left( t \right)} \right]-\dot y\left( t \right)

\quad = \frac{1}{\tau }\left[ {y\left( t \right)+\tau \dot y\left( t \right)+\frac{{{\tau ^2}}}{2}\ddot y\left( \xi \right)-y\left( t \right)} \right]-\dot y\left( t \right),\quad \quad \quad \xi \in \left( {t,t+\tau } \right)

\quad = \frac{{\ddot y\left( \xi \right)}}{2}\tau \quad \leq \quad M{\tau ^1}

Das Eulerverfahren ist also konsistent mit Ordnung 1, falls \ddot y im Intervall \left[ {{t_0},T} \right] beschränkt ist.

1.2.6 Explizite, r-stufige Runge-Kutta-Verfahren, Butcher-Schema

Gebräuchliche Einschrittverfahren höherer Ordnung, die expliziten r-stufigen Runge-Kutta-Verfahren, erhält man mit folgendem Ansatz:

\Phi \left( {t,u,\tau ,f} \right) = \sum\limits_{i = 1}^r {{\beta _i}{f_i}}

Dabei wurden die folgenden Abkürzungen benutzt:

{f_1} = f\left( {t,u} \right)

{f_i} = f\left( {t+{\gamma _i}\tau ,\:\:u+\tau \sum\limits_{s = 1}^{i-1} {{\beta _{i,s}}{f_s}} } \right),\quad \quad i = 2, \ldots ,r

{\beta _i},{\gamma _i},{\beta _{i,s}} sind reelle Parameter. Es ist üblich, diese Parameter in Form einer Tabelle, dem Butcher-Schema, anzugeben:

\begin{array}{*{20}{c}}{{\gamma _2}} &\vline & {{\beta _{2,1}}} & {} & {} & {} & {} \\{{\gamma _3}} &\vline & {{\beta _{3,1}}} & {{\beta _{3,2}}} & {} & {} & {} \\  \vdots &\vline & \vdots & \vdots & \ddots & {} & {} \\{{\gamma _r}} &\vline & {{\beta _{r,1}}} & {{\beta _{r,2}}} & \ldots & {{\beta _{r,r-1}}} & {} \\ \hline{} &\vline & {{\beta _1}} & {{\beta _2}} & \ldots & {{\beta _{r-1}}} & {{\beta _r}} \\   \end{array}

Für r = 1 und {\beta _1} = 1 ergibt sich das Eulerverfahren:

\begin{array}{*{20}{c}}{} &\vline & {} \\ \hline{} &\vline & 1 \\   \end{array}

1.2.7 Konstruktion eines Verfahrens mit Konsistenzordnung 2

Die Taylorentwicklung für y\left( {t+\tau } \right) bis zur dritten Ordnung lautet:

y\left( {t+\tau } \right) = y\left( t \right)+\tau \dot y\left( t \right)+\frac{{{\tau ^2}}}{2}\ddot y\left( t \right)+O\left( {{\tau ^3}} \right)

Für die zweite Ableitung muss das totale Differential benutzt werden:

\ddot y\left( t \right) = \frac{d}{{dt}}f\left( {t,y\left( t \right)} \right) = \frac{{\partial f}}{{\partial t}}+\frac{{\partial f}}{{\partial y}}\frac{{dy}}{{dt}} = \frac{{\partial f}}{{\partial t}}+\frac{{\partial f}}{{\partial y}}f

Damit erhalten wir für den exakten Differenzenquotienten:

\frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau } = \dot y+\frac{\tau }{2}\ddot y+O\left( {{\tau ^2}} \right) = f+\frac{\tau }{2}\left[ {\frac{{\partial f}}{{\partial t}}+\frac{{\partial f}}{{\partial y}}f} \right]+O\left( {{\tau ^2}} \right)

Für die Inkrementfunktion ergibt sich ebenfalls mit der Taylorentwicklung:

\Phi \left( {t,y\left( t \right),\tau ,f} \right) = {\beta _1}f\left( {t,y\left( t \right)} \right)+{\beta _2}f\left( {t+{\gamma _2}\tau ,\:\:y\left( t \right)+{\beta _{21}}\tau f\left( {t,y\left( t \right)} \right)} \right)

\quad \quad = {\beta _1}f+{\beta _2}\left[ {f+{\gamma _2}\tau \frac{{\partial f}}{{\partial t}}+{\beta _{21}}\tau f\frac{{\partial f}}{{\partial y}}} \right]+O\left( {{\tau ^2}} \right)

Daraus folgt für den Konsistenzfehler:

{\ell _\tau }\left( t \right) = \frac{{y\left( {t+\tau } \right)-y\left( t \right)}}{\tau }-\Phi \left( {t,y\left( t \right),\tau ,f} \right)

\quad = f+\frac{\tau }{2}\left[ {\frac{{\partial f}}{{\partial t}}+\frac{{\partial f}}{{\partial y}}f} \right]+O\left( {{\tau ^2}} \right)-{\beta _1}f-{\beta _2}\left[ {f+{\gamma _2}\tau \frac{{\partial f}}{{\partial t}}+{\beta _{21}}\tau f\frac{{\partial f}}{{\partial y}}} \right]+O\left( {{\tau ^2}} \right)

\quad = \left( {1-{\beta _1}-{\beta _2}} \right)f+\left( {\frac{\tau }{2}-{\beta _2}{\gamma _2}\tau } \right)\frac{{\partial f}}{{\partial t}}+\left( {\frac{\tau }{2}-{\beta _2}\tau {\beta _{21}}} \right)f\frac{{\partial f}}{{\partial y}}+O\left( {{\tau ^2}} \right)

Wir erreichen also Konsistenzordnung 2, wenn die Klammerausdrücke jeweils zu 0 werden:

{\beta _1}+{\beta _2} = 1

2{\beta _2}{\gamma _2} = 1

2{\beta _2}{\beta _{21}} = 1

Für beliebiges {\beta _2} folgt daraus:

{\beta _1} = 1-{\beta _2},\quad {\gamma _2} = {\beta _{21}} = \frac{1}{{2{\beta _2}}}

Für {\beta _2} = \frac{1}{2} ergibt sich daraus das Verfahren nach Heun (1900):

\begin{array}{*{20}{c}}  1 &\vline & 1 & {} \\ \hline{} &\vline & {\frac{1}{2}} & {\frac{1}{2}} \\   \end{array}

Ausgeschrieben bedeutet das:

{u_{k+1}} = {u_k}+\tau \left( {\frac{1}{2}{f_{1,k}}+\frac{1}{2}{f_{2,k}}} \right)

{f_{1,k}} = f\left( {{t_k},{u_k}} \right)

{f_{2,k}} = f\left( {{t_k}+\tau ,\:\:{u_k}+\tau {f_{1,k}}} \right)

Für {\beta _2} = 1 ergibt sich das Halbschrittverfahren nach Collatz (1969):

\begin{array}{*{20}{c}}{\frac{1}{2}} &\vline & {\frac{1}{2}} & {} \\ \hline{} &\vline & 0 & 1 \\   \end{array}

Ausgeschrieben bedeutet das:

{u_{k+1}} = {u_k}+\tau {f_{2,k}}

{f_{1,k}} = f\left( {{t_k},{u_k}} \right)

{f_{2,k}} = f\left( {{t_k}+\frac{\tau }{2},\:\:{u_k}+\frac{\tau }{2}{f_{1,k}}} \right)

1.2.8 Beispiele für 2-, 3- und 4-stufige Verfahren

3-stufiges Verfahren von Heun (Konsistenzordnung 3):

\begin{array}{*{20}{c}}{\frac{1}{3}} &\vline & {\frac{1}{3}} & {} & {} \\{\frac{2}{3}} &\vline & 0 & {\frac{2}{3}} & {} \\ \hline{} &\vline & {\frac{1}{4}} & 0 & {\frac{3}{4}} \\   \end{array}

Klassisches Runge-Kutta-Verfahren (Konsistenzordnung 4):

\begin{array}{*{20}{c}}{\frac{1}{2}} &\vline & {\frac{1}{2}} & {} & {} & {} \\{\frac{1}{2}} &\vline & 0 & {\frac{1}{2}} & {} & {} \\  1 &\vline & 0 & 0 & 1 & {} \\ \hline{} &\vline & {\frac{1}{6}} & {\frac{1}{3}} & {\frac{1}{3}} & {\frac{1}{6}} \\   \end{array}

1.2.9 Zusammenhang zwischen Stufenzahl und erreichbarer Ordnung

Resultat von Butcher:

p\quad \leq \quad \left\{ {\begin{array}{*{20}{c}}  r & {r = 1,2,3,4} \\{r-1} & {r = 5,6,7} \\{r-2} & {r = 8,9} \\   \end{array} } \right.

1.2.10 Globaler Diskretisierungsfehler, Konvergenz

In der Regel werden in jedem Zeitschritt {\tau _n} lokale Fehler erzeugt, die sich in komplizierter Weise übertragen und den globalen Diskretisierungsfehler

{e_\tau }\left( {{t_i}} \right) = y\left( {{t_i}} \right)-{u_\tau }\left( {{t_i}} \right)

bilden. Wie vorn ist die Gitterfunktion {u_\tau }:{\omega _\tau } \to {\mathbb{R}^n} durch {u_\tau }\left( {{t_i}} \right) = {u_i} definiert. Gilt für \tau \leq {\tau _0} und für alle {t_i} \in {\omega _\tau }

\left| {{e_\tau }\left( {{t_i}} \right)} \right| = O\left( {{\tau ^p}} \right)

dann heißt das entsprechende Diskretisierungsverfahren konvergent mit der Ordnung p.

1.2.11 Bedingung für Konvergenz

Ist y \in {C^{p+1}}\left[ {{t_0},T} \right] und besitzt das Einschrittverfahren die Konsistenzordnung p, dann ist das Verfahren auch konvergent mit der Ordnung p. Es gilt:

\left| {{e_\tau }\left( {{t_i}} \right)} \right| \leq C{e^{L\left( {{t_i}-{t_0}} \right)}}\max \limits_{t \in \left[ {{t_0},{t_i}} \right]} \left| {{\ell _\tau }\left( t \right)} \right|

wobei die Lipschitzkonstante von f bezüglich y ist. Sie kommt aus der Bedingung

\left\| {f\left( {t,{y^{\left( 1 \right)}}} \right)-f\left( {t,{y^{\left( 2 \right)}}} \right)} \right\| \leq L\left\| {{y^{\left( 1 \right)}}-{y^{\left( 2 \right)}}} \right\|\quad \forall t \in \left[ {{t_0},T} \right]\quad \forall {y^{\left( 1 \right)}},{y^{\left( 2 \right)}} \in {\mathbb{R}^n}

Ist y von geringerer Glattheit, also y \notin {C^{p+1}}\left[ {{t_0},T} \right], so ist die Anwendung eines Einschrittverfahrens der Konsistenzordnung p nicht zu empfehlen. Die Genauigkeit ist dann nicht von Ordnung p. Ein Einschrittverfahren niedrigerer Konsistenzordnung hat dann die gleiche Konvergenzordnung und ist billiger.

1.2.12 Vergleich des Fehlers mit dem physikalischen Modell

Wenn L und {t_i}-{t_0} beide nicht sehr groß sind, dann sollten wir einen kleinen globalen Fehler erhalten, wenn der lokale Fehler auch klein ist. Falls die Lipschitzkonstante L groß ist, so ist die Differentialgleichung instabil. Dann wird ein Fehler im Anfangswert, aber auch jeder lokale Diskretisierungsfehler, deutlich verstärkt. Die Konstante L ist allerdings im Allgemeinen nicht bekannt.
Ist L groß, so ist das eine schlechte Eigenschaft des mathematischen Modells, und es ist dringend notwendig zu überprüfen, ob diese Eigenschaft im realen physikalischen Problem wieder auftritt. Wenn ja, so ist dieses Anwendungsproblem mit numerischen Methoden nur sehr schlecht zugänglich. Wenn nein, so sollte die Modellierung überprüft werden.
Ein weiterer Verstärkungsfaktor für den Fehler tritt auf, wenn wir über sehr lange Zeiten rechnen, d.h. wenn {t_i}-{t_0} sehr groß ist. In diesem Fall addieren sich auch viele kleine, in jedem Schritt gemachte Fehler. Man sollte überprüfen, ob man wirklich über so lange Zeiten auf diese Weise rechnen will.

Ähnliche Artikel

Kommentar verfassen