为了规避柯西原理值，我尝试将使用小位移 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()
主要问题是被积函数在两个位置都有极点x=a and
x=-a. ev-br's post展示如何
处理一根杆子x=a。那么所需要的就是找到一种方法
将积分调整为避免通过另一个极点积分的形式
在x=-a。利用均匀性，我们可以“折叠积分”，
因此，我们不需要处理两个极点，而只需要处理一个极点x=a.
的真实部分
np.exp(-1j*x) / (x**2 - a**2) = (np.cos(x) - 1j * np.sin(x)) / (x**2 - a**2)
是一个偶函数x所以整合实部x = -infinity to
infinity等于积分的两倍x = 0 to infinity。这
被积函数的虚部是奇函数x。积分由x = -infinity to infinity等于积分x = -infinity to 0，加上
积分来自x = 0 to infinity。这两部分相互抵消
因为（虚数）被积函数是奇数。所以虚部的积分等于0。
最后，使用ev-br的建议, since
1 / (x**2 - a**2) = 1 / ((x - a)(x + a))
using weight='cauchy', wvar=a隐式对被积函数进行加权1 / (x - a)因此我们可以将显式被积函数简化为
np.cos(x) / (x + a)
由于被积函数是偶函数a，我们可以不失一般性地假设a是积极的：
a = abs(a)
现在整合来自x = 0 to 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()
