g-rh@cca.UUCP (04/02/87)
The routine given below is a single file containing a package for
finding a real root of a function, using the secant algorithm.
Only one function may be evaluated at a time, since the secant
algorithm is an iteration with memory (i.e. previous values are
used.) This is not a shar file. There is a cut line at the
beginning and at the end.
-------------- Cut here ---------------
/*
* Copyright (c) 1987 by Richard Harter. Permission to use or
* copy is granted, including commercial distributions, provided
* that this notice is included.
*
* Synopsis:
*
* This is a secant root finder. It searches for a root of a
* real function in the neighbourhood of two initial guesses of
* the root.
*
* Author: Richard Harter
* Date: March, 1987
*
* Entry points:
*
* int firstguess(x,y) presents first guess to algorithm
* double x -- first guess at root
* double y -- value of function at first guess
*
* double nextguess(x,y) returns next guess
* double x -- current guess at root
* double y -- value of function at current guess
* value returned is next guess to process
*
* int rootstatus() Status of iteration
* = -1, root not bracketed
* = 0, root has been found
* = 1, root bracketed, still converging
*
* double rooterror() estimated error of current iterate
* returns -1. if in hunt mode
* returns 0. if root is exact
* returns bracket interval otherwise
*
*
* Usage:
*
* Initialize the algorithm by calling firstguess with an
* approximation to the desired root, and the value of the
* function at that point. Call nextguess with the second
* approximation and the value of the function at that point.
* Nextguess will return the next iterate. Evaluate the
* function at the iterate and repeat in a loop. Terminate
* the loop if (a) function rootstatus returns a 0 or (b)
* function rootstatus returns a -1 and and your iteration
* count limit has been exceeded. Call function rooterror
* to get error interval containing the root.
*
* Note that the best root estimate is the value returned
* by nextguess after rootstatus indicates convergence.
* See the example.
*
* Abstract:
*
* The traditional secant iteration formula is
*
* x[n+1] = x[n] - y[n]*(x[n]-x[n-1])/(y[n]-y[n-1])
*
* with a convergence power of 1.6. The method used here is
* to bracket the root and use the end point with the smaller
* absolute function value rather than the previous iterate.
* This reduces to the secant method when the iterates are
* sufficiently close to the root and is superior when the
* two methods differ. (Note: this is not the same as interpolation
* with the end point with opposite sign which yields geometric
* convergence and can be tediously slow.)
*
* It is not assumed that the first two iterates bracket the
* desired root -- if they do not a hunt mode is initiated.
* The hunt mode is primitive; linear interpolation is used
* for the last two iterates with the delta being overestimated
* by 50% to attempt to force bracketing of the root.
*
* The code contains defensive measures against floating point
* underflow in the delta calculations and against out of range
* extrapolations.
*
*
* Example:
*
* main(){
* double x,y,nextguess(),rooterror(),func();
* int i,rootstatus();
* x = 0.;
* y = func(x);
* firstguess(x,y);
* x = 1.;
* for (i=0;i<25 && rootstatus();i++) {
* y = func(x);
* printf("%20.16f %20.16f\n",x,y);
* x = nextguess(x,y);
* }
* printf("root = %20.16f, interval = %20.16f\n",x,rooterror());
* }
* double func(x) double x; {return x*x -2.;}
*/
#define BEGIN {
#define END }
#define true (1)
#define false (0)
/* Bracket values when root is bracketed */
static double xl; /* Low bracket value for x */
static double yl; /* Func value at xlo */
static double xh; /* High bracket value for x */
static double yh; /* Func value at xhi */
/* Previous iterates */
static double x1; /* Previous x value */
static double y1; /* Previous y value */
/* Logic control */
static int mode; /* Processing case */
/* 0: Normal second guess */
/* 1: Hunt mode */
/* -1: Normal bracket mode */
static int found; /* True is root already found */
static int increasing; /* True is increasing func */
static int exact; /* Root is exact */
static double xf; /* Value of found root */
firstguess(x,y) /* Processes first guess */
double x; /* Argument value for 1st guess */
double y; /* Function value for 1st guess */
BEGIN
x1 = x; /* Set previous arg value */
y1 = y; /* Set previous func value */
mode = 0; /* Initialize mode */
if (y==0.0) BEGIN /* First guess is good */
xf = x; /* Set found arg */
found = true; /* Set flag */
exact = true; /* Set exact root flag */
END
else BEGIN /* Not lucky on first guess */
found = false; /* Not found, say so */
exact = false; /* Not an exact root */
END
END
double nextguess(x,y) /* Yields nextguess from current*/
double x; /* Current arg value */
double y; /* current func value */
BEGIN
double xret; /* Next arg value to try */
double delta; /* Change in iterate */
double tdelta; /* L.B. for delta in y= case */
double setfinal(); /* Sets final x on convergence */
int uselow; /* Flag signifying smallest y */
if (found) return xf; /* Root already found */
if (y==0.0) BEGIN /* Root has just been found */
xf = x; /* Set found value */
found = true; /* Set flag */
exact = true; /* Set exact flag */
return xf; /* Return root */
END
if (!mode) BEGIN /* This is 2nd guess, startup */
if (y<0 && y1>0) mode = -1; /* Bracket established */
else if (y>0 && y1<0) mode = -1; /* Likewise */
else mode = 1; /* Must be hunt mode */
if (mode<0) BEGIN /* Bracketed, load bracket */
if (x<x1) BEGIN /* Order is {x,x1} */
xl = x; xh = x1; /* Load x values */
yl = y; yh = y1; /* Load y values */
END
else BEGIN /* Order is {x1,x} */
xl = x1; xh = x; /* Load x values */
yl = y1; yh = y; /* Load y values */
END
if (yl<yh) increasing = true; /* Function is increasing */
else increasing = false; /* Function is decreasing */
xret = xl - yl*(xh-xl)/(yh-yl); /* Straight interpolation */
return xret; /* Return it */
END
END
if (mode<0) BEGIN /* Normal bracket mode */
if (x<xl || x>xh) BEGIN /* Arg is out of range */
xret = (xl + xh)/2.0; /* Get bracket midrange */
return xret; /* Return it */
END
if (increasing) BEGIN /* Cose to find smallest y */
if (yh > -yl) uselow = true; /* Smallest is |yl| */
else uselow = false; /* Smallest is |yh| */
END
else BEGIN /* Function is decreasing */
if (yl < -yh) uselow = true; /* Smallest is |yl| */
else uselow = false; /* Smallest is |yh| */
END
if (uselow) BEGIN /* Use low end x,y */
if (y==yl) BEGIN /* y=yl is special case */
delta = 2.0*(x-xl); /* Explode small end delta */
tdelta= .1*(xh-xl); /* Lower bound for delta size */
if (delta<tdelta) delta=tdelta; /* Use L.B. if needful */
END
else delta = -y*(x-xl)/(y-yl); /* Normal secant delta */
xret = x + delta; /* Compute expected next */
if (delta>0.0) BEGIN /* Increasing delta */
if (xret>xh) BEGIN /* Next x out of range */
delta = (xh-x)/2.0; /* Split the range */
xret = x + delta; /* Recompute next x */
END
while (xret<=x) BEGIN /* Loop to catch underflow */
delta *= 2.0; /* Double delta to try again */
xret = x + delta; /* Recompute next x */
END
if (xret>=xh) BEGIN /* Can't make good delta */
found = true; /* Must be final iterate */
xf = setfinal(x,y); /* Pick best final iterate */
xret = xf; /* Return found root */
END
END
else if (delta<0.0) BEGIN /* Decreasing delta */
if (xret<xl) BEGIN /* Next x out of range */
delta = (xl-x)/2.0; /* Split the range */
xret = x + delta; /* Recompute next x */
END
while (xret>=x) BEGIN /* Loop to catch underflow */
delta *= 2.0; /* Double delta to try again */
xret = x + delta; /* Recompute next x */
END
if (xret<=xl) BEGIN /* Can't make good delta */
found = true; /* Must be final iterate */
xf = setfinal(x,y); /* Pick best final iterate */
xret = xf; /* Return found root */
END
END
else BEGIN /* Delta exactly 0.0, done */
xf = x; /* Accept root */
found = true; /* Record as found */
xret = x; /* Return the root */
END
END
else BEGIN /* Use high end x,y */
if (y==yh) BEGIN /* y=yh, special case */
delta = 2.*(x-xh); /* Explode small delta */
tdelta = .1*(xl-x); /* L.B. for delta */
if (delta>tdelta) delta=tdelta; /* Use LB if needful */
END
else delta = -y*(xh-x)/(yh-y); /* Normal secant delta */
xret = x + delta; /* Compute expected next */
if (delta>0.0) BEGIN /* Increasing delta */
if (xret>xh) BEGIN /* Next x out of range */
delta = (xh-x)/2.0; /* Split the range */
xret = x + delta; /* Recompute next x */
END
while (xret<=x) BEGIN /* Loop to catch underflow */
delta *= 2.0; /* Double delta to try again */
xret = x + delta; /* Recompute next x */
END
if (xret>=xh) BEGIN /* Can't make good delta */
found = true; /* Must be final iterate */
xf = setfinal(x,y); /* Pick best final iterate */
xret = x; /* Return found root */
END
END
else if (delta<0.0) BEGIN /* Decreasing delta */
if (xret<xl) BEGIN /* Next x out of range */
delta = (xl-x)/2.0; /* Split the range */
xret = x + delta; /* Recompute next x */
END
while (xret>=x) BEGIN /* Loop to catch underflow */
delta *= 2.0; /* Double delta to try again */
xret = x + delta; /* Recompute next x */
END
if (xret<=xl) BEGIN /* Can't make good delta */
found = true; /* Must be final iterate */
xf = setfinal(x,y); /* Pick best final iterate */
xret = x; /* Return found root */
END
END
else BEGIN /* Delta exactly 0.0, done */
xf = x; /* Accept root */
found = true; /* Record as found */
xret = x; /* Return the root */
END
END
if ((increasing && y>0.0) || (!increasing && y<0.0)) BEGIN
xh = x; /* Move high end down */
yh = y; /* Same for func value */
END
else BEGIN /* Moving low end forward */
xl = x; /* Moving x */
yl = y; /* Moving y */
END
return xret; /* Return accepted root */
END
if (mode>0) BEGIN /* Hunt mode */
if ((y<0. && y1>0.) || (y>0. && y1<0.)) BEGIN
mode = 0; /* Reset as though 2nd guess */
return nextguess(x,y); /* Reuse switch to bracket code */
END
xret = x - 1.5*y*(x-x1)/(y-y1); /* Overestimate delta */
x1 = x; /* Push x */
y1 = y; /* Push y */
END
return xret; /* Return next guess */
END
double setfinal(x,y) /* Picks best final iterate */
double x; /* Last x given us */
double y; /* Last y given us */
BEGIN
if (increasing) BEGIN /* Function is increasing */
if (y>0.) BEGIN /* Choose between xl and x */
if (y > -yl) return xl; /* |y|>|yl| use xl */
else return x; /* |yl|>=|y| use x */
END
else BEGIN /* Choose between x and xh */
if (-y > yh) return xh; /* |y|>|yh| use xh */
else return x; /* |yh|>=|y| use x */
END
END
else BEGIN /* Function is decreasing */
if (y>0.) BEGIN /* Choose between x and xh */
if (y > -yh) return xh; /* |y|>|yh| use xh */
else return x; /* |y|<=|yh| use x */
END
else BEGIN /* Choose between x and xl */
if (-y > yl) return xl; /* |y|>|yl| use xl */
else return x; /* |y|<=|yl| use x */
END
END
END
int rootstatus() BEGIN /* Returns status of iteration */
if (found) return 0; /* 0 is root is found */
if (mode<1) return 1; /* 1 is bracketed, converging */
return -1; /* Still in hunt mode */
END
double rooterror() BEGIN /* Returns root interval */
if (exact) return 0.; /* Root is exact */
if (mode>0) return -1.; /* Still in hunt mode */
return (xh-xl); /* Return bracket size */
END
----------- Cut here __________________
--
Richard Harter, SMDS Inc. [Disclaimers not permitted by company policy.]