Python求解ODE方程组遇失败：sol.success返回False，时间点长度为1

算力阿峰

2026-6-25

求解自定义ODE方程组遇到的问题

我在求解自定义ODE方程组时遇到了问题：sol.success返回False，且print("Length of time points:", len(solution.t))输出为1。我尝试修改t_span和常数参数均无效，且确认初始条件是合理的（用于本次测试），请问该如何解决？

import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from scipy.integrate import odeint

t_span = [0, 1000]
#dt = 0.1
t_eval = np.arange(0, 1000, 1)

# constants
Mp =1
n = 1
Lambda = 1
Cupsilon =0.014

# 初始条件
phi0 = 12.327
dphi0 =-7.956* Lambda**(1/2)
rad0 =36.9622 * Lambda
a0 = 1
H0 = np.sqrt(Mp**2/2*(dphi0**2/2+phi0**(2*n-1)+rad0))
k=100*H0


J11_0=0
J12_0=0
J13_0=0
J14_0 = 0
J15_0=0
J21_0=0
J22_0=0
J23_0=0
J24_0=0
J25_0= 0
J31_0=0
J32_0=0
J33_0=0
J34_0=0
J35_0= 0
J41_0=0
J42_0=0
J43_0=0
J44_0=0
J45_0= 0
J51_0=0
J52_0=0
J53_0=0
J54_0 = 0
J55_0 = 0

# ODE方程组
def sistema_matricial_m(t, w):
    phi, dphi, rad, a, J11, J12, J13, J14, J15, J21, J22, J23, J24, J25, J31, J32, J33, J34, J35, J41, J42, J43, J44, J45, J51, J52, J53, J54,J55 = w

    dpot = Lambda * phi**(2 *n - 1)
    ddpot = Lambda * (2 * n - 1) * phi**(2 * n - 2)
    dpot0 = Lambda * phi0**(2 * n - 1)

    H = np.sqrt(Mp**2 / 2 * (dphi**2 / 2 + dpot + rad))
    H0 = np.sqrt(Mp**2 / 2 * (dphi0**2 / 2 + dpot0 + rad0))
    k = 100 * H0
    gstar = 12.5
    Cr = gstar * np.pi**2 / 30
    T = (rad / Cr)**(1 / 4)

    Alpha = 0
    Beta = 1

    gamma = Cupsilon * phi**(Alpha) * T**Beta
    gammaT = Beta*Cupsilon*(T)**(-1+Beta)
    gammaPhi = 0
    Q = gamma / (3 * H)


    fpsi = 1 + (k**2 / (3 * a**2 * H**2)) - (dphi / H)**2 / (6 * Mp**2)
    frho = 1 / (6 * Mp**2 * H**2)
    fdphi = (dphi / H) / (6 * Mp**2)
    fphi =dpot / (6 * Mp**2 * H**2)

    gpsi = gamma * H * (dphi / H)**2 - k**2 / (3 * a**2) * ((2 * Mp**2 * k**2 / (a**2 * H**2)) - (dphi / H)**2)
    grho = 4 - gammaT * H * T * (dphi / H)**2 / (4 * rad) - k**2 / (3 * a**2 * H**2)
    gdphi = -(k**2 / (3 * a**2) + 2 * gamma * H) * (dphi / H)
    gphi = -k**2 / (3 * a**2 * H**2) * (3 * (dphi / H) * H**2 + (dpot)) - H * gammaPhi * (dphi / H)**2

    hpsi = 2 * (dpot) / H**2 + gamma / H * (dphi / H)
    hdphi = 3 + gamma / H - 1 / (6 * Mp**2 * H**2) * (3 * dphi**2 + 4 * rad)
    hrho = (T * gammaT / (4 * rad * H)) * (dphi / H)
    hphi = (k**2 / (H**2 * a**2)) + ((ddpot) / (H**2)) + gammaPhi / (H) * (dphi / H)

    Grho = grho + k**2 / (3 * a**2 * H**2)
    Gpsi =gpsi + k**2 / (3 * a**2) * (2 * Mp**2 * k**2 / (a**2 * H**2) - (dphi / H)**2)
    Gphi = gphi + k**2 / (3 * a**2 * H**2) * (3 * H**2 * (dphi / H) + (dpot))
    Gdphi =gdphi + k**2 / (3 * a**2) * (dphi / H)

    A = np.array([[Grho + 4 * rad * frho, -H * k**2 / (a**2 * H**2), Gpsi + 4 * rad * fpsi, Gphi + 4 * rad * fphi, Gdphi + 4 * rad * fdphi],
                  [1 / (3 * H), 3, 4 * rad / (3 * H), gamma * dphi, 0],
                  [frho, 0, fpsi, fphi, fdphi],
                  [0, 0, 0, 0, -1],
                  [hrho + 4 * (dphi / H) * frho, 0, hpsi + 4 * (dphi / H) * fpsi, hphi + 4 * (dphi / H) * fphi, hdphi + 4 * (dphi / H) * fdphi]])

    B = np.array([[-(dphi / H) * np.sqrt(2 * gamma * T * H / a**3)], [0], [0], [0], [np.sqrt(2 * gamma * T / (a * H)**3)]])

    J = np.array([[J11, J12, J13, J14, J15], [J21, J22, J23, J24, J25], [J31, J32, J33, J34, J35], [J41, J42, J43, J44, J45], [J51, J52, J53, J54, J54]])

    matrix = ((-np.dot(A, J) - np.dot(np.transpose(A), J)) + np.dot(np.transpose(B),B))

    dphidt = dphi / H
    ddphidt = -3 * (1 + Q) * dphi - dpot / H
    draddt = -4 * rad + 3 * Q * dphi**2
    dadt = a
    eq1 = matrix[0,0]
    eq2 = matrix[0,1]
    eq3 = matrix[0,2]
    eq4 = matrix[0,3]
    eq5 = matrix[0,4]
    eq6 = matrix[1,0]
    eq7 = matrix[1,1]
    eq8 =matrix[1,2]
    eq9 = matrix[1,3]
    eq10 =matrix[1,4]
    eq11 = matrix[2,0]
    eq12 = matrix[2,1]
    eq13 = matrix[2,2]
    eq14 = matrix[2,3]
    eq15 = matrix[2,4]
    eq16 =matrix[3,0]
    eq17 = matrix[3,1]
    eq18 = matrix[3,2]
    eq19 = matrix[3,3]
    eq20 = matrix[3,4]
    eq21 =matrix[4,0]
    eq22 = matrix[4,1]
    eq23 = matrix[4,2]
    eq24 =matrix[4,3]
    eq25 = matrix[4,4]
    
    dwdt = [dphidt, ddphidt, draddt, dadt, eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16, eq17, eq18, eq19, eq20, eq21, eq22, eq23, eq24,eq25]
    return dwdt

# 初始条件数组
w0 = [phi0, dphi0, rad0, a0, J11_0, J12_0, J13_0, J14_0, J15_0, J21_0, J22_0, J23_0, J24_0,
      J25_0, J31_0, J32_0, J33_0, J34_0, J35_0, J41_0, J42_0, J43_0, J44_0, J45_0, J51_0, J52_0, J53_0, J54_0, J55_0]

# 求解ODE
sol = solve_ivp(sistema_matricial_m, t_span, w0, t_eval=t_eval, method='RK45',rtol=1e-6)
#sol=odeint(sistema_matricial_m, w0,t)

sol.success

解决步骤：

先查看求解器错误详情：执行print(sol.message)，这会直接给出失败原因，比如数值发散、出现NaN/Inf、刚性问题等，是最直接的排查手段。

修正矩阵J的笔误：代码中J矩阵最后一行的最后一个元素写成了J54，应改为J55，否则矩阵元素与初始条件不匹配，后续计算必然出错。修正后：

J = np.array([[J11, J12, J13, J14, J15], 
              [J21, J22, J23, J24, J25], 
              [J31, J32, J33, J34, J35], 
              [J41, J42, J43, J44, J45], 
              [J51, J52, J53, J54, J55]])

检查初始时刻导数合法性：手动调用sistema_matricial_m(0, w0)，查看返回的dwdt列表中是否有NaN、Inf或异常极值，重点检查除法运算的变量（如H、T、Q）是否为0或负数。
更换求解器方法：该方程组包含矩阵运算，大概率属于刚性ODE，RK45对刚性问题处理不佳，尝试使用method='Radau'或BDF，这两种方法专门针对刚性系统。
缩小初始求解范围：先将t_span改为[0,10]，验证短时间范围内是否能正常求解，排除长时间跨度导致的数值不稳定后，再逐步增大时间范围。
调整误差容忍度：暂时放宽rtol至1e-4、atol=1e-6，先让求解器运行起来，再逐步优化精度。
验证ODE方程正确性：核对每个状态变量的导数方程（如ddphidt、draddt）是否符合物理模型，检查符号和系数是否有误。

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

Python求解ODE方程组遇失败：sol.success返回False，时间点长度为1

解决步骤： #

解决步骤：