Differential Equations and Dynamical Systems in Julia

Examples for Biomathematics

Show the code

# import packages
using Plots, LaTeXStrings, DifferentialEquations, OrdinaryDiffEq, ModelingToolkit

Example: A Single Species Population Model

Consider the logistic growth model for a single population:

\[ \frac{dN}{dt} = rN\left(1-\frac{N}{K}\right) \]

where \(N\) is the population size, \(r\) is the intrinsic growth rate, and \(K\) is the carrying capacity. It can be shown that one can choose a scaling so that \(r=1\) and \(K=1\). Then the model becomes

\[ \frac{dN}{dt} = N\left(1-N\right) \]

which is a scalar autonomous differential equation. We will use Julia to obtain a numerical solution to an initial value problem for this equation and then plot the solution.

Numerical Solution for IVP

We will use the DifferentialEquations.jl package to obtain a numerical solution to the initial value problem

\[ \frac{dN}{dt} = N\left(1-N\right), \quad N(0)=N_0 \]

for \(N_0=0.1\). We first define the right-hand side of the differential equation as a function of \(t\) and \(N\).

# define the right-hand side of the differential equation
function logistic_growth(du,u,p,t)
    du[1] = u[1]*(1-u[1])
end

logistic_growth (generic function with 1 method)

Next, we define the initial condition and the time interval for the solution.

# define the initial condition
u0 = [0.1]

# define the time interval
tspan = (0.0, 10.0)

(0.0, 10.0)

Finally, we solve the initial value problem using the solve function from DifferentialEquations.jl.

# solve the initial value problem
prob = ODEProblem(logistic_growth,u0,tspan)
sol = solve(prob);

Show the code

# plot the solution
plot(sol,linewidth=2,xlabel=L"t",ylabel=L"N(t)",legend=false)

Figure 1: A numerical solution to the logistic growth model.

It is also useful to plot the slope field for the differential equation in the \(t,y\) - plane. We can do this as follows:

Show the code

# Define the differential equation dN/dt = N(1-N)
function dNdt(t, N)
    N * (1 - N)
end

# Create a grid of t and y values
t_values = range(0, stop=10, length=20)
N_values = range(-0.5, stop=1.5, length=20)

meshgrid(x, y) = (repeat(x, outer=length(y)), repeat(y, inner=length(x)));
# Create a meshgrid
time, N = meshgrid(t_values, N_values)

scale = 0.2;
# Calculate the slope at each point in the grid
slopes = scale.*dNdt.(time, N)


# Create a quiver plot for the slope field
p_sf = quiver(time, N, quiver=(scale.*ones(size(slopes)), slopes), color="#808080", lw=1.5, legend=:none)

p_sf = hline!(p_sf, [0.0], color="red", lw=1.5,linestyle=:dash)
p_sf = hline!(p_sf, [1.0], color="blue", lw=1.5,linestyle=:dash)

# Customize the plot
title!(L"Slope Field for $\frac{dN}{dt} = N(1-N)$")
xlabel!(L"t")
ylabel!(L"N")

Figure 2: Slope field for the logistic growth model

Let’s add solutions corresponding to \(N(0) = 0.1\) and \(N(0) = 0.5\) to the vector field plot:

Show the code

# define the initial condition
u0_05 = [0.5]

# solve the initial value problem
prob = ODEProblem(logistic_growth,u0_05,tspan)
sol_05 = solve(prob)

p_sf = plot!(p_sf,sol,linewidth=2.5,legend=false,color="purple")
p_sf = plot!(p_sf,sol_05,linewidth=2.5,legend=false,color="purple")

Figure 3: Slope field for the logistic growth model with two solutions corresponding to \(N_0 = 0.1\) and \(N_0 = 0.5\).

Finally, since the logistic growth model is an autonomous differential equation, we can plot the phase line for the differential equation. We can do this as follows:

Show the code

# Define the differential equation dN/dt = N(1-N)
function dNdt(N)
    N * (1 - N)
end

# Create a grid of t and y values
N_values = range(-0.5, stop=1.5, length=101)

# Create a meshgrid
N = N_values

# Calculate the slope at each point in the grid
slopes = dNdt.(N)

# Create a quiver plot for the slope field
p_pl = plot(N, slopes, lw=1.5, legend=:none)

p_pl = hline!(p_pl, [0.0], color="black", lw=1.5,linestyle=:dash)

# add equilibrium points
p_pl = scatter!(p_pl,[0.0],[0.0],color="red",legend=false)
p_pl = scatter!(p_pl,[1.0],[0.0],color="blue",legend=false)

# add arrows for direction of flow
N_values_0 = range(-0.5, stop=1.5, length=10)
slopes_0 = dNdt.(N_values_0)
quiver!(p_pl,N_values_0, zeros(length(N_values_0)), quiver=(scale.*slopes_0, zeros(size(slopes_0))),lw=1.5, legend=:none,color="purple")

# Customize the plot
title!(L"Phase Line for $\frac{dN}{dt} = N(1-N)$")
xlabel!(L"N")
ylabel!(L"\frac{dN}{dt}")

Figure 4: Phase line for the logistic growth model.

Example: Interacting Species Population Model

The Lotka-Volterra model is a model of predator-prey interactions. It is a system of two differential equations:

\[ \begin{aligned} \frac{dN}{dt} &= rN - \alpha N P \\ \frac{dP}{dt} &= \beta N P - \gamma P \end{aligned} \]

where \(N\) is the population of prey, \(P\) is the population of predators, \(r\) is the growth rate of the prey population, \(\alpha\) is the rate at which predators kill prey, \(\beta\) is the rate at which predators reproduce, and \(\gamma\) is the rate at which predators die.

Expand for Session Info

Show the code

VERSION

v"1.9.4"

Show the code

# import Pkg.jl package
import Pkg
# list packages installed in environment
Pkg.status()

Status `~/Documents/GitHub/biomath_des/Project.toml`
⌃ [13f3f980] CairoMakie v0.11.4
  [a93c6f00] DataFrames v1.6.1
  [0c46a032] DifferentialEquations v7.12.0
⌃ [31c24e10] Distributions v0.25.104
⌅ [5b8099bc] DomainSets v0.6.7
  [634d3b9d] DrWatson v2.13.0
⌃ [61744808] DynamicalSystems v3.2.3
  [b964fa9f] LaTeXStrings v1.3.1
  [94925ecb] MethodOfLines v0.10.4
  [961ee093] ModelingToolkit v8.75.0
  [8913a72c] NonlinearSolve v3.4.0
  [429524aa] Optim v1.7.8
⌃ [7f7a1694] Optimization v3.20.2
⌅ [1dea7af3] OrdinaryDiffEq v6.66.0
  [91a5bcdd] Plots v1.39.0
⌃ [f2b01f46] Roots v2.0.22
⌅ [0c5d862f] Symbolics v5.14.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

Reuse

CC BY-NC-SA 4.0