含条件语句的ODE组求解问题：Python模拟结果不符

问题背景 #

需要模拟营养盐(Am)与浮游植物(Phy)的浓度随时间变化，核心需求是每隔20天收获95%的体系生物量，但原代码的条件判断逻辑无法触发收获操作，导致结果不符合预期。

问题根源 #

1. 自适应步长的精度问题：scipy.integrate.solve_ivp采用自适应步长算法，几乎不会精确落在t = 20,40,...这类整数时间点，因此原代码中t % har_f == 0的条件基本无法触发。
2. 错误的突变模拟方式：原代码试图用har_pc/dt将离散收获转化为连续稀释项，但这种方式依赖固定步长，与solve_ivp的自适应步长不兼容，无法正确模拟瞬间收获的状态突变。

解决策略：使用事件函数捕捉离散突变 #

针对带离散状态突变的ODE求解，正确的做法是用solve_ivp的事件机制精确捕捉收获时刻，然后直接更新状态变量，再继续后续积分。

修正后的代码 #

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

# 参数定义
dil = 0.2            # 日常稀释率（天⁻¹）
har_f = 20           # 收获周期（天）
har_pc = 0.95        # 收获比例
ext_Am = 100         # 新鲜培养基营养盐浓度
umax_Phy = 0.693     # 浮游植物最大生长率（天⁻¹）
kAm_Phy = 14         # 营养盐半饱和常数
t0 = 0           
tf = 100              

# 定义连续ODE方程组（仅处理日常稀释和生长）
def odes(t, x):
    Am, Phy = x
    # 日常稀释项
    in_Am = ext_Am * dil
    out_Am = Am * dil
    out_Phy = Phy * dil
    # 浮游植物生长速率
    u_Phy = umax_Phy * Am / (Am + kAm_Phy)
    gro_Phy = u_Phy * Phy
    # 微分方程
    dAm = in_Am - out_Am - gro_Phy
    dPhy = gro_Phy - out_Phy
    return [dAm, dPhy]

# 定义收获事件函数：精确捕捉收获时刻
def harvest_event(t, x):
    # 加入小偏移量避免浮点精度误差，确保触发事件
    return (t % har_f) - 1e-8
harvest_event.direction = 0  # 捕捉时间点（双向穿越均触发）
harvest_event.terminal = False  # 触发事件后不终止积分

# 初始条件
x0 = [99, 1]
tspan = (t0, tf)

# 首次求解，捕捉所有收获事件
sol = solve_ivp(odes, tspan, x0, events=harvest_event, dense_output=True)

# 处理收获事件：更新状态并分段积分
event_times = sol.t_events[0]
# 初始化结果数组
t = sol.t.copy()
Am = sol.y[0].copy()
Phy = sol.y[1].copy()

for event_t in event_times:
    # 找到事件时刻在结果数组中的位置
    idx = np.searchsorted(t, event_t)
    # 计算收获后的状态：移除har_pc比例的培养液，补充新鲜培养基
    new_Am = Am[idx-1] * (1 - har_pc) + ext_Am * har_pc
    new_Phy = Phy[idx-1] * (1 - har_pc)
    # 重新计算事件时刻到终点的解
    sol_post = solve_ivp(odes, (event_t, tf), [new_Am, new_Phy], dense_output=True)
    # 合并结果（移除原数组中事件时刻之后的部分，避免重复）
    t = t[:idx]
    Am = Am[:idx]
    Phy = Phy[:idx]
    t = np.concatenate([t, sol_post.t[1:]])
    Am = np.concatenate([Am, sol_post.y[0][1:]])
    Phy = np.concatenate([Phy, sol_post.y[1][1:]])

# 绘制结果
fig, ax = plt.subplots()
ax.plot(t, Am, '-r', label='Am')
ax.plot(t, Phy, '-g', label='Phy')
ax.legend(loc='best', frameon=True, edgecolor='black', fancybox=False)
ax.set_xlabel('Time [day]', size=12)
ax.set_ylabel('Concentration [$g\\:biomass\\:m^{-3}$]', size=12)
ax.tick_params(axis='both', which='both', direction='in')
ax.grid(True)
plt.show()

关键修改说明 #

1. 分离连续与离散过程：日常稀释和生物生长放在连续ODE中计算，收获作为离散状态突变单独处理，更符合实际物理过程。
2. 事件函数精确触发：通过harvest_event捕捉收获时刻，加入1e-8偏移量避免浮点精度问题，确保不会错过事件。
3. 状态突变更新：收获时直接计算新的营养盐和浮游植物浓度（剩余部分+新鲜培养基补充），然后分段积分后续时间段，保证突变后的状态正确传递。

内容的提问来源于stack exchange，提问作者geo.freitas