# Is there a way to deal with this singularity in numerical integration?

I would like to compute numerically, e.g., using the standard method of trapezes the following definite integral over the interface $$[0,1]$$: $$I = \int_0^1 \frac{f(x)}{\sqrt{1-x^2}} \, \mathrm{d} x \, ,$$
where $$f(x)$$ is continuous in the interval [0,1]. It can readily be shown analytically that the integral is well convergent. However, when proceeding numerically, difficulties arise since the integrand diverges at the upper limit of integration for $$x=1$$.

I was wondering whether there exists a procedure that can help to remove the singularity in this integral. Any help is highly appreciated

Thank you,

hartmut

• $I=\int^{\pi/2}_0f(\sin\theta)\,d\theta$ might help. – user436658 Sep 11 '17 at 10:24
• as soon as $f$ is also diffentiable, integration by parts looks feasible – tired Sep 11 '17 at 12:52
• There is a lot about this sort of question in Forman Acton's nice little book "Real computing made real". – Gareth McCaughan Sep 12 '17 at 10:17
• @GarethMcCaughan Errors of the third kind indeed! I will get my hands on a copy; wish I had done so 20 years ago! – uhoh Sep 12 '17 at 12:57

It's easy to convert your integral into one without any singularities via the substitution $x=\sin\theta$: $$I=\int^{\pi/2}_0f(\sin\theta)\,d\theta.$$ If you're lucky to have a function $f$ that is even (so that $I$ is just a quarter of the integral over the full period, from $-\pi$ to $\pi$), you might be pleasantly surprised by the speed of convergence of the trapezoidal rule in this special situation.

• ...and if indeed certain conditions apply, and you do find the trapezoidal rule to be performing exceptionally well, see this for an explanation. – J. M. isn't a mathematician Sep 11 '17 at 12:15
• The numerical integration indeed performs powerfully well. Only a few number of discretization points are needed to compute the integral – Daddy Sep 11 '17 at 12:24
• Following up on your great answer, $f$ is in fact even, meaning that the integral above = 1/4 of that over a full period. But why does that speeds the convergence up? It is easier to perform a numerical integration over [0, $\pi/2$] than to integrate numerically over [0, $\pi$] and then take the quarter.. your clarification is highly appreciated. Thanks! – Daddy Mar 23 '18 at 14:25

A general technique is to remove the singularity by moving it to an analytically integrable function:

$$\int_0^1\frac{f(x)}{\sqrt{1-x^2}}dx=\int_0^1\frac{f(x)-f(1)+f(1)}{\sqrt{1-x^2}}dx=-\int_0^1\frac{f(x)-f(1)}{\sqrt{1-x^2}}dx+f(1)\left.\arcsin(x)\right|_0^1.$$

The limit of the integrand at $x=1$ should be finite.

For instance with $f(x):=x$,

$$\frac{x-1}{\sqrt{1-x^2}}=-\sqrt{\frac{1-x}{1+x}},$$ which is well-behaved.

• Nice trick, for sure ! – Claude Leibovici Sep 11 '17 at 10:52
• It should be $$\frac{x-1}{\sqrt{1-x^2}}=-\sqrt{\frac{1-x}{1+x}},$$ and "well-behaved" isn't entirely true: the integrand isn't differentiable for $x=1$, so the composite trapezoidal rule with a partition into $N$ equal subintervals has an error of about $\frac1{2N}$ in this case, while it would be $O(N^{-2})$ for a regular ($f''$ bounded) integrand. – user436658 Sep 11 '17 at 12:21
• @ProfessorVector: that's right, you have a vertical tangent, which is not ideal for the Newton-Cotes method. My point was that the unbounded behavior is removed, and the method allows to remove several singuarities, of different types. – Yves Daoust Sep 11 '17 at 13:58
• The $\arcsin(x)$ outside the integral lacks context. I think you mean $\frac\pi2$. (+1) anyway, because this was the answer I was thinking of, but overthought in my answer. – robjohn Sep 12 '17 at 15:54
• @robjohn: you are right, I fixed. – Yves Daoust Sep 12 '17 at 16:29

The other answers so far focus on how to transform the integral so that you don't have a singularity. But that's not always viable so it's important to know how to deal with a singularity when one exists.

The trapezoid rule is a poor choice for this case, because it uses the values of the function at the endpoints, which can be infinite. Simposon's rule will have a similar problem.

More appropriate is the more primitive rectangle rule. That is, divide the domain to segments, and multiply each segment's width by the value of the function at the midpoint of the segment. This will work since you never try to sample the function at the infinite endpoints.

Alternatively, if you want faster convergence, you can use the more sophisticated and extremely powerful Gaussian Quadrature. More complicated but well worth knowing. Here, too, you divide the domain to intervals. For each interval, you sample the function at a few carefully chosen points to obtain surprisingly high accuracy.

There are different methods depending on the number of points in each interval, but the simplest is the 2-point method. And again, it will work with singularities, since the endpoints of intervals aren't sampled.

There is even a variant of Gaussian quadrature tailored specifically for a well-behaved function multiplied by $\frac{1}{\sqrt{1-x^2}}$, called Chebyshev–Gauss quadrature. To use it, you don't even have to divide to intervals - you just choose $n$ (the higher $n$, the better the accuracy) and you have

$$\int_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}\approx\frac{\pi}{n}\sum_{i=1}^nf\left(\cos\left(\frac{2i-1}{2n}\pi\right)\right)$$

Which is very similar to transforming the integral as some other answers suggested. It will maximize the accuracy you can have for a given number of samples from the function $f$, and convergence is exponential.

• For reference: the last formula in this answer is Gauss-Chebyshev quadrature (since the underlying orthogonal polynomials are the Chebyshev polynomials.) – J. M. isn't a mathematician Sep 11 '17 at 14:55
• @J.M.isnotamathematician - Yep, I intended to mention this but apparently forgot to. I've added the name now. – Meni Rosenfeld Sep 11 '17 at 15:02
• FTR, Gauß quadrature per se doesn't buy you anything over simpler methods if there are singularities in the domain. Its only benefit is excellent accuracy for integrating on a compact interval $I$ a function that's smooth on $I$ and in some environment beyond. – leftaroundabout Sep 12 '17 at 14:11

Since your integral is over a subset of the interval $(-1,1)$, your $f$ is continuous on that interval, and your integrand only has singularities on the boundary of that interval, tanh-sinh quadrature is a candidate method for you.

This will change your interval of integration to $[0,\infty]$, but the integrand will decay doubly-exponentially, (e.g., like $\mathrm{e}^{-\mathrm{e}^n}$), so truncation at finite right endpoint need not be significant.

You can continue your plan of using the trapezoidal method with step size $h$ by using the sample points $x_k = \tanh(\frac{\pi}{2} \sinh kh)$ and weights $w_k = \frac{h \pi}{2} \cdot \frac{\cosh kh}{\cosh^2( \frac{\pi}{2} \sinh kh)}$. For instance, the area of your first trapezoid would be $$A_0 = \frac{x_1 - x_0}{2} \left(\frac{w_0 f(x_0)}{\sqrt{1 - x_0^2}} + \frac{w_1 f(x_1)}{\sqrt{1 - x_1^2}} \right) \text{.}$$ Just using $x_k$ and $w_k$ in a left Riemann sum, (and halving $w_0$ when summing this way) convergence is exponential: halving $h$ doubles the number of correct digits. Error analysis suggests the trapezoidal rule applied here should improve convergence slightly (although perhaps is not worth the "effort" compared to just computing $\sum_{k=0}^\infty w_k \frac{f(x_k)}{\sqrt{1-x_k^2}}$, stopping when the terms are small enough to meet your precision and accuracy requirements or when either $1-x_k$ or $w_k$ is so small you cannot easily represent it).

Example: $f(x) = 1$ :

The value of the integral is $\frac{\pi}{2} = 1.57079\dots$. Using left-Riemann summation (and halving $w_0$), with $h = 1/10$, and stopping when $k=20$ (because $1-x_k$ and $w_k$ are about $10^{-5}$), the sum is $1.5659 \dots$. Continuing on to $k = 27$ obtains $1.57079\dots$.

• With $h = 1/20$, stopping at $k = 40$ ($1 - x^k$ around $10^{-5}$ and $w_k$ around $10^{-36}$), obtains $1.56504\dots$.
• With $h = 1/20$, stopping at $k = 54$ ($1 - x^k < 10^{-6}$ and $w_k$ around $10^{-150}$), obtains $1.57078\dots$. (An arbitrary precision calculation gives an error of about $5 \times 10^{-12}$ when we sum $70$ terms.)

We're only summing tens of terms to get these results...

• The double-exponential quadrature of Takahasi and Mori is indeed a good way to flatten singularities at either end of the integration interval, but one needs to be careful with evaluations near the endpoints, as pointed out in their paper. – J. M. isn't a mathematician Sep 11 '17 at 14:56
• @J.M.isnotamathematician : I agree. I don't think I've hidden that. I repeatedly comment on the difficulty of representing $1-x_k$ and $w_k$ for large $k$... – Eric Towers Sep 11 '17 at 15:19

I like Professor Vector's method, but here is a way to take advantage of the integrability of $\frac1{\sqrt{1-x^2}}$ on $(0,1)$.

Using the integral $$\int_0^1\frac{a+bx}{\sqrt{1-x^2}}\,\mathrm{d}x=\frac{\pi a}2+b$$ we can write $$\int_0^1\frac{f(x)}{\sqrt{1-x^2}}\,\mathrm{d}x=\int_0^1\frac{f(x)-(1-x)f(0)-xf(1)}{\sqrt{1-x^2}}\,\mathrm{d}x+\left(\frac\pi2-1\right)f(0)+f(1)$$ If $f$ is smooth at $0$ and $1$, the integrand on the right vanishes at $0$ and $1$

A Lesson in Overthinking

I was wondering why $f(0)$ and $f(1)$ had different weights. Then I realized that $\frac1{\sqrt{1-x^2}}$ doesn't have a singularity at $0$. We don't need to subtract $(1-x)f(0)$ at all. While the formula I give above is correct, I was thinking of $\sqrt{x(1-x)}$ in the denominator, which does have a singularity at $0$ and $1$. Doh!

The standard way is $$x = 1 - t^{2}$$: $$\int_{0}^{1}{\operatorname{f}\left(x\right) \over \,\sqrt{\,{1 - x^{2}}\,}\,}\,\mathrm{d}x = 2\int_{0}^{1}{\operatorname{f}\left(1 - t^{2}\right) \over \,\sqrt{\,{2 - t^{2}}\,}\,}\,\mathrm{d}t$$