henon-heiles.ipynb

Kernel: Unknown Kernel

Solution to Henon-Heiles system

#Importing libraries
import numpy as np
%pylab inline
import matplotlib.pylab as plt

Populating the interactive namespace from numpy and matplotlib

#Runge-Kutta 4 integrator
def RK4_step( f, y, r, h ):
    #Creating solutions
    K0 = h*f(r, y)
    K1 = h*f(r + 0.5*h, y + 0.5*K0)
    K2 = h*f(r + 0.5*h, y + 0.5*K1)
    K3 = h*f(r + h, y + K2)
    y1 = y + 1/6.0*(K0 + 2.0*K1 + 2.0*K2 + K3 )
    #Returning solution
    return y1

#============================================================
#PARAMETERS
#============================================================
#Lambda
lamb = 1.0
#Xmax
Xmax = 0
#Ymax
Ymax = 0.5
#Pymax
Pymax = 0.5

#Tolerance of Poincare section
eps = 0.001

#Dynamic function of the system
def f( t, y ):
    #Extracting values
    X0 = y[0]
    Px0 = y[1]
    Y0 = y[2]
    Py0 = y[3]
    #Derivatives
    dx = Px0
    dpx = -X0 - 2*lamb*X0*Y0
    dy = Py0
    dpy = -Y0 - lamb*( X0**2 - Y0**2 )
    
    return np.array([dx,dpx,dy,dpy])

#Potential
def V(x,y):
    return 0.5*(x**2+y**2) + lamb*(x**2*y - y**3/3.0)

#Energy
def Energy(x,px,y,py):
    return 0.5*(px**2+py**2) + V(x,y)

#Plot of potential
#Xmax
Xmaxp = 2.0
#Ymax
Ymaxp = 2.0
#Meshgrids X
xar = np.linspace(-Xmaxp,Xmaxp,100)
yar = np.linspace(-Ymaxp,Ymaxp,100)
X, Y = np.meshgrid( xar, yar )
Z = V(X,Y)
plt.figure( figsize=(7,7) )
plt.contourf(Z,40,extent = (-Xmaxp,Xmaxp,-Ymaxp,Ymaxp))
plt.xlabel( "x", fontsize=14 )
plt.ylabel( "y", fontsize=14 )
plt.title("Contours of Henon-Heiles potential", fontsize=14)

<matplotlib.text.Text at 0x7fa8b9f54e90>

#Poincare diagram
def Poincare( energy, Ncond, tmax, h ):
    #Specific solution
    ys = []
    pys = []
    #Finding initial conditions
    for n in xrange(Ncond):
        #Finding initial conditions
        while True:
            x0 = -Xmax + 2*np.random.random()*Xmax
            y0 = -Ymax + 2*np.random.random()*Ymax
            py0 = -Pxmax + 2*np.random.random()*Pxmax
            px02 = 2*energy - py0**2 - (x0**2+y0**2) - 2*lamb*(x0**2*y0 - y0**3/3.0)
            if px02>=0:
                px0 = -np.sqrt(px02)
                break
        #Solutions
        X = [x0,]
        PX = [px0,]
        Y = [y0,]
        PY = [py0,]
        #Integration
        i = 0
        for t in np.arange(0,tmax,h):
            y = [ X[i], PX[i], Y[i], PY[i] ]
            Xi, Pxi, Yi, Pyi = RK4_step( f, y, t, h )
            #print Energy(Xi,Pxi,Yi,Pyi)
            #Detecting a passing orbit through plane Y-Py
            if abs(Xi)<eps:
                ys.append( Yi )
                pys.append( Pyi )
            #Adding new solutions
            X.append(Xi)
            PX.append(Pxi)
            Y.append(Yi)
            PY.append(Pyi)
            i += 1
    plt.figure( figsize=(6,6) )
    plt.plot( ys, pys, ".", ms=2.5 )
    plt.title( "Poincare Map for $e =$%1.3f"%(energy), fontsize=14 )
    plt.grid(1)
    plt.xlabel("$y$", fontsize=14)
    plt.ylabel("$dy/dt$", fontsize=14)
    return None

Poincare( 1/100., 60, 400.0, 0.01 )

Poincare( 1/12., 100, 400.0, 0.01 )

Poincare( 1/10., 100, 400.0, 0.01 )

Poincare( 1/8., 100, 400.0, 0.01 )

Poincare( 1/6., 100, 400.0, 0.01 )