HP-71B Program to Find The Root Of A Function

by Namir Shammas

The following program calculates the root of a function using Householder's algorithm. This algorithm uses the following equation to update the guess for the root:

x = x – f(x){[f’(x)² – f(x) f’’(x) / 2]/[f’(x)³ – f(x) f’(x) f’’(x) – f³(x) f(x)² / 6]}

where f(x) is the function whose root is sought. The terms f '(x), f''(x), f³(x) are the first, second, and third derivatives of function f(x). The program approximate the derivatives using the following difference approximations:

f '(x) ≈ (f(x + h) - f(x-h)) / 2h

f''(x) ≈ (f(x+h) - 2 f(x) + f(x-h))/ h²

f³(x) ≈ (f(x+2h) - 2 f(x+h) + 2 f(x-h) - f(x-2h)) / (2 h ^ 3)

where h = 0.01 * (1 + |x|)

The pseudo-code for the Householder algorithm is:

Obtain initial guess x and Tolerance

Obtain or set value to MaxIterations

Clear Diverge flag

Set IterCount = 0

Loop

Increment IterCount

If IterCount > MaxIterations Then Set Diverge flag and exit loop

h = 0.01 * [1 + |x|]

Let f0 = f(x), fp1 = f(x+h), fp2 = f(x+2h), fm1 = f(x-h), and fm2 = f(x-2h)

Let D1 = (fp – fm) / (2h), and D2 = (fp – 2f0 + fm) / (h*h)

D3 = (fp2 - 2 * fp1 + 2 * fm1 - fm2) / 2 / h ^ 3

diff = f0 * (f0^2 – f0 * D1 / 2)/(f0^3 – f0 * D1 * D2 – D3 * f0^2 / 6)

x = x – diff

Until |diff| < Tolerance

Return root as x and Diverge flag

The program prompts you to enter:

1. Guess for the root.

2. Tolerance for the root. The default value is 1E-7.

3. The maximum number of iterations. The default value is 55.

The program displays the following results:

1. The root value.

2. The number of iterations.

If the number of iterations exceeds the maximum limit, the program displays the text SOLUTION FAILED before displaying the above results.

Here is a sample session to find the root of f(x) = e^x - 3*x^2 near x = 5 and using the default values for the tolerance and maximum number of iterations:

PROMPT/DISPLAY	ENTER/PRESS
>	[RUN]
GUESS?	5[END LINE]
TOLER? 1E-7	[END LINE]
MAX ITERS? 55	[END LINE]
(Audio beep)
ROOT= 3.73307902872	[CONT]
ITERS= 4

Here is the BASIC listing:

10 DEF FNF(X) = EXP(X) - 3 * X^2
20 INPUT "GUESS? ";X
30 INPUT "TOLER? ","1E-7";A$
40 T = VAL(A$)
50 INPUT "MAX ITER? ","55";A$
60 M = VAL(A$)
70 I = 0
80 REM START
90 I = I + 1
100 IF I > M THEN 230
110 H = 0.01 * (1 + ABS(X))
120 F0 = FNF(X)
130 F1 = FNF(X+H)
140 F2 = FNF(X+2*H)
150 F3 = FNF(X-H)
160 F4 = FNF(X-2*H)
170 D1 = (F1 - F3) / 2 / H
180 D2 = (F1 - 2 * F0 + F3) / H^2
190 D3 = (F2 - 2 * F1 + 2 * F3 - F4) / 2 / H ^ 3
200 D = F0 * (D1 ^ 2 - F0 * D2 / 2) / (D1 ^ 3 - F0 * D1 * D2 + D3 * F0 ^ 2 / 6)
210 X = X - D
220 IF ABS(D) > T THEN 80
230 REM STOP
240 BEEP
250 If I > M THEN DISP "SOLUTION FAILED" @ PAUSE
260 DISP "ROOT = ";X
270 PAUSE
280 DISP "ITERS = ";I
290 END

The program uses the variables shown in the following table:

Variable Name	Contents
X	Guess for root
T	Tolerance
M	Maximum number of iterations
I	Iteration counter
H	Increment h
F0	Value of the function at X
F1	Value of the function at X+h
F2	Value of the function at X+2h
F3	Value of the function at X-h
F4	Value of the function at X-2h
D1	Value of the first derivative f'(X)
D2	Value of the second derivative f''(X)
D3	Value of the third derivative f³(X)
D	Root refinement
A$	Temporary input for T and M

You can customize the above listing by changing the definition of function FNF in line 10. The current listing is set to solve the following equation:

f(x) = e^x - 3 * x^ 2

The above equation has roots near -0.45, 0.91, and 3.73.

BACK