How to define and integrate a HamiltonianProblem?

I think you have p and q reversed in the definition of H. The documentation is a little unclear on this, but H actually takes 3 arguments H(p, q, params). The params argument defaults to nothing, but can be useful if you want to run multiple integrations of the same equations with different masses for example.

I think you need

function cart_hamiltonian(p, q, params; g = 9.8, m1 = 2.0 , m2 = 1.0, L = 0.5) 
    (x, θ) = q
    (p_x, p_θ) = p
    h =  g * L * m2 * cos(θ) + 
        (L^2 * m2 * p_x^2 + (m1 + m2)* p_θ^2 - 
        2 *L * m2 * p_x * p_θ * cos(θ) ) / 
        ( 2 * L^2*m2*(m1 + m2 * sin(θ)^2))
    return h
end

function initialize()
    a = HamiltonianProblem(cart_hamiltonian,[0,0.3],[0,.1],[0,.2])
    integrator = init(a, Tsit5())
end

Then you can step like this

integrator = initialize()
step!(integrator, dt, true)
2 Likes