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)
17.2 Shooting to a Fitting Point 761 In the program below,the user-supplied function score(xf,y,f)is supposed to map an input N-vector y into an output N-vector F.In most cases,you can dummy this function as the identity mapping Shooting to a fitting point uses globally convergent Newton-Raphson exactly as in 817.1.Comparing closely with the routine shoot of the previous section,you should have no difficulty in understanding the following routine shootf.The main differences in use are that you have to supply both load1 and load2.Also,in the calling program you must supply initial guesses for v1 [1..n2]and v2[1..n1]. Once again a sample program illustrating shooting to a fitting point is given in $17.4. #include "nrutil.h" #define EPS 1.0e-6 granted for 19881992 extern int nn2,nvar; Variables that you must define and set in your main pro- extern float x1,x2,xf; gram. 1.800 int kmax,kount; Communicates with odeint. float *xp,**yp,dxsavi to any /Cambridge from NUMERICAL RECIPESI void shootf(intn,f1oatv[☐，f1oatf0) Routine for use with newt to solve a two point boundary value problem for nvar coupled (Nort to make ODEs by shooting from x1 and x2 to a fitting point xf.Initial values for the nvar ODEs at x1 (x2)are generated from the n2 (n1)coefficients v1 (v2),using the user-supplied routine server computer, e University Press. THE load1 (load2).The coefficients v1 and v2 should be stored in a single array v[1..n1+n2] America in the main program by statements of the form v1=v;and v2 &v [n2];.The input param- one paper ART eter n =n1 +n2 nvar.The routine integrates the ODEs to xf using the Runge-Kutta method with tolerance EPS,initial stepsize h1,and minimum stepsize hmin.At xf it calls the user-supplied routine score to evaluate the nvar functions f1 and f2 that ought to match Programs at xf.The differences f are returned on output.newt uses a globally convergent Newton's for their 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.Set nn2 =n2 in the main program. f void derivs(f1oatx,f1oaty[],f1 oat dydx▣)； OF SCIENTIFIC COMPUTING(ISBN void load1(float x1,float vi[],float y]); 18881920 void1oad2(f1oatx2,f1oatv2[],f1oaty☐)； void odeint(float ystart[],int nvar,float x1,float x2. float eps, float hi,float hmin,int *nok,int *nbad, 10-621 void(*derivs)(f1oat,f1oat],f1oat[☐)， void(*rkgs)(f1oat☐，float☐，int,f1oat*,f1oat,f1oat, .Further reproduction, float []float *float *void (*)(float,float []float []))) Numerical Recipes -43108 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 [])) (outside void score(float xf,float y[],float f☐)； int i,nbad,nok; North Software. float h1,hmin=0.0,*f1,*f2,*yi ying of fi-vector(1,nvar); f2=vector(1,nvar); y=vector(1,nvar); kmax=0; h1=(x2-x1)/100.0: loadi(x1,v,y); Path from x1 to xf with best trial values v1. odeint(y,nvar,x1,xf,EPS,h1,hmin,&nok,&nbad,derivs,rkqs); score(xf,y,f1); load2(x2,&v [nn2],y); Path from x2 to xf with best trial values v2. odeint(y,nvar,x2,xf,EPS,h1,hmin,&nok,&nbad,derivs,rkqs); score(xf,y,f2);
17.2 Shooting to a Fitting Point 761 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). In the program below, the user-supplied function score(xf,y,f) is supposed to map an input N-vector y into an output N-vector F. In most cases, you can dummy this function as the identity mapping. Shooting to a fitting point uses globally convergent Newton-Raphson exactly as in §17.1. Comparing closely with the routine shoot of the previous section, you should have no difficulty in understanding the following routine shootf. The main differences in use are that you have to supply both load1 and load2. Also, in the calling program you must supply initial guesses for v1[1..n2] and v2[1..n1]. Once again a sample program illustrating shooting to a fitting point is given in §17.4. #include "nrutil.h" #define EPS 1.0e-6 extern int nn2,nvar; Variables that you must define and set in your main proextern float x1,x2,xf; gram. int kmax,kount; Communicates with odeint. float *xp,**yp,dxsav; void shootf(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 and x2 to a fitting point xf. Initial values for the nvar ODEs at x1 (x2) are generated from the n2 (n1) coefficients v1 (v2), using the user-supplied routine load1 (load2). The coefficients v1 and v2 should be stored in a single array v[1..n1+n2] in the main program by statements of the form v1=v; and v2 = &v[n2];. The input parameter n = n1 + n2 = nvar. The routine integrates the ODEs to xf using the Runge-Kutta method with tolerance EPS, initial stepsize h1, and minimum stepsize hmin. At xf it calls the user-supplied routine score to evaluate the nvar functions f1 and f2 that ought to match at xf. The differences 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. Set nn2 = n2 in the main program. { void derivs(float x, float y[], float dydx[]); void load1(float x1, float v1[], float y[]); void load2(float x2, float v2[], 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, 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 i,nbad,nok; float h1,hmin=0.0,*f1,*f2,*y; f1=vector(1,nvar); f2=vector(1,nvar); y=vector(1,nvar); kmax=0; h1=(x2-x1)/100.0; load1(x1,v,y); Path from x1 to xf with best trial values v1. odeint(y,nvar,x1,xf,EPS,h1,hmin,&nok,&nbad,derivs,rkqs); score(xf,y,f1); load2(x2,&v[nn2],y); Path from x2 to xf with best trial values v2. odeint(y,nvar,x2,xf,EPS,h1,hmin,&nok,&nbad,derivs,rkqs); score(xf,y,f2);
762 Chapter 17.Two Point Boundary Value Problems for(i=1;i<=n;i+)f[i]=f1[i]-f2[1]: free_vector(y,1,nvar); free_vector(f2,1,nvar); free_vector(f1,1,nvar); There are boundary value 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,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). Stoer.J.,and Bulirsch.R.1980.Introduction to Numerical Analysis(New York:Springer-Verlag). 397.3.5-7.3.6.[1] RECIPES 17.3 Relaxation Methods 6 Press. In relaxation 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 IENTIFIC dy =g(x,y) (17.3.1) dx 6 with an algebraic equation relating function values at two points k,k-1: 张-k-1-(xk-xk-1)g[2(xk+xk-1,(十张-1】=0 (17.3.2) The form of the FDE in (17.3.2)illustrates the idea,but not uniquely:There are many ways to turn the ODE into an FDE.When the problem involves N coupled first-order ODEs represented by FDEs on a mesh of M points,a solution consists of values for N dependent 10621 functions given at each of the M mesh points,or Nx M variables in all.The relaxation Numerica method determines the solution by starting with a guess and improving it,iteratively.As the 43106 iterations improve the solution,the result is said to relax to the true solution. While several iteration schemes are possible,for most problems our old standby,multi- Recipes dimensional Newton's method,works well.The method produces a matrix equation that must be solved,but the matrix takes a special,"block diagonal"form,that allows it to be inverted far more economically both in time and storage than would be possible for a general North Software. matrix of size (MN)x (MN).Since MN can easily be several thousand,this is crucial for the feasibility of the method Our implementation couples at most pairs of points,as in equation (17.3.2).More points can be coupled,but then the method becomes more complex. We will provide enough background so that you can write a more general scheme if you have the patience to do so. Let us develop a general set ofalgebraic equations that represent the ODEs by FDEs.The ODE problem is exactly identical to that expressed in equations(17.0.1)(17.0.3)where we had N coupled first-order equations that satisfy n boundary conditions at and n2=N-n boundary conditions at 2.We first define a mesh or grid by a set ofk=1,2,...,M points at which we supply values for the independent variable r.In particular,1 is the initial boundary,and is the final boundary.We use the notation y to refer to the entire set of
762 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). for (i=1;i<=n;i++) f[i]=f1[i]-f2[i]; free_vector(y,1,nvar); free_vector(f2,1,nvar); free_vector(f1,1,nvar); } There are boundary value 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: Mathematical Association of America). Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell). Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §§7.3.5–7.3.6. [1] 17.3 Relaxation Methods In relaxation 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 1 2 (xk + xk−1), 1 2 (yk + yk−1) =0 (17.3.2) The form of the FDE in (17.3.2) illustrates the idea, but not uniquely: There are many ways to turn the ODE into an FDE. When the problem involves N coupled first-order ODEs represented by FDEs on a mesh of M points, a solution consists of values for N dependent functions given at each of the M mesh points, or N × M variables in all. The relaxation method determines the solution by starting with a guess and improving it, iteratively. As the iterations improve the solution, the result is said to relax to the true solution. While several iteration schemes are possible, for most problems our old standby, multidimensional Newton’s method, works well. The method produces a matrix equation that must be solved, but the matrix takes a special, “block diagonal” form, that allows it to be inverted far more economically both in time and storage than would be possible for a general matrix of size (MN) × (MN). Since MN can easily be several thousand, this is crucial for the feasibility of the method. Our implementation couples at most pairs of points, as in equation (17.3.2). More points can be coupled, but then the method becomes more complex. We will provide enough background so that you can write a more general scheme if you have the patience to do so. Let us develop a general set of algebraic equations that represent the ODEs by FDEs. The ODE problem is exactly identical to that expressed in equations (17.0.1)–(17.0.3) where we had N coupled first-order equations that satisfy n1 boundary conditions at x1 and n2 = N − n1 boundary conditions at x2. We first define a mesh or grid by a set of k = 1, 2, ..., M points at which we supply values for the independent variable xk. In particular, x1 is the initial boundary, and xM is the final boundary. We use the notation yk to refer to the entire set of