使用SciPy solve_ivp求解微分方程组时出现停滞问题求助

一、模拟停滞的原因分析与优化 #

你的模拟在t>95时停滞，核心原因大概率是系统进入刚性状态：
当t接近100时，k_t(t)中(p2*t)^p3项会变得极大，导致k_t(t)趋近于Rmin。如果Rmin数值极小，会让dCT_dt中的k_t(t)*F_t(T)*CT项趋近于0，系统方程的特征值差异变得极大，默认的RK45求解器处理这类刚性问题会不断缩小步长，导致t增长极慢甚至停滞。

针对性优化措施 #

1. 切换刚性求解器
   使用solve_ivp专门针对刚性问题的求解器（如Radau或BDF），这类求解器能高效处理刚性系统：

   solution = solve_ivp(model, t_range, y0, t_eval=t_eval, 
                        args=(p_values[i][0], ..., p_values[i][14]),
                        method='Radau')
   

2. 限制最小步长
   通过dt_min参数强制设置求解器的最小步长，避免步长无限缩小：

   solution = solve_ivp(model, t_range, y0, t_eval=t_eval, 
                        args=(...), dt_min=1e-6)  # 可根据实际调整阈值
   

3. 截断k(t)的极小值
   给k_t(t)设置最小阈值，防止数值过小触发刚性：

   def k_t(t):
       val = Rmin + (p1 / (1 + (p2 * t) ** p3))
       return max(val, 1e-10)  # 替换为适合你模型的最小阈值
   

二、超时跳过当前模拟的实现 #

如果优化后仍有个别参数组停滞，可以用多进程配合超时机制强制跳过，单线程下无法中断正在运行的solve_ivp，因此采用ProcessPoolExecutor实现：

完整代码示例 #

import concurrent.futures
import numpy as np
from scipy.integrate import solve_ivp

# 原模型参数与定义
Upsilon = 0.305
Alpha = 5.5e-7
r = 0.176
b = 5e-13

def model(t, y0, Beta, Eta, Rmin, p1, p2, p3, A, a, Xi, Epsilon, Lambda, Theta, Mu, Delta, Gamma):
    CD, CT, CM, CE, T, TOTAL = y0

    def k_t(t):
        val = Rmin + (p1 / (1 + (p2 * t) ** p3))
        return max(val, 1e-10)  # 可选：添加极小值截断

    def F_t(T):
        return T / (A + T)

    def f_t(CF, T):
        return (CF / T) / (Upsilon + ((a + CF) / T))

    dCD_dt = -(Beta + Eta) * CD
    dCT_dt = (Eta * CD) + (k_t(t) * F_t(T) * CT) - ((Xi + Epsilon + Lambda) * CT) + (Theta * T * CM) - (Alpha * T * CT)
    dCM_dt = (Epsilon * CT) - (Theta * T * CM) - (Mu * CM)
    dCE_dt = (Lambda * CT) - (Delta * CE)
    dT_dt = ((r * T) * (1 - (b * T))) - (Gamma * f_t((CD + CT), T) * T)

    TOTAL = dCD_dt + dCT_dt + dCM_dt + dCE_dt
    return [dCD_dt, dCT_dt, dCM_dt, dCE_dt, dT_dt, TOTAL]

# 封装单组参数的模拟函数
def run_single(params):
    CD = params[15]
    CT, CM, CE = 0, 0, 0
    T = 10 ** 7
    TOTAL = CD
    y0 = [CD, CT, CM, CE, T, TOTAL]
    t_range = (0, 100)
    t_eval = np.arange(0, 100, 0.1)
    sol = solve_ivp(model, t_range, y0, t_eval=t_eval,
                    args=(params[0], params[1], params[2], params[3], params[4],
                          params[5], params[6], params[7], params[8], params[9],
                          params[10], params[11], params[12], params[13], params[14]),
                    method='Radau')
    return sol.y

# 批量运行并处理超时
final_results = []
timeout_duration = 15  # 每组模拟最多允许15秒，根据实际调整

with concurrent.futures.ProcessPoolExecutor() as executor:
    # 提交所有参数组任务
    future_map = {executor.submit(run_single, p): idx for idx, p in enumerate(p_values)}
    
    for future in concurrent.futures.as_completed(future_map):
        idx = future_map[future]
        try:
            result = future.result(timeout=timeout_duration)
            final_results.append((idx, result))
        except concurrent.futures.TimeoutError:
            print(f"参数组 {idx} 模拟超时，已跳过")
            final_results.append((idx, None))

# 按原参数组顺序整理结果（可选）
final_results.sort(key=lambda x: x[0])
final_results = [res for _, res in final_results]

注意事项 #

• 用ProcessPoolExecutor而非ThreadPoolExecutor，避免Python全局解释器锁（GIL）对CPU密集型任务的限制。
• 超时时间需根据你的硬件性能和模型复杂度调整，避免误杀正常运行的参数组。

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