Ordinary Differential Equations (ODEs)

• DifferentialEquations.jl docs

General steps to solve ODEs (in the old way)

• Define a model function representing the right-hand-side (RHS) of the sysstem.

  • Out-of-place form: f(u, p, t) where u is the state variable(s), p is the parameter(s), and t is the independent variable (usually time). The output is the right hand side (RHS) of the differential equation system.

  • In-place form: f!(du, u, p, t), where the output is saved to du. The rest is the same as the out-of-place form. The in-place function may be faster since it does not allocate a new array in each call.

• Initial conditions (u0) for the state variable(s).

• (Optional) parameter(s) p.

• Define a problem (e.g. ODEProblem) using the model function (f), initial condition(s) (u0), simulation time span (tspan == (tstart, tend)), and parameter(s) p.

• Solve the problem by calling solve(prob).

Radioactive decay

The decaying rate of a nuclear isotope is proprotional to its concentration (number):

\[ \frac{d}{dt}C(t) = - \lambda C(t) \]

State variable(s)

• \(C(t)\): The concentration (number) of a decaying nuclear isotope.

Parameter(s)

• \(\lambda\): The rate constant of nuclear decay. The half-life: \(t_{\frac{1}{2}} = \frac{ln2}{\lambda}\).

using DifferentialEquations
using Plots

[ Info: Precompiling IJuliaExt [2f4121a4-3b3a-5ce6-9c5e-1f2673ce168a]

The exponential decay ODE model, out-of-place (3-parameter) form

expdecay(u, p, t) = p * u

expdecay (generic function with 1 method)

p = -1.0 ## Parameter
u0 = 1.0 ## Initial condition
tspan = (0.0, 2.0) ## Simulation start and end time points
prob = ODEProblem(expdecay, u0, tspan, p) ## Define the problem
sol = solve(prob) ## Solve the problem

retcode: Success
Interpolation: 3rd order Hermite
t: 8-element Vector{Float64}:
 0.0
 0.10001999200479662
 0.34208427873632274
 0.6553980290285384
 1.0312652902321524
 1.4709406498424789
 1.9659577002710475
 2.0
u: 8-element Vector{Float64}:
 1.0
 0.9048193287657775
 0.7102883564034526
 0.5192354320104405
 0.35655575232768816
 0.2297097760377979
 0.14002246806154703
 0.13533600284000216

Visualize the solution with Plots.jl.

plot(sol, label="Exp decay")

Solution handling: https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/ The mostly used feature is sol(t), the state variables at time t. t could be a scalar or a vector-like sequence. The result state variables are calculated with interpolation.

sol(1.0)

0.3713271431453589

sol(0.0:0.1:2.0)

t: 0.0:0.1:2.0
u: 21-element Vector{Float64}:
 1.0
 0.9048389983836875
 0.8225844481735297
 0.7453905798423603
 0.6716690052872237
 0.6126431137597175
 0.554633689939333
 0.49719361928210143
 0.45380162043196953
 0.41384825601313535
 0.3713271431453589
 0.3337844113123039
 0.30522377530459677
 0.2789957673333738
 0.2516591524294101
 0.22324380097718674
 0.20365547977113405
 0.18672102774211405
 0.17051216093498966
 0.15310059539462373
 0.13533600284000216

sol.t: time points by the solver.

sol.t

8-element Vector{Float64}:
 0.0
 0.10001999200479662
 0.34208427873632274
 0.6553980290285384
 1.0312652902321524
 1.4709406498424789
 1.9659577002710475
 2.0

sol.u: state variables at sol.t.

sol.u

8-element Vector{Float64}:
 1.0
 0.9048193287657775
 0.7102883564034526
 0.5192354320104405
 0.35655575232768816
 0.2297097760377979
 0.14002246806154703
 0.13533600284000216

The SIR model

A more complicated example is the SIR model describing infectious disease spreading. There are 3 state variables and 2 parameters.

\[\begin{split} \begin{align} \frac{d}{dt}S(t) &= - \beta S(t)I(t) \\ \frac{d}{dt}I(t) &= \beta S(t)I(t) - \gamma I(t) \\ \frac{d}{dt}R(t) &= \gamma I(t) \end{align} \end{split}\]

State variable(s)

• \(S(t)\) : the fraction of susceptible people

• \(I(t)\) : the fraction of infectious people

• \(R(t)\) : the fraction of recovered (or removed) people

Parameter(s)

• \(\beta\) : the rate of infection when susceptible and infectious people meet

• \(\gamma\) : the rate of recovery of infectious people

using DifferentialEquations
using Plots

Here we use the in-place form: f!(du, u, p ,t) for the SIR model. The output is written to the first argument du, without allocating a new array in each function call.

function sir!(du, u, p, t)
    s, i, r = u
    β, γ = p
    v1 = β * s * i
    v2 = γ * i
    du[1] = -v1
    du[2] = v1 - v2
    du[3] = v2
    return nothing
end

sir! (generic function with 1 method)

Setup parameters, initial conditions, time span, and the ODE problem.

p = (β=1.0, γ=0.3)
u0 = [0.99, 0.01, 0.00]
tspan = (0.0, 20.0)
prob = ODEProblem(sir!, u0, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 20.0)
u0: 3-element Vector{Float64}:
 0.99
 0.01
 0.0

Solve the problem

sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 17-element Vector{Float64}:
  0.0
  0.08921318693905476
  0.3702862839612021
  0.798425739836198
  1.323727172235353
  1.9918418559528543
  2.7923706904146552
  3.7547815921022205
  4.90190421552381
  6.260476579481841
  7.764891224501251
  9.390409694084836
 11.483860862519641
 13.372369625696377
 15.961356868207503
 18.681426190579412
 20.0
u: 17-element Vector{Vector{Float64}}:
 [0.99, 0.01, 0.0]
 [0.9890894703413342, 0.010634484617786016, 0.00027604504087978485]
 [0.985833159331863, 0.012901496935960089, 0.0012653437321769186]
 [0.9795270525087807, 0.01728242130889301, 0.0031905261823262586]
 [0.9689082161766913, 0.024631267424412892, 0.006460516398895782]
 [0.9490552303913907, 0.038273388553893645, 0.012671381054715575]
 [0.9118629477834835, 0.06347250081443477, 0.024664551402081673]
 [0.8398871109908902, 0.11078175898927324, 0.04933113001983658]
 [0.7075842208271164, 0.19166147074418158, 0.10075430842870206]
 [0.5081460371757582, 0.29177418995172044, 0.2000797728725213]
 [0.3121322220064624, 0.34158791193414223, 0.34627986605939537]
 [0.18215683750565154, 0.3099983176471322, 0.5078448448472163]
 [0.104272058384535, 0.2206111470381464, 0.6751167945773185]
 [0.07386737657579129, 0.14760143815780902, 0.7785311852663996]
 [0.055450290457826595, 0.07997077517076541, 0.864578934371408]
 [0.047334991613510013, 0.040605658108532304, 0.9120593502779577]
 [0.04522885458990312, 0.029057416110837522, 0.9257137292992594]

Visualize the solution

plot(sol, labels=["S" "I" "R"], legend=:right)

sol[i]: all components at timestep i

sol[2]

3-element Vector{Float64}:
 0.9890894703413342
 0.010634484617786016
 0.00027604504087978485

sol[i, j]: ith component at timestep j

sol[1, 2]

0.9890894703413342

sol[i, :]: the timeseries for the ith component.

sol[1, :]

17-element Vector{Float64}:
 0.99
 0.9890894703413342
 0.985833159331863
 0.9795270525087807
 0.9689082161766913
 0.9490552303913907
 0.9118629477834835
 0.8398871109908902
 0.7075842208271164
 0.5081460371757582
 0.3121322220064624
 0.18215683750565154
 0.104272058384535
 0.07386737657579129
 0.055450290457826595
 0.047334991613510013
 0.04522885458990312

sol(t,idxs=1): the 1st element in time point(s) t with interpolation. t can be a scalar or a vector.

sol(10, idxs=2)

0.28561290260161176

sol(0.0:0.1:20.0, idxs=2)

t: 0.0:0.1:20.0
u: 201-element Vector{Float64}:
 0.01
 0.010714335106962786
 0.011511657165648073
 0.012345202527705887
 0.013171799248016437
 0.01416833789987141
 0.015239587721803413
 0.01630823631440291
 0.017300937892400663
 0.018580292807160968
 0.02000659659162972
 0.021496287236294518
 0.022965802731642962
 ⋮
 0.038492182874362846
 0.03759686349556628
 0.03673462162339544
 0.035896938976267445
 0.035075297272599344
 0.034261178230808306
 0.03344606356931138
 0.03262143500652569
 0.03177877426086837
 0.030909563050756483
 0.030005283094607187
 0.029057416110837522

Lorenz system

The Lorenz system is a system of ordinary differential equations having chaotic solutions for certain parameter values and initial conditions. (Wikipedia)

\[\begin{split} \begin{align} \frac{dx}{dt} &=& \sigma(y-x) \\ \frac{dy}{dt} &=& x(\rho - z) -y \\ \frac{dz}{dt} &=& xy - \beta z \end{align} \end{split}\]

In this example, we will use NamedTuple for the parameters and ComponentArrays.jl for the state variaBLE to access elements by their names.

using ComponentArrays
using DifferentialEquations
using Plots

function lorenz!(du, u, p, t)
    du.x = p.σ * (u.y - u.x)
    du.y = u.x * (p.ρ - u.z) - u.y
    du.z = u.x * u.y - p.β * u.z
    return nothing
end

lorenz! (generic function with 1 method)

u0 = ComponentArray(x=1.0, y=0.0, z=0.0)
p = (σ=10.0, ρ=28.0, β=8 / 3)
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan, p)
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 1281-element Vector{Float64}:
   0.0
   3.5678604836301404e-5
   0.0003924646531993154
   0.0032624087100077666
   0.009058076582749423
   0.016956470605311864
   0.027689959227781235
   0.04185635103821218
   0.060240410627700816
   0.0836854113984534
   0.11336499269451543
   0.14862181409827
   0.18703978025370946
   ⋮
  99.19664352825902
  99.28542618302045
  99.3530549740687
  99.4386918268877
  99.53430077629385
  99.62315909447312
  99.70299380266312
  99.7727905290976
  99.8468719463341
  99.92909223837675
  99.99676433660723
 100.0
u: 1281-element Vector{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x = 1, y = 2, z = 3)}}}}:
 ComponentVector{Float64}(x = 1.0, y = 0.0, z = 0.0)
 ComponentVector{Float64}(x = 0.9996434557625105, y = 0.0009988049817849058, z = 1.781434788799189e-8)
 ComponentVector{Float64}(x = 0.9961045497425811, y = 0.010965399721242457, z = 2.1469553658389193e-6)
 ComponentVector{Float64}(x = 0.9693591550149778, y = 0.08977063252764937, z = 0.0001438019170127846)
 ComponentVector{Float64}(x = 0.924204355043198, y = 0.242289149116772, z = 0.0010461625397616113)
 ComponentVector{Float64}(x = 0.8800455796215916, y = 0.4387364900041282, z = 0.0034242599668253874)
 ComponentVector{Float64}(x = 0.8483309836977301, y = 0.6915629475762161, z = 0.008487624968655677)
 ComponentVector{Float64}(x = 0.8495036679850485, y = 1.0145426495980538, z = 0.018212090123888365)
 ComponentVector{Float64}(x = 0.9139069520070257, y = 1.4425599557988966, z = 0.03669382071284047)
 ComponentVector{Float64}(x = 1.0888638157267199, y = 2.052326562845511, z = 0.07402573450924263)
 ComponentVector{Float64}(x = 1.4608626762755585, y = 3.0206719763684764, z = 0.160039355072216)
 ComponentVector{Float64}(x = 2.162723385463873, y = 4.6333636144497445, z = 0.3771173678677311)
 ComponentVector{Float64}(x = 3.3684642325435363, y = 7.267693727122867, z = 0.9363555414291395)
 ⋮
 ComponentVector{Float64}(x = 2.4154751319743304, y = 0.9019051839955854, z = 23.10904521543824)
 ComponentVector{Float64}(x = 1.9010431518181314, y = 2.075631191607503, z = 18.469762607052317)
 ComponentVector{Float64}(x = 2.3313747500610864, y = 3.427114121484939, z = 15.773720459660732)
 ComponentVector{Float64}(x = 3.8662182925863604, y = 6.477675362450072, z = 13.7026651297458)
 ComponentVector{Float64}(x = 7.596364784839832, y = 12.931698345552782, z = 15.295013031866944)
 ComponentVector{Float64}(x = 12.999465900688321, y = 18.69084638217891, z = 25.77003643982645)
 ComponentVector{Float64}(x = 14.725104821648667, y = 12.127682492556511, z = 37.85997432096073)
 ComponentVector{Float64}(x = 10.45813774612691, y = 2.2355517835584817, z = 37.102951902956555)
 ComponentVector{Float64}(x = 4.858569386282909, y = -1.238820833292597, z = 30.50492838887987)
 ComponentVector{Float64}(x = 1.4267543691218378, y = -1.0945246877110832, z = 24.22832564643038)
 ComponentVector{Float64}(x = 0.2911991490440339, y = -0.7513530324542533, z = 20.18373424388144)
 ComponentVector{Float64}(x = 0.25815851151351493, y = -0.7419175275518352, z = 20.00966894969765)

idxs=(1, 2, 3) makes a phase plot with 1st, 2nd, and the 3rd state variable.

plot(sol, idxs=(1, 2, 3), size=(400, 400))

The plot recipe is using interpolation for smoothing. You can turn off denseplot to see the difference.

plot(sol, idxs=(1, 2, 3), denseplot=false, size=(400, 400))

The zeroth variable in idxs is the independent variable (usually time). The FOLLOWING command plots the time series of the second state variable (y).

plot(sol, idxs=(0, 2))

Non-autonomous ODEs

In non-autonomous ODEs, term(s) in the right-hadn-side (RHS) maybe time-dependent. For example, in this pendulum model has an external, time-dependent, force.

\[\begin{split} \begin{aligned} \dot{\theta} &= \omega(t) \\ \dot{\omega} &= -1.5\frac{g}{l}sin(\theta(t)) + \frac{3}{ml^2}M(t) \end{aligned} \end{split}\]

• \(\theta\): pendulum angle

• \(\omega\): angular rate

• M: time-dependent external torgue

• \(l\): pendulum length

• \(g\): gravitional acceleration

using DifferentialEquations
using Plots

function pendulum!(du, u, p, t)
    l = 1.0                             ## length [m]
    m = 1.0                             ## mass [kg]
    g = 9.81                            ## gravitational acceleration [m/s²]
    du[1] = u[2]                        ## θ'(t) = ω(t)
    du[2] = -3g / (2l) * sin(u[1]) + 3 / (m * l^2) * p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end

u0 = [0.01, 0.0]                     ## initial angular deflection [rad] and angular velocity [rad/s]
tspan = (0.0, 10.0)                  ## time interval

M = t -> 0.1 * sin(t)                   ## external torque [Nm] as the parameter for the pendulum model

prob = ODEProblem(pendulum!, u0, tspan, M)
sol = solve(prob)

plot(sol, linewidth=2, xaxis="t", label=["θ [rad]" "ω [rad/s]"])
plot!(M, tspan..., label="Ext. force")

Linear ODE system

The ODE system could be anything as long as it returns the derivatives of state variables. In this example, the ODE system is described by a matrix differential operator.

\(\dot{u} = Au\)

using DifferentialEquations
using Plots

A = [
    1.0 0 0 -5
    4 -2 4 -3
    -4 0 0 1
    5 -2 2 3
]

u0 = rand(4, 2)

tspan = (0.0, 1.0)
f = (u, p, t) -> A * u
prob = ODEProblem(f, u0, tspan)
sol = solve(prob)
plot(sol)

---

This notebook was generated using Literate.jl.