Some Julia Code

julia
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

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