f17-2

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 the Ncoupled ODEs: one integration to evaluate the current degree of mismatch, and n2for 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 752 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).effort to obtain expansions of “wrong” solutions near the singular point (those not satisfying the desired boundary condition). The solution to the above mentioned difficulties is shooting to a fitting point . Insteadofintegratingfrom x1tox2,weintegratefirstfrom x1tosomepoint xfthat isbetween x1and x2; and second from x2(in the opposite direction) to xf. If (as before) the number of boundary conditions imposed at x1isn1, and the number imposed at x2isn2, then there are n2freely specifiable starting values at x1and n1freely specifiable starting values at x2. (If you are confused by this, go back to §17.1.) We can therefore define an n2-vectorV(1)of starting parameters atx1, and a prescription load1(x1,v1,y) for mapping V(1)into aythat satisfies the boundary conditions at x1, yi(x1)= yi(x1;V(1)1 ,...,V (1) n2) i=1 ,...,N (17.2.1 ) Likewise we can define an n1-vectorV(2)of starting parameters at x2, and a prescription load2(x2,v2,y) formapping V(2)intoaythat satisfies the boundary conditions at x2, yi(x2)= yi(x2;V(2)1 ,...,V (2) n1) i=1 ,...,N (17.2.2 ) We thus have a total of Nfreely adjustable parameters in the combination of V(1)andV(2). The Nconditions that must be satisfied are that there be agreement inNcomponents of yatxfbetween the values obtained integrating from one side and from the other, yi(xf;V(1))= yi(xf;V(2)) i=1 ,...,N (17.2.3 ) In some problems, the Nmatching conditions can be better described (physically, mathematically,ornumerically)byusing Ndifferentfunctions Fi,i=1 ...N,each possiblydependingonthe Ncomponents yi. Inthosecases, (17.2.3)is replacedby Fi[y(xf;V(1))] = Fi[y(xf;V(2))] i=1 ,...,N (17.2.4 ) Intheprogrambelow,theuser-suppliedsubroutine score(xf,y,f) issupposed to map an input N-vectoryinto an output N-vectorF. In most cases, you can dummy this subroutine as the identity mapping. Shooting to a fitting point uses globally convergent Newton-Raphson exactly as in §17.1. Comparingclosely with the routine shootof the previoussection, you should have no difficultyin understandingthe followingroutine shootf. The main differences in use are that you have to supply both load1andload2. Also, in the callingprogramyoumustsupplyinitialguessesfor v1(1:n2) andv2(1:n1) . Once againa sample programillustrating shootingto a fitting point is given in §17.4. C SUBROUTINE shootf(n,v,f) is named "funcv" for use with "newt" SUBROUTINE funcv(n,v,f) INTEGER n,nvar,nn2,kmax,kount,KMAXX,NMAX REAL f(n),v(n),x1,x2,xf,dxsav,xp,yp,EPSPARAMETER (NMAX=50,KMAXX=200,EPS=1.e-6) Atmost NMAXequations. COMMON /caller/ x1,x2,xf,nvar,nn2 COMMON /path/ kmax,kount,dxsav,xp(KMAXX),yp(NMAX,KMAXX) C USES derivs,load1,load2,odeint,rkqs,score Routine for use with newtto solve a two point boundary value problem for nvarcou- pledODEsbyshooting from x1andx2toafittingpoint xf. Initialvaluesforthe nvar 17.3RelaxationMethods 753Sample 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).ODEsat x1 (x2)aregeneratedfromthe n2 (n1)coefficients v1 (v2),usingtheuser- supplied routine load1 (load2) . The coefficients v1andv2shouldbestoredinasin- gle array v(1:n1+n2)inthe mainprogram byan EQUIVALENCE statement of the form (v1(1),v(1)),(v2(1),v(n2+1)) . Theinputparameter n=n1+n2 =nvar. Therou- tineintegratestheODEsto xfusingtheRunge-Kuttamethodwithtolerance EPS,initial stepsize h1,andminimumstepsize hmin.A txfitcallstheuser-suppliedsubroutine score toevaluatethe nvarfunctions f1andf2thatoughttomatchat xf. Thedifferences fare returned onoutput. newtusesagloballyconvergent Newton’smethodtoadjusttheval- uesof vuntilthefunctions farezero. Theuser-suppliedsubroutine derivs(x,y,dydx) suppliesderivativeinformationtotheODEintegrator(seeChapter16). Thecommonblock callerreceives its values from the main program so that funcvcan have the syntax required by newt.S e t nn2 =n2in the main program. The common block pathis for compatibility with odeint. INTEGER i,nbad,nok REAL h1,hmin,f1(NMAX),f2(NMAX),y(NMAX)EXTERNAL derivs,rkqs kmax=0 h1=(x2-x1)/100.hmin=0.call load1(x1,v,y) Pathfrom x1toxfwithbesttrialvalues v1. call odeint(y,nvar,x1,xf,EPS,h1,hmin,nok,nbad,derivs,rkqs) call score(xf,y,f1)call load2(x2,v(nn2+1),y) Pathfrom x2toxfwithbesttrialvalues v2. call odeint(y,nvar,x2,xf,EPS,h1,hmin,nok,nbad,derivs,rkqs) call score(xf,y,f2)do 11i=1,n f(i)=f1(i)-f2(i) enddo 11 return END There are boundaryvalue problems where even shooting to a fitting point fails — the integration interval has to be partitioned by several fitting points with the solution being matched at each such point. For more details see [1]. 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). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §§7.3.5–7.3.6. [1] 17.3 Relaxation Methods Inrelaxation methods we replace ODEs by approximate finite-difference equations (FDEs) on a grid or mesh of points that spans the domain of interest. As a typical example,we could replace a general first-order differential equation dy dx=g(x, y )( 17.3.1 ) with an algebraic equation relating function values at two points k, k−1: yk−yk−1−(xk−xk−1)g/bracketleftbig1 2(xk+xk−1),1 2(yk+yk−1)/bracketrightbig =0 ( 17.3.2 )