392 Chapter 9. Root Finding and Nonlinear Sets of Equations for (den=0.0,i=1;i<=n;i++)den +SQR(s[i]); for (i=1;i<-n;i++)s[i]/=den; qrupdt(r,qt,n,t,s); 5pa2报0 for(i=1;1<=n;1++){ if (r[i][i]==0.0)nrerror("r singular in broydn"); d[i]=r[i][i]; Diagonal of R stored in d. 2 for(i=1;i<=n;i++)[ Right-hand side for linear equations is-QT.F. for (sum=0.0,j=1;j<=n;j++)sum qt[i][j]*fvec[j]; p[i]=-sum; http://www.nr Permission is read able files for (i=n;i>=1;i--){ Computef(Q.R)T.F for the line search. for(sum=0.0,j=1:j<=1;j++) sum -r[j][i]*p[j]; 83g g[i]=sum; granted for 19881992 for (i=1;i<-n;i++){ Store x and F. 1600 (including this one) xold[i]=x[i]; fvcold[i]=fvec[i]; 872 fold=f; Store f. 7423 to any server computer, by Cambridge University Press. users to make from NUMERICAL RECIPES IN rsolv(r,n,d,p); Solve linear equations Insrch(n,xold,fold,g,p,x,&f,stpmax,check,fmin); Insrch returns new x and f.It also calculates fvec at the new x when it calls fmin. (North testso0: Test for convergence on function values. for(i=1;1<=n;i++) America THE if (fabs(fvec[i])>test) test=fabs(fvec[i]); one paper 是 ART if (test TOLF){ *check=0; FREERETURN Programs if (*check){ True if line search failed to find a new x. if (restrt)FREERETURN Failure:already tried reinitializing the Jaco- strictly prohibited else bian. test=0.0; Check for gradient of f zero,i.e. spurious to dir den=FMAX(f,0.5*n); convergence. for(i=1;1<=n;i+){ temp=fabs(g[i])*FMAX(fabs(x [i]),1.0)/den; OF SCIENTIFIC COMPUTING(ISBN if (temp test)test=temp; if (test TOLMIN)FREERETURN else restrt=1; Try reinitializing the Jacobian 10-621 else Successful step:will use Broyden update for restrt=0; next step. Numerical Recipes books or .Further reproduction, 1988-1992 by Numerical Recipes test=0.0; Test for convergence on x. 431085 for(i=1;i<=n;1++)C temp=(fabs(x [i]-xold[i]))/FMAX(fabs(x[i]),1.0); if (temp test)test=temp; (outside North Amer Software. if (test TOLX)FREERETURN nrerror("MAXITS exceeded in broydn"); FREERETURN392 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). for (den=0.0,i=1;i<=n;i++) den += SQR(s[i]); for (i=1;i<=n;i++) s[i] /= den; Store s/(s · s) in s. qrupdt(r,qt,n,t,s); Update R and QT . for (i=1;i<=n;i++) { if (r[i][i] == 0.0) nrerror("r singular in broydn"); d[i]=r[i][i]; Diagonal of R stored in d. } } } for (i=1;i<=n;i++) { Right-hand side for linear equations is −QT · F. for (sum=0.0,j=1;j<=n;j++) sum += qt[i][j]*fvec[j]; p[i] = -sum; } for (i=n;i>=1;i--) { Compute ∇f ≈ (Q · R)T · F for the line search. for (sum=0.0,j=1;j<=i;j++) sum -= r[j][i]*p[j]; g[i]=sum; } for (i=1;i<=n;i++) { Store x and F. xold[i]=x[i]; fvcold[i]=fvec[i]; } fold=f; Store f. rsolv(r,n,d,p); Solve linear equations. lnsrch(n,xold,fold,g,p,x,&f,stpmax,check,fmin); lnsrch returns new x and f. It also calculates fvec at the new x when it calls fmin. test=0.0; Test for convergence on function values. for (i=1;i<=n;i++) if (fabs(fvec[i]) > test) test=fabs(fvec[i]); if (test < TOLF) { *check=0; FREERETURN } if (*check) { True if line search failed to find a new x. if (restrt) FREERETURN Failure; already tried reinitializing the Jaco￾else { bian. test=0.0; Check for gradient of f zero, i.e., spurious den=FMAX(f,0.5*n); convergence. for (i=1;i<=n;i++) { temp=fabs(g[i])*FMAX(fabs(x[i]),1.0)/den; if (temp > test) test=temp; } if (test < TOLMIN) FREERETURN else restrt=1; Try reinitializing the Jacobian. } } else { Successful step; will use Broyden update for restrt=0; next step. test=0.0; Test for convergence on δx. for (i=1;i<=n;i++) { temp=(fabs(x[i]-xold[i]))/FMAX(fabs(x[i]),1.0); if (temp > test) test=temp; } if (test < TOLX) FREERETURN } } nrerror("MAXITS exceeded in broydn"); FREERETURN }