Problem with HamiltonianProblem in Differential

Specific Domains Modelling & Simulations
,
1

Hi there!

Im playing around with some easy physics problems. Here I tried to simulate a simple pendulum. I wanted to solve it via the HamiltonianProblem functionality from DifferentialEquations.jl.
Im not sure what’s the problem with the code below…
Thanks for help and advice!

using DifferentialEquations, DiffEqPhysics

struct Pendulum
    m
    g
    l
    ϕ₀
    dϕ₀
end

function (pen::Pendulum)(q, p) #hamiltonian of a pendulum 
    p^2 / (2 * pen.m * (pen.l)^2) - pen.m * pen.g * pen.l * cos(q)
end

function startvalues(p::Pendulum)
    q₀ = p.ϕ₀
    p₀ = p.dϕ₀ * p.m * (p.l)^2
    q₀, p₀
end

function simulate(p::Pendulum, tspan)
    q₀, p₀ = startvalues(p)
    prob = HamiltonianProblem(p, p₀, q₀, tspan)
    sol = solve(prob, Tsit5())

    sol
end

p = Pendulum(1.0, 9.81, 1.0, pi/2, 0.0)
tspan = (0.0, 1.0)
simulate(p, tspan)
2

That’s not one of the allowed forms.

Define a physical system by its Hamiltonian function `H(p, q, param)` or the function pair
`dp = -∂H/∂q` and `dq = ∂H/∂p`.

You’re missing param. The following works:

using DifferentialEquations, DiffEqPhysics

struct Pendulum
    m
    g
    l
    ϕ₀
    dϕ₀
end

function (pen::Pendulum)(p, q, param) #hamiltonian of a pendulum
    p^2 / (2 * pen.m * (pen.l)^2) - pen.m * pen.g * pen.l * cos(q)
end

function startvalues(p::Pendulum)
    q₀ = p.ϕ₀
    p₀ = p.dϕ₀ * p.m * (p.l)^2
    q₀, p₀
end

function simulate(p::Pendulum, tspan)
    q₀, p₀ = startvalues(p)
    prob = HamiltonianProblem(p, p₀, q₀, tspan)
    sol = solve(prob, Tsit5())

    sol
end

p = Pendulum(1.0, 9.81, 1.0, pi/2, 0.0)
tspan = (0.0, 1.0)
simulate(p, tspan)

(also notice the fix to the p,q order)

Note that this now has a better error message:

ERROR: All methods for the Hamiltonian `H` had too few arguments. A
Hamiltonian `H` must define either `H(p, q, param)` or `H(p, q, param, t)`. This error
can be thrown if you define an Hamiltonian for example as `H(p, q)`.
Note that `param` must be in the arguments list even if it's not used.
For more information on the required number of arguments for the function
you were defining, consult the documentation for the `HamiltonianProblem`.

Offending function: f
Methods:
# 1 method for callable object:
[1] (pen::Pendulum)(q, p) in Main at d:\OneDrive\Computer\Desktop\test.jl:120

Stacktrace:
 [1] HamiltonianProblem{false}(H::Pendulum, p0::Float64, q0::Float64, tspan::Tuple{Float64, Float64}, param::Nothing; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ DiffEqPhysics C:\Users\accou\.julia\dev\DiffEqPhysics\src\hamiltonian.jl:122
 [2] HamiltonianProblem
   @ C:\Users\accou\.julia\dev\DiffEqPhysics\src\hamiltonian.jl:115 [inlined]
 [3] #HamiltonianProblem#1
   @ C:\Users\accou\.julia\dev\DiffEqPhysics\src\hamiltonian.jl:95 [inlined]
 [4] HamiltonianProblem (repeats 2 times)
   @ C:\Users\accou\.julia\dev\DiffEqPhysics\src\hamiltonian.jl:93 [inlined]
 [5] simulate(p::Pendulum, tspan::Tuple{Float64, Float64})
   @ Main d:\OneDrive\Computer\Desktop\test.jl:132
 [6] top-level scope
   @ d:\OneDrive\Computer\Desktop\test.jl:140
3

Thank you for the quick response Chris! Much appreciate your work! Nice that you also changed the error message, with that will be nice for someone in the future who does the same mistake as I did! The error message I received was a bit cryptic :smiley:

1 Like