Singularity Integration

LinearAlgebra

---

from numpy import *
from matplotlib.pyplot import *
from numpy.linalg import *
from scipy.integrate import *

Here is a simple example illustrating how to do numerical integration for singular functions.

We will consider

$$\int_{0}^1\frac{1}{\sqrt{\sin x}}dx$$

A naive idea is to do

$$\int_0^1\frac{1}{\sqrt{\sin x}}\approx \int_\varepsilon^1\frac{1}{\sqrt{\sin x}}$$

where $\varepsilon$ is a small number.

def f1(x):
    if isinstance(x,ndarray):
        y = zeros_like(x)
        for i in range(x.shape[0]):
            y[i] = 1./sqrt(sin(x[i]))
    else:
        y =  1./sqrt(sin(x))
    return y

We use simpson rule to do numerical integration.

$$\int_{x_i}^{x_{i+1}} f(x) dx \approx \frac{x_{i+1}-x_i}{6}\left[f(x_i)+4f\left(\frac{x_i+x_{i+1}}{2} \right)+f(x_{i+1}) \right]$$

def simpson_quad(func,a,b,tol=1e-6):
    n = int((b-a)/tol)
    s = 0.0
    x = a+tol*arange(n+1)
    f = func(x)
    m = a+tol*(0.5+arange(n))
    fm = func(m)
    for i in range(n):
        s +=tol/6.0*(f[i]+4*fm[i]+f[i+1])
    return s

For example, we will have

$$\int_0^1 x dx = 0.5$$

simpson_quad(lambda x:x, 0, 1)

0.49999999999999994

Now we compute the singular integration.

print(simpson_quad(f1, 1e-5,1,1e-6))
print(simpson_quad(f1, 1e-6,1,1e-6))
print(simpson_quad(f1, 1e-7,1,1e-6))

2.02848076409
2.0328057929
2.03425369423

We see even if we use steps $10^{-6}$ and only disgard $10^{-5}\sim 10^{-7}$ we may have error up to $0.002$. It is huge indeed.

We will now split

$$\int_0^1\frac{1}{\sqrt{\sin x}}dx = \int_0^1\left(\frac{1}{\sqrt{\sin x}}-\frac{1}{\sqrt{x}}\right) dx + \int_0^1\frac{1}{\sqrt{x}}dx$$ .

The latter can be computed as

$$\int_0^1\frac{1}{\sqrt{x}}dx=2$$

and for the former, we compute with simpson's rule.

def f2(x):
    y = zeros_like(x)
    for i in range(x.shape[0]):
        t = sqrt(sin(x[i]))
        s = sqrt(x[i])
        if allclose(t,0) or allclose(s,0):
            y[i] = 0.0
        else:
            y[i] = 1./t-1./s
    return y

simpson_quad(f2,0,1,1e-4)+2

2.0348053192075901

simpson_quad(f2,0,1,1e-5)+2

2.0348044178614586

simpson_quad(f2,0,1,1e-6)+2

2.0348053192075706

We see it converges rather quickly here.