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.]