Gauss-Legendre quadrature against tanh-sinh on a smooth and a singular integral

A quadrature rule replaces an integral by a weighted sum of the integrand's values: ${\int }_{-1}^{1}f(x)dx\approx {\sum }_i{w}_{i}f(x_i)$. With $n$ points there are $2n$ numbers to choose, the nodes $x_i$ and the weights $w_i$. The Gauss-Legendre rule chooses them so that the sum is exact for every polynomial of degree up to $2n-1$, which is as high as $2n$ free numbers can reach. The nodes turn out to be the roots of the Legendre polynomial $P_n$.

The cells below compute the nodes and weights, check the degree claim, and then try the rule on two integrands: $e^x$, which is smooth, and $1/\sqrt{1-x^2}$, whose integral is $\pi$ although the integrand is infinite at $x=\pm 1$. On the second one a different rule, tanh-sinh, does far better.

The reference is posts/refs/gauss-legendre-quadrature-and-where-tanh-sinh-wins.py (output beside it): numpy's leggauss, which gets the nodes from the eigenvalues of a tridiagonal matrix rather than by Newton's method as here, closed-form integrals, and the tanh-sinh rule written in numpy.

Exact through degree 9

The cell integrates $(1+x{)}^k$ for $k=0,1,\dots ,13$, whose exact integral is $2^{k+1}/(k+1)$, with the 5-point rule and with the 3-point rule. (A power of $x$ alone would be a poor test: for odd $k$ the symmetric nodes give 0, which is right, for any rule.) The 5-point rule is exact to rounding through degree 9 and then misses. The 3-point rule misses from degree 6 on, as $2n-1=5$ says.

At degree $2n$ the miss is known exactly. The error of the $n$-point rule is

for some $\xi$ in the interval, and for $f=(1+x{)}^{10}$ the tenth derivative is the constant $10!$, so with $n=5$ the error is $2^{11}(5!{)}^4/(11(10!{)}^2)=2.9318\times 10^{-3}$.

In the figure an error of exactly 0 is drawn at $10^{-17}$, since a log axis has no 0.

A smooth integrand and a singular one

$e^x$ is not a polynomial, but it is smooth everywhere, including off the real line, and polynomials approximate it well. The Gauss-Legendre error falls faster than any power of $n$: $8.2\times 10^{-10}$ at 5 points, rounding at 7.

$1/\sqrt{1-x^2}$ is a different case. Its integral is $\pi$, finite, but the function goes to infinity at both ends, and no polynomial follows it there. The rule never evaluates at $\pm 1$, which is good, since the value there is infinite; but its outermost node sits a distance of order $1/n^2$ from the end, and the integral over that last stretch, $\int dx/\sqrt{2(1-x)}$ over a length of order $1/n^2$, is of order $1/n$. The error behaves like that: it halves each time $n$ doubles, and $n$ times the error settles near 1.73. At that rate each further digit costs ten times as many points.

Tanh-sinh

The tanh-sinh rule, from Takahasi and Mori (1974), changes variables first:

As $t$ runs over the whole line, $x$ runs over $(-1,1)$, and the endpoints move off to $t=\pm \infty$. The factor $dx/dt$ falls like $\exp (-\frac{\pi }{2}{e}^{|t|})$, a double exponential, which is faster than $1/\sqrt{1-x^2}$ grows. So the new integrand is smooth and dies off quickly in both directions, and for such a function the plain trapezoid rule, equal steps $h$ in $t$, converges very fast. Here the nodes are $t=kh$ with $|t|\le 4.5$. The last term kept on the singular integrand is $2.1\times 10^{-17}$ at $h=1$, whose last node is $t=4$, and under $10^{-28}$ for every smaller step; the terms left out are smaller still.

There is one trap. Once $\frac{\pi }{2}\sinh t$ passes about 19, $\tanh$ rounds to exactly 1, and $1-x^2$ computed from $x$ is 0: the singular integrand becomes $1/0$. With $h=1/4$ that happens at 12 of the 37 nodes, and the sum is infinite. The fix is to compute $1-x^2$ directly as ${\mathrm{sech}}^2(\frac{\pi }{2}\sinh t)$, which is tiny there but not 0, and to write the integrand in terms of it.

Halving $h$ doubles the points. On the singular integral the error goes $1.9\times 10^{-3}$, $2.0\times 10^{-8}$ and about $10^{-15}$, which is rounding, at 9, 19 and 37 points: the number of correct digits more than doubles with each halving, until rounding stops it. Gauss-Legendre with 128 points was at $1.4\times 10^{-2}$. On $e^x$ the roles reverse: tanh-sinh needs 73 points for rounding where Gauss-Legendre needed 7. The substitution spends most of its points near the ends, and a smooth integrand has nothing there that needs them.