1

我正在尝试计算 [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 越来越小，应该可以恢复上述结果。当我使用 quad 在 scipi 中实现这一点时，我的结果会收敛到错误的结果 6.04685e-11，正如我在 eps 与我所包含的积分结果的图中所示。
为什么quad会这样做？即使我有 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()

每股收益与积分结果

4

1 回答 1

1

通过将积分一分为二，积分到 q2 - eps 然后从 q2 + eps 开始，然后将两个结果相加，也应该能够得到相同的结果

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

我注意到您忽略了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)

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

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

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

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

于 2017-04-11T00:08:22.433 回答