Python中的奇异数值积分（主值）

Question

我正在尝试使用scipy.integrate中的quad函数集成一个具有奇点的函数，但我没有得到期望的答案。以下是代码：

from scipy.integrate import quad
import numpy as np
def fun(x):
    return 1./(1-x**2)
quad(fun, -2, 2, points=[-1, 1])

这会导致IntegrationWarning和关于0.4的返回值。

该函数的极点为[-1,1]。答案应该大约为1.09（使用笔和纸计算得出）。

- user10289025

---

有趣的是，2*quad(fun, 0, 2, points=[1])[0]) 给出了正确的结果。 - AGN Gazer

Answer 1

选项weight='cauchy'可用于有效计算发散积分的主值，例如此类积分。这意味着提供给quad函数的函数将被隐式乘以1/(x-wvar)，因此请相应地调整该函数（通过乘以x-wvar其中wvar是奇异点的位置）。

i1 = quad(lambda x: -1./(x+1), 0, 2, weight='cauchy', wvar=1)[0]
i2 = quad(lambda x: -1./(x-1), -2, 0, weight='cauchy', wvar=-1)[0]
result = i1 + i2

结果为1.0986122886681091。

像这样的简单函数，您还可以使用SymPy进行符号积分：

from sympy import symbols, integrate
x = symbols('x')
f = integrate(1/(1-x**2), x)
result = (f.subs(x, 2) - f.subs(x, -2)).evalf()

结果：1.09861228866811。如果没有使用evalf()，它将是log(3)。

- user6655984

---

我认为你想说的是“发散积分项”而不是“发散积分”。 - AGN Gazer

关于将积分拆分为两个部分：由于被积函数是偶函数，它关于Y轴对称。这意味着您只能在 [0,2] 范围内积分-i1，并将 [-2,0] 范围内的积分 -i2-设置为等于 i1。也就是说，您可以执行 result = 2 * i1，而不是现在计算两个积分。 - AGN Gazer

如果我想要一个权重为1/(x-wvar)**2，我应该如何处理？ - user10289025

@Murali 这些不是支持的权重。带有这样因子的积分可能发散。 - user6655984

@Murali 可以将因子包含在被积函数中，例如可以通过使用Cauchy权重对(x-2)*sin(1/(x-2))进行积分来积分sin(1/(x-2))。此外，另一个答案建议采用去除极点的epsilon邻域的合理方法来提高收敛性。 - user6655984

Answer 2

我也无法使用原始函数使其正常工作。我设计了以下代码，使用scipy计算主值：

def principal_value(func, a, b, poles, eps=10**(-6)):
    #edges
    res = quad(func,a,poles[0]-eps)[0]+quad(func,poles[-1]+eps,b)[0]
    #inner part
    for i in range(len(poles)-1):
        res += quad(func, poles[i]+eps, poles[i+1]-eps)[0]
    return res

当您的函数句柄为func时，a和b是极限值，poles是极点列表，eps是您想要接近极点的程度。您可以将eps越来越小以获得更好的结果，但对于这样的问题，sympy可能会更好。

使用此函数和标准的eps，我得到1.0986112886023367的结果，这几乎与wolframalpha给出的结果相同。

- user8408080