Let's solve the equation:
$$ \frac{du}{dt} = f(u,p,t) = \alpha u $$We know the solution would be:
$$ u(t) = u_0\exp(\alpha t) $$using DifferentialEquations
┌ Info: Precompiling DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa] └ @ Base loading.jl:1278
f(u,p,t) = 1.01*u #
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # flexible choice of algorithms
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 17-element Array{Float64,1}:
0.0
0.012407826196308189
0.042501278333560696
0.0817804940926822
0.12887384498704435
0.18409796286152927
0.24627457447456758
0.31479297816557983
0.3885963515160237
0.4668617724420117
0.5487161305960653
0.6334346972152323
0.7203630000154827
0.808957991167541
0.8987655040395068
0.9894161889652783
1.0
u: 17-element Array{Float64,1}:
0.5
0.5063053789114713
0.5219304750950854
0.5430527156531716
0.5695067765051011
0.6021743618001767
0.6412025634645868
0.6871475244605603
0.7403258398933249
0.8012223416244076
0.8702768595124217
0.9480214633315531
1.0350186537489083
1.1319031171663114
1.2393734610662304
1.3582039071131542
1.3728005076225747
using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
sol(0.45)
0.7876927465687832
plot(sol, lw = 2)
gui()
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
sol(0.45)
0.7876927465687832
Let's solve the following system of differential equations: \begin{align} \frac{dx}{dt} & = \sigma (y-x) \\ \frac{dy}{dx} & = x(\rho - z) - y \\ \frac{dz}{dt} & = xy - \beta z \end{align}
function lorenz!(du,u,p,t)
du[1] = 10.0*(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
lorenz! (generic function with 1 method)
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob);
plot(sol, vars = (1, 2, 3))
plot(sol, vars = (0,2))
plot(sol, vars = (3))
function parameterized_lorenz!(du,u,p,t)
du[1] = p[1]*(u[2]-u[1])
du[2] = u[1]*(p[2]-u[3]) - u[2]
du[3] = u[1]*u[2] - p[3]*u[3]
end
parameterized_lorenz! (generic function with 1 method)
function parameterized_lorenz!(du,u,p,t)
x,y,z = u
σ,ρ,β = p
du[1] = dx = σ*(y-x)
du[2] = dy = x*(ρ-z) - y
du[3] = dz = x*y - β*z
end
parameterized_lorenz! (generic function with 1 method)
u0 = [1.0,0.0,0.0]
tspan = (0.0,10.0)
p = [10.0, 28.0, 8/3]
prob = ODEProblem(parameterized_lorenz!,u0,tspan,p)
sol = solve(prob);
plot(sol, vars = (1,2,3))
Consider the following non-homogeneous equations:
\begin{align} \frac{d\theta(t)}{dt} & = \omega(t) \\ \frac{d\omega(t)}{dt} & = -\frac{3}{2}\frac{g}{l} sin \theta(t) +\frac{3}{ml^2}M(t) \end{align}l = 1.0
m = 1.0
g = 9.81
function pendulum!(du,u,p,t)
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
θ₀ = 0.01 # initial angular deflection [rad]
ω₀ = 0.0 # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀] # initial state vector
tspan = (0.0,10.0) # time interval
M = t->0.1sin(t) # external torque [Nm]
prob = ODEProblem(pendulum!,u₀,tspan,M)
sol = solve(prob)
plot(sol,linewidth=2,xaxis="t",label=["θ [rad]" "ω [rad/s]"],layout=(2,1))
A = [1. 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)
ODEProblem with uType Array{Float64,2} and tType Float64. In-place: false timespan: (0.0, 1.0) u0: 4×2 Array{Float64,2}: 0.613649 0.00364638 0.714742 0.743978 0.50977 0.950863 0.806768 0.212486
plot(sol, lw = 3)
The model is based on Andrew Atkeson (NBER w26876, 2020), Atkeson et al. (NBER w27335, 2020), Fernandez-Villaverde and Jones (NBER w27128, 2020)
using LinearAlgebra, Statistics, Random, SparseArrays
using OrdinaryDiffEq
using Parameters, Plots
gr(fmt=:svg);
In the version of the SEIR model, all individuals in the population are assumed to be in a finite number of states.
The states are: susceptible (S), exposed (E), infected (I) and removed (R).
This type of compartmentalized model has many extensions (e.g. SEIRS relaxes lifetime immunity and allow transitions from R→S).
Comments:
Those in state R have been infected and either recovered or died. Note that in some variations, R may refer only to recovered agents. Those who have recovered, and live, are assumed to have acquired immunity. Those in the exposed group are not yet infectious.
In the SEIR model, the flow across states follows the path $S-> E -> I -> R$ The number of individuals is of size $N$; we ignore birth and non-covid death. This means that
$$ S(t)+E(t)+I(t)+R(t) = N $$for all $t$.
We assume that $N$ is large, so we use continuum approximation for the number of individuals for each state.
Parameters governing the transition rates:
The full SEIR model is characterized by following equations:
\begin{align} \frac{dS}{dt} &= -\beta S \frac{I}{N} \\ \tag{1} \frac{dE}{dt} &= \beta S \frac{I}{N} -\sigma E \\ \frac{dI}{dt} & = \sigma E -\gamma I \\ \frac{dR}{dt} & = \gamma I \end{align}We could define $R_0(t):=\beta (t)/\gamma$. This is the famous basic reproduction number for the SEIR model:
Prior to solving the model, we make a few changes:
The "removed" fraction of the population is $r = 1-s -e -i$
We begin by implementing a simple version of this model with a constant $R_0$ and some baseline parameter values (which we discuss later).
function F_simple(x, p, t; γ = 1/18, R₀ = 3.0, σ = 1/5.2)
s, e, i, r = x
return [-γ*R₀*s*i; # ds/dt = -γR₀si
γ*R₀*s*i - σ*e;# de/dt = γR₀si -σe
σ*e - γ*i; # di/dt = σe -γi
γ*i; # dr/dt = γi
]
end
F_simple (generic function with 1 method)
function sir_model(du, u, p, t)
s,e,i,r = u
γ,R₀,σ = p
du[1] = ds = -γ*R₀*s*i
du[2] = de = γ*R₀*s*i - σ*e
du[3] = di = σ*e - γ*i
du[4] = dr = γ*i
end
sir_model (generic function with 1 method)
i_0 = 1E-7 # 33 = 1E-7 * 330 million population = initially infected
e_0 = 4.0 * i_0 # 132 = 1E-7 *330 million = initially exposed
s_0 = 1.0 - i_0 - e_0
r_0 = 0.0
x_0 = [s_0, e_0, i_0, r_0] # initial condition
tspan = (0.0, 350.0) # ≈ 350 days
prob = ODEProblem(F_simple, x_0, tspan)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: false timespan: (0.0, 350.0) u0: 4-element Array{Float64,1}: 0.9999995 4.0e-7 1.0e-7 0.0
p = [1/18, 3.0, 1/5.2]
prob1 = ODEProblem(sir_model, x_0, tspan, p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 350.0) u0: 4-element Array{Float64,1}: 0.9999995 4.0e-7 1.0e-7 0.0
sol = solve(prob, Tsit5())
plot(sol, labels = ["s" "e" "i" "r"], title = "SEIR Dynamics", lw = 3)
sol1 = solve(prob1, Tsit5())
plot(sol1, labels = ["s" "e" "i" "r"], title = "SEIR Dynamics", lw = 3)
#visualize the proportion in each state
areaplot(sol.t, sol', labels = ["s" "e" "i" "r"], title = "SEIR Proportions")
Cumulative case load is defined as $c(t) = i(t) + r(t)$ and
$$ \frac{dc}{dt} = \sigma e $$We assume that the transmission rate $R_0(t)$ follows:
$$ \frac{dR_0}{dt} = \eta (\bar{R}_0-R_0) $$Define the cumulative number of deaths as $D(t)$ and the proportion of death is $d(t):=D(t)/N$
$$ \frac{d}{dt}d(t) = \delta \gamma i $$The system of equations for the model can be re-written accordingly with 3 additional states.
function F(u, p, t)
s, e, i, r, R₀, c, d = u
@unpack σ, γ, R̄₀, η, δ = p
return [-γ*R₀*s*i; # ds/dt
γ*R₀*s*i - σ*e; # de/dt
σ*e - γ*i; # di/dt
γ*i; # dr/dt
η*(R̄₀(t, p) - R₀);# dR₀/dt
σ*e; # dc/dt (cumulative cases)
δ*γ*i; # dd/dt
]
end;
p_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20,
R₀_n = 1.6, δ = 0.01, N = 3.3E8,
R̄₀ = (t, p) -> p.R₀_n); # use named tuples
p = p_gen() # use all default parameters
i_0 = 1E-7
e_0 = 4.0 * i_0
s_0 = 1.0 - i_0 - e_0
x_0 = [s_0, e_0, i_0, 0.0, p.R₀_n, 0.0, 0.0]
tspan = (0.0, p.T)
prob = ODEProblem(F, x_0, tspan, p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: false timespan: (0.0, 550.0) u0: 7-element Array{Float64,1}: 0.9999995 4.0e-7 1.0e-7 0.0 1.6 0.0 0.0
sol = solve(prob, Tsit5())
@show length(sol.t);
length(sol.t) = 45
While it may seem that 45 time intervals is extremely small for that range, for much of the t, the functions are very flat - and hence adaptive time-stepping algorithms can move quickly and interpolate accurately.
A few more comments:
If you want to ensure that there are specific points that the adaptive-time stepping must include (e.g. at known discontinuities) use tstops.
The built-in plots for the solutions provide all of the attributes in Plots.jl.
See here for details on analyzing the solution and extracting the output.
plot(sol, vars = [6, 7], label = ["c(t)" "d(t)"], lw = 2,
title = ["Cumulative Infected" "Death Proportion"],
layout = (1,2), size = (900, 300))
$\bar{R}_0(t) = R_{0n}$ is constant.
R₀_n_vals = range(1.6, 3.0, length = 6)
sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R₀_n = R₀_n)),
Tsit5(), saveat=0.5) for R₀_n in R₀_n_vals]; #saveat = 0.5 to get evenly spaced every 0.5
labels = permutedims(["R_0 = $r" for r in R₀_n_vals])
infecteds = [sol[3,:] for sol in sols]
plot(infecteds, label=labels, legend=:topleft, lw = 3, xlabel = "t",
ylabel = "i(t)", title = "Current Cases")
cumulative_infected = [sol[6,:] for sol in sols]
plot(cumulative_infected, label=labels ,legend=:topleft, lw = 2, xlabel = "t",
ylabel = "c(t)", title = "Cumulative Cases")
Let’s look at a scenario where mitigation (e.g., social distancing) is successively imposed, but the target (maintaining $R_{0n}$) is fixed.
We will examine the case where $R_0(0)=3$ and then it falls to $R_{0n}=1.6$ due to the progressive adoption of stricter mitigation measures.
η_vals = [1/5, 1/10, 1/20, 1/50, 1/100]
labels = permutedims(["eta = $η" for η in η_vals]);
x_0 = [s_0, e_0, i_0, 0.0, 3.0, 0.0, 0.0]
sols = [solve(ODEProblem(F, x_0, tspan, p_gen(η=η)), Tsit5(), saveat=0.5) for η in η_vals];
Rs = [sol[5,:] for sol in sols]
plot(Rs, label=labels, legend=:topright, lw = 2, xlabel = "t",
ylabel = "R_0(t)", title = "Basic Reproduction Number")
infecteds = [sol[3,:] for sol in sols]
plot(infecteds, label=labels, legend=:topleft, lw = 2, xlabel = "t",
ylabel = "i(t)", title = "Current Infected")
cumulative_infected = [sol[6,:] for sol in sols]
plot(cumulative_infected, label=labels ,legend=:topleft, lw = 2, xlabel = "t",
ylabel = "c(t)", title = "Cumulative Infected")
Consider these two mitigation scenarios:
For both of these, we will choose a large η to focus on the case where rapid changes in the lockdown policy remain feasible.
R₀_L = 0.5 # lockdown
R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0 #step function
R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0 #step function
p_early = p_gen(R̄₀= R̄₀_lift_early, η = 10.0)
p_late = p_gen(R̄₀= R̄₀_lift_late, η = 10.0)
# initial conditions
i_0 = 25000 / p_early.N
e_0 = 75000 / p_early.N
s_0 = 1.0 - i_0 - e_0
x_0 = [s_0, e_0, i_0, 0.0, R₀_L, 0.0, 0.0] # start in lockdown
# create two problems, with rapid movement of R₀(t) towards R̄₀(t)
prob_early = ODEProblem(F, x_0, tspan, p_early)
prob_late = ODEProblem(F, x_0, tspan, p_late)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: false timespan: (0.0, 550.0) u0: 7-element Array{Float64,1}: 0.9996969696969696 0.00022727272727272727 7.575757575757576e-5 0.0 0.5 0.0 0.0
sol_early = solve(prob_early, Tsit5(), tstops = [30.0, 120.0])
sol_late = solve(prob_late, Tsit5(), tstops = [30.0, 120.0])
plot(sol_early, vars = [7], title = "Total Mortality", label = "Lift Early", legend = :topleft)
plot!(sol_late, vars = [7], label = "Lift Late")
# Daily Deaths
flow_deaths(sol, p) = p.N * p.δ * p.γ * sol[3,:]
plot(sol_early.t, flow_deaths(sol_early, p_early), title = "Flow Deaths", label = "Lift Early")
plot!(sol_late.t, flow_deaths(sol_late, p_late), label = "Lift Late")