Principal Value Integration with Python

I am trying to integrate a function with python. It has a singularity and using scipy.integrate I am not able to get satisfying results (I also get an IntegrationWarning).

This is my code thus far. This is my first time handling such functions and I am a bit lost.

from scipy.integrate import quad
import numpy as np

def integrand(x):
    return -(x * np.sin(2*x)) / (x * np.cos(x) + np.sin(x))

def cauchy_principal_value(f, a, b, singular_point):
    integral = quad(f, a, b, points=[a, singular_point, b],limit=1000)
    return integral

result = cauchy_principal_value(integrand, 0, np.pi, 2.0287665755744362)

print("Principal Value: ", result[0])

Can anyone help me?

In general, integrator routines need a bit of help when dealing with singularities.

Here if you know you want a principal value, tell it to quad by using the weight="cauchy" (look up the docs for quad). Here's a quick-and-dirty demo:

In [32]: x0 = brentq(lambda x: x * np.cos(x) + np.sin(x), 1e-2, np.pi)

In [33]: quad(lambda x: -(x * np.sin(2*x)) / (x * np.cos(x) + np.sin(x) * (x - x0)), 1e-4, np.pi, weight='cauchy', wvar=x0)
Out[33]: (-4.14269240855162, 1.4398780224416236e-08)

• The integrator needs a position of the singularity, so find it with a root-finder (brentq below)
• with weight="cauchy" the integrand it interpreted as f(x) in f(x) / (x - c). Therefore, work out the expansion and redefine your integrand around x0. Below I multiply by x - x0 instead.
• there is a 0/0 at x->0, so ideally redefine your integrand function at zero. Below I just step away out of x=0 instead.

The scipy.integrate.quad function is a general-purpose integrator and may encounter difficulties when handling singularities or highly oscillatory functions. In your case, you are trying to integrate a function with a singularity, which can lead to IntegrationWarning or inaccurate results.

To handle the singularity in your integrand, you can use the Cauchy principal value method. Here's an updated version of your code:

from scipy.integrate import quad
import numpy as np

def integrand(x):
    return -(x * np.sin(2*x)) / (x * np.cos(x) + np.sin(x))

def cauchy_principal_value(f, a, b, singular_point, **kwargs):
    def wrapped(x):
        if x == singular_point:
            return np.nan
        return f(x)
    
    integral = quad(wrapped, a, b, **kwargs)
    return integral

result = cauchy_principal_value(integrand, 0, np.pi, 2.0287665755744362)

print("Principal Value:", result[0])

In the updated code, a wrapper function wrapped is introduced to handle the singularity. The wrapper function returns np.nan for the singular point to indicate that the integrand is not defined at that point. The **kwargs parameter allows you to pass additional arguments to the quad function if needed.

By using this approach, the integration should proceed without raising an IntegrationWarning. Let me know, if it works or not.