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)^{2} – f(x) f’’(x) / 2]/[f’(x)^{3} –
f(x) f’(x) f’’(x) – f^{3}(x) f(x)^{2} / 6]}

where f(x) is the function whose root is sought. The terms f '(x),
f''(x), f^{3}(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^{2}

f^{3}(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 *F*x. 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.*

**Copyright (c) Namir Shammas. All rights reserved.**