/* ---------------------- MODULE fpegasus.cpp ----------------------- */ #include #define ITERMAX 300 /* Maximal number of iterations */ #define ABSERR ZERO /* Admissable absolute error */ #define RELERR ((REAL)4.0 * MACH_EPS) /* Admissable relative error */ #define FCTERR EPSQUAD /* Maximal function errror */ int pegasus /* Pegasus Method .........................*/ ( REALFCT fct, /* Function ........................*/ REAL * x1, /* Starting value 1 ................*/ REAL * x2, /* Starting value 2 / solution .....*/ REAL * f2, /* Function value at x2 ............*/ int * iter /* Number of iterations ............*/ ) /*====================================================================* * * * pegasus computes one zero of the continuous function fct, * * provided that the two starting values x1 and x2 satisfy: * * fct(x1) * fct(x2) <= 0 . * * * *====================================================================* * * * Applications: * * ============ * * Determine one root of the continuous function fct, if an * * inclusion interval [x1, x2] is known for the root. * * * *====================================================================* * * * Input parameters: * * ================ * * fct REAL fct (REAL); * * Function, whose root is to be found. * * fct has the form: * * REAL fct (REAL x) * * { REAL f; * * f = ...; * * return(f); * * } * * x1,x2 REAL *x1, *x2; * * Starting values with fct(x1) * fct(x2) <= 0. * * * * Output parameters: * * ================== * * x2 REAL *x2; * * Computed approximation for the root of fct * * * * f2 REAL *f2; * * Functional value at the approximate root, this * * must be nearly zero. * * iter int *iter; * * Number of iterations that were performed. * * * * Return values: * * ============= * * = -1 No inclusion: fct(x2) * fct(x1) > 0 * * = 0 Root has been found with ABS(f2) < FCTERR * * = 1 Break-off with * * ABS(xnew-xold) < ABSERR + xnew * RELERR, * * check functional value * * = 2 Iteration limit reached * * = 3 Improper input parameters * * * *====================================================================* * * * Constants used: ABSERR, RELERR, MACH_EPS, EPSROOT, ITERMAX * * * * Reference: [BIBLI 11] * *====================================================================*/ { REAL f1, x3, f3, s12; int rc = 2; *iter = 0; /* Initialize iteration counter */ if (fct == NULL) return (3); f1 = (*fct)(*x1); /* Function values at *x1, *x2 */ *f2 = (*fct)(*x2); if (f1 * (*f2) > ZERO) return (-1); /* No inclusion -> Error */ if (f1 * (*f2) == ZERO) /* One starting value is a root */ { if (f1 == ZERO) { *x2 = *x1; *f2 = ZERO; } return (0); } while (*iter <= ITERMAX) /* Pegasus iteration */ { (*iter)++; s12 = (*f2 - f1) / (*x2 - *x1); /* Secant slope */ x3 = *x2 - *f2 / s12; /* new approximation */ f3 = (*fct)(x3); if (*f2 * f3 <= ZERO) /* new inclusion interval */ { *x1 = *x2; f1 = *f2; } else f1 *= *f2 / ( *f2 + f3 ); *x2 = x3; *f2 = f3; if ( ABS(*f2) < FCTERR ) /* Root found */ { rc = 0; break; } /* Break-off with small step size */ if ( ABS(*x2 - *x1) <= ABS(*x2) * RELERR + ABSERR ) { rc = 1; break; } } if ( ABS(f1) < ABS(*f2) ) /* Choose approximate root with */ { /* least magnitude function value */ *x2 = *x1; *f2 = f1; } return (rc); } /* ---------------------- END fpegasus.cpp ------------------------ */