scipy 将Python中的指示器函数与quad集成

neskvpey 于 12个月前发布在 Python

关注(0)|答案(1)|浏览(217)

我正在用scipy.integrate.quad集成python中的一个指示器函数。然而，我导出的值似乎不正确。

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

字符串
这给我的是零。然而，如果我像代码中注解的那样从a到b积分，就可以获得正确的值。
它似乎与求积方法有关：在link中也有类似的发现。
任何提示将真的很感激！

scipy

来源：https://stackoverflow.com/questions/77602581/integrate-an-indicator-function-in-python-with-quad

1条答案

按热度按时间

dfty9e191#

注意，quad不会在所有可能的点或甚至 * 许多 * 点处计算被积函数。作为自适应求积方法，它近似积分并估计自己的误差，（为了效率）当它的误差估计福尔斯要求的容差时，它就终止了。它的误差估计策略对于许多表现良好的函数来说都是极好的，但是你的函数并不满足它的假设。在仅仅21个函数计算返回0.0之后，它的 error estimate 为零，所以它以0.0的积分终止。
在更紧的间隔下，它会看到被积函数为1的点，因此它会返回正确的结果。
考虑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

字符串

赞(0）回复(0）举报 12个月前

我来回答

scipy 将Python中的指示器函数与quad集成

1条答案

相关问题

热门标签

最新问答