Classical Physics Models

If you're getting some cold feet to jump in to DiffEq land, here are some handcrafted differential equations mini problems to hold your hand along the beginning of your journey.

First order linear ODE

Radioactive Decay of Carbon-14

\[f(t,u) = \frac{du}{dt}\]

The Radioactive decay problem is the first order linear ODE problem of an exponential with a negative coefficient, which represents the half-life of the process in question. Should the coefficient be positive, this would represent a population growth equation.

using OrdinaryDiffEq, Plots
gr()

#Half-life of Carbon-14 is 5,730 years.
t½ = 5.730

#Setup
u₀ = 1.0
tspan = (0.0, 30.0)

#Define the problem
radioactivedecay(u, p, t) = -log(2) / t½ * u

#Pass to solver
prob = ODEProblem(radioactivedecay, u₀, tspan)
sol = solve(prob, Tsit5())

#Plot
plot(sol, linewidth = 2, title = "Carbon-14 half-life",
    xaxis = "Time in thousands of years", yaxis = "Ratio left",
    label = "Numerical Solution")
plot!(sol.t, t -> 2^(-t / t½), lw = 3, ls = :dash, label = "Analytical Solution")

Second Order Linear ODE

Simple Harmonic Oscillator

Another classical example is the harmonic oscillator, given by:

\[\ddot{x} + \omega^2 x = 0\]

with the known analytical solution

\[\begin{align*} x(t) &= A\cos(\omega t - \phi) \\ v(t) &= -A\omega\sin(\omega t - \phi), \end{align*}\]

where

\[A = \sqrt{c_1 + c_2} \qquad\text{and}\qquad \tan \phi = \frac{c_2}{c_1}\]

with $c_1, c_2$ constants determined by the initial conditions such that $c_1$ is the initial position and $\omega c_2$ is the initial velocity.

Instead of transforming this to a system of ODEs to solve with ODEProblem, we can use SecondOrderODEProblem as follows.

# Simple Harmonic Oscillator Problem
using OrdinaryDiffEq, Plots

#Parameters
ω = 1

#Initial Conditions
x₀ = [0.0]
dx₀ = [π / 2]
tspan = (0.0, 2π)

ϕ = atan((dx₀[1] / ω) / x₀[1])
A = √(x₀[1]^2 + dx₀[1]^2)

#Define the problem
function harmonicoscillator(ddu, du, u, ω, t)
    ddu .= -ω^2 * u
end

#Pass to solvers
prob = SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
sol = solve(prob, DPRKN6())

#Plot
plot(sol, idxs = [2, 1], linewidth = 2, title = "Simple Harmonic Oscillator",
    xaxis = "Time", yaxis = "Elongation", label = ["x" "dx"])
plot!(t -> A * cos(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution x")
plot!(t -> -A * ω * sin(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution dx")

Note that the order of the variables (and initial conditions) is dx, x. Thus, if we want the first series to be x, we have to flip the order with vars=[2,1].

Second Order Non-linear ODE

Simple Pendulum

We will start by solving the pendulum problem. In the physics class, we often solve this problem by small angle approximation, i.e. $ sin(\theta) \approx \theta$, because otherwise, we get an elliptic integral which doesn't have an analytic solution. The linearized form is

\[\ddot{\theta} + \frac{g}{L}{\theta} = 0\]

But we have numerical ODE solvers! Why not solve the real pendulum?

\[\ddot{\theta} + \frac{g}{L}{\sin(\theta)} = 0\]

Notice that now we have a second order ODE. In order to use the same method as above, we need to transform it into a system of first order ODEs by employing the notation $d\theta = \dot{\theta}$.

\[\begin{align*} &\dot{\theta} = d{\theta} \\ &\dot{d\theta} = - \frac{g}{L}{\sin(\theta)} \end{align*}\]

# Simple Pendulum Problem
using OrdinaryDiffEq, Plots

#Constants
const g = 9.81
L = 1.0

#Initial Conditions
u₀ = [0, π / 2]
tspan = (0.0, 6.3)

#Define the problem
function simplependulum(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g / L) * sin(θ)
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5())

#Plot
plot(sol, linewidth = 2, title = "Simple Pendulum Problem", xaxis = "Time",
    yaxis = "Height", label = ["\\theta" "d\\theta"])

So now we know that behaviour of the position versus time. However, it will be useful to us to look at the phase space of the pendulum, i.e., and representation of all possible states of the system in question (the pendulum) by looking at its velocity and position. Phase space analysis is ubiquitous in the analysis of dynamical systems, and thus we will provide a few facilities for it.

p = plot(sol, vars = (1, 2), xlims = (-9, 9), title = "Phase Space Plot",
    xaxis = "Angular position", yaxis = "Angular velocity", leg = false)
function phase_plot(prob, u0, p, tspan = 2pi)
    _prob = ODEProblem(prob.f, u0, (0.0, tspan))
    sol = solve(_prob, Vern9()) # Use Vern9 solver for higher accuracy
    plot!(p, sol, idxs = (1, 2), xlims = nothing, ylims = nothing)
end
for i in (-4pi):(pi / 2):(4π)
    for j in (-4pi):(pi / 2):(4π)
        phase_plot(prob, [j, i], p)
    end
end
plot(p, xlims = (-9, 9))

Double Pendulum

A more complicated example is given by the double pendulum. The equations governing its motion are given by the following (taken from this Stack Overflow question)

\[\frac{d}{dt} \begin{pmatrix} \alpha \\ l_\alpha \\ \beta \\ l_\beta \end{pmatrix}= \begin{pmatrix} 2\frac{l_\alpha - (1+\cos\beta)l_\beta}{3-\cos 2\beta} \\ -2\sin\alpha - \sin(\alpha + \beta) \\ 2\frac{-(1+\cos\beta)l_\alpha + (3+2\cos\beta)l_\beta}{3-\cos2\beta}\\ -\sin(\alpha+\beta) - 2\sin(\beta)\frac{(l_\alpha-l_\beta)l_\beta}{3-\cos2\beta} + 2\sin(2\beta)\frac{l_\alpha^2-2(1+\cos\beta)l_\alpha l_\beta + (3+2\cos\beta)l_\beta^2}{(3-\cos2\beta)^2} \end{pmatrix}\]

#Double Pendulum Problem
using OrdinaryDiffEq, Plots

#Constants and setup
const m₁, m₂, L₁, L₂ = 1, 2, 1, 2
initial = [0, π / 3, 0, 3pi / 5]
tspan = (0.0, 50.0)

#Convenience function for transforming from polar to Cartesian coordinates
function polar2cart(sol; dt = 0.02, l1 = L₁, l2 = L₂, vars = (2, 4))
    u = sol.t[1]:dt:sol.t[end]

    p1 = l1 * map(x -> x[vars[1]], sol.(u))
    p2 = l2 * map(y -> y[vars[2]], sol.(u))

    x1 = l1 * sin.(p1)
    y1 = l1 * -cos.(p1)
    (u, (x1 + l2 * sin.(p2),
        y1 - l2 * cos.(p2)))
end

#Define the Problem
function double_pendulum(xdot, x, p, t)
    xdot[1] = x[2]
    xdot[2] = -((g * (2 * m₁ + m₂) * sin(x[1]) +
                 m₂ * (g * sin(x[1] - 2 * x[3]) +
                  2 * (L₂ * x[4]^2 + L₁ * x[2]^2 * cos(x[1] - x[3])) * sin(x[1] - x[3]))) /
                (2 * L₁ * (m₁ + m₂ - m₂ * cos(x[1] - x[3])^2)))
    xdot[3] = x[4]
    xdot[4] = (((m₁ + m₂) * (L₁ * x[2]^2 + g * cos(x[1])) +
                L₂ * m₂ * x[4]^2 * cos(x[1] - x[3])) * sin(x[1] - x[3])) /
              (L₂ * (m₁ + m₂ - m₂ * cos(x[1] - x[3])^2))
end

#Pass to Solvers
double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan)
sol = solve(double_pendulum_problem, Vern7(), abstol = 1e-10, dt = 0.05);

retcode: Success
Interpolation: specialized 7th order lazy interpolation
t: 300-element Vector{Float64}:
  0.0
  0.05
  0.12217751622058524
  0.21737046905775836
  0.32557870861756055
  0.45380019052324166
  0.6093195606053129
  0.7728295835307077
  0.9531773015934449
  1.1789364520938335
  ⋮
 48.569380798168076
 48.756377534074545
 48.927478887477086
 49.12151506174158
 49.30166659605193
 49.46644518867691
 49.644275424816705
 49.78679810755379
 50.0
u: 300-element Vector{Vector{Float64}}:
 [0.0, 1.0471975511965976, 0.0, 1.8849555921538759]
 [0.05276815671595484, 1.0714389570723506, 0.09384176137401032, 1.8607344940613868]
 [0.13346405022464972, 1.1746261372772198, 0.22461544271215322, 1.7498682922370123]
 [0.2534001539025762, 1.33970622012998, 0.38068162944805334, 1.5194590529384873]
 [0.40362066930016344, 1.402594495269176, 0.5297354554896815, 1.240443705598302]
 [0.5722968387882038, 1.171450176324233, 0.6713718220157865, 0.9861901962522713]
 [0.7085423339309261, 0.539903531551119, 0.8079790486139692, 0.7794821222233647]
 [0.7355031474073365, -0.1908205610455411, 0.9166518850457744, 0.5304216702747947]
 [0.6478012453021919, -0.713845531236494, 0.9733991099427944, 0.057112403958845384]
 [0.4700665697655869, -0.734304363798229, 0.888698434825054, -0.8582765324454241]
 ⋮
 [-0.6555475350480671, -0.43413585939726496, -0.8218749873609228, 1.3383995922630012]
 [-0.6545393233636478, 0.6393519978284158, -0.555165884738273, 1.3845271567723367]
 [-0.4224791976270007, 2.018551723548183, -0.3471917114189364, 1.0423940791239785]
 [0.01664985021925263, 2.10806421640763, -0.14266655415801852, 1.2529804077715108]
 [0.2827442220914206, 0.7861283125754889, 0.14782014738717986, 1.9484415252454967]
 [0.34628967140864897, 0.23457916223370065, 0.48488462162712476, 1.97438594775981]
 [0.4277722810740152, 0.7106352082687376, 0.7731951730946816, 1.2198285330824574]
 [0.5465172554749359, 0.8994327527768179, 0.9011967077171154, 0.5915426120472353]
 [0.7136904261911289, 0.5470830780822246, 0.9471656849456983, -0.10172229763994178]

#Obtain coordinates in Cartesian Geometry
ts, ps = polar2cart(sol, l1 = L₁, l2 = L₂, dt = 0.01)
plot(ps...)

Poincaré section

In this case, the phase space is 4 dimensional, and it cannot be easily visualized. Instead of looking at the full phase space, we can look at Poincaré sections, which are sections through a higher-dimensional phase space diagram. This helps to understand the dynamics of interactions and is wonderfully pretty.

The Poincaré section in this is given by the collection of $(β,l_β)$ when $α=0$ and $\frac{dα}{dt}>0$.

#Constants and setup
using OrdinaryDiffEq
initial2 = [0.01, 0.005, 0.01, 0.01]
tspan2 = (0.0, 500.0)

#Define the problem
function double_pendulum_hamiltonian(udot, u, p, t)
    α = u[1]
    lα = u[2]
    β = u[3]
    lβ = u[4]
    udot .= [2(lα - (1 + cos(β))lβ) / (3 - cos(2β)),
        -2sin(α) - sin(α + β),
        2(-(1 + cos(β))lα + (3 + 2cos(β))lβ) / (3 - cos(2β)),
        -sin(α + β) - 2sin(β) * (((lα - lβ)lβ) / (3 - cos(2β))) +
        2sin(2β) * ((lα^2 - 2(1 + cos(β))lα * lβ + (3 + 2cos(β))lβ^2) / (3 - cos(2β))^2)]
end

# Construct a ContiunousCallback
condition(u, t, integrator) = u[1]
affect!(integrator) = nothing
cb = ContinuousCallback(condition, affect!, nothing,
    save_positions = (true, false))

# Construct Problem
poincare = ODEProblem(double_pendulum_hamiltonian, initial2, tspan2)
sol2 = solve(poincare, Vern9(), save_everystep = false, save_start = false,
    save_end = false, callback = cb, abstol = 1e-16, reltol = 1e-16)

function poincare_map(prob, u₀, p; callback = cb)
    _prob = ODEProblem(prob.f, u₀, prob.tspan)
    sol = solve(_prob, Vern9(), save_everystep = false, save_start = false,
        save_end = false, callback = cb, abstol = 1e-16, reltol = 1e-16)
    scatter!(p, sol, idxs = (3, 4), markersize = 3, msw = 0)
end

poincare_map (generic function with 1 method)

lβrange = -0.02:0.0025:0.02
p = scatter(sol2, idxs = (3, 4), leg = false, markersize = 3, msw = 0)
for lβ in lβrange
    poincare_map(poincare, [0.01, 0.01, 0.01, lβ], p)
end
plot(p, xlabel = "\\beta", ylabel = "l_\\beta", ylims = (0, 0.03))

Hénon-Heiles System

The Hénon-Heiles potential occurs when non-linear motion of a star around a galactic center, with the motion restricted to a plane.

\[\begin{align} \frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x}\\ \frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y} \end{align}\]

where

\[V(x,y)={\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).\]

We pick $\lambda=1$ in this case, so

\[V(x,y) = \frac{1}{2}(x^2+y^2+2x^2y-\frac{2}{3}y^3).\]

Then the total energy of the system can be expressed by

\[E = T+V = V(x,y)+\frac{1}{2}(\dot{x}^2+\dot{y}^2).\]

The total energy should conserve as this system evolves.

using OrdinaryDiffEq, Plots

#Setup
initial = [0.0, 0.1, 0.5, 0]
tspan = (0, 100.0)

#Remember, V is the potential of the system and T is the Total Kinetic Energy, thus E will
#the total energy of the system.
V(x, y) = 1 // 2 * (x^2 + y^2 + 2x^2 * y - 2 // 3 * y^3)
E(x, y, dx, dy) = V(x, y) + 1 // 2 * (dx^2 + dy^2);

#Define the function
function Hénon_Heiles(du, u, p, t)
    x = u[1]
    y = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x * y
    du[4] = y^2 - y - x^2
end

#Pass to solvers
prob = ODEProblem(Hénon_Heiles, initial, tspan)
sol = solve(prob, Vern9(), abstol = 1e-16, reltol = 1e-16);

retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 1845-element Vector{Float64}:
   0.0
   0.011644644151844885
   0.02264023714416834
   0.03374388443401237
   0.04788954072316613
   0.061013377791614205
   0.07550280255229364
   0.09098576673379825
   0.10692621797677565
   0.12281185948039652
   ⋮
  99.57698660332744
  99.63288223314693
  99.69002695380586
  99.74295860063948
  99.79568139744299
  99.85043316365166
  99.90440973346377
  99.95982644217831
 100.0
u: 1845-element Vector{Vector{Float64}}:
 [0.0, 0.1, 0.5, 0.0]
 [0.005822164178948815, 0.09999389777392866, 0.49995932143722593, -0.0010481306031208865]
 [0.011318958127118273, 0.09997692919980947, 0.49984623673656015, -0.0020384490237147167]
 [0.016868100331427598, 0.09994873764623995, 0.4996584585603103, -0.00303968938416814]
 [0.023933789413513467, 0.09989670288341591, 0.4993121945281112, -0.004317887366900626]
 [0.03048398283438544, 0.09983223461360576, 0.49888378195193545, -0.0055073863109945005]
 [0.03770837983827229, 0.09974289088219766, 0.4982911385372982, -0.006825898266222505]
 [0.04541761306212422, 0.09962625160368373, 0.4975193039185924, -0.008242303637579974]
 [0.05334097324069581, 0.09948317991794156, 0.49657542702808843, -0.009710240734282243]
 [0.061220926128609424, 0.09931723251011267, 0.4954845574373833, -0.01118455239229768]
 ⋮
 [0.11422633985261582, 0.42680565269768317, 0.32275163730669065, -0.025465965445907804]
 [0.13191894824167866, 0.42497757894196253, 0.3100060322909237, -0.039982355130911675]
 [0.14921844338165957, 0.42226299237891685, 0.295159422633226, -0.055067468788503944]
 [0.16444331096031972, 0.4189732564497299, 0.27986756414426317, -0.0692713568722036]
 [0.17876687626283724, 0.4149431836821891, 0.26326751565019746, -0.08364378323702079]
 [0.1926788582043158, 0.409949595072511, 0.24470157212591676, -0.09880372526618875]
 [0.20536582405914175, 0.40420808465311453, 0.22520350419186239, -0.11397296670416476]
 [0.2172657238232399, 0.3974555823335333, 0.2040945571150857, -0.12976028393356268]
 [0.22514698918543646, 0.3920106535852793, 0.18818797399244314, -0.1413256709617241]

# Plot the orbit
plot(sol, idxs = (1, 2), title = "The orbit of the Hénon-Heiles system", xaxis = "x",
    yaxis = "y", leg = false)

#Optional Sanity check - what do you think this returns and why?
@show sol.retcode

#Plot -
plot(sol, idxs = (1, 3), title = "Phase space for the Hénon-Heiles system",
    xaxis = "Position", yaxis = "Velocity")
plot!(sol, idxs = (2, 4), leg = false)

#We map the Total energies during the time intervals of the solution (sol.u here) to a new vector
#pass it to the plotter a bit more conveniently
energy = map(x -> E(x...), sol.u)

#We use @show here to easily spot erratic behavior in our system by seeing if the loss in energy was too great.
@show ΔE = energy[1] - energy[end]

#Plot
plot(sol.t, energy .- energy[1], title = "Change in Energy over Time",
    xaxis = "Time in iterations", yaxis = "Change in Energy")

Symplectic Integration

To prevent energy drift, we can instead use a symplectic integrator. We can directly define and solve the SecondOrderODEProblem:

function HH_acceleration!(dv, v, u, p, t)
    x, y = u
    dx, dy = dv
    dv[1] = -x - 2x * y
    dv[2] = y^2 - y - x^2
end
initial_positions = [0.0, 0.1]
initial_velocities = [0.5, 0.0]
prob = SecondOrderODEProblem(HH_acceleration!, initial_velocities, initial_positions, tspan)
sol2 = solve(prob, KahanLi8(), dt = 1 / 10);

retcode: Success
Interpolation: 3rd order Hermite
t: 1001-element Vector{Float64}:
   0.0
   0.1
   0.2
   0.30000000000000004
   0.4
   0.5
   0.6
   0.7
   0.7999999999999999
   0.8999999999999999
   ⋮
  99.19999999999864
  99.29999999999863
  99.39999999999863
  99.49999999999862
  99.59999999999862
  99.69999999999861
  99.7999999999986
  99.8999999999986
 100.0
u: 1001-element Vector{RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
 ([0.5, 0.0], [0.0, 0.1])
 ([0.497004124813899, -0.009071101031878595], [0.049900082497367014, 0.09954822053953173])
 ([0.488065986503409, -0.01856325999777532], [0.09920263962777168, 0.0981717138550656])
 ([0.4733339415930383, -0.028870094978230895], [0.1473200401467506, 0.09580835287880228])
 ([0.45305474121099776, -0.040331734000937175], [0.1936844146808317, 0.09235911550101675])
 ([0.4275722362076315, -0.0532114683862374], [0.2377574398051348, 0.08769462857458447])
 ([0.39732459499168415, -0.06767627591286327], [0.279039924450448, 0.08166386202131776])
 ([0.36283926518715554, -0.08378216859669584], [0.3170810117622551, 0.07410453631898412])
 ([0.3247249463611115, -0.10146505675930295], [0.35148673325356056, 0.06485472395637552])
 ([0.2836600192614762, -0.12053752272622162], [0.3819275865822665, 0.05376507097805092])
 ⋮
 ([0.35754893663581877, 0.06817150704505637], [-0.01689940464435974, 0.41858730953689555])
 ([0.3573608346071073, 0.043775799566172786], [0.018901094264629065, 0.42418546702794857])
 ([0.35056546788922516, 0.01917894635581023], [0.054352325047751296, 0.4273357603171658])
 ([0.3372619546371046, -0.00582469608876845], [0.08879705585509527, 0.4280076678341491])
 ([0.317723101892544, -0.031415069954177], [0.12159668542307678, 0.42615121224761354])
 ([0.2923883191107272, -0.057726485993246215], [0.15214830714764738, 0.42170054930538214])
 ([0.2618505981512214, -0.0848309469334548], [0.17990076750513168, 0.4145793965132324])
 ([0.22683770400643935, -0.11272570532234003], [0.2043691308538718, 0.40470792450133325])
 ([0.1881879739856503, -0.1413256709667938], [0.22514698919068493, 0.39201065357989734])

Notice that we get the same results:

# Plot the orbit
plot(sol2, idxs = (3, 4), title = "The orbit of the Hénon-Heiles system", xaxis = "x",
    yaxis = "y", leg = false)

plot(sol2, idxs = (3, 1), title = "Phase space for the Hénon-Heiles system",
    xaxis = "Position", yaxis = "Velocity")
plot!(sol2, idxs = (4, 2), leg = false)

but now the energy change is essentially zero:

energy = map(x -> E(x[3], x[4], x[1], x[2]), sol2.u)
#We use @show here to easily spot erratic behaviour in our system by seeing if the loss in energy was too great.
@show ΔE = energy[1] - energy[end]

#Plot
plot(sol2.t, energy .- energy[1], title = "Change in Energy over Time",
    xaxis = "Time in iterations", yaxis = "Change in Energy")

And let's try to use a Runge-Kutta-Nyström solver to solve this. Note that Runge-Kutta-Nyström isn't symplectic.

sol3 = solve(prob, DPRKN6());
energy = map(x -> E(x[3], x[4], x[1], x[2]), sol3.u)
@show ΔE = energy[1] - energy[end]
gr()
plot(sol3.t, energy .- energy[1], title = "Change in Energy over Time",
    xaxis = "Time in iterations", yaxis = "Change in Energy")

Note that we are using the DPRKN6 solver at reltol=1e-3 (the default), yet it has a smaller energy variation than Vern9 at abstol=1e-16, reltol=1e-16. Therefore, using specialized solvers to solve its particular problem is very efficient.