Show the code
# import packages
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spiDifferential Equations and Dynamical Systems in Python
Examples for Biomathematics
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 solve_ivp function from the scipy.integrate module 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
def rhs(t, N):
return N * (1 - N)
# Define initial value
N0 = [0.1]
# Define time interval
t0 = 0
tend = 10
# Define time points for output
times = np.arange(t0, tend + 0.1, 0.1)
# Solve the IVP
sol = spi.solve_ivp(fun=lambda t, N: rhs(t, N), t_span=(t0, tend), y0=N0, t_eval=times)Now, let’s plot the solution. To do so, we will use the ggplot2 package. Note that this requires us to format out solution as a data frame
Show the code
# Plot the solution
plt.plot(sol.t, sol.y[0])
plt.xlabel('Time')
plt.ylabel('N')
plt.title('Solution of the ODE')
plt.show()
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)
def dNdt(t, N):
return N * (1 - N)
# Create a grid of t and y values
t_values = np.linspace(0, 10, 20)
N_values = np.linspace(-0.5, 1.5, 20)
# Create a meshgrid
time, N = np.meshgrid(t_values, N_values)
# Calculate the slope at each point in the grid
slopes = dNdt(time, N)
# Create a quiver plot for the slope field
plt.quiver(time, N, np.ones_like(slopes), slopes, color="#808080", linewidth=1.5, angles='xy', scale_units='xy', scale=5)
# Add horizontal lines at N=0 and N=1
plt.axhline(0, color="red", linestyle='--', linewidth=1.5)
plt.axhline(1, color="blue", linestyle='--', linewidth=1.5)
# Customize the plot
plt.title("Slope Field for $\\frac{dN}{dt} = N(1-N)$")
plt.xlabel("t")
plt.ylabel("N")
plt.show()
Let’s add solutions corresponding to \(N(0) = 0.1\) and \(N(0) = 0.5\) to the vector field plot:
Show the code
# Create a grid of t and y values
t_values = np.linspace(0, 10, 20)
N_values = np.linspace(-0.5, 1.5, 20)
N0_05 = [0.5]
# Solve the IVP
sol_05 = spi.solve_ivp(fun=lambda t, N: rhs(t, N), t_span=(t0, tend), y0=N0_05, t_eval=times)
# Create a quiver plot for the slope field
plt.quiver(time, N, np.ones_like(slopes), slopes, color="#808080", linewidth=1.5, angles='xy', scale_units='xy', scale=5)
# Add horizontal lines at N=0 and N=1
plt.axhline(0, color="red", linestyle='--', linewidth=1.5)
plt.axhline(1, color="blue", linestyle='--', linewidth=1.5)
# Add solutions
plt.plot(sol.t, sol.y[0], color="black", linewidth=2)
plt.plot(sol.t, sol_05.y[0], color="black", linewidth=2)
# Customize the plot
plt.title("Slope Field for $\\frac{dN}{dt} = N(1-N)$")
plt.xlabel("t")
plt.ylabel("N")
plt.show()
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)
def dNdt(N):
return N * (1 - N)
# Create a grid of N values
N_values = np.linspace(-0.5, 1.5, 101)
# Calculate the slope at each point in the grid
slopes = dNdt(N_values)
# Create a plot for the phase line
plt.plot(N_values, slopes, linewidth=1.5)
# Add horizontal dashed line at y=0
plt.axhline(0, color="black", linestyle='--', linewidth=1.5)
# Add equilibrium points
plt.scatter([0.0, 1.0], [0.0, 0.0], color=["red", "blue"])
# Add arrows for direction of flow
N_values_0 = np.linspace(-0.5, 1.5, 10)
slopes_0 = dNdt(N_values_0)
plt.quiver(N_values_0, np.zeros_like(slopes_0), slopes_0, np.zeros_like(slopes_0), scale=5, color="purple", width=0.01)
# Customize the plot
plt.title("Phase Line for $\\frac{dN}{dt} = N(1-N)$")
plt.xlabel("N")
plt.ylabel("$\\frac{dN}{dt}$")
plt.show()
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.
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
import sys
import numpy as np
import scipy as sp
import matplotlib as mpl
print("Python version:", sys.version)
print('\n'.join(f'{m.__name__}=={m.__version__}' for m in globals().values() if getattr(m, '__version__', None)))Python version: 3.8.18 | packaged by conda-forge | (default, Oct 10 2023, 15:46:56)
[Clang 16.0.6 ]
numpy==1.24.4
scipy==1.10.1
matplotlib==3.7.2