2

为了规避柯西原理值，我尝试将使用小位移 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()
4

2 回答 2

3

x=a主要问题是被积函数在和 处都有极点x=-a。ev-br的帖子展示了如何处理电线杆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)

是一个偶函数，因此积分 from to x的实部将等于积分 from to的两倍。被积函数的虚部是 的奇函数。从 到 的积分等于从到的积分加上从到的积分。这两个部分相互抵消，因为（假想的）被积函数是奇数。所以虚部的积分等于0。x = -infinityinfinityx = 0infinityxx = -infinityinfinityx = -infinity0x = 0infinity

最后，使用ev-br的建议，因为

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)

现在积分 fromx = 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()

在此处输入图像描述

于 2019-08-04T11:26:47.800 回答
-2

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

于 2019-08-02T13:00:36.497 回答