The tanh-sinh quadrature

[emece67]

Hi all,

This post is aimed to those most "mathematical" of us. On it I roughly describe a numerical integration method I was just aware of, and that I think is very promising to solve integrals in programmable calculators.

Since the high school days I am fascinated with the quadrature methods. For me it is a kind of magic that, given some integral:
\[I = \int_a^bf(x)dx\]
you can approximate it with:
\[I \approx \sum_iw_if(x_i)\]
How's that than you can condense all the "information" carried by function \(f(x)\) over interval \((a, b)\) with a discrete and finite sum? How's that some set of weight and abscissa pairs \(\{w_i, x_i\}\) can gave a much better approximation than another set?

Still remember me manually (with the help of a Casio fx-180p) solving the system of non-linear equations defining the abscissas for the Chebyshev formula just to get such abscissas with more digits than those in my, then, textbook (Piskunov).

Well, back to the present, a few days ago I was aware of the existence of another quadrature method, the so called tanh-sinh quadrature, that seemed quite elegant to me. Some trials convinced me that such method is also precise, fast, robust &, surely, programmable in the little calculating machines we love. As this method (proposed on 1974) seems not to be described in many books, but only on some technical papers, I want to describe it here hoping some of you find it as interesting as I did (provided you still do not know it!)

Let's go, suppose the some integral needs to be computed:
\[I = \int_a^bf(x)dx\]
After the usual variable change \(x={b+a\over2}+{b-a\over2}r\) you get:
\[I = {b-a\over2}\int_{-1}^1f\left({b+a\over2}+{b-a\over2}r\right)dr\]
Now, the next variable change is performed: \(r=g(t)=\tanh\left({\pi\over2}\sinh t\right)\). Thus:
\[I={b-a\over2}\int_{-\infty}^\infty f\left({b+a\over2}+{b-a\over2}g(t)\right)g'(t)dt\]
Now, we compute this last integral as an infinite sum of rectangles, each of width \(h\):
\[I\approx{b+a\over2}h\sum_{i=-\infty}^\infty f\left({b+a\over2}+{b-a\over2}g(hi)\right)g'(hi)={b+a\over2}h\sum_{i=-\infty}^\infty f\left({b+a\over2}+{b-a\over2}r_i\right)w_i\]
with the abscissas \(r_i = \tanh\left({\pi\over2}\sinh(hi)\right)\) and the weights \(w_i = {\pi\over2}{\cosh(hi)\over\cosh^2\left({\pi\over2}\sinh(hi)\right)}\)

But, as \(i\) departs from 0, the weights \(w_i\) behave as:
\[w_i\rightarrow{\pi\over e^{{\pi\over2}e^{|hi|}-|hi|}}\]
that vanishes really fast because of the double exponential on the denominator. Thus, the infinite summation above can be approximated by the finite one:
\[I\approx{b+a\over2}h\sum_{i=-N}^N f\left({b+a\over2}+{b-a\over2}r_i\right)w_i\]
where \(N\) is (hopefully) a small integer.

To implement this algorithm when working with \(d\) decimal digits, you start with \(h = 1\) (the rectangles in the \(t\) domain are of width 1). Then compute the required \(N\) as the largest integer such that:

1. \(w_N\) is below \(\epsilon=10^{-d}\)
2. and \(r_N\), expressed with \(d\) decimal digits, is strictly below 1
3. and \({b + a\over2}+{b - a\over2}r_N\), expressed with \(d\) decimal digits, is strictly below \(b\)
   

These last two requirements ensure that the function \(f\) is not evaluated at the ends of the interval \((a, b)\). As \(hi\) increases, \(w_i\) rapidly vanishes to 0, but the abscissas \(r_i\) do also rapidly approach the interval ends \((-1, 1)\). Not evaluating \(f\) at the ends allows this method to manage integrals where the integrand has discontinuities (or infinite derivatives) at the ends (well, only if the double exponential decay of the weights "overweights" the integrand).
Of the three previous conditions, I have found that the first is slightly less strict than the others, thus from the others, \(N\) can be computed by:
\[N = \lfloor{1\over h}\cosh^{-1}\left({1\over \pi}\ln\left(2\cdot10^d\min\left(1, {b - a\over2}\right)\right)\right)\rfloor\]
With this value of \(N\) the summation is performed and you get a first approximation for \(I\).

Now you repeat this process, but this time halving the rectangle width \(h\). You will halve such width more times, so we define the "level" of the algorithm (or the iteration) as \(k\) and let \(h=2^{-k}\) (the initial \(k\) is 0). At each level you compute an approximation \(I_k\) for \(I\). When computing \(I_k\) you only need to compute half the abscissas and weights, the odd ones, as the even ones where included in the summation for \(I_{k-1}\), So \(I_k\) is \(I_{k-1}\) plus the sum for the new, odd, abscissas and weights. Also, the weights are even functions, so \(w_{-i}=w_i\).

The algorithm stops when the relative difference between \(I_k\) and \(I_{k-1}\) is below \(10^{-d\over2}\). The algorithm doubles the number of correct digits at each level.

In my tests, when working with 16 decimal digits, the algorithm stops at level \(k=3\) having performed 51 function evaluations. At 32 digits it stops at level 4 with 123 function evaluations. The result is usually correct to 15 or 31 digits, but sometimes (in my tests, this happens sometimes when the integrand goes to infinite at one of the interval ends) it gives half the number of correct digits.

There is an obvious drawback in this algorithm: the number of hyperbolic functions that must be computed to get the abscissas and weights. It amounts to twice the number of function evaluations (102 hyperbolics in my 16 digits tests). I have not compared yet the overall performance of this method against others (specifically, against the Romberg method).

I firstly heard of this algorithm when reading the documentation of the Python package mpmath, that I was using to perform comparisons between different root solvers at different digit counts.

There is below some Python code you can reuse to play a bit with this method. As you can see, the final code is compact & straightforward, so I think it can be translated easily to many programmable calculators (the wp34s is my target).

Further reading:
http://crd-legacy.lbl.gov/~dhbailey/dhbp...h-sinh.pdf

http://www.davidhbailey.com/dhbpapers/quadrature-em.pdf

Borwein, Bailey & Girgensohn, “Experimentation in Mathematics - Computational Paths to Discovery”, A K Peters, 2003, pages 312-313.

H. Takahasi and M. Mori. "Double Exponential Formulas for Numerical Integration." Publications of RIMS, Kyoto University 9 (1974), 721–741. (Available online)

Regards.

Code:


from mpmath import *
import time

# number of digits
mp.dps = 16
# repeat this # of times (to get better time estimations)
nloops = 100

# test equations
# equation =  0 => x*log(1 + x); 0; 1; 1/4
#             1 => x**2*atan(x); 0; 1; (pi - 2 + 2ln2)/12
#             2 => exp(x)cos(x); 0; pi/2; (exp(pi/2) - 1)/2
#             3 => atan(sqrt(2 + x**2))/(1 + x**2)/sqrt(2 + x**2); 0; 1; 5pi**2/96
#             4 => sqrt(x)ln(x); 0; 1; -4/9
#             5 => sqrt(1 - x**2); 0; 1; pi/4
#             6 => sqrt(x)/sqrt(1 - x**2); 0; 1; 2sqrt(pi)Gamma(3/4)/Gamma(1/4)
#             7 => ln(x)**2; 0; 1; 2
#             8 => ln(cos(x)); 0; pi/2; -pi*ln(2)/2
#             9 => sqrt(tan(x)); 0; pi/2; pi*sqrt(2)/2

print

# try all equations
for equation in range(10):
  if equation == 0:
    # limits
    a = mpf('0')
    b = mpf('1')
    # function
    def f(x):
      return x*ln(1 + x)
    # expected result (not used during computation)
    s = mpf('1/4')
  if equation == 1:
    a = mpf('0')
    b = mpf('1')
    def f(x):
      return x**2*atan(x)
    s = (pi() - 2 + 2*ln(mpf('2')))/12
  if equation == 2:
    a = mpf('0')
    b = pi()/2
    def f(x):
      return cos(x)*exp(x)
    s = (exp(pi()/2) - 1)/2
  if equation == 3:
    a = mpf('0')
    b = mpf('1')
    def f(x):
      return atan(sqrt(2 + x**2))/(1 + x**2)/sqrt(2 + x**2)
    s = 5*pi()**2/96
  if equation == 4:
    a = mpf('0')
    b = mpf('1')
    def f(x):
      return sqrt(x)*ln(x)
    s = mpf('-4/9')
  if equation == 5:
    a = mpf('0')
    b = mpf('1')
    def f(x):
      return sqrt(1 - x**2)
    s = pi()/4
  if equation == 6:
    a = mpf('0')
    b = mpf('1')
    def f(x):
      return sqrt(x)/sqrt(1 - x**2)
    s = 2*sqrt(pi())*gamma(mpf('3/4'))/gamma(mpf('1/4'))
  if equation == 7:
    a = mpf('0')
    b = mpf('1')
    def f(x):
      return ln(x)**2
    s = mpf('2')
  if equation == 8:
    a = mpf('0')
    b = pi()/2
    def f(x):
      return ln(cos(x))
    s = -pi()*ln(mpf('2'))/2
  if equation == 9:
    a = mpf('0')
    b = pi()/2
    def f(x):
      return sqrt(tan(x))
    s = pi()*sqrt(mpf('2'))/2

  # to measure algorithm execution time
  tt0 = time.time()
  # repeat nloops times to get better time estimations
  for m in range(nloops):
    # counters
    tnfe = 0  # counts function evaluations
    hyp = 0   # counts hyperbolics

    # we need a < b
    a, b = (a, b) if b > a else (b, a)
    # x = bpa2 + bma2*r
    bpa2 = (b + a)/2
    bma2 = (b - a)/2

    # epsilon
    eps = mpf('10')**(-mp.dps)
    # convergence threshold
    thr = mpf('10')**(-mp.dps/2)

    pi2 = pi()/2

    # (approx) maximum t that yields
    # the maximum representable r & x
    # values strictly below the upper
    # limits +1 (for r) and b (for x)
    tm = asinh(ln(2*min(1, bma2)/eps)/(2*pi2))
    hyp += 2

    # level
    k = 0
    # maximum allowed level
    maxlevel = int(ceil(log(mp.dps, 2))) + 2
    # ss is the final result
    # 1st addend at order 0
    ss = f(bpa2)
    tnfe += 1
    # "initial" previous computed result, used
    # in the convergence test
    sshp = 2*ss + 1
    # progress thru levels
    while k <= maxlevel:
      h = mpf('2')**-k
      N = int(floor(tm/h))
      j = 1
      while j <= N:
        t = j*h
        csh = pi2*sinh(t)
        ch = cosh(t)
        w = ch / cosh(csh)**2
        r = tanh(csh)
        hyp += 4
        fp = f(bpa2 + bma2*r)
        fm = f(bpa2 - bma2*r)
        p = w*(fp + fm)
        tnfe += 2
        ss += p
        # at level 0 must sweep all j values,
        # at other levels only the odd ones
        j += 2 if k > 0 else 1
      # converged?
      if abs(2*sshp/ss - 1) < thr: break
      # no, advance to next level
      k += 1
      # store the just computed approximation
      # for the next convergence test
      sshp = ss
    # apply constant coefficients
    ss *= bma2*pi2*h
    # done, measure time
    tt1 = time.time()

  # print results
  print equation, tnfe, hyp, k, (tt1 - tt0)/nloops, ss, ss - s, ss/s - 1