HamiltonianProblem gives inconsistent results

[Vlmbd]

Hello,
I'm trying to integrate this Hamiltonian

H=(p1^2+p2^2)/(2mu)+V(exp(-aq1)-1)^2+V(exp(-aq2)-1)^2-p1p2/m2

but it looks like HamiltonianProblem is giving me inconsistent results

m1,m2,a,V=1.,1.,1.,1.
mu=m1*m2/(m1+m2)
param = [a,V,mu,m2]
tspan = (0., 1000.)
jul_H = Main.eval("""
function H(p,q,param)
    p1, p2 = p
    q1, q2 = q
    a, V, mu, m2 = param
    (p1^2+p2^2)/(2*mu)+V*(exp(-a*q1)-1)^2+V*(exp(-a*q2)-1)^2-p1*p2/m2
    end""")
q1,p1,q2,p2,=0.5,0.,0.,0.    
p0=[p1,p2]    
q0=[q1,q2]
prob=de.HamiltonianProblem(jul_H,p0,q0,tspan,param)
sol = de.solve(prob,de.VelocityVerlet(),dt=0.01)

For example in this case with the VelocityVerlet integrator p2 and q2 remain zero, but the other integrators don't seem to work either.

On the other hand, by applying the Hamilton equations and integrating with ODEProblem, I get good results

jul_f = Main.eval("""
function f(dX,X,p,t)
  q1, p1, q2, p2 = X
  a, V, mu, m2 = p
  dX[1] = p1/mu - p2/m2
  dX[2] = 2*V*(exp(-a*q1)-1)*a*exp(-a*q1)
  dX[3] = p2/mu - p1/m2
  dX[4] = 2*V*(exp(-a*q2)-1)*a*exp(-a*q2)
end""")

X0=[q1,p1,q2,p2]
prob = de.ODEProblem(jul_f, X0, tspan, param)
sol = de.solve(prob,de.ImplicitMidpoint(),dt=0.01)

I would like to make HamiltonianProblem work because it allows to use several symplectic integrators!

Thanks,
Valentin

[ChrisRackauckas]

Hey,
Thanks for the report. It looks like it was just a typo in the docs. It should be H(q,p,params) instead of H(p,q,params). Flipping that to

jul_H = Main.eval("""
function H(q,p,param)
    p1, p2 = p
    q1, q2 = q
    a, V, mu, m2 = param
    (p1^2+p2^2)/(2*mu)+V*(exp(-a*q1)-1)^2+V*(exp(-a*q2)-1)^2-p1*p2/m2
    end""")

fixed it. The docs have now been updated.