Programs: 
ROMBERG, 1ROMBERG 

Version: 
2.0, 25 November 2001 

Description: 
These programs provide Romberg integration of a specified function. The program 1ROMBERG can be called by any program. The program ROMBERG is intended to be used from the keyboard. 

Compatibility: 

Numerical integration of an arbitrary function over a finite interval divides the interval into a number of smaller regions, estimates the integral over each region, and adds the estimates. The “trapezoidal rule” approximates the function in the region by a straight line between the end points; “Simpson’s rule” approximates the function in two adjacent regions by a quadratic function connecting the two outer points and the common midpoint.
In addition to these geometric interpretations, it can be shown that a Simpson’s rule approximation with region width h can be derived from two trapezoidal rule approximations, one with region width h and one with region width 2h. The trapezoidal rule approximation has an error of O(h^{2}) while their combination into a Simpson’s rule approximation has an error of O(h^{4}). By combining a series of finer trapezoidal rule approximations in the same way, the approximation error can be made O(h^{2k}) for as large a value of k as desired. This combination of trapezoidal rule approximations is known as Romberg integration.
Although Romberg integration requires no more function evaluations than Simpson’s rule on an equally fine spacing, its error is much smaller. The calculator’s built in integration function is a Simpson’s rule calculation. Program ROMBERG provides noticeably higher accuracy than the calculator’s built in integration function.
Program ROMBERG interacts with the calculator keyboard and display, and is intended to be the user interface for program 1ROMBERG. You must program the function to be integrated into program F0. This program must accept the value at which the function is evaluated in X, and return the function value in Ans. It must not destroy any of the variables used by program 1ROMBERG, as shown in the “Resources Used” section.
To integrate the function from the keyboard, run program ROMBERG. Program ROMBERG prompts with the lines:
Min? Max? Rows?
Enter the lower bound of the definite integral, the upper bound of the definite integral, and the base 2 logarithm of the number of divisions to be used in the finest trapezoidal rule, respectively. The value for “rows” is interpreted as follows: 0 indicates 1 division, 4 indicates 16 divisions, 8 indicates 256 divisions, and so forth. This is analogous to the calculator’s integration function number of divisions parameter. The program will compute and display the Romberg integral of the function in F0 between the values Min and Max.
To find the definite integral of the function F0 from a program, place the lower bound in variable A, the upper bound in variable B, and the base 2 logarithm of the number of divisions in variable C, and call program 1ROMBERG. Program 1ROMBERG will return the value of the integral in Ans, while destroying variable contents as shown in the “Resources Used” section.
As an example, consider the integral
This value was computed with the calculator’s built in (Simpson’s rule) integration function and with program ROMBERG. The results, where n is the base 2 logarithm of the number of divisions used, are given in the table below.
n 
Calculator value 
Calculator error 
ROMBERG value 
ROMBERG error 

0 
Arg ERROR 
— 
0.34657 35903 
0.03972 77084 
1 
Ma ERROR 
— 
0.38583 46022 
0.00045 97589 
2 
0.3863 
0.00000 56389 
0.38628 78935 
0.00000 64676 
3 
0.38629 
0.00000 43611 
0.38629 43091 
0.00000 00520 
4 
0.38629 4 
0.00000 03611 
0.38629 43609 
0.00000 00002 
5 
0.38629 44 
0.00000 00389 
0.38629 43611 
0 
6 
0.38629 4361 
0.00000 00001 
Unchanged 
0 
7 
0.38629 43611 
0 
Unchanged 
0 
8 
Unchanged 
0 
Unchanged 
0 
The program is available as a text file with .CAT contents, or may be entered as shown below. A semicolon (“;”) marks the beginning of a comment, and is not to be entered into the calculator. Remember that these programs are copyrighted; see the copyright issues page for limitations on redistribution.
Program ROMBERG (85 bytes):
; Variables: ; A is lower bound of integration interval ; B is upper bound of integration interval ; C is number of rows in Romberg table ; Symbols: ; > is assignment arrow ; > is greater than relational operator "Min"?>A ; Get integration interval lower bound "Max"?>B ; Get integration interval upper bound If A>B ; Test if bounds are out of order Then A>C ; If so, B>A ; reverse the bounds C>B ; so they are in increasing order IfEnd ; [End of test if bounds are out of order] "Rows"?>C ; Get number of Romberg rows Int Abs C>C ; Ensure is nonnegative integer Prog "1ROMBERG" ; Do the integration
Program 1ROMBERG (294 bytes):
; Variables: ; A is lower bound of integration interval ; B is upper bound of integration interval ; C is number of rows in Romberg table ; D, E, F, G are work ; X is the argument to F0 ; List 1 is the column of the Romberg table ; Symbols: ; > is assignment arrow ; ^ is exponentiation operator ; / is division operator ; = is equality relational operator ; <= is less than or equal to relational operator ; > is greater than relational operator ; ; Phase 1: Accumulate the sums of function values ; Seq(0,D,1,C+1,1)>List 1 ; Get space for Romberg table column 2^C>D ; Get number of subdivisions For 0>E To D ; Loop through subdivision boundary points 1+Dim List 1>F ; Figure out the Romberg row index into which to D+E>G ; accumulate the function value at this point. Do ; The appropriate row is the last one if the G/2>G ; point index mod 2 is 1, the next to last if F1>F ; the point index mod 4 is 2, the second to LpWhile Frac G=0 ; last if the point index mod 8 is 4, etc. The If F<=0 ; point index is added to the number of indices Then 1>F ; to ensure something is found in this loop, IfEnd ; but the computed index may be 1 or 0, which ; must be corrected to 1. ((DE)A+EB)/D>X ; Get the boundary point Prog "F0" ; Evaluate the function List 1[F]+Ans>List 1[F] ; Add the function value into the appropriate ; Romberg row Next ; [End of loop through subdivision boundaries] ; ; Phase 2: Turn the sums of function values into trapezoidal rule ; approximations to the integral ; BA>D ; Get the integration interval size DList 1[1]/2>List 1[1] ; Apply the trapezoidal rule with one interval to ; the first row of the Romberg column For 2>E To Dim List 1 ; Loop through the Romberg column rows D/2>D ; The new interval is half as small as before List 1[E1]/2+DList 1[E]>List 1[E] ; Turn the sum of function values into the ; trapezoidal rule value Next ; [End of loop through Romberg column rows] ; ; Phase 3: Apply the Romberg recursion to refine the estimates ; C>D ; Get the number of rows 1>E ; Get the initial multiplier While D>0 ; Iterate once per row 4E>E ; Get this iteration's multiplier For 1>F To D ; Loop through first rows (EList 1[F+1]List 1[F])/(E1)>List 1[F] ; Apply the correction Next ; [End of loop through rows] D1>D ; Decrement the iteration index WhileEnd ; [End of iterate once per row] ; ; The Romberg integral is in List 1[1] ; List 1[1] ; Return the integral value
[ Previous page  Top of page  Next page ]
Copyright © 2001 Brian Hetrick
Page last updated 25 November 2001.