CODEWITHSUNDEEP
pythonscipy
user10289025
user10289025

Reputation:

Numerical integration with singularities in python (principal value)

I'm trying to integrate a function with singularities using the quad function in scipy.integrate but I'm not getting the desired answer. Here is the code:

from scipy.integrate import quad
import numpy as np
def fun(x):
    return 1./(1-x**2)
quad(fun, -2, 2, points=[-1, 1])

This results in IntegrationWarning and return value about 0.4.

The poles for the function are [-1,1]. The answer should be of roughly 1.09 (calculated using pen and paper).

Upvotes: 7

Views: 2968

Answers (2)

user6655984
user6655984

Reputation:

The option weight='cauchy' can be used to efficiently compute the principal value of divergent integrals like this one. It means that the function provided to quad will be implicitly multiplied by 1/(x-wvar), so adjust that function accordingly (multiply it by x-wvar where wvar is the point of singularity).

i1 = quad(lambda x: -1./(x+1), 0, 2, weight='cauchy', wvar=1)[0]
i2 = quad(lambda x: -1./(x-1), -2, 0, weight='cauchy', wvar=-1)[0]
result = i1 + i2

The result is 1.0986122886681091.

With a simple function like this, you can also do symbolic integration with SymPy:

from sympy import symbols, integrate
x = symbols('x')
f = integrate(1/(1-x**2), x)
result = (f.subs(x, 2) - f.subs(x, -2)).evalf()

Result: 1.09861228866811. Without evalf() it would be log(3).

Upvotes: 5

user8408080
user8408080

Reputation: 2478

I also couldn't get it to work with the original function. I came up with this to evaluate the principal value in scipy:

def principal_value(func, a, b, poles, eps=10**(-6)):
    #edges
    res = quad(func,a,poles[0]-eps)[0]+quad(func,poles[-1]+eps,b)[0]
    #inner part
    for i in range(len(poles)-1):
        res += quad(func, poles[i]+eps, poles[i+1]-eps)[0]
    return res

Where func is your function handle, a and b are the limits, poles is a list of poles and eps is how near you want to approach the poles. You can make eps smaller and smaller to get a better result, but maybe sympy will be better for a problem like this.

With this function and the standard eps I get 1.0986112886023367 as a result, which is almost the same as wolframalpha gives.

Upvotes: 2

Related Questions