Solving differential equations in Julia

Solving differential equations in Julia#

Standard procedures

  • 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 form has potential performance benefits since it allocates less than the out-of-place (f(u, p, t)) counterpart.

    • Using ModelingToolkit.jl : define equations and build an ODE system.

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

  • (Optional) define parameter(s) p.

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

  • Solve the problem by calling solve(prob).

Solve ODEs using DifferentialEquations.jl#

Documentation: https://docs.sciml.ai/DiffEqDocs/stable/

Exponential decay model#

The concentration of a decaying nuclear isotope could be described as an exponential decay:

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

State variable

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

Parameter

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

using DifferentialEquations
using Plots
Plots.default(linewidth=2)

The model function is the 3-argument out-of-place form, f(u, p, t).

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

p = -1.0            ## Rate of exponential decay
u0 = 1.0            ## Initial condition
tspan = (0.0, 2.0)  ## Start time and end time

prob = ODEProblem(expdecay, u0, tspan, p)
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 8-element Vector{Float64}:
 0.0
 0.10001999200479662
 0.34208427066999536
 0.6553980136343391
 1.0312652525315806
 1.4709405856363595
 1.9659576669700232
 2.0
u: 8-element Vector{Float64}:
 1.0
 0.9048193287657775
 0.7102883621328676
 0.5192354400036404
 0.35655576576996556
 0.2297097907863828
 0.14002247272452764
 0.1353360028400881

Solution at time t=1.0 (with interpolation)

sol(1.0)

0.3678796381978344

Time points

sol.t

8-element Vector{Float64}:
 0.0
 0.10001999200479662
 0.34208427066999536
 0.6553980136343391
 1.0312652525315806
 1.4709405856363595
 1.9659576669700232
 2.0

Solutions at corresponding time points

sol.u

8-element Vector{Float64}:
 1.0
 0.9048193287657775
 0.7102883621328676
 0.5192354400036404
 0.35655576576996556
 0.2297097907863828
 0.14002247272452764
 0.1353360028400881

Visualize the solution

plot(sol)
_images/44cc398a05ab37fc9e7550294b532de1260fb823ca01f5fd06997563ee9deac5.png

SIR model#

This 3-variable model describes the spreading of an contagious disease can be described by the SIR model:

\[\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 variables

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

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

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

Parameters

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

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

using DifferentialEquations
using Plots
Plots.default(linewidth=2)

SIR model (in-place form)

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)

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

retcode: Success
Interpolation: 3rd order Hermite
t: 17-element Vector{Float64}:
  0.0
  0.08921318693905476
  0.3702862715172094
  0.7984257132319627
  1.3237271485666187
  1.991841832691831
  2.7923706947355837
  3.754781614278828
  4.901904318934307
  6.260476636498209
  7.7648912410433075
  9.39040980993922
 11.483861023017885
 13.372369854616487
 15.961357172044833
 18.681426667664056
 20.0
u: 17-element Vector{Vector{Float64}}:
 [0.99, 0.01, 0.0]
 [0.9890894703413342, 0.010634484617786016, 0.00027604504087978485]
 [0.9858331594901347, 0.012901496825852227, 0.0012653436840130785]
 [0.9795270529591532, 0.017282420996456258, 0.003190526044390597]
 [0.9689082167415561, 0.02463126703444545, 0.006460516223998508]
 [0.9490552312363142, 0.03827338797605378, 0.012671380787632141]
 [0.9118629475333939, 0.06347250098224964, 0.024664551484356558]
 [0.8398871089274511, 0.11078176031568547, 0.049331130756863524]
 [0.7075842068024722, 0.19166147882272844, 0.1007543143747994]
 [0.508146028721987, 0.29177419341470584, 0.20007977786330722]
 [0.31213222024413995, 0.3415879120018046, 0.34627986775405545]
 [0.18215683096365565, 0.3099983134156389, 0.5078448556207055]
 [0.10427205468919205, 0.22061114011133276, 0.6751168051994751]
 [0.07386737407725845, 0.14760143051851143, 0.7785311954042301]
 [0.05545028910907714, 0.07997076922865315, 0.8645789416622697]
 [0.047334990695892025, 0.04060565321383335, 0.9120593560902746]
 [0.04522885458929332, 0.029057416110814603, 0.925713729299892]

Visualize the solution

plot(sol, label=["S" "I" "R"], legend=:right)
_images/f086ae656fc169db7a012223b28e28e315f02c2a50684c56ba9e34949bd3c837.png

Saving simulation results#

using DataFrames
using CSV

df = DataFrame(sol)
CSV.write("lorenz.csv", df)
rm("lorenz.csv")

SIR model (MTK)#

using DifferentialEquations
using ModelingToolkit
using Plots
Plots.default(linewidth=2)

@independent_variables t
@parameters β γ
@variables s(t) i(t) r(t)
D = Differential(t)

eqs = [
    D(s) ~ -β * s * i,
    D(i) ~ β * s * i - γ * i,
    D(r) ~ γ * i
]

@mtkbuild sirSys = ODESystem(eqs, t)
\[\begin{split} \begin{align} \frac{\mathrm{d} s\left( t \right)}{\mathrm{d}t} &= - i\left( t \right) s\left( t \right) \beta \\ \frac{\mathrm{d} i\left( t \right)}{\mathrm{d}t} &= - i\left( t \right) \gamma + i\left( t \right) s\left( t \right) \beta \\ \frac{\mathrm{d} r\left( t \right)}{\mathrm{d}t} &= i\left( t \right) \gamma \end{align} \end{split}\]
p = [β => 1.0, γ => 0.3]
u0 = [s => 0.99, i => 0.01, r => 0.00]
tspan = (0.0, 20.0)

prob = ODEProblem(sirSys, u0, tspan, p)
sol = solve(prob)

plot(sol)
_images/ad0caf374533912e67ce5ffa361c806af8d60b23eadbb3fab9bdba577e7e19e6.png

Catalyst.jl#

Catalyst.jl is a domain-specific language (DSL) package to simulate chemical reaction networks.

using DifferentialEquations
using Catalyst
using Plots
Plots.default(linewidth=2)

Exponential decay model#

decay_rn = @reaction_network begin
    λ, C --> 0
end
\[ \begin{align*} \mathrm{C} &\xrightarrow{\lambda} \varnothing \end{align*} \]
p = [:λ => 1.0]
u0 = [:C => 1.0]
tspan = (0.0, 2.0)

prob = ODEProblem(decay_rn, u0, tspan, p)
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 8-element Vector{Float64}:
 0.0
 0.10001999200479662
 0.34208427066999536
 0.6553980136343391
 1.0312652525315806
 1.4709405856363595
 1.9659576669700232
 2.0
u: 8-element Vector{Vector{Float64}}:
 [1.0]
 [0.9048193287657775]
 [0.7102883621328676]
 [0.5192354400036404]
 [0.35655576576996556]
 [0.2297097907863828]
 [0.14002247272452764]
 [0.1353360028400881]

plot(sol, title="Exponential Decay")
_images/ecde3e543fe1a6f98adc79b6254d886b1fa0c4f353da69e6f4836053b469d261.png

SIR model#

sir_rn = @reaction_network begin
    β, S + I --> 2I
    γ, I --> R
end
\[\begin{split} \begin{align*} \mathrm{S} + \mathrm{I} &\xrightarrow{\beta} 2 \mathrm{I} \\ \mathrm{I} &\xrightarrow{\gamma} \mathrm{R} \end{align*} \end{split}\]

Extract the symbols for later use

@unpack β, γ, S, I, R = sir_rn

p = [β => 1.0, γ => 0.3]
u0 = [S => 0.99, I => 0.01, R => 0.00]
tspan = (0.0, 20.0)

prob = ODEProblem(sir_rn, u0, tspan, p)
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 17-element Vector{Float64}:
  0.0
  0.08921318693905476
  0.3702862715172094
  0.7984257132319627
  1.3237271485666187
  1.991841832691831
  2.7923706947355837
  3.754781614278828
  4.901904318934307
  6.260476636498209
  7.7648912410433075
  9.39040980993922
 11.483861023017885
 13.372369854616487
 15.961357172044833
 18.681426667664056
 20.0
u: 17-element Vector{Vector{Float64}}:
 [0.99, 0.01, 0.0]
 [0.9890894703413342, 0.010634484617786016, 0.00027604504087978485]
 [0.9858331594901347, 0.012901496825852227, 0.0012653436840130785]
 [0.9795270529591532, 0.017282420996456258, 0.003190526044390597]
 [0.9689082167415561, 0.02463126703444545, 0.006460516223998508]
 [0.9490552312363142, 0.03827338797605378, 0.012671380787632141]
 [0.9118629475333939, 0.06347250098224964, 0.024664551484356558]
 [0.8398871089274511, 0.11078176031568547, 0.049331130756863524]
 [0.7075842068024722, 0.19166147882272844, 0.1007543143747994]
 [0.508146028721987, 0.29177419341470584, 0.20007977786330722]
 [0.31213222024413995, 0.3415879120018046, 0.34627986775405545]
 [0.18215683096365565, 0.3099983134156389, 0.5078448556207055]
 [0.10427205468919205, 0.22061114011133276, 0.6751168051994751]
 [0.07386737407725845, 0.14760143051851143, 0.7785311954042301]
 [0.05545028910907714, 0.07997076922865315, 0.8645789416622697]
 [0.047334990695892025, 0.04060565321383335, 0.9120593560902746]
 [0.04522885458929332, 0.029057416110814603, 0.925713729299892]

plot(sol, legend=:right, title="SIR Model")
_images/21028590b3075a72ae9addf3a9df95641ea23608e05125c3a1f8861b7b4750fd.png

This notebook was generated using Literate.jl.