Preface: This article is an assignment from the author’s Mathematical Modeling course. It explores the properties of the solutions to a system of ordinary differential equations (ODEs) modeling the competition between two biological species, demonstrating the basic ideas of the qualitative theory of differential equations. Of course, the most important purpose of this article is to showcase the perfect combination of LaTeX and Python. (All figures in this article were generated using Python’s Matplotlib module, while the document was edited using LaTeX.)
Problem Statement
We study the competitive relationship between two species living in the same natural environment. It is assumed that when each species lives alone in this environment, its population evolution follows the Logistic law. Furthermore, it is assumed that when they compete with each other, they slow down each other’s population growth, and the decrease in growth rate is proportional to the product of their populations. The ODE model established according to these assumptions is:
\label{eq:jingzhengfangcheng} \left\{\begin{aligned} \frac{dx_1}{dt}&=r_1 x_1\left(1-\frac{x_1}{N_1}\right)-a_1 x_1 x_2 \\ \frac{dx_2}{dt}&=r_2 x_2\left(1-\frac{x_2}{N_2}\right)-a_2 x_1 x_2 \end{aligned}\right.
This article analyzes the properties of this equation from both quantitative and qualitative perspectives.
Model Analysis
In the analysis process, the model parameters we adopt are:
\label{eq:fangchengcanshu} \left\{\begin{aligned} &r_1 = 0.1, \quad \frac{1}{N_1}=0.002, \quad a_1=0.0001\\ &r_2=0.3, \quad \frac{1}{N_2}=0.003, \quad a_2=0.0002\\ &x_1(0)=100, \quad x_2(0)=150 \end{aligned}\right.
These parameters indicate that x_1 is a species with relatively slow growth, but its intraspecific competition is small, and it is not easily affected by the x_2 species. Therefore, in the long run, x_1 will be dominant in this natural environment. In contrast, the x_2 species has a faster growth rate but faces greater intraspecific competition and is more susceptible to the influence of x_1. Consequently, in long-term competition, x_2 will be in a slightly weaker position. The following numerical solution results support this conjecture.
Numerical Solution
Combining [eq:jingzhengfangcheng] and [eq:fangchengcanshu], we performed a numerical solution of the model using Python. The results are shown in the figure below:
The numerical solution shows that the x_2 species initially grows rapidly, starts to decline slowly after reaching a maximum value, and finally stabilizes at 250. Meanwhile, x_1 maintains a relatively slow growth rate until it reaches the equilibrium point of 375. Next, we analyze the properties of the equilibrium points of [eq:jingzhengfangcheng].
Equilibrium Points
The determining equations for the equilibrium points of equation [eq:jingzhengfangcheng] are:
\label{eq:pinghengdianfangcheng} \left\{\begin{aligned} r_1 x_1\left(1-\frac{x_1}{N_1}\right)-a_1 x_1 x_2=0 \\ r_2 x_2\left(1-\frac{x_2}{N_2}\right)-a_2 x_1 x_2=0 \end{aligned}\right.
This system has 4 solutions:
\label{eq:pinghengdian} \left\{\begin{aligned}&x_1 =0\\&x_2 =0\end{aligned}\right., \quad \left\{\begin{aligned}&x_1 =0\\&x_2 =N_2\end{aligned}\right.,\quad \left\{\begin{aligned}&x_1 = N_1\\&x_2 =0\end{aligned}\right.,\quad \left\{\begin{aligned}&x_1 = \frac{N_1 r_2 (a_1 N_2 -r_1 )}{a_1 a_2 N_1 N_2-r_1 r_2}\\&x_2 =\frac{N_2 r_1 (a_2 N_1- r_2)}{a_1 a_2 N_1 N_2-r_1 r_2} \end{aligned}\right.
The first three equilibrium points are trivial, while the last one is non-trivial. It can be observed that when the parameters are as in [eq:fangchengcanshu], equation [eq:jingzhengfangcheng] eventually converges to the fourth equilibrium point, indicating that the fourth equilibrium point has strong stability. Below, we analyze the stability of the equilibrium points and the overall properties of the equation through the direction field and phase portrait.
Direction Field
The figure below shows the direction field for [eq:jingzhengfangcheng] and [eq:fangchengcanshu], where the equilibrium points are marked with green dots. The intensity of the arrow color represents the speed of change at that point.
From the direction field, it is relatively easy to judge the stability of the four equilibrium points:
The first equilibrium point, (0,0) at the bottom left, representing the case where the populations of both species are zero, is an unstable equilibrium point. All arrows point away from this point, meaning that a tiny perturbation near this point will lead to a final state far from it. In other words, as long as a small amount of the species exists here, they will eventually multiply and reach another equilibrium.
The second equilibrium point, (0, N_2) at the top left, representing the case where only species x_2 exists, is an unstable equilibrium point (saddle point). Only arrows on the vertical line passing through this point point toward it; all other arrows point away. This indicates that while this point is more stable than (0,0), it is still unstable; the presence of even a small amount of x_1 will lead away from this equilibrium due to competition.
The third equilibrium point, (N_1, 0) at the bottom right, representing the case where only species x_1 exists, is an unstable equilibrium point (saddle point). Only arrows on the horizontal line passing through this point point toward it; all other arrows point away. Similar to the second point, it is unstable; the presence of a small amount of x_2 will lead away from the current equilibrium.
The fourth equilibrium point, \left(\frac{N_1 r_2 (a_1 N_2 -r_1 )}{a_1 a_2 N_1 N_2-r_1 r_2},\frac{N_2 r_1 (a_2 N_1- r_2)}{a_1 a_2 N_1 N_2-r_1 r_2}\right) at the top right, where both species coexist and reach equilibrium, is a stable equilibrium point. All arrows in the vicinity point toward this point, indicating that even if the system deviates from this equilibrium, it will eventually return to it.
Phase Trajectories
From the family of phase trajectories shown below, the flow near the equilibrium points can be seen more clearly. It is evident that almost all streamlines in the first quadrant (corresponding to physically meaningful solutions) converge to the fourth equilibrium point, demonstrating its very strong stability. This equilibrium point is the final stable state for any initial condition.
It can be observed that there are no closed phase trajectories in the plot; therefore, the equation does not have periodic oscillating solutions.
Code Listing
The following code requires a Python 3 environment with Numpy, Scipy, and Matplotlib.
Numerical Solution
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
# Time array
t = np.arange(0, 100, 0.1)
def deriv(w, t, a, b, c, d, e, f):
x, y = w
return np.array([a*(1-b*x)*x - c*y*x, d*(1-e*y)*y - f*x*y])
# Parameters
p = [0.1, 0.002, 0.0001, 0.3, 0.003, 0.0002, 100, 150]
a, b, c, d, e, f, x0, y0 = p
yinit = np.array([x0, y0]) # Initial values
# Solving ODE
yyy = odeint(deriv, yinit, t, args=(a, b, c, d, e, f))
plt.figure(figsize=(7, 5))
plt.plot(t, yyy[:, 0], "b-", label="$x_1$ curve")
plt.plot(t, yyy[:, 1], "r-", label="$x_2$ curve")
plt.plot([0, 100], [250, 250], "g--")
plt.plot([0, 100], [375, 375], "g--")
plt.xlabel('Time t')
plt.ylabel('Population')
plt.title('Evolution of Two Competing Species')
plt.legend(loc=4)
plt.show()Direction Field
import numpy as np
import matplotlib.pyplot as plt
p = [0.1, 0.002, 0.0001, 0.3, 0.003, 0.0002]
x0 = np.array([0, 0, 1/0.002, 375])
y0 = np.array([0, 1/0.003, 0, 250])
a, b, c, d, e, f = p
x, y = np.mgrid[-100:601:25, -100:501:25]
s = a*(1-b*x)*x - c*y*x
t = d*(1-e*y)*y - f*x*y
r = np.sqrt(s**2 + t**2)
plt.figure(figsize=(8.4, 6))
plt.plot([-100, 600], [1/0.003, 1/0.003], 'r--', alpha=0.5)
plt.plot([-100, 600], [250, 250], 'r--', alpha=0.5)
plt.plot([500, 500], [-100, 500], 'r--', alpha=0.5)
plt.plot([375, 375], [-100, 500], 'r--', alpha=0.5)
plt.scatter(x0, y0, s=75, color='green', label='Equilibrium Point')
plt.quiver(x, y, s/r, t/r, r)
plt.colorbar()
plt.xlim(-100, 601)
plt.ylim(-100, 501)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Direction Field of the Biological Competition Model')
plt.legend()
plt.show()Phase Portrait
import numpy as np
import matplotlib.pyplot as plt
x0 = np.array([0, 0, 1/0.002, 375])
y0 = np.array([0, 1/0.003, 0, 250])
p = [0.1, 0.002, 0.0001, 0.3, 0.003, 0.0002]
a, b, c, d, e, f = p
y, x = np.mgrid[-100:501:5, -100:601:5]
s = a*(1-b*x)*x - c*y*x
t = d*(1-e*y)*y - f*x*y
r = np.sqrt(s**2 + t**2)
plt.figure(figsize=(8.4, 6))
plt.plot([-100, 600], [1/0.003, 1/0.003], 'b--', alpha=0.5)
plt.plot([-100, 600], [250, 250], 'b--', alpha=0.5)
plt.plot([500, 500], [-100, 500], 'b--', alpha=0.5)
plt.plot([375, 375], [-100, 500], 'b--', alpha=0.5)
plt.plot([0, 0], [-100, 500], 'b--', alpha=0.5)
plt.plot([-100, 600], [0, 0], 'b--', alpha=0.5)
plt.scatter(x0, y0, s=75, color='green', label='Equilibrium Point')
plt.streamplot(x, y, s, t, color=r, density=1.5, linewidth=1, cmap=plt.cm.autumn)
plt.colorbar()
plt.xlim(-100, 601)
plt.ylim(-100, 501)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Phase Trajectories of the Biological Competition Model')
plt.legend()
plt.show()Original Address: https://kexue.fm/archives/3120
For more details on reprinting, please refer to: Scientific Space FAQ