mset of the Henon map

[1]A. Shudo and K. S. Ikeda, Complex Trajectory Description for Chaotic Tunneling, Progress of Theoretical Physics Supplement 139, 246 (2000).

plot mset henon
import numpy as np
import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings('ignore')

#np.seterr(invalid='ignore')

twopi = 2*np.pi

def get_Mset(evolve, pinit, tmax, xr=[0,twopi], yr =[-3,3],sample=500):
    qr = np.linspace(xr[0], xr[1], sample)
    qi = np.linspace(yr[0], yr[1], sample)
    if yr[0]*yr[1]<0:
        qi = np.linspace(yr[0], 0, sample//2,endpoint=False)
        qi = np.append(qi, np.linspace(0, yr[1], sample//2))

    QR,QI = np.meshgrid(qr,qi)
    q0 = QR + 1.j*QI
    p0 = np.zeros(q0.shape, dtype=np.complex128) + pinit
    x = np.array([np.copy(q0), np.copy(p0)])
    initial = np.copy(q0)

    x = evolves(x, tmax, evolve)
    index = (np.abs(x[0]) > 1e100) & (np.abs(x[1]) > 1e100)
    #x[0][index] = np.nan + 1.j*np.nan
    #x[1][index] = np.nan + 1.j*np.nan
    x[0][index] = np.inf
    x[1][index] = np.inf
    prod = x[1].imag*np.roll(x[1].imag,1)
    index = (prod > 0)
    return [initial, q0[~index], x]

def plot_contour(ax,levels, init, lset):
    #vmax = lset[1][~np.isnan(lset[1])].imag.max()
    ax.contour(init.real, init.imag, lset[1].imag, levels, colors='k',alpha=0.5,lw=1)



def fmap(x, c):
    q,p = x
    p = p - (c  - q**2)
    q = q + p #+ np.sin(p)
    return np.array([q, p])

def evolves(x, tmax, evolve):
    for i in range(tmax):
        x = evolve(x)
    return x

tmax = 6
c=0.8
evolve = lambda x: fmap(x, c)
init, mset, lset = get_Mset(evolve, 0, tmax, sample=2000, xr=[-2,2], yr=[-2,2])
fig, ax = plt.subplots(1,1)
ax.plot(mset.real, mset.imag, ',k')

levels = np.linspace(-2,2,20)
plot_contour(ax, levels, init, lset)

levels = [-np.inf, 0, np.inf]
colors = ["#00ff7f", "#00ffff"]
ax.contour(init.real, init.imag, lset[1].imag, levels, colors='k',linewidths=1)
ax.contourf(init.real, init.imag, lset[1].imag, levels, colors=colors,alpha=0.5, linewidths=3)

ax.set_xlabel(r"$\mathrm{Re}[x_0]$")
ax.set_ylabel(r"$\mathrm{Im}[x_0]$")
plt.show()

Total running time of the script: ( 0 minutes 1.905 seconds)

Gallery generated by Sphinx-Gallery