Einführung in die Programmierung für Studierende der Physik

(Introduction to Programming for Physicists)

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main

(Sommersemester 2025)

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske

Frankfurt am Main 20.03.2025

Das Doppelpendel

Das Doppelpendel ist ein Anwendungsfall aus der klassischen Mechanik (siehe Walter Greiner, 'Klassische Mechanik II' [8. Auflage, 2008, Kapitel V15. Seite 269]). Das System besteht aus zwei miteinander verbundenen Pendeln, wobei wir die Luftreibung zunächst vernachlässigen. Die Herleitung der Bewegungsgleichungen des Doppelpendels erfolgt am elegantesten mittels der Euler-Lagrange Gleichungen, bzw. mittels der Hamilton Theorie. In der Euler-Lagrange Theorie beschreibt man das Doppelpendel mittels zweier generalisierten Koordinaten $\theta_1(t)$ und $\theta_2(t)$, welche die Ausschläge der beiden Pendel als Pendelwinkel zur Vertikalen beschreiben. Die generalisierten Geschwindigkeiten $\dot{\theta_1}(t)$ und $\dot{\theta_2}(t)$ stellen die Winkelgeschwindigkeiten der beiden Pendel dar. Die Lagrangefunktion ergibt sich, indem man die potentielle Energie $V$ des Systems von seiner kinetischen Energie $T$ subtrahiert ($L = T - V$). Die potentielle Energie definieren wir dabei relativ zu einer Ebene im Abstand $l_1 + l_2$ unterhalb des Aufhängepunktes des Doppelpendels. Betrachtet man die Bewegung des Doppelpendels in der x-y-Ebene, so ergibt sich die folgende Lagrangefunktion des Systems

$$ \begin{eqnarray} L = T - V &=& \frac{1}{2} m_1 \left[ \dot{x_1}(t)^2 + \dot{y_1}(t)^2 \right] + \frac{1}{2} m_2 \left[ \dot{x_2}(t)^2 + \dot{y_2}(t)^2 \right] - V \\ \hbox{mit:} && x_1 = l_1 \, \hbox{sin}(\theta_1) \, , \,\, y_1 = - l_1 \, \hbox{cos}(\theta_1) \quad , \,\, x_2 = l_1 \, \hbox{sin}(\theta_1) + l_2 \, \hbox{sin}(\theta_2) \, , \,\, y_2 = - l_1 \, \hbox{cos}(\theta_1) - l_2 \, \hbox{cos}(\theta_2) \\ && V = V(\theta_1,\theta_2) = m_1 \, g \left[ l_1 - l_1 \, \hbox{cos}(\theta_1) \right] + m_2 \, g \left[ l_1 + l_2 - \left( l_1 \, \hbox{cos}(\theta_1) + l_2 \, \hbox{cos}(\theta_2) \right) \right] \quad , \end{eqnarray} $$

wobei $m_1$ und $m_2$ die Massen und $l_1$ und $l_2$ die Längen der beiden Pendel darstellen. </p>

Mittels der Lagrange-Gleichungen

$$ \begin{equation} \frac{d }{dt} \left( \frac{\partial L}{\partial \dot{\theta}_1} \right) - \frac{\partial L}{\partial \theta_1} \, =\, 0 \quad \hbox{und} \quad \frac{d }{dt} \left( \frac{\partial L}{\partial \dot{\theta}_2} \right) - \frac{\partial L}{\partial \theta_2} \, =\, 0 \end{equation} $$

gelangt man zu zwei gekoppelten Differentialgleichungen zweiter Ordnung, die man dann in ein System von vier gekoppelten DGLs erster Ordung umschreiben muss um sie numerisch lösen zu können.

Von den Euler-Lagrange Gleichungen zur Bewegungsgleichung

In diesem Jupyter Notebook berechnen wir zunächst die Bewegungsgleichungen des Doppelpendels auf analytischem Weg, indem wir unter Verwendung der SymPy Bibliothek die Euler-Lagrange Gleichungen des Doppelpendels bestimmen. Dann schreiben wir die Bewegungsgleichung in ein System von erster Ordnung Differentialgleichungen um.

from sympy import *
init_printing()

Wir definieren die nötigen sympy Symbole und Funktionen um die analytischen Ausdrücke der kinetischen Energie $T$, potentiellen Energie $V$ und Gesamtenergie $E$ festzulegen. Dann definieren wir die Lagrangefunktion $L$ und geben diese aus.

t = symbols('t',real=True)
m_1, m_2, l_1, l_2, g = symbols('m_1, m_2, l_1, l_2, g', positive = True, real = True)
theta_1 = Function('theta_1')(t)
theta_2 = Function('theta_2')(t)

x_1 = l_1 * sin(theta_1)
y_1 = - l_1 * cos(theta_1)
x_2 = l_1 * sin(theta_1) + l_2 * sin(theta_2)
y_2 = - l_1 * cos(theta_1) - l_2 * cos(theta_2)

T = 1/2 * m_1 * (diff(x_1,t)**2 + diff(y_1,t)**2) + 1/2 * m_2 * (diff(x_2,t)**2 + diff(y_2,t)**2)
#V = m_1 * g * (l_1 + l_2 - l_1 * cos(theta_1)) + m_2 * g * (l_1 + l_2 - (l_1 * cos(theta_1) + l_2 * cos(theta_2)))
V = m_1 * g * l_1 *(1 - cos(theta_1)) + m_2 * g * (l_1 + l_2 - (l_1 * cos(theta_1) + l_2 * cos(theta_2)))
L = T - V
E = T + V
simplify(L)

$\displaystyle g l_{1} m_{1} \left(\cos{\left(\theta_{1}{\left(t \right)} \right)} - 1\right) + g m_{2} \left(l_{1} \cos{\left(\theta_{1}{\left(t \right)} \right)} - l_{1} + l_{2} \cos{\left(\theta_{2}{\left(t \right)} \right)} - l_{2}\right) + 0.5 l_{1}^{2} m_{1} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 0.5 m_{2} \left(l_{1}^{2} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 2 l_{1} l_{2} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d}{d t} \theta_{1}{\left(t \right)} \frac{d}{d t} \theta_{2}{\left(t \right)} + l_{2}^{2} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}\right)$

Mittels der so definierten Lagrangefunktion $L$ und den Euler-Lagrange Gleichungen bestimmen wir die Bewegungsgleichungen des Doppelpendels. Die beiden Euler-Lagrange Gleichungen lauten wie folgt:

ELG_1 = diff(L.diff(theta_1.diff(t)),t) - L.diff(theta_1)
Eq(0,simplify(ELG_1))

$\displaystyle 0 = l_{1} \left(g m_{1} \sin{\left(\theta_{1}{\left(t \right)} \right)} + g m_{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} + l_{1} m_{1} \frac{d^{2}}{d t^{2}} \theta_{1}{\left(t \right)} + l_{1} m_{2} \frac{d^{2}}{d t^{2}} \theta_{1}{\left(t \right)} + l_{2} m_{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + l_{2} m_{2} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \theta_{2}{\left(t \right)}\right)$

ELG_2 = diff(L.diff(theta_2.diff(t)),t) - L.diff(theta_2)
Eq(0,simplify(ELG_2))

$\displaystyle 0 = l_{2} m_{2} \left(g \sin{\left(\theta_{2}{\left(t \right)} \right)} - l_{1} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + l_{1} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \theta_{1}{\left(t \right)} + l_{2} \frac{d^{2}}{d t^{2}} \theta_{2}{\left(t \right)}\right)$

Dies ist ein gekoppeltes System von zwei Differentialgleichungen zweiter Ordnung in der Zeit. Wir formen die erste Gleichung nach $\frac{d^2 \theta_1}{dt^2}$ und die zweite Gleichung nach $\frac{d^2 \theta_2}{dt^2}$ um.

ELG_1a = Eq(diff(theta_1,t,t),solve(ELG_1,diff(theta_1,t,t))[0].simplify())
ELG_1a

$\displaystyle \frac{d^{2}}{d t^{2}} \theta_{1}{\left(t \right)} = - \frac{g m_{1} \sin{\left(\theta_{1}{\left(t \right)} \right)} + g m_{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} + l_{2} m_{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + l_{2} m_{2} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \theta_{2}{\left(t \right)}}{l_{1} \left(m_{1} + m_{2}\right)}$

ELG_2a = Eq(diff(theta_2,t,t),solve(ELG_2,diff(theta_2,t,t))[0].simplify())
ELG_2a

$\displaystyle \frac{d^{2}}{d t^{2}} \theta_{2}{\left(t \right)} = \frac{- g \sin{\left(\theta_{2}{\left(t \right)} \right)} + l_{1} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - l_{1} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \theta_{1}{\left(t \right)}}{l_{2}}$

Wir substituieren nun in der ersten Gleichung den Ausdruck von $\frac{d^2 \theta_2}{dt^2}$ der zweiten Gleichung und formen wieder nach $\frac{d^2 \theta_1}{dt^2}$ um.

ELG_1b = Eq(diff(theta_1,t,t),solve(ELG_1a.subs(diff(theta_2,t,t),ELG_2a.rhs),diff(theta_1,t,t))[0].simplify())
ELG_1b

$\displaystyle \frac{d^{2}}{d t^{2}} \theta_{1}{\left(t \right)} = \frac{g m_{1} \sin{\left(\theta_{1}{\left(t \right)} \right)} + \frac{g m_{2} \sin{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)}}{2} + \frac{g m_{2} \sin{\left(\theta_{1}{\left(t \right)} \right)}}{2} + \frac{l_{1} m_{2} \sin{\left(2 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2}}{2} + l_{2} m_{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}}{l_{1} \left(- m_{1} + m_{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - m_{2}\right)}$

In gleicher Weise substituieren wir auch in der zweiten Gleichung den Ausdruck von $\frac{d^2 \theta_1}{dt^2}$ der ersten Gleichung und formen wieder nach $\frac{d^2 \theta_2}{dt^2}$ um.

ELG_2b = Eq(diff(theta_2,t,t),solve(ELG_2a.subs(diff(theta_1,t,t),ELG_1a.rhs),diff(theta_2,t,t))[0].simplify())
ELG_2b

$\displaystyle \frac{d^{2}}{d t^{2}} \theta_{2}{\left(t \right)} = \frac{g m_{1} \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - g m_{1} \sin{\left(\theta_{2}{\left(t \right)} \right)} + g m_{2} \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - g m_{2} \sin{\left(\theta_{2}{\left(t \right)} \right)} + 2 l_{1} m_{1} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 2 l_{1} m_{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + l_{2} m_{2} \sin{\left(2 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}}{2 l_{2} \left(m_{1} - m_{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + m_{2}\right)}$

Die beiden oberen Gleichungen stellen die Bewegungsgleichungen des Doppelpendels dar.

Umschreiben in ein System von Differentialgleichungen erster Ordnung

Um das gekoppelte System von DGLs numerisch lösen zu können, schreiben wir es in ein System von DGLs erster Ordung um. Wir betrachten nun die numerische Lösung und machen die folgenden Variablenumbenennung: $u_1(t)=\theta_1(t)$, $u_2(t)=\theta_2(t)$, $u_3(t)=\dot{\theta}_1(t)$, $u_4(t)=\dot{\theta}_2(t)$ und definieren das System von DGLs wie folgt:

$$ \begin{eqnarray} \dot{u_1}(t) &=& \frac{d u_1}{dt} = \frac{d \theta_1}{dt} = u_3(t) \\ \dot{u_2}(t) &=& \frac{d u_2}{dt} = \frac{d \theta_2}{dt} = u_4(t) \\ \dot{u_3}(t) &=& \frac{d u_3}{dt} = \frac{d \dot{\theta_1}}{dt} = \ddot{\theta}_1(t) := f_1(u_1,u_2,u_3,u_4)\\ \dot{u_4}(t) &=& \frac{d u_4}{dt} = \frac{d \dot{\theta_2}}{dt} = \ddot{\theta}_2(t) := f_2(u_1,u_2,u_3,u_4) \end{eqnarray} $$

u_1 = Function('u_1')(t)
u_2 = Function('u_2')(t)
u_3 = Function('u_3')(t)
u_4 = Function('u_4')(t)
f_1 = ELG_1b.rhs.subs(diff(theta_1,t),u_3).subs(diff(theta_2,t),u_4).subs(theta_1,u_1).subs(theta_2,u_2)
f_2 = ELG_2b.rhs.subs(diff(theta_1,t),u_3).subs(diff(theta_2,t),u_4).subs(theta_1,u_1).subs(theta_2,u_2)

Eq(diff(u_3,t),f_1)

$\displaystyle \frac{d}{d t} u_{3}{\left(t \right)} = \frac{g m_{1} \sin{\left(u_{1}{\left(t \right)} \right)} + \frac{g m_{2} \sin{\left(u_{1}{\left(t \right)} - 2 u_{2}{\left(t \right)} \right)}}{2} + \frac{g m_{2} \sin{\left(u_{1}{\left(t \right)} \right)}}{2} + \frac{l_{1} m_{2} u_{3}^{2}{\left(t \right)} \sin{\left(2 u_{1}{\left(t \right)} - 2 u_{2}{\left(t \right)} \right)}}{2} + l_{2} m_{2} u_{4}^{2}{\left(t \right)} \sin{\left(u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)}}{l_{1} \left(- m_{1} + m_{2} \cos^{2}{\left(u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)} - m_{2}\right)}$

Eq(diff(u_4,t),f_2)

$\displaystyle \frac{d}{d t} u_{4}{\left(t \right)} = \frac{g m_{1} \sin{\left(2 u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)} - g m_{1} \sin{\left(u_{2}{\left(t \right)} \right)} + g m_{2} \sin{\left(2 u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)} - g m_{2} \sin{\left(u_{2}{\left(t \right)} \right)} + 2 l_{1} m_{1} u_{3}^{2}{\left(t \right)} \sin{\left(u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)} + 2 l_{1} m_{2} u_{3}^{2}{\left(t \right)} \sin{\left(u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)} + l_{2} m_{2} u_{4}^{2}{\left(t \right)} \sin{\left(2 u_{1}{\left(t \right)} - 2 u_{2}{\left(t \right)} \right)}}{2 l_{2} \left(m_{1} - m_{2} \cos^{2}{\left(u_{1}{\left(t \right)} - u_{2}{\left(t \right)} \right)} + m_{2}\right)}$

Um die mittels sympy gewonnenen analytische Ausdrücke direkt beim numerischen Lösen zu verwenden, schreiben wir die Funktionen $f_1(u_1,u_2,u_3,u_4)$ und $f_2(u_1,u_2,u_3,u_4)$ in numerische Lambda-Funktionen um. Zusätzlich definieren wir auch eine Lambda-Funktion der kinetischen Energie, potentiellen Energie und Gesamtenergie. Auch legen wir ab jetzt den Wert der Erdbeschleunigung fest $g = 9.81$.

set_g = 9.81
f_1L = lambdify((u_1,u_2,u_3,u_4,m_1,m_2,l_1,l_2), f_1.subs(g,set_g))
f_2L = lambdify((u_1,u_2,u_3,u_4,m_1,m_2,l_1,l_2), f_2.subs(g,set_g))
E_L = lambdify((u_1,u_2,u_3,u_4,m_1,m_2,l_1,l_2), E.subs(diff(theta_1,t),u_3).subs(diff(theta_2,t),u_4).subs(theta_1,u_1).subs(theta_2,u_2).subs(g,set_g))
T_L = lambdify((u_1,u_2,u_3,u_4,m_1,m_2,l_1,l_2), T.subs(diff(theta_1,t),u_3).subs(diff(theta_2,t),u_4).subs(theta_1,u_1).subs(theta_2,u_2).subs(g,set_g))
V_L = lambdify((u_1,u_2,u_3,u_4,m_1,m_2,l_1,l_2), V.subs(diff(theta_1,t),u_3).subs(diff(theta_2,t),u_4).subs(theta_1,u_1).subs(theta_2,u_2).subs(g,set_g))

Wir definieren nun das System der gekoppelten Differentialgleichung des Doppelpendels. Die Werte der Massenparameter $m_1$ und $m_2$ setzen wir auf Eins und die Längen der Pendel $l_1$ und $l_2$ legen wir noch nicht fest; diese werden als zusätzliche Argumente in der DGL-Funktion definiert.

set_m_1 = 1
set_m_2 = 1

def DGLsys(t,u_vec,l_1,l_2):
    u1, u2, u3, u4 = u_vec
    du1_dt = u3
    du2_dt = u4
    du3_dt = f_1L(u1,u2,u3,u4,set_m_1,set_m_2,l_1,l_2)
    du4_dt = f_2L(u1,u2,u3,u4,set_m_1,set_m_2,l_1,l_2)
    return np.array([du1_dt,du2_dt,du3_dt,du4_dt])

Alternativ kann man auch die gewonnenen Ausdrücke für $f_1(u_1,u_2,u_3,u_4)$ und $f_2(u_1,u_2,u_3,u_4)$ explizit in der DGL-Funktion angeben. Diese Gleichungen werden wir auch in dem C++ Programm implementieren, dessen Ergebnisse wir uns auch später ansehen werden.

import numpy as np

def DGLsysa(t,u_vec,l_1,l_2):
    u_1, u_2, u_3, u_4 = u_vec
    g = set_g
    m_1 = set_m_1
    m_2 = set_m_2
    du1_dt = u_3
    du2_dt = u_4
    du3_dt = (g*m_1*np.sin(u_1) + g*m_2*np.sin(u_1)/2 + g*m_2*np.sin(u_1 - 2*u_2)/2 + l_1*m_2*u_3**2*np.sin(2*u_1 - 2*u_2)/2 + l_2*m_2*u_4**2*np.sin(u_1 - u_2))/(l_1*(-m_1 + m_2*np.cos(u_1 - u_2)**2 - m_2))
    du4_dt = (-g*m_1*np.sin(u_2) + g*m_1*np.sin(2*u_1 - u_2) - g*m_2*np.sin(u_2) + g*m_2*np.sin(2*u_1 - u_2) + 2*l_1*m_1*u_3**2*np.sin(u_1 - u_2) + 2*l_1*m_2*u_3**2*np.sin(u_1 - u_2) + l_2*m_2*u_4**2*np.sin(2*u_1 - 2*u_2))/(2*l_2*(m_1 - m_2*np.cos(u_1 - u_2)**2 + m_2))
    return np.array([du1_dt,du2_dt,du3_dt,du4_dt])

Numerische Lösung des Problems

Zum Lösen des gekoppelten Systems von Differenzialgleichungen benutzen wir das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Semester die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt. (Bemerkung: In den vergangenen Semestern benutzten wir noch die Funktion 'solve_odeint(...)', welche auf einem älteren Fortran-basierten Solver basiert.)

import matplotlib.pyplot as plt 
import matplotlib
from scipy.integrate import solve_ivp

Bei dem numerischen Lösen von Differentialgleichungen in Python (mittels solve_ivp(...)) kann der Benutzer zusätzlich die relativen und absoluten Fehler-Toleranzen der berechneten numerischen Werte festlegen, wobei der relative Fehler mit der Zusatzoption 'rtol' und der absolute Fehler mit der Zusatzoption 'atol' kontrolliert wird (näheres siehe scipy.integrate.solve_ivp). In Folgenden setzen wir die relativen und absoluten Fehler-Toleranzen (rtol und atol) auf $10^{-13}$.

Die Bewegungsgleichung des Doppelpendels ist nicht-linear und, abhängig von den Parametern $m_1$, $m_2$, $l_1$, $l_2$ und Anfangsbedingungen können sehr unterschiedliche und teilweise chaotische Bewegungsformen auftreten. Im Folgenden werden wir nur einen kleinen Ausschnitt aus dem Parameterraum untersuchen und die Abhängigkeit von der Anfangsbedingungen und Pendellängen $l_1$ und $l_2$ bei festgehaltenen Masse $m_1=1$ und $m_2=1$ betrachten.

Reguläre Bewegungen bei geringer Energie

Zunächst legen wir die Längen des Doppelpendels fest (gleiche Längen der Pendel $l_1=1$ und $l_2=1$). Als Anfangswert zur Zeit $t=0$ betrachten wir das Pendel als ruhend, wobei wir die beiden Massen nur geringfügig aus der Nulllage auslenken:

$\theta_1(0) = \pi/14 \,\, , \,\, \theta_2(0) = \sqrt{2} \cdot \theta_1(0) \,\, , \,\, \dot{\theta}_1(0) = 0 \,\, , \,\, \dot{\theta}_2(0) = 0 \longrightarrow E = V = 0.9817662576217473 \, , \quad T=0$

Wir lösen das Anfangswertproblem im Bereich $t \in [0,15]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points).

t_end = 15
N = 10000
t_val = np.linspace(0, t_end, N)
set_l_1 = 1
set_l_2 = 1
theta1_0 = np.pi/14
theta2_0 = np.sqrt(2)*theta1_0
dtheta1_0 = 0
dtheta2_0 = 0
u_init_a = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]
fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init_a, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  0.9817662576217473
Potentielle Energie V=  0.9817662576217473
Kinetische Energie T =  0.0

Wir veranschaulichen uns nun die numerischen Ergebnisse nun in unterschiedlichen Diagrammen.

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.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \theta_1, \, \theta_2$")
plt.plot(Loes.t, Loes.y[0],c="blue", label=r"$\rm \theta_1(t)$")
plt.plot(Loes.t, Loes.y[1],c="red", label=r"$\rm \theta_2(t)$")
plt.legend(loc='lower right',fontsize=16);

Die obere Abbildung zeigt, wie sich die Winkel $\theta_1$ und $\theta_2$ des Doppelpendels im Laufe der Zeit verändern und die untere Abbildung zeigt das zeitliche Verhalten der Winkelgeschwindigkeiten $\frac{d \theta_1(t)}{dt}$ und $\frac{d \theta_2(t)}{dt}$. Man erkennt eine sehr gleichförmige Bewegung der beiden Pendel.

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \dot{\theta_1}, \, \dot{\theta_2}$")
plt.plot(Loes.t, Loes.y[2],c="blue", label=r"$\rm \dot{\theta}_1(t)$")
plt.plot(Loes.t, Loes.y[3],c="red", label=r"$\rm \dot{\theta}_2(t)$")
plt.legend(loc='lower right',fontsize=16);

In der unteren Abbildung erkennt man, dass sich die potentielle und kinetische Energie periodisch austauschen und die Gesamtenergie des Systems erhalten bleibt.

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm Energie \, E \, , \, \, T  \, , \, \, V$")
plt.plot(Loes.t, E_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="black", label=r"$\rm Gesamtenergie \, E$");
plt.plot(Loes.t, T_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="green", label=r"$\rm kinetische \, Energie \, T$");
plt.plot(Loes.t, V_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="orange", label=r"$\rm potentielle \, Energie \, V$");
plt.legend(loc='upper right',fontsize=16);

Die Bewegung im Konfigurationsraum $(x,y)$ sieht für die beiden Pendel wie folgt aus.

import matplotlib.gridspec as gridspec

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Pendel 1")
ax2.set_title("Pendel 2")
ax1.set_xlabel(r"$\rm x_1$")
ax1.set_ylabel(r"$\rm y_1$")
ax2.set_xlabel(r"$\rm x_2$")
ax2.set_ylabel(r"$\rm y_2$")

ax1.plot(set_l_1 * np.sin(Loes.y[0]), - set_l_1 * np.cos(Loes.y[0]),c="blue")
ax2.plot(set_l_1 * np.sin(Loes.y[0]) + set_l_2 * np.sin(Loes.y[1]), - set_l_1 * np.cos(Loes.y[0]) - set_l_2 * np.cos(Loes.y[1]),c="red");

Die untere Abbildung veranschaulicht, dass die Phasenraumtrajektorien des Doppelpendels bei dieser Gesamtenergie und Anfangswert stabil auf einem Attraktor-Grenzzyklus bleiben (links: Phasenraumtrajektorie der Masse 1 $\frac{d \theta_1(t)}{dt}$ vs. $\theta_1(t)$, rechts : Phasenraumtrajektorie der Masse 2 $\frac{d \theta_2(t)}{dt}$ vs. $\theta_2(t)$).

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

Wir möchten nun die Bewegung des Pendels zusammen mit seiner Phasenraumtrajektorie in einer Animation veranschaulichen. Wir erstellen diese Animation mittels des Pythonmoduls Matplotlib und benutzen hierbei die Funktion "FuncAnimation(...)" (siehe Animations using Matplotlib).

import matplotlib.animation as animation
from IPython.display import HTML

plot_max=set_l_1 + set_l_2 + 0.2
set_frames=150
stepi=int(N/set_frames)
fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_1(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[0], Loes.y[2],c="blue",alpha=0.2)
    ax2.plot(Loes.y[0][0:stepi*i], Loes.y[2][0:stepi*i],c="blue")
    ax2.scatter(Loes.y[0][stepi*i], Loes.y[2][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[1], Loes.y[3],c="red",alpha=0.2)
    ax2.plot(Loes.y[1][0:stepi*i], Loes.y[3][0:stepi*i],c="red")
    ax2.scatter(Loes.y[1][stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

In den beiden oberen Animationen erkennt man gut, dass die beiden Pendel in gleicher Richtung schwingen.

Wir ändern jetzt die Anfangsbedingungen (nun $\theta_1(0) = \pi/14 \,\, , \,\, \theta_2(0) = - \sqrt{2} \cdot \theta_1(0)$) bei gleicher Gesamtenergie.

theta2_0 = - np.sqrt(2)*theta1_0
dtheta1_0 = 0
dtheta2_0 = 0
u_init_b = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]
fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init_b, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  0.9817662576217473
Potentielle Energie V=  0.9817662576217473
Kinetische Energie T =  0.0

Die untere Abbildung zeigt, wie sich die Winkel $\theta_1$ und $\theta_2$ des Doppelpendels im Laufe der Zeit verändern. Man erkennt wieder eine sehr gleichförmige Bewegung der beiden Pendel, wobei die Pendel nun entgegengesetzt schwingen.

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \theta_1, \, \theta_2$")
plt.plot(Loes.t, Loes.y[0],c="blue", label=r"$\rm \theta_1(t)$")
plt.plot(Loes.t, Loes.y[1],c="red", label=r"$\rm \theta_2(t)$")
plt.legend(loc='lower right',fontsize=16);

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \dot{\theta_1}, \, \dot{\theta_2}$")
plt.plot(Loes.t, Loes.y[2],c="blue", label=r"$\rm \dot{\theta}_1(t)$")
plt.plot(Loes.t, Loes.y[3],c="red", label=r"$\rm \dot{\theta}_2(t)$")
plt.legend(loc='lower right',fontsize=16);

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm Energie \, E \, , \, \, T  \, , \, \, V$")
plt.plot(Loes.t, E_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="black", label=r"$\rm Gesamtenergie \, E$");
plt.plot(Loes.t, T_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="green", label=r"$\rm kinetische \, Energie \, T$");
plt.plot(Loes.t, V_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="orange", label=r"$\rm potentielle \, Energie \, V$");
plt.legend(loc='upper right',fontsize=16);

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Pendel 1")
ax2.set_title("Pendel 2")
ax1.set_xlabel(r"$\rm x_1$")
ax1.set_ylabel(r"$\rm y_1$")
ax2.set_xlabel(r"$\rm x_2$")
ax2.set_ylabel(r"$\rm y_2$")

ax1.plot(set_l_1 * np.sin(Loes.y[0]), - set_l_1 * np.cos(Loes.y[0]),c="blue")
ax2.plot(set_l_1 * np.sin(Loes.y[0]) + set_l_2 * np.sin(Loes.y[1]), - set_l_1 * np.cos(Loes.y[0]) - set_l_2 * np.cos(Loes.y[1]),c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_1(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[0], Loes.y[2],c="blue",alpha=0.2)
    ax2.plot(Loes.y[0][0:stepi*i], Loes.y[2][0:stepi*i],c="blue")
    ax2.scatter(Loes.y[0][stepi*i], Loes.y[2][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[1], Loes.y[3],c="red",alpha=0.2)
    ax2.plot(Loes.y[1][0:stepi*i], Loes.y[3][0:stepi*i],c="red")
    ax2.scatter(Loes.y[1][stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Wir möchten nun noch weitere Simulationen bei gleicher Energie $E_0 = 0.9817662576217473$ machen. Dazu betrachten wir uns den Ausdruck der Gesamtenergie $E_0 = E(\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2)$:

theta_1_0, theta_2_0 = symbols("theta_1, theta_2")
dtheta_1_0 = symbols(r"\dot{\theta}_1")
dtheta_2_0 = symbols(r"\dot{\theta}_2")
E_0 = symbols(r"E_0")

GL_E = Eq(E_0,simplify(E.subs({(diff(theta_1,t),dtheta_1_0),(diff(theta_2,t),dtheta_2_0),(theta_1,theta_1_0),(theta_2,theta_2_0)})))
GL_E

$\displaystyle E_{0} = 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{1} - g l_{1} m_{1} \left(\cos{\left(\theta_{1} \right)} - 1\right) - g m_{2} \left(l_{1} \cos{\left(\theta_{1} \right)} - l_{1} + l_{2} \cos{\left(\theta_{2} \right)} - l_{2}\right) + 0.5 m_{2} \left(\dot{\theta}_1^{2} l_{1}^{2} + 2 \dot{\theta}_1 \dot{\theta}_2 l_{1} l_{2} \cos{\left(\theta_{1} - \theta_{2} \right)} + \dot{\theta}_2^{2} l_{2}^{2}\right)$

Bei gegebenen Parameterwerten $m_1, m_2, l_1, l_2$ und $g$ muss man nun die Anfangsbedingungen $\theta_1(t_0),\theta_2(t_0),\dot{\theta}_1(t_0)$ und $\dot{\theta}_2(t_0)$ so wählen, dass der vorgegebene Energiewert $E_0$ eingehalten wird. Um dies zu gewährleisten, formen wir den oberen Ausdruck nach $\dot{\theta}_2$ um und benutzen dann diese Funktion $\dot{\theta}_2(\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2,E_0)$ um unterschiedliche Anfangswerte gleicher Gesamtenergie zu generieren. Man erhält formal zwei Lösungen:

GL_E_dt2 = Eq(dtheta_2_0,solve(GL_E,dtheta_2_0)[1])
GL_E_dt2

$\displaystyle \dot{\theta}_2 = - \frac{\dot{\theta}_1 l_{1} \cos{\left(\theta_{1} - \theta_{2} \right)}}{l_{2}} + \frac{1.4142135623731 \sqrt{E_{0} - 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{1} + 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{2} \cos^{2}{\left(\theta_{1} - \theta_{2} \right)} - 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{2} + g l_{1} m_{1} \cos{\left(\theta_{1} \right)} - g l_{1} m_{1} + g l_{1} m_{2} \cos{\left(\theta_{1} \right)} - g l_{1} m_{2} + g l_{2} m_{2} \cos{\left(\theta_{2} \right)} - g l_{2} m_{2}}}{l_{2} \sqrt{m_{2}}}$

Eq(dtheta_2_0,solve(GL_E,dtheta_2_0)[0])

$\displaystyle \dot{\theta}_2 = - \frac{\dot{\theta}_1 l_{1} \cos{\left(\theta_{1} - \theta_{2} \right)}}{l_{2}} - \frac{1.4142135623731 \sqrt{E_{0} - 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{1} + 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{2} \cos^{2}{\left(\theta_{1} - \theta_{2} \right)} - 0.5 \dot{\theta}_1^{2} l_{1}^{2} m_{2} + g l_{1} m_{1} \cos{\left(\theta_{1} \right)} - g l_{1} m_{1} + g l_{1} m_{2} \cos{\left(\theta_{1} \right)} - g l_{1} m_{2} + g l_{2} m_{2} \cos{\left(\theta_{2} \right)} - g l_{2} m_{2}}}{l_{2} \sqrt{m_{2}}}$

Beim Generieren von unterschiedlichen Lösungen muss man darauf achten, dass der Wert unter der Wurzel positiv bleibt; dies ist nicht automatisch bei jeder Energie erfüllt.

Nimmt man nun z.B. noch die folgende Nebenbedingung an: $\theta_1=\theta_2=0$ (beide Massen werden aus der Ruhelage angestoßen) und legt die Parameter fest ($m_1=1, m_2=1, l_1=1, l_2=1$ und $g=9.81$), so erhält man die folgende Gleichung.

GL_E_2 = Eq(dtheta_2_0,GL_E_dt2.rhs.subs({(g,set_g),(m_1,set_m_1),(m_2,set_m_2),(l_1,set_l_1),(l_2,set_l_2),(theta_1_0,0),(theta_2_0,0)}))
GL_E_2

$\displaystyle \dot{\theta}_2 = - \dot{\theta}_1 + 1.4142135623731 \sqrt{E_{0} - 0.5 \dot{\theta}_1^{2}}$

Man hätte jedoch auch annehmen können, dass beide Pendel zur Zeit $t=0$ ruhen ($\dot{\theta}_1(t_0)=\dot{\theta}_2(t_0)=0$) und die dann entstehende Gleichung nach $\theta_2$ aufzulösen. Wieder erhält man zwei Lösungen

GL_E_t2 = Eq(theta_2_0,simplify(solve(GL_E.subs(dtheta_1_0,0).subs(dtheta_2_0,0),theta_2_0)[1]))
GL_E_t2

$\displaystyle \theta_{2} = \operatorname{acos}{\left(\frac{- E_{0} + g l_{1} m_{1} \cdot \left(1 - \cos{\left(\theta_{1} \right)}\right) + g l_{1} m_{2} \cdot \left(1 - \cos{\left(\theta_{1} \right)}\right) + g l_{2} m_{2}}{g l_{2} m_{2}} \right)}$

Eq(theta_2_0,simplify(solve(GL_E.subs(dtheta_1_0,0).subs(dtheta_2_0,0),theta_2_0)[0]))

$\displaystyle \theta_{2} = - \operatorname{acos}{\left(- \frac{E_{0} + g l_{1} m_{1} \left(\cos{\left(\theta_{1} \right)} - 1\right) + g l_{1} m_{2} \left(\cos{\left(\theta_{1} \right)} - 1\right) - g l_{2} m_{2}}{g l_{2} m_{2}} \right)} + 2 \pi$

Nach Festlegung der Parameter ($m_1=1, m_2=1, l_1=1, l_2=1$ und $g=9.81$) erhält man dann:

GL_E_1 = Eq(theta_2_0,GL_E_t2.rhs.subs({(g,set_g),(m_1,set_m_1),(m_2,set_m_2),(l_1,set_l_1),(l_2,set_l_2)}))
GL_E_1

$\displaystyle \theta_{2} = \operatorname{acos}{\left(- 0.101936799184506 E_{0} - 2.0 \cos{\left(\theta_{1} \right)} + 3.0 \right)}$

Wir betrachten nun die entstehenden Bewegungen zu jeweils einer der beiden Gleichungen; speziell für $E_0 = 0.9817662576217473$:

$\dot{\theta}_1=\dot{\theta}_2=\theta_1=0 \, , \quad \theta_2=0.451205974942728$

$\theta_1=\theta_2=\dot{\theta}_1=0 \, , \quad \dot{\theta}_2=1.40126104464639$

GL_E_1.subs(theta_1_0,0).subs(E_0,0.9817662576217473)

$\displaystyle \theta_{2} = 0.451205974942728$

GL_E_2.subs(dtheta_1_0,0).subs(E_0,0.9817662576217473)

$\displaystyle \dot{\theta}_2 = 1.40126104464639$

Wir stellen im folgenden nur die Lösung mit den Anfangsbedingungen $\theta_1=\theta_2=\dot{\theta}_1=0 $ und $ \dot{\theta}_2=1.40126104464639$ dar.

theta1_0 = 0
dtheta1_0 = 0
dtheta2_0 = 0
theta2_0 = float(GL_E_1.rhs.subs(theta_1_0,theta1_0).subs(E_0,0.9817662576217473))
u_init_c = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]

theta1_0 = 0
theta2_0 = 0
dtheta1_0 = 0
dtheta2_0 = float(GL_E_2.rhs.subs(dtheta_1_0,dtheta1_0).subs(E_0,0.9817662576217473))
u_init_d = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]

fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init_d, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  0.9817662576217473
Potentielle Energie V=  0.0
Kinetische Energie T =  0.9817662576217473

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \theta_1, \, \theta_2$")
plt.plot(Loes.t, Loes.y[0],c="blue", label=r"$\rm \theta_1(t)$")
plt.plot(Loes.t, Loes.y[1],c="red", label=r"$\rm \theta_2(t)$")
plt.legend(loc='lower right',fontsize=16);

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \dot{\theta_1}, \, \dot{\theta_2}$")
plt.plot(Loes.t, Loes.y[2],c="blue", label=r"$\rm \dot{\theta}_1(t)$")
plt.plot(Loes.t, Loes.y[3],c="red", label=r"$\rm \dot{\theta}_2(t)$")
plt.legend(loc='lower right',fontsize=16);

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm Energie \, E \, , \, \, T  \, , \, \, V$")
plt.plot(Loes.t, E_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="black", label=r"$\rm Gesamtenergie \, E$");
plt.plot(Loes.t, T_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="green", label=r"$\rm kinetische \, Energie \, T$");
plt.plot(Loes.t, V_L(Loes.y[0],Loes.y[1],Loes.y[2],Loes.y[3],1,1,set_l_1,set_l_2),c="orange", label=r"$\rm potentielle \, Energie \, V$");
plt.legend(loc='upper right',fontsize=16);

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Pendel 1")
ax2.set_title("Pendel 2")
ax1.set_xlabel(r"$\rm x_1$")
ax1.set_ylabel(r"$\rm y_1$")
ax2.set_xlabel(r"$\rm x_2$")
ax2.set_ylabel(r"$\rm y_2$")

ax1.plot(set_l_1 * np.sin(Loes.y[0]), - set_l_1 * np.cos(Loes.y[0]),c="blue")
ax2.plot(set_l_1 * np.sin(Loes.y[0]) + set_l_2 * np.sin(Loes.y[1]), - set_l_1 * np.cos(Loes.y[0]) - set_l_2 * np.cos(Loes.y[1]),c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.32, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_1(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[0], Loes.y[2],c="blue",alpha=0.2)
    ax2.plot(Loes.y[0][0:stepi*i], Loes.y[2][0:stepi*i],c="blue")
    ax2.scatter(Loes.y[0][stepi*i], Loes.y[2][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.32, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[1], Loes.y[3],c="red",alpha=0.2)
    ax2.plot(Loes.y[1][0:stepi*i], Loes.y[3][0:stepi*i],c="red")
    ax2.scatter(Loes.y[1][stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Man erkennt, dass diese Bewegung zwar komplizierter erscheint, aber anscheinend keinen chaotischen Charakter besitzt. Um die Frage nach einem möglichen chaotischen Verhalten der Bewegung näher zu betrachten, werden wir uns im Folgenden die Trajektorie der Simulation im Phasenraum ansehen und die Schnittpunkte zur Hyperfläche $\theta_2=0$ (bzw. $\theta_1=0$) erstellen (Poincaré Abbildung oder Poincaré-Schnitt).

Die Poincaré Abbildung

Wir werden nun eine Poincaré Abbildung durchführen und speziell den Poincaré-Schnitt bei $\theta_2=0$ betrachten. Zusätzlich wollen wir nur die Poincaré-Punkte berücksichtigen, bei welchen die Trajektorie des zweiten Pendels die $\theta_2=0$-Ebene der Hypersphäre des Phasenraumes mit positiver Geschwindigkeit in x-Richtung durchstösst. Bei $\theta_2=0$ gilt für die Geschwindigkeit des zweiten Pendels in x-Richtung.

diff(x_2,t)

$\displaystyle l_{1} \cos{\left(\theta_{1}{\left(t \right)} \right)} \frac{d}{d t} \theta_{1}{\left(t \right)} + l_{2} \cos{\left(\theta_{2}{\left(t \right)} \right)} \frac{d}{d t} \theta_{2}{\left(t \right)}$

Bei $\theta_2=0$ gilt somit für die Geschwindigkeit des zweiten Pendels in x-Richtung $\dot{x}(t)= l_1\,\hbox{cos}(\theta_1) \, \dot{\theta}_1 + l_1\,\dot{\theta}_2$ und wir haben somit die folgenden Bedingungen bei unserem Poincaré-Schnitt

$$ \begin{equation} \theta_2(t)=0 \quad \hbox{und} \quad l_1\,\hbox{cos}(\theta_1(t)) \, \dot{\theta}_1(t) + l_1\,\dot{\theta}_2(t) > 0 \quad \forall \, t \in {\cal T} = [t_1,t_2,...,t_n] \quad, \end{equation} $$

wobei ${\cal T}$ die Menge der Zeitpunkte ist, bei denen die Trajektorie die ($\theta_2=0$) - Ebene in positiver Richtung durchläuft.

Um die gesamte Vorgehensweise der Poincaré-Abbindung zu verdeutlichen, berechnen wir nun diese Poincaré-Punkte und stellen sie zusammen mit der Phasenraumtrajektorie dar.

t1 = Loes.y[0]
t2 = Loes.y[1]
dt1 = Loes.y[2]
dt2 = Loes.y[3]

t_n_index = []
for i in range(len(t1)-1):
    if t2[i]*t2[i+1] < 0 and  set_l_1*np.cos(t1[i])*dt1[i] + set_l_2*dt2[i] > 0:
        t_n_index.append(i)

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

point=1
i_end=t_n_index[point]
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(Loes.y[0][0], Loes.y[2][0], -np.max(np.abs(Loes.y[1])), marker='o', s=30, c = "yellow")
ax.scatter3D(Loes.y[0][0], Loes.y[2][0], Loes.y[1][0], marker='o', s=30, c = "yellow")
ax.scatter3D(Loes.y[0][t_n_index[0:point+1]], Loes.y[2][t_n_index[0:point+1]], [-np.max(np.abs(Loes.y[1]))]*(point+1), marker='o', s=30, c = "black")
ax.scatter3D(Loes.y[0][t_n_index[0:point+1]], Loes.y[2][t_n_index[0:point+1]], Loes.y[1][t_n_index[0:point+1]], marker='o', s=30, c = "black")
bild = ax.scatter3D(Loes.y[0][0:i_end], Loes.y[2][0:i_end], Loes.y[1][0:i_end], marker='o', s=2, alpha=0.3, c = Loes.y[1][0:i_end], cmap=cm.seismic, vmin=-np.max(np.abs(Loes.y[1])), vmax=np.max(np.abs(Loes.y[1])))
ax.view_init(azim=20, elev=20)
cbar = fig.colorbar(bild, shrink=0.5, aspect=10)
cbar.set_label(r"$\rm Winkel \, \theta_2$", rotation=270)
plt.title(r"$\rm Schnitt \, \theta_2=0 \,\, und \,\, \frac{d x_2}{dt}>0$")
ax.set_xlabel(r"$\rm \theta_1$")
ax.set_ylabel(r"$\rm \frac{d \theta_1}{dt}$")
ax.set_zlabel(r"$\rm \theta_2$");

Die obere Abbildung zeigt die Phasenraumtrajektorie in einem dreidimensionalen Koordinatensystem $(x,y,z) \rightarrow (\theta_1 , \frac{d \theta_1}{dt} , \theta_2)$. Die Farbskala der Trajektorie ist so gewählt, dass positive $\theta_2$-Werte rot und negative blau erscheinen, sodass der Übergang von blau zu rot gerade die Poincaré-Bedingung ( $\theta_2(t)=0 \, , \, l_1\,\hbox{cos}(\theta_1(t)) \, \dot{\theta}_1(t) + l_1\,\dot{\theta}_2(t) > 0$) repräsentiert. Die Trajektorie ist bis zum zweiten Poincaré-Punkt gezeichnet, wobei der Anfangswert der Trajektorie als gelber Punkt und die beiden Poincaré-Punkte als schwarze Punkte auf der Trajektorie markiert sind. Zusätzlich ist auch die Projektion der Punkte auf die untere $(\theta_1 , \frac{d \theta_1}{dt})$-Ebene in der Grafik visualisiert.

Betrachtet man sich die Trajektorie zu längeren Zeiten, dann entsteht in der $(\theta_1 , \frac{d \theta_1}{dt})$-Ebene eine geometrische Form der Punkte (ein einzelner Punkt, eine Kurve oder ein begrenzter Flächenbereich), wobei reguläre (nicht-chaotische Bewegungen) durch einen Punkt oder eine Kurve gekennzeichnet sind.

Die folgende Abbildung zeigt den Poincaré-Schnitt für einige weitere Poincaré-Punkte.

point=len(t_n_index)-1
i_end=t_n_index[point]
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(Loes.y[0][0], Loes.y[2][0], -np.max(np.abs(Loes.y[1])), marker='o', s=30, c = "yellow")
ax.scatter3D(Loes.y[0][0], Loes.y[2][0], Loes.y[1][0], marker='o', s=30, c = "yellow")
ax.scatter3D(Loes.y[0][t_n_index[0:point+1]], Loes.y[2][t_n_index[0:point+1]], [-np.max(np.abs(Loes.y[1]))]*(point+1), marker='o', s=30, c = "black")
ax.scatter3D(Loes.y[0][t_n_index[0:point+1]], Loes.y[2][t_n_index[0:point+1]], Loes.y[1][t_n_index[0:point+1]], marker='o', s=30, c = "black")
bild = ax.scatter3D(Loes.y[0][0:i_end], Loes.y[2][0:i_end], Loes.y[1][0:i_end], marker='o', s=2, alpha=0.3, c = Loes.y[1][0:i_end], cmap=cm.seismic, vmin=-np.max(np.abs(Loes.y[1])), vmax=np.max(np.abs(Loes.y[1])))
ax.view_init(azim=20, elev=20)
cbar = fig.colorbar(bild, shrink=0.5, aspect=10)
cbar.set_label(r"$\rm Winkel \, \theta_2$", rotation=270)
plt.title(r"$\rm Schnitt \, \theta_2=0 \,\, und \,\, \frac{d x_2}{dt}>0$")
ax.set_xlabel(r"$\rm \theta_1$")
ax.set_ylabel(r"$\rm \frac{d \theta_1}{dt}$")
ax.set_zlabel(r"$\rm \theta_2$");

Um eine aussagekräftige Poincaré-Abbildung zu erhalten, muss man jedoch eine längere Simulation durchführen. Wir erhöhen deshalb den Endwert auf $t_{end}=500$.

t_end_long = 500
N_long = 300000
t_val_long = np.linspace(0, t_end_long, N_long )
fehler_long = 10**(-10)

Loes_long = solve_ivp(DGLsysa, [0, t_end_long], u_init_d, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)

t1 = Loes_long.y[0]
t2 = Loes_long.y[1]
dt1 = Loes_long.y[2]
dt2 = Loes_long.y[3]

t_n_index = []
for i in range(len(t1)-1):
    if t2[i]*t2[i+1] < 0 and  set_l_1*np.cos(t1[i])*dt1[i] + set_l_2*dt2[i] > 0:
        t_n_index.append(i)

point=len(t_n_index)-1
i_end=t_n_index[point]
i_end1=t_n_index[10]
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(Loes_long.y[0][0], Loes_long.y[2][0], -np.max(np.abs(Loes_long.y[1])), marker='o', s=30, c = "yellow")
ax.scatter3D(Loes_long.y[0][0], Loes_long.y[2][0], Loes_long.y[1][0], marker='o', s=30, c = "yellow")
ax.scatter3D(Loes_long.y[0][t_n_index[0:point+1]], Loes_long.y[2][t_n_index[0:point+1]], [-np.max(np.abs(Loes_long.y[1]))]*(point+1), marker='o', s=30, c = "black")
ax.scatter3D(Loes_long.y[0][t_n_index[0:point+1]], Loes_long.y[2][t_n_index[0:point+1]], Loes_long.y[1][t_n_index[0:point+1]], marker='o', s=30, c = "black")
bild = ax.scatter3D(Loes_long.y[0][0:i_end1], Loes_long.y[2][0:i_end1], Loes_long.y[1][0:i_end1], marker='o', s=2, alpha=0.3, c = Loes_long.y[1][0:i_end1], cmap=cm.seismic, vmin=-np.max(np.abs(Loes_long.y[1])), vmax=np.max(np.abs(Loes_long.y[1])))
ax.view_init(azim=20, elev=20)
cbar = fig.colorbar(bild, shrink=0.5, aspect=10)
cbar.set_label(r"$\rm Winkel \, \theta_2$", rotation=270)
plt.title(r"$\rm Schnitt \, \theta_2=0 \,\, und \,\, \frac{d x_2}{dt}>0$")
ax.set_xlabel(r"$\rm \theta_1$")
ax.set_ylabel(r"$\rm \frac{d \theta_1}{dt}$")
ax.set_zlabel(r"$\rm \theta_2$");

In der obigen Abbildung haben wir, der Übersichtlichkeit halber, die Trajektorie nur bis zum 10-ten Poincaré-Punkt dargestellt. Man erkennt gut, dass die Poincaré-Abbildung der Simulation eine Kurve liefert, was einer nicht-chaotischen Bewegung entspricht.

Wir wollen nun noch einen weiteren Poincaré-Schnitt durchführen und nehmen die folgenden Bedingungen an

$$ \begin{equation} \theta_1(t)=0 \quad \hbox{und} \quad \dot{\theta}_1(t) > 0 \quad \forall \, t \in {\cal T} = [t_1,t_2,...,t_n] \quad, \end{equation} $$

wobei ${\cal T}$ nun die Menge der Zeitpunkte ist, bei denen die Trajektorie die ($\theta_1=0$) - Ebene in positiver Richtung durchläuft.

t_n_indexa = []
for i in range(len(t1)-1):
    if t1[i]*t1[i+1] < 0 and dt1[i] > 0:
        t_n_indexa.append(i)

point=len(t_n_indexa)-1
i_end=t_n_indexa[point]
i_end1=t_n_indexa[10]
fig = plt.figure(figsize = (16, 9))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(Loes_long.y[1][0], Loes_long.y[3][0], -np.max(np.abs(Loes_long.y[0])), marker='o', s=30, c = "yellow")
ax.scatter3D(Loes_long.y[1][0], Loes_long.y[3][0], Loes_long.y[0][0], marker='o', s=30, c = "yellow")
ax.scatter3D(Loes_long.y[1][t_n_indexa[0:point+1]], Loes_long.y[3][t_n_indexa[0:point+1]], [-np.max(np.abs(Loes_long.y[0]))]*(point+1), marker='o', s=30, c = "black")
ax.scatter3D(Loes_long.y[1][t_n_indexa[0:point+1]], Loes_long.y[3][t_n_indexa[0:point+1]], Loes_long.y[0][t_n_indexa[0:point+1]], marker='o', s=30, c = "black")
bild = ax.scatter3D(Loes_long.y[1][0:i_end1], Loes_long.y[3][0:i_end1], Loes_long.y[0][0:i_end1], marker='o', s=2, alpha=0.3, c = Loes_long.y[1][0:i_end1], cmap=cm.seismic, vmin=-np.max(np.abs(Loes_long.y[0])), vmax=np.max(np.abs(Loes_long.y[0])))
ax.view_init(azim=20, elev=20)
cbar = fig.colorbar(bild, shrink=0.5, aspect=10)
cbar.set_label(r"$\rm Winkel \, \theta_1$", rotation=270)
plt.title(r"$\rm Schnitt \, \theta_1=0 \,\, und \,\, \frac{d x_1}{dt}>0$")
ax.set_xlabel(r"$\rm \theta_2$")
ax.set_ylabel(r"$\rm \frac{d \theta_2}{dt}$")
ax.set_zlabel(r"$\rm \theta_1$");

In der obigen Abbildung haben wir wieder, der Übersichtlichkeit halber, die Trajektorie nur bis zum 10-ten Poincaré-Punkt dargestellt. Man erkennt auch hier gut, dass die Poincaré-Abbildung der Simulation eine Kurve liefert, was einer nicht-chaotischen Bewegung entspricht.

Zusammenfassend erhalten wir somit folgende Poincaré-Schnitte für die Simulation.

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title(r"$\rm Schnitt \, \theta_2=0 \,\, , \,\, \dot{x}_2>0$")
ax2.set_title(r"$\rm Schnitt \, \theta_1=0 \,\, , \,\, \dot{x}_1>0$")
ax1.set_xlabel(r"$\rm \theta_1(t_n)$")
ax1.set_ylabel(r"$\rm \dot{\theta_1}(t_n)$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta_2}(t)$")

ax1.scatter(Loes_long.y[0][t_n_index], Loes_long.y[2][t_n_index],c = "blue",s=10)


ax2.scatter(Loes_long.y[1][t_n_indexa], Loes_long.y[3][t_n_indexa],c = "red",s=10);

Wir berechnen diese Schnitte nun auch für die anderen drei Anfangsbedingungen.

t_end_long = 1000
N_long = 500000
t_val_long = np.linspace(0, t_end_long, N_long )
fehler_long = 10**(-10)

u1 = []
u2 = []
u3 = []
u4 = []
u_init_long = [u_init_a,u_init_b,u_init_c,u_init_d]
for u_init in u_init_long:
    Loes_long = solve_ivp(DGLsysa, [0, t_end_long], u_init, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)

    t1 = Loes_long.y[0]
    t2 = Loes_long.y[1]
    dt1 = Loes_long.y[2]
    dt2 = Loes_long.y[3]

    t_n_index = []
    t_n_indexa = []
    for i in range(len(t1)-1):
        if t2[i]*t2[i+1] < 0 and  set_l_1*np.cos(t1[i])*dt1[i] + set_l_2*dt2[i] > 0:
            t_n_index.append(i)
        if t1[i]*t1[i+1] < 0 and dt1[i] > 0:
            t_n_indexa.append(i)
    u1.append([t1[t_n_index],t1[t_n_indexa]])
    u2.append([t2[t_n_index],t2[t_n_indexa]])
    u3.append([dt1[t_n_index],dt1[t_n_indexa]])
    u4.append([dt2[t_n_index],dt2[t_n_indexa]])

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title(r"$\rm Schnitt \, \theta_2=0 \,\, , \,\, \dot{x}_2>0$")
ax2.set_title(r"$\rm Schnitt \, \theta_1=0 \,\, , \,\, \dot{x}_1>0$")
ax1.set_xlabel(r"$\rm \theta_1(t_n)$")
ax1.set_ylabel(r"$\rm \dot{\theta_1}(t_n)$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta_2}(t)$")

ax1.scatter(u1[0][0], u3[0][0],c = "blue",s=10, label=r"$\rm a$")
ax2.scatter(u2[0][1], u4[0][1],c = "blue",s=10, label=r"$\rm a$")
ax1.scatter(u1[1][0], u3[1][0],c = "red",s=10, label=r"$\rm b$")
ax2.scatter(u2[1][1], u4[1][1],c = "red",s=10, label=r"$\rm b$")
ax1.scatter(u1[2][0], u3[2][0],c = "orange",s=10, label=r"$\rm c$")
ax2.scatter(u2[2][1], u4[2][1],c = "orange",s=10, label=r"$\rm c$")
ax1.scatter(u1[3][0], u3[3][0],c = "black",s=10, label=r"$\rm d$")
ax2.scatter(u2[3][1], u4[3][1],c = "black",s=10, label=r"$\rm d$")

ax1.legend(loc='lower right',fontsize=16)
ax2.legend(loc='upper right',fontsize=16);

Die obere Abbildung zeigt die beiden Poincaré-Schnitte für alle vier Simulationen, wobei die Farben den unterschiedlichen Anfangsbedingungen der Simulationen entsprechen:

Anfangsbedingung a (blau): $\theta_1(0) = \pi/14 \,\, , \,\, \theta_2(0) = \sqrt{2} \cdot \theta_1(0) \,\, , \,\, \dot{\theta}_1(0)=\dot{\theta}_2(0)=0$

Anfangsbedingung b (rot): $\theta_1(0) = \pi/14 \,\, , \,\, \theta_2(0) = - \sqrt{2} \cdot \theta_1(0) \,\, , \,\, \dot{\theta}_1(0)=\dot{\theta}_2(0)=0$

Anfangsbedingung c (orange): $\dot{\theta}_1=\dot{\theta}_2=\theta_1=0 \, , \quad \theta_2=0.451205974942728$

Anfangsbedingung d (schwarz): $\theta_1=\theta_2=\dot{\theta}_1=0 \, , \quad \dot{\theta}_2=1.40126104464639$

Die Anfangsbedingungen a und b liefern nur einen Punkt in ihren Poincaré-Schnitten, wobei die Simulationen mit den Anfangsbedingungen c und d nicht durchgezogene Kurven liefern. Zusammengefasst sind bei dieser geringen Energie alle untersuchten Simulationen nicht chaotisch.

Für eine wissenschaftliche Analyse muss man jedoch noch viele weitere Anfangsbedingungen untersuchen, was bei unserer Vorgehensweise hier mittels Python sehr zeitintensiv wäre. Die folgende Abbildung wurde mit einem entsprechenden parallelisierten C++ Programm erzeugt. Sie zeigt den Poincaré-Schnitt bei der Gesamtenergie von $E_0 = 0.9817662576217473$ für sehr viele unterschiedliche Anfangsbedingungen.

Reguläre und chaotische Bewegungen bei mittlerer Energie

Wir werden nun die Bewegung des Doppelpendels bei höheren Gesamtenergien studieren. Wir starten wieder mit einem gleichförmigen Pendeln mit $E_0=11.200769559519784$.

theta1_0 = np.pi/4
theta2_0 = np.sqrt(2)*theta1_0
dtheta1_0 = 0
dtheta2_0 = 0
u_init1_a = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]
fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init1_a, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  11.200769559519784
Potentielle Energie V=  11.200769559519784
Kinetische Energie T =  0.0

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[1], Loes.y[3],c="red",alpha=0.2)
    ax2.plot(Loes.y[1][0:stepi*i], Loes.y[3][0:stepi*i],c="red")
    ax2.scatter(Loes.y[1][stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Die obere Animation zeigt wiederum eine sehr gleichförmige Bewegung des Doppelpendels. Wir berechnen nun das gegenentsetzte Pendeln bei gleicher Energie.

theta1_0 = np.pi/4
theta2_0 = - np.sqrt(2)*theta1_0
dtheta1_0 = 0
dtheta2_0 = 0
u_init1_b = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]
fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init1_b, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  11.200769559519784
Potentielle Energie V=  11.200769559519784
Kinetische Energie T =  0.0

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[1], Loes.y[3],c="red",alpha=0.2)
    ax2.plot(Loes.y[1][0:stepi*i], Loes.y[3][0:stepi*i],c="red")
    ax2.scatter(Loes.y[1][stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Zusätzlich berechnen wir noch zwei weitere Bewegungen bei gleicher Energie und stellen eine davon dar.

theta1_0 = 0
dtheta1_0 = 0
dtheta2_0 = 0
theta2_0 = float(GL_E_1.rhs.subs(theta_1_0,theta1_0).subs(E_0,11.200769559519784))
u_init1_c = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]

theta1_0 = 0
theta2_0 = 0
dtheta1_0 = 0
dtheta2_0 = float(GL_E_2.rhs.subs(dtheta_1_0,dtheta1_0).subs(E_0,11.200769559519784))
u_init1_d = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]

fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init1_c, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  11.200769559519784
Potentielle Energie V=  0.0
Kinetische Energie T =  11.200769559519784

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)
    ax2.plot(Loes.y[1], Loes.y[3],c="red",alpha=0.2)
    ax2.plot(Loes.y[1][0:stepi*i], Loes.y[3][0:stepi*i],c="red")
    ax2.scatter(Loes.y[1][stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Man erkennt in der oberen Abbildung (Anfangsbedingung c1), dass die Bewegung der beiden Pendel schon ein wenig chaotisch erscheint. Um dies näher zu beleuchten, betrachten wir im Folgenden die Poincaré Abbildungen der Bewegungen.

Die Poincaré Abbildung

Wir berechnen nun die Poincaré-Schnitte für alle vier Anfangsbedingungen mit $E_0=11.200769559519784$.

t_end_long = 1000
N_long = 500000
t_val_long = np.linspace(0, t_end_long, N_long )
fehler_long = 10**(-10)

u1 = []
u2 = []
u3 = []
u4 = []
u_init_long = [u_init1_a,u_init1_b,u_init1_c,u_init1_d]
for u_init in u_init_long:
    Loes_long = solve_ivp(DGLsysa, [0, t_end_long], u_init, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)

    t1 = Loes_long.y[0]
    t2 = Loes_long.y[1]
    dt1 = Loes_long.y[2]
    dt2 = Loes_long.y[3]

    t_n_index = []
    t_n_indexa = []
    for i in range(len(t1)-1):
        if t2[i]*t2[i+1] < 0 and  set_l_1*np.cos(t1[i])*dt1[i] + set_l_2*dt2[i] > 0:
            t_n_index.append(i)
        if t1[i]*t1[i+1] < 0 and dt1[i] > 0:
            t_n_indexa.append(i)
    u1.append([t1[t_n_index],t1[t_n_indexa]])
    u2.append([t2[t_n_index],t2[t_n_indexa]])
    u3.append([dt1[t_n_index],dt1[t_n_indexa]])
    u4.append([dt2[t_n_index],dt2[t_n_indexa]])

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title(r"$\rm Schnitt \, \theta_2=0 \,\, , \,\, \dot{x}_2>0$")
ax2.set_title(r"$\rm Schnitt \, \theta_1=0 \,\, , \,\, \dot{x}_1>0$")
ax1.set_xlabel(r"$\rm \theta_1(t_n)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t_n)$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t)$")

ax1.scatter(u1[0][0], u3[0][0],c = "blue",s=5, label=r"$\rm a1$")
ax2.scatter(u2[0][1], u4[0][1],c = "blue",s=5, label=r"$\rm a1$")
ax1.scatter(u1[1][0], u3[1][0],c = "red",s=5, label=r"$\rm b1$")
ax2.scatter(u2[1][1], u4[1][1],c = "red",s=5, label=r"$\rm b1$")
ax1.scatter(u1[2][0], u3[2][0],c = "orange",s=5, label=r"$\rm c1$")
ax2.scatter(u2[2][1], u4[2][1],c = "orange",s=5, label=r"$\rm c1$")
ax1.scatter(u1[3][0], u3[3][0],c = "black",s=5, label=r"$\rm d1$")
ax2.scatter(u2[3][1], u4[3][1],c = "black",s=5, label=r"$\rm d1$")

ax1.legend(loc='lower right',fontsize=16)
ax2.legend(loc='upper right',fontsize=16);

Die obere Abbildung zeigt die beiden Poincaré-Schnitte für alle vier Simulationen bei $E_0=11.200769559519784$, wobei die Farben wieder den unterschiedlichen Anfangsbedingungen (a1,b1,c1,d1) der Simulationen entsprechen

Die Anfangsbedingungen a1 und b1 zeigen durchgezogene Kurven und Anfangsbedingungen d1 erzeugt eine nicht durchgezogene Kurve. Zusammengefasst sind diese Simulationen somit regulär und nicht chaotisch.

Unsere Vermutung, dass die Anfangsbedingung c1 ein chaotisches Verhalten zeigt, bestätigt sich auch beim Betrachten ihrer Poincaré-Schnitte. Dies zeigt auch die nebenstehende Abbildung, welche mit einem entsprechenden parallelisierten C++ Programm erzeugt wurde. Die orangen Poincaré-Punkte zeigen die Simulationsergebnisse mit der Anfangsbedingung c1. Man erkennt, dass diese eine begrenzte Fläche ausfüllen, was bei chaotischen Bewegungen und seltsamen Attraktoren üblich ist.

Die nebenstehende Abbildung wurde auch mit dem parallelisierten C++ Programm erzeugt. Sie zeigt den Poincaré-Schnitt bei der Gesamtenergie von $E_0=11.200769559519784$ für sehr viele unterschiedliche Anfangsbedingungen. Man erkennt eine komplizierte Struktur von regulären und chaotischen Bereichen.

Das chaotische Verhalten der Trajektorie mit der Anfangsbedingung 1c bedeutet, dass sich kleine $\epsilon$-ähnliche Störungen nach einiger Zeit so aufschaukeln, dass die Bewegung nicht mehr adäquat vorhergesagt werden kann. Dies wollen wir im Folgenden numerisch betrachten, indem wir eine Vergleichrechnung machen, bei der der $\theta_2(0)$-Wert der Anfangsbedingungen um ein $\epsilon$ (hier $\epsilon=0.00001$) verschoben ist.

t_end_long = 500
N_long = 100000
t_val_long = np.linspace(0, t_end_long, N_long )
fehler_long = 10**(-13)
epsilon=0.00001

u_init1_c_eps = [u_init1_c[0],u_init1_c[1]+epsilon,u_init1_c[2],u_init1_c[3]]
u_init1_d_eps = [u_init1_d[0],u_init1_d[1]+epsilon,u_init1_d[2],u_init1_d[3]]

Loes_long_c = solve_ivp(DGLsysa, [0, t_end_long], u_init1_c, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)
Loes_long_eps_c = solve_ivp(DGLsysa, [0, t_end_long], u_init1_c_eps, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)

Die untere Abbildung zeigt den Unterschied der berechneten zeitlichen Verläufe der Winkelkoordinaten des ersten und zweiten Pendels ($\Delta\theta_1(t)$ und $\Delta\theta_2(t)$).

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta\theta_1, \, \Delta\theta_2$")
plt.plot(Loes_long_c.t, Loes_long_c.y[0]-Loes_long_eps_c.y[0],c="blue", label=r"$\rm \Delta\theta_1(t)$")
plt.plot(Loes_long_c.t, Loes_long_c.y[1]-Loes_long_eps_c.y[1],c="red", label=r"$\rm \Delta\theta_2(t)$")
plt.legend(loc='lower right',fontsize=16);

Man erkennt, einen merklichen Unterschied der Simulationen für $t>150$. Dies wird noch deutlicher, wenn man sich die Winkelgeschwindigkeiten $\Delta\dot{\theta}_1(t)$ und $\Delta\dot{\theta}_2(t)$ betrachtet (siehe untere Abbildung). Eine genaue Berechnung der Bewegung ist bei dieser Anfangsbedingung und Energie somit nicht mehr möglich.

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta\dot{\theta}_1, \, \Delta\dot{\theta}_2$")
plt.plot(Loes_long_c.t, Loes_long_c.y[2]-Loes_long_eps_c.y[2],c="blue", label=r"$\rm \Delta\dot{\theta}_1(t)$")
plt.plot(Loes_long_c.t, Loes_long_c.y[3]-Loes_long_eps_c.y[3],c="red", label=r"$\rm \Delta\dot{\theta}_2(t)$")
plt.legend(loc='lower right',fontsize=16);

Zum Vergleich zeigen die beiden folgenden Abbildungen die Abweichungen für die Anfangsbedingung a1, welche eine reguläre, nicht-chaotische Bewegung darstellt. Man erkennt, dass die Abweichungen hier auch für große Zeiten klein bleiben.

Loes_long_d = solve_ivp(DGLsysa, [0, t_end_long], u_init1_d, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)
Loes_long_eps_d = solve_ivp(DGLsysa, [0, t_end_long], u_init1_d_eps, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta\theta_1, \, \Delta\theta_2$")
plt.plot(Loes_long_d.t, Loes_long_d.y[0]-Loes_long_eps_d.y[0],c="blue", label=r"$\rm \Delta\theta_1(t)$")
plt.plot(Loes_long_d.t, Loes_long_d.y[1]-Loes_long_eps_d.y[1],c="red", label=r"$\rm \Delta\theta_2(t)$")
plt.legend(loc='lower right',fontsize=16);

plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta\dot{\theta}_1, \, \Delta\dot{\theta}_2$")
plt.plot(Loes_long_d.t, Loes_long_d.y[2]-Loes_long_eps_d.y[2],c="blue", label=r"$\rm \Delta\dot{\theta}_1(t)$")
plt.plot(Loes_long_d.t, Loes_long_d.y[3]-Loes_long_eps_d.y[3],c="red", label=r"$\rm \Delta\dot{\theta}_2(t)$")
plt.legend(loc='lower right',fontsize=16);

Chaotische Bewegungen bei hoher Energie und Pendelüberschläge

Wir werden nun die Bewegung des Doppelpendels bei noch höheren Gesamtenergien studieren ($E_0=48.89208445507095$). Wir starten mit einer Simulation, die am Anfang nur aus potentieller Energie besteht.

theta1_0 = np.pi/1.5
theta2_0 = np.sqrt(2)*theta1_0
dtheta1_0 = 0
dtheta2_0 = 0
u_init2_a = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]
fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init2_a, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  48.89208445507095
Potentielle Energie V=  48.89208445507095
Kinetische Energie T =  0.0

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

Die rechte obere Abbildung zeigt, dass die Bewegung des zweiten Pendels mehrere Überschläge erzeugt, da der Winkel $\theta_2$ teilweise über $\pi$ bzw. unter $-\pi$ ist. Bei solchen Bewegungen werden die Phasenraumtrajektorien oft Modulo $2 \pi$ dargestellt. Mittels der unteren Funktion rechnen wir die ursprünglichen Lösungen der Winkelveränderungen Modulo $2 \pi$ um. Wir stellen dann die neuen Abbildungen der Phasenraumtrajektorien des Pendels dar.

def Mod_Winkel(liste_winkel):
    shift=0
    lwn=[]
    for w in liste_winkel:
        if(w+shift>=np.pi):
            shift=shift-2*np.pi
        if(w+shift<=-np.pi):
            shift=shift+2*np.pi
        lwn.append(w+shift)
    return np.array(lwn)

wn1 = Mod_Winkel(Loes.y[0])
wn2 = Mod_Winkel(Loes.y[1])

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.scatter(wn1, Loes.y[2],c="blue", s=1, marker='.')
ax2.scatter(wn2, Loes.y[3],c="red", s=1, marker='.');

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)    
    ax2.scatter(wn2, Loes.y[3],c="red", s=1, marker='.',alpha=0.1)
    ax2.scatter(wn2[0:stepi*i], Loes.y[3][0:stepi*i],c="red", s=1, marker='.')
    ax2.scatter(wn2[stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Wir betrachten eine weitere Simulation bei $E_0=48.89208445507095$. Auch hier sind die Anfangsbedingungen so gewählt, dass das System bei $t=0$ keine kinetische Energie besitzt.

theta1_0 = np.pi/1.5
theta2_0 = -np.sqrt(2)*theta1_0
dtheta1_0 = 0
dtheta2_0 = 0
u_init2_b = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]
fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init2_b, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  48.89208445507095
Potentielle Energie V=  48.89208445507095
Kinetische Energie T =  0.0

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

wn1 = Mod_Winkel(Loes.y[0])
wn2 = Mod_Winkel(Loes.y[1])

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.scatter(wn1, Loes.y[2],c="blue", s=1, marker='.')
ax2.scatter(wn2, Loes.y[3],c="red", s=1, marker='.');

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)    
    ax2.scatter(wn2, Loes.y[3],c="red", s=1, marker='.',alpha=0.1)
    ax2.scatter(wn2[0:stepi*i], Loes.y[3][0:stepi*i],c="red", s=1, marker='.')
    ax2.scatter(wn2[stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Wir betrachten zwei weitere Simulation bei $E_0=48.89208445507095$. Hier sind die Anfangsbedingungen so gewählt, dass das System bei $t=0$ aus der Ruhelage angestoßen wird und somit nur kinetische Energie besitzt.

theta1_0 = 0
theta2_0 = 0
dtheta1_0 = 5
dtheta2_0 = float(GL_E_2.rhs.subs(dtheta_1_0,dtheta1_0).subs(E_0,48.89208445507095))
u_init2_c = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]

theta1_0 = 0
theta2_0 = 0
dtheta1_0 = 0
dtheta2_0 = float(GL_E_2.rhs.subs(dtheta_1_0,dtheta1_0).subs(E_0,48.89208445507095))
u_init2_d = [theta1_0,theta2_0,dtheta1_0,dtheta2_0]

fehler = 10**(-13)
Loes = solve_ivp(DGLsysa, [0, t_end], u_init2_c, t_eval=t_val, args=(set_l_1,set_l_2, ), rtol=fehler, atol=fehler)
print("Gesamtenergie E = ",E_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Potentielle Energie V= ",V_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))
print("Kinetische Energie T = ",T_L(theta1_0,theta2_0,dtheta1_0,dtheta2_0,set_m_1,set_m_2,set_l_1,set_l_2))

Gesamtenergie E =  48.89208445507095
Potentielle Energie V=  0.0
Kinetische Energie T =  48.89208445507095

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.plot(Loes.y[0], Loes.y[2],c="blue")
ax2.plot(Loes.y[1], Loes.y[3],c="red");

wn1 = Mod_Winkel(Loes.y[0])
wn2 = Mod_Winkel(Loes.y[1])

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title("Phasenraumtrajektorie Pendel 1")
ax2.set_title("Phasenraumtrajektorie Pendel 2")
ax1.set_xlabel(r"$\rm \theta_1(t)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t) = \frac{d \theta_1(t)}{dt}$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")

ax1.scatter(wn1, Loes.y[2],c="blue", s=1, marker='.')
ax2.scatter(wn2, Loes.y[3],c="red", s=1, marker='.');

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28, hspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def animate(i):
    ax1.cla()
    ax2.cla()
    ax1.set_xlabel(r"$\rm x $")
    ax1.set_ylabel(r"$\rm y $")
    ax1.set_aspect('equal')
    ax2.set_xlabel(r"$\rm \theta_2(t)$")
    ax2.set_ylabel(r"$\rm \dot{\theta}_2(t) = \frac{d \theta_2(t)}{dt}$")
    ax1.set_xlim(-plot_max,plot_max)
    ax1.set_ylim(-plot_max,plot_max)
    ax1.scatter(0, 0, s=30, marker='o', c="black")
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i]), c="blue", alpha=0.4)
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i]), s=50, marker='o', c="blue")
    ax1.plot([0,set_l_1*sin(Loes.y[0][stepi*i])],[0,-set_l_1*cos(Loes.y[0][stepi*i])] ,c="black",linewidth=1)
    ax1.plot(set_l_1*np.sin(Loes.y[0][0:stepi*i])+set_l_2*np.sin(Loes.y[1][0:stepi*i]), -set_l_1*np.cos(Loes.y[0][0:stepi*i])-set_l_2*np.cos(Loes.y[1][0:stepi*i]), c="red")
    ax1.scatter(set_l_1*sin(Loes.y[0][stepi*i])+set_l_2*sin(Loes.y[1][stepi*i]), -set_l_1*cos(Loes.y[0][stepi*i])-set_l_2*cos(Loes.y[1][stepi*i]), s=50, marker='o', c="red")
    ax1.plot([set_l_1*sin(Loes.y[0][stepi*i]),set_l_1*sin(Loes.y[0][stepi*i]) + set_l_2*sin(Loes.y[1][stepi*i])],[-set_l_1*cos(Loes.y[0][stepi*i]),-set_l_1*cos(Loes.y[0][stepi*i]) - set_l_2*cos(Loes.y[1][stepi*i])] ,c="black",linewidth=1)    
    ax2.scatter(wn2, Loes.y[3],c="red", s=1, marker='.',alpha=0.1)
    ax2.scatter(wn2[0:stepi*i], Loes.y[3][0:stepi*i],c="red", s=1, marker='.')
    ax2.scatter(wn2[stepi*i], Loes.y[3][stepi*i], s=80, marker='o', c="black")
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=set_frames,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Die Poincaré Abbildung

Wir berechnen nun die Poincaré-Schnitte für alle vier Anfangsbedingungen mit $E_0=48.89208445507095$.

t_end_long = 1000
N_long = 500000
t_val_long = np.linspace(0, t_end_long, N_long )
fehler_long = 10**(-10)

u1 = []
u2 = []
u3 = []
u4 = []
u_init_long = [u_init2_a,u_init2_b,u_init2_c,u_init2_d]
for u_init in u_init_long:
    Loes_long = solve_ivp(DGLsysa, [0, t_end_long], u_init, t_eval=t_val_long, args=(set_l_1,set_l_2, ), rtol=fehler_long, atol=fehler_long)

    t1 = Mod_Winkel(Loes_long.y[0])
    t2 = Mod_Winkel(Loes_long.y[1])
    dt1 = Loes_long.y[2]
    dt2 = Loes_long.y[3]

    t_n_index = []
    t_n_indexa = []
    for i in range(len(t1)-1):
        if t2[i]*t2[i+1] < 0 and abs(t2[i]) < 1 and set_l_1*np.cos(t1[i])*dt1[i] + set_l_2*dt2[i] > 0:
            t_n_index.append(i)
        if t1[i]*t1[i+1] < 0 and abs(t1[i]) < 1 and dt1[i] > 0:
            t_n_indexa.append(i)
    u1.append([t1[t_n_index],t1[t_n_indexa]])
    u2.append([t2[t_n_index],t2[t_n_indexa]])
    u3.append([dt1[t_n_index],dt1[t_n_indexa]])
    u4.append([dt2[t_n_index],dt2[t_n_indexa]])

fig = plt.figure()
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.28)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_title(r"$\rm Schnitt \, \theta_2=0 \,\, , \,\, \dot{x}_2>0$")
ax2.set_title(r"$\rm Schnitt \, \theta_1=0 \,\, , \,\, \dot{x}_1>0$")
ax1.set_xlabel(r"$\rm \theta_1(t_n)$")
ax1.set_ylabel(r"$\rm \dot{\theta}_1(t_n)$")
ax2.set_xlabel(r"$\rm \theta_2(t)$")
ax2.set_ylabel(r"$\rm \dot{\theta}_2(t)$")

ax1.scatter(u1[0][0], u3[0][0],c = "blue",s=5, label=r"$\rm a2$")
ax2.scatter(u2[0][1], u4[0][1],c = "blue",s=5, label=r"$\rm a2$")
ax1.scatter(u1[1][0], u3[1][0],c = "red",s=5, label=r"$\rm b2$")
ax2.scatter(u2[1][1], u4[1][1],c = "red",s=5, label=r"$\rm b2$")
ax1.scatter(u1[2][0], u3[2][0],c = "orange",s=5, label=r"$\rm c2$")
ax2.scatter(u2[2][1], u4[2][1],c = "orange",s=5, label=r"$\rm c2$")
ax1.scatter(u1[3][0], u3[3][0],c = "black",s=5, label=r"$\rm d2$")
ax2.scatter(u2[3][1], u4[3][1],c = "black",s=5, label=r"$\rm d2$")

ax1.legend(loc='lower right',fontsize=16)
ax2.legend(loc='upper right',fontsize=16);

Die obere Abbildung zeigt die beiden Poincaré-Schnitte für alle vier Simulationen bei $E_0=48.89208445507095$, wobei die Farben wieder den unterschiedlichen Anfangsbedingungen (a2,b2,c2,d2) der Simulationen entsprechen. Alle vier Anfangsbedingungen erzeugen chaotische Bewegungen.

Unsere Vermutung, dass bei dieser Gesamtenergie stets ein chaotisches Verhalten erzeugt wird, bestätigt sich auch beim Betrachten der Poincaré-Schnitte von vielen weiteren Anfangsbedingungen. Dies zeigt die nebenstehende Abbildung welche mit einem entsprechenden parallelisierten C++ Programm erzeugt wurde.

Erhöht man die Energie weiter, so entstehen bei sehr hohen Energien wieder Bereiche, die reguläre (nicht chaotische) Pendelbewegungen zeigen. Dies zeigen die weiteren Abbildungen