使用 scipy.integrate.ode 检索额外的返回值

如何解决使用 scipy.integrate.ode 检索额外的返回值

我正在设置一个相当复杂的函数，我希望将其作为 ODE 求解。最好我想使用在大多数 python 环境中都可以访问的求解器，所以我尝试使用 scipy.integrate.ode。问题是我的函数产生了很多我想在每次迭代时提取的值。

这是一个例子：我想在 someValues 的每次完成迭代中获得 r.integrate() 并将它们保存在某处：

from scipy.integrate import ode

def f(t,y):
    someValues = [10,5,4]
    return -1

y0 = 1

r = ode(f)
r.set_initial_value(y0,0)

dt = 0.1
while r.successful() and r.y[0] >= 0:
    r.integrate(r.t + dt)
    # print(someValues)

然而，求解器禁止 f() 返回额外的值。由于不一定清楚求解器将如何访问该函数，因此从 f() 内部将其推送到全局变量也不是一个好主意。

虽然可以重新计算 f() 的每一步，返回所有值，但在求解器完成后，最好一次性完成以提高性能。

解决方法

我也遇到过这个问题。使用 Scipy ODE 求解器，我认为不可能按照您的要求执行，您必须“重新计算 f() 的每一步”。这是我的做法：

import numpy as np
from scipy.integrate import solve_ivp

def f_all(t,y):
    someValue = 8.0*y[0]
    dydt = y
    return dydt,someValue

def f(t,y):
    dydt,someValue = f_all(t,y)
    return dydt

def main():
    y0 = np.array([1.0])
    tspan = [0,1]
    t_eval = np.linspace(tspan[0],tspan[1],100)
    sol = solve_ivp(f,tspan,y0,t_eval=t_eval,method='LSODA')
    ysol = sol.y.T

    someValue = np.empty((len(sol.t),),np.float64)
    for i in range(len(sol.t)):
        dydt,someValue[i] = f_all(sol.t[i],ysol[i])  
    return ysol,someValue

ysol,someValue = main()

几乎所有时间，ODE 求解都比积分后执行数百或数千个函数调用花费的时间更多。但是如果你关心性能，你可以使用带有 numba 的 NumbaLSODA ode 求解器来让你的所有代码都超快：

import numpy as np
from NumbaLSODA import lsoda_sig,lsoda
import numba as nb

@nb.njit
def f_all(t,someValue

@nb.cfunc(lsoda_sig)
def f(t,y_,dydt,p):
    y = nb.carray(y_,(1,))
    dydt_,y)
    dydt[0] = dydt_[0]
funcptr = f.address

@nb.njit
def main():
    y0 = np.array([1.0])
    tspan = [0,100)
    ysol,success = lsoda(funcptr,t_eval)
    someValue = np.empty((len(t_eval),np.float64)
    for i in range(len(t_eval)):
        dydt,someValue[i] = f_all(t_eval[i],someValue = main()

上面的 NumbaLSODA 版本比 Scipy 版本快 100 倍。