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

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

@time prob, t, x = sis_model()

 69.001844 seconds (89.91 M allocations: 5.064 GiB, 1.76% gc time, 92.34% compilation time: 9% of which was recompilation)

(prob = SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ModelingToolkitBase.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xd34fbd62, 0x7e1fa66b, 0xbeb7e4c5, 0x4c3a5514, 0x9aee684e), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xb4cac23a, 0x3112a02d, 0x41ec0037, 0x99b3038d, 0x1a8f5029), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Union{Nothing, SciMLBase.OverrideInitData}, Union{Nothing, SciMLBase.ODENLStepData}}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, ModelingToolkitBase.PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}(SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xd34fbd62, 0x7e1fa66b, 0xbeb7e4c5, 0x4c3a5514, 0x9aee684e), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xb4cac23a, 0x3112a02d, 0x41ec0037, 0x99b3038d, 0x1a8f5029), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Union{Nothing, SciMLBase.OverrideInitData}, Union{Nothing, SciMLBase.ODENLStepData}}(ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xd34fbd62, 0x7e1fa66b, 0xbeb7e4c5, 0x4c3a5514, 0x9aee684e), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xb4cac23a, 0x3112a02d, 0x41ec0037, 0x99b3038d, 0x1a8f5029), Nothing}}(RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xd34fbd62, 0x7e1fa66b, 0xbeb7e4c5, 0x4c3a5514, 0x9aee684e), Nothing}(nothing), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xb4cac23a, 0x3112a02d, 0x41ec0037, 0x99b3038d, 0x1a8f5029), Nothing}(nothing)), LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}(Model pdesys:
Equations (198):
  198 standard: see equations(pdesys)
Unknowns (198): see unknowns(pdesys)
  (I(t))[100]
  (I(t))[99]
  (I(t))[98]
  (I(t))[97]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ
Observed (6): see observed(pdesys), Dict{Any, Any}(), false, false, ModelingToolkitBase, false, true, nothing), nothing, Model pdesys:
Equations (198):
  198 standard: see equations(pdesys)
Unknowns (198): see unknowns(pdesys)
  (I(t))[100]
  (I(t))[99]
  (I(t))[98]
  (I(t))[97]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ
Observed (6): see observed(pdesys), nothing, nothing), [0.19980267284282716, 0.1992114701314478, 0.19822872507286887, 0.1968583161128631, 0.19510565162951538, 0.19297764858882513, 0.190482705246602, 0.18763066800438638, 0.18443279255020154, 0.18090169943749476  …  0.9587785252292473, 0.9535826794978997, 0.9481753674101716, 0.9425779291565073, 0.9368124552684678, 0.9309016994374948, 0.9248689887164855, 0.9187381314585725, 0.9125333233564304, 0.9062790519529313], (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, Union{}, Tuple{}, @NamedTuple{}}(), MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, ModelingToolkitBase.PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}(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}}(x => (0.0, 1.0), t => (0.0, 10.0)), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, ReadOnlyArrays.ReadOnlyVector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.SmallVec{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{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, MethodOfLines.ScalarizedDiscretization}(Dict{Symbolics.Num, Float64}(x => 0.01), t, 2, MethodOfLines.UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}()), PDESystem
Equations: Symbolics.Equation[Differential(t, 1)(S(t, x)) - ((-(1 + x)*(brn + sinpi(2x)*ϵ)*I(t, x)*S(t, x)) / (I(t, x) + S(t, x))) - dS*Differential(x, 2)(S(t, x)) - (1 + x)*I(t, x) ~ 0, Differential(t, 1)(I(t, x)) - (((1 + x)*(brn + sinpi(2x)*ϵ)*I(t, x)*S(t, x)) / (I(t, x) + S(t, x))) - dI*Differential(x, 2)(I(t, x)) + (1 + x)*I(t, 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}}}(), false, Base.RefValue{Any}(Model pdesys:
Equations (202):
  202 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, Pair{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}[(S(t))[1] => 0.9 + 0.1sinpi(0.0), (S(t))[2] => 0.9 + 0.1sinpi(0.02), (S(t))[3] => 0.9 + 0.1sinpi(0.04), (S(t))[4] => 0.9 + 0.1sinpi(0.06), (S(t))[5] => 0.9 + 0.1sinpi(0.08), (S(t))[6] => 0.9 + 0.1sinpi(0.1), (S(t))[7] => 0.9 + 0.1sinpi(0.12), (S(t))[8] => 0.9 + 0.1sinpi(0.14), (S(t))[9] => 0.9 + 0.1sinpi(0.16), (S(t))[10] => 0.9 + 0.1sinpi(0.18)  …  (I(t))[92] => 0.1 + 0.1cospi(1.82), (I(t))[93] => 0.1 + 0.1cospi(1.84), (I(t))[94] => 0.1 + 0.1cospi(1.86), (I(t))[95] => 0.1 + 0.1cospi(1.88), (I(t))[96] => 0.1 + 0.1cospi(1.9), (I(t))[97] => 0.1 + 0.1cospi(1.92), (I(t))[98] => 0.1 + 0.1cospi(1.94), (I(t))[99] => 0.1 + 0.1cospi(1.96), (I(t))[100] => 0.1 + 0.1cospi(1.98), (I(t))[101] => 0.1 + 0.1cospi(2.0)])), t = t, x = x)

Solving time-dependent SIS epidemic model¶

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

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

 10.009446 seconds (12.77 M allocations: 803.065 MiB, 1.48% gc time, 99.30% compilation time: 3% 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.879621 0.879606 … 0.7900…
  I(t, x) => [0.2 0.199803 … 0.199803 0.2; 0.204072 0.203977 … 0.245383 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.203691  0.203211     0.2449    0.245383  0.245544
 0.239567  0.239532  0.239427  0.239248     0.320235  0.320608  0.320732
 0.300266  0.300276  0.300306  0.300352     0.413704  0.413992  0.414088
 0.376036  0.376073  0.376184  0.376363     0.503742  0.503967  0.504042
 0.452112  0.452158  0.452297  0.452525  …  0.573952  0.574128  0.574187
 0.516857  0.516901  0.517035  0.517254     0.620401  0.620539  0.620584
 0.566217  0.566254  0.566367  0.566552     0.647111  0.647219  0.647254
 0.601153  0.601183  0.601272  0.601419     0.660266  0.660352  0.66038
 0.624956  0.624979  0.625049  0.625163     0.665592  0.665662  0.665685
 ⋮                                       ⋱                      ⋮
 0.67636   0.676364  0.676374  0.676391     0.656473  0.656498  0.656507
 0.676357  0.676361  0.676372  0.676388     0.656475  0.6565    0.656509
 0.676354  0.676358  0.676369  0.676385     0.656478  0.656503  0.656511
 0.676351  0.676354  0.676365  0.676382  …  0.65648   0.656506  0.656514
 0.676347  0.67635   0.676361  0.676378     0.656483  0.656509  0.656517
 0.676343  0.676346  0.676357  0.676374     0.656487  0.656512  0.656521
 0.676338  0.676342  0.676352  0.676369     0.656491  0.656516  0.656524
 0.676332  0.676336  0.676347  0.676363     0.656495  0.656521  0.656529
 0.676325  0.676329  0.67634   0.676357  …  0.656501  0.656526  0.656534

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.