当 scipy.integrate.solve_ivp 中使用的变量是另一个微分方程的解时，如何求解微分方程？

Question

我想求解定义为的微分方程

def eqn1(t, y):
    dydt = (A + B @ u) @ x
    return dydt
y = solve_ivp(eqn1, (t0, tf), y0, method='LSODA', t_eval=np.linspace(0, tf, 1000)) # this will not work cause u is unkonwn

其中

A

和

B

是已知的常数矩阵，但变量

u

是另一个方程的解，求解为

def eqn2(t, u):
    dudt = Omega + u @ S @ u
    return dudt

和

u

应该通过

获得

u = solve_ivp(eqn2, (t0, tf), u0, method='LSODA', t_eval=np.linspace(0, tf, 1000))

u

可以很容易地解决，因为eqn2中的变量都是已知的，但是如何通过解决

y

来获得

eqn1？

一旦调用

solve_ivp

，我尝试调用 eqn2 的

u

来获取

eqn1

：

def eqn2(t, u):
    dudt = Omega + u @ S @ u
    return dudt

def eqn1(t, y):
    t_eval = np.linspace(0, tf, 1000)
    u = solve_ivp(eqn2, (t0, tf), u0, method='LSODA', t_eval=np.sort(np.append(t_eval, t)))
    u = u.y[np.where(u.t == t)[0][0]]
    dydt = (A + B @ u) @ x
    return dydt
y = solve_ivp(eqn1, (t0, tf), y0, method='LSODA', t_eval=np.linspace(0, tf, 1000))

但是速度太慢并且输出很多警告。

Answer 1

您的代码的问题是您为同一个变量 u 调用了两次solve_ivp。这效率低下，还可能导致错误。

更好的方法是使用包装函数。包装函数是调用另一个函数然后返回结果的函数。在这种情况下，您可以创建一个包装函数，为 eqn2 调用solve_ivp，然后返回解。

以下代码展示了如何执行此操作：

    def wrapper(t, y):
        u = solve_ivp(eqn2, (t0, tf), u0, method='LSODA', t_eval=np.sort(np.append(t_eval, t)))
        u = u.y[np.where(u.t == t)[0][0]]
        dydt = (A + B @ u) @ x
        return dydt

y = solve_ivp(wrapper, (t0, tf), y0, method='LSODA', t_eval=np.linspace(0, tf, 1000))

此代码的工作原理是首先使用当前时间 t 和初始条件 u0 调用 eqn2 的solve_ivp。该调用的结果存储在变量 u 中。然后从 u 中提取 t 时刻的 u 值并用于计算导数 dydt。最后，包装函数返回 dydt。

此代码比原始代码更高效，因为它只对每个 t 值调用一次solve_ivp。它也不太可能导致错误，因为 u 的值仅针对每个 t 值计算一次。

我希望这有帮助！如果您还有其他问题，请告诉我。