import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Konfiguriere matplotlib um schönere Plots zu erhalten
plt.rcParams.update({
"interactive": True,
"figure.dpi": 200,
"figure.figsize": (10, 6),
"font.family": "serif",
"mathtext.fontset": "cm",
"axes.grid": True,
"axes.labelsize": 14,
"grid.linewidth": .5,
"grid.alpha": .5,
"lines.linewidth": 1,
})
def euler(f, h, t, y):
return y + h * f(t, y)
def rk2(f, h, t, y):
k1 = f(t, y)
k2 = f(t + h/2, y + h/2 * k1)
return y + h * k2
def rk4(f, h, t, y):
k1 = f(t, y)
k2 = f(t + h/2, y + h/2 * k1)
k3 = f(t + h/2, y + h/2 * k2)
k4 = f(t + h, y + h * k3)
return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)
def simulate(f, y0, h, steps, method):
"""
Simulates a given function f(t, y) = dy/dt starting from (t=0, y=y0) for the
given number of steps and step size h using the specified integration method.
Returns: The calculated times and y-values.
"""
t = np.empty(steps)
t[0] = 0
y = np.empty(steps)
y[0] = y0
for i in range(steps-1):
t[i+1] = t[i] + h
y[i+1] = method(f, h, t[i], y[i])
return t, y
def compare_methods(f, y0, h, steps, analytic=None):
methods = {"Euler": euler, "RK2": rk2, "RK4": rk4}
# Plot der Werte
plt.figure(dpi=200)
plt.title("Approximationen im Vergleich zur analytischen Lösung")
for label, method in methods.items():
t, y = simulate(f, y0, h, steps, method)
plt.plot(t, y, ".--", label=label)
if analytic is not None:
t = np.linspace(0, h * (steps-1), 200)
y = analytic(t)
plt.plot(t, y, label="Analytische Lösung")
plt.xlabel("$t$")
plt.ylabel("$y$")
plt.legend()
plt.show()
# Plot der Fehler im Vergleich zur analytischen Lösung, aufgetragen gegen die Zeit
if analytic is not None:
plt.figure(dpi=200)
plt.title("Entwicklung des Approximationsfehlers mit der Zeit")
for label, method in methods.items():
t, y = simulate(f, y0, h, steps, method)
delta_y = np.abs(y - analytic(t))
# Nutze [1:] um den Anfangswert zu ignorieren, da dieser nicht fehlerbehaftet ist
plt.plot(t[1:], delta_y[1:], ".--", label=label)
plt.xlabel("$t$")
plt.ylabel("$\\Delta y$")
plt.yscale("log")
plt.legend()
plt.show()
# Plot der Fehler im Vergleich zur analytischen Lösung, aufgetragen gegen die Schrittweite
# (Konvergenz der Approximationen für h -> 0)
plt.figure(dpi=200)
plt.title("Konvergenz der Approximationen")
h_arr = np.logspace(-3, -1, 20) # Verteile 20 Werte für h logarithmisch zwischen 10^-3 und 10^0 = 1
for label, method in methods.items():
delta_y_list = []
for h in h_arr:
t, y = simulate(f, y0, h, steps, method)
end_delta_y = np.abs(y[-1] - analytic(t[-1]))
delta_y_list.append(end_delta_y)
plt.plot(h_arr, delta_y_list, ".--", label=label)
plt.xlabel("$h$")
plt.ylabel("$\\Delta y$")
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.show()
def f1(t, y):
return y
compare_methods(f1, y0=1, h=.5, steps=16, analytic=np.exp)
def f2(t, y):
return -t * y
# Hinweis: Hier wird mithilfe des lambda-Befehls eine Funktion erstellt.
# Mehr dazu: https://www.programiz.com/python-programming/anonymous-function
compare_methods(f2, y0=1, h=.5, steps=10, analytic=lambda t: np.exp(-t**2 / 2))