HP-71B Program to Find Root

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 a guess for the root. The program uses internally defined values of 1e-8 and 55 for the tolerance and maximum iterations, respectively.

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 the BASIC listing:

OPTION TYPO
OPTION NOLET

! Find root of nonlinear function using Householder's method

! User function
DEF Fx(X) = EXP(X) - 3 * X^2

Declare Numeric Tolerance, X, h, Diff, Iter, MaxIter
Declare Numeric F0, FP1, FM1, FP2, FM2, D1, D2, D3

Tolerance = 1e-8
MaxIter = 55

CLEAR
PRINT "HOUSEHOLDER'S METHOD TO SOLVE FOR THE ROOT OF A NONLINEAR EQUATION"
PRINT

INPUT PROMPT "Enter guess for root: ": X

Iter = 0
Do
  Iter = Iter + 1
  If Iter > MaxIter Then Exit Do
  h = 0.01 * (1 + ABS(X))
  F0 = Fx(X)
  FP1 = Fx(X+h)
  FM1 = Fx(X-h)
  FP2 = Fx(X+2*h)
  FM2 = Fx(X-2*h)
  D1 = (FP1 - FM1) / 2 / h
  D2 = (FP1 - 2 * F0 + FM1) / h^2
  D3 = (FP2 - 2 * FP1 + 2 * FM1 - FM2) / 2 / H ^ 3
  Diff = F0 * (D1^2 - F0 * D2 / 2) / (D1^3 - F0 * D1 * D2 + D3 * F0^2 / 6)
  X = X - Diff
Loop While Abs(Diff) > Tolerance

If Iter > MaxIter Then PRINT "SOLUTION DIVERGED!"
PRINT "Root = ";X
PRINT "Iterations = ";Iter

END

You can customize the above listing by changing the definition of function Fx. 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.

You can also change the tolerance for the root by assigning a different value to the variable Tolerance. You can also change the maximum number of iterations by assigning a different value to the variable MaxIter.

BACK

True BASIC Program to Find The Root Of A Function

by Namir Shammas