我们可以使用 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。

ode python scipy scipy