1D PDE: SIS diffusion model

1D PDE: SIS diffusion model#

Source

\[\begin{split} \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} \end{split}\]

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

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

Setup parameters, variables, and differential operators

@parameters t x
@parameters dS dI brn ϵ
@variables S(..) I(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

Differential(x) ∘ Differential(x)

Helper functions

γ(x) = x + 1
ratio(x, brn, ϵ) =  brn + ϵ * sinpi(2x)

ratio (generic function with 1 method)

1D PDE and boundary conditions

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)
]
\[\begin{split} \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} S\left( t, x \right) =& \frac{ - \left( 1 + x \right) S\left( t, x \right) I\left( t, x \right) \left( brn + \mathrm{sinpi}\left( 2 x \right) \epsilon \right)}{S\left( t, x \right) + I\left( t, x \right)} + dS \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} S\left( t, x \right) + \left( 1 + x \right) I\left( t, x \right) \\ \frac{\mathrm{d}}{\mathrm{d}t} I\left( t, x \right) =& \frac{\left( 1 + x \right) S\left( t, x \right) I\left( t, x \right) \left( brn + \mathrm{sinpi}\left( 2 x \right) \epsilon \right)}{S\left( t, x \right) + I\left( t, x \right)} + dI \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} I\left( t, x \right) - \left( 1 + x \right) I\left( t, x \right) \end{align} \end{split}\]

Boundary 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
]
\[\begin{split} \begin{align} S\left( 0, x \right) =& 0.9 + 0.1 \mathrm{sinpi}\left( 2 x \right) \\ I\left( 0, x \right) =& 0.1 + 0.1 \mathrm{cospi}\left( 2 x \right) \\ \frac{\mathrm{d}}{\mathrm{d}x} S\left( t, 0 \right) =& 0 \\ \frac{\mathrm{d}}{\mathrm{d}x} S\left( t, 1 \right) =& 0 \\ \frac{\mathrm{d}}{\mathrm{d}x} I\left( t, 0 \right) =& 0 \\ \frac{\mathrm{d}}{\mathrm{d}x} I\left( t, 1 \right) =& 0 \end{align} \end{split}\]

Space and time domains

domains = [
    t ∈ Interval(0.0, 10.0),
    x ∈ Interval(0.0, 1.0)
]

2-element Vector{Symbolics.VarDomainPairing}:
 Symbolics.VarDomainPairing(t, 0.0 .. 10.0)
 Symbolics.VarDomainPairing(x, 0.0 .. 1.0)

Define the PDE system

@named pdesys = PDESystem(eqs, bcs, domains,
    [t, x], ## Independent variables
    [S(t, x), I(t, x)],  ## Dependent variables
    [dS => 0.5, dI => 0.1, brn => 3, ϵ => 0.1] ## Initial conditions
)
\[\begin{split} \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} S\left( t, x \right) =& \frac{ - \left( 1 + x \right) S\left( t, x \right) I\left( t, x \right) \left( brn + \mathrm{sinpi}\left( 2 x \right) \epsilon \right)}{S\left( t, x \right) + I\left( t, x \right)} + dS \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} S\left( t, x \right) + \left( 1 + x \right) I\left( t, x \right) \\ \frac{\mathrm{d}}{\mathrm{d}t} I\left( t, x \right) =& \frac{\left( 1 + x \right) S\left( t, x \right) I\left( t, x \right) \left( brn + \mathrm{sinpi}\left( 2 x \right) \epsilon \right)}{S\left( t, x \right) + I\left( t, x \right)} + dI \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} I\left( t, x \right) - \left( 1 + x \right) I\left( t, x \right) \end{align} \end{split}\]

Finite difference method (FDM) converts the PDE system into an ODE problem

dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx], t, approx_order=order)
prob = discretize(pdesys, discretization)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 198-element Vector{Float64}:
 0.9062790519529313
 0.9125333233564304
 0.9187381314585725
 0.9248689887164855
 0.9309016994374948
 0.9368124552684678
 0.9425779291565073
 0.9481753674101716
 0.9535826794978997
 0.9587785252292473
 0.963742398974869
 0.9684547105928689
 0.9728968627421412
 ⋮
 0.17289686274214117
 0.17705132427757894
 0.18090169943749476
 0.18443279255020154
 0.18763066800438638
 0.190482705246602
 0.19297764858882513
 0.19510565162951538
 0.1968583161128631
 0.19822872507286887
 0.1992114701314478
 0.19980267284282716

Solving time-dependent SIS epidemic model#

KenCarp4 is good at solving semilinear problems (like this one).

sol = solve(prob, KenCarp4(), saveat=0.2)

retcode: Success
Interpolation: Dict{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
  2.0
  2.2
  2.4
  ⋮
  7.8
  8.0
  8.2
  8.4
  8.6
  8.8
  9.0
  9.2
  9.4
  9.6
  9.8
 10.0ivs: 2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 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{Num, Matrix{Float64}} with 2 entries:
  S(t, x) => [0.904194 0.906279 … 0.893721 0.895806; 0.879603 0.879585 … 0.7899…
  I(t, x) => [0.2 0.199803 … 0.199803 0.2; 0.204066 0.20397 … 0.245396 0.245557…

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
  2.0
  2.2
  2.4
  ⋮
  7.8
  8.0
  8.2
  8.4
  8.6
  8.8
  9.0
  9.2
  9.4
  9.6
  9.8
 10.0

Results (Matrices)

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.204066  0.20397   0.203684  0.203205     0.244912  0.245396  0.245557
 0.2396    0.239565  0.23946   0.239282     0.320343  0.320715  0.32084
 0.300271  0.300281  0.300311  0.300356     0.413564  0.413852  0.413948
 0.37595   0.375986  0.376097  0.376277     0.503725  0.50395   0.504025
 0.451575  0.451621  0.45176   0.451988  …  0.573497  0.573673  0.573731
 0.516821  0.516865  0.516998  0.517215     0.620087  0.620225  0.62027
 0.566316  0.566354  0.566466  0.56665      0.646898  0.647006  0.647042
 0.601147  0.601177  0.601266  0.601413     0.660184  0.66027   0.660299
 0.624959  0.624983  0.625052  0.625166     0.665474  0.665543  0.665567
 0.64117   0.641188  0.641243  0.641331  …  0.666954  0.667012  0.667032
 0.651984  0.651998  0.652041  0.65211      0.666264  0.666313  0.666329
 0.65916   0.659171  0.659205  0.65926      0.664656  0.664698  0.664713
 ⋮                                       ⋱                      ⋮
 0.676301  0.676305  0.676315  0.676332     0.656521  0.656547  0.656555
 0.676303  0.676306  0.676317  0.676334  …  0.65652   0.656545  0.656554
 0.676303  0.676307  0.676318  0.676335     0.656519  0.656545  0.656553
 0.676304  0.676307  0.676318  0.676335     0.656519  0.656544  0.656553
 0.676304  0.676307  0.676318  0.676335     0.656519  0.656544  0.656553
 0.676303  0.676307  0.676317  0.676334     0.656519  0.656544  0.656553
 0.676302  0.676306  0.676317  0.676333  …  0.656519  0.656545  0.656553
 0.676301  0.676305  0.676316  0.676333     0.65652   0.656545  0.656554
 0.676301  0.676304  0.676315  0.676332     0.65652   0.656546  0.656554
 0.6763    0.676304  0.676314  0.676331     0.656521  0.656546  0.656554
 0.6763    0.676303  0.676314  0.676331     0.656521  0.656546  0.656555
 0.6763    0.676304  0.676315  0.676331  …  0.656521  0.656546  0.656554

Visualize

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")

Runtime environment#

using Pkg
Pkg.status()

Status `/tmp/docs/pde/Project.toml`
⌅ [5b8099bc] DomainSets v0.6.7
  [94925ecb] MethodOfLines v0.11.0
  [961ee093] ModelingToolkit v9.15.0
  [8913a72c] NonlinearSolve v3.12.4
  [1dea7af3] OrdinaryDiffEq v6.80.1
  [91a5bcdd] Plots v1.40.4
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

using InteractiveUtils
InteractiveUtils.versioninfo()

Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 4 default, 0 interactive, 2 GC (on 4 virtual cores)
Environment:
  JULIA_CI = true
  JULIA_NUM_THREADS = auto
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_PATH = /usr/local/julia/
  JULIA_DEPOT_PATH = /srv/juliapkg/
Contents