与 c++ 相比，与 solve_ivp 集成非常慢答案

【问题标题】：Integrating with solve_ivp is very slow in comparison with c++与 c++ 相比，与 solve_ivp 集成非常慢
【发布时间】：2019-07-26 18:52:06
【问题描述】：

使用 scipy.integrate 中的 solve_ivp 来集成刚性系统大约需要 3 分钟才能完成。积分器是：

sol = solve_ivp(DC_model,[t0,tf],y0,method='LSODA')

使用 c++ 和 Boost 库求解同一系统需要 2.12742 秒。这是一个巨大的差异。有没有办法改善我的 python 脚本的时间？

完整代码：

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp 
import time

def DC_model(t,y):
    T,H,CTL,Den,IL2,FBeta,FGamma,Ml = y

    dT = 0.002 * T * np.log(1e12 / T) - (0.1136 * T * CTL * Ml / (50 + Ml)) * ((0.69 * FBeta + 1e4) / (FBeta + 1e4))
    dH = 1e-4 - 0.005 * H + 10e-2 * Den*(H*(1 - H / 1))
    dC = 1e-4 - 0.01925 * CTL + 0.00004e-2 * IL2*(CTL*(1 - CTL / 1))
    dDen = -0.009625 * Den * CTL
    dI = 1e-2 * H * Den - 1e-7 * CTL * IL2 - 1e-2 * IL2
    dFbeta = 5.57e-6 * T - 6.93 * FBeta
    dFgamma = 1.02e-4 * CTL - 0.102 * FGamma
    dMl = 1.44 + (2.89 * FGamma) / (FGamma + 3.38e5) - 0.0144 * Ml

    return np.array([dT,dH,dC,dDen,dI,dFbeta,dFgamma,dMl])


y0 = [6e4, 0, 0, 0, 0, 0, 0, 0]
t0 = 0.0 
delay = 232
tf = 168 + delay
ef=0.05
start = time.time()
sol = solve_ivp(DC_model,[t0,tf],y0,method='LSODA')

y0 = sol.y[:,-1] + [0, 0, 0, 1e6 * ef, 0, 0, 0, 0]
t0 = tf 
tf = tf + 168   
sol = solve_ivp(DC_model,[t0,tf],y0,method='LSODA')

y0 = sol.y[:,-1] + [0, 0, 0, 1e6 * ef, 0, 0, 0, 0]
t0 = tf 
tf = tf + 168   
sol = solve_ivp(DC_model,[t0,tf],y0,method='LSODA')

y0 = sol.y[:,-1] + [0, 0, 0, 1e6 * ef, 0, 0, 0, 0]
t0 = tf 
tf = 1400   
sol = solve_ivp(DC_model,[t0,tf],y0,method='LSODA')
end = time.time()
print('Solving took: ' + str((end-start)/60) + ' min')

【问题讨论】：

可悲的是，这是一个合理的减速，人们经常观察到编写良好的 Python 与编写良好的 C++ 相比。 numpy 有一个 ODE 求解器，它使用 C 扩展模块进行计算，这可能会对您有所帮助。
不幸的是，这是一个未解决的问题 (github.com/scipy/scipy/issues/8257)。你最好使用 ode 或 odeint，这会更快，尽管当我用你的代码尝试解决方案时，解决方案有很大不同。
尝试 cythonise 你的 DC_model 函数。

标签： python scipy numerical-integration

【解决方案1】：

我为 LSODA 编写了一个 python 包装器，它应该和 C/C++ 一样快：https://github.com/Nicholaswogan/NumbaLSODA。

不过我有点困惑，因为当我在上面运行您的代码时，它使用solve_ivp 执行得非常快。无论如何，这是您使用 NumbaLSODA 的代码：

import matplotlib.pyplot as plt
import numpy as np 
import time

from NumbaLSODA import lsoda_sig, lsoda
import numba as nb


@nb.cfunc(lsoda_sig)
def DC_model(t, y_, dy, p):
    y = nb.carray(y_, (8,))

    
    T,H,CTL,Den,IL2,FBeta,FGamma,Ml = y

    dT = 0.002 * T * np.log(1e12 / T) - (0.1136 * T * CTL * Ml / (50 + Ml)) * ((0.69 * FBeta + 1e4) / (FBeta + 1e4))
    dH = 1e-4 - 0.005 * H + 10e-2 * Den*(H*(1 - H / 1))
    dC = 1e-4 - 0.01925 * CTL + 0.00004e-2 * IL2*(CTL*(1 - CTL / 1))
    dDen = -0.009625 * Den * CTL
    dI = 1e-2 * H * Den - 1e-7 * CTL * IL2 - 1e-2 * IL2
    dFbeta = 5.57e-6 * T - 6.93 * FBeta
    dFgamma = 1.02e-4 * CTL - 0.102 * FGamma
    dMl = 1.44 + (2.89 * FGamma) / (FGamma + 3.38e5) - 0.0144 * Ml

    dy_ = np.array([dT,dH,dC,dDen,dI,dFbeta,dFgamma,dMl])
    for i in range(len(dy_)):
        dy[i] = dy_[i]

funcptr = DC_model.address

@nb.njit
def main():
    y0 = np.array([6e4, 0, 0, 0, 0, 0, 0, 0],np.float64)
    t0 = 0.0 
    delay = 232
    tf = 168 + delay
    ef=0.05
    t_eval = np.linspace(t0,tf,100)
    sol, success = lsoda(funcptr, y0, t_eval)

    y0 = sol[-1] + np.array([0, 0, 0, 1e6 * ef, 0, 0, 0, 0],np.float64)
    t0 = tf 
    tf = tf + 168   
    t_eval = np.linspace(t0,tf,100)
    sol, success = lsoda(funcptr, y0, t_eval)

    y0 = sol[-1] + np.array([0, 0, 0, 1e6 * ef, 0, 0, 0, 0],np.float64)
    t0 = tf 
    tf = tf + 168 
    t_eval = np.linspace(t0,tf,100)
    sol, success = lsoda(funcptr, y0, t_eval)

    y0 = sol[-1] + np.array([0, 0, 0, 1e6 * ef, 0, 0, 0, 0],np.float64)
    t0 = tf 
    tf = 1400   
    t_eval = np.linspace(t0,tf,100)
    sol, success = lsoda(funcptr, y0, t_eval)
   
main() # to compile the code
start = time.time()
main()
end = time.time()
print('Solving took: ' + str((end-start)) + ' sec')

结果

Solving took: 0.00039005279541015625 sec

【讨论】：