Numerical Principal Value integration

Question

I want to calculate the Principal Value Integral of

$$ \int_a^\infty dx \frac{f(x)}{(x-b)^n}, $$

where $b$ is in the integration region. For the case $ n=1 $ I do know the trick

$$ \int_a^\infty dx\frac{f(x)}{x-b} = \int_a^\infty dx \frac{f(x)-f(b)}{x-b} + \int_a^\infty dx \frac{f(b)}{x-b}, $$ where after the trick the first part is numerically stable and the second part can be integrated analytically.

How can I extend this trick for higher $n$?

Answer 1

It might be difficult to come up with a solution for a general $f(x)$ and $n$ without implementing a QAWC adaptive integration for Cauchy principal values, using $$ I =\int_a^b dx \frac{f(x)}{x-c} = \lim_{\epsilon\rightarrow0}\bigg(\int_a^{c-\epsilon} dx \frac{f(x)}{x-c}+\int_{c+\epsilon}^{b} dx \frac{f(x)}{x-c} \bigg) $$

You can also experiment with the ‘cauchy’ weight for the quadrature. As an example consider the following for $f(x)=\cos(x)$. I used this example since I know the analytical solution.

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

def complex_quadrature(func, a, b, **kwargs):
    def real_func(x):
        return np.real(func(x))
    def imag_func(x):
        return np.imag(func(x))
    real_integral = quad(real_func, a, b, **kwargs)
    imag_integral = quad(imag_func, a, b, **kwargs)
    return real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]

def I(a,n):
    a = abs(a)
    return 2*(complex_quadrature(lambda x: np.cos(x)/(x + a)**n, 0, 1e6, weight='cauchy', wvar=a)[0])

I_vec = np.vectorize(I)

fig, ax = plt.subplots()

a = np.linspace(-15, 15, 150)
ax.plot(a, -np.pi*np.sin(a)/(a), label="analytical result", lw=2)
ax.plot(a, (I_vec(a,1)), "--", label="numerical result (cauchy)")
ax.set_ylim(-5, 5)
ax.set_ylabel("f(x)")
ax.set_xlabel("x")
ax.legend()
plt.title(r'$\int_{-\infty}^\infty \frac{\cos(x)}{x+a} \, dx$')
plt.show()

In principle this should work for any $n$ but you have to tweak the numerical parameters.

plot