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

x' = x – 2 f(x)f'(x)/{2 [f’(x)]^{2} – f(x) f''(x)}

where f(x) is the function whose root is sought, f '(x) is the first derivative of function f(x), and f ''(x) is the second derivative 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}

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

The pseudo-code for the Richmond 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), fp = f(x+h), and fm =
f(x-h)

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

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

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.73307902863 | [CONT] |

ITERS= 5 |

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 200

110 H = 0.01 * (1 + ABS(X))

120 F0 = FNF(X)

130 F1 = FNF(X+H)

140 F2 = FNF(X-H)

150 D1 = (F1 - F2) / 2 / H

160 D2 = (F1 - 2 * F0 + F2) / H^2

170 D = 2 * F0 * D1/(2 * D1^2 - F0 * D2)

180 X = X - D

190 IF ABS(D) > T THEN 80

200 REM STOP

210 BEEP

220 If I > M THEN DISP "SOLUTION FAILED" @ PAUSE

230 DISP "ROOT = ";X

240 PAUSE

250 DISP "ITERS = ";I

260 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-h |

D1 | Value of the first derivative f'(X) |

D2 | Value of the second 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.

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