使用Scipy solve_ivp求解简单IVP时陷入无限运行的问题求助

你的代码陷入无限运行的核心原因是这个IVP存在有限时间奇点：当t=6.25时，h(t)会变为0，之后方程在实数域没有定义。solve_ivp默认会持续尝试逼近奇点，导致步长不断缩小，无法完成求解。

1. 先确认解析解与奇点位置 #

通过分离变量法可以得到该IVP的解析解：

dh/dt = -2/h
分离变量积分：$\frac{1}{2}h^2 = -2t + C$
代入初始条件h(0)=5，得$C=\frac{25}{2}$
最终解析解：$h(t) = \sqrt{25 - 4t}$

显然，当$25-4t=0$即$t=6.25$时，h(t)=0，此时方程的右端项$-2/h$趋于无穷，求解器无法继续推进时间。

2. 调整求解区间（推荐方案） #

将求解区间限制在解析解的定义域内，即t不超过6.25，修改后的代码如下：

import numpy as np
from scipy.integrate import solve_ivp

def dh_dt_zerodriver(t, h):
    return -2 / h

t_max = 6.25  # 奇点位置，解析解的最大有效时间
steps = 10
dt = t_max / steps
t_span = [0, t_max]
t_eval = np.arange(0, t_max + dt, dt)

h0 = 5

sol = solve_ivp(dh_dt_zerodriver, t_span=t_span, y0=[h0], t_eval=t_eval)
# 查看结果
print("求解时间点：", sol.t)
print("对应h值：", sol.y[0])

3. 强制终止求解（若需覆盖到t=50） #

如果必须设置t_span为[0,50]，可以通过定义终止事件让求解器在h接近0时自动停止，避免无限运行：

import numpy as np
from scipy.integrate import solve_ivp

def dh_dt_zerodriver(t, h):
    return -2 / h

# 定义终止事件：当h=0时触发停止
def stop_at_zero(t, h):
    return h[0]
stop_at_zero.terminal = True  # 设置为触发后终止求解
stop_at_zero.direction = -1   # 仅当h从正变负时触发（可选）

t = 50
steps = 10
dt = t / steps
t_span = [0, t]
t_eval = np.arange(0, t + dt, dt)

h0 = 5

sol = solve_ivp(dh_dt_zerodriver, t_span=t_span, y0=[h0], t_eval=t_eval, events=stop_at_zero)
print(f"求解自动终止于t={sol.t[-1]:.4f}，此时h={sol.y[0][-1]:.4f}")

关于数值方法的说明 #

你的论文基于RK45方法，solve_ivp默认使用的就是RK45（method='RK45'），无需更换方法。只要调整求解区间或添加终止事件，即可正常运行。

内容的提问来源于stack exchange，提问作者nicholas-t-nguyen