2025年美赛E题超详细解题思路+开端示例模子代码参考+代码开源 - ToB企服应用市场:ToB评测及商务社交产业平台

% 示例主脚本：演示基于简化的农业生态系统ODE模型进行仿真
clc; clear; close all;
%% 1. 模拟时间设置
tspan = [0 100]; % 模拟 0 ~ 100（单位：天）可根据需求自行调整
%% 2. 初始条件 (y0)
% y = [P; I; B; A]
% P: 植物(农作物)，I: 害虫，B: 蝙蝠，A: 鸟类
P0 = 40;
I0 = 5;
B0 = 2;
A0 = 1;
y0 = [P0; I0; B0; A0];
%% 3. 模型参数设置
% 下面给出一组示例参数，仅供演示，实际需根据实验数据来修正
params.r_p = 0.5; % 植物的固有增长率
params.K_p = 100; % 植物环境承载量
params.alpha_PI = 0.01; % 害虫对植物的取食系数
params.mu_p = 0.01; % 除草剂或其他负面影响的基数
params.r_i = 1.0; % 害虫固有增长率
params.K_i = 50; % 害虫环境承载量
params.beta_IP = 0.01; % 植物对害虫繁殖的促进系数
params.gamma_IB = 0.01; % 蝙蝠捕食害虫的效率
params.gamma_IA = 0.01; % 鸟类捕食害虫的效率
params.mu_i = 0.02; % 杀虫剂或其他负面影响的基数
params.r_b = 0.3; % 蝙蝠固有增长率
params.K_b = 20; % 蝙蝠环境承载量
params.delta_BI = 0.001; % 蝙蝠捕食害虫获益对自身繁殖的促进
params.mu_b = 0.1; % 蝙蝠的自然死亡率
params.r_a = 0.3; % 鸟类固有增长率
params.K_a = 15; % 鸟类环境承载量
params.delta_AI = 0.001; % 鸟类捕食害虫获益对自身繁殖的促进
params.mu_a = 0.1; % 鸟类的自然死亡率
% f_herbicide(C_h) = theta1 * C_h
% f_pesticide(C_p) = theta2 * C_p
% 这里给出常数系数作为示例
params.theta1 = 0.01; % 除草剂影响系数
params.theta2 = 0.02; % 杀虫剂影响系数
%% 4. 人为决策场景设置
% 场景1: 高度依赖化学投入（C_h=1, C_p=1）
% 场景2: 逐步减少化学投入（C_h=0.3, C_p=0.3）
% 场景3: 完全停止化学投入（C_h=0, C_p=0）
scenarioNames = {'HighChem','LowChem','NoChem'};
Ch_values = [1, 0.3, 0]; % 除草剂使用强度
Cp_values = [1, 0.3, 0]; % 杀虫剂使用强度
%% 5. 循环跑不同场景并绘图
figure('Position',[100 100 1200 400]);
for iScen = 1:length(scenarioNames)
params.C_h = Ch_values(iScen);
params.C_p = Cp_values(iScen);
% 调用 ode45 求解
[tSol, ySol] = ode45(@(t,y) ecosystemODE(t,y,params), tspan, y0);
% ySol(:,1) = P, ySol(:,2) = I, ySol(:,3) = B, ySol(:,4) = A
subplot(1,length(scenarioNames),iScen);
plot(tSol, ySol(:,1), 'g-', 'LineWidth',1.5); hold on;
plot(tSol, ySol(:,2), 'r-', 'LineWidth',1.5);
plot(tSol, ySol(:,3), 'b-', 'LineWidth',1.5);
plot(tSol, ySol(:,4), 'm-', 'LineWidth',1.5);
xlabel('Time');
ylabel('Population (arbitrary unit)');
legend({'P(植物)','I(害虫)','B(蝙蝠)','A(鸟类)'}, 'Location','best');
title(['Scenario: ' scenarioNames{iScen} ...
' | C_h=' num2str(params.C_h) ' | C_p=' num2str(params.C_p)]);
grid on;
end
function dydt = ecosystemODE(t, y, params)
%ecosystemODE 定义四元常微分方程组
% 输入：
% t: 时间（标量），但在该模型中并未显式使用
% y: 状态变量向量 [P; I; B; A]
% params: 参数结构体，包含以下字段示例：
% r_p, K_p, alpha_PI, mu_p,
% r_i, K_i, beta_IP, gamma_IB, gamma_IA, mu_i,
% r_b, K_b, delta_BI, mu_b,
% r_a, K_a, delta_AI, mu_a,
% C_h, C_p, theta1, theta2
%
% 输出：
% dydt: 长度为4的列向量，对应 [dP/dt; dI/dt; dB/dt; dA/dt]
% 提取当前状态变量
P = y(1); % 植物
I = y(2); % 害虫
B = y(3); % 蝙蝠
A = y(4); % 鸟类
% 提取参数
r_p = params.r_p; K_p = params.K_p; alpha_PI = params.alpha_PI; mu_p = params.mu_p;
r_i = params.r_i; K_i = params.K_i; beta_IP = params.beta_IP; mu_i = params.mu_i;
gamma_IB = params.gamma_IB; gamma_IA = params.gamma_IA;
r_b = params.r_b; K_b = params.K_b; delta_BI = params.delta_BI; mu_b = params.mu_b;
r_a = params.r_a; K_a = params.K_a; delta_AI = params.delta_AI; mu_a = params.mu_a;
C_h = params.C_h; C_p = params.C_p; theta1 = params.theta1; theta2 = params.theta2;
% 定义除草剂和杀虫剂影响
% f_herbicide(C_h) = theta1 * C_h
% f_pesticide(C_p) = theta2 * C_p
f_herb = theta1 * C_h;
f_pest = theta2 * C_p;
%% 计算微分方程
% dP/dt: 植物
% 植物的逻辑斯蒂增长: r_p * P (1 - P/K_p)
% 被害虫取食: - alpha_PI * P * I
% 除草剂的负面影响(简化表示): - mu_p * f_herb
dPdt = r_p * P * (1 - P/K_p) ...
- alpha_PI * P * I ...
- mu_p * f_herb;
% dI/dt: 害虫
% 害虫逻辑斯蒂增长: r_i * I (1 - I/K_i)
% 从植物获取食物繁殖: + beta_IP * I * P (示意)
% 被蝙蝠、鸟类捕食: - gamma_IB * I * B - gamma_IA * I * A
% 杀虫剂的负面影响(简化表示): - mu_i * f_pest
dIdt = r_i * I * (1 - I/K_i) ...
+ beta_IP * I * P ...
- gamma_IB * I * B ...
- gamma_IA * I * A ...
- mu_i * f_pest;
% dB/dt: 蝙蝠
% 蝙蝠逻辑斯蒂增长: r_b * B (1 - B/K_b)
% 捕食害虫获益: + delta_BI * B * I
% 自然死亡率: - mu_b
dBdt = r_b * B * (1 - B/K_b) ...
+ delta_BI * B * I ...
- mu_b;
% dA/dt: 鸟类
% 鸟类逻辑斯蒂增长: r_a * A (1 - A/K_a)
% 捕食害虫获益: + delta_AI * A * I
% 自然死亡率: - mu_a
dAdt = r_a * A * (1 - A/K_a) ...
+ delta_AI * A * I ...
- mu_a;
% 输出
dydt = [dPdt; dIdt; dBdt; dAdt];
end

复制代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
plt.rcParams['font.sans-serif'] = ['SimHei'] # 使用SimHei字体
plt.rcParams['axes.unicode_minus'] = False # 正确显示负号
# ========== 1. 定义微分方程函数 ==========
def ecosystem_ode(t, y, params):
"""
ecosystem_ode: 定义四元常微分方程组
y = [P, I, B, A], 分别代表植物(作物)、害虫、蝙蝠、鸟类
params: 字典，包含各种参数
返回: [dPdt, dIdt, dBdt, dAdt]
"""
# 提取当前状态变量
P, I, B, A = y
# 提取参数
r_p = params['r_p']
K_p = params['K_p']
alpha_PI = params['alpha_PI']
mu_p = params['mu_p']
r_i = params['r_i']
K_i = params['K_i']
beta_IP = params['beta_IP']
mu_i = params['mu_i']
gamma_IB = params['gamma_IB']
gamma_IA = params['gamma_IA']
r_b = params['r_b']
K_b = params['K_b']
delta_BI = params['delta_BI']
mu_b = params['mu_b']
r_a = params['r_a']
K_a = params['K_a']
delta_AI = params['delta_AI']
mu_a = params['mu_a']
# 化学投入强度及其影响系数
C_h = params['C_h'] # 除草剂使用强度
C_p = params['C_p'] # 杀虫剂使用强度
theta1 = params['theta1'] # 除草剂影响系数
theta2 = params['theta2'] # 杀虫剂影响系数
# 定义化学品的影响函数
# f_herbicide(C_h) = theta1 * C_h
# f_pesticide(C_p) = theta2 * C_p
f_herb = theta1 * C_h
f_pest = theta2 * C_p
# --------- dP/dt: 植物 ---------
# 1) 植物逻辑斯蒂增长: r_p * P(1 - P/K_p)
# 2) 被害虫取食: - alpha_PI * P * I
# 3) 除草剂负面影响(简化): - mu_p * f_herb
dPdt = r_p * P * (1 - P/K_p) \
- alpha_PI * P * I \
- mu_p * f_herb
# --------- dI/dt: 害虫 ---------
# 1) 害虫逻辑斯蒂增长: r_i * I(1 - I/K_i)
# 2) 从植物获取食物的繁殖促进: + beta_IP * I * P
# 3) 被蝙蝠、鸟类捕食: - gamma_IB * I * B, - gamma_IA * I * A
# 4) 杀虫剂负面影响(简化): - mu_i * f_pest
dIdt = r_i * I * (1 - I/K_i) \
+ beta_IP * I * P \
- gamma_IB * I * B \
- gamma_IA * I * A \
- mu_i * f_pest
# --------- dB/dt: 蝙蝠 ---------
# 1) 蝙蝠逻辑斯蒂增长: r_b * B(1 - B/K_b)
# 2) 捕食害虫获益(简化): + delta_BI * B * I
# 3) 自然死亡率: - mu_b
dBdt = r_b * B * (1 - B/K_b) \
+ delta_BI * B * I \
- mu_b
# --------- dA/dt: 鸟类 ---------
# 1) 鸟类逻辑斯蒂增长: r_a * A(1 - A/K_a)
# 2) 捕食害虫获益(简化): + delta_AI * A * I
# 3) 自然死亡率: - mu_a
dAdt = r_a * A * (1 - A/K_a) \
+ delta_AI * A * I \
- mu_a
return [dPdt, dIdt, dBdt, dAdt]
def main():
# ========== 2. 设置仿真时间 ==========
t_start = 0
t_end = 100
# 为了在绘图时有较好的分辨率，定义一个时间采样点
t_eval = np.linspace(t_start, t_end, 300)
# ========== 3. 初始条件 ==========
# P(植物), I(害虫), B(蝙蝠), A(鸟类)
P0 = 40
I0 = 5
B0 = 2
A0 = 1
y0 = [P0, I0, B0, A0]
# ========== 4. 设置模型基本参数 ==========
params = {
# 植物相关
'r_p': 0.5, # 植物固有增长率
'K_p': 100, # 植物环境承载量
'alpha_PI': 0.01, # 害虫对植物取食系数
'mu_p': 0.01, # 植物受除草剂的影响基数
# 害虫相关
'r_i': 1.0, # 害虫固有增长率
'K_i': 50, # 害虫环境承载量
'beta_IP': 0.01, # 植物对害虫繁殖的促进系数
'mu_i': 0.02, # 害虫受杀虫剂的影响基数
'gamma_IB': 0.01, # 蝙蝠捕食害虫效率
'gamma_IA': 0.01, # 鸟类捕食害虫效率
# 蝙蝠相关
'r_b': 0.3, # 蝙蝠固有增长率
'K_b': 20, # 蝙蝠环境承载量
'delta_BI': 0.001, # 蝙蝠捕食害虫对自身繁殖的促进
'mu_b': 0.1, # 蝙蝠自然死亡率
# 鸟类相关
'r_a': 0.3, # 鸟类固有增长率
'K_a': 15, # 鸟类环境承载量
'delta_AI': 0.001, # 鸟类捕食害虫对自身繁殖的促进
'mu_a': 0.1, # 鸟类自然死亡率
# 化学影响系数
'theta1': 0.01, # 除草剂影响系数
'theta2': 0.02, # 杀虫剂影响系数
# 初始先给个默认值, 后面会覆盖
'C_h': 1.0, # 除草剂使用强度
'C_p': 1.0 # 杀虫剂使用强度
}
# ========== 5. 定义场景 ==========
scenario_names = ['HighChem', 'LowChem', 'NoChem']
Ch_values = [1.0, 0.3, 0.0] # 除草剂强度
Cp_values = [1.0, 0.3, 0.0] # 杀虫剂强度
# ========== 6. 循环仿真并绘图 ==========
fig, axes = plt.subplots(1, 3, figsize=(12,4), sharey=True)
for i, scen_name in enumerate(scenario_names):
params['C_h'] = Ch_values[i]
params['C_p'] = Cp_values[i]
# 调用 solve_ivp 进行数值求解
sol = solve_ivp(
fun=lambda t, y: ecosystem_ode(t, y, params),
t_span=(t_start, t_end),
y0=y0,
t_eval=t_eval
)
t_sol = sol.t
y_sol = sol.y # shape=(4, len(t_eval))
# y_sol[0,:] = P, y_sol[1,:] = I, y_sol[2,:] = B, y_sol[3,:] = A
ax = axes[i]
ax.plot(t_sol, y_sol[0, :], 'g-', label='P(植物)', linewidth=1.5)
ax.plot(t_sol, y_sol[1, :], 'r-', label='I(害虫)', linewidth=1.5)
ax.plot(t_sol, y_sol[2, :], 'b-', label='B(蝙蝠)', linewidth=1.5)
ax.plot(t_sol, y_sol[3, :], 'm-', label='A(鸟类)', linewidth=1.5)
ax.set_title(f'Scenario: {scen_name}\nC_h={params["C_h"]}, C_p={params["C_p"]}')
ax.set_xlabel('Time')
if i == 0:
ax.set_ylabel('Population')
ax.legend(loc='best')
ax.grid(True)
plt.tight_layout()
plt.show()
if __name__ == '__main__':
main()

复制代码