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)
1.5707963267_9489661923_1321692195_2032049740_0882640697_2503877808_704745288
2_7101748171_3942893866_2316742183_6796150502_5937316895_0417577408_475800854
4_5628142836_8673040694_1249390695_6344814583_2236941053_851601884

(12) -> abs(int3-%pi/2)
0.5554517628_7542412671_9419593390_3364085913_8006787437_7214628876_453560615
6_8437688746_1769296163_5714948354_3607228896_9582948438_7726267453_442704425
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)
1.5707963267_9489661923_1321691639_7514420985_8469968755_2910487472_296153908
2_0314310449_9314017412_6710585339_9107404325_6641153323_5469223238_784499800
4_5287150025_7861373847_7334689949_7174954259_0751407012_954089269

(18) -> abs(int5-%pi/2)
0.1910315886_4182607445_9615455786_5963129554_8527914401_4792528292_268179437
3_5963818590_2038180838_1938683423_3585421354_3187584226_8658998324_329887020
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:
f(x)==-0.5*log(log(2/(x+1)))
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)

0.5772156649_0153286060_6512090082_4024310421_5933593992_3598805767_2348848677
_2677766467_0936947063_2917467495_1463144724_9807082480_9605040144_8654283622
_4173997644_9235362535_0033374293_7337737673_9427925952_5824709492
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 *