本文探讨了在使用 SciPy 的 quad 函数积分包含指示函数的复杂函数时可能遇到的问题，即由于 quad 的自适应性质，可能导致积分结果不准确。文章提供了问题分析，并介绍了使用 qmc_quad 函数的解决方案，该函数通过增加采样点数量来提高积分精度，并给出了详细的代码示例。
在使用 SciPy 的 quad 函数进行数值积分时，如果被积函数包含指示函数（indicator function），可能会遇到一些问题。指示函数在特定区间内取值为 1，在其他区间取值为 0，这可能导致 quad 函数的自适应积分算法失效，从而得到不准确的结果。
问题分析
scipy.integrate.quad 采用自适应求积方法，它会根据被积函数的性质动态调整采样点。如果指示函数定义的区间非常小，并且初始采样点恰好都落在指示函数值为 0 的区域，quad 可能会过早地认为积分已经收敛，从而返回不正确的结果（例如 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, np.pi] 区间上积分，得到的结果是 0。但是，如果我们指定积分区间为指示函数非零的区间 [a, b]，就能得到正确的结果。
解决方案：使用 qmc_quad
为了解决这个问题，可以使用 scipy.integrate.qmc_quad 函数。该函数使用拟蒙特卡洛（Quasi-Monte Carlo, QMC）方法进行积分，它通过在积分区间内均匀分布采样点来更准确地估计积分值。
以下是使用 qmc_quad 函数的示例代码：
import numpy as np from scipy import integrate def indac(x, xc, rad): return (xc - rad <= x) & (x <= xc + rad) phi = lambda ii, x: np.sin(ii * x) xc = 0.1586663 rad = 0.01 * np.pi # The integrand callable needs to be vectorized to evaluate # the integrand at `n_points` points in a single call. # Increase `n_points` for more accurate results. res = integrate.qmc_quad(lambda x: phi(1, x) * indac(x, xc, rad), 0., np.pi, n_points=10000) print(res) # QMCQuadResult(integral=0.009904273812591187, standard_error=1.5619537172522532e-05
代码解释
- indac 函数：定义指示函数，当 x 落在 [xc – rad, xc + rad] 区间内时返回 True (等价于1)，否则返回 False (等价于0)。需要注意的是，这里返回的是布尔值，可以被 numpy 直接转换为数值。
- phi 函数：定义被积函数的另一部分，这里是一个正弦函数。
- integrate.qmc_quad 函数：使用 QMC 方法进行积分。
- 第一个参数是被积函数，这里使用 lambda 表达式将 phi 和 indac 函数组合在一起。
- 第二个和第三个参数是积分区间的下限和上限。
- n_points 参数指定采样点的数量。增加 n_points 可以提高积分的精度，但也会增加计算时间。
注意事项
- qmc_quad 函数要求被积函数是向量化的，即可以接受一个数组作为输入，并返回一个包含每个元素积分值的数组。在上面的例子中，indac 函数已经满足这个要求。
- 选择合适的 n_points 值非常重要。如果 n_points 太小，积分结果可能不准确。如果 n_points 太大，计算时间会很长。通常需要根据具体问题进行调整。
总结
当使用 scipy.integrate.quad 积分包含指示函数的复杂函数时，可能会遇到精度问题。scipy.integrate.qmc_quad 函数提供了一种更可靠的解决方案，通过增加采样点数量来提高积分精度。在使用 qmc_quad 函数时，需要注意被积函数的向量化以及 n_points 参数的选择。