# Numerical Principal Value integration

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

• Can you maybe use Integration-By-Parts? For example for $n=2$ we have $\frac{1}{(x-b)^2} = -\frac{d}{dx} \frac{1}{x-b}$ and using IBP you end up with the $n=1$ computation just with $f'(x)$ instead. Commented Aug 30, 2023 at 12:27
• I'm not sure the PV exists for even $n$, unless $f$ has a suitable zero at $b$
– Sal
Commented Aug 30, 2023 at 13:11

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
import matplotlib.pyplot as plt

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