*This class was translated by a statistician from the FORTRAN *version of dfzero. It is NOT an official translation. When *public domain Java optimization routines become available *from professional numerical analysts, then THE CODE PRODUCED *BY THE NUMERICAL ANALYSTS SHOULD BE USED. * *
*Meanwhile, if you have suggestions for improving this *code, please contact Steve Verrill at steve@www1.fpl.fs.fed.us. * *@author Steve Verrill *@version .5 --- April 17, 2001 * */ public class Fzero extends Object { /** * *
*This method searches for a zero of a function f(x) between *the given values b and c until the width of the interval *(b,c) has collapsed to within a tolerance specified by *the stopping criterion, Math.abs(b-c) <= 2.0*(rw*Math.abs(b)+ae). *The method used is an efficient combination of bisection *and the secant rule. *The introductory comments from *the FORTRAN version are provided below. * *This method is a translation from FORTRAN to Java of the Netlib *(actually it is in the SLATEC library which is available at Netlib) *function dfzero. In the FORTRAN code, L.F. Shampine and H.A. Watts *are listed as the authors. * *Translated by Steve Verrill, April 17, 2001. * *@param zeroclass A class that defines a method, f_to_zero, * that returns a value that is to be zeroed. * The class must implement * the Fzero_methods interface (see the definition * in Fzero_methods.java). See FzeroTest.java * for an example of such a class. * f_to_zero must have one * double valued argument. *@param b[1] One endpoint of the initial interval. The value returned * for b[1] is usually the better approximation to a * zero of f_to_zero. *@param c[1] The other endpoint of the initial interval. *@param r Initial guess of a zero. A (better) guess of a zero * of f_to_zero which could help in * speeding up convergence. If f_to_zero(b) and f_to_zero(r) have * opposite signs, a root will be found in the interval * (b,r); if not, but f_to_zero(r) and f_to_zero(c) have opposite * signs, a root will be found in the interval (r,c); * otherwise, the interval (b,c) will be searched for a * possible root. When no better guess is known, it is * recommended that r be set to b or c; because if r is * not interior to the interval (b,c), it will be ignored. * *@param re Relative error used for rw in the stopping criterion. * If the requested re is less than machine precision, * then rw is set to approximately machine precision. * *@param ae Absolute error used in the stopping criterion. If the * given interval (b,c) contains the origin, then a * nonzero value should be chosen for AE. * *@param iflag[1] A status code. User must check iflag[1] after each call. * Control returns to the user from dfzero in all cases. * * 1: b is within the requested tolerance of a zero. * The interval (b,c) collapsed to the requested * tolerance, the function changes sign in (b,c), and * f_to_zero(x) decreased in magnitude as (b,c) collapsed. * * 2: f_to_zero(b) = 0. However, the interval (b,c) may not have * collapsed to the requested tolerance. * * 3: b may be near a singular point of f_to_zero(x). * The interval (b,c) collapsed to the requested tol- * erance and the function changes sign in (b,c), but * f_to_zero(x) increased in magnitude as (b,c) collapsed,i.e. * Math.abs(f_to_zero(b out)) > * Math.max(Math.abs(f_to_zero(b in)), * Math.abs(f_to_zero(c in))) * * 4: No change in sign of f_to_zero(x) was found although the * interval (b,c) collapsed to the requested tolerance. * The user must examine this case and decide whether * b is near a local minimum of f_to_zero(x), or b is near a * zero of even multiplicity, or neither of these. * * 5: Too many (> 500) function evaluations used. * */ public static void fzero (Fzero_methods zeroclass, double b[], double c[], double r, double re, double ae, int iflag[]) { /* Here is a copy of the Netlib documentation: SUBROUTINE DFZERO(F,B,C,R,RE,AE,IFLAG) C***BEGIN PROLOGUE DFZERO C***DATE WRITTEN 700901 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. F1B C***KEYWORDS LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(FZERO-S DFZERO-D), C BISECTION,NONLINEAR,NONLINEAR EQUATIONS,ROOTS,ZEROES, C ZEROS C***AUTHOR SHAMPINE,L.F.,SNLA C WATTS,H.A.,SNLA C***PURPOSE Search for a zero of a function F(X) in a given C interval (B,C). It is designed primarily for problems C where F(B) and F(C) have opposite signs. C***DESCRIPTION C C **** Double Precision version of FZERO **** C C Based on a method by T J Dekker C written by L F Shampine and H A Watts C C DFZERO searches for a zero of a function F(X) between C the given values B and C until the width of the interval C (B,C) has collapsed to within a tolerance specified by C the stopping criterion, DABS(B-C) .LE. 2.*(RW*DABS(B)+AE). C The method used is an efficient combination of bisection C and the secant rule. C C Description Of Arguments C C F,B,C,R,RE and AE are DOUBLE PRECISION input parameters. C B and C are DOUBLE PRECISION output parameters and IFLAG (flagged C by an * below). C C F - Name of the DOUBLE PRECISION valued external function. C This name must be in an EXTERNAL statement in the C calling program. F must be a function of one double C precision argument. C C *B - One end of the interval (B,C). The value returned for C B usually is the better approximation to a zero of F. C C *C - The other end of the interval (B,C) C C R - A (better) guess of a zero of F which could help in C speeding up convergence. If F(B) and F(R) have C opposite signs, a root will be found in the interval C (B,R); if not, but F(R) and F(C) have opposite C signs, a root will be found in the interval (R,C); C otherwise, the interval (B,C) will be searched for a C possible root. When no better guess is known, it is C recommended that r be set to B or C; because if R is C not interior to the interval (B,C), it will be ignored. C C RE - Relative error used for RW in the stopping criterion. C If the requested RE is less than machine precision, C then RW is set to approximately machine precision. C C AE - Absolute error used in the stopping criterion. If the C given interval (B,C) contains the origin, then a C nonzero value should be chosen for AE. C C *IFLAG - A status code. User must check IFLAG after each call. C Control returns to the user from FZERO in all cases. C XERROR does not process diagnostics in these cases. C C 1 B is within the requested tolerance of a zero. C The interval (B,C) collapsed to the requested C tolerance, the function changes sign in (B,C), and C F(X) decreased in magnitude as (B,C) collapsed. C C 2 F(B) = 0. However, the interval (B,C) may not have C collapsed to the requested tolerance. C C 3 B may be near a singular point of F(X). C The interval (B,C) collapsed to the requested tol- C erance and the function changes sign in (B,C), but C F(X) increased in magnitude as (B,C) collapsed,i.e. C abs(F(B out)) .GT. max(abs(F(B in)),abs(F(C in))) C C 4 No change in sign of F(X) was found although the C interval (B,C) collapsed to the requested tolerance. C The user must examine this case and decide whether C B is near a local minimum of F(X), or B is near a C zero of even multiplicity, or neither of these. C C 5 Too many (.GT. 500) function evaluations used. C***REFERENCES L. F. SHAMPINE AND H. A. WATTS, *FZERO, A ROOT-SOLVING C CODE*, SC-TM-70-631, SEPTEMBER 1970. C T. J. DEKKER, *FINDING A ZERO BY MEANS OF SUCCESSIVE C LINEAR INTERPOLATION*, 'CONSTRUCTIVE ASPECTS OF THE C FUNDAMENTAL THEOREM OF ALGEBRA', EDITED BY B. DEJON C P. HENRICI, 1969. C***ROUTINES CALLED D1MACH C***END PROLOGUE DFZERO c */ int i,ic,icnt,ierr,ipass,ipss,j,klus,kount,kprint,lun,ndeg; int itest[] = new int[37]; int itmp[] = new int[16]; double a,acbs,acmb,aw,cmb,er,fa,fb,fc,fx,fz,p,q,rel,rw, t,tol,wi,work,wr,z; // 1.1102e-16 is machine precision er = 2.0*1.2e-16; // Initialize z = r; if (r <= Math.min(b[1],c[1]) || r >= Math.max(b[1],c[1])) { z = c[1]; } rw = Math.max(re,er); aw = Math.max(ae,0.0); ic = 0; t = z; fz = zeroclass.f_to_zero(t); fc = fz; t = b[1]; fb = zeroclass.f_to_zero(t); kount = 2; if (Blas_f77.sign_f77(1.0,fz) != Blas_f77.sign_f77(1.0,fb)) { c[1] = z; } else { if (z != c[1]) { t = c[1]; fc = zeroclass.f_to_zero(t); kount = 3; if (Blas_f77.sign_f77(1.0,fz) != Blas_f77.sign_f77(1.0,fc)) { b[1] = z; fb = fz; } } } a = c[1]; fa = fc; acbs = Math.abs(b[1] - c[1]); fx = Math.max(Math.abs(fb),Math.abs(fc)); // Main loop while (true) { if (Math.abs(fc) < Math.abs(fb)) { // Perform interchange a = b[1]; fa = fb; b[1] = c[1]; fb = fc; c[1] = a; fc = fa; } cmb = 0.5*(c[1] - b[1]); acmb = Math.abs(cmb); tol = rw*Math.abs(b[1]) + aw; // Test stopping criterion and function count. if (acmb <= tol) { // Finished. Process the results for the proper setting // of the flag. if (Blas_f77.sign_f77(1.0,fb) == Blas_f77.sign_f77(1.0,fc)) { iflag[1] = 4; return; } if (Math.abs(fb) > fx) { iflag[1] = 3; return; } iflag[1] = 1; return; } if (fb == 0.0) { iflag[1] = 2; return; } if (kount >= 500) { iflag[1] = 5; return; } /* C CALCULATE NEW ITERATE IMPLICITLY AS B+P/Q C WHERE WE ARRANGE P .GE. 0. C THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW. */ p = (b[1] - a)*fb; q = fa - fb; if (p < 0.0) { p = -p; q = -q; } /* C UPDATE A AND CHECK FOR SATISFACTORY REDUCTION C IN THE SIZE OF THE BRACKETING INTERVAL. C IF NOT, PERFORM BISECTION. */ a = b[1]; fa = fb; ic++; if (ic >= 4 && 8.0*acmb >= acbs) { // Use bisection b[1] = .5*(c[1] + b[1]); } else { if (ic >= 4) { ic = 0; acbs = acmb; } // Test for too small a change. if (p <= Math.abs(q)*tol) { // Increment by the tolerance. b[1] += Blas_f77.sign_f77(tol,cmb); } else { // Root ought to be between b[1] and (c[1] + b[1])/2. if (p >= cmb*q) { // Use bisection. b[1] = 0.5*(c[1] + b[1]); } else { // Use the secant rule. b[1] += p/q; } } } // Have completed the computation for a new iterate b[1]. t = b[1]; fb = zeroclass.f_to_zero(t); kount++; // Decide whether the next step is an interpolation // or an extrapolation. if (Blas_f77.sign_f77(1.0,fb) == Blas_f77.sign_f77(1.0,fc)) { c[1] = a; fc = fa; } } } }