0

I am trying to plot the poincaré surface sections for the Henon- Heiles Hamiltonian using the symplectic integrator Ruth of 3rd order, but I am having some problems because after several steps the variables start to take huge values resulting in an error. I do not understand why this is happening as the integrator I am using is supposed to be symplectic.

I have tried to write this code several times but I always get some kind of error, so I would appreciate a lot if someone can explain me what I am doing wrong.

The code I wrote in python looks like this:

import math as m
import numpy as np
from matplotlib import pyplot as plt

def xpunto(vx0): #time derivative of x
    return vx0
def pxpunto(x0,y0,lambda0): #time derivative of px
    return -x0-(2.*lambda0*x0*y0)
def ypunto(vy0): #time derivative of y
    return vy0
def pypunto(x0,y0,lambda0): #time derivative of py
    return -y0-(lambda0*((x0*x0)-(y0*y0)))
def energia_poinc(y0,py0, H0): #initial condition for px using the constraints
    return m.sqrt((2.*H0)-(py0*py0)-(y0*y0)+((2./3.)*(y0*y0*y0)))


iteraciones=100000
dt=0.01
t=[0.]
x=[]
y=[]
px=[]
py=[]

def ruth3rd(xx,yy,ppx,ppy,lambda00):#define function with symplectic integrator

    aa=[1,-1./24.,3./4.,7./24.] #constants for position variable of the integrator
    bb=[1,1.,-2./3., 2./3.] #constants of velocity of the integrator

    x.append(xx)
    y.append(yy)
    px.append(ppx)
    py.append(ppy)



    xi , yi , pxi , pyi= xx, yy, ppx , ppy
    for i in range (iteraciones):

        vvx = pxi+bb[1]*pxpunto(xi,yi,lambda00)*dt
        VV1y = pyi+bb[1]*pypunto(xi,yi,lambda00)*dt
        xx = xi + aa[1]*vvx*dt
        yy =yi + aa[1]*VV1y*dt
        for j in range(2,len(aa)):
            vvx += bb[j]*pxpunto(xx,yy,lambda00)*dt
        VV1y += bb[j]*pxpunto(xx,yy,lambda00)*dt
            xx += aa[j]*vvx*dt
            yy += aa[j]*VV1y*dt

        xi = xx
        pxi = vvx
        yi = yy
        pyi = VV1y

        x.append(xi)
        px.append(pxi) 
        y.append(yy)
        py.append(pyi) 
        t.append(t[-1] + dt)

setsrange=10
sets=np.array(np.linspace(-0.4, 0.4, 10),dtype=np.float64)


for I in range(1):

    xpoinc=[]
    ppoinc=[]

    for ii in range (1):
        for jj in range (setsrange):

            Hinicial=1./12.
            x=[0.] 
            y=[0.4]
            py=[sets[jj]] 
            print (py[0],y[0],x[0],Hinicial)
            px=[energia_poinc(y[0],py[0],Hinicial)] 

            ruth3rd(x[0],y[0],px[0],py[0],2.)       


        crossings=[]
            crossingsindex=[]    
            for idx, item in enumerate(x[:-1]): #intersections
                if x[idx] < 0 and x[idx+1] > 0:
                    crossings.append(x[idx])
                    crossingsindex.append(idx)


            for index in crossingsindex:
                xpoinc.append(y[index])
                ppoinc.append(py[index])
Dean
  • 1
  • 1
  • Your function `energia_poinc` is wrong, you need the square root of the whole expression. Also, the missing terms are only zero for $x0=0$, but then you would not need that as parameter. – Lutz Lehmann Apr 07 '19 at 08:23
  • Thank you for pointing that out, I fixed it on the code. Now I get a ValueError: math domain error on that `energia_poinc` function on the first iteration but I cant figure out why because the numbers I am using are not that small. – Dean Apr 07 '19 at 13:53
  • Did you check that the energy is large enough that the radicant is positive? – Lutz Lehmann Apr 07 '19 at 14:41
  • For `x0, y0, py0 = 0,0.4,0.2` I start getting deviations in the energy from `t=3.6..4.0` on. Note that some small deviation is normal, as the conserved quantity of the numerical method is a perturbation of the Hamiltonian. And that the Hamiltonian due to the last term does not have bounded level sets, so that a divergence to infinity can be a correct behavior of a solution. – Lutz Lehmann Apr 07 '19 at 14:56

0 Answers0