GENERIC INTERFACE FindZero(R);

Arithmetic for Modula-3, see doc for details

Abstract: Root finders

FROM Arithmetic IMPORT Error;

TYPE
  Ftn = PROCEDURE (x: R.T; ): R.T;

  DerivativeArray = ARRAY OF R.T;
  DerivativeArray2 = ARRAY [0 .. 1] OF R.T;
  DifFtn = PROCEDURE (x: R.T; ): DerivativeArray2;

  Bracket = RECORD l, r: R.T END;

PROCEDURE AreBracketing (y1, y2: R.T; ): BOOLEAN;

Test if two function values have different signs. If they have and the function is continuous then between the function points is a zero.

PROCEDURE BracketOut
  (func: Ftn;                    (* find brackets for this function *)
   VAR (*INOUT*) x      : Bracket;  (* starting with these points *)
                 maxiter: CARDINAL  := 55 (* growing maxiter times *)
  ): BOOLEAN;                    (* true if successful *)

Given x.l,x.r, search for points (returned in x.l, x.r) for which func(x.l) is opposite sign from func(x.r). Grow outward from the original x.l,x.r by golden ratio, for geometric growth. Return true if a good x.l,x.r can be found before getting to maxiter, else return false.

requires: x.l<x.r.

PROCEDURE BracketIn
  (func: Ftn;                    (* find brackets for this function *)
   READONLY x: Bracket;          (* starting with these points *)
   n: [1 .. LAST(CARDINAL)];     (* using n equi-sized segments *)
  ): REF ARRAY OF Bracket;       (* returning pairs here *)

Break up the x.l..x.r range into n equi-sized segments. Select all pairs which allow bracketing.

requires: x.l<x.r.

PROCEDURE Bisection (func: Ftn;  (* find root of this function *)
                     READONLY x  : Bracket;  (* between these brackets *)
                              tol: R.T;      (* to within +/- tolerance *)
                     maxiter := 45 (* but no more than maxiter cuts *)
  ): R.T RAISES {Error};         (* returning the root *)

Given brackets x.l,x.r, find a root via bisection, and refine it to within +/- tol

PROCEDURE Brent (func: Ftn;      (* find a root of this function *)
                 READONLY x  : Bracket;  (* between these bracket points *)
                          tol: R.T;      (* to this tolerance *)
                 maxiter := 100  (* with <= maxiter iterations *)
  ): R.T RAISES {Error};

Use Brent's algorithm to find the real root between the bracket points. x.l and x.r must be of opposite signs.

PROCEDURE NewtonRaphson
  (         func   : DifFtn;     (* this ftn *)
   READONLY x      : Bracket;    (* bracketed by these points *)
            xtol   : R.T;        (* find root to this precision of x *)
            maxiter            := 25 (* with no more than maxiter loops *)
  ): R.T RAISES {Error};         (* returning root *)

Given a function which returns both f(x) and df(x), and brackets x.l and x.r, find the root to xtol precision. Works via newton-raphson and bisection.

END FindZero.