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
758 Chapter 17.Two Point Boundary Value Problems Now,as far as Newton-Raphson is concerned,we are nearly in business.We want to find a vector value of V that zeros the vector value of F.We do this by invoking the globally convergent Newton's method implemented in the routine newt of 89.7.Recall that the heart of Newton's method involves solving the set of n2 linear equations J.6V=-F (17.1.3) and then adding the correction back, ynew vold+v (17.1.4) In(17.1.3),the Jacobian matrix J has components given by F J=0 (17.1.5) It is not feasible to compute these partial derivatives analytically.Rather,each requires a separate integration of the N ODEs,followed by the evaluation of ≈E(y+△）-FB(,,y山 av; (17.1.6) △V This is done automatically for you in the routine fdjac that comes with newt.The 9 only input to newt that you have to provide is the routine vecfunc that calculates F by integrating the ODEs.Here is the appropriate routine,called shoot,that is to be passed as the actual argument in newt: #include "nrutil.h" #define EPS 1.0e-6 extern int nvar; Variables that you must define and set in your main pro- extern float x1,x2; gram. int kmax,kount; Communicates with odeint 色 float *xp,**yp,dxsav; Numerical uctio 431086 void shoot(intn,f1oatv▣，f1oatf▣) Routine for use with newt to solve a two point boundary value problem for nvar coupled ODEs Recipes by shooting from x1 to x2.Initial values for the nvar ODEs at x1 are generated from the n2 input coefficients v[1..n2],using the user-supplied routine load.The routine integrates the ODEs to x2 using the Runge-Kutta method with tolerance EPS,initial stepsize h1,and minimum stepsize hmin.At x2 it calls the user-supplied routine score to evaluate the n2 functions North Software. f[1..n2]that ought to be zero to satisfy the boundary conditions at x2.The functions f are returned on output.newt uses a globally convergent Newton's method to adjust the values of v until the functions f are zero.The user-supplied routine derivs(x,y,dydx)supplies derivative information to the ODE integrator (see Chapter 16).The first set of global variables above receives its values from the main program so that shoot can have the syntax required for it to be the argument vecfunc of newt. void derivs(float x,float y[],float dydx[]); void1oad(f1oatx1,f1oatv[0,f1oaty[☐)； void odeint(float ystart],int nvar,float xi,float x2 float eps,float hi,float hmin,int *nok,int *nbad, void(*derivs.)(f1oat,f1oat□，f1oatJ), void (*rkqs)(float []float []int,float *float,float
758 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). Now, as far as Newton-Raphson is concerned, we are nearly in business. We want to find a vector value of V that zeros the vector value of F. We do this by invoking the globally convergent Newton’s method implemented in the routine newt of §9.7. Recall that the heart of Newton’s method involves solving the set of n2 linear 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 J has components given by Jij = ∂Fi ∂Vj (17.1.5) It is not feasible to compute these partial derivatives analytically. Rather, each requires a separate integration of the N ODEs, followed by the evaluation of ∂Fi ∂Vj ≈ Fi(V1,...,Vj + ∆Vj ,...) − Fi(V1,...,Vj ,...) ∆Vj (17.1.6) This is done automatically for you in the routine fdjac that comes with newt. The only input to newt that you have to provide is the routine vecfunc that calculates F by integrating the ODEs. Here is the appropriate routine, called shoot, that is to be passed as the actual argument in newt: #include "nrutil.h" #define EPS 1.0e-6 extern int nvar; Variables that you must define and set in your main proextern float x1,x2; gram. int kmax,kount; Communicates with odeint. float *xp,**yp,dxsav; void shoot(int n, float v[], float f[]) Routine for use with newt to solve a two point boundary value problem for nvar coupled ODEs by shooting from x1 to x2. Initial values for the nvar ODEs at x1 are generated from the n2 input coefficients v[1..n2], using the user-supplied routine load. The routine integrates the ODEs to x2 using the Runge-Kutta method with tolerance EPS, initial stepsize h1, and minimum stepsize hmin. At x2 it calls the user-supplied routine score to evaluate the n2 functions f[1..n2] that ought to be zero to satisfy the boundary conditions at x2. The functions f are returned on output. newt uses a globally convergent Newton’s method to adjust the values of v until the functions f are zero. The user-supplied routine derivs(x,y,dydx) supplies derivative information to the ODE integrator (see Chapter 16). The first set of global variables above receives its values from the main program so that shoot can have the syntax required for it to be the argument vecfunc of newt. { void derivs(float x, float y[], float dydx[]); void load(float x1, float v[], float y[]); void odeint(float ystart[], int nvar, float x1, float x2, float eps, float h1, float hmin, int *nok, int *nbad, void (*derivs)(float, float [], float []), void (*rkqs)(float [], float [], int, float *, float, float
17.1 The Shooting Method 759 float []float *float *void (*)(float,float []float []))) void rkgs(float y[],float dydx],int n,float *x, float htry,float eps,float yscal[],float *hdid,float *hnext, void(*derivs)(f1oat,f1oat[],f1oat☐))： void score(float xf,float y[],float f[]); int nbad,nok; float h1,hmin=0.0,*y; y=vector(1,nvar); kmax=0; h1=(x2-x1)/100.0: load(x1,v,y); odeint(y,nvar,x1,x2,EPS,h1,hmin,&nok,&nbad,derivs,rkqs); score(x2,y,f); free_vector(y,1,nvar); 83g granted for 19881992 11800 For some problems the initial stepsize AV might depend sensitively upon the initial conditions.It is straightforward to alter load to include a suggested stepsize from NUMERICAL RECIPESI h1 as another output variable and feed it to fdjac via a global variable. A complete cycle of the shooting method thus requires n2+1 integrations of the N coupled ODEs:one integration to evaluate the current degree of mismatch, and n2 for the partial derivatives.Each new cycle requires a new round of n2+1 integrations.This illustrates the enormous extra effort involved in solving two point 2a分旧 America computer, boundary value problems compared with initial value problems. If the differential equations are linear,then only one complete cycle is required, since (17.1.3)(17.1.4)should take us right to the solution.A second round can be useful,however,in mopping up some(never all)of the roundoff error. As given here,shoot uses the quality controlled Runge-Kutta method of 816.2 OF SCIENTIFIC to integrate the ODEs,but any of the other methods of Chapter 16 could just as well be used. You,the user,must supply shoot with:(i)a function load(x1,v,y)which calculates the n-vector y[1..n](satisfying the starting boundary conditions,of course),given the freely specifiable variables of v[1..n2]at the initial point x1; (ii)a function score(x2,y,f)which calculates the discrepancy vector f [1..n2] of the ending boundary conditions,given the vector y[1..n]at the endpoint x2: 10621 (iii)a starting vector v[1..n2];(iv)a function derivs for the ODE integration;and other obvious parameters as described in the header comment above. Fuunrggoirioh Numerical Recipes 43106 In $17.4 we give a sample program illustrating how to use shoot. (outside Software. CITED REFERENCES AND FURTHER READING: Ame ying of Acton,F.S.1970,Numerica/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.1 The Shooting Method 759 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). float [], float *, float *, void (*)(float, float [], float []))); void rkqs(float y[], float dydx[], int n, float *x, float htry, float eps, float yscal[], float *hdid, float *hnext, void (*derivs)(float, float [], float [])); void score(float xf, float y[], float f[]); int nbad,nok; float h1,hmin=0.0,*y; y=vector(1,nvar); kmax=0; h1=(x2-x1)/100.0; load(x1,v,y); odeint(y,nvar,x1,x2,EPS,h1,hmin,&nok,&nbad,derivs,rkqs); score(x2,y,f); free_vector(y,1,nvar); } For some problems the initial stepsize ∆V might depend sensitively upon the initial conditions. It is straightforward to alter load to include a suggested stepsize h1 as another output variable and feed it to fdjac via a global variable. A complete cycle of the shooting method thus requires n2 + 1 integrations of the N coupled ODEs: one integration to evaluate the current degree of mismatch, and n2 for the partial derivatives. Each new cycle requires a new round of n 2 + 1 integrations. This illustrates the enormous extra effort involved in solving two point boundary value problems compared with initial value problems. If the differential equations are linear, then only one complete cycle is required, since (17.1.3)–(17.1.4) should take us right to the solution. A second round can be useful, however, in mopping up some (never all) of the roundoff error. As given here, shoot uses the quality controlled Runge-Kutta method of §16.2 to integrate the ODEs, but any of the other methods of Chapter 16 could just as well be used. You, the user, must supply shoot with: (i) a function load(x1,v,y) which calculates the n-vector y[1..n] (satisfying the starting boundary conditions, of course), given the freely specifiable variables of v[1..n2] at the initial point x1; (ii) a function score(x2,y,f) which calculates 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 function derivs for 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: Mathematical Association of America). Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell)
760 Chapter 17.Two Point Boundary Value Problems 17.2 Shooting to a Fitting Point The shooting method described in 817.1 tacitly assumed that the "shots"would be able to traverse the entire domain of integration,even at the early stages of convergence to a correct solution.In some problems it can happen that,for very wrong starting conditions,an initial solution can't even get from I to x2 without encountering some 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 is feasible to integrate in the direction away from a singular point,using the special w 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 effort to obtain expansions of"wrong"solutions near the singular point(those not satisfying the desired boundary condition). 9 The solution to the above mentioned difficulties is shooting to a fitting point. Instead of integrating from to z2,we integrate first from i to some point rf that is between i and x2;and second from x2 (in the opposite direction)to xf. If(as before)the number of boundary conditions imposed at x1 is n,and the number imposed at 2 is n2,then there are n2 freely specifiable starting values at 至飞么 9 1 and nI freely specifiable starting values at x2.(If you are confused by this,go back to $17.1.)We can therefore define an n2-vector V()of starting parameters at 1,and a prescription load1(x1,v1,y)for mapping V(1)into a y that satisfies 61 the boundary conditions at 1, (c1)=(z1;Va)1,,V1)n2）i=1,,N (17.2.1) Likewise we can define an ni-vector V(2)of starting parameters at x2,and a prescription load2(x2,v2,y)for mapping V(2)into a y that satisfies the boundary Numerica 10621 conditions at x2, i(r2)=h(r2;V21,,V2)n) i=1,.,N (17.2.2) (outside We thus have a total of N freely adjustable parameters in the combination of V()and V(2).The N conditions that must be satisfied are that there be agreement in N components of y at xf between the values obtained integrating from one side and from the other, (xf;V(1)）=h(xf;V(2)i=1,,N (17.2.3) In some problems,the N matching conditions can be better described(physically mathematically,or numerically)by using N different functions F,i=1...N,each possibly depending on the N components y:.In those cases,(17.2.3)is replaced by F(xf;V(a)】=Fy(xf;V(2】i=1,,N (17.2.4)
760 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). 17.2 Shooting to a Fitting Point The shooting method described in §17.1 tacitly assumed that the “shots” would be able to traverse the entire domain of integration, even at the early stages of convergence to a correct solution. In some problems it can happen that, for very wrong starting conditions, an initial solution can’t even get from x1 to x2 without encountering some 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 is feasible to integrate in the direction away from 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 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. Instead of integrating from x1 to x2, we integrate first from x1 to some point xf that is between x1 and x2; and second from x2 (in the opposite direction) to xf . If (as before) the number of boundary conditions imposed at x 1 is n1, and the number imposed at x2 is n2, then there are n2 freely specifiable starting values at x1 and n1 freely specifiable starting values at x2. (If you are confused by this, go back to §17.1.) We can therefore define an n2-vector V(1) of starting parameters at x1, and a prescription load1(x1,v1,y) for mapping V(1) into a y that 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-vector V(2) of starting parameters at x2, and a prescription load2(x2,v2,y) for mapping V(2) into a y that 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 N freely adjustable parameters in the combination of V(1) and V(2). The N conditions that must be satisfied are that there be agreement in N components of y at xf between 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 N matching conditions can be better described (physically, mathematically, or numerically) by using N different functions Fi, i = 1 ...N, each possibly depending on the N components yi. In those cases, (17.2.3) is replaced by Fi[y(xf ; V(1))] = Fi[y(xf ; V(2))] i = 1,...,N (17.2.4)