使用 SciPy 的四边形通过积分到奇异点的正下方和正上方来获得积分的主值

Question

我正在尝试计算 [Ecut, inf] 上 1/((s - q02)*(s - q2)) 的积分（在 s 上）的主值，其中 q02 < Ecut < q2。手工（或 Mathematica）计算 Principle 值得到一般结果

ln((q2-Ecut)/(Ecut-q02)) / (q02 -q2)

在下面的具体示例中，这给出了结果 -1.58637*10^-11。人们也应该能够通过将积分一分为二，积分到 q2 - eps 然后从 q2 + eps 开始，然后将两个结果相加（分歧应该抵消）来获得相同的结果。通过使 eps 越来越小，应该可以恢复上面的结果。当我在 scipi 中使用 quad 实现它时，我的结果收敛到错误的结果 6.04685e-11，正如我在 eps 与我包含的积分结果的图中所示。
为什么四方这样做？即使我有 eps = 0 它也会给我这个错误的结果，当我期望它会在事情爆炸时给我一个错误...

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


q02  = 485124412.
Ecut = 17909665929.
q2   = 90000000000.

def integrand(s):
    return 1/((s - q02)*(s - q2))

xx=[1.,0.1,0.01,0.001,0.0001,0.00001,0.000001,0.0000001,0.00000001,
    0.000000001,0.0000000001,0.00000000001,0.]

integral = [0*y for y in xx]
i=0
for eps in xx:

    ans1,err = quad(integrand, Ecut, q2 -eps )
    ans2,err= quad(integrand, q2 + eps, np.inf)

    integral[i] = ans1 + ans2
    i=i+1

plt.semilogx(xx,integral,marker='.')
plt.show()

Answer 1

One should also be able to get the same result by splitting the integral in two, integrating up to q2 - eps and then starting from q2 + eps, and then adding the two results

只有在计算完全准确的情况下。在数值实践中，你所描述的基本上是一个人能做的最糟糕的事情。你得到两个符号相反的大积分，它们相加时几乎相互抵消；剩下的更多是与积分误差有关，而不是与积分的实际值有关。

我注意到您忽略了脚本中的错误值 err，甚至没有将它们打印出来。 坏主意：它们的大小为 1e-10，这已经告诉您 "something e-11" 的最终结果是垃圾。

计算科学问题 Numerical Principal Value Integration - Hilbert like 解决了这个问题。他们指出的一种方法是在尝试对奇点求积分之前，在关于奇点对称的点处添加被积函数的值。这需要对以奇点q2为中心的对称区间（即从Ecut到2*q2-Ecut）进行积分，然后加上积分从2*q2的贡献-Ecut 到无穷大。无论如何，这种分裂是有道理的，因为 quad 对待无限极限的方式非常不同（使用傅里叶积分），这是影响奇点抵消方式的另一件事。

因此，这种方法的实现是

ans1, err = quad(lambda s: integrand(s) + integrand(2*q2-s), Ecut, q2)
ans2, err = quad(integrand, 2*q2-Ecut, np.inf)

不需要每股收益。然而，结果仍然不对：大约 -2.5e-11。事实证明，第二个积分是罪魁祸首。不幸的是，傅里叶积分方法在这里似乎并不有效（或者我没有找到让它起作用的方法）。事实证明，提供一个大但有限的值作为上限会导致更好的结果，特别是如果还使用选项 epsabs，例如epsabs=1e-20。

更好的是，仔细阅读 documentation of quad 并注意它直接支持柯西权重为 1/(s-q2) 的积分，为它们选择合适的数值方法。这仍然需要一个有限的上限，以及 epsabs 的一个小值，但结果非常准确：

quad(lambda s: 1/(s - q02), Ecut, 1e9*q2, weight='cauchy', wvar=q2, epsabs=1e-20)

returns -1.5863735715967363e-11，与精确值 -1.5863735704856253e-11 相比。请注意，因子 1/(s-q2) 没有出现在上面的被积函数中，被降级为权重选项。

使用 SciPy 的四边形通过积分到奇异点的正下方和正上方来获得积分的主值

Using SciPy's quad to get the principal value of an integral by integrating to just below and from just above the singular point

scipy

numerical-integration

quad