如何用obeint库求解23年美赛A题基础模型

在这里我们求解的是23年美赛的A题中的两个基本模型为例子，用obeint函数来求解改题里面的的微分方程和偏微分方程。在这一年的美赛里面的最基本的模型是人口理论模型：
$$\frac{dN}{dt}=r\ast N\ast (1-\frac{N}{K})$$
$N$表示人口数量，$t$表示时间，$K$表示物种的环境容纳量。我们尝试用python求解这个微分方程，并且画出函数图像。代码如下：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义微分方程
def logistic_growth(t, N, r, K):
    return r * N * (1 - N / K)

# 参数
r = 0.1  # 固有增长率
K = 100  # 环境承载量
N0 = 10  # 初始种群大小

# 时间跨度
t_span = (0, 100)
t_eval = np.linspace(*t_span, 1000)

# 解微分方程
sol = solve_ivp(logistic_growth, t_span, [N0], args=(r, K), t_eval=t_eval)

# 绘图
plt.plot(sol.t, sol.y[0], label=f'Logistic Growth\nr={r}, K={K}, N0={N0}')
plt.xlabel('Time')
plt.ylabel('Population Size')
plt.legend()
plt.title('Logistic Growth Model')
plt.show()

最后画出来的函数图像如下图：

​当考虑到两个物种的情况时，模型公式如下：
$$\begin{gathered} \frac{d{N}_1}{dt}=r_1\ast N_1\ast (1-\frac{N_1}{K_1}-{\delta }_{12}\frac{N_2}{K_2}) \\ \frac{d{N}_2}{dt}=r_2\ast N_2\ast (1-\frac{N_2}{K_2}-{\delta }_{21}\frac{N_1}{K_1}) \end{gathered}$$
$N_1，N_2$：表示物种的种群数量
$K_1，K_2$：物种的环境容纳量
$r_1，r_2$：物种的种群增长率s
求解的代码如下：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义微分方程组
def competition_model(t, N, params):
    N1, N2 = N
    r1, K1, delta12, r2, K2, delta21 = params
    dN1dt = r1 * N1 * (1 - N1 / K1 - delta12 * N2 / K2)
    dN2dt = r2 * N2 * (1 - N2 / K2 - delta21 * N1 / K1)
    return [dN1dt, dN2dt]

# 参数
r1, K1, delta12 = 0.1, 100, 0.02  # 物种1的参数
r2, K2, delta21 = 0.1, 90, 0.02  # 物种2的参数
params = (r1, K1, delta12, r2, K2, delta21)
N0 = [10, 10]  # 初始种群大小

# 时间跨度
t_span = (0, 300)
t_eval = np.linspace(*t_span, 300)

# 解微分方程组
sol = solve_ivp(competition_model, t_span, N0, args=(params,), t_eval=t_eval, method='RK45')

# 绘图
plt.plot(sol.t, sol.y[0], label='Species 1')
plt.plot(sol.t, sol.y[1], label='Species 2')
plt.xlabel('Time')
plt.ylabel('Population Size')
plt.legend()
plt.title('Two-Species Competition Model')
plt.show()

画出来的函数图像为：