using MethodOfLines
using ModelingToolkit
using DomainSets
using OrdinaryDiffEq
using Plots
using LaTeXStrings
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:
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:
= Differential(t)
Dt = Differential(x)
Dx = Differential(x)^2 Dxx
We can now define the model symbolically through:
= Dt(u(t, x)) ~ r * u(t,x) * (1-u(t,x)) + D * Dxx(u(t,x)) eq
\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} u\left( t, x \right) = D \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( t, x \right) + r u\left( t, x \right) \left( 1 - u\left( t, x \right) \right) \end{equation}
Domains of integration
Let us introduce some parameters for space and time domains:
= 30.0
x_max = 14.0 t_max
And the domains of integration:
= [x ∈ Interval(0.0, x_max),
domain ∈ Interval(0.0, t_max)] t
We also introduce (initial and) boundary conditions:
= [u(0.0, x) ~ 0.0,
ic_bc 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
:
= 0.1
dx = MOLFiniteDifference([x => dx], t) discretization
And we set up the (ODE) problem to be integrated:
= discretize(sys, discretization) prob
And we finally integrate it through the OrdinaryDiffEq.jl
solver with Tsit5
algorithm.
= solve(prob, Tsit5(), saveat = .1) sol
Graphical representation
We retrieve the components of the solution for easier manipulation:
= sol[x]
gridx = sol[t]
gridt = sol[u(t,x)] solu
And we plot the animation of the solution through time:
= @animate for i in eachindex(gridt)
anim plot(
gridx,:];
solu[i, = "position "*L"$x$",
xlabel ylabel = "population density "*L"$u$",
= L"$u(x,t)$",
label title = "t=$(gridt[i])",
)end
gif(anim, "fisherKPP.gif", fps = 10)
And that’s it !