Julia has a great ODE package, providing a vast number of algorithms. The DiffEqPhys.jl package allows to use them for the numerical solution of Hamiltonian systems.

Published

September 16, 2022

Lissajous curves

usingPlotsa = [1,3,5,3]b = [1,5,7,4]p = []δ =π/2t =LinRange(-π, π, 300)for i ∈1:4 x =sin.(a[i] .* t .+ δ) y =sin.(b[i] .* t)push!(p, plot(x, y, legend =false, aspect_ratio=:equal))endplot(p..., layout=4)

┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
└ @ Base loading.jl:1664

┌ Info: Precompiling GR_jll [d2c73de3-f751-5644-a686-071e5b155ba9]
└ @ Base loading.jl:1664

Numerical Solution of ODEs

Julia has a great ODE package, providing a vast number of algorithms. The DiffEqPhys.jl package allows to use them for the numerical solution of Hamiltonian systems:

Provide:

a Hamiltonian function H(p,q) where p and q can be scalars, vectors or arrays of arbitrary (but equal) dimension,

initial values for p and q and

a timespan

and the function HamiltonianProblem(H, p₀, q₀, timespan,..) uses automatic differentiation to generate the system of first order Hamiltonian equations \[\begin{align}
\dot q &=\frac{\partial H}{\partial p} \\
\dot p &=-\frac{\partial H}{\partial q}.
\end{align}
\] in a form suitable for the solvers of the DifferentialEquations.jl package.

Simple Example: A charge moving in an homogeneous magnetic field

Let \(\mathbf{B} = (0,0,B)\) be a field in \(z\)-direction. A charge with some initial momentum is coerced by the Lorentz force into a circular motion in the perpendicular \(xy\) plane.

(In quantum mechanics the energy and therefore the radii of these orbits are quantized. This is the very first step to understand the quantum Hall effect.)

If the initial momentum has a component along the \(z\) axis, the charge will move on a helix.

The Hamiltonian for the problem is \[
H = \frac{(\mathbf{p}-e\mathbf{A})^2}{2m}\qquad\text{with the potential}\qquad \mathbf{A}=(0, Bx, 0)
\]

and that’s the Julia code:

usingDifferentialEquations, DiffEqPhysics, PlotsfunctionH(p, q, parms) (; m, B) = parms h = p[1]^2+ (p[2] - B*q[1])^2+ p[3]^2return h/2mendmyparms = (m=1.0, B=4.0)p₀ = [0, .1, .04]q₀ = [0.06, 0, 0 ]timespan = (0, 9)a =HamiltonianProblem(H, p₀, q₀, timespan, myparms)sol =solve(a, alg=OwrenZen4()) # Owren-Zennaro optimized 4/3 method plot(sol, legend =false, vars=(4,5,6)) # vars 1,2,3 are the p_i and 4,5,6 the q_i