11

I'm solving the integral numerically using python:

enter image description here

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: enter image description here

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

pv.
  • 33,875
  • 8
  • 55
  • 49
Rasmus Gross Søgaard
  • 291
  • 4
  • 14
  • How to implement it in MATLAB? – Ka Wa Yip Jun 22 '16 at 08:18
  • In MATLAB the symbolic integration "int" can handle principal values: http://se.mathworks.com/help/symbolic/int.html Otherwise the numerical integrator "integral" can also handle singularities at endpoints. So you can split the integral in two exactly add the singularity and then add the two results: http://se.mathworks.com/help/matlab/ref/integral.html?searchHighlight=integral – Rasmus Gross Søgaard Jun 22 '16 at 14:43

1 Answers1

8

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.

pv.
  • 33,875
  • 8
  • 55
  • 49