# How to define and integrate a HamiltonianProblem?

**URL:** <https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303>\
**Category:** General Usage\
**Tags:** package, diffeq, forwarddiff\
**Created:** [May 19, 2022, 9:07am UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303 "2022-05-19T09:07:55Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [May 19, 2022, 9:07am UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303/1 "2022-05-19T09:07:55Z")

</div>

Dear friends. I have a hamiltonian for the cart and pole system which is a standard problem in RL. [Lagrangian version shown here](https://www-robotics.cs.umass.edu/~grupen/503/SLIDES/cart-pole.pdf), I just legendre transformed it.

I don’t need a full solution. I just want to step it forward in time in a symplectic, energy preserving way. Here’s what I got so far

```julia
using DiffEqPhysics
using DifferentialEquations

function cart_hamiltonian(q, p; 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

a = HamiltonianProblem(cart_hamiltonian,[0,0.3],[0,.1],[0,.2]) # gives me ODE problem?
sol = init(a, Tsit5())

```

which yields

`MethodError: no method matching cart_hamiltonian(::Vector{Float64}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqPhysics.PhysicsTag, Float64}, Float64, 2}}, ::Nothing)`

when cart\_hamiltonian sees dual numbers.

ultimately I just want a function

```julia
function step(p,q,δt)
return (p_new, q_new) #new p q at slightly later time

```

I could write an euler integrator, but I feel like I’m close. I see `step!` in the documentation. Part of the problem with Julia is that it’s so easy to write blazing fast integrators I never needed to learn the details. 😍 But autodiffed hamiltonians would be so freaking awesome. Thanks

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [May 19, 2022, 11:05pm UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303/2 "2022-05-19T23:05:24Z")

</div>

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

```julia
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`](https://diffeq.sciml.ai/stable/basics/integrator/#Initialization-and-Stepping) like this

```julia
integrator = initialize()
step!(integrator, dt, true)

```

---

<div class="post-metadata">

**Author:** ![genkuroki](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/genkuroki/32/18030_2.png) [@genkuroki](https://discourse.julialang.org/u/genkuroki)\
**Post date:** [May 19, 2022, 11:51pm UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303/3 "2022-05-19T23:51:40Z")

</div>

Working Example: [https://github.com/genkuroki/public/blob/main/0030/DiffEqPhysics%20HamiltonianProblem.ipynb](https://github.com/genkuroki/public/blob/main/0030/DiffEqPhysics%20HamiltonianProblem.ipynb)

```julia
using DiffEqPhysics
using DifferentialEquations
using Plots
using StaticArrays

params = (g = 9.8, m1 = 2.0 , m2 = 1.0, L = 0.5)

function cart_hamiltonian(p, q, params, t)
    (; g, m1, m2, L) = params
    (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

p0 = SVector(0.0, 0.3)
q0 = SVector(0.0, 0.1)
tspan = (0.0, 0.2)
a = HamiltonianProblem(cart_hamiltonian, p0, q0, tspan, params)
sol = solve(a, Tsit5())
plot(sol; legend=:outertopright)

```

Output:

 ![2022-05-20](https://global.discourse-cdn.com/julialang/original/3X/8/5/857150b300ce6fd57b3015d2413dd96450e944fe.png)

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [May 20, 2022, 4:36pm UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303/4 "2022-05-20T16:36:27Z")

</div>

Thanks for the great replies. I think this is fantastic and I’m marking as solved. But now I think I should have disclosed my full use case:

I need really this function with the reinforcement learning form: update(state, action) → state

Now I see this “param” field, and it seems like a great place to put in the action (i.e. the control choice as an added linear potential.) I have

```julia
 
function cart_hamiltonian(p, q, control; 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)) 
        - control * x # this is the control, an added linear potential that kicks the system
    return h
end

```

So the following code works for my purposes, but I feel like it could be made more elegant than just rebuilding an integrator each time step. Can I change the parameter “on the fly” somehow with the same integrator?

```julia

function initialize(p,q,c)
    a = HamiltonianProblem(cart_hamiltonian,p,q,[0.0,1.0],c)
    return init(a, Tsit5())
end

function update(s,c; dt = .2)
   (t,p,q) = s
   integrator = initialize(p,q,c)
   step!(integrator, dt, true)
   (p1,q1) = integrator.u
   return (t+dt,p1,q1)
end

update((1.0,[0.1,0.2],[.3,.4]),.2) # the propagator in general control form

```

Thanks for the great help!

---

<div class="post-metadata">

**Author:** ![genkuroki](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/genkuroki/32/18030_2.png) [@genkuroki](https://discourse.julialang.org/u/genkuroki)\
**Post date:** [May 21, 2022, 12:28am UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303/5 "2022-05-21T00:28:25Z")

</div>

Second simple example: [https://github.com/genkuroki/public/blob/main/0030/example%20of%20init%20and%20step!%20of%20integrator.ipynb](https://github.com/genkuroki/public/blob/main/0030/example%20of%20init%20and%20step!%20of%20integrator.ipynb)

```julia
using DiffEqPhysics
using DifferentialEquations
using Plots
default(fmt = :png, size = (450, 250))
using StaticArrays

module My
mutable struct ParamSinglePendulum{Tg, Tl}
    g::Tg
    l::Tl
end
end

function singlependulum(p, q, param, t = nothing)
    (; g, l) = param
    dθ = p[1]
    θ = q[1]
    (1/2)*(dθ/l)^2 - g*l*cos(θ)
end

plot_singlependulum(sol) = plot(sol; label=["dθ/dt" "θ"], legend=:outertopright)

param = My.ParamSinglePendulum(9.8, 1.0) # on the earth
q0 = SVector(-0.9π)
p0 = SVector(0.0)
tspan = (0.0, 10.0)
prob = HamiltonianProblem(singlependulum, p0, q0, tspan, param)

```

First, a single pendulum is shaken on the earth.

```julia
integrator = init(prob, Tsit5())
param.g = 9.8 # on the earth
dt = 0.1
for _ in 1:40
    step!(integrator, dt, true)
end
plot_singlependulum(integrator.sol)

```

![1](https://global.discourse-cdn.com/julialang/original/3X/b/d/bd82f9195146c5e010b5c8ad2313b2ce1a9624c0.png)

Then, suddenly, the gravitational acceleration becomes equal to the gravitational acceleration of the moon. 😜

```julia
param.g = 1.6 # on the moon
dt = 0.1
for _ in 1:210
    step!(integrator, dt, true)
end
plot_singlependulum(integrator.sol)

```

![2](https://global.discourse-cdn.com/julialang/original/3X/9/7/977de36ee15236ad4fa0e73bcf24a89d37830676.png)

The gravitational acceleration returns to its original level.

```julia
param.g = 9.8 # on the earth
dt = 0.1
for _ in 1:200
    step!(integrator, dt, true)
end
plot_singlependulum(integrator.sol)

```

![3](https://global.discourse-cdn.com/julialang/original/3X/c/4/c4e9bc0b6d7755a2c15935066387a36e6503dfa5.png)

---

<div class="post-metadata">

**Author:** ![cdawg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdawg/32/9811_2.png) [@cdawg](https://discourse.julialang.org/u/cdawg)\
**Post date:** [May 21, 2022, 6:51pm UTC](https://discourse.julialang.org/t/how-to-define-and-integrate-a-hamiltonianproblem/81303/6 "2022-05-21T18:51:49Z")

</div>

very cool. @genkuroki. just read your post about the [**Problem-Algorithm-Solver pattern**](https://discourse.julialang.org/t/function-depending-on-the-global-variable-inside-module/64322/10) As someone who never uses globals, or modules outside of the parent source file in a package, I am amazed at your awareness and control of Julia scope.

However, I wonder if there is an irremovable disadvantage to your pattern in the case where you want to have multiple actors acting in parallel, each with different state, each changing this potential parameter 😅. I think I need this “update(state, action) → state” mapping as a closed, stateless function.

Sadly for my use case, the transparency of a hard coded Euler solver may ultimately be unbeatable here. 😶 But I still learned a lot and hope this post will be helpful for HamiltonianProblem users in the future. Thanks very much!