Python: Find principal value of an integral numerically

Question

I'm solving the integral numerically using python:

where a(x) can take on any value; positive, negative, inside or outside the the [-1;1] and eta is an infinitesimal positive quantity. There is a second outer integral of which changes the value of a(x)

I'm trying to solve this using the Sokhotski–Plemelj theorem:

However this involves determining the principle value, which I can't find any method to in python. I know it's implemented in Matlab, but does anyone know of either a library or some other way of the determining the principal value in python (if a principle value exists)?

Answer 1

You can use sympy to evaluate the integral directly. Its real part with eta->0 is the principal value:

from sympy import *
x, y, eta = symbols('x y eta', real=True)
re(integrate(1/(x - y + I*eta), (x, -1, 1))).simplify().subs({eta: 0})
# -> log(Abs(-y + 1)/Abs(y + 1))

Matlab's symbolic toolbox int gives you the same result, of course (I'm not aware of other relevant tools in Matlab for this --- please specify if you know a specific one).

You asked about numerical computation of a principal value. The answer there is that if you only have a function f(y) whose analytical form or behavior you don't know, it's in general impossible to compute them numerically. You need to know things such as where the poles of the integrand are and what order they are.

If you on the other hand know your integral is of the form f(y) / (y - y_0), scipy.integrate.quad can compute the principal value for you, for example:

import numpy as np
from scipy import integrate, special

# P \int_{-1}^1 dx 1/(x - wvar) * (1 + sin(x))
print(integrate.quad(lambda x: 1 + np.sin(x), -1, 1, weight='cauchy', wvar=0))
# -> (1.8921661407343657, 2.426947531830592e-13)

# Check against known result
print(2*special.sici(1)[0])
# -> 1.89216614073

See here for details.