import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
npts = 101
time = np.linspace(0, 10, npts)
z0 = np.array([0, 0.05])
def rhs(t, z):
x = z[0]
y = z[1]
zdot = [x - y - x*(x**2+2*y**2), x + y - y*(x**2+y**2)]
return zdot
sol = solve_ivp(rhs, [time[0], time[-1]], z0, t_eval=time)
z0 = sol.y[:, -1]
time = np.linspace(0, 6.5, npts)
sol_cyc = solve_ivp(rhs, [time[0], time[-1]], z0, t_eval=time)
theta = np.linspace(0, 2*np.pi, npts)
anni = np.column_stack((np.cos(theta), np.sin(theta)))
anno = 2/np.sqrt(5) * anni
traj = np.column_stack((sol.y.transpose(), sol_cyc.y.transpose(), anni, anno))
plt.figure()
plt.plot(time, sol.y.transpose())
plt.figure()
plt.plot(sol.y[0], sol.y[1], '-', sol_cyc.y[0], sol_cyc.y[1], 'o', anni[:, 0], anni[:, 1], anno[:, 0], anno[:, 1])
plt.show(block=False)
with open("cycle.dat", "w") as f:
np.savetxt(f, traj, fmt='%f')