c9-7 (779537), страница 3

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 3)

The disadvantage of this method is thatit cannot easily be embedded in a globally convergent strategy, for which the gradient ofequation (9.7.4) requires B, not B−1 ,∇( 12 F · F) ' BT · F(9.7.20)Accordingly, we implement the update formula in the form (9.7.17).However, we can still preserve the O(N 2 ) solution of (9.7.3) by using QR decomposition(§2.10) instead of LU decomposition. The reason is that because of the special form of equation(9.7.17), the QR decomposition of Bi can be updated into the QR decomposition of Bi+1 inO(N 2 ) operations (§2.10). All we need is an initial approximation B0 to start the ball rolling.It is often acceptable to start simply with the identity matrix, and then allow O(N ) updates toproduce a reasonable approximation to the Jacobian.

We prefer to spend the first N functionevaluations on a finite-difference approximation to initialize B via a call to fdjac.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Newton’s method as implemented above is quite powerful, but it still has severaldisadvantages. One drawback is that the Jacobian matrix is needed. In many problemsanalytic derivatives are unavailable. If function evaluation is expensive, then the cost offinite-difference determination of the Jacobian can be prohibitive.Just as the quasi-Newton methods to be discussed in §10.7 provide cheap approximationsfor the Hessian matrix in minimization algorithms, there are quasi-Newton methods thatprovide cheap approximations to the Jacobian for zero finding. These methods are often calledsecant methods, since they reduce to the secant method (§9.2) in one dimension (see, e.g., [1]).The best of these methods still seems to be the first one introduced, Broyden’s method [2].Let us denote the approximate Jacobian by B.

Then the ith quasi-Newton step δx iis the solution of390Chapter 9.Root Finding and Nonlinear Sets of Equations#include <math.h>#include "nrutil.h"#define MAXITS 200#define EPS 1.0e-7#define TOLF 1.0e-4#define TOLX EPS#define STPMX 100.0#define TOLMIN 1.0e-6Here MAXITS is the maximum number of iterations; EPS is a number close to the machineprecision; TOLF is the convergence criterion on function values; TOLX is the convergence criterionon δx; STPMX is the scaled maximum step length allowed in line searches; TOLMIN is used todecide whether spurious convergence to a minimum of fmin has occurred.#define FREERETURN {free_vector(fvec,1,n);free_vector(xold,1,n);\free_vector(w,1,n);free_vector(t,1,n);free_vector(s,1,n);\free_matrix(r,1,n,1,n);free_matrix(qt,1,n,1,n);free_vector(p,1,n);\free_vector(g,1,n);free_vector(fvcold,1,n);free_vector(d,1,n);\free_vector(c,1,n);return;}int nn;Global variables to communicate with fmin.float *fvec;void (*nrfuncv)(int n, float v[], float f[]);void broydn(float x[], int n, int *check,void (*vecfunc)(int, float [], float []))Given an initial guess x[1..n] for a root in n dimensions, find the root by Broyden’s methodembedded in a globally convergent strategy.

The vector of functions to be zeroed, calledfvec[1..n] in the routine below, is returned by the user-supplied routine vecfunc(n,x,fvec).The routine fdjac and the function fmin from newt are used. The output quantity checkis false (0) on a normal return and true (1) if the routine has converged to a local minimumof the function fmin or if Broyden’s method can make no further progress. In this case tryrestarting from a different initial guess.{void fdjac(int n, float x[], float fvec[], float **df,void (*vecfunc)(int, float [], float []));float fmin(float x[]);void lnsrch(int n, float xold[], float fold, float g[], float p[], float x[],float *f, float stpmax, int *check, float (*func)(float []));void qrdcmp(float **a, int n, float *c, float *d, int *sing);void qrupdt(float **r, float **qt, int n, float u[], float v[]);void rsolv(float **a, int n, float d[], float b[]);int i,its,j,k,restrt,sing,skip;float den,f,fold,stpmax,sum,temp,test,*c,*d,*fvcold;float *g,*p,**qt,**r,*s,*t,*w,*xold;c=vector(1,n);d=vector(1,n);fvcold=vector(1,n);g=vector(1,n);p=vector(1,n);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Since B is not the exact Jacobian, we are not guaranteed that δx is a descent direction forf = 12 F · F (cf. equation 9.7.5). Thus the line search algorithm can fail to return a suitable stepif B wanders far from the true Jacobian.

In this case, we reinitialize B by another call to fdjac.Like the secant method in one dimension, Broyden’s method converges superlinearlyonce you get close enough to the root. Embedded in a global strategy, it is almost as robustas Newton’s method, and often needs far fewer function evaluations to determine a zero.Note that the final value of B is not always close to the true Jacobian at the root, evenwhen the method converges.The routine broydn given below is very similar to newt in organization.

The principaldifferences are the use of QR decomposition instead of LU , and the updating formula insteadof directly determining the Jacobian. The remarks at the end of newt about scaling thevariables apply equally to broydn.9.7 Globally Convergent Methods for Nonlinear Systems of Equations391Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).qt=matrix(1,n,1,n);r=matrix(1,n,1,n);s=vector(1,n);t=vector(1,n);w=vector(1,n);xold=vector(1,n);fvec=vector(1,n);Define global variables.nn=n;nrfuncv=vecfunc;f=fmin(x);The vector fvec is also computed by thistest=0.0;call.for (i=1;i<=n;i++)Test for initial guess being a root.

Use moreif (fabs(fvec[i]) > test)test=fabs(fvec[i]);stringent test than simif (test < 0.01*TOLF) {ply TOLF.*check=0;FREERETURN}for (sum=0.0,i=1;i<=n;i++) sum += SQR(x[i]);Calculate stpmax for line searches.stpmax=STPMX*FMAX(sqrt(sum),(float)n);restrt=1;Ensure initial Jacobian gets computed.for (its=1;its<=MAXITS;its++) {Start of iteration loop.if (restrt) {fdjac(n,x,fvec,r,vecfunc);Initialize or reinitialize Jacobian in r.qrdcmp(r,n,c,d,&sing);QR decomposition of Jacobian.if (sing) nrerror("singular Jacobian in broydn");for (i=1;i<=n;i++) {Form QT explicitly.for (j=1;j<=n;j++) qt[i][j]=0.0;qt[i][i]=1.0;}for (k=1;k<n;k++) {if (c[k]) {for (j=1;j<=n;j++) {sum=0.0;for (i=k;i<=n;i++)sum += r[i][k]*qt[i][j];sum /= c[k];for (i=k;i<=n;i++)qt[i][j] -= sum*r[i][k];}}}for (i=1;i<=n;i++) {Form R explicitly.r[i][i]=d[i];for (j=1;j<i;j++) r[i][j]=0.0;}} else {Carry out Broyden update.for (i=1;i<=n;i++) s[i]=x[i]-xold[i];s = δx.for (i=1;i<=n;i++) {t = R · s.for (sum=0.0,j=i;j<=n;j++) sum += r[i][j]*s[j];t[i]=sum;}skip=1;for (i=1;i<=n;i++) {w = δF − B · s.for (sum=0.0,j=1;j<=n;j++) sum += qt[j][i]*t[j];w[i]=fvec[i]-fvcold[i]-sum;if (fabs(w[i]) >= EPS*(fabs(fvec[i])+fabs(fvcold[i]))) skip=0;Don’t update with noisy components of w.else w[i]=0.0;}if (!skip) {for (i=1;i<=n;i++) {t = QT · w.for (sum=0.0,j=1;j<=n;j++) sum += qt[i][j]*w[j];t[i]=sum;}392Chapter 9.Root Finding and Nonlinear Sets of Equationsfor (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.}}nrerror("MAXITS exceeded in broydn");FREERETURN}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Характеристики

Список файлов книги