area-preserving Hénon map

\[\begin{split}p_{n+1} &= p_{n} - V'(q_{n}) \\ q_{n+1} &= q_{n} + p_{n}\end{split}\]

where \(V(q) = cq + q^3/3\)

import numpy as np
import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings('ignore')


def fmap(q,p,c):
    p = p - (c - q**2)
    q = q + p
    return q, p

tmax = 1000
sample=1000
c = 0.8
q = np.linspace(-5,5,sample)
p = np.array([0.0]*sample)
traj = [np.array([])]*2
for i in range(tmax):
    q,p = fmap(q,p,c)


    traj[0] = np.append(traj[0], q)
    traj[1] = np.append(traj[1], p)

fig, ax = plt.subplots(1,1)
ax.plot(traj[0], traj[1], ",k")
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
plt.show()