Hamiltonian Systems
Hamilton's equations of motion are given in terms of the Hamiltonian $H(q,p)$ by
\[\begin{align*} \frac{dq}{dt} &= \frac{\partial H}{\partial p} , & \frac{dp}{dt} &= - \frac{\partial H}{\partial q} . \end{align*}\]
In the following, we show how these equations can be obtained for the example of a harmonic oscillator.
Harmonic Oscillator
Before any use, we need to load EulerLagrange
:
using EulerLagrange
Next, we generate symbolic variables for a one-dimensional Hamiltonian system:
t, q, p = hamiltonian_variables(1)
(t, (q(t))[1:1], (p(t))[1:1])
We define a named tuple with typical values for the parameters, e.g.,
params = (k=0.5, ω=√0.5)
(k = 0.5, ω = 0.7071067811865476)
We use the function symbolize
to generate a symbolic version of the parameters:
sparams = symbolize(params)
(k = kₚ, ω = ωₚ)
Now we can define the Hamiltonian function:
using LinearAlgebra
H(t, q, p, params) = p ⋅ p / 2 + params.k * (q ⋅ q) / 2
H (generic function with 1 method)
The Hamiltonian, evaluated on and together with the symbolic variables and parameters is used to construct a HamiltonianSystem
:
ham_sys = HamiltonianSystem(H(t, q, p, sparams), t, q, p, sparams)
Hamiltonian system with
H = (1//2)*((p(t))[1]^2) + (1//2)*kₚ*((q(t))[1]^2)
and equations of motion
-(p(t))[1] + Differential(t)((q(t))[1])
Differential(t)((p(t))[1]) + kₚ*(q(t))[1]
The constructor computes Hamilton's equations and generates the corresponding Julia code. In the last step, we can now construct a HODEProblem
from the HamiltonianSystem
and some appropriate initial conditions, a time span to integrate over and a time step:
tspan = (0.0, 10.0)
tstep = 0.01
q₀, p₀ = [0.5], [0.0]
hprob = HODEProblem(ham_sys, tspan, tstep, q₀, p₀; parameters = params)
Geometric Equation Problem for Hamiltonian Ordinary Differential Equation (HODE)
with vector fields
v = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :Q, :P, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x468c1e7b, 0x820745e4, 0x03626fee, 0x51a7af52, 0x2f5ada10), Expr}(quote
#= none:1 =#
#= none:5 =#
begin
#= none:7 =#
#= none:7 =# @inbounds begin
#= none:9 =#
ˍ₋out[1] = getindex(P, 1)
#= none:11 =#
nothing
end
end
end)
f = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :t, :Q, :P, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x228e1638, 0x59733f35, 0xc90cfba5, 0xec3190e8, 0xfa0063d3), Expr}(quote
#= none:1 =#
#= none:5 =#
begin
#= none:7 =#
#= none:7 =# @inbounds begin
#= none:9 =#
ˍ₋out[1] = (-1 // 1 * getindex(Q, 1)) * params.k
#= none:11 =#
nothing
end
end
end)
Hamiltonian: H = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:t, :Q, :P, :params), EulerLagrange.var"#_RGF_ModTag", EulerLagrange.var"#_RGF_ModTag", (0x9beadfb6, 0x42cd9174, 0x39dfa402, 0xcc13b23c, 0x822ba1e9), Expr}(quote
#= none:1 =#
#= none:5 =#
1 // 2 * getindex(P, 1) ^ 2 + (1 // 2 * getindex(Q, 1) ^ 2) * params.k
end)
Invariants:
GeometricBase.NullInvariants()
Timespan: (0.0, 10.0)
Timestep: 0.01
Initial conditions:
(t = 0.0, q = [0.5], p = [0.0])
Parameters:
(k = 0.5, ω = 0.7071067811865476)
We can integrate this system using GeometricIntegrators:
using GeometricIntegrators
sol = integrate(hprob, Gauss(1))
using CairoMakie
fig = lines(parent(sol.q[:,1]), parent(sol.p[:,1]);
axis = (; xlabel = "q₁", ylabel = "p₁", title = "Harmonic Oscillator"),
figure = (; size = (800,600), fontsize = 22))