386 Chapter 9.Root Finding and Nonlinear Sets of Equations temp=fabs(p[i])/FMAX(fabs(xold[i]),1.0); if (temp test)test=temp; alamin=TOLX/test; alam=1.0: Always try full Newton step first. for(;;){ Start of iteration loop. for (i=1;i<=n;i++)x[i]=xold[i]talam*p[i]; *f=(*func)(x); if (alam alamin){ Convergence on Ar.For zero find- for (i=1;i<=n;i++)x[i]=xold[i]; ing,the calling program should *check=1; verify the convergence return: Sufficient function decrease. 5常 else if (*f <fold+ALF*alam*slope)return; 8 else Backtrack. 1f(a1am==1.0) tmplam =-slope/(2.0*(*f-fold-slope)); First time. 1881992 else Subsequent backtracks. rhs1 *f-fold-alam*slope; 11.00 rhs2=f2-fold-alam2+slope; a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2); b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2); if (a ==0.0)tmplam =-slope/(2.0*b); else 6 disc=b*b-3.0*a*slope; (Nort serve if (disc 0.0)tmplam=0.5*alam; else if (b <0.0)tmplam=(-b+sgrt(disc))/(3.0*a); else tmplam=-slope/(b+sqrt(disc)); America computer one paper University Press. ART if(tmplam>0.5*alam） tmplam=0.5*alam; 入≤0.5λ1 Programs alam2=alam: for their f2=*f; alam=FMAX(tmplam,0.1*alam); λ20.1入1- Try again. to dir Here now is the globally convergent Newton routine newt that uses Insrch.A feature 1788-1982 SCIENTIFIC COMPUTING(ISBN of newt is that you need not supply the Jacobian matrix analytically;the routine will attempt to compute the necessary partial derivatives of F by finite differences in the routine fdjac.This routine uses some of the techniques described in $5.7 for computing numerical derivatives.Of course,you can always replace fdjac with a routine that calculates the Jacobian analytically 10-621 if this is easy for you to do. #include <math.h> Numerical Recipes -43108 #include "nrutil.h" #define MAXITS 200 #define TOLF 1.0e-4 (outside #define TOLMIN 1.0e-6 Software. #define TOLX 1.0e-7 North #define STPMX 100.0 Here MAXITS is the maximum number of iterations;TOLF sets the convergence criterion on function values:TOLMIN sets the criterion for deciding whether spurious convergence to a minimum of fmin has occurred;TOLX is the convergence criterion on ox;STPMX is the scaled maximum step length allowed in line searches. int nni Global variables to communicate with fmin. float *fvec; void (*nrfuncv)(int n,float v[],float f[]); #define FREERETURN {free_vector(fvec,1,n);free_vector(xold,1,n);\ free_vector(p,1,n);free_vector(g,1,n);free_matrix(fjac,1,n,1,n);\ free_ivector(indx,1,n);return;386 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). temp=fabs(p[i])/FMAX(fabs(xold[i]),1.0); if (temp > test) test=temp; } alamin=TOLX/test; alam=1.0; Always try full Newton step first. for (;;) { Start of iteration loop. for (i=1;i<=n;i++) x[i]=xold[i]+alam*p[i]; *f=(*func)(x); if (alam < alamin) { Convergence on ∆x. For zero find￾ing, the calling program should verify the convergence. for (i=1;i<=n;i++) x[i]=xold[i]; *check=1; return; } else if (*f <= fold+ALF*alam*slope) return; Sufficient function decrease. else { Backtrack. if (alam == 1.0) tmplam = -slope/(2.0*(*f-fold-slope)); First time. else { Subsequent backtracks. rhs1 = *f-fold-alam*slope; rhs2=f2-fold-alam2*slope; a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2); b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2); if (a == 0.0) tmplam = -slope/(2.0*b); else { disc=b*b-3.0*a*slope; if (disc < 0.0) tmplam=0.5*alam; else if (b <= 0.0) tmplam=(-b+sqrt(disc))/(3.0*a); else tmplam=-slope/(b+sqrt(disc)); } if (tmplam > 0.5*alam) tmplam=0.5*alam; λ ≤ 0.5λ1. } } alam2=alam; f2 = *f; alam=FMAX(tmplam,0.1*alam); λ ≥ 0.1λ1. } Try again. } Here now is the globally convergent Newton routine newt that uses lnsrch. A feature of newt is that you need not supply the Jacobian matrix analytically; the routine will attempt to compute the necessary partial derivatives of F by finite differences in the routine fdjac. This routine uses some of the techniques described in §5.7 for computing numerical derivatives. Of course, you can always replace fdjac with a routine that calculates the Jacobian analytically if this is easy for you to do. #include <math.h> #include "nrutil.h" #define MAXITS 200 #define TOLF 1.0e-4 #define TOLMIN 1.0e-6 #define TOLX 1.0e-7 #define STPMX 100.0 Here MAXITS is the maximum number of iterations; TOLF sets the convergence criterion on function values; TOLMIN sets the criterion for deciding whether spurious convergence to a minimum of fmin has occurred; TOLX is the convergence criterion on δx; STPMX is the scaled maximum step length allowed in line searches. int nn; Global variables to communicate with fmin. float *fvec; void (*nrfuncv)(int n, float v[], float f[]); #define FREERETURN {free_vector(fvec,1,n);free_vector(xold,1,n);\ free_vector(p,1,n);free_vector(g,1,n);free_matrix(fjac,1,n,1,n);\ free_ivector(indx,1,n);return;}