python中带奇异点的数值积分（主值）

Question

我试图用scipy.integrate中的四元函数来积分一个具有奇点的函数，但是我没有得到想要的答案。代码如下：

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

这将导致关于0.4的集成警告和返回值。在

函数的极点是[-1,1]。答案应该是1.09左右（用笔和纸计算）。在

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进行符号集成：

^{pr2}$

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

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和{}是极限，poles是极点列表，eps是接近极点的距离。你可以让每股收益越来越小，以获得更好的结果，但也许sympy会更好地解决这样的问题。在

有了这个函数和标准eps我得到了1.0986112886023367，这与wolframalpha给出的结果几乎相同。在