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:

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^{3}(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.

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