Python求解耦合ODE：振荡异常与虚假过早稳态问题

问题描述 #

我正在复现一篇论文中的生物耦合一阶微分方程组模型，该系统具有生物特性，负解无实际意义。已知论文预期结果图，但当前求解结果完全不符合预期，存在两种异常情况：

• 默认求解器：no与np出现剧烈振荡（甚至产生负解），而nc与ns保持恒定，不符合模型预期；
• BDF求解器：0秒后立即进入非系统真实状态的“稳态”，无参考价值。

为满足ntitin总量恒定的约束，我仅使用四个微分方程中的三个进行求解，场景1的代码如下：

import numpy as nump

# 常数定义
kp = 0.07   # M^-1 s^-1
kr = 6  # s^-1
ks = (10**(-8)) - (10**(-6))  # s^-1
kdns = 0.002  # s^-1
ATP = 7.2*10**(-3)  # M
w0 = 10**6  # s^-1
kb = 1.3806452*10**(-23)  # J/K
T = 310  # K
u0 = 1*10**(-9)  # m
umax = 28*10**(-9)  # m
delta_G = 35*kb*T
ntitin = 4.7 * 10**(-6)  # 单位：摩尔(M)  

# 微分方程组
def equations(t, y):
    nc, no, np = y
    
    # 通过约束条件计算ns
    ns = ntitin - (no + nc + np)
    
    # 根据当前力值计算Eplus和Eminus
    f = force_input(t)  # 调用已定义的力输入函数
    Eplus = delta_G - f * u0
    Eminus = f * (umax - u0)
    
    # 根据当前力值计算速率常数
    kplus = w0 * nump.exp((-Eplus) / (kb * T))
    kminus = w0 * nump.exp((-Eminus) / (kb * T))
    
    # 微分方程组（因ns由约束条件给出，故不包含dns_dt）
    dnc_dt = -kplus * nc + kminus * no
    dno_dt = kplus * nc - kminus * no - kp * no * ATP + kr * np + kdns * ns
    dnp_dt = kp * no * ATP - kr * np - ks * np
    #dns_dt = ks * np - kdns * ns

    return [dnc_dt, dno_dt, dnp_dt]

# 初始条件
no0 = 0.0006*ntitin
np0 = 0.0025*ntitin
ns0 = 10e-12
nc0 = ntitin-no0-np0-ns0 # 假设ns初始值为10e-12，且ntitin总量保持恒定
y0 = [nc0, no0, np0]
# 若NS初始值设为10E-6，no会在时间范围内急剧增加，而NC、NP和NS保持稳定

# 求解时间范围
t_span = [-1, 3 * 260]  # 3个完整输入周期
t_values = nump.linspace(t_span[0], t_span[1], 1000)

# 求解方程组
solution = solve_ivp(
    equations, t_span, y0, t_eval = t_values)

# 根据约束条件计算每个时间点的ns值
ns_values = ntitin - (solution.y[0] + solution.y[1] + solution.y[2])
sum_values = (ns_values + solution.y[1] + solution.y[1])/ntitin

# 绘制结果
#plt.plot(solution.t, solution.y[0]/ntitin, label='nc')
plt.plot(solution.t, solution.y[1]/ntitin, label='no')
plt.plot(solution.t, solution.y[2]/ntitin, label='np')
plt.plot(solution.t, sum_values, label='ns+np+no', linestyle='--')  # ns用虚线表示
plt.legend()
plt.xlabel('时间 (s)')
plt.ylabel('开放/磷酸化肌联蛋白的比例')
plt.title('带约束条件的系统随时间动态变化')
plt.show()

已验证力输入函数force_input的正确性，其代码如下：

def force_input(t):
    """
    生成随时间变化的力输入：先重复10次开关循环，然后进入休息状态。
    
    :param t: 时间（秒）
    :return: 对应时间点的力（单位：N）
    """
    high_force = 20 * 10**(-12)  # "开"状态的力值
    low_force = 4.5 * 10**(-12) # "关"状态和休息状态的力值
    
    # 一次"开关"周期为20秒，整个序列周期为260秒
    on_off_period = 20
    full_set_period = 260
    
    # 计算当前在序列中的位置
    position_in_set = t % full_set_period
    
    if t<0:
        return low_force
    # 判断是否处于休息周期
    if position_in_set > 10 * on_off_period:
        return low_force
    else:
        # 若不在休息周期，判断当前处于"开关"周期的"开"还是"关"阶段
        position_in_on_off = position_in_set % on_off_period
        return high_force if position_in_on_off < 10 else low_force

前10个时间点的浓度计算结果显示，no与np相互补偿振荡，甚至进入负值区间：

nc:  [0.99689787 0.99689787 0.99689785 0.99689779 0.99689774 0.99689769
 0.99689764 0.99689758 0.99689753 0.99689748]
no:  [ 0.0006      0.19058381 -0.33508352 -0.02302532 -0.01392478 -0.14079521
  0.01207641 -0.07641563  0.0918557   0.04086486]
np:  [ 0.0025     -0.18748378  0.3381835   0.0261254   0.01702492  0.14389538
 -0.00897615  0.07951593 -0.08875532 -0.03776443]
ns:  [1.00000000e-11 9.83700948e-12 1.02291693e-11 9.97153525e-12
 9.94892783e-12 1.00318165e-11 9.89773534e-12 9.95090262e-12
 9.80492439e-12 9.82905055e-12]

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