当然有办法解决这个问题！你遇到的u超出0-1范围的情况，在solve_ivp里可以通过几种不同的思路来处理，我给你详细拆解一下：

方法一：用事件函数实现分段强制截断 #

solve_ivp支持事件函数，可以在变量触及边界时触发积分停止，之后我们手动修正变量值，再继续积分直到结束时间。这种方法能严格保证u始终在0-1区间内。

实现代码 #

import numpy as np
from scipy.integrate import solve_ivp

# 保留你原有的模型函数
def model_case_3(t, z, Kmax, k, b, list_Kmax, sigma):
    a=1.5
    x, u= z
    m=1
    r = list_Kmax[0][0]
    dxdt = (x)*(r*(1-a*u**2)*(1-x/(Kmax*(1-0.999*u**2)))-m/(k+b*u)-0.05)
    dudt = sigma*((-2*a*(b**2)*r*(u**3)+4*a*b*k*r*(u**2)+2*a*(k**2)*r*u-b*m)/((b*u+k)**2))
    return [dxdt, dudt]

# 定义触发边界事件的函数
def u_hit_lower(t, z, *args):
    return z[1] - 0  # u=0时返回0，触发停止
u_hit_lower.terminal = True  # 触发事件时终止当前积分段
u_hit_lower.direction = -1  # 仅在u从正值降到0时触发

def u_hit_upper(t, z, *args):
    return z[1] - 1  # u=1时返回0，触发停止
u_hit_upper.terminal = True
u_hit_upper.direction = 1  # 仅在u从小于1升到1时触发

# 封装分段积分逻辑
def solve_with_bounded_u(t_span, y0, args, t_eval):
    t_start, t_end = t_span
    current_t = t_start
    current_y = y0.copy()
    
    # 存储最终结果
    ts, xs, us = [], [], []
    
    while current_t < t_end:
        # 筛选当前积分段需要评估的时间点
        current_eval_points = [t for t in t_eval if t >= current_t]
        if not current_eval_points:
            break
        next_segment_end = min(current_eval_points[-1], t_end)
        
        # 求解当前积分段
        sol_segment = solve_ivp(
            fun=model_case_3,
            t_span=[current_t, next_segment_end],
            y0=current_y,
            t_eval=current_eval_points,
            args=args,
            events=[u_hit_lower, u_hit_upper]
        )
        
        # 保存当前段结果
        ts.extend(sol_segment.t)
        xs.extend(sol_segment.y[0])
        us.extend(sol_segment.y[1])
        
        # 如果触发了边界事件，修正u并继续积分
        if sol_segment.status == 1:
            # 获取事件触发的时间点
            current_t = sol_segment.t_events[0][0] if sol_segment.t_events[0].size > 0 else sol_segment.t_events[1][0]
            # 将u修正到边界值
            current_y[1] = 0.0 if sol_segment.y_events[0].size > 0 else 1.0
            current_y[0] = sol_segment.y[0][-1]
        else:
            current_t = next_segment_end
            current_y = sol_segment.y[:, -1]
    
    # 整理成类似solve_ivp的返回格式
    from collections import namedtuple
    BoundedResult = namedtuple('BoundedResult', ['t', 'y'])
    return BoundedResult(t=np.array(ts), y=np.array([xs, us]))

# 使用封装后的函数求解
sol = solve_with_bounded_u(
    t_span=[scaled_days[i][j][0], scaled_days[i][j][-1]],
    y0=[scaled_pop[i][j][0], list_u[0][0][0]],
    args=(list_Kmax[0][0], k0, b0, list_b, sigma0),
    t_eval=scaled_days[i][j]
)

方法二：修改微分方程，添加边界惩罚项 #

这种方法更简单：直接在dudt中添加惩罚项，当u超出0-1范围时，强制将它拉回区间内。适合对约束严格性要求不是极高的场景。

实现代码 #

def model_case_3_with_bounds(t, z, Kmax, k, b, list_Kmax, sigma):
    a=1.5
    x, u= z
    m=1
    r = list_Kmax[0][0]
    
    # 计算原始的dudt
    dudt_original = sigma*((-2*a*(b**2)*r*(u**3)+4*a*b*k*r*(u**2)+2*a*(k**2)*r*u-b*m)/((b*u+k)**2))
    
    # 添加边界约束惩罚项：超出范围时拉回区间
    penalty_strength = 10  # 可根据实际情况调整强度
    if u < 0:
        dudt = dudt_original + penalty_strength*(0 - u)
    elif u > 1:
        dudt = dudt_original + penalty_strength*(1 - u)
    else:
        dudt = dudt_original
    
    dxdt = (x)*(r*(1-a*u**2)*(1-x/(Kmax*(1-0.999*u**2)))-m/(k+b*u)-0.05)
    return [dxdt, dudt]

# 直接用solve_ivp求解修改后的模型
sol = solve_ivp(
    fun=model_case_3_with_bounds,
    t_span=[scaled_days[i][j][0], scaled_days[i][j][-1]],
    y0=[scaled_pop[i][j][0], list_u[0][0][0]],
    t_eval=scaled_days[i][j],
    args=(list_Kmax[0][0], k0, b0, list_b, sigma0)
)

方法三：切换到支持内置约束的求解器 #

solve_ivp本身没有内置变量约束功能，如果你需要更专业的带约束微分方程求解，可以尝试CasADi或Pyomo这类库。以CasADi为例，它可以直接设置变量的上下界：

实现代码示例 #

from casadi import *

# 定义变量与参数
t = MX.sym('t')
x = MX.sym('x')
u = MX.sym('u')
Kmax = MX.sym('Kmax')
k = MX.sym('k')
b = MX.sym('b')
r = MX.sym('r')
sigma = MX.sym('sigma')
a = 1.5
m = 1

# 定义微分方程
dxdt = (x)*(r*(1-a*u**2)*(1-x/(Kmax*(1-0.999*u**2)))-m/(k+b*u)-0.05)
dudt = sigma*((-2*a*(b**2)*r*(u**3)+4*a*b*k*r*(u**2)+2*a*(k**2)*r*u-b*m)/((b*u+k)**2))

# 创建ODE对象
ode = {'x': vertcat(x, u), 't': t, 'p': vertcat(Kmax, k, b, r, sigma), 'ode': vertcat(dxdt, dudt)}

# 创建带边界约束的求解器：设置u的上下界为0-1
opts = {'ipopt.print_level': 0, 'print_time': 0}
solver = integrator(
    'solver', 
    'cvodes', 
    ode, 
    {'t0': scaled_days[i][j][0], 'tf': scaled_days[i][j][-1], 
     'bound_x': [vertcat(-inf, 0), vertcat(inf, 1)]},  # x无约束，u约束在0-1
    opts
)

# 传入初始值与参数
params = vertcat(list_Kmax[0][0], k0, b0, list_Kmax[0][0], sigma0)
x0 = vertcat(scaled_pop[i][j][0], list_u[0][0][0])

# 求解并插值到目标时间点
res = solver(x0=x0, p=params)
sol_t = np.array(scaled_days[i][j])
sol_x = np.interp(sol_t, [scaled_days[i][j][0], scaled_days[i][j][-1]], res['xf'][[0]].full().flatten())
sol_u = np.interp(sol_t, [scaled_days[i][j][0], scaled_days[i][j][-1]], res['xf'][[1]].full().flatten())

方法对比 #

• 方法一：严格保证变量在约束范围内，但需要手动处理分段逻辑；
• 方法二：实现最简单，适合快速验证，但约束不是绝对严格；
• 方法三：适合复杂约束场景，但需要学习新的库。

内容的提问来源于stack exchange，提问作者user606273