scipy 的solve_ivp 无法解决简单的ivp

Question

我正在尝试用Python解决一个非常简单的IVP，以便进行错误分析：

import numpy as np
from scipy.integrate import solve_ivp


def dh_dt_zerodriver(t, h):
    return -2 / 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)

但是，该解决方案不会计算，因为它会无限期地运行。我已经使用solve_ivp在没有解析解的情况下解决了更复杂的IVP，但这似乎让我难住了，尽管它看起来是一个相当简单的问题。

我的论文基于使用 RK45 方法，但如果我必须使用某种不同的数值方法，那么也可以。

Answer 1

问题是您想要解决方案的领域。您的 ODE 是

h'(t) = -2/h(t)

，其中

h(0) = 5

。解析解为

h(t) = sqrt(25 - 4t)

。对于

t < 6.25

，解是实数，但对于

t >= 6.25

，解变得复杂。如果您改为对

t = [0, 6]

进行积分，您很快就会得到正确的结果。