使用scipy.integrate.quad积分指示函数结果异常求助

算力阿峰

2026-7-4

问题：scipy.integrate.quad积分指示函数结果异常

使用scipy.integrate.quad对指示函数与正弦函数的乘积积分时，出现区间[0, π]积分结果为0，但在指示函数非零区间内积分得到正确值的异常情况。

复现代码

import numpy as np
from scipy.integrate import quad

def indac(x, xc, rad):
    if xc - rad <= x <= xc + rad:
        return 1
    else:
        return 0

phi = lambda ii, x: np.sin(ii * x)

xc = 0.1586663
rad = 0.01 * np.pi

# 全区间积分结果异常
result, _ = quad(lambda x: phi(1, x) * indac(x, xc, rad), 0., np.pi)
print(result)  # 0.0

# 缩小到非零区间积分结果正确
a, b = xc - rad, xc + rad
result, _ = quad(lambda x: phi(1, x) * indac(x, xc, rad), a, b)
print(result)  # 0.009925887836572549

异常现象

积分区间为[0, π]时，输出结果为0.0，与实际不符
积分区间缩小到指示函数非零的[a, b]时，输出正确结果

解决提示

问题根源：quad默认的自适应求积算法对阶跃不连续函数的检测能力有限，当非零区间相对于整个积分区间极小时，算法的采样点可能完全错过非零区域，导致积分结果为0。
针对性解决方案：
- 手动指定间断点：调用quad时通过points参数传入指示函数的间断点，强制算法在这些位置细化采样：
```
result, _ = quad(lambda x: phi(1, x) * indac(x, xc, rad), 0., np.pi, points=[a, b])
```
- 向量化指示函数：替换Python原生if-else分支为NumPy向量化逻辑，提升算法对函数特征的识别能力：
```
def indac(x, xc, rad):
    return np.where((xc - rad <= x) & (x <= xc + rad), 1, 0)
```
- 拆分积分区间：手动将[0, π]拆分为[0, a]、[a, b]、[b, π]三个区间分别积分后求和，前、后区间积分结果为0，仅需计算中间有效区间。