使用scipy.integrate.odeint覆盖ODE集成中的值

编程问答 2022-08-15

问题描述

我正在使用scipy.integrate.odeint对ODE系统进行数值积分。 抽象地讲，我希望能够给求解器一个任意条件以应用于其输出值。例如：如果一项的值（例如X[0]）小于1，将值覆盖为0，而不是使用通过数值积分计算的值。然后进行整合。

这是我如何进行整合的一个例子。

# Toy ODE system to be integrated
def dXdt(X,t,args):
    return [args[0]*X[0],t*args[1]*X[1]]

# Setup initial conditions and integrate the system.
from scipy.integrate import odeint
import numpy as np
t = np.linspace(0,100,100) # timepoints to integrate
x0 = [100,100] # initial conditions
args = [-.01,-.02]
X = odeint(dXdt,x0,args=(args,))

一种大致实现我想要的行为的解决方法是按如下方式修改dXdt：

def dXdt(X,args):
    if X[0] > 1: # make dX[0]/dt extremely negative to force X[0] toward 0
        return [-1e200,t*args[1]*X[1]] 
    return [args[0]*X[0],t*args[1]*X[1]]

这显然不能推广到可能需要的其他条件，也不能精确解决问题（精确的解决方案将需要知道dXdt中求解器的步长。

解决方法

这是不可能的。您必须将c(10,25)视为黑盒，在这里您无权访问内部状态和集成期间执行的步骤。 在其他上下文（语言，库）中，您可以使用事件操作机制，您可以在其中定义向下穿越某个值的事件，并采取操作来修改状态。

在odeint中，您实现了事件机制的一半。您将必须将事件定义为终端，提取事件的最后状态，进行修改，然后以修改后的状态重新开始集成。最后，您可以连接解决方​​案细分。

numpyodeodeintpythonscipyscipy