我们可以使用 VODE 修改具有 scipy.integrate.ode 的积分步骤之间的解决方案向量吗？

Question

我正在尝试为刚性 ODE 问题找到一个解决方案，在每个积分步骤中，我必须在继续积分之前修改解向量。为此，我在 bdf 模式下使用 scipy.integrate.ode 和积分器 VODE。这是我正在使用的代码的简化版本。功能比那个复杂很多，涉及到CANTERA的使用。

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

def yprime(t,y):
    return y

vode = ode(yprime)
vode.set_integrator('vode', method='bdf', with_jacobian=True)

y0 = np.array([1.0])
vode.set_initial_value(y0, 0.0)
y_list = np.array([])
t_list = np.array([])
while vode.t<5.0 and vode.successful:
    vode.integrate(vode.t+1e-3,step=True)
    y_list = np.append(y_list,vode.y)
    t_list = np.append(t_list,vode.t)

plt.plot(t_list,y_list)

输出：

到目前为止一切顺利。

现在的问题是，在每个步骤中，我想在 y 被 VODE 集成后修改它。当然，我希望 VODE 继续与修改后的解决方案集成。这是我到目前为止所尝试的：

while vode.t<5.0 and vode.successful:
    vode.integrate(vode.t+1e-3,step=True)
    vode.y[0] += 1  # Will change the solution until vode.integrate is called again
    vode._y[0] += 1 # Same here.

我也试过查看 vode._integrator，但似乎一切都保存在求解器的 fortran 实例中。作为快速参考，here is the source code of scipy.integrate.ode, and here 是 pyf 接口 scipy 用于 VODE。

有没有人试过类似的东西？我也可以更改我正在使用的求解器和/或包装器，但我想为此继续使用 python。

非常感谢！

Answer 1

对于遇到同样问题的人，问题出在 Scipy 的 Fortran 包装器中。

我的解决方案是更改使用的包，从 ode 到 solve_ivp。不同之处在于 solve_ivp 完全是用 Python 制作的，您将能够通过实施破解您的方式。请注意，与其他包使用的 vode link 相比，代码将运行变慢，即使代码编写得很好并且使用了 numpy（基本上，C 级性能尽可能).

以下是您必须遵循的几个步骤。

首先，重现已经工作的代码：

from scipy.integrate import _ivp # Not supposed to be used directly. Be careful.
import numpy as np
import matplotlib.pyplot as plt

def yprime(t,y):
    return y

y0 = np.array([1.0])
t0 = 0.0
t1 = 5.0

# WITHOUT IN-BETWEEN MODIFICATION
bdf = _ivp.BDF(yprime,t0,y0,t1)

y_list = np.array([])
t_list = np.array([])
while bdf.t<t1:
    bdf.step()
    y_list = np.append(y_list,bdf.y)
    t_list = np.append(t_list,bdf.t)

plt.plot(t_list,y_list)

输出：

现在，实现一种在积分步骤之间修改 y 值的方法。

# WITH IN-BETWEEN MODIFICATION
bdf = _ivp.BDF(yprime,t0,y0,t1)

y_list = np.array([])
t_list = np.array([])
while bdf.t<t1:
    bdf.step()
    bdf.D[0] -= 0.1 # The first column of the D matrix is the value of your vector y.
    # By modifying the first column, you modify the solution at this instant.
    y_list = np.append(y_list,bdf.y)
    t_list = np.append(t_list,bdf.t)

plt.plot(t_list,y_list)

给出剧情：

不幸的是，这对这个问题没有任何物理意义，但暂时有效。

注意：求解器完全有可能变得不稳定。它与 Jacobian 没有在正确的时间更新有关，因此必须再次重新计算它，这在大多数情况下都是性能很高的。好的解决方案是重写 class BDF 以在更新雅可比矩阵之前实施修改。源代码 here.

我们可以使用 VODE 修改具有 scipy.integrate.ode 的积分步骤之间的解决方案向量吗？

Can we modify the solution vector between integrations steps with scipy.integrate.ode, using VODE?

python

scipy

ode