Voraussetzungen

  • Schwingungen

  • Numerische Integration analytischer Funktionen

Lerninhalte

  • Explizite Integrationsverfahren von gewöhnlichen Differentialgleichungen

Integration von ODEs mit expliziten Verfahren#

Gewöhnliche Differentialgleichungen (Ordinary Differential Equations, ODEs) können oft nicht analytisch gelöst werden.

Während die harmonische Schwingungsgleichung

\[\ddot x + \omega_0^2\cdot x = 0\]

eine analytische Lösung besitzt, hat das inverse Tangentenproblem (auch Traktrix genannt)

\[y' \cdot \sqrt{d^2-y^2} = - y\]

schon keine klare analytische Lösung mehr. In der Praxis begegnen uns viele Fälle, in denen sich zwar eine Differentialgleichung aufstellen, aber nicht analytisch lösen lässt.

Wichtige Begriffe#

Explizite Formulierungen#

Die Traktrix können wir nach \(y'\) umformen zu

\[y'(x) = - \frac{y(x)}{\sqrt{d^2-y^2(x)}}.\]

Wenn wir die ODE in der Form

\[y'(x) = f(x,y(x))\]

schreiben, sprechen wir von einer expliziten Formulierung. Diese Form ist besonders interessant für uns, denn wir können unsere Kenntnisse aus der numerischen Integration auf sie anwenden. Das ist besonders hilfreich, wenn wir die Funktion, wie die Traktrix, nicht analytisch integrieren können.

Schrittweise Integration#

Wie die numerische Integration bekannter Funktionen können wir \(y'(x)\) über das Teilintervall \([x_i,x_{i+1}]\) integrieren.

\[y(x_{i+1})-y(x_i) = \int_{x_i}^{x_{i+1}}f(x,y(x))dx\]

Durch Umstellen nach \(y_{i+1} = y(x_{i+1})\) erhalten wir nach jedem Teilintervall einen neuen Funktionswert für \(y_{i+1}\).

\[y_{i+1} = y_i + \int_{x_i}^{x_{i+1}} f(x,y(x))dx\]

Für die schrittweise Integration führen wir außerdem ein

\[h = x_{i+1} - x_i.\]

Anfangswertprobleme#

Mit der Annahme, dass \(y_{i+1} = f(y_i)\) ist, müssen wir ein \(y_0\) festlegen, von dem aus wir starten. Für das Taschenuhrenbeispiel ist es einfach \((0,d)\). Für andere ODEs, wie die Schwingungsgleichung, müssen wir Anfangswerte sinnvoll festlegen. Für das Pendel heißt das, wie hoch es anfangs hängt und wie schnell es sich bewegt. Bei ODEs sprechen wir daher auch von Anfangswertproblemen (AWPs oder IVPs, Initial Value Problems).

Jetzt müssen wir noch das Integral geeignet approximieren.

Aufgabe 1: Linke Rechteckregel#

Für die Integration bekannter Funktionen entspricht die linke Rechteckregel einem Rechteck unter dem Punkt \(x_i\):

../../_images/quadratur_rechteck_x.png

Fig. 2 Integral einer bekannten Funktion mit der linken Rechteckregel.#

Bei der Integration einer ODE entspricht die linke Rechteckregel dem Vektorpfeil (der Ableitung) im Punkt \((x_i,y_i)\):

../../_images/quadratur_rechteck_xdot.png

Fig. 3 Integration einer Differentialgleichung mit der linken Rechteckregel. Das Vektorfeld (rote Pfeile), das die Differentialgleichung an beliebigen Punkten \((x,y)\) beschreibt, wird immer am linken Rand des diskretisierten Intervalls (\(x_0, x_1, ...\)) ausgewertet, um auf den nächsten Startwert zu schließen (blaue Pfeile). Dabei wird hier die ideale Lösung (blau gestrichelt) unterschätzt.#

Wenden Sie die linke Rechteckregel

\[\int_{x_i}^{x_{i+1}}f(x,y(x)) dx \approx h \cdot f(x_i,y_i)\]

auf die Traktrix an. In dem Code ist in Zeile 8 bereits x=xspan(i) definiert. Das kommt in der Traktrix nicht vor, ist aber für andere Gleichungen \(f(x,y)\) sehr relevant, wie Sie weiter unten sehen werden.

Was passiert bei \(x=0\)? Wie können Sie dem Problem ausweichen?

Hinweis

Dieses Vorgehen entspricht dem expliziten Eulerverfahren, das in dem Kapitel Stabilität des expliziten versus impliziten Eulerverfahrens weiter dikutiert wird. Die rechte Rechteckregel entspricht dem impliziten Eulerverfahren, das in dem Kapitel Integration von ODEs mit impliziten Verfahren behandelt wird.

% your code here
y0    = 1;
n     = 20; % number of timesteps
xmax  = 10;
xspan = linspace(0,xmax,n);
y     = zeros(1,n);
for i = 1:n-1
    x    = xspan(i);
    y(i) = ...
end

plot(xspan,y)

Mehrstufige Verfahren#

Aufgabe 2: Mittelpunktsregel#

Um eine größere Genauigkeit zu erreichen, wenden Sie die Mittelpunktsregel an.

\[\int_{x_i}^{x_{i+1}}f\left(x,y(x)\right)dx \approx h \cdot f\left(x_{i+\frac{1}{2}},y\left(x_i\right)\right)\]
../../_images/quadratur_mittelpunkt_x.png

Fig. 4 Integral einer bekannten Funktion mit der Mittelpunktregel.#

Hieraus ergibt sich für die schrittweise Integration

\[y_{i+1} \approx y_i + h \cdot f\left(\frac{x_i+x_{i+1}}{2},y\left(x_i\right)\right).\]

Hierbei ist \(\frac{x_i+x_{i+1}}{2}\) gleichbedeutend mit \(x_{i+\frac{1}{2}}\), beziehungsweise \(x_i + \frac{1}{2}h\). Wir verwenden weiter unten vor allem die letzte Form.

../../_images/quadratur_mittelpunkt_xdot.png

Fig. 5 Integration einer Differentialgleichung mit der Mittelpunktregel. Das Vektorfeld (rote Pfeile), das die Differentialgleichung an beliebigen Punkten \((x,y)\) beschreibt, wird zwischen zwei Punkten des diskretisierten Intervalls (\(x_{1/2}, x_{3/2}, ...\)) ausgewertet (dicke rote Pfeile), um auf die Steigung zu schließen (blaue Pfeile). Dabei wird die Lösung schon deutlich besser angenähert.#

Hinweis

Kopieren Sie Ihren Code in die neue Codebox oder auf Ihrem Rechner in eine neue Datei, um die Lösungen der verschiedenen Verfahren später vergleichen zu können. Übernehmen Sie die Anpassung für x.

% your code here

Streng genommen können Sie für die Approximation nicht \(y(x_i)\) verwenden, sondern müssen \(y\left(x_{i+1/2}\right)\) schätzen. Verwenden Sie dazu die linke Rechteckregel von oben.

\[y\left(x_{i+\frac{1}{2}}\right) = y_{i+\frac{1}{2}} = y_i + \frac{h}{2} f(x_i,y_i)\]

Wir führen zur besseren Übersicht \(K_1\) für die linke und \(K_2\) für die mittlere Stützstelle ein und kommen zu

\[K_1 := f(x_i,y_i)\]
\[K_2 := f\left(x_i+\frac{1}{2}h,y_i + \frac{h}{2} K_1\right)\]
\[y_{i+1} \approx y_i + h \cdot K_2.\]

Übernehmen Sie die verbesserte Stützstelle \(K_2\) in Ihrem Code. In der Abbildung sehen Sie, wie sich die verbesserte Ableitung auswirkt. Die kleinen blauen Pfeile beschreiben die Ableitung \(f(x,y(x))\) an ausgewählten Stellen im Raum ausgewertet, da die ODE von der Zeit und dem Funktionswert abhängen kann.

Name of image

Abbildung 6: Lösung der Traktrix (orange) mit explizitem Eulerverfahren (dunkelblau) und einem Runke-Kutta-Verfahren zweiter Ordnung (lila). Die Differentialgleichung kann an beliebigen Punkten \((x,y)\) ausgewertet werden (hellblaue Pfeile). Für die RK Lösung wird zuerst die Steigung am linken Rand ausgewertet (roter Pfeil), anschließend an dessen Hälfte die Steigung ermittelt (grüner Pfeil) und zuletzt diese Steigung am Startpunkt angesetzt.

% your code here

Anwendung auf höherstufige Verfahren#

Formulieren wir die Simpsonregel

\[\int_{x_i}^{x_{i+1}} f(x,y(x)) dx \approx \frac{x_{i+1} - x_i}{6}\left(f(x_i,y_i)+4 f\left(x_{i+\frac{1}{2}},y_i\right)+f(x_{i+1},y_i)\right)\]

um zu

\[y_{i+1} \approx y_i + h \cdot \left(\frac{1}{6} \cdot f\left(x_i,y_i\right) + \frac{4}{6} \cdot f\left(x_i+\frac{1}{2}h,y_i\right) + \frac{1}{6} \cdot f\left(x_i+h,y_i\right)\right).\]

Außerdem benennen wir die drei verwendeten Stufenwerte als \(K_1, K_2, K_3\), also

\[\begin{split}\begin{align}K_1 &= f\left(x_i,y_i\right)\\ K_2 &= f\left(x_i+\frac{1}{2}h,y_i\right)\\ K_3 &= f\left(x_i+h,y_i\right).\end{align}\end{split}\]

Wir können das aber noch besser: Bisher haben wir angenommen, dass wir \(y_i\) auch für die mittlere Stützstelle \(K_2 = f\left(x_i+\frac{1}{2}h,y_i\right)\) und die rechte Stützstelle \(K_3 = f\left(x_i+h,y_i\right)\) verwenden können. Wir haben aber schon eine Approximation, für die Stelle \(x_{i+\frac{1}{2}}\). Mit \(\frac{h}{2}\) multipliziert, können wir \(K_1\) nutzen, um mit der linken Rechteckregel \(y_{i+\frac{1}{2}}\) zu approximieren:

\[K_2 = f\left(x_i+\frac{1}{2}h,y_i + \frac{h}{2}\cdot f(x_i,y_i)\right) = f\left(x_i+\frac{1}{2}h,y_i + h\cdot\frac{1}{2}K_1\right).\]

Genauso verwenden wir für \(K_3\) zusätzlich auch \(K_2\) als Stützstelle:

\[\begin{split}\begin{align} & \\ & \end{align} K_3 = f\left(x_i + h, y_i + h \cdot \left(-1 \cdot K_1 + 2 \cdot K_2 \right)\right).\end{split}\]

Setzen wir das zurück in die Simpsonregel ein, ergibt sich ein Runge-Kutta-Verfahren der Ordnung 3:

\[y_{i+1} = y_i + h \cdot \left(\frac{1}{6}K_1 + \frac{4}{6}K_2 + \frac{1}{6}K_3\right)\]

In der Abbildung sehen Sie, wie die einzelnen Faktoren das Ergebnis für eine Traktrix ergeben. Da die Auswertungsstelle für \(K_3\) unschön ist, finden Sie darunter auch die Aufstellung der Faktoren für eine Kosinus-Funktion.

../../_images/vectorfield_animated_traktrix.gif

Fig. 7 Lösung der Traktrix (orange) mit explizitem Eulerverfahren (dunkelblau) und einem Runke-Kutta-Verfahren dritter Ordnung (lila). Die Differentialgleichung kann an beliebigen Punkten \((x,y)\) ausgewertet werden (hellblaue Pfeile). Für die RK Lösung werden in verschiedenen Schritten die Faktoren \(K_1, K_2, K_3\) ermittelt und zuletzt gewichtet mit \(\frac{1}{6},\frac{4}{6},\frac{1}{6}\) aneinandergereiht, um zum ersten Schritt zu gelangen.#

../../_images/vectorfield_animated_cos.gif

Fig. 8 Lösung der Funktion \(y' = f(x,y) = \sin(x)\) mit der Lösung \(y = \cos(x)\) (orange) mit explizitem Eulerverfahren (dunkelblau) und einem Runke-Kutta-Verfahren dritter Ordnung (lila). Die Differentialgleichung kann an beliebigen Punkten \((x,y)\) ausgewertet werden (hellblaue Pfeile). Für die RK Lösung werden in verschiedenen Schritten die Faktoren \(K_1, K_2, K_3\) ermittelt und zuletzt gewichtet mit \(\frac{1}{6},\frac{4}{6},\frac{1}{6}\) aneinandergereiht, um zum ersten Schritt zu gelangen.#

Hinweis

Verwechseln Sie nicht die mehrstufigen Verfahren mit Mehrschrittverfahren. Letztere verwenden zusätzlich zu \(y_i\) auch vorherige Auswertungen \(y_{i-c\cdot h}\).

Aufgabe 3: Implementierung des Runge-Kutta-Verfahrens 3. Ordnung#

Passen Sie Ihren Code so an, dass zur Lösung anstelle des expliziten Eulerverfahrens das RK-Verfahren 3. Ordnung verwendet wird.

% your code here

Aufgabe 4: Ordnungen der Verfahren#

Die Ordnung \(\mathcal{O(h)}\) gibt an, wie sich ein Fehler bei Verkleinerung der Schrittweite ändern. Ein Verfahren 2. Ordnung, \(\mathcal{O}(h^2)\) hat also bei Halbierung der Schrittweite nur noch ein Viertel des Fehlers.

Ermitteln Sie nun experimentell die Ordnungen der von Ihnen entwickelten Verfahren anhand der Traktrix:

  • explizites Eulerverfahren (Linke Rechteckregel)

  • Mittelpunktsregel ohne Schätzung von \(y_{i+\frac{1}{2}}\)

  • Runge-Kutta-Verfahren 2. Ordnung (Mittelpunktsregel mit Schätzung von \(y_{i+\frac{1}{2}})\)

  • Runge-Kutta-Verfahren 3. Ordnung (Simpsonregel)

% your code here