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, seg in enumerate(contours.collections):
for j, path in enumerate(seg.get_paths()):
header = f"Segment {i + 1}, Path {j + 1}"
np.savetxt(myfile, path.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)