振荡函数的积分在python中不收敛

Question

我想使用 scipy.integrate.quad 评估函数的积分。以下是被积函数的样子： 我们可以注意到，该被积函数的大部分贡献将来自 0.1 到 4 或 5 ish。对于 x = 10 及以上，虽然函数是振荡的（在图片上很难分辨），但它非常小并且不断变小。

这是从 0 到某个上限的积分结果的样子。积分的上限在 x 轴上。 在这里，虽然前一百个 x 的结果似乎是稳定的，但之后就不再是这样了，即使我期待的是一条直线......

我对 python 很陌生，我不知道最好的做法是什么。我现在最好的猜测是取一些小于 100 的值作为我的积分的上限，并丢弃其他值，因为它只是从 integration.quad 收敛的不好。

编辑： 为了绘制第二张图，我使用了 scipy.integrate.quad 函数。但是，如果我只使用我生成的点来绘制被积函数（第一个图）并在 scipy.integrate.simps 中使用它并改变我积分的最大 x，我会得到一致的结果。

Answer 1

当被积函数具有远小于积分范围的重要特征时，它可能会被自适应quad 例程“忽略”。相反，simps 如果您使用足够细的网格，则不会错过它们，但可能需要更长的时间来评估。我将描述两种处理方法，第二种更实用。

积分参数

您可以使用quad 的points 参数来确保不会发生这种情况。这是一个示例，我在区间 [-1000, 5000] 上集成了高斯函数 exp(-x**2)。该函数定位在 0 附近；几乎所有这些都在 [-5,5] 区间内，所以我包括 points=[-5, 5] 以确保不会忽略此范围。 （积分要求在积分范围内，所以出现if）。

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

f = lambda x: np.exp(-x**2) 
numpoints = 1000
t = np.linspace(-1000, 5000, numpoints)
y = np.zeros((numpoints,))

for i in range(numpoints):
    y[i] = quad(f, -1000, t[i])[0]      # integration without points
plt.plot(t, y, 'r') 

for i in range(numpoints):            # using points if upper bound is above 5 
    y[i] = quad(f, -1000, t[i], points=[-5,5])[0] if t[i] > 5 else quad(f, -1000, t[i])[0]
plt.plot(t, y, 'b') 
plt.show()

红色曲线是没有points的输出，蓝色曲线是带有points的输出。后者的行为应有尽有：从 0 上升到 pi/2 并保持在那里。

重用之前的计算

在多个点计算反导数的一种更有效的方法是使用先前计算的值，将未积分的区间的贡献添加到该值中。

y = np.zeros((numpoints,))
for i in range(1, numpoints):
    y[i] = y[i-1] + quad(f, t[i-1], t[i])[0]
plt.plot(t, y, 'g') 

这与上面的蓝色曲线具有相同的输出。