1D PDE: SIS diffusion model

Source

\begin{align} \frac{\partial S}{\partial t} &= d_S S_{xx} - \beta(x)\frac{SI}{S+I} + \gamma(x)I \\ \frac{\partial I}{\partial t} &= d_I I_{xx} + \beta(x)\frac{SI}{S+I} - \gamma(x)I \end{align}

(1)

where $x \in (0, 1)$

Solve the steady-state problem $\frac{\partial S}{\partial t} = \frac{\partial I}{\partial t} = 0$

The boundary condition: $\frac{\partial S}{\partial x} = \frac{\partial I}{\partial x} = 0$ for x = 0, 1

The conservative relationship: $\int^{1}_{0} (S(x) + I(x) ) dx = 1$

Notations:

$x$ : location
$t$ : time
$S(x, t)$ : the density of susceptible populations
$I(x, t)$ : the density of infected populations
$d_S$ / $d_I$ : the diffusion coefficients for susceptible and infected individuals
$\beta(x)$ : transmission rates
$\gamma(x)$ : recovery rates

using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots
using OrdinaryDiffEqBDF

Setup parameters, variables, and differential operators

function sis_model(;dx = 0.01, order = 2)
    @independent_variables t x
    @parameters dS=0.5 dI=0.1 brn=3 ϵ=0.1
    @variables S(..) I(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2
    # Helper functions
    γ(x) = x + 1
    ratio(x, brn, ϵ) = brn + ϵ * sinpi(2x)

    # 1D PDE for disease spreading
    eqs = [
        Dt(S(t, x)) ~ dS * Dxx(S(t, x)) - ratio(x, brn, ϵ) * γ(x) * S(t, x) * I(t, x) / (S(t, x) + I(t, x)) + γ(x) * I(t, x),
        Dt(I(t, x)) ~ dI * Dxx(I(t, x)) + ratio(x, brn, ϵ) * γ(x) * S(t, x) * I(t, x) / (S(t, x) + I(t, x)) - γ(x) * I(t, x)
    ]

    # Boundary conditions (including initial conditions)
    bcs = [
        S(0, x) ~ 0.9 + 0.1 * sinpi(2x),
        I(0, x) ~ 0.1 + 0.1 * cospi(2x),
        Dx(S(t, 0)) ~ 0.0,
        Dx(S(t, 1)) ~ 0.0,
        Dx(I(t, 0)) ~ 0.0,
        Dx(I(t, 1)) ~ 0.0
    ]

    # Space and time domains
    domains = [
        t ∈ Interval(0.0, 10.0),
        x ∈ Interval(0.0, 1.0)
    ]

    # Build the PDE system
    @named pdesys = PDESystem(eqs, bcs, domains,
        [t, x], ## Independent variables
        [S(t, x), I(t, x)],  ## Dependent variables
        [dS, dI, brn, ϵ],    ## parameters
    )
    # Finite difference method (FDM) converts the PDE system into an ODE problem
    discretization = MOLFiniteDifference([x => dx], t, approx_order=order)
    prob = discretize(pdesys, discretization)
    return (; prob, t, x)
end

sis_model (generic function with 1 method)

There are 202 ODEs in the system.

@time "Build problem" prob, t, x = sis_model()

Build problem: 73.719513 seconds (122.90 M allocations: 6.635 GiB, 2.42% gc time, 99.78% compilation time: 7% of which was recompilation)

(prob = SciMLBase.DAEProblem{Vector{Float64}, Vector{Float64}, Tuple{Float64, Float64}, true, ModelingToolkitBase.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.DAEFunction{true, SciMLBase.AutoDespecialize, ModelingToolkitBase.GeneratedFunctionWrapper{Tuple{3, 4, true}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x6813f0e4, 0xfa94abce, 0xf5d7826a, 0x6e2e00c9, 0x93afcc25), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1187374493630987813, :__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xdc45730f, 0xdbd0a6bf, 0x0717304a, 0x9b048fa2, 0x3409b121), Nothing}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Nothing, Nothing}, Base.Pairs{Symbol, DiffEqBase.BrownFullBasicInit{Float64, Nothing}, Nothing, @NamedTuple{initializealg::DiffEqBase.BrownFullBasicInit{Float64, Nothing}}}, Vector{Bool}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}, ModelingToolkitBase.PDESystem, Base.RefValue{Any}, Nothing}}(SciMLBase.DAEFunction{true, SciMLBase.AutoDespecialize, ModelingToolkitBase.GeneratedFunctionWrapper{Tuple{3, 4, true}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x6813f0e4, 0xfa94abce, 0xf5d7826a, 0x6e2e00c9, 0x93afcc25), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1187374493630987813, :__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xdc45730f, 0xdbd0a6bf, 0x0717304a, 0x9b048fa2, 0x3409b121), Nothing}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Nothing, Nothing}(ModelingToolkitBase.GeneratedFunctionWrapper{Tuple{3, 4, true}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x6813f0e4, 0xfa94abce, 0xf5d7826a, 0x6e2e00c9, 0x93afcc25), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1187374493630987813, :__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xdc45730f, 0xdbd0a6bf, 0x0717304a, 0x9b048fa2, 0x3409b121), Nothing}}(RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x6813f0e4, 0xfa94abce, 0xf5d7826a, 0x6e2e00c9, 0x93afcc25), Nothing}(nothing), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1187374493630987813, :__mtk_arg_1, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1565065830599138339), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xdc45730f, 0xdbd0a6bf, 0x0717304a, 0x9b048fa2, 0x3409b121), Nothing}(nothing)), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}(Model pdesys:
Equations (6):
  6 standard: see equations(pdesys)
Unknowns (202): see unknowns(pdesys)
  (S(t))[1]
  (S(t))[2]
  (S(t))[3]
  (S(t))[4]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ, Dict{Any, Any}(), false, false, ModelingToolkitBase, false, nothing), nothing, Model pdesys:
Equations (6):
  6 standard: see equations(pdesys)
Unknowns (202): see unknowns(pdesys)
  (S(t))[1]
  (S(t))[2]
  (S(t))[3]
  (S(t))[4]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ, nothing, nothing), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.9, 0.9062790519529313, 0.9125333233564304, 0.9187381314585725, 0.9248689887164855, 0.9309016994374948, 0.9368124552684678, 0.9425779291565073, 0.9481753674101716, 0.9535826794978997  …  0.18443279255020154, 0.18763066800438638, 0.190482705246602, 0.19297764858882513, 0.19510565162951538, 0.1968583161128631, 0.19822872507286887, 0.1992114701314478, 0.19980267284282716, 0.2], (0.0, 10.0), ModelingToolkitBase.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}([0.5, 0.1, 3.0, 0.1], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], (), (), (), ()), Base.Pairs{Symbol, DiffEqBase.BrownFullBasicInit{Float64, Nothing}, Nothing, @NamedTuple{initializealg::DiffEqBase.BrownFullBasicInit{Float64, Nothing}}}(:initializealg => DiffEqBase.BrownFullBasicInit{Float64, Nothing}(1.0e-10, nothing)), Bool[0, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1, 1, 1, 1, 0], MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}, ModelingToolkitBase.PDESystem, Base.RefValue{Any}, Nothing}(MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}(PDEBase.VariableMap(Any[S(t, x), I(t, x)], SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}[x], Symbolics.Num[dS, dI, brn, ϵ], t, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Tuple{Float64, Float64}}(t => (0.0, 10.0), x => (0.0, 1.0)), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, ReadOnlyArrays.ReadOnlyVector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.ArgsT{SymbolicUtils.SymReal}}}(I => [t, x], S => [t, x]), SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}[S, I], Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Int64}(x => 1), Dict{Int64, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}(1 => x), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Symbolics.Num}()), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}}(S(t, x) => [(S(t))[1], (S(t))[2], (S(t))[3], (S(t))[4], (S(t))[5], (S(t))[6], (S(t))[7], (S(t))[8], (S(t))[9], (S(t))[10]  …  (S(t))[92], (S(t))[93], (S(t))[94], (S(t))[95], (S(t))[96], (S(t))[97], (S(t))[98], (S(t))[99], (S(t))[100], (S(t))[101]], I(t, x) => [(I(t))[1], (I(t))[2], (I(t))[3], (I(t))[4], (I(t))[5], (I(t))[6], (I(t))[7], (I(t))[8], (I(t))[9], (I(t))[10]  …  (I(t))[92], (I(t))[93], (I(t))[94], (I(t))[95], (I(t))[96], (I(t))[97], (I(t))[98], (I(t))[99], (I(t))[100], (I(t))[101]]), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(x => 0.0:0.01:1.0), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(x => 0.0:0.01:1.0), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Float64}(x => 0.01), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}(S(t, x) => CartesianIndices((101,)), I(t, x) => CartesianIndices((101,))), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}(S(t, x) => CartesianIndices((101,)), I(t, x) => CartesianIndices((101,))), nothing), MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}(Dict{Symbolics.Num, Float64}(x => 0.01), t, 2, MethodOfLines.UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, true, Any[], Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}()), PDESystem
Equations: Symbolics.Equation[Differential(t, 1)(S(t, x)) - ((-I(t, x)*S(t, x)*(brn + sinpi(2x)*ϵ)*(1 + x)) / (I(t, x) + S(t, x))) - Differential(x, 2)(S(t, x))*dS - I(t, x)*(1 + x) ~ 0, Differential(t, 1)(I(t, x)) - ((I(t, x)*S(t, x)*(brn + sinpi(2x)*ϵ)*(1 + x)) / (I(t, x) + S(t, x))) - Differential(x, 2)(I(t, x))*dI + I(t, x)*(1 + x) ~ 0]
Boundary Conditions: Symbolics.Equation[S(0, x) ~ 0.9 + 0.1sinpi(2x), I(0, x) ~ 0.1 + 0.1cospi(2x), Differential(x, 1)(S(t, 0)) ~ 0.0, Differential(x, 1)(S(t, 1)) ~ 0.0, Differential(x, 1)(I(t, 0)) ~ 0.0, Differential(x, 1)(I(t, 1)) ~ 0.0]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0 .. 10.0), Symbolics.VarDomainPairing(x, 0.0 .. 1.0)]
Dependent Variables: Symbolics.Num[S(t, x), I(t, x)]
Independent Variables: Symbolics.Num[t, x]
Parameters: Symbolics.Num[dS, dI, brn, ϵ]
Default Parameter ValuesModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}}(), Base.RefValue{Any}(nothing), nothing, Pair{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Float64}[(S(t))[1] => 0.9, (S(t))[2] => 0.9062790519529313, (S(t))[3] => 0.9125333233564304, (S(t))[4] => 0.9187381314585725, (S(t))[5] => 0.9248689887164855, (S(t))[6] => 0.9309016994374948, (S(t))[7] => 0.9368124552684678, (S(t))[8] => 0.9425779291565073, (S(t))[9] => 0.9481753674101716, (S(t))[10] => 0.9535826794978997  …  (I(t))[92] => 0.18443279255020154, (I(t))[93] => 0.18763066800438638, (I(t))[94] => 0.190482705246602, (I(t))[95] => 0.19297764858882513, (I(t))[96] => 0.19510565162951538, (I(t))[97] => 0.1968583161128631, (I(t))[98] => 0.19822872507286887, (I(t))[99] => 0.1992114701314478, (I(t))[100] => 0.19980267284282716, (I(t))[101] => 0.2])), t = t, x = x)

Solving time-dependent SIS epidemic model¶

KenCarp47 and FBDF are good at solving reaction-diffusion problems.

@time sol = solve(prob, DFBDF(), saveat=0.2)

 14.258168 seconds (18.49 M allocations: 1.033 GiB, 1.62% gc time, 98.94% compilation time: <1% of which was recompilation)

retcode: Success
Interpolation: Dict{Symbolics.Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}
t: 51-element Vector{Float64}:
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0
  1.2
  1.4
  1.6
  1.8
  ⋮
  8.4
  8.6
  8.8
  9.0
  9.2
  9.4
  9.6
  9.8
 10.0ivs: 2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}:
 t
 xdomain:([0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8  …  8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, 10.0], 0.0:0.01:1.0)
u: Dict{Symbolics.Num, Matrix{Float64}} with 2 entries:
  S(t, x) => [0.904194 0.906279 … 0.893721 0.895806; 0.879549 0.879535 … 0.7900…
  I(t, x) => [0.2 0.199803 … 0.199803 0.2; 0.204072 0.203977 … 0.245385 0.24554…

Grid points

discrete_x = sol[x]
discrete_t = sol[t]

51-element Vector{Float64}:
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0
  1.2
  1.4
  1.6
  1.8
  ⋮
  8.4
  8.6
  8.8
  9.0
  9.2
  9.4
  9.6
  9.8
 10.0

Results (Matrices)

@variables S(..) I(..)
sol
S_solution = sol[S(t, x)]
I_solution = sol[I(t, x)]

51×101 Matrix{Float64}:
 0.2       0.199803  0.199211  0.198229  …  0.199211  0.199803  0.2
 0.204072  0.203977  0.20369   0.203211     0.244902  0.245385  0.245547
 0.23955   0.239515  0.23941   0.239232     0.32025   0.320622  0.320746
 0.300221  0.300231  0.300261  0.300307     0.413676  0.413965  0.414061
 0.375973  0.37601   0.37612   0.376301     0.503818  0.504044  0.504119
 0.452152  0.452198  0.452337  0.452565  …  0.574007  0.574183  0.574242
 0.517042  0.517086  0.517219  0.517437     0.620276  0.620413  0.620459
 0.566331  0.566369  0.566481  0.566665     0.646839  0.646946  0.646982
 0.601197  0.601226  0.601316  0.601462     0.66008   0.660165  0.660194
 0.624949  0.624972  0.625042  0.625156     0.665493  0.665562  0.665585
 ⋮                                       ⋱                      ⋮
 0.676329  0.676333  0.676344  0.67636      0.656498  0.656523  0.656531
 0.676329  0.676333  0.676344  0.676361     0.656498  0.656523  0.656531
 0.676329  0.676333  0.676343  0.67636      0.656498  0.656523  0.656532
 0.676328  0.676332  0.676343  0.676359  …  0.656499  0.656524  0.656532
 0.676327  0.676331  0.676341  0.676358     0.656499  0.656525  0.656533
 0.676325  0.676329  0.67634   0.676356     0.656501  0.656526  0.656535
 0.676323  0.676327  0.676338  0.676354     0.656502  0.656528  0.656536
 0.676321  0.676324  0.676335  0.676352     0.656505  0.65653   0.656538
 0.676318  0.676321  0.676332  0.676349  …  0.656507  0.656532  0.656541

Visualize the solution

surface(discrete_x, discrete_t, S_solution, xlabel="Location", ylabel="Time", title="Susceptible")

surface(discrete_x, discrete_t, I_solution, xlabel="Location", ylabel="Time", title="Infectious")

This notebook was generated using Literate.jl.