# Numerical Integration for integrable singularity

Till this time i have learned three numerical technique to find the definite integration. They are Simpson, Trapezoidal and Gauss-legendre formula. The sad thing is that I can't apply these theorem directly of my integration has any integrable singularity within the interval.

Can you give me any special technique so i can use these theorem for that type of integrations?

With an open quadrature method such as Gauss-Legendre you may not need to evaluate the integrand at the point of singularity. However, proceeding this way will most likely result in a severe loss of accuracy,

If you are lucky, you may be able to remove the singularity with a change of variables. For example if $f \in C([0,1])$ and $0 < \alpha < 1/2$ then the improper integral

$$\int_0^1 \frac{f(x)}{x^\alpha} \, dx,$$ can be transformed under the change of variables $t = x^\alpha$ to

$$\frac{1}{\alpha}\int_0^1 f(t^{1/\alpha})t^{(1 - 2\alpha)/\alpha}\, dt,$$

which is a proper integral and can be handled efficiently by the techniques you mentioned.

More generally, an understanding of the asymptotic behavior of the integrand near the singularity is important.

• What if the singularity is of type $\ln x$? You'd have to substitute $t=\ln x$, which would require to calculate an integral on infinite interval, which seems not much better. – Ruslan Mar 15 '17 at 7:09
• @Ruslan: This is just an example of a method that can handle one type of singularity. Note the word "lucky" above. In the case of a logarithmic singularity there are special forms of Gaussian quadrature that can be applied. – RRL Mar 15 '17 at 7:57
• @Ruslan Being specific, $\int_0^1 f(\ln(x)) dx = \int_{-\infty}^0 e^u f(u) du$ which is now conducive to Gauss-Laguerre quadrature. – Ian Mar 15 '17 at 10:13

The standard numerical methods are designed for smooth (or at least sufficiently differentiable) functions, and don't perform very well in the presence of a singularity. One thing that may work well is to subtract out the singularity. That is, if you want $\int_a^b f(x)\; dx$, try to write $f(x) = g(x) + h(x)$ where $g(x)$ has the singularity but is integrable in closed form, while $h(x)$ has no singularity. See e.g. these notes of mine.

• Thank you! Just used this advice for quite a complicated integral and it worked very well, as it just had an $1/\sqrt{x}$ singularity to be removed – Yuriy S Nov 4 '18 at 11:33

There are a few different methods to deal with this problem. Gaussian quadrature can work in some cases, you then absorb the singularity in a suitably chosen weight function and then you proceed in the usual way with constructing the orthogonal polynomials where the inner product is defined using that weight function.

Another method is to do a change of variables to eliminate the singularity as mentioned in RRL's answer, or you can split the integrand in a singular part and a non-singular part as mentioned in Robert Israel's answer.

A general purpose method that doesn't require much effort to implement is the Tanh-sinh quadrature method, here you transform an integral from minus 1 to 1 to one over the entire real line, singularities at the endpoints are then not going to cause problems.

Ideally, you would know something about the integrand in a neighborhood of the singularity which would allow you to resolve its impact. For example, suppose you want $\int_0^1 f(x) dx$, but $f$ has a singularity at $x=0$, behaving as $x^{-1/2}+O(1)$ there, and it has no other singularities. Then you can write the integral as $\int_0^\delta x^{-1/2} dx + \int_0^\delta f(x)-x^{-1/2} dx + \int_\delta^1 f(x) dx$ where $\delta>0$ is a tuning parameter. Then the first term is exactly $2\delta^{1/2}$ and the other terms have no singularities. You can make a similar "regularizing" transformation using a change of variables as RRL suggested.

You could also use adaptive quadrature in this situation. Essentially, if the singularity is of order $-1<p<0$, then any reasonable quadrature method on that interval will be of order $p+1$, so that Richardson extrapolation can be applied, provided you know exactly what $p$ is. Note that this requires a much less detailed asymptotic than the first suggestion from the previous paragraph; like the second suggestion from that paragraph, it only requires you to know the order of the singularity. I can go into specifics about this if you like, feel free to comment.

If you are repeatedly dealing with the same "type" of singularity behavior, such as a singularity of order $-1/2$ at both endpoints, then it may be worth developing a Gaussian quadrature method (or you may be in a case where a standard one applies). This simplifies the problem of handling the singularity because it allows you to deal with much simpler integrands, but it does still require you to be able to compute some singular integrals in order to compute the nodes and weights.