使用solve_ivp时有没有办法存储中间值？

Question

我正在研究一个用 scipy.integrate.solve_ivp 求解的微分方程系统。对于这个系统的每个变化项，我结合了不同的中间变量。我正在寻找一种方法来存储这些中间变量，而不仅仅是结果变量。

我正在考虑使用列表或其他类型的可迭代将中间值直接存储在传递给 solve_ivp 的 func 中，但这似乎不是最好的方法，并且不允许我使用插值，当使用密集输出时。 我也可以使用solve_ivp的结果重新计算这些中间变量，但是由于计算有点复杂，直接存储值并且只计算一次应该更容易。

我的 func 函数看起来像这样：

def dY(t, Y):
    C1, C2 = Y

    dC1 = A(C1, C2) + B(C1, C2)
    dC2 = C(C1, C2) + B(C1, C2)

    return [dC1, dC2]

我想在每一步存储 A、B、C 和 D 的值，以及 C1 和 C2 的值（对于它所说的人来说，这是一个化学工程问题，A、B、C 和D 是不同的汇和源术语）

感谢您的帮助！

Answer 1

使用解的值并在事后调用函数。不要使用存储在dY 函数中的中间值。僵硬的 ODE 积分器将在一个时间步内多次调用该函数，最终将存储无意义的数据。非刚性 ODE 积分器可能还会在一次尝试中尝试多个时间步长，但会遇到相同的问题。

以下是事后如何重构 A 和 B 的示例：

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

tspan = (0, 10)
y0 = [1, 2]

def A(C1, C2):
    return 2*C1 + 3*C2 - C1**2
def B(C1, C2):
    return 1*C1 + 4*C2
def C(C1, C2):
    return 1*C1 + 3*C2 - C2**2
def D(C1, C2):
    return 2*C1 + 5*C2
def dY(t, Y):
    C1, C2 = Y

    dC1 = A(C1, C2) + B(C1, C2)
    dC2 = C(C1, C2) + B(C1, C2)

    return [dC1, dC2]
sol = solve_ivp(dY, tspan, y0)

fig,ax = plt.subplots()
ax.plot(sol.t, sol.y[0], label='C1')
ax.plot(sol.t, sol.y[1], label='C2')
ax.legend()
fig2,ax2 = plt.subplots()
ax2.plot(sol.t, A(sol.y[0], sol.y[1]), label='A')
ax2.plot(sol.t, B(sol.y[0], sol.y[1]), label='B')
ax2.legend()