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
In [4]:
def euler_step(f, dt, tn, xn):
return xn + dt * f(tn, xn)
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
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
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()
In [10]:
def poincare_plot(px, x, py, y)