DGLshow

Nachdem ich so vielen Differenzialgleichungssystemen [1, 2, 3, 4] begegnet bin, die sich nicht analytisch lösen lassen, habe ich mir ein Programm zur numerischen Lösung und Visualisierung derselben geschrieben.

Die grundlegende Idee zur numerischen Lösung von Differentialgleichungen ist es, die Zeit in diskreten Schritten \(\tau\) vergehen zu lassen. Nach jedem Schritt wird der Zustand so geändert, als ob sich während des Zeitschrittes nichts geändert hätte und die “Kräfte” werden entsprechend der Bewegungsgleichungen neu berechnet. Für infinitesimal kleine \(\tau \to \mathrm{d}t\) ist diese Methode schließlich exakt.

Im einfachsten Fall, dem Euler Verfahren, sähe das für ein einfaches Fadenpendel nach \(k\) Zeitschritten so aus

\begin{align*} \ddot\vartheta_{k+1} &= - mgl \sin(\vartheta_k)\\ \dot\vartheta_{k+1} &= \tau \ddot\vartheta_{k} + \dot\vartheta_{k}\\ \vartheta_{k+1} &= \tau \dot\vartheta_{k} + \vartheta_{k} \end{align*}

Unglücklicherweise hat dieses Verfahren ernsthafte Probleme mit der Energieerhaltung und braucht sehr kleine \(\tau\) für brauchbare Ergebnisse. Es gibt deutlich ausgefeiltere Methoden, wie den klassischen Runge-Kutta Algorithmus. Es gibt Methoden, den Zeitschritt \(\tau\) während der Simulation adaptiv anzupassen, um nur wenig Rechenaufwand in den wenig fehleranfälligen Phasen zu verbringen. Es gibt spezialisierte Methoden, die sehr gut für bestimmte Bewegungsgleichungen funktionieren, wie Velocity-Verlet, der oft für Molekulardynamiksimulationen eingesetzt wird.

Chaotische Systeme haben in der Regel etwas kompliziertere Bewegungsgleichungen. Das oben abgebildete Doppelpendel etwa wird, wie ich in einem anderen Post beschrieben habe durch folgendes Ungetüm beschrieben.

\begin{align*} \ddot\theta_1 &= \frac{m_2 \cos(\theta_1 - \theta_2) (l_1 \sin(\theta_1 - \theta_2) \dot\theta_1^2 - g \sin(\theta_2)) + m_2 l_2 \sin(\theta_1 - \theta_2) \dot\theta_2^2 + (m_1 + m_2) g \sin(\theta_1)}{m_2 l_1 \cos^2(\theta_1 - \theta_2) - (m_1+m_2) l_1} \\ \ddot\theta_2 &= \frac{m_2 l_2 \cos(\theta_1 - \theta_2) \sin(\theta_1 - \theta_2) \dot\theta_2^2 + (m_1+m_2) l_1 \sin(\theta_1 - \theta_2) \dot\theta_1^2 + (m_1+m_2) g \cos(\theta_1 - \theta_2) \sin(\theta_1) - (m_1+m_2) g \sin(\theta_2)}{(m_1+m_2) l_2 - m_2 l_2 \cos^2(\theta_1 - \theta_2)} \end{align*}

Anfangs empfiehlt es sich also etwas einfacheres und vertrauteres zu lösen, wie den Lorenz-Attraktor

\begin{align*} \dot{X} &= a(Y - X) \\ \dot{Y} &= X(b - Z) - Y \\ \dot{Z} &= XY - cZ \\ \end{align*}

Oder das Dreikörperproblem

\begin{align*} \ddot{\vec{x}_1} &= -\frac{Gm_2}{\left(x_1 - x_2\right)^3} (\vec{x}_1 - \vec{x}_2) - \frac{Gm_3}{\left(x_1 - x_3\right)^3} (\vec{x}_1 - \vec{x}_3)\\ \ddot{\vec{x}_2} &= -\frac{Gm_1}{\left(x_2 - x_1\right)^3} (\vec{x}_2 - \vec{x}_1) - \frac{Gm_3}{\left(x_2 - x_3\right)^3} (\vec{x}_2 - \vec{x}_3)\\ \ddot{\vec{x}_3} &= -\frac{Gm_1}{\left(x_3 - x_1\right)^3} (\vec{x}_3 - \vec{x}_1) - \frac{Gm_2}{\left(x_3 - x_2\right)^3} (\vec{x}_3 - \vec{x}_2)\\ \end{align*}

Da man das 3-Körperproblem trivial auf ein \(N\)-Körperproblem erweitern kann, habe ich hier ein “Sonnensystem” bzw. Bohrsches “Atom”-modell simuliert.

Sonnensystem

Um die obigen (bewegten) Bilder zu erzeugen und um ein bewegtes Doppelpendel für meinen Schreibtisch zu haben, — wennauch nur auf einem Bildschirm — habe ich in C++ einen adaptiven Runge-Kutta-4 Löser geschrieben, der mit den Qt Zeichenprimitiven animiert wird.

Auch wenn der Code nicht sehr aufgeräumt ist und Startwerte im Quellcode angepasst werden müssen, sind die Quellen auf GitHub: github.com/surt91/DGLshow.

Doppelpendel

Das ist ein Doppelpendel. Ein Doppelpendel ist neben dem Dreikörperproblem und dem Lorenz-Attraktor [1, 2] das Paradebeispiel für analytisch unlösbare Bewegungsgleichungen und chaotisches Verhalten. Aus diesem Grund sollte ein Doppelpendel auf keinem Schreibtisch fehlen und bietet sich als grandiose Geschenkidee für Physiker an.

Dass es analytisch unlösbar ist, lässt sich mit einem nicht rigorosen Argument anschaulich machen: Ein Blick auf die Bewegungsgleichungen:

\begin{align*} (m_1 + m_2) l_1 \ddot\vartheta_1 + m_2 l_2 \ddot\vartheta_2 \cos(\vartheta_1 - \vartheta_2) + m_2 l_2 \dot\vartheta_2^2 \sin(\vartheta_1 - \vartheta_2) + g(m_1 + m_2) \sin(\vartheta_1) &= 0\\ m_2 l_2 \ddot\vartheta_2 + m_2 l_1 \ddot\vartheta_1 \cos(\vartheta_1 - \vartheta_2) - m_2 l_1 \dot\vartheta_1^2 \sin(\vartheta_1 - \vartheta_2) + m_2 g \sin(\vartheta_2) &= 0 \end{align*}

Das sind die Differentialgleichungen für die beiden Winkel \(\vartheta_1\) und \(\vartheta_2\) des Doppelpendels. \(m_i\) sind die beiden Massen und \(l_i\) die Fadenlängen.

Unser Ziel ist es das obige Video zu erstellen, dazu müssen wir die Bahnkurve, also \(\vartheta_1(t)\) und \(\vartheta_2(t)\) bestimmen. Dazu müssen wir die obigen Gleichungen, die sich relativ simpel, wenn auch mühsam, per Lagrange-Formalismus herleiten lassen, zunächst nach den Winkelbeschleunigungen aufgelösen.

\begin{align*} \ddot\vartheta_1 &= \frac{m_2 \cos(\vartheta_1 - \vartheta_2) (l_1 \sin(\vartheta_1 - \vartheta_2) \dot\vartheta_1^2 - g \sin(\vartheta_2)) + m_2 l_2 \sin(\vartheta_1 - \vartheta_2) \dot\vartheta_2^2 + (m_1 + m_2) g \sin(\vartheta_1)}{m_2 l_1 \cos^2(\vartheta_1 - \vartheta_2) - (m_1+m_2) l_1} \\ \ddot\vartheta_2 &= \frac{m_2 l_2 \cos(\vartheta_1 - \vartheta_2) \sin(\vartheta_1 - \vartheta_2) \dot\vartheta_2^2 + (m_1+m_2) l_1 \sin(\vartheta_1 - \vartheta_2) \dot\vartheta_1^2 + (m_1+m_2) g \cos(\vartheta_1 - \vartheta_2) \sin(\vartheta_1) - (m_1+m_2) g \sin(\vartheta_2)}{(m_1+m_2) l_2 - m_2 l_2 \cos^2(\vartheta_1 - \vartheta_2)} \end{align*}

Diese Gleichungen sind durchaus sehr unhandlich und können nicht analytisch, gelöst werden — aber numerisch ist es kein Problem.

Seltsamer Attraktor

Zuvor habe ich bereits den Schmetterlingseffekt erwähnt. Um den Zusammenhang mit Chaos zu zeigen, betrachten wir folgendes Video von der Projektion in die y-z-Ebene von 13 Teilchen, die den Attraktor durchlaufen.

Alle Teilchen starten auf fast dem selben Punkt, aber nehmen sehr verschiedene Wege. Nach kurzer Zeit kann man den einzelnen Teilchen nicht mehr ansehen, dass sie fast die gleichen Anfangsbedingungen hatten.

Lorenz war Meteorologe und sein Differentialgleichungssystem

\begin{align} \dot{X} &= a(Y - X) \\ \dot{Y} &= X(b - Z) - Y \\ \dot{Z} &= XY - cZ, \\ \end{align}

das dieses chaotische Verhalten zeigt, sollte die Atmosphäre modellieren.

Jetzt kann man verstehen, was es mit dem Schmetterling aus Jurassic Park auf sich hat.

Er bewegt in Peking die Flügel, und im Central Park gibt’s Regen statt Sonne.

Dr. Ian Malcolm (1993)

Sein Flügelschlag ändert den Zustand eines chaotischen Systems, dem Wetter, ein wenig und nach einiger Zeit hat das System einen grundlegend anderen Weg eingeschlagen, als ohne diesen Flügelschlag.

Dennoch sieht das Video irgendwie geordnet aus. Fast schon vorhersagbar. Seltsam.