The Museum of HP Calculators

Numerical Integration for the HP-67

This program is Copyright © 1976 by Hewlett-Packard and is used here by permission. This program was originally published in the HP-67 Math Pac 1.

This program is supplied without representation or warranty of any kind. Hewlett-Packard 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.

Card Labels
Numerical Integration
Shift a↑b n →∫(a..b)f   i
Label h f(xj) →TRAP →SIMP  
Key A B C D E

Overview

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.

Discrete case

Let x0, x1, . . ., xn be n equally spaced points (xj = x0 + jh, j = 1, 2, . . ., n) at which corresponding values f(x0), f(x1), . . ., f(xn) 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(xj), 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(x0)+2Σ(j=1..n-1)f(xj) + f(xn)]

2. Simpson's rule:

∫(x0..xn)f(x) dx ≈ h/3[f(x0) + 4f(x1) + 2f(x2) + ... + 4f(xn-3) 2f(xn-2) + 4f(xn-1) +f(xn)]

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.

Explicit functions

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 x0 = a, x1 = x0 + jh, j = 1, 2, ..., n-1, and xn = b where h = (b-a)/n.

The ∫(a..b) f(x) dx is approximated by equation (2) above, Simpson's rule.

Up to five different functions fi(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 fi(x), as well as registers R1 through R8, RS0 through RS9, and the four stack registers. The functions should assume x is in the X-register upon entry. Two levels of subroutines are allowed in the functions fi(x), but recall that the only labels available are 1 through 5.

Functions fi(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:

  1. Load side 1 of Numerical Integration.
  2. Press GTO . 112.
  3. Press MERGE.
  4. Load your card with the functions fi(x).

Remarks

Note that the function values for the discrete case f(xj), 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 user-definable keys be pressed during the entry of the f(xj), lest the next f(xj) entered go into the wrong LBL B.

Instructions

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 x-values h A  
 4 Repeat this step for j = 0, 1, ..., n: Key in the function value at xj.  f(xj) 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      

Example 1

Given the values below for f(xj), 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   
xi 0 .25 .5 .75 1 1.25 1.5 1.75 2
f(xi) 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)

Example 2

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 π × 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.

The Program

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   ×
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←(b-a)/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 1-5.
112   RTN

Register Use

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