问题描述
我正在设置一个相当复杂的函数,我希望将其作为 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 倍。