Integrating Hamilton's equations

Specific Domains Numerics
,
1

I want to integrate the equations of motion of a classical physical system.
My Hamiltonian is changing in time, and actually, I am interested in considering the general case where it is not necessarily of the form H(q, p) = \frac{p^2}{2m} + U(q).

My goal is to get high accuracy in the solution, so I guess that using a symplectic integrator is the way-to-go. (Or maybe not? Note that since the Hamiltonian is changing in time, energy is not conserved, but still the symplectic structure holds).

From the documentation of DifferentialEquations.jl I could not understand what is a good approach.

It seems to me like the only options to use symplectic integrators in my case (time changing Hamiltonian, and non-standard form of q, p) are SecondOrderODEProblem or DynamicalODEFunction. However, I could not find an example that showcases DynamicalODEFunction. Is there such an example?

I tried to use SecondOrderODEProblem from the example here, and to modify from the two coordinates x,y to just one coordinate x. However, I get that the velocity does not change in time. I am surely doing something wrong. I would love to get some help.
Here is my modification of the example:

function HH_acceleration!(dv,v,u,p,t)
    x = u
    dx = dv
    dv = -cos.(x)
end
initial_positions = [0.1]
initial_velocities = [0.5]
prob = SecondOrderODEProblem(HH_acceleration!,initial_velocities,initial_positions,(0.0, 10.0))
sol2 = solve(prob, KahanLi8(), dt=1/10)

And what I get back is:

retcode: Success
Interpolation: 3rd order Hermite
t: 101-element Vector{Float64}:
  0.0
  0.1
  0.2
  0.30000000000000004
  0.4
  0.5
  0.6
  0.7
  0.7999999999999999
  0.8999999999999999
  0.9999999999999999
  1.0999999999999999
  1.2
  ⋮
  8.899999999999984
  8.999999999999984
  9.099999999999984
  9.199999999999983
  9.299999999999983
  9.399999999999983
  9.499999999999982
  9.599999999999982
  9.699999999999982
  9.799999999999981
  9.89999999999998
 10.0
u: 101-element Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
 ([0.5], [0.1])
 ([0.5], [0.15])
 ([0.5], [0.20000000000000004])
 ([0.5], [0.25000000000000006])
 ([0.5], [0.3000000000000001])
 ([0.5], [0.35000000000000014])
 ([0.5], [0.4000000000000002])
 ([0.5], [0.45000000000000023])
 ([0.5], [0.5000000000000002])
 ([0.5], [0.5500000000000003])
 ([0.5], [0.6000000000000003])
 ([0.5], [0.6500000000000004])
 ([0.5], [0.7000000000000004])
 ⋮
 ([0.5], [4.550000000000033])
 ([0.5], [4.600000000000033])
 ([0.5], [4.650000000000034])
 ([0.5], [4.700000000000035])
 ([0.5], [4.7500000000000355])
 ([0.5], [4.800000000000036])
 ([0.5], [4.850000000000037])
 ([0.5], [4.900000000000038])
 ([0.5], [4.950000000000038])
 ([0.5], [5.000000000000039])
 ([0.5], [5.05000000000004])
 ([0.5], [5.1000000000000405])
2

You need to write

dx .= dv
dv .= -cos.(x)

to overwrite the input vectors instead of creating new ones.

1 Like
3

The second order ODE equation is defined as follows:

function eqn!(ddu, du, u, p, t)
    ddu[1]=-cos(u[1])
end

function getsol()
    #du0=[0.5]
    #u0=[0.1]
    #tspan=(0, 10.0)
    probl=SecondOrderODEProblem(eqn!, [0.5], [0.1], (0.0,10.0))
    solve(probl, KahanLi8(), dt=0.05)
end
s=getsol()
2 Likes
4

Thanks! Actually, this is what is working for me

function HH_acceleration!(dv,v,u,p,t)
    x = u
    dx = dv
    dv .= -cos.(x)
end
5

Oh, yes, you are right! You can actually remove dx = dv since it doesn’t do anything (for second order problems) :wink:

2 Likes