Figure 2.15:

Contours of an energy function V(x_{1},x_{2}) or H(x_{1},x_{2}).

Code for Figure 2.15

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

n = 201
x = np.linspace(-5, 5, n)
y = np.linspace(-5, 5, n)
kappa = 2

xx, yy = np.meshgrid(x, y)

xx2 = (xx - 2) ** 2
yy2 = yy ** 2

z = 10 * yy2 - kappa * np.cos(xx * np.pi / 2) + 0.2 * xx2

val = [-1, 0, 1, 2, 3, 5, 8, 10]
contours = plt.contour(x, y, z, levels=val)
plt.axis([-1.75, 3, -1, 1])
plt.xlabel("$q$")
plt.text(-5.75, 0, "$p$", rotation=0)

# Write contour data for output
with open('energy.dat', 'w') as myfile:
    for i, level_segs in enumerate(contours.allsegs):
        for j, vertices in enumerate(level_segs):
            header = f"Segment {i + 1}, Path {j + 1}"
            np.savetxt(myfile, vertices, fmt='%f', header=header)
            myfile.write('\n')

def rhs(t, w):
    x, y, a = w
    dxdt = -np.pi * np.sin(np.pi/2 * x) - 0.4 * (x - 2)
    dydt = -20 * y
    dadt = -(dxdt ** 2 + dydt ** 2)
    return [dxdt, dydt, dadt]

w0 = [1.6, 0.36666, 3]
npts = 101
time = np.linspace(0, 0.5, npts)
sol = solve_ivp(rhs, [time[0], time[-1]], w0, t_eval=time)

plt.plot(sol.y[0], sol.y[1], '-o')

plt.show(block=False)