Adaptive Simpson and Romberg integration methods

[robve]

Moved this post to this new thread as suggested by others

I assume most readers are familiar with integration and quadratures. The methods and results reported here should not be surprising to anyone. Nevertheless, I hope this answers some of the questions frequently raised about numerical integration with these methods and how they compare.

As originally intended, this post is to be viewed in the context of vintage pocket computers with limited hardware, to find an approximation fast (within 1E-5 for example) and then we will to take a look what it would take to increase the error tolerance to say 1E-9 and what the pitfalls are of these methods. Many modern and better integration methods exist that I am very familiar with, but those may not run on machines with just 1K~2K of RAM for programs and data. These little machine's could run Simpson's rule that people could enter in BASIC, e.g. from listings in accompanying "application books" and from magazines. Likewise, this post includes listings for Romberg and Adaptive Simpson that will run on these little machines.

Adaptive Simpson is one of my favorite integration methods (on those machines). Adaptive Simpson is usually much more efficient than composite Simpson's rule and Romberg methods, since it uses fewer function evaluations in places where the integrand is well-approximated by a cubic function.

"Adaptive Simpson's method uses an estimate of the error we get from calculating a definite integral using Simpson's rule. If the error exceeds a user-specified tolerance, the algorithm calls for subdividing the interval of integration in two and applying adaptive Simpson's method to each subinterval in a recursive manner." - https://en.wikipedia.org/wiki/Adaptive_S...27s_method

The HP-71B Math Pac uses the Romberg method. One advantage of Romberg over Adaptive Simpson is its error estimate, which is based on the error of a high-degree interpolating polynomial constructed for the integrand. This error estimate can be reliably used, unless the integrand is much "too spiky", that is, has high order large nonzero derivatives, such in step functions and square and sawtooth wave functions that have non-differentiable points. Nevertheless, the final error estimate when convergence is not reached is rarely reported by an integration method. The listings at the end of this post include code to report the final error estimate for non-convergence.

For this post I will use true and tested implementations of Adaptive Simpson (qasi) and Romberg (qromb and qromo) that I wrote based on optimized versions of published algorithms in reputable sources such as Numerical Recipes and NIST FORTRAN programs. I believe these better and safer than those found in places like Wikipedia that have errors in the listings. I've tested these implementations and compared the results and the number of function evaluations to other Romberg integrators to corroborate the results. Nevertheless, there are many other good implementations possible. The listings of the implementations in C and in BASIC are further discussed below if you are interested and want to reproduce the results.

I received some feedback and questions about the Romberg implementations in SHARP BASIC that I showed and used in earlier posts. I will use the same programs here, but written in C for speed. The BASIC implementations are derived from these C versions and also included at the end of this post.

The integration methods take 5 parameters: the function f, the integration range a to b (finite), the maximum number of iterations n, and the error tolerance eps. There are two versions of Romberg, one for closed intervals (qromb) and one for open intervals (qromo).

To get started, we compute \( \int_{-1}^1 p(x)\,dx \) and collect the average number of integrand evaluations over 1000 different polynomials \( p(x) \) of degrees 1 to 30 to reach convergence:

$$
\begin{array}{rrrr}
\textit{eps} & \texttt{qromb} & \texttt{qromo} & \texttt{qasi} \\
\hline
10^{-3} & 36 & 96 & 23 \\
10^{-5} & 72 & 234 & 62 \\
10^{-7} & 127 & 448 & 195 \\
10^{-9} & 198 & 887 & 606 \\
\end{array}
$$

The averages are rounded to the nearest integer. Polynomials are randomly generated with coefficients \( -1\le a_i\le 1 \).

For polynomials (and many other functions), Adaptive Simpson converges quickly if the error threshold eps is not too tight to allow lower precision. Otherwise, Romberg over closed intervals appears to be best to achieve high accuracy for polynomials. Not entirely surprising because the methods essential perform polynomial interpolation and extrapolation to determine the definite integral.

In practice, Adaptive Simpson is really, really good to get a quick approximation, with only two minor drawbacks:
1. Adaptive Simpson may not produce usable error estimates of the result when the function is not well approximated on points by a cubic. (note: modifications are possible that produce an error estimate, but we're talking about the straightforward method).
2. Adaptive Simpson uses closed intervals and may fail to integrate improper integrals.

Improper integrals require integration over an open interval to avoid evaluating the endpoints, for example using Romberg with midpoint quadratures, see e.g. Numerical Recipes 2nd Ed. Ch.4.4 qromo.

If the integrand has singularities (non-evaluable points), then one can split up the interval in parts and use a quadrature suitable for open intervals to evaluate the parts separately. This may sometimes not be necessary because equidistant methods like Romberg can skip over them if you can control the points carefully. Adaptive Simpson will very likely not work, because it may subdivide its intervals until it hits the problematic points.

On the other hand, open intervals may not be necessary for integrands that produce a valid value at the endpoints, such as \( \int_0^1 \arctan(1/x)\,dx \) for example that evaluates to \( \pi/2 \) at \( x=0 \). Evaluation of this function at \( x=0 \) is fine as long as no floating point exception is raised to stop the program, because \( \arctan\infty = \frac{\pi}{2} \). Further, if the integrand's value at an endpoint is zero in the limit, then the floating point exception can be ignored (zero does not contribute).

I should add that Romberg implementations terminate at convergence, which can be too early if the number of points evaluated is insufficient to accurately assess convergence. I made sure that in my implementations the Romberg methods iterate at least 3 times (9 trapezoidal points) and 2 times (9 midpoints). While not always necessary, this increases the confidence we have in the error estimate when reaching convergence.

Next we evaluate \( \int_0^1 f(x)\,dx \) for a variety of functions \( f(x) \) with different characteristics.

In the next test, ten "regular" functions are tested (1 to 10) and six extreme ones (11 to 16) to hit the methods hard. For the hard cases no convergence can be expected. To test the extremes and make the methods work hard, we set the parameters fairly high: the Romberg qromb method tested performs at most N=23 subdivisions to evaluate the function up to 4 194 305 times. The Romberg qromo method tested performs at most N=15 subdivisions to evaluate the function up to 4 782 969 times. The Adaptive Simpson qasi method tested is recursive up to a depth of 21, which when evaluated in full performs 2²¹⁺¹+1 = 4 194 305 function evaluations.

The following table shows the number of function evaluations \( nnn \) performed by each method in the three columns for each of the 16 functions tested with tolerance \( 10^{-9} \), annotated with success or failure of the method to produce the desired result:

- \( nnn \) result with an acceptable error within the given tolerance \( 10^{-9} \)
- \( nnn^= \) exact result within \( 10^{-15} \)
- \( nnn^* \) max iterations or adaptive depth reached, residual error still within \( 10^{-9} \)
- \( nnn^e \) max iterations or adaptive depth reached, high residual error within \( 10^{-5} \)
- \( nnn^\dagger \) failed with incorrect result
- \( \dagger \) fatal: failed with NaN

\begin{array}{cccccc}
& f(x) & \int_0^1 f(x)\,dx & \texttt{qromb} & \texttt{qromo} & \texttt{qasi} \\
\hline
1 & x & \frac12 & 9^= & 9^= & 5^= \\
2 & x^3-2x^2+x & \frac{1}{12} & 9^= & 9^= & 5^= \\
3 & 4/(1+x^2) & \pi & 65 & 243 & 153 \\
4 & 1/(1+x) & \log2 & 33 & 243 & 97 \\
5 & \textrm{e}^x & \textrm{e}-1 & 17 & 81 & 65 \\
6 & \sin(x\pi) & \frac{2}{\pi} & 65 & 243 & 217 \\
7 & \arccos x & 1 & 262\,145 & 531\,441 & 489^* \\
8 & \arctan(1/x) & (2\log2+\pi)/4 & 33 & 81 & 93 \\
9 & \log x & -1 & \dagger & 4\,782\,969^e & \dagger \\
10 & 16x^4\log(2x+\sqrt{4x^2+1}) & 4.0766820599 & 65 & 243 & 449 \\
11 & \lfloor 10x\rfloor \mod 2 & \frac12 & 65 & 4\,782\,969^e & 705^* \\
12 & \lfloor 8x\rceil \mod 2 & \frac12 & 257 & 4\,782\,969^e & 5^\dagger \\
13 & \lfloor 6x\rceil \mod 2 & \frac12 & 65 & 4\,782\,969^e & 441 \\
14 & \lfloor 10x\rfloor & \frac12 & 4\,194\,305^e & 4\,782\,969^e & 441^e \\
15 & x - \lfloor 10x\rfloor & 4\frac12 & 4\,194\,305^e & 4\,782\,969^e & 441^\dagger \\
16 & \textrm{random} & \approx\frac12 & 4\,194\,305^e & 4\,782\,969^e & 4\,194\,305^\dagger \\
\end{array}

Trigs are evaluated in radians, log has base e and \( \lfloor x \rceil \) means round to nearest integer.

Integrals 1 and 2 are exact, as can be expected, with 9 point evaluation for Romberg to ensure accuracy (otherwise only 3 points are evaluated for integrals 1 and 2). Integral 2 over a cubic requires 5 points for Adaptive Simpson which is exact for cubics. Integrals 3 to 7 are proper and pose no significant issues. Integrals 8 and 9 are improper, but 8 can be evaluated on a closed interval. Integral 10 is from Numerical Recipes for reference and verification. Integrals 11 to 15 have spiky integrands and are "integration busters". Integrands 11 to 13 are square wave functions that can fool a method's convergence check. Integrand 14 is a sawtooth function and integrand 15 is a step function.

When we tally the results for this experiment, we see that Romberg wins if we count the fewest function evaluations to obtain a reliable result within 10^-6 relative error. Not surprising because we set the error threshold to 10^-9 to make these methods work hard. On the other hand, even in this more extreme case Adaptive Simpson is still fast and very good most of the time as long as the functions are well behaved, i.e. can be approximated by a cubic at the points evaluated.

While this is not an exhaustive comparison, it shows the successes and pitfalls of integration methods and that comparisons are sometimes apples to oranges. It depends on the integrand and interval what method we should to pick.

More modern adaptive methods based on Gauss–Kronrod quadrature and Clenshaw–Curtis quadrature are now generally preferred.

For the remainder of the post I will comment on the C and BASIC implementations of qromb, qromo, and qasi. You may freely use any of the code here and please ask questions or PM me if something is not clear.

You may ask how it is possible to implement Adaptive Simpson in old BASIC dialects when BASIC has no recursive functions and a limited call stack depth for GOSUB?

It is not too difficult when we use a stack to save and restore the recursive call parameters. First, we write the code in C for clarity:

Code:

double qasi(double (*f)(double), double a, double b, double eps) {
  double fa = f(a);
  double fm = f((a+b)/2);
  double fb = f(b);
  double v  = (fa + 4*fm + fb)*(b-a)/6;
  return qas_(f, a, b, fa, fm, fb, v, eps, 20, 0); // last 0 only for the tail-recursive version!
}

double qas_(double (*f)(double), double a, double b, double fa, double fm, double fb, double v, double eps, int n) {
  double h  = (b-a)/2;
  double f1 = f(a + h/2);
  double f2 = f(b - h/2);
  double sl = h*(fa + 4*f1 + fm)/6;
  double sr = h*(fm + 4*f2 + fb)/6;
  double s  = sl+sr;
  double d  = (s-v)/15;
  double m  = a+h;
  if (n <= 0 || fabs(d) < eps)
    return s + d; // convergence or max depth exceeded
  return qas_(f, a, m, fa, f1, fm, sl, eps/2, n-1) + qas_(f, m, b, fm, f2, fb, sr, eps/2, n-1);
}

Next, we optimize for tail-recursion by passing along the grand total t as a variable to be updated:

Code:

double qas_(double (*f)(double), double a, double b, double fa, double fm, double fb, double v, double eps, int n, double t) {
  double h  = (b-a)/2;
  double f1 = f(a + h/2);
  double f2 = f(b - h/2);
  double sl = h*(fa + 4*f1 + fm)/6;
  double sr = h*(fm + 4*f2 + fb)/6;
  double s  = sl+sr;
  double d  = (s-v)/15;
  double m  = a+h;
  if (n <= 0 || fabs(d) < eps)
    return t + s + d; // convergence or max depth exceeded
  t = qas_(f, a, m, fa, f1, fm, sl, eps/2, n-1, t);
  return qas_(f, m, b, fm, f2, fb, sr, eps/2, n-1, t);
}

Finally, we convert the recursive version to a non-recursive version using a stack stk and stack pointer sp, where we also replaced the check for depth n <= 0 with a comparison on eps that is halved at each level (making one or the other redundant):

Code:

double qas_(double (*f)(double), double a, double b, double fa, double fm, double fb, double v, double eps, int n, double t, int sp, double stk[]) {
ntry:
  double h  = (b-a)/2;
  double f1 = f(a + h/2);
  double f2 = f(b - h/2);
  double sl = h*(fa + 4*f1 + fm)/6;
  double sr = h*(fm + 4*f2 + fb)/6;
  double s  = sl+sr;
  double d  = (s-v)/15;
  double m  = a+h;
  // convergence or max depth exceeded?
  if (eps < 1E-15 || fabs(d) < eps) {
    // exit if stack is empty
    if (sp <= 0)
      return t + s + d;
    // pop parameters for the second tail-recursive call
    eps = stk[--sp]; s = stk[--sp]; fb = stk[--sp]; fm = stk[--sp]; fa = stk[--sp]; b = stk[--sp]; a = stk[--sp];
    goto entry;
  }
  eps /= 2;
  // save the parameters for the second "recursive" call on the stack to pop when "returning"
  stk[sp++] = m; stk[sp++] = b; stk[sp++] = fm; stk[sp++] = f2; stk[sp++] = fb; stk[sp++] = sr; stk[sp++] = eps;
  // the first "recursive" call
  b = m; fb = fm; fm = f1; s = sl;
  goto entry;
}

The Adaptive Simpson BASIC program, here for the SHARP PC-1350, but easily converted to other BASIC dialects such as HP-71B, note that A(27..166) is SHARP array (indexed by stack pointer O) that does not require DIM:

Code:

' A,B        range
' D          (s-v)/15
' E          error 1E-9, changed afterwards
' F$         function label to integrate
' H          interval
' I          fa=f(a)
' J          fm=f((a+b)/2)
' K          fb=f(b)
' L          sl=h*(fa+4*f1+fm)/6
' M          a
' N          b
' O          stack pointer
' P          f1=f(a+h/2)
' Q          f2=f(b-h/2)
' R          sr=h*(fm+4*f2+fb)/6
' S          s=sl+sr
' T          total sum
' U          partial sum v
' X          scratch
' A(27..166) auto-array stack for up to 22 recursive calls with 7 parameters per call (E=2^-20*1E-9<1E-15)

100 "QASI" E=1E-9: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init fa=f(a),fb=f(b),fm=f((a+b)/2)
130 X=A: GOSUB F$: I=Y,X=(A+B)/2: GOSUB F$: J=Y,X=B: GOSUB F$: K=Y
140 U=(B-A)*(I+4*J+K)/6,M=A,N=B,T=0,O=27
' recursive call entry
150 H=(N-M)/2,X=M+H/2: GOSUB F$: P=Y,X=N-H/2: GOSUB F$: Q=Y
160 L=H*(I+4*P+J)/6,R=H*(J+4*Q+K)/6,S=L+R,D=(S-U)/15: IF E<1E-15 OR ABS D<E LET T=T+S+D: GOTO 190
' save the parameters for the second "recursive" call on the stack to pop when "returning"
170 E=E/2,A(O)=M+H,A(O+1)=N,A(O+2)=J,A(O+3)=Q,A(O+4)=K,A(O+5)=R,A(O+6)=E,O=O+7
' the first "recursive" call
180 N=M+H,K=J,J=P,U=L: GOTO 150
' exit if stack is empty
190 IF O<=27 LET Y=T: PRINT Y: END
' if stack not empty then pop parameters for the second tail-recursive call
200 O=O-7,M=A(O),N=A(O+1),I=A(O+2),J=A(O+3),K=A(O+4),U=A(O+5),E=A(O+6): GOTO 150

300 "F1" Y=4/(1+X*X): RETURN
310 "F2" Y=EXP X: RETURN
320 "F3" Y=X-INT X: RETURN
330 "F4" Y=X^4*LN(X+SQR(X*X+1)): RETURN

The Romberg integration methods in C as used in the experiments with n=23 and n=15 for qromb and qromo, respectively, and eps=1e-9:

Code:

double qromb(double (*f)(double), double a, double b, int n, double eps) {
  double R1[n], R2[n];
  double *Ro = &R1[0], *Ru = &R2[0];
  double h = b-a;
  int i, j;
  Ro[0] = (f(a)+f(b))*h/2;
  for (i = 1; i < n; ++i) {
    unsigned long long k = 1UL << i; // k = 2^(i-1)
    unsigned long long s = 1;
    double sum = 0;
    double *Rt;
    h /= 2;
    for (j = 1; j < k; j += 2)
      sum += f(a+j*h);
    Ru[0] = h*sum + Ro[0]/2;
    for (j = 1; j <= i; ++j) {
      s *= 4;
      Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
    }
    if (i > 2 && fabs(Ro[i-1]-Ru[i]) <= eps*fabs(Ru[i])+eps)
      return Ru[i];
    Rt = Ro;
    Ro = Ru;
    Ru = Rt;
  }
  return Ro[n-1]; // no convergence, return best result, error is fabs((Ru[n-2]-Ro[n-1])/Ro[n-1])
}

double qromo(double (*f)(double), double a, double b, int n, double eps) {
  double R1[n], R2[n];
  double *Ro = &R1[0], *Ru = &R2[0];
  double h = b-a;
  int i, j;
  unsigned long long k = 1;
  Ro[0] = f((a+b)/2)*h;
  for (i = 1; i < n; ++i) {
    unsigned long long s = 1;
    double sum = 0;
    double *Rt;
    k *= 3;
    h /= 3;
    for (j = 1; j < k; j += 3)
      sum += f(a+(j-1)*h+h/2) + f(a+(j+1)*h+h/2);
    Ru[0] = h*sum + Ro[0]/3;
    for (j = 1; j <= i; ++j) {
      s *= 9;
      Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
    }
    if (i > 1 && fabs(Ro[i-1]-Ru[i]) <= eps*fabs(Ru[i])+eps)
      return Ru[i];
    Rt = Ro;
    Ro = Ru;
    Ru = Rt;
  }
  return Ro[n-1]; // no convergence, return best result, error is fabs((Ru[n-2]-Ro[n-1])/Ro[n-1])
}

The same Romberg integration methods in BASIC, with lower default N=10 for 512 function evaluations max and N=7 for 729 point function evaluations max to limit the CPU time in case of slow or non-convergence:

Code:

' VARIABLES
'  A,B   range
'  F$    function label to integrate
'  Y     result
'  E     relative error: integral = Y with precision E (attempts E = 1E-10)
'  H     step size
'  N     max number of Romberg steps (=10)
'  I     iteration step
'  U     current row
'  O     previous row
'  J,S,X scratch
'  A(27..26+2*N) scratch auto-array

100 "QROMB" E=1E-9,N=10: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init and first trapezoidal step
130 H=B-A,X=A: GOSUB F$: S=Y,X=B: GOSUB F$: O=27,U=O+N,A(O)=H*(S+Y)/2,I=1
' next trapezoidal step
140 H=H/2,S=0
150 FOR J=1 TO 2^I STEP 2: X=A+J*H: GOSUB F$: S=S+Y: NEXT J
' integrate and combine with previous results
160 A(U)=H*S+A(O)/2,S=1
170 FOR J=1 TO I: S=4*S,A(U+J)=(S*A(U+J-1)-A(O+J-1))/(S-1): NEXT J
' loop until convergence
180 IF I>2 LET Y=A(U+I): IF ABS(Y-A(O+I-1))<=E*ABS Y+E PRINT Y: END
190 J=O,O=U,U=J,I=I+1: IF I<N GOTO 140
' no convergence, output result with error estimate
200 E=ABS((Y-A(U+N-2))/(ABS Y+E)): PRINT Y,E: END

' VARIABLES
'  A,B   range
'  F$    function label to integrate
'  Y     result
'  E     relative error: integral = Y with precision E (attempts E = 1E-10)
'  H     step size
'  N     max number of Romberg steps (=7)
'  I     iteration step
'  U     current row
'  O     previous row
'  J,S,X scratch
'  A(27..26+2*N) scratch auto-array

100 "QROMO" E=1E-9,N=7: INPUT "f=F";F$: F$="F"+F$
110 INPUT "a=";A
120 INPUT "b=";B
' init and first midpoint step
130 H=B-A,X=A+H/2: GOSUB F$: O=27,U=O+N,A(O)=H*Y,I=1
' next midpoint step
140 H=H/3,S=0
150 FOR J=1 TO 3^I STEP 3: X=A+(J-.5)*H: GOSUB F$: S=S+Y,X=A+(J+1.5)*H: GOSUB F$: S=S+Y: NEXT J
' integrate and combine with previous results
160 A(U)=H*S+A(O)/3,S=1
170 FOR J=1 TO I: S=9*S,A(U+J)=(S*A(U+J-1)-A(O+J-1))/(S-1): NEXT J
' loop until convergence
180 IF I>2 LET Y=A(U+I): IF ABS(Y-A(O+I-1))<=E*ABS Y+E PRINT Y: END
190 J=O,O=U,U=J,I=I+1: IF I<N GOTO 140
' no convergence, output result with error estimate
200 E=ABS((Y-A(U+N-2))/(ABS Y+E)): PRINT Y,E: END

Hope this is useful.

Enjoy!

- Rob

Edit: minor additions to clarify n and N values in the implementations. Edit2: the initial post assumes readers are familiar with integration methods; I added some paragraphs and table to help casual readers.