Im Unterpunkt Systeme von gekoppelten Differentialgleichungen und Differentialgleichungen zweiter Ordnung der Vorlesung 9 behandelten wir unter anderem das numerische Lösen von Systemen von gekoppelten Differentialgleichungen erster Ordnung mittels der C++ und Python Programmiersprachen. In dem C++ Programm DGL_2.cpp und dem Jupyter Notebook DGL_2.ipynb hatten wir exemplarisch die Lösungen eines Systems bestehend aus zwei gekoppelten Differentialgleichungen erster Ordnung bestimmt. Erweitern Sie das C++ Programm oder das Jupyter Notebook auf ein System bestehend aus drei gekoppelten Differentialgleichungen erster Ordnung und lösen Sie das folgende Differentialgleichung-System des Lorenz-Modells mittels C++ oder Python. \begin{eqnarray} \dot{x}(t) &=& \frac{d x}{dt} = \, \sigma \cdot \left( y - x \right) \\ \dot{y}(t) &=& \frac{d y}{dt} = \, r \cdot x - y - x \cdot z \\ \dot{z}(t) &=& \frac{d z}{dt} = \, x \cdot y - b \cdot z \end{eqnarray} Das Lorenz-System wurde im Jahre 1962 von dem Meteorologen Edward N. Lorenz als eine Idealisierung der hydrodynamischen Gleichungen der Erdatmosphäre entwickelt und stellt eines der am meisten untersuchten kontinuierlichen Modelle mit chaotischem Verhalten dar. Den Parameter $\sigma$ bezeichnet man als 'Prandtl-Zahl', $r$ als 'relative Rayleigh-Zahl' und $b$ ist eine Konstante mit $b \in (0,4)$.
Betrachten Sie speziell das Lorenz-Modell mit $\sigma = 10$ und $b = \frac{8}{3}$ und berechnen Sie zumindest für einen der folgenden $r$-Werte die numerische Lösung: $r = 0.9, 12, 22, 28$. Betrachten Sie dabei sowohl die zeitliche Entwicklung ausgehend von dem Startwert $x(0) = 0 \, , \,\, y(0) = 1$ und $z(0) = 0$ als auch zumindest einen weiteren Anfangswert. Visualisieren Sie Ihre Ergebnisse, indem sie die zeitliche Entwicklung der Funktionen $x(t) \, , \,\, y(t) \, , \,\, z(t)$ darstellen. Erstellen Sie zusätzlich ein Bild der Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ in einem dreidimensionalen $x,y,z$-Koordinatensystem. Verwenden Sie dazu z.B. das folgende Jupyter Notebook als Vorlage: Vorlage_Raumkurve.ipynb (View Notebook, Download Notebook).
Wir berechnen nun die numerische Lösung des oben beschriebenen Lorenz-Modells mit $\sigma = 10$ und $b = \frac{8}{3}$ in Python und visualisieren auch die numerischen Lösungen des C++ Programms Lorenz.cpp. Wir benutzen hierfür das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Notebook die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.integrate import solve_ivp
Wir definieren uns das System der DGLs als eine Funktion
def DGLsys(t,vec,r):
x, y, z = vec
dx_dt = 10*(y - x)
dy_dt = r*x - y - x*z
dz_dt = x*y - 8/3*z
return np.array([dx_dt,dy_dt,dz_dt])
und lösen das Anfangswertproblem für unterschiedliche Werte der 'Rayleigh-Zahl' $r$.
Aufgrund der nicht-linearen Bewegungsgleichung des Lorenz-Modells können sich in gewissen Parameterbereichen kleine Ungenauigkeiten exponentiell vergrößern, sodass eine genaue Vorhersage der Bewegung nach einiger Zeit nicht mehr möglich ist. Bei $r=28$ entsteht entsteht z.B. ein chaotischer bzw. seltsamen Attraktor (Lorenz-Attraktor). Im folgenden betrachten wir die folgenden unterschiedlichen r-Werte: $r = 0.9, 12, 22, 28$
Wir lösen das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,40]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points). In den folgenden beiden Simulationen setzen wir die relativen und absoluten Fehler-Toleranzen (rtol und atol) auf zwei unterschiedliche Werte: $\Delta_0=10^{-10}$ und $\Delta_1=10^{-13}$ (näheres siehe scipy.integrate.solve_ivp).
t_end = 40
N = 10000
t_val = np.linspace(0, t_end, N+1)
vec_init = [0,1,0]
set_r = 0.9
fehler_0 = 10**(-10)
fehler_1 = 10**(-13)
Loes_0 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_0, atol=fehler_0)
Loes_1 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Wir stellen uns die numerischen Lösungen grafisch dar:
params = {
'figure.figsize' : [11,5],
# 'text.usetex' : True,
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
}
matplotlib.rcParams.update(params)
plt.title(r"$\rm Python-Lösung \; mit \; r=0.9 \quad und \quad \Delta_1=10^{-13}$")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(Loes_1.t, Loes_1.y[0],c="blue", label="x");
plt.plot(Loes_1.t, Loes_1.y[1],c="green", label="y");
plt.plot(Loes_1.t, Loes_1.y[2],c="red", label="z");
plt.legend(loc='upper right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen dar (z.B. $\Delta x = x_{\Delta_0} - x_{\Delta_1}$).
plt.title(r"$\rm Python-Lösung \; mit \; r=0.9 \quad , \quad Abweichungen: \, \Delta \, = \, \Delta_0 \, - \, \Delta_1$")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_0.t, Loes_1.y[0]-Loes_0.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_0.t, Loes_1.y[1]-Loes_0.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_0.t, Loes_1.y[2]-Loes_0.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Die beiden Lösungen stimmen, innerhalb der angegebenen Fehler-Toleranzen, gut miteinander überein. Die zeitliche Entwicklung zeigt somit bei $r=0.9$ ein nicht-chaotisches, reguläres Verhalten.
Unter Verwendung des Jupyter Notebooks Vorlage_Raumkurve.ipynb stellen nun noch die Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ in einem drei-dimensionalen $x,y,z$-Koordinatensystem dar.
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=1, c = Loes_1.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=20)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm Python-Lösung \; mit \; r=0.9 \quad und \quad \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Betrachtet man sich die zeitliche Entwicklung ausgehend von anderen Anfangswerten (hier speziell $[x(0),y(0),z(0)] = [1,0,0.2]$ (blau) , $[0.5,0,0.5]$ (rot), $[0.5,1,0.5]$ (grün)) oder verändert sie nur um einen Epsilon-Wert (hier speziell $[x(0),y(0),z(0)] = [0,0.999,0]$ (gelb)), so entwickelt sich die Trajektorie ebenfalls hin zu ihrem einzigen gemeinsamen Attraktor, dem Ursprung $[0,0,0]$.
vec_init_eps = [0,1-10**(-3),0]
vec_init_a = [1,0,0.2]
vec_init_b = [0.5,0,0.5]
vec_init_c = [0.5,1,0.5]
vec_init = [0,1,0]
Loes_1_eps = solve_ivp(DGLsys, [0, t_end], vec_init_eps, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_a = solve_ivp(DGLsys, [0, t_end], vec_init_a, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_b = solve_ivp(DGLsys, [0, t_end], vec_init_b, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_c = solve_ivp(DGLsys, [0, t_end], vec_init_c, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=3, c = "black")
ax.scatter3D(Loes_1_eps.y[0],Loes_1_eps.y[1],Loes_1_eps.y[2], marker='o', s=1, c = "yellow")
ax.scatter3D(Loes_1_a.y[0],Loes_1_a.y[1],Loes_1_a.y[2], marker='o', s=1, c = "blue")
ax.scatter3D(Loes_1_b.y[0],Loes_1_b.y[1],Loes_1_b.y[2], marker='o', s=1, c = "red")
ax.scatter3D(Loes_1_c.y[0],Loes_1_c.y[1],Loes_1_c.y[2], marker='o', s=1, c = "green")
ax.view_init(azim=-75, elev=20)
plt.title(r"$\rm r=0.9 \quad , \quad unterschiedliche \, Anfangsbedingungen$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,40]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points) mittels des C++ Programms Lorenz.cpp. Die Ergebnisse des C++ Programms wurden in die Datei "Lorenz-r_0.9-init_0.dat" geschrieben. Im Folgenden stellen wir die Ergebnisse der numerischen Euler- und Runge Kutta RK4 Methoden dar.
data = np.genfromtxt("./Lorenz-r_0.9-init_0.dat")
plt.title("C++ Lösung (Euler Methode): r=0.9 , N=10000")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(data[:,1],data[:,2],c="blue", label=r"$\rm x$");
plt.plot(data[:,1],data[:,3],c="green", label=r"$\rm x$");
plt.plot(data[:,1],data[:,4],c="red", label=r"$\rm z$");
plt.legend(loc='lower right',fontsize=16);
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=0.9 , N=10000")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(data[:,1],data[:,5],c="blue", label=r"$\rm x$");
plt.plot(data[:,1],data[:,6],c="green", label=r"$\rm y$");
plt.plot(data[:,1],data[:,7],c="red", label=r"$\rm z$");
plt.legend(loc='lower right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen dar.
plt.title(r"C++ Lösung (Abweichungen $\Delta$ Euler vs. Runge Kutta RK4): r=0.9 , N=10000")
plt.xlabel("t")
plt.ylabel(r"$\rm\Delta x, \Delta y, \Delta z$")
plt.plot(data[:,1],data[:,2]-data[:,5],c="blue", label=r"$\rm \Delta x$");
plt.plot(data[:,1],data[:,3]-data[:,6],c="green", label=r"$\rm \Delta y$");
plt.plot(data[:,1],data[:,4]-data[:,7],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Runge Kutta RK4 Methode im Vergleich zur Python Lösung ($\Delta_1=10^{-13}$) ergibt sich:
plt.title(r"Abweichungen $\Delta$ Runge Kutta RK4 vs. Phython $\Delta_1=10^{-13}$): r=0.9 , N=10000")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm y_1,y_2,y_3$")
plt.plot(data[:,1],Loes_1.y[0]-data[:,5],c="blue", label=r"$\rm y_1$");
plt.plot(data[:,1],Loes_1.y[1]-data[:,6],c="green", label=r"$\rm y_2$");
plt.plot(data[:,1],Loes_1.y[2]-data[:,7],c="red", label=r"$\rm y_3$");
plt.legend(loc='lower right',fontsize=16);
Stellt man sich ausgehend von den simulierten Daten der Runge Kutta RK4 Methode die Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ in einem dreidimensionalen $x,y,z$-Koordinatensystem dar, so ergibt sich wieder das folgende Bild.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(data[:,5],data[:,6],data[:,7], marker='o', s=0.1, c = data[:,1], cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=0.9 , N=10000")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,15]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points) für $r = 12$.
t_end = 15
N = 10000
t_val = np.linspace(0, t_end, N+1)
set_r = 12
Loes_0 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_0, atol=fehler_0)
Loes_1 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
plt.title(r"$\rm Python-Lösung \; mit \; r=12 \quad und \quad \Delta_1=10^{-13}$")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(Loes_1.t, Loes_1.y[0],c="blue", label="x");
plt.plot(Loes_1.t, Loes_1.y[1],c="green", label="y");
plt.plot(Loes_1.t, Loes_1.y[2],c="red", label="z");
plt.legend(loc='upper right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen dar (z.B. $\Delta x = x_{\Delta_0} - x_{\Delta_1}$).
plt.title(r"$\rm Python-Lösung \; mit \; r=12 \quad , \quad Abweichungen: \, \Delta \, = \, \Delta_0 \, - \, \Delta_1$")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_0.t, Loes_1.y[0]-Loes_0.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_0.t, Loes_1.y[1]-Loes_0.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_0.t, Loes_1.y[2]-Loes_0.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Die beiden Lösungen stimmen, innerhalb der angegebenen Fehler-Toleranzen, gut miteinander überein. Die zeitliche Entwicklung zeigt somit auch bei $r=12$ ein nicht-chaotisches, reguläres Verhalten.
Die dreidimensionale Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ sieht wie folgt aus.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=0.1, c = Loes_1.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=20)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm Python-Lösung \; mit \; r=12 \quad und \quad \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Betrachtet man sich die zeitliche Entwicklung ausgehend von einem anderen Anfangswert (hier speziell $[x(0),y(0),z(0)] = [-10,12,30]$), so kann es geschehen, dass sich das System zu einem anderen Fixpunkt hin entwickelt (siehe folgende Abbildung). Verändert man hingegen den Anfangswert nur um einen Epsilon-Wert (hier speziell $[x(0),y(0),z(0)] = [0,0.999,0]$), so entwickelt sich die Trajektorie ebenfalls hin zu ihrem ursprünglichem Attraktor.
vec_init_a = [-10,12,30]
Loes_1_eps = solve_ivp(DGLsys, [0, t_end], vec_init_eps, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_a = solve_ivp(DGLsys, [0, t_end], vec_init_a, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=0.1, c = Loes_1.t, cmap=cm.gnuplot2)
ax.scatter3D(Loes_1_eps.y[0],Loes_1_eps.y[1],Loes_1_eps.y[2], marker='o', s=0.1, c = Loes_1_eps.t, cmap=cm.gnuplot2)
ax.scatter3D(Loes_1_a.y[0],Loes_1_a.y[1],Loes_1_a.y[2], marker='o', s=0.1, c = Loes_1_a.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=12 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,15]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points) mittels des C++ Programms Lorenz.cpp. Die Ergebnisse des C++ Programms wurden in die Datei "Lorenz-r_12-init_0.dat" geschrieben. Im Folgenden stellen wir die Ergebnisse der numerischen Euler- und Runge Kutta RK4 Methoden dar.
data = np.genfromtxt("./Lorenz-r_12-init_0.dat")
plt.title("C++ Lösung (Euler Methode): r=12 , N=10000")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(data[:,1],data[:,2],c="blue", label=r"$\rm x$");
plt.plot(data[:,1],data[:,3],c="green", label=r"$\rm x$");
plt.plot(data[:,1],data[:,4],c="red", label=r"$\rm z$");
plt.legend(loc='lower right',fontsize=16);
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=12 , N=10000")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(data[:,1],data[:,5],c="blue", label=r"$\rm x$");
plt.plot(data[:,1],data[:,6],c="green", label=r"$\rm y$");
plt.plot(data[:,1],data[:,7],c="red", label=r"$\rm z$");
plt.legend(loc='lower right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen dar.
plt.title(r"C++ Lösung (Abweichungen $\Delta$ Euler vs. Runge Kutta RK4): r=12 , N=10000")
plt.xlabel("t")
plt.ylabel(r"$\rm\Delta x, \Delta y, \Delta z$")
plt.plot(data[:,1],data[:,2]-data[:,5],c="blue", label=r"$\rm \Delta x$");
plt.plot(data[:,1],data[:,3]-data[:,6],c="green", label=r"$\rm \Delta y$");
plt.plot(data[:,1],data[:,4]-data[:,7],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Runge Kutta RK4 Methode im Vergleich zur Python Lösung ($\Delta_1=10^{-13}$) ergibt sich:
plt.title(r"Abweichungen $\Delta$ Runge Kutta RK4 vs. Phython $\Delta_1=10^{-13}$): r=12 , N=10000")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm y_1,y_2,y_3$")
plt.plot(data[:,1],Loes_1.y[0]-data[:,5],c="blue", label=r"$\rm y_1$");
plt.plot(data[:,1],Loes_1.y[1]-data[:,6],c="green", label=r"$\rm y_2$");
plt.plot(data[:,1],Loes_1.y[2]-data[:,7],c="red", label=r"$\rm y_3$");
plt.legend(loc='lower right',fontsize=16);
Die Ergebnisse der Runge Kutta RK4 Methode stimmen gut mit der Python Lösung überein, wohingegen die Euler Methode keine guten Ergebnisse lieferte und hier eine Erhöhung der Anzahl der Gitterpunkte nötig wäre um gute Ergebnisse zu erhalten.
Stellt man sich ausgehend von den simulierten Daten der Runge Kutta RK4 Methode die Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ in einem dreidimensionalen $x,y,z$-Koordinatensystem dar, so ergibt sich wieder das folgende Bild.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(data[:,5],data[:,6],data[:,7], marker='o', s=0.1, c = data[:,1], cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=12 , N=10000")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,50]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points) für $r = 22$.
t_end = 50
N = 10000
t_val = np.linspace(0, t_end, N+1)
set_r = 22
Loes_0 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_0, atol=fehler_0)
Loes_1 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Wir stellen uns die numerische Lösungen grafisch dar:
plt.title(r"$\rm Python-Lösung \; mit \; r=22 \quad und \quad \Delta_1=10^{-13}$")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(Loes_1.t, Loes_1.y[0],c="blue", label="x");
plt.plot(Loes_1.t, Loes_1.y[1],c="green", label="y");
plt.plot(Loes_1.t, Loes_1.y[2],c="red", label="z");
plt.legend(loc='upper right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen dar (z.B. $\Delta x = x_{\Delta_0} - x_{\Delta_1}$).
plt.title(r"$\rm Python-Lösung \; mit \; r=22 \quad , \quad Abweichungen: \, \Delta \, = \, \Delta_0 \, - \, \Delta_1$")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_0.t, Loes_1.y[0]-Loes_0.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_0.t, Loes_1.y[1]-Loes_0.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_0.t, Loes_1.y[2]-Loes_0.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Die beiden Lösungen stimmen, innerhalb der angegebenen Fehler-Toleranzen, gut miteinander überein. Man würde vermuten, dass die zeitliche Entwicklung somit bei $r=22$ ein nicht-chaotisches, reguläres Verhalten zeigt. Dies stimmt auch für diese Anfangsbedingung, jedoch bei anderen Anfangsbedingungen kann es in gewissen Zeitabschnitten zu einem vorübergehenden chaotischen Verhalten kommen ('transient chaos', siehe später).
Die dreidimensionale Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ sieht wie folgt aus.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=0.1, c = Loes_1.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=22 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Betrachtet man sich die zeitliche Entwicklung ausgehend von einem anderen Anfangswert (hier speziell $[x(0),y(0),z(0)] = [-10,12,30]$), so kann es geschehen, dass sich das System zu einem anderen Fixpunkt hin entwickelt (siehe folgende Abbildung). Verändert man hingegen den Anfangswert nur um einen Epsilon-Wert (hier speziell $[x(0),y(0),z(0)] = [0,0.999,0]$), so entwickelt sich die Trajektorie ebenfalls hin zu ihrem ursprünglichem Attraktor.
vec_init_a = [-10,12,30]
vec_init_b = [20,12,-10]
Loes_1_eps = solve_ivp(DGLsys, [0, t_end], vec_init_eps, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_a = solve_ivp(DGLsys, [0, t_end], vec_init_a, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=0.1, c = Loes_1.t, cmap=cm.gnuplot2)
ax.scatter3D(Loes_1_eps.y[0],Loes_1_eps.y[1],Loes_1_eps.y[2], marker='o', s=0.1, c = Loes_1_eps.t, cmap=cm.gnuplot2)
ax.scatter3D(Loes_1_a.y[0],Loes_1_a.y[1],Loes_1_a.y[2], marker='o', s=0.1, c = Loes_1_a.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=22 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir wählen jedoch noch die Anfangsbedingung $x(0) = \alpha_1 = 20 \, , \,\, y(0) = \alpha_2 = 12$ und $z(0) = \alpha_3 = -10$ so entsteht ein vorübergehendes chaotisches Verhalten ('transient chaos'), wie wir im Folgenden zeigen werden. Wir lösen das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 20 \, , \,\, y(0) = \alpha_2 = 12$ und $z(0) = \alpha_3 = -10$ im Bereich $t \in [0,70]$ unter Verwendung von $N=100000$ zeitlichen Gitterpunkten (mesh points) für $r = 22$. In den folgenden beiden Simulationen setzen wir die relativen und absoluten Fehler-Toleranzen (rtol und atol) auf zwei unterschiedliche Werte: $\Delta_0=10^{-10}$ und $\Delta_1=10^{-13}$. Zusätzlich zu dieser Simulation machen wir auch eine Simulation, bei der die Anfangsbedingung um einen Epsilon-Wert verschoben ist (hier speziell $[x(0),y(0),z(0)] = [20,11.999,-10]$).
vec_init_b = [20,12,-10]
vec_init_b_eps = [20,11.999,-10]
t_end_b = 70
N = 100000
t_val_b = np.linspace(0, t_end_b, N+1)
Loes_0_b = solve_ivp(DGLsys, [0, t_end_b], vec_init_b, t_eval=t_val_b, args=(set_r, ), rtol=fehler_0, atol=fehler_0)
Loes_1_b = solve_ivp(DGLsys, [0, t_end_b], vec_init_b, t_eval=t_val_b, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_b_eps = solve_ivp(DGLsys, [0, t_end_b], vec_init_b_eps, t_eval=t_val_b, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
plt.title(r"$\rm Python-Lösung \; mit \; r=22 \quad und \quad \Delta_1=10^{-13}$")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(Loes_1_b.t, Loes_1_b.y[0],c="blue", label="x");
plt.plot(Loes_1_b.t, Loes_1_b.y[1],c="green", label="y");
plt.plot(Loes_1_b.t, Loes_1_b.y[2],c="red", label="z");
plt.legend(loc='upper right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen mit unterschiedlichen Fehler-Toleranzen dar (z.B. $\Delta x = x_{\Delta_0} - x_{\Delta_1}$).
plt.title(r"$\rm Python-Lösung \; mit \; r=22 \quad , \quad Abweichungen: \, \Delta \, = \, \Delta_0 \, - \, \Delta_1$")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_0_b.t, Loes_1_b.y[0]-Loes_0_b.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_0_b.t, Loes_1_b.y[1]-Loes_0_b.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_0_b.t, Loes_1_b.y[2]-Loes_0_b.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Simulation zur Epsilon-verschobenen Anfangsbedingung ergibt sich.
plt.title(r"Python-Lösung mit r=22 , Abweichungen zu einer $\epsilon$ verschobenen Anfangsbedingung")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_0_b.t, Loes_1_b.y[0]-Loes_1_b_eps.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_0_b.t, Loes_1_b.y[1]-Loes_1_b_eps.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_0_b.t, Loes_1_b.y[2]-Loes_1_b_eps.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
In der oberen Abbildung erkennt man, dass der zeitliche Abschnitt des vorübergehenden chaotischen Verhaltens dazu geführt hat, dass sich die Simulation der $\epsilon$-verschobene Anfangsbedingung zu einem anderen Fixpunkt entwickelt als die ursprüngliche. Dies können wir auch in den folgenden Trajektorien gut sehen.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1_b.y[0],Loes_1_b.y[1],Loes_1_b.y[2], marker='o', s=0.1, c = Loes_1_b.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=22 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1_b_eps.y[0],Loes_1_b_eps.y[1],Loes_1_b_eps.y[2], marker='o', s=0.1, c = Loes_1_b_eps.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=22 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13} \quad , \quad Anfangsbed. \, \epsilon-verschoben$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Diese Eigenschaft, dass kleine Veränderungen der Anfangsbedingungen die gesamte Entwicklung des Systems nach einiger Zeit verändern können, bezeichnet man als der "Schmetterlingseffekt" in der Chaostheorie.
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 20 \, , \,\, y(0) = \alpha_2 = 12$ und $z(0) = \alpha_3 = -10$ im Bereich $t \in [0,70]$ unter Verwendung von $N=100000$ zeitlichen Gitterpunkten (mesh points) mittels des C++ Programms Lorenz.cpp. Zusätzlich zu dieser Simulation machen wir auch eine Simulation, bei der die Anfangsbedingung um einen $\epsilon$-Wert verschoben ist (hier speziell $[x(0),y(0),z(0)] = [20,11.999,-10]$). Die Ergebnisse des C++ Programms wurden in die Dateien "Lorenz-r_22-init_b.dat" und "Lorenz-r_22-init_b_eps.dat" geschrieben. Im Folgenden stellen wir nur die Ergebnisse der Runge Kutta RK4 Methode dar.
data = np.genfromtxt("./Lorenz-r_22-init_b.dat")
data_eps = np.genfromtxt("./Lorenz-r_22-init_b_eps.dat")
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=22 , N=100000")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(data[:,1],data[:,5],c="blue", label=r"$\rm x$");
plt.plot(data[:,1],data[:,6],c="green", label=r"$\rm y$");
plt.plot(data[:,1],data[:,7],c="red", label=r"$\rm z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Simulation zur Epsilon-verschobenen Anfangsbedingung ergibt sich.
plt.title(r"C++ Lösung (Abweichungen $\Delta$ zur $\epsilon$-verschobenen Anfangsb: r=22 , N=100000")
plt.xlabel("t")
plt.ylabel(r"$\rm\Delta x, \Delta y, \Delta z$")
plt.plot(data[:,1],data[:,5]-data_eps[:,5],c="blue", label=r"$\rm \Delta x$");
plt.plot(data[:,1],data[:,6]-data_eps[:,6],c="green", label=r"$\rm \Delta y$");
plt.plot(data[:,1],data[:,7]-data_eps[:,7],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Runge Kutta RK4 Methode im Vergleich zur Python Lösung ($\Delta_1=10^{-13}$) ergibt sich:
plt.title(r"Abweichungen $\Delta$ Runge Kutta RK4 vs. Phython $\Delta_1=10^{-13}$): r=22 , N=100000")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm y_1,y_2,y_3$")
plt.plot(data[:,1],Loes_1_b.y[0]-data[:,5],c="blue", label=r"$\rm y_1$");
plt.plot(data[:,1],Loes_1_b.y[1]-data[:,6],c="green", label=r"$\rm y_2$");
plt.plot(data[:,1],Loes_1_b.y[2]-data[:,7],c="red", label=r"$\rm y_3$");
plt.legend(loc='lower right',fontsize=16);
Stellt man sich ausgehend von den simulierten Daten der Runge Kutta RK4 Methode die Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ in einem dreidimensionalen $x,y,z$-Koordinatensystem dar, so ergibt sich wieder das folgende Bild.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(data[:,5],data[:,6],data[:,7], marker='o', s=0.1, c = data[:,1], cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=22 , N=100000")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(data_eps[:,5],data_eps[:,6],data_eps[:,7], marker='o', s=0.1, c = data_eps[:,1], cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"C++ Lösung (Runge Kutta RK4 Methode): r=22 , N=100000, Anfangsbed. $\epsilon$-verschoben")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,50]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points) für $r = 28$.
t_end = 50
N = 100000
t_val = np.linspace(0, t_end, N+1)
set_r = 28
Loes_0 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_0, atol=fehler_0)
Loes_1 = solve_ivp(DGLsys, [0, t_end], vec_init, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Wir stellen uns die numerische Lösungen grafisch dar:
plt.title(r"$\rm r=28 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(Loes_1.t, Loes_1.y[0],c="blue", label="x");
plt.plot(Loes_1.t, Loes_1.y[1],c="green", label="y");
plt.plot(Loes_1.t, Loes_1.y[2],c="red", label="z");
plt.legend(loc='upper right',fontsize=16);
Wir stellen uns die Abweichungen $\Delta$ der beiden Simulationen dar (z.B. $\Delta x = x_{\Delta_0} - x_{\Delta_1}$).
plt.title(r"$\rm r=28 \quad , \quad Abweichungen: \, \Delta_1=10^{-13} \, - \, \Delta_0=10^{-10}$")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_0.t, Loes_1.y[0]-Loes_0.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_0.t, Loes_1.y[1]-Loes_0.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_0.t, Loes_1.y[2]-Loes_0.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Die beiden Lösungen stimmen für $t>30$ nichtmehr miteinander überein. Die zeitliche Entwicklung zeigt somit bei $r=28$ ein chaotisches, nicht mehr vorhersagbares Verhalten.
Die dreidimensionale Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ sieht wie folgt aus.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1.y[0],Loes_1.y[1],Loes_1.y[2], marker='o', s=0.1, c = Loes_1.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=28 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Die obere Abbildung zeigt das chaotische Verhalten des "seltsamen Lorenz-Attraktors". Betrachtet man sich die zeitliche Entwicklung ausgehend von anderen Anfangswerten (hier speziell $[x(0),y(0),z(0)] = [0,30,45]$ und $[20,12,-10]$), so erhält man für $r=28$ stets diesen "seltsamen Attraktor" (siehe folgende Abbildungen). Verändert man den Anfangswert nur um einen Epsilon-Wert (hier speziell $[x(0),y(0),z(0)] = [0,0.999,0]$), so kann man die zeitliche Entwicklung der Trajektorie nach einer gewissen Zeit nicht mehr vorhersagen.
vec_init_a = [0,30,45]
vec_init_b = [20,12,-10]
Loes_1_eps = solve_ivp(DGLsys, [0, t_end], vec_init_eps, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_a = solve_ivp(DGLsys, [0, t_end], vec_init_a, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Loes_1_b = solve_ivp(DGLsys, [0, t_end], vec_init_b, t_eval=t_val, args=(set_r, ), rtol=fehler_1, atol=fehler_1)
Für die Abweichungen $\Delta$ der Simulation zur Epsilon-verschobenen Anfangsbedingung ergibt sich.
plt.title(r"Python-Lösung mit r=28 , Abweichungen zu einer $\epsilon$ verschobenen Anfangsbedingung")
plt.xlabel("t")
plt.ylabel(r"$\rm \Delta x, \Delta y, \Delta z$")
plt.plot(Loes_1.t, Loes_1.y[0]-Loes_1_eps.y[0],c="blue", label=r"$\rm \Delta x$");
plt.plot(Loes_1.t, Loes_1.y[1]-Loes_1_eps.y[1],c="green", label=r"$\rm \Delta y$");
plt.plot(Loes_1.t, Loes_1.y[2]-Loes_1_eps.y[2],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Für den $\epsilon$ verschobenen Anfangswert sieht der Lorenz-Attraktor wie folgt aus:
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1_eps.y[0],Loes_1_eps.y[1],Loes_1_eps.y[2], marker='o', s=0.1, c = Loes_1_eps.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=28 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Für den Anfangswert $[x(0),y(0),z(0)] = [0,30,45]$ sieht der Lorenz-Attraktor wie folgt aus:
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1_a.y[0],Loes_1_a.y[1],Loes_1_a.y[2], marker='o', s=0.1, c = Loes_1_a.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=28 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Für den Anfangswert $[x(0),y(0),z(0)] = [20,12,-10]$ sieht der Lorenz-Attraktor wie folgt aus:
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(Loes_1_b.y[0],Loes_1_b.y[1],Loes_1_b.y[2], marker='o', s=0.1, c = Loes_1_b.t, cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"$\rm r=28 \quad , \quad Fehler-Toleranzen: \, \Delta_1=10^{-13}$")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
Wir lösen nun das Anfangswertproblem mit der Anfangsbedingung $x(0) = \alpha_1 = 0 \, , \,\, y(0) = \alpha_2 = 1$ und $z(0) = \alpha_3 = 0$ im Bereich $t \in [0,50]$ unter Verwendung von $N=100000$ zeitlichen Gitterpunkten (mesh points) mittels des C++ Programms Lorenz.cpp. Zusätzlich zu dieser Simulation machen wir auch eine Simulation, bei der die Anfangsbedingung um einen $\epsilon$-Wert verschoben ist (hier speziell $[x(0),y(0),z(0)] = [0,0.999,0]$). Die Ergebnisse des C++ Programms wurden in die Dateien "Lorenz-r_28-init_0.dat" und "Lorenz-r_28-init_0_eps.dat" geschrieben. Im Folgenden stellen wir nur die Ergebnisse der Runge Kutta RK4 Methode dar.
data = np.genfromtxt("./Lorenz-r_28-init_0.dat")
data_eps = np.genfromtxt("./Lorenz-r_28-init_0_eps.dat")
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=28 , N=100000")
plt.xlabel("t")
plt.ylabel("x,y,z")
plt.plot(data[:,1],data[:,5],c="blue", label=r"$\rm x$");
plt.plot(data[:,1],data[:,6],c="green", label=r"$\rm y$");
plt.plot(data[:,1],data[:,7],c="red", label=r"$\rm z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Simulation zur Epsilon-verschobenen Anfangsbedingung ergibt sich.
plt.title(r"C++ Lösung (Abweichungen $\Delta$ zur $\epsilon$-verschobenen Anfangsb: r=28 , N=100000")
plt.xlabel("t")
plt.ylabel(r"$\rm\Delta x, \Delta y, \Delta z$")
plt.plot(data[:,1],data[:,5]-data_eps[:,5],c="blue", label=r"$\rm \Delta x$");
plt.plot(data[:,1],data[:,6]-data_eps[:,6],c="green", label=r"$\rm \Delta y$");
plt.plot(data[:,1],data[:,7]-data_eps[:,7],c="red", label=r"$\rm \Delta z$");
plt.legend(loc='lower right',fontsize=16);
Für die Abweichungen $\Delta$ der Runge Kutta RK4 Methode im Vergleich zur Python Lösung ($\Delta_1=10^{-13}$) ergibt sich:
plt.title(r"Abweichungen $\Delta$ Runge Kutta RK4 vs. Phython $\Delta_1=10^{-13}$): r=28 , N=100000")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm y_1,y_2,y_3$")
plt.plot(data[:,1],Loes_1.y[0]-data[:,5],c="blue", label=r"$\rm y_1$");
plt.plot(data[:,1],Loes_1.y[1]-data[:,6],c="green", label=r"$\rm y_2$");
plt.plot(data[:,1],Loes_1.y[2]-data[:,7],c="red", label=r"$\rm y_3$");
plt.legend(loc='lower right',fontsize=16);
Stellt man sich ausgehend von den simulierten Daten der Runge Kutta RK4 Methode die Trajektorie der Bewegung der Lösungsfunktionen $x(t), y(t)$ und $z(t)$ in einem drei-dimensionalen $x,y,z$-Koordinatensystem dar, so ergibt sich wieder das folgende Bild des seltsamen Lorenz-Attrakors.
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(data[:,5],data[:,6],data[:,7], marker='o', s=0.1, c = data[:,1], cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title("C++ Lösung (Runge Kutta RK4 Methode): r=28 , N=100000")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
bild = ax.scatter3D(data_eps[:,5],data_eps[:,6],data_eps[:,7], marker='o', s=0.1, c = data_eps[:,1], cmap=cm.gnuplot2)
ax.view_init(azim=-75, elev=10)
cbar = fig.colorbar(bild, shrink=0.6, aspect=15)
cbar.set_label(r"$\rm t$", rotation=270)
plt.title(r"C++ Lösung (Runge Kutta RK4 Methode): r=28 , N=100000, Anfangsbed. $\epsilon$-verschoben")
ax.set_xlabel(r"$\rm x(t)$")
ax.set_ylabel(r"$\rm y(t)$")
ax.set_zlabel(r"$\rm z(t)$");
