Scipy solve_ivp求解时变复矩阵刚性ODE的解发散归零与LSODA性能异常问题

我目前正在Python中使用solve_ivp的LSODA方法求解一个刚性ODE问题，控制方程是：

u' = Lu

其中u是复2M维向量，L是随时间变化的复2M×2M矩阵。矩阵的定义比较复杂，但我已经在外部完成了正确计算，现在只需要基于给定数据求解问题。

矩阵计算代码实现 #

def compute_matrices(psi):

    Q = M - N * q(psi)
    
    A = ((2 * np.pi) ** 2) * (
        N @ (G_22(psi) @ N + G_23(psi) @ M)
      + M @ (G_23(psi) @ N + G_33(psi) @ M)
    )
    
    B = -2 * np.pi * 1j * chi_r(psi) * (
        N @ (G_22(psi) + q(psi) * G_23(psi))
      + M @ (G_23(psi) + q(psi) * G_33(psi))
    )
    
    C = (
        (-2 * np.pi * 1j) * (
            chi_rr(psi) * (N @ G_22(psi) + M @ G_23(psi))
          + qchi_r_derivative(psi) * (N @ G_23(psi) + M @ G_33(psi))
        )
        - ((2 * np.pi) ** 2) * chi_r(psi) * (
            (N @ G_12(psi) + M @ G_31(psi)) @ Q
            )
        - np.pi * 1j * (
            Q @ J_theta(psi) + J_theta(psi) @ Q
            + (P_r(psi) / chi_r(psi)) * (N @ J(psi) + J(psi) @ N)
            )
        )
  
    D = (chi_r(psi) ** 2) * (
        G_22(psi) + q(psi) * G_23(psi)
        + q(psi) * (G_23(psi) + q(psi) * G_33(psi))
    )
    
    E = (
        chi_r(psi) * (
            chi_rr(psi) * (G_22(psi) + q(psi) * G_23(psi))
          + qchi_r_derivative(psi) * (G_23(psi) + q(psi) * G_33(psi))
        )
        - 2 * np.pi * 1j * (chi_r(psi) ** 2) * (
            (G_12(psi) + q(psi) * G_31(psi)) @ Q
        )
        + P_r(psi) * J(psi)
    )
    
    H = (
        chi_rr(psi) * (
          chi_rr(psi) * G_22(psi) + qchi_r_derivative(psi) * G_23(psi)
        )
        + qchi_r_derivative(psi) * (chi_rr(psi) * G_23(psi) + qchi_r_derivative(psi) * G_33(psi))
        
        + 2 * np.pi * 1j * chi_r(psi) * (
            chi_rr(psi) * (Q @ G_12(psi) - G_12(psi) @ Q)
          + qchi_r_derivative(psi) * (Q @ G_31(psi) - G_31(psi) @ Q)
        )
        + ((2 * np.pi * chi_r(psi)) ** 2) * (Q @ G_11(psi) @ Q)
        + (
            P_r(psi) * ((J(psi) * chi_rr(psi) / chi_r(psi)) + J_r(psi))
          + q_r(psi) * chi_r(psi) * J_theta(psi)
        )
    )

    A_inv = np.linalg.inv(A)
    F = D - B.conj().T @ A_inv @ B
    K = E - B.conj().T @ A_inv @ C
    G = H - C.conj().T @ A_inv @ C
    F_inv = np.linalg.inv(F)
    
    L_upper = np.hstack([-F_inv @ K, F_inv])
    L_lower = np.hstack([G - K.conj().T @ F_inv @ K, K.conj().T @ F_inv])
    L = np.vstack([L_upper, L_lower])
    
    return L, F, K

我已经检查过矩阵的插值和求逆逻辑，确认这部分没有问题。

ODE系统与求解代码 #

def ode_system(psi, U_flat):

    U = np.ascontiguousarray(U_flat).reshape(2 * size, -1).view(dtype=np.complex128)
    L, _, _ = compute_matrices(psi)
    U_dot = L @ U                            # (2M, N)

    return np.ascontiguousarray(U_dot).view(dtype=np.float64).ravel()

psi_range = np.linspace(1e-1, 0.2, 20)
_, F0, K0 = compute_matrices(psi_range[0])

u0_complex = np.vstack([ # (2 * size, size)
    np.zeros((size, size), dtype=complex),
    F0 @ np.eye(size, dtype=complex)
])  

u0_real = u0_complex.view(dtype=np.float64).ravel()
ucrit = 1e4
current_t = psi_range[0]
final_t = psi_range[-1]
sing_flag = False
test_flag = False
solutions_t = []
solutions_y = []
current_u_real = u0_real.copy()

sol = solve_ivp(
        ode_system,
        [current_t, final_t],
        current_u_real,
        method='LSODA',
        atol=1e-3,
        rtol=1e-5,
        max_step=1e-3,
        min_step=1e-5,
        dense_output=True
    )
    
solutions_t.append(sol.t)
solutions_y.append(sol.y)

遇到的主要问题 #

1. 解的发散与突然归零现象
   计算解的欧几里得范数时，发现它会在某个点突然增长，然后又骤降到零，而且解的复数部分会直接变为零。这种现象的出现位置还会随积分区间变化，导致IVP问题的解不一致。即使切换到RK45、Radau或BDF求解器，这个问题依然存在。

• 积分区间t = [0.1, 0.2]的结果：
  
• 积分区间t = [0.1, 0.15]的结果：
  

2. LSODA性能异常
   LSODA理论上应该在非刚性区域大步推进，在刚性区域自动减小步长，但实际运行起来比RK45慢了近100倍。而RK45似乎能给出合理的解，我怀疑LSODA在非刚性部分并没有高效运行。

我求解的是n耦合的一阶复ODE，感觉问题出在数值计算层面，希望能找到解决办法。

备注：内容来源于stack exchange，提问作者Junyoung Jang