Integration of Singularity within integral

Question

(Kramer Kroenig EQ Huang, et al.,)

Hello, does anyone know how to integrate the function below to solve for the Phase Amp using Kramers Kroenig Relation's in python. The Integral can be definite through the range of 1.2 - 8 however I am unsure on how to properly integrate the singularity point using the quad function. I have tried specifying singularity points with 'points' and weights like with 'cauchys' I am unsure if I am doing it correctly as the integral goes wild. R(E') takes values from a text file with values for the points of the integral.

def reflection_function(energy):
    return np.interp(energy, ref_ener, ref_ref)

def integrand(e_prime, energy):
    return np.log(reflection_function(e_prime))/(e_prime**2 - energy**2)

integral_result = quad(integrand, 1, 8, args=(energy,), weight = 'cauchy', wvar = energy)[0]

print(integral_result)

It returns wild values like -220310228.0672309 and it should be around -0.7655680665661697. I am unsure why this is happening. Also I get this error

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  integral_result = quad(integrand, 1, 8, args=(energy,), weight =

Answer 1

Putting weight="cauchy", wvar=energy will automatically apply a weighting 1/(eprime-energy) associated with this pole, so this doesn't need to be included in the integrand.

Thus, you can keep your statement using quad, but you need to send it a modified integrand:

def integrand(e_prime, energy):
    return np.log(reflection_function(e_prime))/(e_prime + energy)     # the rest

integral_result = quad(integrand, 1, 8, args=(energy,), weight = 'cauchy', wvar = energy)[0]

(This does assume that your reflection_function is positive, but since we can't see it we can't do anything about that.)

I'm also assuming that your energy is positive - use the other pole if not.

As a simpler example, the following program would find the principal value of the integral of 1/(x²-1) from 0 to 3. The integrand has a simple pole at x=1, so is modified in the call to quad by removing the factor (x-1) from the denominator. The principal value at the end should be (1/2)ln(1/2), or -0.34657...

from scipy.integrate import quad

def integrand_modified( x ):
    return 1.0 / ( x + 1.0 )

i = quad(integrand_modified, 0.0, 3.0, weight='cauchy', wvar = 1.0 )[0]
print( i )