我们可以使用 VODE 使用 scipy.integrate.ode 修改集成步骤之间的解向量吗？

如何解决我们可以使用 VODE 使用 scipy.integrate.ode 修改集成步骤之间的解向量吗？

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

输出：

到目前为止一切顺利。

现在的问题是，在每一步中，我想在VODE集成后修改y。自然，我希望 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 是 scipy.integrate.ode 的源代码，here 是 scipy 用于 VODE 的 pyf 接口。

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

非常感谢！

解决方法

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

我的解决方案是将使用的包从 ode 更改为 solve_ivp。不同之处在于 solve_ivp 完全是用 Python 制作的，您将能够通过自己的方式进行实现。请注意，与其他包使用的 vode 链接相比，代码运行速度会很慢，即使代码编写得很好并且使用了 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,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,y_list)

给出情节：

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

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