Tanh-sinh quadrature (2)

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

    \[ \sum_{i=a+1}^b f(i) =     \int^b_a f(x)\,dx + B_1 \left(f(b) - f(a)\right) +     \sum_{k=1}^p\frac{B_{2k}}{(2k)!}\left(f^{(2k - 1)}(b) - f^{(2k - 1)}(a)\right) +      R. \]

Here a and b are natural numbers, B_k are the Bernoulli numbers, and R is the remainder. With a double exponential function, such as \tanh(\pi/2 \sinh x), 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 \psi(x) is odd, and so the weight function \psi'(x) 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 [-1,1]:

    \[ \sum^{N}_{k=-N}\psi'(kh)f(\psi(kh))=\psi'(0)f(0)+\sum^N_{k=1}\psi'(kh)\Bigl(f(\psi(kh))+f(-\psi(kh))\Bigr). \]

Suppose that I_n is the value obtained using h=(1/2)^n and kh being all multiples of h for which 0<kh<M for some value M. We can choose M in advance so that all the derivatives

    \[ \frac{d^k}{dx^k}\psi(x)|_{x=M} \]

are small. We can choose M=10: as we have seen above, derivatives at 10 are not far off 10^{-2000}.

If h is halved, then the new values at which we will have to compute \psi(kh) and \psi'(kh) will be

    \[ h,3h,5h,7h,\ldots (10(2^{n+1})-1)h. \]

Let S_{n+1} be the sum

    \[ \sum h\psi'(kh)\Bigl(f(\psi(kh))+f(-\psi(kh))\Bigr) \]

for the new value h=(1/2)^{n+1} and the odd multiples of h just given. Then:

    \[ I_{n+1}=I_n/2+S_{n+1}. \]

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:

    \[ f(x)=x\sin(2\exp(2\sin(2\exp(2x)))). \]

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.

Leave a Reply

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