将 scipy.quad 与 iε 技巧一起使用：不良结果

Question

为了规避柯西原理值，我尝试用一个小的位移iε将一个积分积分到复平面上以避开极点。但是，从下图可以看出，结果很糟糕。此结果的代码如下所示。您有改进此方法的想法吗？为什么它不起作用？我已经尝试过更改 ε 或积分中的极限。

编辑：我将方法 "cauchy" 包含在原理值中，但似乎根本不起作用。

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

def cquad(func, a, b, **kwargs):
    real_integral = quad(lambda x: np.real(func(x)), a, b, limit = 200,**kwargs)
    imag_integral = quad(lambda x: np.imag(func(x)), a, b, limit = 200,**kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

def k_(a):
    ϵ = 1e-32
    return (cquad(lambda x: np.exp(-1j*x)/(x**2 - a**2 - 1j*ϵ),-np.inf,np.inf)[0])

def k2_(a):
    return (cquad(lambda x: np.exp(-1j*x)/(x**2 - a**2),-1e6,1e6, weight='cauchy', wvar = a)[0])

k  = np.vectorize(k_)
k2 = np.vectorize(k2_)

fig, ax = plt.subplots()

a = np.linspace(-10,10,300)
ax.plot(a,np.real(k(a)),".-",label = "numerical result")
ax.plot(a,np.real(k2(a)),".-",label = "numerical result (cauchy)")
ax.plot(a, - np.pi*np.sin(a)/a,"-",label="analytical result")
ax.set_ylim(-5,5)
ax.set_ylabel("f(x)")
ax.set_xlabel("x")
ax.set_title(r"$\mathcal{P}\int_{-\infty}^{\infty} \frac{e^{-i y}}{y^2 - x^2}\mathrm{d}y = -\frac{\pi\sin(x)}{x}$")
plt.legend()
plt.savefig("./bad_result.png")
plt.show()

Answer 1

您可以使用 weight="cauchy" 参数代替 quad。 https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html

Answer 2

主要问题是被积函数在 x=a 和 x=-a。显示如何处理 x=a 处的极点。那么所需要的就是找到一种方法将积分按摩成避免积分通过另一极的形式在 x=-a。利用均匀性可以让我们 "fold the integral over", 所以我们只需要在 x=a.

处处理一个极点，而不是有两个极点

的实部

np.exp(-1j*x) / (x**2 - a**2) = (np.cos(x) - 1j * np.sin(x)) / (x**2 - a**2)

是 x 的偶函数，所以将 x = -infinity 的实部积分为 infinity 等于从 x = 0 到 infinity 的积分的两倍。这被积函数的虚部是x的奇函数。从x = -infinity到infinity的积分等于从x = -infinity到0的积分，加上从 x = 0 到 infinity 的积分。这两部分相互抵消因为（虚数）被积函数是奇数。所以虚部的积分等于0.

最后，使用，因为

1 / (x**2 - a**2) = 1 / ((x - a)(x + a))

使用 weight='cauchy', wvar=a 隐含地为被积函数加权 1 / (x - a) 从而允许我们将显式被积函数减少到

np.cos(x) / (x + a)

由于被积函数是 a 的偶函数，我们可以不失一般性地假设 a 是正数：

a = abs(a)

现在从 x = 0 积分到 infinity 避开了 x = -a 处的极点。

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


def cquad(func, a, b, **kwargs):
    real_integral = quad(lambda x: np.real(func(x)), a, b, limit=200, **kwargs)
    imag_integral = quad(lambda x: np.imag(func(x)), a, b, limit=200, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])


def k2_(a):
    a = abs(a)
    # return 2*(cquad(lambda x: np.exp(-1j*x)/(x + a), 0, 1e6, weight='cauchy', wvar=a)[0]) # also works
    # return 2*(cquad(lambda x: np.cos(x)/(x + a), 0, 1e6, weight='cauchy', wvar=a)[0]) # also works, but not necessary
    return 2*quad(lambda x: np.cos(x)/(x + a), 0, 1e6, limit=200, weight='cauchy', wvar=a)[0]


k2 = np.vectorize(k2_)

fig, ax = plt.subplots()

a = np.linspace(-10, 10, 300)
ax.plot(a, np.real(k2(a)), ".-", label="numerical result (cauchy)")
ax.plot(a, - np.pi*np.sin(a)/a, "-", label="analytical result")
ax.set_ylim(-5, 5)
ax.set_ylabel("f(x)")
ax.set_xlabel("x")
ax.set_title(
    r"$\mathcal{P}\int_{-\infty}^{\infty} \frac{e^{-i y}}{y^2 - x^2}\mathrm{d}y = -\frac{\pi\sin(x)}{x}$")
plt.legend()
plt.show()

将 scipy.quad 与 iε 技巧一起使用：不良结果

Using scipy.quad with iε trick: Bad results

python

integral

scipy

numerical-methods

quad