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

Question

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

我在考虑使用列表或其他某种可迭代对象将中间值直接存储在传递给 solve_ivp 的函数中，但这似乎不是最好的方法，也不是允许我在使用 dense_output 时使用插值。我也可以使用 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()

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

Is there a way to store intermediate values when using solve_ivp?

python

scipy

numerical-methods

differential-equations