San code – the secant root finderThe 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>>
end
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 CODE FOR 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. |