python - 如何使用 solve_ivp 通过精确点？

Question

我有一个 ODE 系统，我想在求解器到达精确时间点时更改变量的值。我正在尝试做的与此类似 example in Julia

我尝试做的是使用 else 和 if 来检查时间 t 是否已到。然而，这是行不通的，因为包含时间的数组并没有经过那个确切的点。还尝试使用 t_eval 但没有成功。

使用 julia 示例中显示的相同问题，我们有:

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

def f(t, u):
    if t == 4 and u[0] < 4:
        u[0] += 10
    du = np.empty(1)
    du[0] = -u[0]
    return du

u0=[10.0]
V = 1
t = [0,10]
u0 = [10.0]
sol = solve_ivp(f,t,u0)

plt.plot(sol.t, sol.y[0])

有没有办法复制 julia 回调函数？

Answer 1

与其尝试使用回调来精确复制该代码，分两个阶段求解方程似乎更简单:

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


def f(t, u):
    du = -u
    return du


dose = 10.0

u0 = dose
V = 4

t0 = 0
t1 = 4
t2 = 10

t = np.linspace(t0, t1, 100)
sol1 = solve_ivp(f, [t0, t1], [u0], dense_output=True, t_eval=t)

u0 = sol1.y[0, -1]
if u0 / V < 4:
    u0 += dose

t = np.linspace(t1, t2, 150)
sol2 = solve_ivp(f, [t1, t2], [u0], dense_output=True, t_eval=t)

plt.plot(sol1.t, sol1.y[0], 'b')
plt.plot([sol1.t[-1], sol2.t[0]], [sol1.y[0, -1], sol2.y[0, 0]], 'b--')
plt.plot(sol2.t, sol2.y[0], 'b')
plt.grid()
plt.xlabel('t')
plt.show()