Subtraction method for numerical integration with singularity.

Question

In https://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities, the author explains the subtraction method to get rid of singularities when performing numerical integration.

The example he gives is the integral $\int_{0}^{1}\frac{1}{\sqrt{\sin x}}dx$.

To accurately compute this integral numerically, one first finds the approximation of integrand near $x=0$, which is $\frac{1}{\sqrt{x}}$ for this case.

Then, the original integral can be written like $\int_{0}^{1}(\frac{1}{\sqrt{\sin x}}-\frac{1}{\sqrt{x}})dx+\int_{0}^{1}\frac{1}{\sqrt{x}}dx$.

The latter term can be computed analytically by hand, so the author mentions that the work only left is computing the former one with numerical integration method such as trapezoidal, Simpson, etc.

I understand that the singularity is eliminated from the former integrand but I could not understand how one can actually compute it with numerical method.

For example, let's say I want to compute it with Simpson's method:

$\int_{0}^{2h}f(x)dx=\frac{h}{3}[f_{0}+4f_{1}+f_{2}]+\mathcal{O}(h^{5})$.

Then, I must evaluate $f(0)$ with the computer to fully compute the former integral. However, how can one compute $\frac{1}{\sqrt{sin(0)}}-\frac{1}{\sqrt{0}}$ with the computer? Shouldn't it give DivisionByZero error?

I want to understand how can this problem solved or am I misunderstanding something?

Answer 1

Obviously you may not attempt numerical evaluation of the function at the origin. But obviously also the value is $0$ so that you needn't even compute the term !

Answer 2

Obviously, numerical evaluation of the difference will replicate the problem with the singularity. One has thus to first transform the difference algebraically, $$ \frac1{\sqrt{\sin x}}-\frac1{\sqrt{x}}=\frac{\sqrt{x}-\sqrt{\sin x}}{\sqrt{x\sin x}}=\frac{x-\sin x}{x\sqrt{\sin x}+\sqrt{x}\sin x}, $$ and analytically, $$ =\frac{\frac16x^3(1-\frac1{20}x^2+\frac1{840}x^4+...)}{\sqrt{x}^3(2 - \frac14x^2 + \frac{13}{1440}x^4+...) }\approx x^{3/2}\frac{1-\frac1{20}x^2}{12-\frac32x^2}. $$

---

Actually comparing plots of the two, the difference between the original and the de-singularized approximation is only relevant for step sizes $h<10^{-6}$.

plot of difference orig minus approx

Otherwise just setting $f_0=0$ is indeed sufficient. Between $h=10^{-6}$ and the best switchover point at $h=10^{-3}$ the error in the function value $f_1=f(x_1)$ at $x_1=h$ will be dominated by the method error of the quadrature rule.