Использование scipy.quad с трюком iε: плохие результаты

Основная проблема заключается в том, что подынтегральное выражение имеет полюсы как в 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)

является четной функцией от x, поэтому интегрирование действительной части от x = -infinity до infinity будет равняться удвоенному интегралу от x = 0 до infinity. Мнимая часть подынтегрального выражения является нечетной функцией x. Интеграл от x = -infinity до infinity равен интегралу от x = -infinity до 0 плюс интеграл от x = 0 до infinity. Эти две части компенсируют друг друга, поскольку (мнимое) подынтегральное выражение нечетное. Таким образом, интеграл от мнимой части равен 0.

Наконец, используя предложение 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)

Теперь интеграция от 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()

Вместо этого вы можете использовать аргумент weight = "cauchy" для quad. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html