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?

Solution 1:

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.

Solution 2:

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.

Solution 3:

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.

Solution 4:

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.

If you don't know anything quantitative about the singularity, then about all you can do is adaptive quadrature. Unfortunately the usual construction of adaptive quadrature methods using Richardson extrapolation will fall through in the presence of a singularity, because the order of the method is different in the presence of the singularity, and you need to know the exact order to know how to set up your Richardson extrapolation. So you will be stuck with essentially "brute force" adaptive quadrature instead: keep locally refining repeatedly until the difference between successive refinements is less than your tolerance. This method is not mathematically safe (in particular, it may fail to detect non-integrable singularities even when they are present), but it may work in practice anyway.

If you don't even know where the singularity is...good luck.