# 4th order Symplectic integration of a 2 degree of freedom Hamiltonian system

#f,g are the dH/dq and -dH/dp respectively

def f(v1,v2):

alpha = (a*v1 - b*v2)/k
beta =  (a*v2 - c*v1)/k

dot1 = (c-a)/(3*k) +  (a*np.exp(alpha)/k - c*np.exp(beta)/k)/(1 + np.exp(alpha) + np.exp(beta))
dot2 = (b-a)/(3*k) +  (a*np.exp(beta)/k - b*np.exp(alpha)/k)/(1 + np.exp(alpha) + np.exp(beta))
return dot1, dot2

def g(u1,u2):

dot1 = -1/3 + np.exp(u1)/(1 + np.exp(u1)+ np.exp(u2))
dot2 = -1/3 + np.exp(u2)/(1 + np.exp(u1)+ np.exp(u2))
return -dot1, -dot2

#these are the constants for 4th order symplectic rk method where c1 is for canonical position and c2 for canonical momentum cordinates. (As given by Ruth and Forest)

c1 = np.array((1.0/(2.0*(2.0-2.0**(1.0/3.0))), (1.0-2.0**(1.0/3.0))/(2.0*(2.0-2.0**(1.0/3.0))),(1.0-2.0**(1.0/3.0))/(2.0*(2.0-2.0**(1.0/3.0))), 1.0/(2.0*(2.0-2.0**(1.0/3.0)))))

c2 = np.array((1.0/(2.0-2.0**(1.0/3.0)), -2.0**(1.0/3.0)/(2.0-2.0**(1.0/3.0)),1.0/(2.0-2.0**(1.0/3.0)), 0.0))

#this is the function where the integration is happening and returns the 4 arrays which are time series of the phase-space variables

def sol(x0):

u1 = np.zeros(n+1)
u2 = np.zeros(n+1)
v1 = np.zeros(n+1)
v2 = np.zeros(n+1)

u1[0] = x0[0]
u2[0] = x0[1]
v1[0] = x0[2]
v2[0] = x0[3]
for k in xrange (n):
u10 = u1[k]
u20 = u2[k]
v10 = v1[k]
v20 = v2[k]
for i in xrange(len(c2)):
U1 = u10 + c1[i]*h*f(v10,v20)[0]
U2 = u20 + c1[i]*h*f(v10,v20)[1]
V1 = v10 + c2[i]*h*g(U1,U2)[0]
V2 = v20 + c2[i]*h*g(U1,U2)[1]
u10 = U1
u20 = U2
v10 = V1
v20 = V2

U1 = u10 + h*c1[-1]*f(v10,v20)[0]
U2 = u20 + h*c1[-1]*f(v10,v20)[1]
u10 = U1
u20 = U2
u1[k+1] = U1
u2[k+1] = U2
v1[k+1] = V1
v2[k+1] = V2

return u1,u2,v1,v2


When I run the code for some values of parameters (like a,b,c, etc) I get an overflow error in dot2. For other parameters I am getting a phase space which is not bounded, but a Hamiltonian system is supposed to have a bounded trajectory.

Is something wrong with my algorithm? Kindly point the error. I can provide the whole code and the paper, whose results I am trying to reproduce.

• A hamiltonian system does not necessarily have a bounded trajectory. Try plotting the contour of your hamiltonian using the parameter values where you get the overflow error to get a clearer picture of what is going on. – zap Mar 14 '17 at 18:44
• The Hamiltonian is autonomus. Hence it has to be bounded. Is there an error in the algorithm of symplectic Runge Kutta (4th order) code that I have written? – Prakhar Godara Mar 14 '17 at 21:35
• Can you write your equations of motion? Also, which symplectic integrator are you using? – nicoguaro Mar 15 '17 at 16:44
• H = H(u,v) and I have given the form of dH/du and dH/dv (a,b,c,k are constant parameters). The equations of motion are du/dt = dH/dv and dv/dt = -dH/du. Forest, E.; Ruth, Ronald D. (1990). "Fourth-order symplectic integration". Physica D. 43: 105. I am using this integrator – Prakhar Godara Mar 15 '17 at 23:20
• Take for instance the Hamiltonian of a free particle H=p^2/2. It is autonomus. It is bounded in any bounded domain of phase space, although it can take arbitrarily large values. The orbits however are straight lines, i.e. not bounded. – zap Mar 16 '17 at 18:25