The Museum of HP Calculators


Solution To f(x)=0 On An Interval 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

Solution to f(x) = 0

Shift

 

 

 

 

 

Label

b->f(b)

c->f(c)

TOL

->root

x->f(x)

Key

A

B

C

D

E

Overview

This program finds one real root of the equation f(x) = 0 in a finite interval [b,c], where f(x) is a function specified by the user which must be continuous and real on the interval. The program assumes without checking that of the values f(b) and f(c), one will be positive and one negative, i.e., f(b) x f(c) < 0. In this way, b and c will bracket the root. An accuracy tolerance TOL (>=0) must also be specified. This number should be the greatest allowable error in the final approximation for the root. That is, the actual root will be no farther away than TOL from the program's solution for the root.

The function f(x) should be keyed into program memory under LBL E and should assume that x will be in the X-register upon entry. 85 program steps, registers R1 through R7, RS0 through RS9 and the stack are available for defining f(x).

The method used is a combination of bisection (interval-halving) and the secant method. Bisection is often slow but is guaranteed to converge to a root, if one exists in the interval; the secant method is fast but does not always converge. The algorithm employed in this program combines the safety of bisection with some of the speed of the secant method. If the function is known to be well behaved on the interval in question, then the program in Standard Pac, Calculus and Roots of f(x), may lead to a faster and more convenient solution.

Remarks

As each value for b or c is input, its function value will be computed and displayed. If you are in doubt about values for b and c which will satisfy f(b) x f(c) < 0, you may simply keep inputting different values until you strike a good combination. Each new value input overwrites the old.

References

George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler, Computer Methods in Mathematical Computation, Computer Science Department, Stanford University, 1972.

Richard P. Brent, Algorithms for Minimization without Derivatives, Prentice-Hall, 1973.

T. J. Dekker, "Finding a zero by means of successive linear interpolation," in B. Dejon and P. Henrici (editors), Constructive Aspects of the Fundamental Theorem of Algebra, Interscience, 1969.

Instructions

Step

Instructions

Input Data/Units

Keys

Output Data/Units

1

Load side 1 and side 2.

     

2

Prepare to key in function.

 

GTO E

 

3

Switch to PRGM. See line 138.

 

 

 

 4

Key in the function f (x).

 

 

 

5

Switch to RUN.

 

 

 

6

Key in the end points of the interval (remember f(b) x f(c) < 0)

b

A

f(b)

   

c

B

f(c)

7

Key in the accuracy tolerance.

TOL

C

TOL

8

Compute a real root.

 

D

root

9

 To evaluate the function at any point

x

E

f(x)

Examples

Find an angle alpha between 100 and 101 radians such that sin alpha = 0.1. Hence let f(x) = sin x - 0.1. Assume a tolerance of 10-3.

Keystrokes                        Outputs
Load side 1 and side 2.
GTO E
Switch to PRGM. See line 138.
RAD SIN .1 - DEG
Switch to RUN.
l00 A                                -0.61   (f(100))
101 B                                 0.35   (f(101))
EEX CHS 3 C                 1.000000000-03
D                                   100.63   (root) 

Find a root of the equation ln x + 3x - 10.8074 = 0 in the interval [1,5]. An accuracy of 10-4 is acceptable. Store the constant 10.8074 in R1.

Keystrokes                        Outputs
Load side 1 and side 2.
GTO E
Switch to PRGM. See line 138. 
LN LASTx 3 x + RCL 1 -
Switch to RUN.
10.8074 STO 1                        10.81         
1 A                                  -7.81    (f(1))
5 B                                   5.80    (f(5))
EEX CHS 4 C                 1.000000000-04
D                                     3.21    (root)

Check the solution by computing its function value.

E                          -1.901000000-05 (f(3.21))

The Program

LINE  KEYS
001  *LBL A     b
002   STO B
003   GSB E     f(b)
004   STO 8
005   RTN
006  *LBL B
007   STO A     c
008   STO C
009   GSB E
010   STO 0     f(c)
011   STO 9
012   RTN
013  *LBL C     TOL
014   STO E
015   RTN
016  *LBL D
017   RCL 8     If f(b) = 0 exit.
018   X=0?
019   GTO 5
020   RCL B
021   RCL C
022   -
023   ABS
024   RCL E     If TOL > |b-c|, exit.
025   X>Y?
026   GTO 5
027   2
028   ÷
029   EEX
030   CHS
031   9         TOL1 = 10^-9b + .5*TOL.
032   RCL B
033   x
034   +
035   STO I
036   RCL 8
037   RCL 0     Linear interpolation (secant method):
038   RCL 8
039   -
040   RCL A
041   RCL B     b' = b'f(b)/(f(a)-f(b))/(a-b)
042   -
043   ÷
044   ÷         b' may be next b.
045   RCL B
046   X<>Y
047   -
048   STO D     Test b' e[b,c]
049   RCL B
050   -
051   RCL C     s = (b'-b)/(c-b)
052   RCL B
053   -
054   ÷
055   X<0?      If s<0, reject b'.
056   GTO 1
057   1         If s>=1, reject b'.
058   X<=Y?
059   GTO 1
060   RCL B
061   RCL C
062   -
063   ABS
064   4
065   ÷         If b' closer than |b-c|/4 to c
066   RCL D     then reject b'.
067   RCL C
068   -
069   ABS
070   X<=Y?
071   GTO 1
072   RCL I
073   RCL D
074   RCL B
075   -         If |b'-b| <= TOL1, b'<-b+TOL1 x sgn(c-b)
076   ABS
077   X>Y?
078   GTO 2
079   X<>Y
080   RCL C
081   RCL B
082   -
083   ENTER
084   ABS
085   ÷
086   x
087   RCL B
088   +
089   STO D
090   GTO 2
091  *LBL 1
092   RCL B     Reject b', set b'= (b+c)/2, i.e. midpoint of [b,c].
093   RCL C
094   +
095   2
096   ÷
097   STO D
098  *LBL 2     Set new values for next iteration.
099   RCL B
100   STO A     a<-b
101   RCL 8
102   STO 0     f(a)<-f(b)
103   RCL D
104   STO B     b<-b'
105   GSB E
106   STO 8     f(b)<-f(b')
107   RCL 9
108   x
109   X<0?      If f(b) x f(c) < 0, leave c unchanged
110   GTO 3
111   RCL A     Else replace c<-a
112   STO C
113   RCL 0     f(c)<-f(a)
114   STO 9
115  *LBL 3
116   RCL 9
117   ABS
118   RCL 8
119   ABS       If |f(b)| <= |f(c)|, loop again.
120   X<=Y?
121   GTO D
122   RCL B     Else swap b and c and set a<-(new) c.
123   RCL C
124   STO B
125   X<>Y
126   STO C
127   STO A
128   RCL 8
129   RCL 9
130   STO 8
131   X<>Y
132   STO 9
133   STO 0
134   GTO D     Loop again.
135  *LBL 5
136   RCL B     Display root.
137   RTN
138  *LBL E
139   RTN       User defined f(x)

Register Use

R0  f(a)
R8  f(b)
R9  f(c)
A   a
B   b
C   c
D   b'
E   TOL
I   TOL1

Go back to the software library
Go back to the main exhibit hall