Tanh-Sinh quadrature

Tanh-sinh quadrature is one of a sequence of methods known as “double exponential” methods, because the weights used die away at the rate of the double exponential function

    \[ e^{-e^x}. \]

In order to integrate

    \[ I=\int^1_{-1}f(x)\,dx \]

we use the substitution

    \[ x=\tanh\left(\frac{\pi}{2}\sinh t\right) \]

which transforms the integral into

    \[ I=\int^\infty_{-\infty}f\left(\tanh\left(\frac{\pi}{2}\sinh t\right)\right)\frac{\frac{\pi}{2}\cosh t}{\cosh^2\left(\frac{\pi}{2}\sinh t\right)}\,dt. \]

So far this seems like we’re going in the wrong direction – ordinarily we’d want to transform an infinite integral into a nice finite one for computation.

But consider the functions

    \begin{align*} \psi(t)&=\tanh\left(\frac{\pi}{2}\sinh t\right)\\ \psi'(t)&=\frac{\frac{\pi}{2}\cosh t}{\cosh^2\left(\frac{\pi}{2}\sinh t\right)} \end{align*}

At even so small a value as t=6, we have \psi(t)\approx 1-1.225\times 10^{-273} and \psi'(t)\approx 7.767\times 10^{-271}. Suppose we approximate the transformed integral by the trapezoidal rule as

    \[ I=\sum_{k=-N}^Nf(\psi(kh))\psi'(kh) \]

for a small enough h and large enough N. Here the values \psi(kh) are the nodes or abscissae, and \psi'(kh) are the weights.

We only need N large enough so that the die-off we have just seen happens. For example, if h=1/4, then for k such that |kh|>6 – so that |k|>24 we are dealing with weights of the order of 10^{-271}.

The tanh-sinh rule thus approximates an integral of a function f(x) over [-1,1] by choosing a small h and a reasonable N, and computing the above sum. The remarkable thing is how accurate this rule is for relatively large values of h and small values of N. (This can be shown using the Euler-Maclaurin summation formula, and I’ll discuss this is a later post).

The other very nice aspect of this rule is that if we are using real arithmetic to a high enough precision, then we can deal happily with functions with singularities at one or both endpoints, as we never actually calculate the function at those values. For t=6 as above, the closest we get to \pm 1 is about 1.225\times 10^{-273}.

So let’s look at two examples, the first one is

    \[ \int^1_{-1}\sqrt{1-x^2}\dx=\pi/2. \]

We’ll do this in PanAxiom, starting by setting a high precision and output
(1) -> digits(1000);outputGeneral(200);
and entering the functions \psi(x) and \psi'(x):
(2) -> psi:Float->Float
(3) -> psi(x)==tanh(%pi/2*sinh(x))
(4) -> psid:Float->Float
(5) -> psid(x)==%pi/2*cosh(x)/(cosh(%pi/2*sinh(x)))^2
Now we’ll try a few values of h and associated vaues of N:
(6) -> m:=3
(7) -> h:=0.5^m
(8) -> N:=6*2^m-1
(9) -> xs := [psi(h*t) for t in -N..N]
(10) -> ws := [psid(h*t) for t in -N..N]
(11) -> int3 := reduce(+,map((w,x)+->h*w*f(x),ws,xs)

(12) -> abs(int3-%pi/2)
5_5987084261_8034013826_3090769499_7396765187_2670125846_199135277 E -27

So even with this fairly large value of h=1/8 and only a few function calls (93) we have accuracy to 27 decimal places. Clearly we need to try with another m:
(13) -> m:=5
(14) -> h:=0.5^m;N:=6*2^m-1
(15) -> xs := [psi(h*t) for t in -N..N]
(16) -> ws := [psid(h*t) for t in -N..N]
(17) -> int5 := reduce(+,map((w,x)+->h*w*f(x),ws,xs)

(18) -> abs(int5-%pi/2)
8_3521709505_2582597229_5158482247_8277177970_8098898047_2412558066 E -127

and we have 127 decimal place accuracy! In general, working to high enough precision, we should expect an approximate doubling of the number of correct decimal places for every halving of h.

For a second example, consider

    \[ -\frac{1}{2}\int^1_{-1}\log\left\log\left(\frac{2}{x+1}\right)\right)\,dx=\gamma \]

where \gamma is the Euler-Mascheroni constant. This integral has a singularity at x=-1. Here goes:
Put m=6 and compute h,N and the two lists of nodes and weights as above. Then
(23) -> int6:=reduce(+,map((w,x)+->h*w*f(x),ws,xs)

and this can be checked as being correct to all displayed 200 decimal places.

Tanh-sinh quadrature is an impressive method, but seems little used. I’m guessing that’s because it’s quite new (1974) – well, new in comparison to other methods, and also hard to perform a solid error analysis. However, David Bailey, who is one of the champions of this method – see a paper of his here – has done several decades of computations and experiments, and claims that “Overall, the tanh-sinh scheme appears to be the best for integrands of the type most often encountered in experimental math research.”

You can read the original 1974 paper of Hidetosi Takahasi and Masatake Mori here.

One thought on “Tanh-Sinh quadrature

Leave a Reply

Your email address will not be published. Required fields are marked *