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$?
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.
