f17-0

Sample 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).Chapter 17. Two Point Boundary Value Problems 17.0 Introduction Whenordinarydifferentialequationsarerequiredtosatisfyboundaryconditions at morethanonevalueofthe independentvariable,the resultingproblemis called a two pointboundaryvalueproblem . As theterminologyindicates,the mostcommon casebyfariswhereboundaryconditionsaresupposedtobesatisfiedattwopoints— usually the starting and ending values of the integration. However, the phrase “two point boundary value problem” is also used loosely to include more complicated cases, e.g., where some conditions are specified at endpoints, others at interior (usually singular) points. The crucial distinction between initial value problems (Chapter 16) and two point boundary value problems (this chapter) is that in the former case we are able tostartanacceptablesolutionatitsbeginning(initialvalues)andjustmarchitalongby numerical integration to its end (final values); while in the present case, the boundaryconditionsat the starting point do not determinea uniquesolution to start with — and a “random” choice among the solutions that satisfy these (incomplete) startingboundaryconditionsis almostcertain nottosatisfy theboundaryconditions at the other specified point(s). It should not surprise you that iteration is in general required to meld these spatiallyscatteredboundaryconditionsintoasingleglobalsolutionofthedifferential equations. Forthis reason,two pointboundaryvalueproblemsrequireconsiderablymore effort to solve than do initial value problems. You have to integrate your dif- ferentialequationsovertheintervalofinterest,orperformananalogous“relaxation” procedure (see below), at least several, and sometimes very many, times. Only in the special case of linear differential equations can you say in advance just how many such iterations will be required. The “standard”two point boundaryvalue problemhas the followingform: We desire the solution to a set of Ncoupled first-order ordinary differential equations, satisfying n 1boundary conditions at the starting point x1, and a remaining set of n2=N−n1boundaryconditions at the final point x2. (Recall that all differential equations of order higher than first can be written as coupled sets of first-order equations, cf. §16.0.) The differential equations are dy i(x) dx=gi(x, y 1,y2,...,y N) i=1,2,...,N (17.0.1 ) 745 746 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).required boundaryvaluedesiredboundaryvalue1 3 2y x Figure17.0.1. Shootingmethod(schematic). Trialintegrations thatsatisfytheboundarycondition atone endpoint are “launched. ”Thediscrepancies fromthedesired boundary condition attheother endpoint are usedtoadjustthestartingconditions, untilboundaryconditions atbothendpoints areultimately satis fied. Atx1, the solution is supposed to satisfy B1j(x1,y1,y2,...,y N)=0 j=1,...,n 1 (17.0.2 ) while at x2, it is supposed to satisfy B2k(x2,y1,y2,...,y N)=0 k=1,...,n 2 (17.0.3 ) There are two distinct classes of numerical methods for solving two point boundary value problems. In the shooting method (§17.1) we choose values for all of the dependent variables at one boundary. These values must be consistent withany boundary conditions for thatboundary, but otherwise are arranged to depend on arbitrary free parameters whose values we initially “randomly ”guess. We then integratetheODEsbyinitialvaluemethods,arrivingattheotherboundary(and/orany interiorpointswithboundaryconditionsspeci fied). Ingeneral,we finddiscrepancies from the desired boundary values there. Now we have a multidimensional root-finding problem, as was treated in §9.6 and §9.7: Find the adjustment of the free parameters at the starting point that zeros the discrepancies at the other boundary point(s). If we likenintegratingthe differentialequationsto followingthetrajectoryofashotfromguntotarget,thenpickingtheinitialconditionscorrespondstoaiming (see Figure 17.0.1). The shooting method provides a systematic approachto taking a set of“ranging”shots that allow us to improve our “aim”systematically. As anothervariantof the shootingmethod( §17.2),we canguess unknownfree parametersatbothendsofthedomain,integratetheequationstoacommonmidpoint,and seek to adjust the guessed parameters so that the solution joins “smoothly”at thefitting point. In all shooting methods, trial solutions satisfy the differential equations “exactly”(or as exactly as we care to make our numerical integration), but the trial solutions come to satisfy the required boundary conditions only after the iterations are finished. Relaxation methods use a different approach. The differential equations are replaced by finite-difference equations on a mesh of points that covers the range of 17.0 Introduction 747Sample 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).required boundaryvaluerequiredboundaryvalueinitialguess1stiteration 2nditeration true solution Figure 17.0.2. Relaxation method(schematic). Aninitial solution isguessed thatapproximately satis fies the differential equation and boundary conditions. An iterative process adjusts the function to bring it into close agreement with the true solution. theintegration. Atrialsolutionconsistsofvaluesforthedependentvariablesateach meshpoint, notsatisfyingthedesired finite-differenceequations,nornecessarilyeven satisfying the required boundary conditions. The iteration, now called relaxation , consists ofadjustingall thevaluesonthemeshso asto bringthemintosuccessivelycloser agreement with the finite-difference equations and, simultaneously, with the boundaryconditions(see Figure17.0.2). Forexample,if theprobleminvolvesthree coupled equations and a mesh of one hundred points, we must guess and improvethree hundred variables representing the solution. Withallthisadjustment,youmaybesurprisedthatrelaxationiseveranef ficient method, but (for the right problems) it really is! Relaxation works better than shooting when the boundary conditions are especially delicate or subtle, or where they involve complicated algebraic relations that cannot easily be solved in closedform. Relaxationworksbestwhenthesolutionis smoothandnothighlyoscillatory. Such oscillations would require many grid points for accurate representation. The numberandpositionofrequiredpointsmaynotbeknown apriori. Shootingmethods areusuallypreferredinsuchcases,becausetheirvariablestepsizeintegrationsadjust naturally to a solution ’s peculiarities. Relaxation methods are often preferred when the ODEs have extraneous solutions which, while not appearing in the final solution satisfying all boundary conditions, may wreak havoc on the initial value integrations required by shooting.The typical case is that of trying to maintain a dying exponential in the presence of growing exponentials. Good initial guesses are the secret of ef ficient relaxation methods. Often one hastosolveaproblemmanytimes,eachtimewitha slightlydifferentvalueofsome parameter. In that case, the previous solution is usually a good initial guess when the parameter is changed, and relaxation will work well. Untilyouhaveenoughexperiencetomakeyourownjudgmentbetweenthetwo methods, you might wish to follow the advice of your authors, who are notoriouscomputer gunslingers: We always shoot first, and only then relax. 748 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).Problems Reducible tothe Standard Boundary Problem Therearetwoimportantproblemsthatcanbereducedtothestandardboundary valueproblemdescribedbyequations(17.0.1) –(17.0.3). The first is theeigenvalue problem for differential equations . Here the right-hand side of the system of differential equations depends on a parameter λ, dy i(x) dx=gi(x, y 1,...,y N,λ)( 17.0.4 ) and one has to satisfy N+1boundary conditions instead of just N. The problem is overdeterminedand in general there is no solution for arbitrary values of λ.F o r certainspecial valuesof λ, theeigenvalues,equation(17.0.4)doeshavea solution. We reduce this problem to the standard case by introducing a new dependent variable yN+1≡λ (17.0.5 ) and another differential equation dy N+1 dx=0 ( 17.0.6 ) An example of this trick is given in §17.4. Theothercase that canbe putin the standardformis a free boundaryproblem . Here only one boundary abscissa x1is specified, while the other boundary x2is to be determined so that the system (17.0.1)has a solution satisfying a total of N+1 boundaryconditions. Here we againadd an extra constant dependentvariable: yN+1≡x2−x1 (17.0.7 ) dy N+1 dx=0 ( 17.0.8 ) We also de fin ean e w independent variable tby setting x−x1≡ty N+1, 0≤t≤1( 17.0.9 ) The system of N+1differential equations for dy i/dtis now in the standard form, with tvarying between the known limits 0and 1. CITED REFERENCES AND FURTHER READING: Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell). Kippenhan, R., Weigert, A., and Hofmeister, E. 1968, in Methods in Computational Physics , vol. 7 (New York: Academic Press), pp. 129ff. Eggleton, P.P. 1971, Monthly Notices of the RoyalAstronomical Society , vol. 151, pp.351–364. London, R.A., and Flannery, B.P. 1982, Astrophysical Journal , vol. 258, pp. 260–269. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §§7.3–7.4. 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 speci fied, but subjectto n1conditions. Thereforethereare n2=N−n1freelyspecifiable starting values. Let us imagine that these freely speci fiable 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, ratherthanspeci ficvalues 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 ’sodeint). Now, at x2, let us de fine adiscrepancy 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.