The working of this quadrature method can be explained in part by the Euler-Maclaurin summation formula, which elegantly relates integral and sums:

Here and are natural numbers, are the Bernoulli numbers, and is the remainder. With a double exponential function, such as , the higher derivatives approach zero very quickly. For example:

(1) -> psi:=(x:Float):Float+-> tanh(%pi/2*sinh(x)) (2) -> function(D(psi(x),x),psid,x) (3) -> function(D(psi(x),x,3),psi3d,x) (4) -> function(D(psi(x),x,5),psi5d,x) (5) -> digits(2000); outputFloating(10); (6) -> psid(7.5) (6) 0.2135882778 E -1229 (7) -> psi3d(10.0) (7) 0.7961568243 E -1987 (8) -> psi5d(10.0) (8) 0.8548668808 E -1977

This means that we can approximate the integral of a double exponential function using a trapezoidal sum with remarkable precision, as we saw in my previous post.

Also, since the points at which the function is computed are nested, tanh-sinh integration can be computed efficiently by at each stage only computing the “new” points. This is precisely analogous to the computation of the trapezoidal rules needed for Romberg integration.

The function is odd, and so the weight function is even. This means we only have to compute the positive values of the weights, thus halving the amount of work at each stage. This means that, for a function defined over :

Suppose that is the value obtained using and being all multiples of for which for some value . We can choose in advance so that all the derivatives

are small. We can choose : as we have seen above, derivatives at are not far off .

If is halved, then the new values at which we will have to compute and will be

Let be the sum

for the new value and the odd multiples of just given. Then:

We can encapsulate all of this in a simple program, which I’ve called `tsquad.input`:

p2 := %pi/2.0; psi(x:Float):Float == tanh(p2*sinh(x)); psid(x:Float):Float == p2*cosh(x)/(cosh(p2*sinh(x)))^2; tanhsinh(f:Float->Float,a:Float,b:Float,m:NNI):Float== g := (x:Float):Float +-> f((b-a)/2*x+(b+a)/2) h := 0.5; N := 13; hs := [h*t for t in 1..N] int := reduce(+,[h*psid(k)*(g(psi(k))+g(-psi(k))) for k in hs]) int := (int+h*psid(0.0)*g(0.0))*(b-a)/2 for j in 2..m repeat h := 0.5*h N := 2*N+1 hs := [h*t for t in 1..N by 2]; nn := reduce(+,[h*psid(k)*(g(psi(k))+g(-psi(k))) for k in hs]) int := 0.5*int+nn return(int)

For example, we might try to integrate the wacky function given as an example by the eminent numerical analyst Nick Trefethen of his chebfun package:

You can see a picture of the function here.

So here goes:

(9) -> )read tsquad (10) -> digits(100); outputGeneral(20); (11) -> f:=x+->x*sin(2*exp(2*sin(2*exp(2*x)))) (12) -> tanhsinh(f,-1,1,8) (12) 0.3367328347_8172753599

and this certainly agrees with the digits given on the cited page. In fact , we can go a little further and see how accurate we can get by looking at the differences between successive approximations:

(13) -> Is:=[tanhsinh(f,-1,1,k) for k in 1..12]; (14) -> for i in 2..12 repeat print(abs(Is.i-Is.(i-1))) 0.0654260836_6531745520_3 0.0914750030_3956870582_7 0.1216370035_8161735876 0.0444735581_9242514396_5 0.0003994963_7625185334_381 0.0000156965_4202311656_8617 0.9403290530_5548381954 E -15 0.5253804892_4141132667 E -37 0.1420739498_0670464881 E -91 0.1142987391_2822749822 E -99 0.3286088749_9365405739 E -99

and the fact that we are not getting better approximations is simply that we are working at the limit of our accuracy, having set the precision to 100 digits. If we increased the precision, we would get greater accuracy, at the expense of more time required to perform the computations.

At any rate, tanh-sinh quadrature has performed extremely well.