在这里我们求解的是23年美赛的A题中的两个基本模型为例子,用obeint函数来求解改题里面的的微分方程和偏微分方程。在这一年的美赛里面的最基本的模型是人口理论模型:
d
N
d
t
=
r
∗
N
∗
(
1
−
N
K
)
\frac{dN}{dt} = r*N*(1-\frac{N}{K}) \\
dtdN=r∗N∗(1−KN)
N
N
N表示人口数量,
t
t
t表示时间,
K
K
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()
最后画出来的函数图像如下图:
当考虑到两个物种的情况时,模型公式如下:
d
N
1
d
t
=
r
1
∗
N
1
∗
(
1
−
N
1
K
1
−
δ
12
N
2
K
2
)
d
N
2
d
t
=
r
2
∗
N
2
∗
(
1
−
N
2
K
2
−
δ
21
N
1
K
1
)
\frac{dN_1}{dt} = r_1*N_1*(1-\frac{N_1}{K_1}-\delta_{12}\frac{N_2}{K_2}) \\ \frac{dN_2}{dt} = r_2*N_2*(1-\frac{N_2}{K_2}-\delta_{21}\frac{N_1}{K_1}) \\
dtdN1=r1∗N1∗(1−K1N1−δ12K2N2)dtdN2=r2∗N2∗(1−K2N2−δ21K1N1)
N
1
,
N
2
N_1,N_2
N1,N2:表示物种的种群数量
K
1
,
K
2
K_1,K_2
K1,K2:物种的环境容纳量
r
1
,
r
2
r_1,r_2
r1,r2:物种的种群增长率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()
画出来的函数图像为: