This program is Copyright © 1976 by HewlettPackard and is used here by permission. This program was originally published in the HP67 Math Pac 1.
This program is supplied without representation or warranty of any kind. HewlettPackard Company and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.
Numerical Integration 

Shift 
a^b 
n 
>_{(a..b)}f 

i 
Label 
h 
f(x_{j}) 
>TRAP 
>SIMP 

Key 
A 
B 
C 
D 
E 
This program will perform numerical integration whether a function is known explicitly or only at a finite number of equally spaced points (discrete case). The integrals of explicit functions are found using Simpson's rule; discrete case integrals may be approximated by either the trapezoidal rule or Simpson's rule.
Let x_{0}, x_{1}, . . ., x_{n} be n equally spaced points (x_{j} = x_{0} + jh, j = 1, 2, . . ., n) at which corresponding values f(x_{0}), f(x_{1}), . . ., f(x_{n}) of the function f(x) are known. The function itself need not be known explicitly. After input of the step size h and the values of f(x_{j}), j = 0, 1, . . ., n, then the integral
_{(x0..xn)} f(x)dx
may be approximated using
1. The trapezoidal rule:
_{(x0..xn)}f(x) dx ~ h/2[f(x_{0})+2_{(j=1..n1)}f(x_{j}) + f(x_{n})]
2. Simpson's rule:
_{(x0..xn)}f(x) dx ~ h/3[f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + ... + 4f(x_{n3}) 2f(x_{n2}) + 4f(x_{n1}) +f(x_{n})]
In order to apply Simpson's rule, n must be even. If n is not even, the calculator will halt displaying "Error" if D is pressed.
If an explicit formula is known for the function f(x), then the function may be keyed into program memory and numerically integrated by Simpson's rule. The user must specify the endpoints a and b of the interval over which integration is to be performed, and the number of subintervals n into which the interval (a,b) is to be divided. This n must be even; if it is not, Error will be displayed. The program will go on to compute x_{0} = a, x_{1} = x_{0} + jh, j = 1, 2, ..., n1, and x_{n} = b where h = (ba)/n.
The _{(a..b) }f(x) dx is approximated by equation (2) above, Simpson's rule.
Up to five different functions f_{i}(x), i = 1, ..., 5, may be loaded into program memory at one time under labels 1 through 5. The function to be integrated is selected by keying in a digit 1, 2, 3, 4, or 5, and pressing f E. The function under the appropriate label will then be selected. 112 program stops are available for defining the f_{i}(x), as well as registers R1 through R8, RS0 through RS9, and the four stack registers. The functions should assume x is in the Xregister upon entry. Two levels of subroutines are allowed in the functions f_{i}(x), but recall that the only labels available are 1 through 5.
Functions f_{i}(x) may be keyed into program memory after loading side 1 of Numerical Integration, or you may record these functions beforehand on a magnetic card and load them in the following manner:
Note that the function values for the discrete case f(x_{j}), j = 0, 1, ..., n, are keyed into B . There are actually three routines in the program which begin with LBL B, one for j = 0, one for j odd, and one for j even. It is important that no other userdefinable keys be pressed during the entry of the f(x_{j}), lest the next f(x_{j}) entered go into the wrong LBL B.
Step 
Instructions 
Input Data/Units 
Keys 
Output Data/Units 
1 
Load side 1. 

2 
For explicit functions, go to step 8; for discrete case, go to step 3 



DISCRETE  
3 
Key in spacing between xvalues 
h 
A 

4 
Repeat this step for j = 0, 1, ..., n: Key in the function value at x_{j}. 
f(x_{j}) 
B 
j 
5 
Compute the area by trapezoidal rule. 
C 
Trap 

6 
Compute the area by Simpson's rule. 
D 
Simp 

7 
For a new case, go to step 2. 



Explicit Functions 




8 
Load the function(s) into program memory (either key them in with LBL and RTN, or link from step 112.)  
9 
Specify the function i to be selected. 
i(1  5) 
f E 

10 
Key in the beginning and final endpoints of the integration interval. 
a 
ENTER 


b 
f A 

11 
Key in the number of subintervals. (must be even). 
n 
f B 

12 
Compute the area by Simpson's rule. 

f C 
Simp 
13 
To change a, b, or n, go to the appropriate step; for a new case, go to step 2 

Given the values below for f(x_{j}), j = 0, 1, ...,8, compute the approximations to the integral
_{(0..2)}f(x) dx
by the trapezoidal rule and by Simpson's rule. The value for h is 0.25.
i 
0 
1 
2 
3 
4 
5 
6 
7 
8 
x_{i} 
0 
.25 
.5 
.75 
1 
1.25 
1.5 
1.75 
2 
f(x_{i}) 
2 
2.8 
3.8 
5.2 
7 
9.2 
12.1 
15.6 
20 
Keystrokes Outputs .25 A 2 B 2.8 B 3.8 B 5.2 B 7 B 9.2 B 12.1 B 15.6 B 20 B C 16.68 *** (Trapezoidal) D 16.58 *** (Simpson's)
Find the value of
dx [0..2pi]  1  cos x + 0.25
for n = 10 and then for n = 16. Note that x is assumed to be in radians. For safety, if you work mostly in degrees, it is good programming practice to set the angular mode to radians at the beginning of the routine, then back to degrees at the end. Key the function in under LBL 1.
Keystrokes Outputs GTO . 112 Switch to PRGM. LBL 1 RAD COS 1 x<>y  .25 + 1/x DEG RTN Switch to RUN. 0 ENTER 2 pi x f A 10 f B 1 f E f C 8.22 *** (n=10) 16 f B f C 8.36 *** (n=16)
The exact solution is 8/3 = 8.38.
LINE KEYS 001 *LBL A Input h. 002 STO D 003 RTN 004 *LBL B Input f(X0) 005 STO 0 R0 contains TRAP sum. 006 STO 9 R9 contains SIMP sum. 007 0 n = 0 008 STO C 009 *LBL 9 010 RTN 011 *LBL B Input f(xj), j odd. 012 STO A R0<R0 + 2f(Xj) 013 GSB 6 014 ENTER 015 + 016 STO + 9 R9<R9 + 4f(Xj) 017 RCL C 018 1 019 + 020 STO C n<n+1 021 RTN 022 *LBL B 023 STO A Input f(Xj), j even. 024 GSB 6 R0<R0 + 2f(Xj) 025 STO + 9 026 RCL C 027 1 R9<R9 + 2f(Xj) 028 + 029 STO C n<n+1 030 GTO 9 Exit 031 *LBL C Compute trapezoidal area. 032 2 033 RCL 0 034 GTO 0 035 *LBL D Compute Simpson area. 036 RCL C 037 GSB 7 038 3 test n even. 039 RCL 9 040 *LBL 0 041 RCL A 042  2f(Xn) were added, so subtract f(Xn) 043 *LBL d 044 RCL D 045 x 046 X<>Y 047 ÷ 048 PRTX 049 RTN Area 050 *LBL a Input a and b. 051 STO B 052 X<>Y Store b. 053 STO A 054 RTN Store a. 055 *LBL b Input n. 056 STO C 057 *LBL 7 058 2 If n is not even, display 059 ÷ "Error" by call to E. 060 FRAC 061 X!=0? 062 GTO E 063 RTN Compute Simpson area. 064 *LBL c 065 RCL A 066 GSB (i) R0<f(a) 067 STO 0 068 RCL B 069 GSB (i) 070 STO + 0 R0<R0+f(b) 071 RCL B 072 RCL A 073 STO E Set initial x top a. 074  075 RCL C 076 ÷ 077 STO D h<(ba)/n 078 0 079 STO 9 R9 will count to n. 080 *LBL 8 081 RCL D 082 RCL E 083 + 084 STO E x<x+h. 085 GSB(i) 086 GSB 6 087 STO + 0 088 2 R0<R0+4f(X) 089 STO + 9 090 RCL C 091 RCL 9 092 X=Y? If R9=n, exit. 093 GTO 0 094 RCL D 095 RCL E 096 + 097 STO E 098 GSB (i) x<x+h 099 GSB 6 100 GTO 8 R0<R0+2f(X) 101 *LBL 0 102 3 103 RCL 0 104 GTO d Area = h/3 x R0 105 *LBL 6 Subroutine 106 ENTER 107 + 108 STO + 0 109 RTN 110 *LBL e R0<R0+2f(X) 111 STO I Input i to select function 15. 112 RTN
R0 Used R9 Used A f(Xj), a B b C n D h E x I Function i
Go back to the software library
Go back to the main exhibit hall