Integrating with Insight
Vector 9:3, page 113
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 113
This article was triggered by the function for Romberg integration given by C.E. Williams in Education Vector 9.1 with a note expressing hope for an improvement — well let’s try!
First Williams’ comparison of Romberg integration with the Simpson Rule is not one of like with like — Simpson’s Rule and the Trapezium rule on which Williams’ algorithm is based, are base integration techniques. Romberg’s method by contrast is a refinement technique which can be exercised on one or other of the base methods. In fact Romberg with one iteration is Simpson’s Rule, so the former is necessarily slower provided there is more than one iteration.
The refinement technique commences every time a new base integration is performed and consists of taking a scalar s which will in fact be the result of the new base integral evaluation and a vector v which will be the sequence of refinement steps since the last base integration and constructing a result vector r where
and so on, a process somewhat like (but not directly translatable into) that for evaluating s⊥v. The quantities k[i] are a sequence of constants — 1/3, 1/15, 1/63,… if the base technique is the Trapezium rule, 1/15, 1/127, 1/4096, … for the Simpson Rule. We are thus now in a position to define the refinement procedure by a function REFINE:
[0] r←v REFINE s
[1] ⍝ global variable Q = 1(Trap),2(Simpson)
[2] →L1 IF 0=⍴v←,v
[3] r←,(¯1↓v)REFINE s
[4] →0 r←r,(¯1↑r)+(¯1↑r-v)÷¯1+2*2×Q×⍴r
[5] L1:r←s
The function REFINE uses only ISO-APL functions. From this point however the power of APL2 to define operators is invaluable for the purpose of identifying the algorithmic structure.
The next step is to define the base methods themselves. The algorithms are so similar that it is sensible to combine them in a single operator METHOD whose second operand is 1 for Trapezium, 2 for Simpson:
[0] Z←L(P METHOD Q)R;T
[1] ⍝ L=no. intvs, P=integrand, Q=1(Trap),2(Simp), R=range of integr'n
[2] Z←T×+/(÷Q+1)×(1,((L-1)⍴(-Q)↑4 2),1)×P¨(↑R)+(T←-R-.÷L)×0,⍳L
The minimum number of intervals is 1 for Trapezium, 2 for Simpson.
Now the Romberg process can be described in terms of the base method METHOD and the refinement technique REFINE:
[0] Z←N(P ROMB Q)R
[1] ⍝ N=no of iter'ns, P=integrand, Q=1(trap),2(Simp), R=range of int'n
[2] →L1 IF N=0
[3] →0 Z←((N-1)P ROMB Q R)REFINE(2*N)(P METHOD Q)R
[4] L1:Z←1 P METHOD Q R
I prefer to identify a new iteration as starting with each invocation of the base technique which is the sense in which N is used in my ROMB. Williams considers each item in a new "REFINE" vector (what I call a refinement step) as an iteration so that 10 of his iterations are equivalent to 4 of mine or more generally \(\tfrac{1}{2}n.(n+1)\) of his are equivalent to \(n\) of mine. The operator ROMB returns a base integration together with the complete set of refinement steps which follow it, so taking the function which Williams uses as an example:
define the APL equivalent function:
and the numerical equivalence of his ROMB can be confirmed by:
The relative execution performance of my ROMB and Williams’ function depend on the integrand. Over a range of functions those for my ROMB are a small fraction of his — in most cases under 10%.
Rather than compute a given number of iterations, it is in practice more usual to seek an integration procedure which terminates when two iterates agree to a specified degree of precision EPS. An appropriate operator for this purpose is:
[0] Z←EPS(P ROMBERG Q)R;T;N
[1] ⍝ P is integrnd, Q = 1 (trap), 2 (Simpson), R = range of integration
[2] Z←T←(1+N←0)P METHOD Q R
[3] L1:Z←(T←¯1↑Z)REFINE(2*N←N+1)(P METHOD Q)R
[4] →L1 IF EPS<|T-¯1↑Z
[5] Z←¯1↑Z
To count iterations insert a counter-incrementing statement between lines 3 and 4 of ROMB.
ROMBERG is not the only possible refinement technique. In essence it doubles the number of intervals at each integration whereas it may be the case that adequate precision has already been achieved in one half of the interval. An alternative refinement technique is so-called adaptive integration which divides the tolerance EPS equally between the two halves of the range, and tests each separately, proceeding only in that half or halves where further precision is necessary. This refinement technique is described by:
[0] Z←EPS(P ADAPT Q)R
[1] ⍝ P is integrand, Q = 1(trap),2(Simp), R = range of integrn
[2] Z←(2×Q)(P METHOD Q)R
[3] →0 IF EPS>|Z-Q(P METHOD Q)R
[4] Z←+/(0.5×EPS)(P ADAPT Q)¨0 1⌽¨R,¨0.5×+/R
The value of having several alternative techniques at one’s disposal is illustrated by integrating the following "difficult" function:
whose APL rendition is:
The Simpson rule (method 2) with 2 and 4 intervals gives results which agree to 4 decimal places:
which is confirmed by using ADAPT with Simpson which converges immediately:
However integrating to a precision of 5 decimal places using ROMBERG shows that this answer is wrong:
The above integration required 6 iterations. Using the Trapezium rule as base method to obtain the same precision
requires 14 iterations. ADAPT with the Trapezium rule converges to the correct answer following 1023 interval subdivisions, that is 2047 separate calls to ADAPT.
To summarise, compared with Williams’ function the number of APL lines has been reduced from 30 to 12, four techniques are available instead of one, execution times have been drastically cut, and, most important of all in my view, the APL code provides an insightful description of the true nature of the methods by separating them into their essential stages of approximation and refinement. For the practically minded, five minutes typing these 12 lines provides a remarkably comprehensive integration package — worth a year’s subscription to Vector in itself! I think I can modestly claim an improvement.