bounty将在17小时后过期。回答此问题可获得+50声望奖励。Rafnus希望引起更多人对此问题的关注：我并不关心任何具体的集成方法，我只是想把下面给出的代码加速一个数量级。这应该是绝对可能的，我只是不知道怎么做。任何人谁管理的加速水平是有资格的赏金。

我想通过结合Numba的cfunc和Scipy的LowLevelCallable功能来优化基于时间的ODE系统的求解，以大大提高odeint()的速度。但是，我很难找到正确解的确切语法。
我的（未优化的）代码看起来像这样：

import numpy as np
import numba as nb
from numba import types
from matplotlib import pyplot as plt
import scipy.integrate as si

# hours in the simulation
T_sim = 72 
# simulate fluctuating outside temp for 3 days
T_out = np.sin(np.linspace(0, T_sim, T_sim)/3.8) * 8 + 18 
# heat capacity of the air and floor
cap_air, cap_flr = 1, 10

def func(T, t):
    T_air, T_flr= T[0], T[1]
    t_ = int(t)
    H_airOut = 4 * U_vent[t_] * (T_air - T_out[t_])
    H_blowAir = U_blow[t_]
    H_airFlr = (T_air - T_flr)
    dT_air = ((H_blowAir 
               - H_airOut
               - H_airFlr) / cap_air)
    dT_flr = H_airFlr / cap_flr
    return dT_air, dT_flr

T0 = [16, 18] # starting temperature for T_air and T_flr
t_eval = np.linspace(0, T_sim-1, T_sim-1)
U_blow = np.full(T_sim, 0.1) # control for the heating
U_vent = np.full(T_sim, 0.3) # control for the ventilation
result = si.odeint(func,T0,t_eval)

# this piece of code is a standin for the real code, where func() is run 
# tens of thousands of times. All variables are constant, except for the 
# two arrays U_blow and U_vent
for _ in range(10):
    U_blow = np.random.rand(T_sim)
    U_vent = np.random.rand(T_sim)
    result = si.odeint(func,T0,t_eval)

plt.plot(result[:, 0], label='T_air')
plt.plot(result[:, 1], label='T_flr')
plt.plot(T_out[:-1], label='T_out')
plt.legend()

根据其他SO帖子，我非常确定解决方案应该看起来像下面的代码，但我自己不能让它工作。如果只有一些变量在运行之间改变，问题会容易得多，但这不是事实，因为我需要能够在运行之间改变数组U_blow和U_vent。

import scipy as sp
from numba.types import float64, CPointer, intc

def jit_integrand_function(integrand_function):
    jitted_function = nb.jit(integrand_function, nopython=True)

    @nb.cfunc(float64(intc, CPointer(float64)))
    def wrapped(n, xx):
        return jitted_function(xx[0], xx[1])
    return sp.LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def func(T, t):
    T_air, T_flr= T[0], T[1]
    t_ = int(t)
    H_airOut = 4 * U_vent[t_] * (T_air - T_out[t_])
    H_blowAir = U_blow[t_]
    H_airFlr = (T_air - T_flr)
    dT_air = ((H_blowAir 
               - H_airOut
               - H_airFlr) / cap_air)
    dT_flr = H_airFlr / cap_flr
    return dT_air, dT_flr

T0 = [16, 18] # starting temperature for T_air and T_flr
t_eval = np.linspace(0, T_sim-1, T_sim-1)

for _ in range(10):
    U_blow = np.random.rand(T_sim)
    U_vent = np.random.rand(T_sim)
    result = si.odeint(func,T0,t_eval, (U_blow, U_vent))

plt.plot(result[:, 0], label='T_air')
plt.plot(result[:, 1], label='T_flr')
plt.plot(T_out[:-1], label='T_out')
plt.legend()

顺便说一句，我的完整模型在数值上不够稳定，欧拉前向积分法无法工作，所以这不是一个选择。