table of contents
Comp Sci
October 2003

San code - the secant root finder

The following example is the code for a secant root finder routine. The secant root finder algorithm is an algorithm for finding a zero of a function, i.e., given a function f(x), find an x such that f(x) = 0.

The algorithm is initialized with two evaluated points, y0 = f(x0), and y1 = f(x1). Succesive approximations to the root are calculated using the iteration equation:

        x_n+1 = x_n - y_n*(x_n - x_n-1)/(y_n - y_n-1)
This iteration can be transformed into
        beta      = y_n/y_n-1
        delta_n+1 = - delta_n * beta / (1 - beta)
        x_n+1     = x_n + delta_n+1
There are a number of considerations that one would want to attend to in a production implementation. The following implementation handles the core algorithm.

The idea is to define a morph called secant that has defined within it two functions, secant$f and secant$init. The former is the iteration routine, and the latter is the initialization routine. Within a morph all fields using the $ qualifier are values, and all routines within the morph have access to the morph's fields. The usage would run something like this:

    begin loop while secant$delta gt? epsilon
        init x := <secant$init x0, <f x0>, x1, <f x1>>
        x := <secant <f x>>
Note that secant$f can be replaced by secant because $f is the default function for a morph. The various fields in secant, $delta, $x, and $y persist as long as the morph does.

The reader might ask, suppose we have a routine in which we need to concurrently find zeros of two functions. The answer is that we can duplicate (clone) the morph, e.g.,

        set secant1 = secant
        set secant2 = secant
begin morph secant
        begin proc $f
                args: y
                beta := y / $y
                $delta := - $delta * beta / (1 - beta)
                $x += $delta
                $y := y
                return $x
                end $f
        begin proc $init
                args: x0, y0, x1, y1
                beta   := y1 / y0
                delta  := x1 - x0
                $delta := - delta * beta / (1 - beta)
                $y     := y1
                $x     := x1 + $delta
                return $x
                end $init
        end secant                

This page was last updated October 3, 2003.

table of contents
Comp Sci
October 2003