使用 SciPy quad 积分指示函数：问题与解决方案

本文探讨了在使用 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

AI建筑知识问答

用人工智能ChatGPT帮你解答所有建筑问题

22

查看详情

为了解决这个问题，可以使用 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

代码解释

1. indac 函数：定义指示函数，当 x 落在 [xc - rad, xc + rad] 区间内时返回 True (等价于1)，否则返回 False (等价于0)。需要注意的是，这里返回的是布尔值，可以被 numpy 直接转换为数值。
2. phi 函数：定义被积函数的另一部分，这里是一个正弦函数。
3. 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 参数的选择。

以上就是使用 SciPy quad 积分指示函数：问题与解决方案的详细内容，更多请关注php中文网其它相关文章！