Computational Physics – Übungsblatt 10

Hausaufgabe 16: Poincaréschnitte des Hénon-Heiles-Systems

Das Potential des Hénon-Heiles-Systems ist \begin{equation} V(x, y) = \frac{1}{2}\left(x^2 + y^2\right) + x^2y - \frac{1}{3}y^3 \end{equation} wobei $x$ und $y$ die räumlichen 2D-Koordinaten sind.

Bevor das Potential mit Hilfe von Poincaré-Schnitten untersucht wird, werden zunächst verschiedene Verfahren zur Integration von Differentialgleichungen implementiert.


In [1]:
%matplotlib inline

In [2]:
import numpy as np
import matplotlib.pyplot as plt

In [3]:
def integrate(step_method, f, N, dt, t0, x0):
    ts = []
    xs = []
    t = t0
    x = x0
    for _ in range(0, N):
        ts.append(t)
        xs.append(x)
        x = step_method(f, dt, t, x)
        t += dt
    return ts, xs

Eulerverfahren


In [4]:
def euler_step(f, dt, tn, xn):
    return xn + dt * f(tn, xn)

Runge-Kutta 4. Ordnung


In [5]:
def runge_kutta_4_step(f, dt, tn, xn):
    k1 = f(tn, xn) * dt
    k2 = dt * f(tn + dt / 2, xn + k1 / 2)
    k3 = dt * f(tn + dt / 2, xn + k2 / 2)
    k4 = dt * f(tn + dt, xn + k3)
    return 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) + xn

Symplektischer Algorithmus 4. Ordnung


In [6]:
def symplectic_4_step(f, dt, tn, xn):
    c1 = c4 = 1 / (2 * (2 - 2 ** 1/3))
    c2 = c3 = (1 - 2 ** 1/3) / (2 * (2 - 2 ** 1/3))
    d1 = d3 = 1 / (2 - 2 ** 1/3)
    d2 = - 2 ** 1/3 / (2 - 2 ** 1/3)
    d4 = 0

Harmonischer Oszillator


In [7]:
def f_harmosz_2d(t, x, k1=1, k2=5):
    return np.matrix(
        [[0, -1/k1, 0, 0],
         [1, 0, 0, 0],
         [0, 0, 0, -1/k2],
         [0, 0, 1, 0]]
    ) * x

In [8]:
N = 10000
dt = .01
t0 = 0
x0 = np.matrix([[0], [1], [0], [1]])

for method in (euler_step, runge_kutta_4_step):
    ts, xs = integrate(method, f_harmosz_2d, N, dt, t0, x0)
    px = [z[0, 0] for z in xs]
    x = [z[1, 0] for z in xs]
    py = [z[2, 0] for z in xs]
    y = [z[3, 0] for z in xs]
    #true_p = np.sin(ts)
    #true_q = np.cos(ts)
    # Make a nice plot
    plt.plot(ts, px, label='$p_x$')
    plt.plot(ts, x, label='$x$')
    plt.plot(ts, py, label='$p_y$')
    plt.plot(ts, y, label='$y$')
    #plt.plot(ts, true_x, label='x true')
    #plt.plot(ts, true_v, label='v true')
    plt.title('{}'.format(method.__name__.replace('_step', '')))
    plt.legend(bbox_to_anchor=(1.2, 1))
    plt.show()



In [9]:
def f_henon_2d(t, z):
    px = z[0, 0]
    x = z[1, 0]
    py = z[2, 0]
    y = z[3, 0]
    return np.matrix(
        [[-x - 2 * x * y],
         [px],
         [-y - x ** 2 + y ** 2],
         [py]]
    )

In [10]:
N = 1000
dt = .01
t0 = 0
x0 = np.matrix([[0], [1], [0], [1]])

for method in (euler_step, runge_kutta_4_step):
    ts, xs = integrate(method, f_henon_2d, N, dt, t0, x0)
    ps = [x[0, 0] for x in xs]
    qs = [x[1, 0] for x in xs]
    #true_p = np.sin(ts)
    #true_q = np.cos(ts)
    # Make a nice plot
    plt.plot(ts, ps, label='p')
    plt.plot(ts, qs, label='q')
    #plt.plot(ts, true_x, label='x true')
    #plt.plot(ts, true_v, label='v true')
    plt.title('{}'.format(method.__name__.replace('_step', '')))
    plt.legend(loc='best')
    plt.show()


-c:7: RuntimeWarning: overflow encountered in double_scalars
-c:9: RuntimeWarning: overflow encountered in double_scalars
-c:9: RuntimeWarning: invalid value encountered in double_scalars

In [10]:
def poincare_plot(px, x, py, y)