如何用solve_ivp+LSODA求解复值ODE或为odeintw添加事件处理？

一、支持复值ODE + LSODA + 事件处理的替代框架 #

可以试试pyodesys库，它原生支持复值微分方程组，能够调用LSODA求解器（自动处理刚性/非刚性切换），并且内置了事件处理功能，API设计和solve_ivp较为接近，能满足你的需求。

另外，Scipy的旧版积分接口scipy.integrate.ode搭配lsoda方法也支持复值ODE，你可以通过set_solout回调函数手动实现事件检测：在每个积分步后触发回调，检查事件函数的状态，一旦满足触发条件就终止积分，再从事件点继续后续计算。

二、为odeintw添加事件处理的手动实现方式 #

由于odeintw基于Scipy的odeint修改，而odeint本身没有内置事件处理，你可以通过分步积分+事件检测的方式手动实现：

• 步骤1：分步积分并获取中间状态
  调用odeintw时，设置t_eval参数为密集的时间点，或者使用full_output=True获取积分过程中的所有中间时间和状态数据。
• 步骤2：检测事件触发
  遍历中间状态，计算事件函数的值，检查相邻时间点间事件函数是否发生符号变化（或满足其他触发条件）。
• 步骤3：精确定位事件时间
  当检测到事件可能发生时，用插值方法（比如scipy.interpolate.interp1d）在相邻时间点间精确定位事件的精确时间，然后获取该时间点的状态。
• 步骤4：继续积分
  从事件时间点开始，继续调用odeintw完成剩余区间的积分，同时根据事件需求修改状态（如重置、切换参数等）。

示例代码片段：

import numpy as np
from odeintw import odeintw
from scipy.interpolate import interp1d

def ode_system(y, t, params):
    # 你的复值ODE定义
    return ...

def event_func(y, t):
    # 事件函数，返回0时触发事件
    return np.real(y[0])  # 示例：检测第一个状态实部为0

def integrate_with_events(y0, t_span, params):
    t_start, t_end = t_span
    current_t = t_start
    current_y = y0
    results_t = [current_t]
    results_y = [current_y]
    
    while current_t < t_end:
        # 选择下一个积分区间（可根据需求调整步长）
        next_t = min(current_t + 0.1, t_end)
        t_eval = np.linspace(current_t, next_t, 100)
        y_vals = odeintw(ode_system, current_y, t_eval, args=(params,))
        
        # 计算事件函数值
        event_vals = np.array([event_func(y, t) for y, t in zip(y_vals, t_eval)])
        
        # 检查事件触发（符号变化）
        sign_changes = np.where(np.diff(np.sign(event_vals)) != 0)[0]
        if len(sign_changes) > 0:
            # 取第一个触发的事件
            idx = sign_changes[0]
            t_left, t_right = t_eval[idx], t_eval[idx+1]
            e_left, e_right = event_vals[idx], event_vals[idx+1]
            
            # 线性插值找事件时间
            t_event = t_left - e_left * (t_right - t_left)/(e_right - e_left)
            # 插值获取事件点状态
            interp_y = interp1d(t_eval, y_vals, axis=0, kind='linear')
            y_event = interp_y(t_event)
            
            # 记录事件点
            results_t.append(t_event)
            results_y.append(y_event)
            
            # 处理事件（比如重置状态，这里示例直接继续）
            current_t = t_event
            current_y = y_event
        else:
            # 无事件，记录区间终点
            results_t.append(next_t)
            results_y.append(y_vals[-1])
            current_t = next_t
            current_y = y_vals[-1]
    
    return np.array(results_t), np.array(results_y)

三、关于solve_ivp的优化 #

你之前用solve_ivp+BDF耗时过长，也可以尝试将复值ODE拆分为实部和虚部的实值方程组，然后用solve_ivp的LSODA方法（注意Scipy 1.10+版本的solve_ivp已支持LSODA），LSODA会自动切换刚性/非刚性求解模式，可能比BDF更高效。拆分的方式很简单：将每个复值状态变量拆成实部和虚部两个实变量，重新定义ODE方程组即可。

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