Numerische Methoden für Anfangswertprobleme¶

Lernziele¶

Sie können das explizite und implizite Euler-Verfahren zur Lösung skalarer Anfangswertprobleme

\[\begin{split}\begin{split} y'(x) &= f(x,y(x)),\\ y(x_0) &= y_0, \end{split}\end{split}\]

anwenden.
Sie verstehen den Unterschied zwischen den beiden Verfahren.
Sie können modularen Code erstellen, um eine hohe Wiederverwendbarkeit zu erhalten.
Sie können die Methoden systematisch testen.

Theorie¶

Die Theorie zum Praktikum ist im Skript Kapitel 3.1 - 3.2.2 gegeben.

Das Praktikum ist wie folgt aufgebaut:

Implementieren explizites Euler-Verfahren
- Anwenden auf Testproblem mit bekannter analytischer Lösung
- Kontrolle der Konvergenz (1. Ordnung vgl. Satz 3.1 im Skript)
Implementieren implizites Euler-Verfahren
- Anwenden auf Testproblem mit bekannter analytischer Lösung
- Kontrolle der Konvergenz (1. Ordnung vgl. Satz 3.1 im Skript)
Anwendung auf Anfangswertproblem mit unbekannter analytischer Lösung
- Unterschied zwischen den beiden Verfahren

Aufträge¶

Implementieren Sie das explizite Euler-Verfahren

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

(Formel 3.2 im Skript) als Funktion. Übergeben Sie die rechte Seite \(f(x,y)\) der Differentialgleichung als Parameter:
```
def explizitEuler(xend, h, x0, y0, f):
    x = [x0]
    y = [y0]

    while x[-1] < xend-h/2:
        <<snipp>>
    return np.array(x), np.array(y)
```
Testen Sie Ihr Programm mit dem Modellproblem

(1)¶\[\begin{split}y'(x) & = -4\,y(x)\quad x\in [0,2]\\ y(0) & = 1.\end{split}\]

Berechnen Sie dazu die analytische Lösung und den absoluten Fehler, visualisieren Sie diesen.

Ein sehr guter Test für numerische Verfahren ist die Kontrolle der Konvergenzordnung. Gehen Sie dazu wie folgt vor:

n = 10**np.linspace(2,5)
hs = 2/n
err = []
for h in hs:
    x, y = explizitEuler(2,h,0,1,f)
    err.append(np.linalg.norm(y-ya(x),np.inf)) # ya(x) ist die exakte Lösung

plt.loglog(hs,err,'-')
plt.xlabel('h')
plt.ylabel(r'$\max_k \|e(x_k,h)\|$')
plt.grid()
plt.show()

Implementieren Sie das implizite Euler-Verfahren

\[y_{k+1} = y_k + h\cdot f(x_{k+1}, y_{k+1})\]

(Formel 3.4 im Skript) als Funktion. Übergeben Sie die rechte Seite \(f(x,y)\) der Differentialgleichung sowie die partielle Ableitung \(\partial_y f(x,y)\) als Parameter. Ziel ist eine möglichst modulare Implementierung. Dazu implementieren wir das Newton-Verfahren aus der letzten Woche als Funktion einer generischen Abbildung \(G(s)\) deren Nullstelle gesucht wird. Die generische Abbildung \(G(s)\) ist im Fall des impliziten Euler-Verfahren gegeben durch

\[G(s) = s - y_k - h\cdot f(x_{k+1}, s).\]

Sie können sich am folgenden Template orientieren:

def implizitEuler(xend, h, x0 y0, f, df):
    x = [x0]
    y = [y0]

    # Verfahrensfunktion für implizit Euler
    def G(s, xk, yk):
        return <<snipp>>

    # Partielle Ableitung nach s der Verfahrensfunktion
    def dG(s, xk, yk):
        return <<snipp>>

    def newton(s, xk, yk, tol=1e-12, maxIter=20):
        delta = 10*tol
        for k in range(maxIter):
            delta = <<snipp>>
            s -= delta
            if np.abs(delta) < tol:
                break
        return s

    while x[-1] < xend-h/2:
        y.append(newton(y[-1],x[-1],y[-1]))
        x.append(x[-1]+h)
    return np.array(x), np.array(y)

Testen Sie Ihr Programm wieder mit dem Modellproblem (1) und kontrollieren / visualisieren Sie die Konvergenzordnung.
Berechnen Sie mit den beiden Verfahren auf dem Intervall \([0,2]\) numerisch Lösungen für das Anfangswertproblem

\[y'(x) + \frac{x^2}{y(x)} = 0, \quad y(0) = -4.\]
- Da die nichtlineare Differentialgleichung separabel ist, kann eine analytische Lösung berechnet werden. Berechnen Sie mit Hilfe der analytischen Lösung den absoluten Fehler für \(x=2\). Visualisieren Sie diesen in Abhängigkeit verschiedener Schrittweiten, gegeben durch die Anzahl Intervalle \(N=3^j\), \(j \in \{1,2,3,4,5,6,7,8\}\).
  
  (Zur Kontrolle \(y(2) = -4\sqrt{\frac{2}{3}}\))
- Visualisieren Sie Ihre numerische Lösung für \(N=3^8\) im Richtungsfeld der Differentialgleichung.
Verständnisfrage: Erklären Sie mit Hilfe des Richtungsfeldes den Unterschied zwischen dem expliziten und impliziten Euler-Verfahren.

Abgabe¶

Bitte geben Sie Ihre Lösungen bis spätestens vor dem nächsten Praktikum ab.