Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 92
Текст из файла (страница 92)
The vector of functions to be zeroed, called fvec[1..n] in the routinebelow, is returned by the user-supplied routine vecfunc(n,x,fvec). The output quantitycheck is false (0) on a normal return and true (1) if the routine has converged to a localminimum of the function fmin defined below. In this case try restarting from a different initialguess.{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 lubksb(float **a, int n, int *indx, float b[]);void ludcmp(float **a, int n, int *indx, float *d);int i,its,j,*indx;float d,den,f,fold,stpmax,sum,temp,test,**fjac,*g,*p,*xold;indx=ivector(1,n);fjac=matrix(1,n,1,n);g=vector(1,n);p=vector(1,n);xold=vector(1,n);fvec=vector(1,n);Define global variables.nn=n;nrfuncv=vecfunc;f=fmin(x);fvec is also computed by this call.test=0.0;Test for initial guess being a root.
Usefor (i=1;i<=n;i++)more stringent test than simply TOLF.if (fabs(fvec[i]) > test) test=fabs(fvec[i]);if (test < 0.01*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);for (its=1;its<=MAXITS;its++) {Start of iteration loop.fdjac(n,x,fvec,fjac,vecfunc);If analytic Jacobian is available, you can replace the routine fdjac below with yourown routine.for (i=1;i<=n;i++) {Compute ∇f for the line search.for (sum=0.0,j=1;j<=n;j++) sum += fjac[j][i]*fvec[j];g[i]=sum;}for (i=1;i<=n;i++) xold[i]=x[i];Store x,fold=f;and f .for (i=1;i<=n;i++) p[i] = -fvec[i]; Right-hand side for linear equations.ludcmp(fjac,n,indx,&d);Solve linear equations by LU decompolubksb(fjac,n,indx,p);sition.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 valfor (i=1;i<=n;i++)ues.if (fabs(fvec[i]) > test) test=fabs(fvec[i]);if (test < TOLF) {*check=0;FREERETURN}if (*check) {Check for gradient of f zero, i.e., spuritest=0.0;ous convergence.den=FMAX(f,0.5*n);for (i=1;i<=n;i++) {temp=fabs(g[i])*FMAX(fabs(x[i]),1.0)/den;388Chapter 9.Root Finding and Nonlinear Sets of Equationsif (temp > test) test=temp;}*check=(test < TOLMIN ? 1 : 0);FREERETURN}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 newt");}#include <math.h>#include "nrutil.h"#define EPS 1.0e-4Approximate square root of the machine precision.void fdjac(int n, float x[], float fvec[], float **df,void (*vecfunc)(int, float [], float []))Computes forward-difference approximation to Jacobian. On input, x[1..n] is the point atwhich the Jacobian is to be evaluated, fvec[1..n] is the vector of function values at thepoint, and vecfunc(n,x,f) is a user-supplied routine that returns the vector of functions atx.
On output, df[1..n][1..n] is the Jacobian array.{int i,j;float h,temp,*f;f=vector(1,n);for (j=1;j<=n;j++) {temp=x[j];h=EPS*fabs(temp);if (h == 0.0) h=EPS;x[j]=temp+h;Trick to reduce finite precision error.h=x[j]-temp;(*vecfunc)(n,x,f);x[j]=temp;for (i=1;i<=n;i++) df[i][j]=(f[i]-fvec[i])/h;Forward difference for}mula.free_vector(f,1,n);}#include "nrutil.h"extern int nn;extern float *fvec;extern void (*nrfuncv)(int n, float v[], float f[]);float fmin(float x[])Returns f = 12 F · F at x.
The global pointer *nrfuncv points to a routine that returns thevector of functions at x. It is set to point to a user-supplied routine in the calling program.Global variables also communicate the function values back to the calling program.{int i;float sum;(*nrfuncv)(nn,x,fvec);for (sum=0.0,i=1;i<=nn;i++) sum += SQR(fvec[i]);return 0.5*sum;}3899.7 Globally Convergent Methods for Nonlinear Systems of EquationsThe routine newt assumes that typical values of all components of x and of F are of orderunity, and it can fail if this assumption is badly violated.
You should rescale the variables bytheir typical values before invoking newt if this problem occurs.Multidimensional Secant Methods: Broyden’s MethodNewton’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 ofBi · δxi = −Fi(9.7.15)where δxi = xi+1 − xi (cf. equation 9.7.3). The quasi-Newton or secant condition is thatBi+1 satisfyBi+1 · δxi = δFi(9.7.16)where δFi = Fi+1 − Fi . This is the generalization of the one-dimensional secant approximation to the derivative, δF/δx. However, equation (9.7.16) does not determine Bi+1 uniquelyin more than one dimension.Many different auxiliary conditions to pin down Bi+1 have been explored, but thebest-performing algorithm in practice results from Broyden’s formula.
This formula is basedon the idea of getting Bi+1 by making the least change to Bi consistent with the secantequation (9.7.16). Broyden showed that the resulting formula isBi+1 = Bi +(δFi − Bi · δxi ) ⊗ δxiδxi · δxi(9.7.17)You can easily check that Bi+1 satisfies (9.7.16).Early implementations of Broyden’s method used the Sherman-Morrison formula,equation (2.7.2), to invert equation (9.7.17) analytically,−1B−1+i+1 = Bi· δFi ) ⊗ δxi · B−1(δxi − B−1iiδxi · B−1· δFii(9.7.18)Then instead of solving equation (9.7.3) by e.g., LU decomposition, one determinedδxi = −B−1· Fii(9.7.19)by matrix multiplication in O(N ) operations.
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 ,2∇( 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.390Chapter 9.Root Finding and Nonlinear Sets of EquationsSince 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.#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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
