Fisher KPP with Julia

This implementation of numerical solve of a reaction diffusion equation is based on the presentation of the package MethodOfLines.jl at JuliaCon 2022 by A. Jones.

Fisher KPP equation

The Fisher KPP equation (Fisher’s version) reads (Fisher (1937), Kolmogorov, Petrovskii, and Piskunov (1937)):

\frac{\partial u}{\partial t} = ru\left(1-u\right) + D \frac{\partial^2 u}{\partial x^2},

with u(t,x) the population density at time t and position x (scaled to the local carrying capacity K), r the intrinsic growth rate of the population, and D the diffusion coefficient.

Packages

Let us first import the packages used for the simulation:

using MethodOfLines
using ModelingToolkit
using DomainSets
using OrdinaryDiffEq
using Plots
using LaTeXStrings

Model definition

MethodsOfLines.jl makes use of ModelingToolkit.jl to symbolically define the model to integrate.

Let us first define the time and space parameters:

@parameters t x

The model parameters:

@parameters r D

Now the variable u(t,x):

@variables u(..)

And finally the derivatives:

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

We can now define the model symbolically through:

eq = Dt(u(t, x)) ~ r * u(t,x) * (1-u(t,x)) + D * Dxx(u(t,x))

Differential(t)(u(t, x)) ~ D*Differential(x)(Differential(x)(u(t, x))) + r*u(t, x)*(1 - u(t, x))

Domains of integration

Let us introduce some parameters for space and time domains:

x_max = 30.0
t_max = 14.0

And the domains of integration:

domain = [x ∈ Interval(0.0, x_max),
          t ∈ Interval(0.0, t_max)]

We also introduce (initial and) boundary conditions:

ic_bc = [u(0.0, x) ~ 0.0,
         u(t, 0.0) ~ 1.0,
         u(t, x_max) ~ 0.0]

Simulation

We define the model to be integrated as a PDESystem, from the equation eq, the initial and boundary conditions ic_bc, the domains of integration domain, the time and space parameters t and x, the solution we want to retrieve u(t,x), and the model parameters r and D:

@named sys = PDESystem(eq, ic_bc, domain, [t, x], [u(t,x)], [r => 1.0, D => 1.0])

We set up the discretization of space, through MethodOfLines.jl:

dx = 0.1
discretization = MOLFiniteDifference([x => dx], t)

And we set up the (ODE) problem to be integrated:

prob = discretize(sys, discretization)

And we finally integrate it through the OrdinaryDiffEq.jl solver with Tsit5 algorithm.

sol = solve(prob, Tsit5(), saveat = .1)

Graphical representation

We retrieve the components of the solution for easier manipulation:

gridx = sol[x]
gridt = sol[t]
solu = sol[u(t,x)]

And we plot the animation of the solution through time:

anim = @animate for i in eachindex(gridt)
    plot(
        gridx,
        solu[i, :];
        xlabel = "position "*L"$x$",
        ylabel = "population density "*L"$u$",
        label = L"$u(x,t)$",
        title = "t=$(gridt[i])",
    )
end

gif(anim, "fisherKPP.gif", fps = 10)

[ Info: Saved animation to /home/ludo/Nextcloud/MyDrive/Programmes/julia/biomaths_julia_www/fisherKPP.gif

And that’s it !

References

Fisher, Ronald Aylmer. 1937. “The Wave of Advance of Advantageous Genes.” Annals of Eugenics 7 (4): 355–69.

Kolmogorov, A. N., I. G. Petrovskii, and N. S. Piskunov. 1937. “A Study of the Diffusion Equation with Increase in the Amount of Substance, and Its Application to a Biological Problem.” Bull. Moscow Univ. Math. Mech. 1:6: 1–26.

Reuse

CC BY-NC 4.0

--- title: "Fisher KPP with Julia" engine: julia ---  ```{julia} #| include: false #| eval: true # # this cell does not appear in the rendering, but is executed # # for reproducibility purposes, we load the local julia environment/project, containing # [13f3f980] CairoMakie v0.11.6 # [a93c6f00] DataFrames v1.6.1 # [864edb3b] DataStructures v0.18.16 # [0c46a032] DifferentialEquations v7.12.0 # [5b8099bc] DomainSets v0.6.7 # [e9467ef8] GLMakie v0.9.6 # [b964fa9f] LaTeXStrings v1.3.1 # [94925ecb] MethodOfLines v0.10.6 # [961ee093] ModelingToolkit v8.75.0 # [1dea7af3] OrdinaryDiffEq v6.69.0 # [91a5bcdd] Plots v1.40.0 # [f27b6e38] Polynomials v4.0.6 # [90137ffa] StaticArrays v1.9.1 # [0c5d862f] Symbolics v5.16.1 # # to share the environment, copy Project.toml and Manifest.toml files in some directory # `activate` the local environment # if necessary `instantiate` to get the correct package versions # # to check if some PackageName.jl is used from the local environment ## Pkg.status("PackageName") using Pkg Pkg.activate(".") ``` This implementation of numerical solve of a reaction diffusion equation is based on the [presentation](https://www.youtube.com/watch?v=8gLhaWRYvfQ) of the package `MethodOfLines.jl` at JuliaCon 2022 by A. Jones. ## Fisher KPP equation The Fisher KPP equation (Fisher's version) reads (@fisher1937, @Kolmogorov1937): $$ \frac{\partial u}{\partial t} = ru\left(1-u\right) + D \frac{\partial^2 u}{\partial x^2}, $$ with $u(t,x)$ the population density at time $t$ and position $x$ (scaled to the local carrying capacity $K$), $r$ the intrinsic growth rate of the population, and $D$ the diffusion coefficient. ## Packages Let us first import the packages used for the simulation: ```{julia} using MethodOfLines using ModelingToolkit using DomainSets using OrdinaryDiffEq using Plots using LaTeXStrings ``` ## Model definition `MethodsOfLines.jl` makes use of `ModelingToolkit.jl` to symbolically define the model to integrate. Let us first define the time and space parameters: ```{julia} @parameters t x ``` The model parameters: ```{julia} @parameters r D ``` Now the variable $u(t,x)$: ```{julia} @variables u(..) ``` And finally the derivatives: ```{julia} Dt = Differential(t) Dx = Differential(x) Dxx = Differential(x)^2 ``` We can now define the model symbolically through: ```{julia} #| output: true eq = Dt(u(t, x)) ~ r * u(t,x) * (1-u(t,x)) + D * Dxx(u(t,x)) ``` ## Domains of integration Let us introduce some parameters for space and time domains: ```{julia} x_max = 30.0 t_max = 14.0 ``` And the domains of integration: ```{julia} domain = [x ∈ Interval(0.0, x_max), t ∈ Interval(0.0, t_max)] ``` We also introduce (initial and) boundary conditions: ```{julia} ic_bc = [u(0.0, x) ~ 0.0, u(t, 0.0) ~ 1.0, u(t, x_max) ~ 0.0] ``` ## Simulation We define the model to be integrated as a `PDESystem`, from the equation `eq`, the initial and boundary conditions `ic_bc`, the domains of integration `domain`, the time and space parameters `t` and `x`, the solution we want to retrieve `u(t,x)`, and the model parameters $r$ and $D$: ```{julia} @named sys = PDESystem(eq, ic_bc, domain, [t, x], [u(t,x)], [r => 1.0, D => 1.0]) ``` We set up the discretization of space, through `MethodOfLines.jl`: ```{julia} dx = 0.1 discretization = MOLFiniteDifference([x => dx], t) ``` And we set up the (ODE) problem to be integrated: ```{julia} prob = discretize(sys, discretization) ``` And we finally integrate it through the `OrdinaryDiffEq.jl` solver with `Tsit5` algorithm. ```{julia} sol = solve(prob, Tsit5(), saveat = .1) ``` ## Graphical representation We retrieve the components of the solution for easier manipulation: ```{julia} gridx = sol[x] gridt = sol[t] solu = sol[u(t,x)] ``` And we plot the animation of the solution through time: ```{julia} #| output: true anim = @animate for i in eachindex(gridt) plot( gridx, solu[i, :]; xlabel = "position "*L"$x$", ylabel = "population density "*L"$u$", label = L"$u(x,t)$", title = "t=$(gridt[i])", ) end gif(anim, "fisherKPP.gif", fps = 10) ``` And that's it !