Previous page Next page Navigation bar

Programs

Mathematics Programs

Romberg Integration

Summary

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:

Is compatible CFX-9850G

Compatible with CFX-9850G

Is compatible FX-7400G

Compatible with FX-7400G

Detailed Description

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(h2) while their combination into a Simpson’s rule approximation has an error of O(h4). By combining a series of finer trapezoidal rule approximations in the same way, the approximation error can be made O(h2k) 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.

Usage

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

Integral from 0 to 1 of ln (1+x)

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

Program Source

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 non-negative 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
F-1->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.
((D-E)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
;
B-A->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[E-1]/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])/(E-1)->List 1[F]
                   ;     Apply the correction
Next               ;   [End of loop through rows]
D-1->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 ]

Previous page Top of page Next page Navigation bar

Copyright © 2001 Brian Hetrick
Page last updated 25 November 2001.

Brian’s Casio Calculator Corner

Home

Programs

Index

Linkage Conventions

Mathematics

Index

Factorization

Bisection

Newton’s Method

Complex Zeroes

Best Fraction Approximation

Pseudorandom Numbers

Romberg Integration

First Degree ODE

Rounding

Permutations

Finance

Physics

Utility

Tutorial

Puzzles

Site Information

Your Privacy

Site Map

E-mail

Site Technical Data