f17-1

17.1TheShootingMethod 749Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).17.1 The Shooting Method Inthissectionwediscuss“pure”shooting,wheretheintegrationproceedsfrom x1tox2, and we try to match boundary conditions at the end of the integration. In the next section, we describe shooting to an intermediate fitting point, where the solution to the equations and boundary conditions is found by launching “shots” from both sides of the interval and trying to match continuity conditions at some intermediate point. Our implementation of the shooting method exactly implements multidimen- sional, globally convergent Newton-Raphson ( §9.7). It seeks to zero n2functions ofn2variables. The functions are obtained by integrating Ndifferential equations from x1tox2. Let us see how this works: At the starting point x1there are Nstarting values yito be specified, but subjectto n1conditions. Thereforethereare n2=N−n1freelyspecifiable starting values. Let us imagine that these freely specifiable values are the components of a vectorVthat lives in a vector space of dimension n2. Then you, the user, knowing the functionalform of the boundaryconditions (17.0.2),can write a subroutine thatgenerates a complete set of Nstarting values y, satisfying the boundary conditions atx 1,fromanarbitraryvectorvalueof Vinwhichtherearenorestrictionsonthe n2 component values. In other words, (17.0.2) converts to a prescription yi(x1)=yi(x1;V1,...,V n2) i=1,...,N (17.1.1 ) Below, the subroutine that implements (17.1.1) will be called load. Notice that the components of Vmight be exactly the values of certain “free” components of y, with the other components of ydetermined by the boundary conditions. Alternatively,the componentsof Vmightparametrizethe solutions that satisfy the starting boundary conditions in some other convenient way. Boundary conditionsoftenimposealgebraicrelationsamongthe yi, ratherthanspecificvalues for each of them. Using some auxiliary set of parameters often makes it easier to “solve” the boundary relations for a consistent set of yi’s. It makes no difference which way you go, as long as your vector space of V’s generates (through 17.1.1) all allowed starting vectors y. Givena particular V, a particular y(x1)is thusgenerated. Itcan thenbeturned into ay(x2)by integrating the ODEs to x2as an initial value problem (e.g., using Chapter 16’s odeint). Now, at x2, let us define a discrepancy vector F, also of dimension n2, whose components measure how far we are from satisfying the n2 boundary conditions at x2(17.0.3). Simplest of all is just to use the right-hand sides of (17.0.3), Fk=B2k(x2,y) k=1,...,n 2 (17.1.2 ) As in the case of V, however, you can use any other convenient parametrization, as long as your space of F’s spans the space of possible discrepancies from the desired boundary conditions, with all components of Fequal to zero if and only if the boundary conditions at x2are satisfied. Below, you will be asked to supply a user-writtensubroutine scorewhichuses(17.0.3)toconvertan N-vectorofending valuesy(x2)into an n2-vector of discrepancies F. 750 Chapter17. TwoPointBoundaryValueProblemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Now, as far as Newton-Raphson is concerned, we are nearly in business. We want to find a vector value of Vthat zeros the vector value of F. We do this by invoking the globally convergent Newton’s method implemented in the routine newtof§9.7. Recall that the heart of Newton’s method involves solving the set ofn2linear equations J·δV=−F (17.1.3 ) and then adding the correction back, Vnew=Vold+δV (17.1.4 ) In (17.1.3), the Jacobian matrix Jhas components given by Jij=∂F i ∂V j(17.1.5 ) It is not feasible to compute these partial derivatives analytically. Rather, each requires a separateintegrationof the NODEs, followed by the evaluationof ∂F i ∂V j≈Fi(V1,...,V j+∆Vj,...)−Fi(V1,...,V j,...) ∆Vj(17.1.6 ) This is doneautomaticallyfor youin the routine fdjacthat comes with newt. The only input to newtthat you have to provide is the routine funcvthat calculates F by integrating the ODEs. Here is the appropriate routine: C SUBROUTINE shoot(n2,v,f) is named "funcv" for use with "newt" SUBROUTINE funcv(n2,v,f) INTEGER n2,nvar,kmax,kount,KMAXX,NMAX REAL f(n2),v(n2),x1,x2,dxsav,xp,yp,EPSPARAMETER (NMAX=50,KMAXX=200,EPS=1.e-6) At most NMAXcoupled ODEs. COMMON /caller/ x1,x2,nvar COMMON /path/ kmax,kount,dxsav,xp(KMAXX),yp(NMAX,KMAXX) C USES derivs,load,odeint,rkqs,score Routine for use with newtto solve a two point boundary value problem for nvarcoupled ODEs by shooting from x1tox2. Initial values for the nvarODEs at x1are generated fromthe n2inputcoefficients v(1:n2),usingtheuser-suppliedroutine load. Theroutine integratestheODEsto x2usingtheRunge-Kuttamethodwithtolerance EPS,initialstepsize h1, and minimum stepsize hmin.A t x2it calls the user-supplied subroutine scoreto evaluatethe n2functions f(1:n2)thatoughttobezerotosatisfytheboundaryconditions atx2. Thefunctions farereturned onoutput. newtusesagloballyconvergent Newton’s methodtoadjustthevaluesof vuntilthefunctions farezero. Theuser-suppliedsubroutine derivs(x,y,dydx) supplies derivative information to the ODE integrator (see Chapter 16). Thecommonblock callerreceives itsvaluesfromthemainprogram sothat funcv canhavethesyntaxrequiredby newt. Thecommonblock pathisincludedforcompatibility with odeint. INTEGER nbad,nok REAL h1,hmin,y(NMAX)EXTERNAL derivs,rkqskmax=0 h1=(x2-x1)/100. hmin=0.call load(x1,v,y) call odeint(y,nvar,x1,x2,EPS,h1,hmin,nok,nbad,derivs,rkqs) call score(x2,y,f)returnEND 17.2Shootingto aFittingPoint 751Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).For some problems the initial stepsize ∆Vmight depend sensitively upon the initial conditions. It is straightforwardto alter loadto include a suggested stepsize h1as another returned argumentand feed it to fdjacvia a commonblock. A complete cycle of the shooting method thus requires n2+1integrations of theNcoupled ODEs: one integration to evaluate the current degree of mismatch, andn2for the partial derivatives. Each new cycle requires a new round of n2+1 integrations. This illustrates the enormousextraeffortinvolvedin solvingtwo point boundary value problems compared with initial value problems. Ifthedifferentialequationsare linear,thenonlyonecompletecycleisrequired, since (17.1.3)–(17.1.4)should take us right to the solution. A second round can beuseful, however, in mopping up some (never all) of the roundoff error. As givenhere, shootuses thequalitycontrolledRunge-Kuttamethodof §16.2 to integrate the ODEs, but any of the other methods of Chapter 16 could just aswell be used. You,the user,must supply shootwith: (i) a subroutine load(x1,v,y) which returnsthe n-vector y(1:n)(satisfyingthestartingboundaryconditions,ofcourse), given the freely specifiable variables of v(1:n2) at the initial point x1; (ii) a subroutine score(x2,y,f) which returns the discrepancy vector f(1:n2) of the ending boundary conditions, given the vector y(1:n)at the endpoint x2; (iii) a starting vector v(1:n2); (iv) a subroutine derivsfor the ODE integration; and other obvious parameters as described in the header comment above. In§17.4 we give a sample program illustrating how to use shoot. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America). Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell). 17.2 Shooting to a Fitting Point Theshootingmethoddescribedin §17.1tacitlyassumedthatthe“shots”would be able to traverse the entire domain of integration, even at the early stages ofconvergence to a correct solution. In some problems it can happen that, for very wrong starting conditions, an initial solution can’t even get from x 1tox2without encounteringsome incalculable, or catastrophic, result. For example, the argument of a square root might go negative, causing the numerical code to crash. Simple shooting would be stymied. A different, but related, case is where the endpoints are both singular points of the set of ODEs. One frequently needs to use special methods to integrate near the singular points, analytic asymptotic expansions,for example. In such cases it isfeasible to integrate in the direction awayfrom a singular point, using the special method to get through the first little bit and then reading off “initial” values for further numerical integration. However it is usually not feasible to integrate into a singular point, if only because one has not usually expended the same analytic