Some Julia Code

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.

Lissajous curves

using Plots

a = [1,3,5,3]
b = [1,5,7,4]
p = []
δ = π/2
t = 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))
end

plot(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:

using DifferentialEquations, DiffEqPhysics, Plots

function H(p, q, parms)
    (; m, B) = parms
    h = p[1]^2 + (p[2] - B*q[1])^2 + p[3]^2
    return h/2m
end

myparms = (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