  GNU-Darwin Web
i/roots Home Manual Packages Global Index Keywords Quick Reference ``` /* ROOTS.I Collection of root, maximum, and minimum finders. \$Id\$ */ /* Copyright (c) 1996. The Regents of the University of California. All rights reserved. */ local roots ; /* DOCUMENT roots.i defines: nraphson - Newton-Raphson/bisection root solver (scalar) f_inverse - function inverse by Newton-Raphson (vectorized) mnbrent - Brent's method minimizer (scalar) mxbrent - Brent's method maximizer (scalar) */ /* ------------------------------------------------------------------------ */ func nraphson (f_and_dfdx, x0, x1, xerr) /* DOCUMENT nraphson(f_and_dfdx, x0, x1) or nraphson(f_and_dfdx, x0, x1, xerr) Find a root of a function by Newton-Raphson iteration, backed up by bisection if the convergence seems poor. The subroutine F_AND_DFDX must be defined as: func F_AND_DFDX (x, &f, &dfdx) returning both the function value f(x) and derivative dfdx(x). If F_AND_DFDX always returns dfdx==0, nraphson uses bisection. The value of x is constrained to lie within the interval from X0 to X1; the function values at these two points must have opposite sign. The iteration stops when the root is known to within XERR, or to machine precision if XERR is nil or zero. f_inverse is a "vectorized" version of nraphson. Based on rtsafe from Press, et. al. Numerical Recipes, Ch 9. SEE ALSO: mnbrent, mxbrent, f_inverse */ { if (is_void(xerr) || xerr<0.0) xerr= 0.0; x0= double(x0); x1= double(x1); local f, dfdx; /* get function value at endpoints -- derivatives unused */ f_and_dfdx, x0, f, dfdx; if (f == 0.0) return x0; dxo= f; f_and_dfdx, x1, f, dfdx; if (f == 0.0) return x1; if (f*dxo > 0.0) error, "f(x0) and f(x1) have same sign"; if (f > 0.0) { xlo= x0; xhi= x1; } else { xlo= x1; xhi= x0; } /* first guess is midpoint */ x= 0.5*(xhi+xlo); dx= dxo= abs(xhi-xlo); f_and_dfdx, x, f, dfdx; if (f < 0.0) xlo= x; else if (f > 0.0) xhi= x; else return x; for (i=1 ; i<=nr_maxits ; ++i) { xo= x; if (((x-x1)*dfdx-f)*((x-x0)*dfdx-f) >= 0.0 || abs(2.*f) > abs(dxo*dfdx)) { /* take bisection step if N-R step would be out of bounds or if previous step did not converge fast enough */ dxo= dx; dx= 0.5*(xhi-xlo); x= xlo+dx; } else { /* take N-R step */ dxo= dx; dx= f/double(dfdx); x-= dx; } /* quit on either machine precision or requested precision */ if (x==xo || abs(dx) 0.0)) error, "f(x0)-y and f(x1)-y have same sign"; dfdx= x0; mask= double(f > dxo); maskc= 1.0-mask; xlo= x0*mask + x1*maskc; xhi= x1*mask + x0*maskc; /* first guess is midpoint */ x= 0.5*(xhi+xlo); dx= dxo= abs(xhi-xlo); f_and_dfdx, x, f, dfdx; f-= y; mask= (f < 0.0); list= where(mask); if (numberof(list)) xlo(list)= x(list); list= where(!mask); if (numberof(list)) xhi(list)= x(list); for (i=1 ; i<=nr_maxits ; ++i) { xo= x; mask= ((((x-x1)*dfdx-f)*((x-x0)*dfdx-f) >= 0.0) | (abs(2.*f) > abs(dxo*dfdx))); list= where(mask); if (numberof(list)) { /* take bisection step where N-R step would be out of bounds or if previous step did not converge fast enough */ dxob= dx(list); xob= xlo(list); dxb= 0.5*(xhi(list)-xob); xb= xob+dxb; } else { xb= dxb= dxob= []; } list= where(!mask); if (numberof(list)) { /* otherwise take N-R step */ dxon= dx(list); dxn= f(list)/dfdx(list); xon= x(list); xn= xon-dxn; } else { xn= dxn= dxon= []; } x= merge(xb, xn, mask); dx= merge(dxb, dxn, mask); dxo= merge(dxob, dxon, mask); /* check for uniform convergence either to requested precision or to machine precision */ if (allof((x==xo) | (dx abserr) { /* attempt a trial parabolic fit to (w,v,x) -- the three smallest points seen so far */ r= (x-w)*(fx-fv); q= (x-v)*(fx-fw); p= (x-v)*q - (x-w)*r; q= 2.*(q-r); if (q > 0.0) p= -p; q= abs(q); e1= e; e= d; if (abs(p)>=abs(0.5*q*e1) || p<=q*(xlo-x) || p>=q*(xhi-x)) { /* take golden section step -- parabolic step is crazy */ e= (x>=xm? xlo-x : xhi-x); d= golden*e; } else { /* take inverse parabolic step */ d= p/q; u= x+d; if ((u-xlo)=xm? xlo-x : xhi-x); d= golden*e; } /* evaluate at next point, avoiding meaninglessly small step size */ u= x + (abs(d)>=abserr? d : sign(d)*abserr); fu= (from_mxbrent? -f(u) : f(u)); if (fu <= fx) { /* new point is best so far */ if (u >= x) xlo= x; else xhi= x; v= w; fv= fw; w= x; fw= fx; x= u; fx= fu; } else { if (u < x) xlo= u; else xhi= u; if (fu<=fw || w==x) { /* new point is second best so far */ v= w; fv= fw; w= u; fw= fu; } else if (fu<=fw || v==x || v==w) { /* new point is third best so far */ v= u; fv= fu; } } } error, "br_maxits iteration count exceeded"; } br_maxits= 100; /* ------------------------------------------------------------------------ */ ```