# Integrating Hamilton's equations

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])


You need to write

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


to overwrite the input vectors instead of creating new ones.

1 Like

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

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


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

2 Likes