Chapter 17.Two Point Boundary Value Problems .com or call 装 17.0 Introduction 11-800-72 When ordinary differential equations are required to satisfy boundary conditions Cambridge at more than one value of the independent variable,the resulting problem is called a NUMERICAL RECIPES IN two point boundary value problem.As the terminology indicates,the most common server case by far is where boundary conditions are supposed to be satisfied at two points- usually the starting and ending values of the integration.However,the phrase"two compu point boundary value problem"is also used loosely to include more complicated Press. cases,e.g.,where some conditions are specified at endpoints,others at interior C:THE ART OF (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 SCIENTIFIC to start an acceptable solution at its beginning(initial values)and just march it along by numerical integration to its end(final values);while in the present case,the 6 boundary conditions at the starting point do not determine a unique solution to start with-and a"random"choice among the solutions that satisfy these (incomplete) COMPUTING starting boundary conditions is almost certain not to satisfy the boundary conditions at the other specified point(s). r Numerical 18-1992bm It should not surprise you that iteration is in general required to meld these Further spatially scattered boundary conditions into a single global solution of the differential Recipes equations.For this reason,two point boundary value problems require considerably more effort to solve than do initial value problems.You have to integrate your dif- ferential equations over the interval of interest,or perform an analogous"relaxation" procedure(see below),at least several,and sometimes very many,times.Only in (outside the special case of linear differential equations can you say in advance just how Software. many such iterations will be required. North The "standard"two point boundary value problem has the following form:We desire the solution to a set of N coupled first-order ordinary differential equations, America) satisfying n boundary conditions at the starting point 1,and a remaining set of visit website machine n2 =N-n boundary conditions at the final point z2.(Recall that all differential equations of order higher than first can be written as coupled sets of first-order equations,cf.816.0.) The differential equations are dyi(x) =9(c,1,2,·,N) i=1,2,.,N (17.0.1) d红 753
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Chapter 17. Two Point Boundary Value Problems 17.0 Introduction When ordinary differential equations are required to satisfy boundary conditions at more than one value of the independent variable, the resulting problem is called a two point boundary value problem. As the terminology indicates, the most common case by far is where boundary conditions are supposed to be satisfied at two points — 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 to start an acceptable solution at its beginning (initial values) and just march it along by numerical integration to its end (final values); while in the present case, the boundary conditions at the starting point do not determine a unique solution to start with — and a “random” choice among the solutions that satisfy these (incomplete) starting boundary conditions is almost certain not to satisfy the boundary conditions at the other specified point(s). It should not surprise you that iteration is in general required to meld these spatially scattered boundary conditions into a single global solution of the differential equations. For this reason, two point boundary value problems require considerably more effort to solve than do initial value problems. You have to integrate your differential equations over the interval of interest, or perform an analogous “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 boundary value problem has the following form: We desire the solution to a set of N coupled first-order ordinary differential equations, satisfying n1 boundary conditions at the starting point x1, and a remaining set of n2 = N − n1 boundary conditions 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 dyi(x) dx = gi(x, y1, y2,...,yN ) i = 1, 2,...,N (17.0.1) 753
754 Chapter 17.Two Point Boundary Value Problems ① ③ y desired boundary value ② required 83g boundary value r Figure 17.0.1.Shooting method(schematic).Trial integrations that satisfy the boundary condition at one endpoint are"launched."The discrepancies from the desired boundary condition at the other endpoint are used to adjust the starting conditions,until boundary conditions at both endpoints are ultimately satisfied. RECIPES I At z1,the solution is supposed to satisfy B1(1,1,2,·,N)=0 j=1,...,1 (17.0.2) N旧的 Press. while at r2,it is supposed to satisfy B2k(t2,h1,2,,N）=0k=1,,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 6 of the dependent variables at one boundary.These values must be consistent with any boundary conditions for that boundary,but otherwise are arranged to depend on arbitrary free parameters whose values we initially "randomly"guess.We then integrate the ODEs by initial value methods,arriving at the other boundary(and/or any interior points with boundary conditions specified).In general,we find discrepancies 10.621 from the desired boundary values there.Now we have a multidimensional root- Recipes Numerica finding problem,as was treated in 89.6 and 89.7:Find the adjustment of the free 43106 parameters at the starting point that zeros the discrepancies at the other boundary (outside Recipes point(s).If we liken integrating the differential equations to following the trajectory of a shot from gun to target,then picking the initial conditions corresponds to aiming (see Figure 17.0.1).The shooting method provides a systematic approach to taking North a set of"ranging shots that allow us to improve our "aim systematically. As another variant of the shooting method (817.2),we can guess unknown free parameters at both ends of the domain,integrate the equations to a common midpoint, and seek to adjust the guessed parameters so that the solution joins"smoothly"at the fitting 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
754 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). required boundary value desired boundary value 1 3 2 y x Figure 17.0.1. Shooting method (schematic). Trial integrations that satisfy the boundary condition at one endpoint are “launched.” The discrepancies from the desired boundary condition at the other endpoint are used to adjust the starting conditions, until boundary conditions at both endpoints are ultimately satisfied. At x1, the solution is supposed to satisfy B1j (x1, y1, y2,...,yN )=0 j = 1,...,n1 (17.0.2) while at x2, it is supposed to satisfy B2k(x2, y1, y2,...,yN )=0 k = 1,...,n2 (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 with any boundary conditions for that boundary, but otherwise are arranged to depend on arbitrary free parameters whose values we initially “randomly” guess. We then integrate the ODEs by initial value methods,arriving at the other boundary (and/or any interior points with boundary conditions specified). In general, we find discrepancies 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 liken integrating the differential equations to following the trajectory of a shot from gun to target, then picking the initial conditions corresponds to aiming (see Figure 17.0.1). The shooting method provides a systematic approach to taking a set of “ranging” shots that allow us to improve our “aim” systematically. As another variant of the shooting method (§17.2), we can guess unknown free parameters at both ends of the domain, integrate the equations to a common midpoint, and seek to adjust the guessed parameters so that the solution joins “smoothly” at the fitting 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 755 Ist iteration initial guess 2nd iteration required boundary true solution value required 83 boundary 鱼 18881892 value 三g Figure 17.0.2.Relaxation method(schematic).An initial solution is guessed that approximately satisfies the differential equation and boundary conditions.An iterative process adjusts the function to bring it into close agreement with the true solution. from NUMERICAL RECIPESI the integration.A trial solution consists of values for the dependent variables at each 令 mesh point,not satisfying the desired finite-difference equations,nor necessarily even satisfying the required boundary conditions.The iteration,now called relaxation, consists of adjusting all the values on the mesh so as to bring them into successively America Press. closer agreement with the finite-difference equations and,simultaneously,with the boundary conditions(see Figure 17.0.2).For example,if the problem involves three 9 coupled equations and a mesh of one hundred points,we must guess and improve three hundred variables representing the solution. With all this adjustment,you may be surprised that relaxation is ever an efficient 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 closed form.Relaxation works best when the solution is smooth and not highly oscillatory. Such oscillations would require many grid points for accurate representation.The number and position of required points may not be known apriori.Shooting methods are usually preferred in such cases,because their variable stepsize integrations adjust naturally to a solution's peculiarities. Numerical 10-521 Relaxation methods are often preferred when the ODEs have extraneous 43106 solutions which,while not appearing in the final solution satisfying all boundary (outside Recipes 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. North Good initial guesses are the secret of efficient relaxation methods.Often one has to solve a problem many times,each time with a slightly different value of some parameter.In that case,the previous solution is usually a good initial guess when the parameter is changed,and relaxation will work well. Until you have enough experience to make your own judgment between the two methods,you might wish to follow the advice of your authors,who are notorious computer gunslingers:We always shoot first,and only then relax
17.0 Introduction 755 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). required boundary value required boundary value initial guess 1st iteration 2nd iteration true solution Figure 17.0.2. Relaxation method (schematic). An initial solution is guessed that approximately satisfies the differential equation and boundary conditions. An iterative process adjusts the function to bring it into close agreement with the true solution. the integration. A trial solution consists of values for the dependent variables at each mesh point, notsatisfying the desired finite-difference equations, nor necessarily even satisfying the required boundary conditions. The iteration, now called relaxation, consists of adjusting all the values on the mesh so as to bring them into successively closer agreement with the finite-difference equations and, simultaneously, with the boundary conditions (see Figure 17.0.2). For example, if the problem involves three coupled equations and a mesh of one hundred points, we must guess and improve three hundred variables representing the solution. With all this adjustment, you may be surprised that relaxation is ever an efficient 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 closed form. Relaxation works best when the solution is smooth and not highly oscillatory. Such oscillations would require many grid points for accurate representation. The number and position of required points may not be known a priori. Shooting methods are usually preferred in such cases, because their variable stepsize integrations adjust 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 efficient relaxation methods. Often one has to solve a problem many times, each time with a slightly different value of some parameter. In that case, the previous solution is usually a good initial guess when the parameter is changed, and relaxation will work well. Until you have enough experience to make your own judgment between the two methods, you might wish to follow the advice of your authors, who are notorious computer gunslingers: We always shoot first, and only then relax
756 Chapter 17.Two Point Boundary Value Problems Problems Reducible to the Standard Boundary Problem There are two important problems that can be reduced to the standard boundary value problem described by equations(17.0.1)-(17.0.3).The first is the eigemalue problem for differential equations.Here the right-hand side of the system of differential equations depends on a parameterA, d@=gc,h,,N,月 (17.0.4) d and one has to satisfy N+1 boundary conditions instead of just N.The problem 81 is overdetermined and in general there is no solution for arbitrary values of A.For certain special values of A,the eigenvalues,equation(17.0.4)does have a solution. We reduce this problem to the standard case by introducing a new dependent 容 variable yN+1三入 (17.0.5) and another differential equation dyN+10 3 dr (17.0.6) An example of this trick is given in 817.4. 9 The other case that can be put in the standard form is a free boundary problem. Here only one boundary abscissa z is specified,while the other boundary 2 is to be determined so that the system (17.0.1)has a solution satisfying a total of N+1 boundary conditions.Here we again add an extra constant dependent variable: 6 yN+1三x2-x1 (17.0.7) dyN+i=0 (17.0.8) dx 3方豆 Numerica 10.521 We also define a new independent variable t by setting 43126 x-x1三tyN+1, 0≤t≤1 (17.0.9) The system of N+1 differential equations for dyi/dt is now in the standard form, with t varying between the known limits 0 and 1. North 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 Royal Astronomical Society,vol.151,pp.351-364 London,R.A.,and Flannery.B.P.1982.Astrophysical Journal,vol.258,pp.260-269. Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). 887.3-7.4
756 Chapter 17. Two Point Boundary Value Problems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Problems Reducible to the Standard Boundary Problem There are two important problems that can be reduced to the standard boundary value problem described by equations (17.0.1) – (17.0.3). The first is the eigenvalue problem for differential equations. Here the right-hand side of the system of differential equations depends on a parameter λ, dyi(x) dx = gi(x, y1,...,yN , λ) (17.0.4) and one has to satisfy N + 1 boundary conditions instead of just N. The problem is overdetermined and in general there is no solution for arbitrary values of λ. For certain special values of λ, the eigenvalues, equation (17.0.4) does have a solution. We reduce this problem to the standard case by introducing a new dependent variable yN+1 ≡ λ (17.0.5) and another differential equation dyN+1 dx =0 (17.0.6) An example of this trick is given in §17.4. The other case that can be put in the standard form is a free boundary problem. Here only one boundary abscissa x1 is specified, while the other boundary x2 is to be determined so that the system (17.0.1) has a solution satisfying a total of N + 1 boundary conditions. Here we again add an extra constant dependent variable: yN+1 ≡ x2 − x1 (17.0.7) dyN+1 dx =0 (17.0.8) We also define a new independent variable t by setting x − x1 ≡ t yN+1, 0 ≤ t ≤ 1 (17.0.9) The system of N + 1 differential equations for dyi/dt is now in the standard form, with t varying between the known limits 0 and 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 Royal Astronomical Society, vol. 151, pp. 351–364. London, R.A., and Flannery, B.P. 1982, Astrophysical Journal, vol. 258, pp. 260–269. Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §§7.3–7.4
17.1 The Shooting Method 757 17.1 The Shooting Method In this section we discuss"pure"shooting,where the integration proceeds from I to 2,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(89.7).It seeks to zero n2 functions 81 of n2 variables.The functions are obtained by integrating N differential equations from to x2.Let us see how this works: At the starting point 1 there are N starting values yi to be specified,but 、 subject to n conditions.Therefore there are n2=N-n1 freely specifiable starting values.Let us imagine that these freely specifable values are the components of a vector V that lives in a vector space of dimension n2.Then you,the user,knowing the functional form of the boundary conditions(17.0.2),can write a function that generates a complete set of N starting values y,satisfying the boundary conditions at 1,from an arbitrary vector value of V in which there are no restrictions on the n2 component values.In other words,(17.0.2)converts to a prescription (c1)=班(1；，，V2)i=1,,N (17.1.1) 经6 Below,the function that implements(17.1.1)will be called load. Notice that the components of V might be exactly the values of certain"free" b云P components of y,with the other components of y determined by the boundary 61 conditions.Alternatively,the components of V might parametrize the solutions that satisfy the starting boundary conditions in some other convenient way.Boundary conditions often impose algebraic relations among the yi,rather than specific values 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) Numerica 10621 all allowed starting vectors y Given a particular V,a particular y()is thus generated.It can then be turned 431 into a y(2)by integrating the ODEs to x2 as an initial value problem(e.g,using Recipes 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 2(17.0.3).Simplest of all is just to use the right-hand North sides of(17.0.3), f=B2k(c2,y）k=1,,n2 (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 F equal to zero if and only if the boundary conditions at x2 are satisfied.Below,you will be asked to supply a user-written function score which uses(17.0.3)to convert an N-vector of ending values y(2)into an n2-vector of discrepancies F
17.1 The Shooting Method 757 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 17.1 The Shooting Method In this section we discuss “pure” shooting, where the integration proceeds from x1 to x2, 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 multidimensional, globally convergent Newton-Raphson (§9.7). It seeks to zero n 2 functions of n2 variables. The functions are obtained by integrating N differential equations from x1 to x2. Let us see how this works: At the starting point x1 there are N starting values yi to be specified, but subject to n1 conditions. Therefore there are n2 = N − n1 freely specifiable starting values. Let us imagine that these freely specifiable values are the components of a vector V that lives in a vector space of dimension n2. Then you, the user, knowing the functional form of the boundary conditions (17.0.2), can write a function that generates a complete set of N starting values y, satisfying the boundary conditions at x1, from an arbitrary vector value of V in which there are no restrictions on the n 2 component values. In other words, (17.0.2) converts to a prescription yi(x1) = yi(x1; V1,...,Vn2 ) i = 1,...,N (17.1.1) Below, the function that implements (17.1.1) will be called load. Notice that the components of V might be exactly the values of certain “free” components of y, with the other components of y determined by the boundary conditions. Alternatively, the components of V might parametrize the solutions that satisfy the starting boundary conditions in some other convenient way. Boundary conditions often impose algebraic relations among the y i, rather than specific values 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. Given a particular V, a particular y(x1) is thus generated. It can then be turned into a y(x2) by integrating the ODEs to x2 as 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,...,n2 (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 F equal to zero if and only if the boundary conditions at x2 are satisfied. Below, you will be asked to supply a user-written function score which uses (17.0.3) to convert an N-vector of ending values y(x2) into an n2-vector of discrepancies F