The Tip Of The Romberg

| 0 Comments

Last time[1] we took a first look at the numerical approximation of integrals with the trapezium rule, a rather simplistic algorithm with the unfortunate property that its accuracy was implicitly, rather than explicitly, specified.
As a user of numerical libraries, including my own, I would much rather have an algorithm that works to a given accuracy, or at the very least tries to, than have to figure it out for myself and I suggested that we might apply the lessons learned from our attempts to numerically approximate derivatives to do just that[2]-[6].

Polynomial Approximation Will Alleviate Your Calculus Blues

One of the most accurate approximations we came up with was to use polynomial approximation[4], not of the function we were differentiating, but of the symmetric finite difference as a function of \(\delta\)
\[ g_{f,x}(\delta) = \frac{f(x+\delta) - f(x-\delta)}{2\delta} \]
By extrapolating the approximating polynomial to zero we were able to reduce the approximation error to the least significant decimal floating point digit, at least for our, admittedly rather simple, test case.
Given that the more accurate approximations required the implementation of the rules of differentiation as transformations of data structures, something we simply can't do for integration, it is extremely tempting to use this kind of polynomial approximation for integrals and it should come as no surprise that it's already been done. The approach is known as Romberg's method[7] and approximates the trapezium rule as a function of the width \(h\)
\[ g_{f,x_1,x_{n+1}}(h) = \sum_{i=1}^n \tfrac{1}{2} \left(f\left(x_i\right) + f\left(x_{i+1}\right)\right) \times h \approx \int_{x_1}^{x_{n+1}} f(x) \; \mathrm{d}x \]
where \(\sum\) is the summation sign and the wavy equals sign means approximately equals.
Note that whilst we can only evaluate this function when
Derivation 1: The Taylor Series Of A Trapezium's Area
Given a trapezium from \(x_i\) to \(x_{i+1}\), define
\[ \begin{align*} \bar{x} &= \tfrac{1}{2}\left(x_i + x_{i+1}\right)\\ x_i &= \bar{x} - \tfrac{1}{2}h\\ x_{i+1} &= \bar{x} + \tfrac{1}{2}h \end{align*} \]
yielding Taylor series of the function at \(x_i\) and \(x_{i+1}\) of
\[ \begin{align*} f\left(x_i\right) &= f\left(\bar{x}\right) - \tfrac{1}{2} h f^\prime\left(\bar{x}\right) + \tfrac{1}{8} h^2 f^{\prime\prime}\left(\bar{x}\right) - \dots\\ f\left(x_{i+1}\right) &= f\left(\bar{x}\right) + \tfrac{1}{2} h f^\prime\left(\bar{x}\right) + \tfrac{1}{8} h^2 f^{\prime\prime}\left(\bar{x}\right) + \dots \end{align*} \]
and consequently the trapezium's area is given by
\[ \begin{align*} \tfrac{1}{2}\left(f\left(x_i\right) + f\left(x_{i+1}\right)\right) \times h &= \left(f\left(\bar{x}\right) + \tfrac{1}{8} h^2 f^{\prime\prime}\left(\bar{x}\right) + \dots\right) \times h\\ &= h f\left(\bar{x}\right) + \tfrac{1}{8} h^3 f^{\prime\prime}\left(\bar{x}\right) + \dots \end{align*} \]
Figure 1: Doubling The Number Of Trapeziums
\[ h = x_{i+1} - x_i = \frac{x_{n+1} - x_1}{n} \]
it doesn't stop us approximating it with a polynomial and extrapolating to the limit of \(h\) equal to zero.

A rather convenient property of the trapezium rule is that the Taylor series for the area of each trapezium about its midpoint on the \(x\) axis only contains odd powers of \(h\), as shown by derivation 1. Rearranging the above formula for \(h\), we have
\[ n = \frac{x_{n+1} - x_1}{h} \]
implying that their sum can be represented using a polynomial having only even powers of \(h\), halving the number of terms that we need to calculate.

Furthermore, if we double the number of trapeziums each time we add a new term to the polynomial then half of the function evaluations we'll need will have been made in the previous step, saving us yet more computational expense. This is illustrated by figure 1 in which the red lines represent a four trapezium approximation of an integral of the normal PDF and the green lines the additional function evaluations required for an eight trapezium approximation.

Putting It All Together

The implementation of ak.rombergIntegral, giving in listing 1, follows our now familiar approach of first checking the arguments, providing default values if necessary, and then returning a function that forwards to the actual implementation.

Listing 1: ak.rombergIntegral
ak.rombergIntegral = function(f, threshold, steps) {
 if(ak.nativeType(f)!==ak.FUNCTION_T) {
  throw new Error('invalid function in ak.rombergIntegral');
 }

 threshold = ak.nativeType(threshold)===ak.UNDEFINED_T ? Math.pow(ak.EPSILON, 0.75)
                                                       : Math.abs(threshold);
 if(isNaN(threshold)) {
  throw new Error('non-numeric threshold passed to ak.rombergIntegral');
 }

 steps = ak.nativeType(steps)===ak.UNDEFINED_T ? 18
                                               : ak.floor(Math.abs(steps));
 if(isNaN(steps)) {
  throw new Error('non-numeric steps passed to ak.rombergIntegral');
 }

 return function(x0, x1) {
  return rombergIntegral(f, x0, x1, threshold, steps);
 };
};

Here f is the function to be integrated, threshold is the convergence threshold and steps is the maximum number of iterations. Note that the convergence threshold defaults to three quarters of the digits of a double if it's not specified and the maximum number of iterations defaults to 18.

Listing 2 gives the implementation of rombergIntegral.

Listing 2: rombergIntegral
function rombergIntegral(f, x0, x1, eps, steps) {
 var terms = 5;
 var a = [];
 var n = 1;
 var i, j, h0, s0, s1;

 if(ak.nativeType(x1)===ak.UNDEFINED_T) {
  x1 = x0;
  x0 = 0;
 }

 if(!isFinite(x0) || !isFinite(x1)) {
  throw new Error('invalid range of integration in ak.rombergIntegral');
 }

 h0 = x1-x0;
 a.push(0.5*(f(x0)+f(x1))*h0);
 for(i=1;i<terms;++i) a.push(rombergStep(f, x0, x1, a[i-1], n*=2));
 s0 = rombergPredict(a, h0);

 a.shift();
 a.push(rombergStep(f, x0, x1, a[terms-2], n*=2));
 s1 = rombergPredict(a, h0/=2);

 for(i=0;i<steps && ak.diff(s0, s1)>eps;++i) {
  a.shift();
  a.push(rombergStep(f, x0, x1, a[terms-2], n*=2));
  s0 = s1;
  s1 = rombergPredict(a, h0/=2);
 }
 if(!(ak.diff(s0, s1)<=eps)) {
  throw new Error('failure to converge in ak.rombergIntegral');
 }

 return s1;
}

Note that, as we did in ak.trapeziumIntegral, if a single argument is passed to the integral we compute the integral from zero to it.

This implementation follows the algorithm as described in Numerical Recipes in C[8], using a fifth order polynomial to extrapolate the trapezium width to zero. The array a is initially filled with five successive trapezium rule approximations, calculated by the helper function rombergStep, and is kept at that length at each step by shifting off the first value before pushing on the next approximation for which we double the number of trapeziums, n. The extrapolation is calculated by rombergPredict, given together with rombergStep in listing 3, and we terminate the algorithm once the normalised difference of successive extrapolations falls below the convergence threshold eps or we exceed the maximum number of steps.

Listing 3: rombergStep And rombergPredict
function rombergStep(f, x0, x1, s0, n) {
 var h = (x1-x0)/n;
 var s1 = 0;
 var i;

 if(h===0) return s0;
 for(i=1;i<n;i+=2) s1 += f(x0 + i*h);
 return s0/2 + s1*h;
}

function rombergPredict(a, h) {
 var i = ak.nevilleInterpolate(0);
 var terms = a.length;
 var j;

 h *= h;
 for(j=1;j<terms && h>0;++j,h/=4) i.refine(h, a[j-1]);
 return h>0 ? i.refine(h, a[terms-1]) : a[j-1];
}

As noted above, rombergStep exploits the fact that we double the number of trapeziums at each step by only adding up the function's values at odd multiples of the current width, assuming that it is not zero, that is. The product of the width and the sum of the function's values at even multiples is simply half the value of the previous approximation, s0, since it was made with double the current width.
We use ak.nevilleInterpolate to extrapolate the successive trapezium rule approximations to a width of zero in rombergPredict from the squares of the trapezium widths, again as noted above. In the extremely unlikely event that h underflows, or in the rather less unlikely event that it starts at zero because the bounds of integration are equal, we bail out of the loop and return the trapezium rule approximation at which it equalled zero.

The implementation of ak.rombergIntegral can be found in RombergIntegral.js.

Using ak.rombergIntegral

The great advantage of ak.rombergIntegral over ak.trapeziumIntegral is that we specify the convergence threshold rather than the trapezium width giving us direct control over the approximation error, or at least a reasonable proxy for it.
To demonstrate how much easier this makes this algorithm to use we can simply apply it to the same quadratic function that we integrated with ak.trapeziumIntegral, which we know has an exact integral of
\[ \begin{align*} \int_{0}^{z} \left(ax^2 + bx + c\right) \; \mathrm{d}x &= \left[\tfrac{1}{3}ax^3 + \tfrac{1}{2}bx^2 + cx\right]_{0}^{z}\\ &= \tfrac{1}{3}az^3 + \tfrac{1}{2}bz^2 + cz \end{align*} \]
Program 1 follows the implementation of that of the previous post, with the first column of the output showing the approximate integral, the second the exact and the third the absolute difference between them, simply replacing ak.trapeziumIntegral with ak.rombergIntegral.

Program 1: Using ak.rombergIntegral

This represents a tremendous improvement in both the control of the approximation error and the computational expense of achieving accurate results, reaching more than ten decimal orders of magnitude better than those in the last post's example almost as quickly.

Now you might very well accuse me of having stacked the deck by using a polynomial integrand to demonstrate the accuracy of polynomial extrapolation of the integral, so program 2 does the same for the exponential function which has an exact integral of
\[ \int_0^z e^x \; \mathrm{d}x = \left[e^x\right]_0^z = e^z - e^0 = e^z - 1 \]
Program 2: Integrating The Exponential Function

Hopefully, this will have convinced you that the effectiveness of this algorithm isn't restricted to integrals of polynomials. However, I'm sorry to report that it is not universally applicable and there are many functions for which it fails to converge, the simplest of which are discontinuous functions. For example
\[ f(x) = \begin{cases} -1 & x < 0\\ +1 & x \geqslant 0 \end{cases} \]
Program 3 shows how Romberg's method fails to integrate this function.

Program 3: A Discontinuous Failure

Now this isn't entirely a bad thing since other algorithms will return inaccurate results rather than explicitly fail, an arguably much worse outcome.
Fortunately, if we know that a function is discontinuous and, crucially, where the discontinuities are, we can easily fix the problem by adding together integrals of the continuous regions of the function, as illustrated by program 4.

Program 4: A Discontinuous Success

Unfortunately, it's not really very much more difficult to find continuous functions that numerical integration algorithms struggle with. A trivial example is

\(\displaystyle{ f(x) = \sin(e^{x^2}) }\)

which, as \(x\) increases, moves from positive to negative so rapidly that it requires a huge number of function evaluations to accurately calculate the difference between the positive and negative areas between the curve and the \(x\) axis which, as demonstrated by program 5, means that the integral will terminate before it converges.

Program 5: A Continuous Failure

Badly behaved continuous functions like this are much, much harder to deal with than discontinuous functions and, rather irritatingly, it won't be long before we're faced with another one.

References

[1] Trapezium Artistry, www.thusspakeak.com, 2015.

[2] You're Going To Have To Think! Why [Insert Algorithm Here] Won't Cure Your Calculus Blues, www.thusspakeak.com, 2013.

[3] You're Going To Have To Think! Why Finite Differences Won't Cure Your Calculus Blues, www.thusspakeak.com, 2014.

[4] You're Going To Have To Think! Why Polynomial Approximation Won't Cure Your Calculus Blues, www.thusspakeak.com, 2014.

[5] You're Going To Have To Think! Why Computer Algebra Won't Cure Your Calculus Blues, www.thusspakeak.com, 2014.

[6] You're Going To Have To Think! Why Automatic Differentiation Won't Cure Your Calculus Blues., www.thusspakeak.com, 2014.

[7] Romberg, W., Vereinfachte Numerische Integration, Det Kongelige Norske Videnskabers Selskab Forhandlinger, Vol. 28, No. 7, Trondheim, 1955.

[8] Press, W.H. et al, Numerical Recipes in C (2nd ed.), Cambridge University Press, 1992.

Leave a comment