The Museum of HP Calculators The Museum of HP Calculators
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.
| Solution to f(x) = 0 | |||||
| Shift | |||||
| Label | b→f(b) | c→f(c) | TOL | →root | x→f(x) |
| Key | A | B | C | D | E |
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.
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.
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.
| 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) |
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))
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 × 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 × 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 093 RCL C of [b,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 × 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)
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