如何向量化Scipy.integrate.quad()实现区间积分采样？

我定义了一个在区间[t0, t1]上的一维函数f(t)，想要以时间步长delta_t在该区间上均匀采样，得到函数的积分结果。

目前我用的常规方法是在Python循环中调用scipy.integrate.quad：

import numpy as np
import scipy.integrate

rez = np.zeros(nStep-1)
for i in range(1, nStep):
    rez[i-1] = scipy.integrate.quad(my_func, t0, t0 + i * delta_t)[0]  # 注：quad返回(积分值, 误差)，需取第一个元素

但我觉得这不是最快的实现方式，想知道有没有向量化的实现方法？另外我自己想到一种优化思路：分段积分再累加，因为每次积分的区间更短，减少冗余计算？实现代码如下：

rez = np.zeros(nStep-1)
for i in range(1, nStep):
    rez[i-1] = scipy.integrate.quad(my_func, t0 + (i - 1) * delta_t, t0 + i * delta_t)[0]
rez = np.cumsum(rez)

首先明确：scipy.integrate.quad本身是非向量化的单区间积分函数，无法直接传入一组区间批量计算，因此不存在完全原生的向量化调用方式，但可以通过技巧减少开销，或用其他工具实现近似向量化的积分效果。

1. 分段累加思路的核心优势 #

你提出的分段积分再累加的方式，确实比原方法高效很多：

• 原方法每次都从t0积分到t0+i*delta_t，相当于重复计算了前i-1段的积分，冗余计算量随nStep增大呈指数级增长；
• 分段积分仅计算每一小段的独立积分，再用np.cumsum累加，完全避免重复计算，nStep越大，效率提升越明显。

2. 并行化加速分段积分 #

因为每一段的积分计算完全独立，可通过并行化工具利用多核CPU资源，替代串行循环：

from joblib import Parallel, delayed
import numpy as np
import scipy.integrate

# 生成所有分段的区间对
intervals = [(t0 + k*delta_t, t0 + (k+1)*delta_t) for k in range(nStep-1)]

# 并行计算每一段的积分
segment_integrals = Parallel(n_jobs=-1)(
    delayed(scipy.integrate.quad)(my_func, start, end)[0] for start, end in intervals
)

# 累加得到累计积分
rez = np.cumsum(segment_integrals)

这种方式适合nStep极大的场景，能大幅缩短计算时间；若my_func本身计算开销大，并行化的收益会更显著。

3. 完全向量化的近似积分（针对可向量化函数） #

如果你的my_func是用numpy实现的可向量化函数（能直接接受数组输入），可以用数值积分的向量化实现，比如scipy.integrate.cumtrapz或scipy.integrate.simpson：

import numpy as np
from scipy.integrate import cumtrapz

# 生成所有采样点
t = np.linspace(t0, t1, nStep)
# 计算函数在所有采样点上的值
y = my_func(t)
# 计算累计积分（cumtrapz返回nStep个值，[1:]对应原需求的nStep-1个结果）
rez = cumtrapz(y, t, initial=0)[1:]

这种方法完全无循环，速度极快，但属于数值近似积分，精度取决于采样点密度（即delta_t的大小）；而quad是自适应精度的解析积分方法，精度更高。需根据需求选择：

• 精度要求高：优先用并行化的分段quad；
• 精度要求适中且函数可向量化：用cumtrapz/simpson的向量化方法速度优势极大。

4. 关键注意点 #

• 调用scipy.integrate.quad时必须取返回值的第一个元素（积分结果），原代码遗漏这一步会导致数组存入(积分值, 误差)元组，后续计算报错；
• 若my_func计算开销极小，并行化的额外调度开销可能抵消收益，此时分段循环+cumsum就足够高效。

内容的提问来源于stack exchange，提问作者Aleksejs Fomins