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

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.
C₁ = 5.730

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

#Define the problem
radioactivedecay(u, p, t) = -C₁ * 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 = "Percentage left",
     label = "Numerical Solution")
plot!(sol.t, t -> exp(-C₁ * t), lw = 3, ls = :dash, label = "Analytical Solution")

Second Order Linear ODE

Simple Harmonic Oscillator

Another classical example is the harmonic oscillator, given by:

with the known analytical solution

where

with constants determined by the initial conditions such that is the initial position and 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, vars = [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. , because otherwise, we get an elliptic integral which doesn’t have an analytic solution. The linearized form is

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

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 .

# 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 = "Velocity", yaxis = "Position", 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, vars = (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)

#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: 303-element Vector{Float64}:
  0.0
  0.05
  0.12217751559731858
  0.21737046504966892
  0.32557871102475433
  0.4538001831296513
  0.6093195612067883
  0.7728295888824992
  0.9531773177034342
  1.1789364907473971
  ⋮
 48.69306260963273
 48.91021204184887
 49.07673144456275
 49.26950145313383
 49.413714336117856
 49.58639229759933
 49.735034200347606
 49.92709346084078
 50.0
u: 303-element Vector{Vector{Float64}}:
 [0.0, 1.0471975511965976, 0.0, 1.8849555921538759]
 [0.05276815671595484, 1.071438957072351, 0.09384176137401032, 1.8607344940613866]
 [0.13346404949254442, 1.1746261361732084, 0.22461544162151867, 1.7498682935067824]
 [0.25340014853291515, 1.3397062143942808, 0.38068162335792477, 1.5194590635443397]
 [0.40362067267647705, 1.4025944942000763, 0.5297354584756716, 1.240443699864503]
 [0.5722968301269858, 1.171450199082293, 0.6713718147242973, 0.9861902079902164]
 [0.7085423342556784, 0.5399035287557558, 0.8079790490828, 0.7794821214727354]
 [0.7355031463861064, -0.1908205822316285, 0.9166518878844778, 0.5304216598453477]
 [0.6478012338021546, -0.7138455590058932, 0.9733991108628701, 0.057112350467565534]
 [0.4700665413819645, -0.734304300371504, 0.888698401649624, -0.8582767128313853]
 ⋮
 [-0.6790142112405692, 0.1666713249565996, -0.6449607059987224, 1.445707398332814]
 [-0.453711223749631, 1.907772695621923, -0.3662711842230573, 1.0713657099932998]
 [-0.08088755404381312, 2.283085230936697, -0.19537310870031996, 1.1164019590612861]
 [0.25211550091772744, 1.0205834408401855, 0.08848962912628953, 1.8493662815063088]
 [0.3327108424404666, 0.2519694839950239, 0.37935712531815796, 2.073023391326819]
 [0.39200740461082473, 0.5635940086140278, 0.6950083099731763, 1.4813704445243048]
 [0.5030292288143082, 0.8775432091086321, 0.8641350711861532, 0.8027494923114267]
 [0.6687647784296428, 0.7437942081020751, 0.9462230447555693, 0.0961700207091559]
 [0.7159266100692314, 0.5372537705329169, 0.9459035223663264, -0.09804066261876586]

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