Using Scipy.quad With Iε Trick: Bad Results

In order to circumvent the cauchy principle value, I tried to integrate an integral using a small shift iε into the complex plane to evade the pole. However, as can be inferred fr

Solution 1:

The main problem is that the integrand has poles at both x=a and x=-a. ev-br's post show how to deal with a pole at x=a. All that's needed then is to find a way to massage the integral into a form that avoids integrating through the other pole at x=-a. Taking advantage of evenness allows us to "fold the integral over", so instead of having two poles we just need to deal with one pole at x=a.

---

The real part of

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

is an even function of x so integrating the real part from x = -infinity to infinity would equal twice the integral from x = 0 to infinity. The imaginary part of the integrand is an odd function of x. The integral from x = -infinity to infinity equals the integral from x = -infinity to 0, plus the integral from x = 0 to infinity. These two parts cancel each other out since the (imaginary) integrand is odd. So the integral of the imaginary part equals 0.

Finally, using ev-br's suggestion, since

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

using weight='cauchy', wvar=a implicitly weights the integrand by 1 / (x - a) thus allowing us to reduce the explicit integrand to

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

Since the integrand is an even function of a, we can assume without loss of generality that a is positive:

a = abs(a)

Now integrating from x = 0 to infinity avoids the pole at 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()

Solution 2:

You can use instead the weight="cauchy" argument to quad. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html