#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.

  • 1
    $\begingroup$ 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. $\endgroup$ – zap Mar 14 '17 at 18:44
  • $\begingroup$ 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? $\endgroup$ – Prakhar Godara Mar 14 '17 at 21:35
  • 2
    $\begingroup$ Can you write your equations of motion? Also, which symplectic integrator are you using? $\endgroup$ – nicoguaro Mar 15 '17 at 16:44
  • $\begingroup$ 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 $\endgroup$ – Prakhar Godara Mar 15 '17 at 23:20
  • $\begingroup$ 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. $\endgroup$ – zap Mar 16 '17 at 18:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.