subroutine aptnewt (a, np, xinit, xmin, xmax, nitmax, tol, & xroot, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTNEWT c c call aptnewt (a, np, xinit, xmin, xmax, nitmax, tol, xroot, nerr) c c Version: aptnewt Updated 1997 September 3 13:20. c aptnewt Originated 1997 July 30 15:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find a root xroot of the polynomial equation in x, c P(x) = SUM {a(n)*x**n}, (n = 0, np), with an initial guess of c xinit, a search range from xmin to xmax, allowing no more then c nitmax Newtonian iterations, and a relative precision of tol. c Flag nerr, if non-zero, indicates an input error, or a failure c to find a root. c c Method: A modified Newtonian iteration method is used, with an initial c x of xinit. The n'th iterate of x is found as follows: c x(n) = x(n-1) - P{x(n-1)} / P'{x(n-1)}, where P'(x) is the c derivative of P(x) with respect to x. A root has been found c when the relative error in calculating P(x(n)) is no more than c tol. A maximum of nitmax iterations are allowed. If P'(x) is c zero, or if x gets stuck at xmin or xmax, a new trial value of c x is randomly selected between xmin and xmax. c c Note: Roots at a minimum or maximum of P(x) can not be found by this c method. Try finding the roots of P'(x), and seeing if P(x) is c zero at any of those roots. Repeat for higher derivatives, c if necessary. c If an interval (xmin, xmax) has been found for which P(xmin) c and P(xmax) have opposite signs, a good initial guess xinit is c (xmin * P(xmax) - xmax * P(xmin)) / (P(xmax) - P(xmin)). c Up to np roots may exist between xmin and xmax. c c Input: a, np, xinit, xmin, xmax, nitmax, tol. c c Output: xroot, nerr c c Glossary: c c a Input The array of coefficients of the polynomial equation c in x of degree np: c P(x) = SUM {a(n) * x**n} (n = 0, np). c Must be dimensioned a(0:np) or bigger. c c nitmax Input The maximum number of iterations to be tried. c Should be at least -3.32 log (tol), or about 50 c for 64-bit floating point numbers. c c nerr Output Indicates no input errors were found, and a root was c found, if zero. Otherwise: c 1 np was not positive. c 2 xmin was not less than xmax. c 3 xinit was not between xmin and xmax. c 4 nitmax was not positive. c 5 no convergence after nitmax iterations. c c np Input The largest exponent of x in the polynomial equation c P(x), for which the coefficient is non-zero. c Must be positive. If 4, use aptquar. If 3, use c aptcubs. If 2, use aptqrts or aptqrtv. c If 1, xroot = -a(0) / a(1). c c tol Input Numerical iteration convergence control. c If zero, a value of 1.e-11 is used instead. c If using 64-bit floating point numbers, a value c between 1.e-11 and 1.e-14 is recommended. c c xinit Input The initial value of x to try. Must be between c xmin and xmax. c c xmax Input The maximum allowed value of x to try. c Must be more than xmin. c c xmin Input The minimum allowed value of x to try. c Must be less than xmax. c c xroot Output A root of the polynomial P(x), for which the following c conditions are satisfied: c The initial guess is xinit; c the root is between xmin and xmax; c no more than nitmax iterations are required; c the relative error in calculating P(xroot) is no c more than tol. c If no root was found, -1.e99 is returned. c Up to np roots may exist between xmin and xmax. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. dimension a (0:1) ! The coefficients of x**n. cbugc***DEBUG begins. cbug 9901 format (/ 'aptnewt using Newton''s iteration to find a root', cbug & ' of a polynomial.') cbug 9902 format ('tol=',1pe22.14,' np=',i5,' nitmax=',i6) cbug 9903 format ('xinit,min,max=',1p3e22.14) cbug 9904 format (i5,2x,1pe22.14) cbug write ( 3, 9901) cbug write ( 3, 9902) tol, np, nitmax cbug write ( 3, 9903) xinit, xmin, xmax cbug write ( 3, 9904) (n, a(n), n = 0, np) cbugc***DEBUG ends. c.... Initialize. nerr = 0 xroot = -1.e99 niter = 0 tolx = tol if (tolx .eq. 0.0) tolx = 1.e-11 c.... Test for input errors. if (np .le. 0) then nerr = 1 go to 210 endif xerr = 2.0 * tolx * (abs (xmax) + abs (xmin)) if ((xmin - xmax) .gt. xerr) then nerr = 2 go to 210 endif xerr1 = 2.0 * tolx * (abs (xinit) + abs (xmin)) xerr2 = 2.0 * tolx * (abs (xinit) + abs (xmax)) if (((xmin - xinit) .gt. xerr1) .or. & ((xinit - xmax) .gt. xerr2)) then nerr = 3 go to 210 endif if (nitmax .le. 0) then nerr = 4 go to 210 endif c.... Initialize iterate. xold = xinit xnew = xinit c.... Iterate using Newtonian iteration. 110 niter = niter + 1 c.... Find value of polynomial. poly = a(0) perr = 0.0 ! Error in polynomial values. xperr = 0.0 xpow = 1.0 do 120 n = 1, np xperr = abs (xold) * (xperr + tolx * abs (xpow)) xpow = xpow * xold perr = perr + abs (a(n)) * (xperr + tolx * abs (xpow)) poly = poly + a(n) * xpow 120 continue cbugcbugc***DEBUG begins. cbugcbug 9911 format ('n=',i5,' x,p,e=',1p3e20.12) cbugcbug write ( 3, 9911) niter, xold, poly, perr cbugcbugc***DEBUG ends. c.... Use iterate if polynomial value is less than allowed error. if (abs (poly) .le. perr) then xroot = xold cbugcbugc***DEBUG begins. cbugcbug 9912 format ('CONVERGED.') cbugcbug write ( 3, 9912) cbugcbugc***DEBUG ends. go to 210 endif c.... Find value of slope of polynomial. slope = a(1) serr = 0.0 ! Error in slope value. xperr = 0.0 xpow = 1.0 do 130 n = 2, np xperr = abs (xold) * (xperr + tolx * abs (xpow)) xpow = xpow * xold serr = serr + n * abs (a(n)) * (xperr + tolx * abs (xpow)) slope = slope + n * a(n) * xpow 130 continue c.... Pick an x at random in interval if slope is zero. if (abs (slope) .le. serr) then cbugcbugc***DEBUG begins. cbugcbug 9913 format ('PICKING NEW XINT. SLOPE =',1pe22.14) cbugcbug write ( 3, 9913) slope cbugcbugc***DEBUG ends. xnew = xmin + ranf () * (xmax - xmin) go to 200 endif c.... Find next iterate. xnew = xold - poly / slope if (xnew .lt. xmin) xnew = xmin if (xnew .gt. xmax) xnew = xmax c.... Pick an x at random in interval if stuck at xmin or xmax. if ((xold .eq. xmin) .and. (xnew .eq. xmin)) then cbugcbugc***DEBUG begins. cbugcbug 9914 format ('STUCK AT XMIN. SLOPE =',1pe22.14) cbugcbug write ( 3, 9914) slope cbugcbugc***DEBUG ends. xnew = xmin + ranf () * (xmax - xmin) endif if ((xold .eq. xmax) .and. (xnew .eq. xmax)) then cbugcbugc***DEBUG begins. cbugcbug 9915 format ('STUCK AT XMAX. SLOPE =',1pe22.14) cbugcbug write ( 3, 9915) slope cbugcbugc***DEBUG ends. xnew = xmin + ranf () * (xmax - xmin) endif 200 xold = xnew if (niter .le. nitmax) go to 110 cbugcbugc***DEBUG begins. cbugcbug 9916 format ('TOO MANY ITERATIONS.') cbugcbug write ( 3, 9916) cbugcbugc***DEBUG ends. nerr = 5 210 continue cbugc***DEBUG begins. cbug 9940 format (/ 'aptnewt results. nerr=',i3,' nit=',i5, cbug & ' xroot=',1pe22.14) cbug write ( 3, 9940) nerr, niter, xroot cbugc***DEBUG ends. return c.... End of subroutine aptnewt. (+1 line.) end UCRL-WEB-209832