Euler-Verfahren - Ingenieurmathematik (i939)

Das explizite und das implizite Euler-Verfahren setzen die naheliegende Idee um, den Ableitung $y'$ durch einen Differenzenquotient auf einen äquidistanten Gitter zu ersetzen (vgl. numerische Differentiation). Damit gehören sie zu den einfachsten Verfahren für das numerische Lösen von gewöhnlichen Differentialgleichungen und Systemen von Differentialgleichungen.

Das explizite Euler-Verfahren ist praktisch kaum relevant, da die Lösungen zu ungenau sind. Das implizite Euler-Verfahren ist hingegen praktisch einsetzbar.

6.5.1Explizites Euler-Verfahren¶

Führen das Verfahren zunächst für eine einzelne Gleichung

y'=f(x,y),\qquad y(x_0)=y_0

(1)

einschließlich Anfangswert ein.

6.5.1.1Formulierung¶

Wählen im Intervall $[x_0,\infty)$ ein äquidistantes Gitter mit Stützstellenabstand $h>0$ :

x_k:=x_0+k\,h,\qquad k=0,1,2,\ldots.

(2)

Für $k=0,1,2,\ldots$ können damit $y'(x_k)$ durch den üblichen Vorwärtsdifferenzenquotient ersetzen und erhalten damit für die Differentialgleichung an den Stellen $x_0,x_1,x_2,\ldots$ :

\frac{y(x_{k+1})-y(x_k)}{h}=f(x_k,y(x_k)).

(3)

Mit

y_{k}:=y(x_k),\qquad k=1,2,\ldots,

(4)

liefert Umstellen nach $y(x_{k+1})$ die Iterationsvorschrift für das explizite Euler-Verfahren:

y_{k+1}=y_k+h\,f(x_k,y_k),

(5)

wobei $y_0$ durch die Anfangsbedingung $y(x_0)=y_0$ der Differentialgleichung gegeben ist.

Dieses Verfahren folgt ausgehend vom Punkt $(x_0,y_0)$ schrittweise dem Richtungsfeld (IDV 630).

6.5.1.2Eigenschaften¶

Das explizite Euler-Verfahren liefert nur für kleine Schrittweiten $h$ und auch nur für wenige Schritte brauchbare Lösungen. Je größer $h$ ist oder je mehr Schritte ausgeführt werden, desto größer wird die Differenz zwischen exakter und numerisch berechneter Lösung. Man kann zeigen, dass stets

|y(x_k)-y_k|\leq c_1\,h\,\bigl(\rme^{c_2\,h\,k}-1\bigr)

(6)

gilt mit Konstanten $c_1,c_2>0$ die nur von $f$ abhängen. Diese Abschätzung liefert zwei Erkenntnisse:

Für feste Schrittweite $h$ kann die Differenz zwischen berechnetem und exaktem Funktionswert exponentiell bezüglich der Schrittanzahl wachsen.
Betrachtet man immer kleinere Schrittweiten $h\to 0$ , wählt dabei aber $k$ in Abhängigkeit von $h$ so, dass $x_k$ stets die selbe Stelle $x$ beschreibt (also $h\,k$ konstant ist), so geht der Fehler bei $|y(x)-y_k|$ gegen Null. Genauer gilt
$|y(x)-y_k|\leq c_3\,h\qquad\text{bei }h\to 0.$
(7)
Halbieren der Schrittweite führt also beispielsweise zum Halbieren des Fehlers.

Beispiel

Wollen das Anfangswertproblem

y'=\sin x-\frac{y}{t},\qquad y(1)=1

(8)

numerisch mit dem explizten Euler-Verfahren lösen. Die exakte Lösung (zu Vergleichszwecken) ist

y(x)=\frac{\sin x-x\,\cos x+1-\sin 1+\cos 1}{x}.

(9)

import matplotlib.pyplot as plt
import numpy as np

f = lambda x, y: np.sin(x) - y / x
x0 = 1
y0 = 1

x_max = x0 + 8  # Intervall, auf dem die Lösung berechnet werden soll
h = 0.5  # Schrittweite

# Listen für berechnete Werte
x = [x0]  # zunächst nur Startwert
y = [y0]

# Berechnung
while x[-1] < x_max:
    x_new = x[-1] + h
    y_new = y[-1] + h * f(x[-1], y[-1])
    x.append(x_new)
    y.append(y_new)

# Plot
fig, ax = plt.subplots(figsize=(10, 5))

sol_x = np.linspace(x0 - 0.5, x_max + 0.5, 100)
sol_y = (np.sin(sol_x) - sol_x * np.cos(sol_x) + 1 - np.sin(1) + np.cos(1)) / sol_x
ax.plot(sol_x, sol_y, '-r', label='exakte Lösung')
ax.plot(x0, y0, 'or', label='Anfangsbedingung')

ax.plot(x, y, '-b', label='numerische Lösung')

ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
plt.show()

Das Abdriften der numerischen Lösung ist deutlich erkennbar.

Das Abdriften kann extreme Züge annehmen bis hin zur völligen Unbrauchbarkeit der Ergebnisse.

Beispiel

Wollen das Anfangswertproblem

y'=-\frac{5}{2}\,y,\qquad y(0)=1

(10)

numerisch lösen. Die exakte Lösung (zu Vergleichszwecken) ist

y(x)=\rme^{-\frac{5}{2}\,x}.

(11)

f = lambda x, y: -5/2 * y
x0 = 0
y0 = 1

x_max = x0 + 8  # Intervall, auf dem die Lösung berechnet werden soll
h = 0.9  # Schrittweite

# Listen für berechnete Werte
x = [x0]  # zunächst nur Startwert
y = [y0]

# Berechnung
while x[-1] < x_max:
    x_new = x[-1] + h
    y_new = y[-1] + h * f(x[-1], y[-1])
    x.append(x_new)
    y.append(y_new)

# Plot
fig, ax = plt.subplots(figsize=(10, 5))

sol_x = np.linspace(x0 - 0.5, x_max + 0.5, 100)
sol_y = np.exp(-5/2 * sol_x)
ax.plot(sol_x, sol_y, '-r', label='exakte Lösung')
ax.plot(x0, y0, 'or', label='Anfangsbedingung')

ax.plot(x, y, '-b', label='numerische Lösung')

ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
plt.show()

Die numerische Lösung oszilliert mit zunehmender Schrittanzahl immer stärker und hat mit der tatsächlichen Lösung nichts mehr zu tun.

6.5.1.3Systeme¶

Für Systeme

{\begin{align} y_1'&=f_1(x,y_1,\ldots,y_n),\\ &\vdots\\ y_n'&=f_n(x,y_1,\ldots,y_n).\\ \end{align}}

(12)

von Differentialgleichungen mit Anfangswerten

y_1(x_0)=y_{1,0},\quad\ldots,\quad y_n(x_0)=y_{n,0}

(13)

kann das explizite Euler-Verfahren auf jede Gleichung angewendet werden. Man erhält die Iterationsvorschrift

{\begin{bmatrix}y_{1,k+1}\\\vdots\\y_{n,k+1}\end{bmatrix}} ={\begin{bmatrix}y_{1,k}\\\vdots\\y_{n,k}\end{bmatrix}} +h\,{\begin{bmatrix}f_1(x_k,y_{1,k},\ldots,y_{n,k})\\\vdots\\f_n(x_k,y_{1,k},\ldots,y_{n,k})\end{bmatrix}}

(14)

für $k=1,2,\ldots$ .

Gleichungen höherer Ordnung können stets als Systeme erster Ordnung geschrieben werden, sodass das explizite Euler-Verfahren auch für Gleichungen höherer Ordnung einsetzbar ist.

6.5.2Implizites Euler-Verfahren¶

6.5.2.1Formulierung¶

Nutzt man statt Vorwärtsdifferenzen Rückwärtsdifferenzen, so erhält man

\frac{y(x_{k+1})-y(x_k)}{h}=f(x_{k+1},y(x_{k+1})),\quad k=0,1,2,\ldots

(15)

aus der Differentialgleichung; also den Zusammenhang

y_{k+1}=y_k+h\,f(x_{k+1},y_{k+1}).

(16)

Zur Berechnung von $y_{k+1}$ aus $y_k$ ist somit ein (nichtlineares) Gleichungssystem zu lösen, da die gesuchte Größe $y_{k+1}$ nur implizit durch eine Gleichung beschrieben wird.

Das implizite Euler-Verfahren wählt den Wert $y_{k+1}$ gerade so, dass man bei einem Schritt entgegen des Richtungsfeldes der Differentialgleichung gerade den vorhergehenden Punkt $(x_k,y_k)$ erreicht (IDV 635).

6.5.2.2Eigenschaften¶

Beim Fehlerverhalten des impliziten Euler-Verfahrens gelten grundsätzlich die gleichen Aussagen wie beim expliziten Euler-Verfahren. Allerdings sind die Situationen, in den des exponentielle Fehlerwachstum tatsächlich auftritt deutlich seltener. Diese Situationen kann man recht konkret angeben, was wir hier aber nicht tun.

Beispiel

Lösen das Anfangswertproblem

y'=\sin x-\frac{y}{t},\qquad y(1)=1

(17)

aus dem obigen Beispiel numerisch mit dem impliziten Euler-Verfahren.

import scipy

f = lambda x, y: np.sin(x) - y / x
x0 = 1
y0 = 1

x_max = x0 + 8  # Intervall, auf dem die Lösung berechnet werden soll
h = 0.5  # Schrittweite

# Listen für berechnete Werte
x = [x0]  # zunächst nur Startwert
y = [y0]

# Berechnung
while x[-1] < x_max:
    x_new = x[-1] + h
    y_new = scipy.optimize.newton(
        lambda yk1: yk1 - y[-1] - h * f(x_new, yk1),
        y[-1],  # Startwert
        lambda yk1: 1 - h * (-1 / x_new)  # Ableitung
    )
    x.append(x_new)
    y.append(y_new)

# Plot
fig, ax = plt.subplots(figsize=(10, 5))

sol_x = np.linspace(x0 - 0.5, x_max + 0.5, 100)
sol_y = (np.sin(sol_x) - sol_x * np.cos(sol_x) + 1 - np.sin(1) + np.cos(1)) / sol_x
ax.plot(sol_x, sol_y, '-r', label='exakte Lösung')
ax.plot(x0, y0, 'or', label='Anfangsbedingung')

ax.plot(x, y, '-b', label='numerische Lösung')

ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
plt.show()

Das Abdriften der numerischen Lösung ist deutlich erkennbar.

Beispiel

Lösen das Anfangswertproblem

y'=-\frac{5}{2}\,y,\qquad y(0)=1.

(18)

f = lambda x, y: -5/2 * y
x0 = 0
y0 = 1

x_max = x0 + 8  # Intervall, auf dem die Lösung berechnet werden soll
h = 0.9  # Schrittweite

# Listen für berechnete Werte
x = [x0]  # zunächst nur Startwert
y = [y0]

# Berechnung
while x[-1] < x_max:
    x_new = x[-1] + h
    y_new = scipy.optimize.newton(
        lambda yk1: yk1 - y[-1] - h * f(x_new, yk1),
        y[-1],  # Startwert
        lambda yk1: 1 - h * (-5/2)  # Ableitung
    )
    x.append(x_new)
    y.append(y_new)

# Plot
fig, ax = plt.subplots(figsize=(10, 5))

sol_x = np.linspace(x0 - 0.5, x_max + 0.5, 100)
sol_y = np.exp(-5/2 * sol_x)
ax.plot(sol_x, sol_y, '-r', label='exakte Lösung')
ax.plot(x0, y0, 'or', label='Anfangsbedingung')

ax.plot(x, y, '-b', label='numerische Lösung')

ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
plt.show()

Im Gegensatz zum expliziten Euler-Verfahren treten hier keine immer stärker werdenden Oszillationen auf.

Der Preis für das bessere Fehlerverhalten ist die Notwendigkeit, in jedem Schritt eine nichtlineare Gleichung zu Lösen. Da allerdings mit dem vorhergehenden Funktionswert $y_k$ ein guter Startwert zur Verfügung steht, sind meist nur sehr wenige Newton-Iterationen nötig.

6.5.2.3Systeme¶

Analog zum expliziten Euler-Verfahren kann auch das implizite Euler-Verfahren zum Lösen von Differentialgleichungssystemen (und damit von Gleichungen höherer Ordnung) eingesetzt werden. In jedem Schritt ist dann ein nichtlineares Gleichungssystem zu lösen.

6.5 Euler-Verfahren

6.5.1Explizites Euler-Verfahren¶

6.5.1.1Formulierung¶

6.5.1.2Eigenschaften¶

6.5.1.3Systeme¶

6.5.2Implizites Euler-Verfahren¶

6.5.2.1Formulierung¶

6.5.2.2Eigenschaften¶

6.5.2.3Systeme¶