f17-3

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).ODEs at x1 (x2)are generated fromthe n2 (n1)coefficients v1 (v2),using the user- supplied routine load1 (load2) . The coefficients v1andv2should be stored in a sin- gle array v(1:n1+n2) in the main program by an EQUIVALENCE statement of the form (v1(1),v(1)),(v2(1),v(n2+1)) . Theinputparameter n=n1+n2 =nvar. Therou- tine integrates the ODEs to xfusing the Runge-Kutta method with tolerance EPS,i n i t i a l stepsize h1,andminimumstepsize hmin.Atxfitcallstheuser-suppliedsubroutine score toevaluatethe nvarfunctions f1andf2thatoughttomatchat xf. Thedifferences fare returned on output. newtuses a globally convergent Newton’s method to adjust the val- ues of vuntil the functions fare zero. The user-supplied subroutine 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) Path from x1toxfwith best trial values 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) Path from x2toxfwith best trial values 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 ) 754 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).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 Ncoupled first-order ODEs represented by FDEs on a mesh of Mpoints, a solution consists of values for Ndependent functions given at each of the Mmesh points, or N×Mvariables in all. The relaxation method determines the solution by starting with a guess and improving it, iteratively. As theiterations improve the solution, the result is said to relaxto the true solution. While several iteration schemes are possible, for most problems our old standby, multi- dimensional Newton’s method, works well. The method produces a matrix equation thatmust be solved, but the matrix takes a special, “block diagonal” form, that allows it to beinverted far more economically both in time and storage than would be possible for a generalmatrix of size (MN)×(MN). Since MNcan 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. LetusdevelopageneralsetofalgebraicequationsthatrepresenttheODEsbyFDEs. The ODEproblemisexactlyidenticaltothatexpressedinequations(17.0.1)–(17.0.3)wherewehadNcoupled first-order equations that satisfy n 1boundary conditions at x1andn2=N−n1 boundary conditions at x2. We first define a mesh or grid by a set of k=1,2, ..., Mpoints at which we supply values for the independent variable xk. In particular, x1is the initial boundary, and xMis the final boundary. We use the notation ykto refer to the entire set of dependent variables y1,y2,...,y Nat point xk. At an arbitrary point kin the middle of the mesh, we approximate the set of Nfirst-order ODEs by algebraic relations of the form 0=Ek≡yk−yk−1−(xk−xk−1)gk(xk,x k−1,yk,yk−1),k =2,3,...,M (17.3.3 ) The notation signifies that gkcan be evaluated using information from both points k,k−1. TheFDEslabeledby Ekprovide Nequationscoupling 2Nvariablesatpoints k, k−1. There areM−1points, k=2,3,...,M,atwhich difference equations oftheform (17.3.3) apply. ThustheFDEsprovideatotalof (M−1)Nequationsforthe MNunknowns. Theremaining Nequations come from the boundary conditions. At the first boundary we have 0=E1≡B(x1,y1)( 17.3.4 ) while at the second boundary 0=EM+1≡C(xM,yM)( 17.3.5 ) The vectors E1andBhave only n1nonzero components, corresponding to the n1boundary conditions at x1. It will turn out to be useful to take these nonzero components to be the lastn1components. In other words, Ej,1/negationslash=0only for j=n2+1,n2+2,...,N.A t the other boundary, only the first n2components of EM+1andCare nonzero: Ej,M +1/negationslash=0 only for j=1,2,...,n 2. The “solution” of the FDE problem in (17.3.3)–(17.3.5) consists of a set of variables yj,k, the values of the Nvariables yjat the Mpoints xk. The algorithm we describe below requires an initial guess for the yj,k. We then determine increments ∆yj,ksuch that yj,k+∆yj,kis an improved approximation to the solution. Equations for the increments are developed by expanding the FDEsin first-orderTaylor series with respect to small changes ∆yk. At an interior point, k=2,3,...,Mthis gives: Ek(yk+∆yk,yk−1+∆yk−1)≈Ek(yk,yk−1) +N/summationdisplay n=1∂Ek ∂y n,k−1∆yn,k−1+N/summationdisplay n=1∂Ek ∂y n,k∆yn,k(17.3.6 ) Forasolutionwewanttheupdatedvalue E(y+∆y)tobezero,sothegeneralsetofequations at an interior point can be written in matrix form as N/summationdisplay n=1Sj,n∆yn,k−1+2N/summationdisplay n=N+1Sj,n∆yn−N,k=−Ej,k,j=1,2,...,N (17.3.7 ) 17.3RelaxationMethods 755Sample 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).where Sj,n=∂E j,k ∂y n,k−1,S j,n+N=∂E j,k ∂y n,k,n =1,2,...,N (17.3.8 ) The quantity Sj,nis an N×2Nmatrix at each point k. Each interior point thus supplies a block of Nequations coupling 2Ncorrections tothe solution variablesatthe points k, k−1. Similarly, the algebraic relations at the boundaries can be expanded in a first-order Taylor series for increments that improve the solution. Since E1depends only on y1,w e find at the first boundary: N/summationdisplay n=1Sj,n∆yn,1=−Ej,1,j=n2+1,n2+2,...,N (17.3.9 ) where Sj,n=∂E j,1 ∂y n,1,n =1,2,...,N (17.3.10 ) At the second boundary, N/summationdisplay n=1Sj,n∆yn,M=−Ej,M +1,j=1,2,...,n 2 (17.3.11 ) where Sj,n=∂E j,M +1 ∂y n,M,n =1,2,...,N (17.3.12 ) We thus have in equations (17.3.7)–(17.3.12) a set of linear equations to be solved for the corrections ∆y, iterating until the corrections are sufficiently small. The equations have a special structure, because each Sj,ncouples only points k,k−1. Figure 17.3.1 illustrates the typical structure of the complete matrix equation for the case of 5 variables and 4 meshpoints, with 3 boundary conditions at the first boundary and 2 at the second. The 3×5 block of nonzero entries in the top left-hand corner of the matrix comes from the boundarycondition S j,nat point k=1. The next three 5×10blocks are the Sj,nat the interior points, coupling variables at mesh points (2,1), (3,2), and (4,3). Finally we have the blockcorresponding to the second boundary condition. We can solve equations (17.3.7)–(17.3.12) for the increments ∆yusing a form of Gaussian elimination that exploits the special structure of the matrix to minimize the totalnumber of operations, and that minimizes storage of matrix coefficients by packing theelements in a special blocked structure. (You might wish to review Chapter 2, especially§2.2, if you are unfamiliar with the steps involved in Gaussian elimination.) Recall that Gaussian elimination consists of manipulating the equations by elementary operations suchas dividing rows of coefficients by a common factor to produce unity in diagonal elements,and adding appropriate multiples of other rows to produce zeros below the diagonal. Herewe take advantage of the block structure by performing a bit more reduction than in pureGaussian elimination, so that the storage of coefficients is minimized. Figure 17.3.2 showstheformthatwewishtoachievebyelimination,justpriortothebacksubstitutionstep. Onlyasmallsubsetofthereduced MN×MNmatrixelementsneedstobestoredastheelimination progresses. Once the matrix elements reach the stage in Figure 17.3.2, the solution followsquickly by a backsubstitution procedure. Furthermore, the entire procedure, except the backsubstitution step, operates only on one block of the matrix at a time. The procedure contains four types of operations: (1)partial reduction to zero of certain elements of a block using results from a previous step,(2) elimination of the square structure of the remaining block elements such that the squaresectioncontainsunityalongthediagonal,andzeroinoff-diagonalelements,(3)storageoftheremaining nonzero coefficients for use in later steps, and (4) backsubstitution. We illustratethe steps schematically by figures. Considertheblockofequationsdescribingcorrectionsavailablefromtheinitialboundary conditions. Wehave n 1equations for Nunknown corrections. We wishto transform the first 756 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).X XXXXXXXX XXXXXXXX XXXXXXXX XXXXXXXX XXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX XXXXXXX XXXXXXX XXXXXXX XXXXXXX XXXXXXV VVVVVVVVVVVVVVVVVVVB BBBBBBBBBBBBBBBBBBB Figure 17.3.1. Matrix structure of a set of linear finite-difference equations (FDEs) with boundary conditions imposed at both endpoints. Here Xrepresents a coef ficient of the FDEs, Vrepresents a component of the unknown solution vector, and Bis a component of the known right-hand side. Empty spaces represent zeros. The matrix equation is to be solved by a special form of Gaussian elimination. (See text for details.) 1 1 1X XX 1X XX 1 1 1 1X XXXX 1X XXXX 1 1 1 1X XXXX 1 1 1 1X XXXX 1X XXXX 1V VVVVVVVVVVVVVVVVVVVB BBBBBBBBBBBBBBBBBBBX XXXX 1 Figure 17.3.2. Target structure of the Gaussian elimination. Once the matrix of Figure 17.3.1 has been reduced to this form, the solution follows quickly by backsubstitution. 17.3RelaxationMethods 757Sample 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).block so that itsleft-hand n1×n1square section becomes unity along the diagonal, and zero in off-diagonal elements. Figure 17.3.3 shows the original and final form of the first block of the matrix. In the figure we designate matrix elements that are subject to diagonalization by“D”, and elements that will be altered by “A”;i nt h efinal block, elements that are stored are labeled by “S”. We get from start to finish by selecting in turn n1“pivot”elements from amongthe firstn1columns,normalizingthepivotrowsothatthevalueofthe “pivot”element is unity, and adding appropriate multiples of this row to the remaining rows so that theycontain zerosinthepivotcolumn. Inits finalform,the reduced block expresses valuesforthe corrections to the firstn 1variables at mesh point 1in terms of values for the remaining n2 unknown corrections at point 1, i.e., we now know what the firstn1elements are in terms of theremaining n2elements. Westoreonly the finalset of n2nonzero columns from theinitial block, plus the column for the altered right-hand side of the matrix equation. We must emphasize here an important detail of the method. To exploit the reduced storage allowed by operating on blocks, it is essential that the ordering of columns in the s matrix of derivatives be such that pivot elements can be found among the firstn1rows of the matrix. This means that the n1boundary conditions at the first point must contain some dependence onthe firstj=1,2,..., n1dependent variables, y(j,1). Ifnot,thentheoriginal square n1×n1subsection of the first block will appear to be singular, and the method will fail. Alternatively, we would have to allow the search for pivot elements to involve all N columns of the block, and this would require column swapping and far more bookkeeping.The code provides a simple method of reordering the variables, i.e., the columns of the s matrix, so that this can be done easily. End of important detail. Next consider the block of Nequations representing the FDEsthatdescribe therelation betweenthe 2Ncorrectionsatpoints2and1. Theelementsofthatblock,togetherwithresults fromthe previous step,areillustrated inFigure17.3.4. Note thatby adding suitablemultiplesof rows from the first block we can reduce to zero the firstn 1columns of the block (labeled by“Z”), and, to do so, we will need to alter only the columns from n1+1toNand the vector elementon theright-hand side. Ofthe remaining columns wecan diagonalize a squaresubsection of N×Nelements, labeled by “D”in thefigure. In the process we alter the final set of n 2+1columns, denoted “A”in thefigure. The second half of the figure shows the block when we finish operating on it,with the stored (n2+1 )×Nelements labeled by “S.” IfweoperateonthenextsetofequationscorrespondingtotheFDEscouplingcorrections atpoints3and2,weseethatthestateofavailableresultsandnewequationsexactlyreproducesthe situation described in the previous paragraph. Thus, we can carry out those steps againforeach block in turnthrough block M. Finallyon block M+1weencounter the remaining boundary conditions. Figure 17.3.5 shows the final block of n 2FDEs relating the Ncorrections for variables at mesh point M, together with the result of reducing the previous block. Again, we can first use the prior results to zero the firstn1columns of the block. Now, when we diagonalize the remaining square section, we strike gold: We get values for the finaln2corrections at mesh point M. With the final block reduced, the matrix has the desired form shown previously in Figure 17.3.2, and the matrix is ripe for backsubstitution. Starting with the bottom row andworking up towards the top, at each stage we can simply determine one unknown correctionin terms of known quantities. The subroutine solvdeorganizes the steps described above. The principal procedures used in the algorithm are performed by subroutines called internally by solvde. The subroutine redeliminates leading columns of the smatrix using results from prior blocks. pinvsdiagonalizesthesquaresubsectionof sandstoresunreducedcoef ficients. bksubcarries out the backsubstitution step. The user of solvdemust understand the calling arguments, as described below, and supply a subroutine difeq, called by solvde, that evaluates the smatrix for each block. Most of the arguments in the call to solvdehave already been described, but some require discussion. Array y(j,k)contains the initial guess for the solution, with jlabeling the dependent variables at mesh points k. The problem involves neFDEs spanning points k=1,..., m .nbboundary conditions apply at the first point k=1. The array indexv(j) establishesthecorrespondencebetweencolumnsofthe smatrix,equations(17.3.8),(17.3.10), 758 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).(a) (b)D DD 1 00D DD 0 10D DD 0 01A AA S SSA AA S SSV VV V VVA AA S SS Figure 17.3.3. Reduction process for the first (upper left) block of the matrix in Figure 17.3.1. (a) Original form of the block, (b) final form. (See text for explanation.) (a) 1 00 ZZZZZV VVVVVVVS SS AAAAA (b) 1 00000000 0100V VVVVVVVS SSSSSSS0 10 ZZZZZ0 01 ZZZZZS SS DDDDDS SS DDDDDD DDDDD DDDDD DDDDA AAAAA AAAA 0 10000000 0100000S SS 10000S SS 010000 00100 0001S SSSSS SSSS Figure 17.3.4. Reduction process for intermediate blocks of the matrix in Figure 17.3.1. (a) Original form, (b) final form. (See text for explanation.) (a) 0 00000 00000 00001 00000 10000 0100 ZZ0 0010 ZZ0 0001 ZZS SSSS DDS SSSS DDV VVVVVVS SSSS AA (b) 0 00000 00000 00001 00000 10000 0100000 0010000 000100S SSSS 10S SSSS 01V VVVVVVS SSSSSS Figure 17.3.5. Reduction process for the last (lower right) block of the matrix in Figure 17.3.1. (a) Original form, (b) final form. (See text for explanation.) 17.3RelaxationMethods 759Sample 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).and (17.3.12), and the dependent variables. As described above it is essential that the nb boundaryconditionsat k=1involvethedependentvariablesreferencedbythe firstnbcolumns of the smatrix. Thus, columns jof the smatrix can be ordered by the user in difeqto refer to derivatives with respect to the dependent variable indexv(j) . The subroutine only attempts itmaxcorrection cycles before returning, even if the solution has not converged. The parameters conv, slowc, scalv relate to convergence. Each inversion of thematrix produces corrections for nevariables at mmesh points. Wewant these tobecome vanishingly smallas theiterationsproceed, butwemustde fine ameasure for the size of corrections. This error “norm”is very problem speci fic, so the user might wish to rewrite this section of the code as appropriate. In the program below we compute a valuefor the average correction errby summing the absolute value of all corrections, weighted by a scale factor appropriate to each type of variable: err=1 m×nem/summationdisplay k=1ne/summationdisplay j=1|∆Y(j,k)| scalv(j)(17.3.13 ) When err≤conv,themethodhasconverged. Notethattheusergetstosupplyanarray scalv which measures the typical size of each variable. Obviously, if erris large, we are far from a solution, and perhaps it is a bad idea to believe that the corrections generated from a first-order Taylor series are accurate. The number slowcmodulates application of corrections. After each iteration we apply only a fraction of the corrections found by matrix inversion: Y(j,k) →Y(j,k) +slowc max( slowc,err )∆Y(j,k) (17.3.14 ) Thus, when err>slowconly a fraction of the corrections are used, but when err≤slowc the entire correction gets applied. The call statement also supplies solvdewith the array y(1:nyj,1:nyk) containing the initial trial solution, and workspace arrays c(1:nci,1:ncj,1:nck) ,s(1:nsi,1:nsj) . The array cis the blockbuster: It stores the unreduced elements of the matrix built up for the backsubstitution step. If there are mmesh points, then there will be nck=m+1 blocks, each requiring nci=nerows and ncj=ne-nb+1 columns. Although large, this is small compared with (ne×m)2elements required for the whole matrix if we did not break it into blocks. We now describe the workings of the user-supplied subroutine difeq. The parameters of the subroutine are given by SUBROUTINE difeq(k,k1,k2,jsf,is1,isf,indexv,ne,s,nsi,nsj,y,nyj,nyk) The only information returned from difeqtosolvdeis the matrix of derivatives s(i,j); all other arguments are input to difeqand should not be altered. kindicates the currentmeshpoint,orblocknumber. k1,k2labelthefirstandlastpointinthemesh. If k=k1 ork>k2, the block involves the boundary conditions at the first orfinal points; otherwise the block acts on FDEs coupling variables at points k-1,k. The convention on storing information into the array s(i,j)follows that used in equations (17.3.8), (17.3.10), and (17.3.12): Rows ilabel equations, columns jrefer to derivatives with respect to dependent variables in the solution. Recall that each equation willdepend on the nedependent variables at either one or two points. Thus, jruns from 1to either neor2*ne. The column ordering for dependent variables at each point must agree with the list supplied in indexv(j) . Thus, for a block not at a boundary, the first column multiplies ∆Y(l=indexv(1),k-1 ),andthecolumn ne+1multiplies ∆Y(l=indexv(1),k ). is1,isf give the numbers of the starting and finalrowsthat need to be filled in the smatrix for this block. jsflabels the column in which the difference equations E j,kof equations (17.3.3)–(17.3.5) are stored. Thus, −s(i,jsf) is the vector on the right-hand side of the matrix. The reason for the minus sign is that difeqsupplies the actual difference equation, Ej,k, not its negative. Note that solvdesupplies a value for jsfsuch that the difference equation isput inthe column just after allderivatives inthe smatrix. Thus, difeqexpects to find values entered into s(i,j)for rows is1≤i≤isfand1≤j≤jsf. 760 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).Finally, s(1:nsi,1:nsj) andy(1:nyj,1:nyk) supply difeqwith storage for sand the solution variables yfor this iteration. An example of how to use this routine is given in the next section. Many ideas in the following code are due to Eggleton [1]. SUBROUTINE solvde(itmax,conv,slowc,scalv,indexv,ne,nb,m, * y,nyj,nyk,c,nci,ncj,nck,s,nsi,nsj) INTEGER itmax,m,nb,nci,ncj,nck,ne,nsi,nsj, * nyj,nyk,indexv(nyj),NMAX REAL conv,slowc,c(nci,ncj,nck),s(nsi,nsj), * scalv(nyj),y(nyj,nyk) PARAMETER (NMAX=10) Largest expected value of ne. C USES bksub,difeq,pinvs,red Driverroutineforsolutionoftwopointboundaryvalueproblemsbyrelaxation. itmaxisthe maximum number of iterations. convis the convergence criterion (see text). slowccon- trolsthefractionofcorrections actuallyusedaftereachiteration. scalv(1:nyj) contains typical sizes for each dependent variable, used to weight errors. indexv(1:nyj) lists the columnorderingofvariablesusedtoconstructthematrix sofderivatives. (The nbboundary conditions at the first mesh point must contain some dependence on the first nbvariables listedin indexv.) Theprobleminvolves neequationsfor neadjustabledependentvariables at each point. At the first mesh point there are nbboundary conditions. There are atotal ofmmesh points. y(1:nyj,1:nyk) is the two-dimensional array that contains the initial guessforallthedependentvariablesateachmeshpoint. Oneachiteration,itisupdatedby the calculated correction. The arrays c(1:nci,1:ncj,1:nck) ,s(1:nsi,1:nsj) sup- ply dummy storage used by the relaxation code; the minimum dimensions must satisfy: nci=ne,ncj=ne-nb+1 ,nck=m+1,nsi=ne,nsj=2*ne+1 . INTEGER ic1,ic2,ic3,ic4,it,j,j1,j2,j3,j4,j5,j6,j7,j8, * j9,jc1,jcf,jv,k,k1,k2,km,kp,nvars,kmax(NMAX) REAL err,errj,fac,vmax,vz,ermax(NMAX)k1=1 S e tu pr o wa n dc o l u m nm a r k e r s . k2=m nvars=ne*mj1=1j2=nb j3=nb+1 j4=nej5=j4+j1 j6=j4+j2 j7=j4+j3j8=j4+j4j9=j8+j1 ic1=1 ic2=ne-nbic3=ic2+1ic4=ne jc1=1 jcf=ic3do 16it=1,itmax Primary iteration loop. k=k1 Boundary conditions at first point. call difeq(k,k1,k2,j9,ic3,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk)call pinvs(ic3,ic4,j5,j9,jc1,k1,c,nci,ncj,nck,s,nsi,nsj)do 11k=k1+1,k2 Finite difference equations at all point pairs. kp=k-1 call difeq(k,k1,k2,j9,ic1,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk)call red(ic1,ic4,j1,j2,j3,j4,j9,ic3,jc1,jcf,kp, * c,nci,ncj,nck,s,nsi,nsj) call pinvs(ic1,ic4,j3,j9,jc1,k,c,nci,ncj,nck,s,nsi,nsj) enddo 11 k=k2+1 Final boundary conditions. call difeq(k,k1,k2,j9,ic1,ic2,indexv,ne,s,nsi,nsj,y,nyj,nyk) call red(ic1,ic2,j5,j6,j7,j8,j9,ic3,jc1,jcf,k2, * c,nci,ncj,nck,s,nsi,nsj) call pinvs(ic1,ic2,j7,j9,jcf,k2+1,c,nci,ncj,nck,s,nsi,nsj) 17.3RelaxationMethods 761Sample 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).call bksub(ne,nb,jcf,k1,k2,c,nci,ncj,nck) Backsubstitution. err=0. do13j=1,ne Convergence check, accumulate average error. jv=indexv(j)errj=0. km=0 vmax=0.do 12k=k1,k2 Find point with largest error, foreach dependent variable. vz=abs(c(jv,1,k)) if(vz.gt.vmax) then vmax=vzkm=k endif errj=errj+vz enddo 12 err=err+errj/scalv(j) Note weighting for each dependent variable. ermax(j)=c(jv,1,km)/scalv(j) kmax(j)=km enddo 13 err=err/nvars fac=slowc/max(slowc,err) Reduce correction applied when error is large. do15j=1,ne Apply corrections. jv=indexv(j) do14k=k1,k2 y(j,k)=y(j,k)-fac*c(jv,1,k) enddo 14 enddo 15 write(*,100) it,err,fac Summaryofcorrectionsforthisstep. Pointwithlargest errorforeachvariablecanbemonitored bywrit-ing out kmaxand ermax.if(err.lt.conv) return enddo 16 pause ’itmax exceeded in solvde’ Convergence failed. 100 format(1x,i4,2f12.6) returnEND SUBROUTINE bksub(ne,nb,jf,k1,k2,c,nci,ncj,nck) INTEGER jf,k1,k2,nb,nci,ncj,nck,ne REAL c(nci,ncj,nck) Backsubstitution, used internally by solvde. INTEGER i,im,j,k,kp,nbfREAL xxnbf=ne-nb im=1 do 13k=k2,k1,-1 Userecurrence relationstoeliminateremainingdependences. if (k.eq.k1) im=nbf+1 Special handling of first point. kp=k+1 do12j=1,nbf xx=c(j,jf,kp)do 11i=im,ne c(i,jf,k)=c(i,jf,k)-c(i,j,k)*xx enddo 11 enddo 12 enddo 13 do16k=k1,k2 Reorder corrections to be in column 1. kp=k+1do 14i=1,nb c(i,1,k)=c(i+nbf,jf,k) enddo 14 do15i=1,nbf c(i+nb,1,k)=c(i,jf,kp) enddo 15 enddo 16 returnEND 762 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).SUBROUTINE pinvs(ie1,ie2,je1,jsf,jc1,k,c,nci,ncj,nck,s,nsi,nsj) INTEGER ie1,ie2,jc1,je1,jsf,k,nci,ncj,nck,nsi,nsj,NMAX REAL c(nci,ncj,nck),s(nsi,nsj) PARAMETER (NMAX=10) Diagonalize the square subsection of the smatrix, and store the recursion coefficients in c; used internally by solvde. INTEGER i,icoff,id,ipiv,irow,j,jcoff,je2,jp,jpiv,js1,indxr(NMAX)REAL big,dum,piv,pivinv,pscl(NMAX)je2=je1+ie2-ie1 js1=je2+1 do 12i=ie1,ie2 Implicit pivoting, as in §2.1. big=0.do 11j=je1,je2 if(abs(s(i,j)).gt.big) big=abs(s(i,j)) enddo 11 if(big.eq.0.) pause ’singular matrix, row all 0 in pinvs’ pscl(i)=1./big indxr(i)=0 enddo 12 do18id=ie1,ie2 piv=0. do14i=ie1,ie2 Find pivot element. if(indxr(i).eq.0) then big=0. do13j=je1,je2 if(abs(s(i,j)).gt.big) then jp=j big=abs(s(i,j)) endif enddo 13 if(big*pscl(i).gt.piv) then ipiv=ijpiv=jppiv=big*pscl(i) endif endif enddo 14 if(s(ipiv,jpiv).eq.0.) pause ’singular matrix in pinvs’ indxr(ipiv)=jpiv In place reduction. Save column ordering. pivinv=1./s(ipiv,jpiv)do 15j=je1,jsf Normalize pivot row. s(ipiv,j)=s(ipiv,j)*pivinv enddo 15 s(ipiv,jpiv)=1.do 17i=ie1,ie2 Reduce nonpivot elements in column. if(indxr(i).ne.jpiv) then if(s(i,jpiv).ne.0.) then dum=s(i,jpiv) do16j=je1,jsf s(i,j)=s(i,j)-dum*s(ipiv,j) enddo 16 s(i,jpiv)=0. endif endif enddo 17 enddo 18 jcoff=jc1-js1 Sort and store unreduced coefficients. icoff=ie1-je1do 21i=ie1,ie2 irow=indxr(i)+icoff do19j=js1,jsf c(irow,j+jcoff,k)=s(i,j) enddo 19 enddo 21 17.3RelaxationMethods 763Sample 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).return END SUBROUTINE red(iz1,iz2,jz1,jz2,jm1,jm2,jmf,ic1,jc1,jcf,kc, * c,nci,ncj,nck,s,nsi,nsj) INTEGER ic1,iz1,iz2,jc1,jcf,jm1,jm2,jmf,jz1,jz2,kc,nci,ncj, * nck,nsi,nsj REAL c(nci,ncj,nck),s(nsi,nsj) Reduce columns jz1-jz2of the smatrix,using previousresults asstored inthe cmatrix. Only columns jm1-jm2,jmf are affected by the prior results. redis used internally by solvde. INTEGER i,ic,j,l,loff REAL vxloff=jc1-jm1 ic=ic1 do 14j=jz1,jz2 Loop over columns to be zeroed. do12l=jm1,jm2 Loop over columns altered. vx=c(ic,l+loff,kc) do11i=iz1,iz2 Loop over rows. s(i,l)=s(i,l)-s(i,j)*vx enddo 11 enddo 12 vx=c(ic,jcf,kc)do 13i=iz1,iz2 Plus final element. s(i,jmf)=s(i,jmf)-s(i,j)*vx enddo 13 ic=ic+1 enddo 14 return END “AlgebraicallyDifficult” Sets ofDifferential Equations Relaxationmethods allowyoutotakeadvantage ofanadditional opportunity that,while not obvious, can speed up some calculations enormously. It is not necessary that the setof variables y j,kcorrespond exactly with the dependent variables of the original differential equations. They can be related to those variables through algebraic equations. Obviously, itis necessary only that the solution variables allow us to evaluatethe functions y,g,B,Cthat are used to construct the FDEs from the ODEs. In some problems gdepends on functions of ythatareknown only implicitly,so thatiterativesolutions arenecessary toevaluate functions in the ODEs. Often one can dispense with this “internal”nonlinear problem by de fining a new set of variables from which both y,gand the boundary conditions can be obtained directly. A typical example occurs in physical problems where the equations require solutionofacomplex equation ofstatethatcan be expressed inmore convenient termsusing variablesother than the original dependent variables in the ODE. While this approach is analogous toperforming an analyticchange of variables directly on the original ODEs, such an analytic transformation might be prohibitively complicated. The change of variables in the relaxationmethod is easy and requires no analytic manipulations. CITED REFERENCES AND FURTHER READING: Eggleton,P.P.1971, MonthlyNoticesoftheRoyalAstronomicalSociety ,vol.151,pp.351–364.[1] 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. 764 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).17.4 AWorkedExample: SpheroidalHarmonics The best way to understand the algorithms of the previous sections is to see them employedto solve an actual problem. As a sample problem, we have selectedthe computation of spheroidal harmonics. (The more common name is spheroidal angle functions, but we prefer the explicit reminder of the kinship with spherical harmonics.) We will show how to find spheroidal harmonics, first by the method of relaxation ( §17.3), and then by the methods of shooting ( §17.1) and shooting to afitting point ( §17.2). Spheroidal harmonics typically arise when certain partial differential equations are solved by separation of variables in spheroidal coordinates. They satisfy the following differential equation on the interval −1≤x≤1: d dx/bracketleftbigg (1−x2)dS dx/bracketrightbigg +/parenleftbigg λ−c2x2−m2 1−x2/parenrightbigg S=0 ( 17.4.1 ) Heremisaninteger, cisthe“oblatenessparameter, ”andλistheeigenvalue. Despite the notation, c2can be positive or negative. For c2>0the functions are called “prolate,”while if c2<0they are called “oblate.”The equationhas singular points atx=±1andistobesolvedsubjecttotheboundaryconditionsthatthesolutionbe regularat x=±1. Onlyforcertainvaluesof λ,theeigenvalues,willthisbepossible. Ifweconsider firstthesphericalcase,where c=0,werecognizethedifferential equation for Legendre functions Pm n(x). In this case the eigenvalues are λmn = n(n+1 ),n=m, m +1,.... The integer nlabels successive eigenvalues for fixedm: When n=mwe have the lowest eigenvalue, and the corresponding eigenfunctionhas no nodes in the interval −1<x< 1; when n=m+1we have the nexteigenvalue,andthe eigenfunctionhas onenodeinside (−1,1); andso on. A similar situation holdsforthe generalcase c2/negationslash=0. We write the eigenvalues of (17.4.1) as λmn(c)and the eigenfunctions as Smn(x;c).F o rfixedm,n= m, m +1,...labels the successive eigenvalues. Thecomputationof λmn(c)andSmn(x;c)traditionallyhasbeenquitedif ficult. Complicated recurrence relations, power series expansions, etc., can be found in references [1-3]. Cheap computing makes evaluation by direct solution of the differential equation quite feasible. Thefirst step is to investigate the behavior of the solution near the singular points x=±1. Substituting a power series expansion of the form S=( 1±x)α∞/summationdisplay k=0ak(1±x)k(17.4.2 ) in equation (17.4.1), we find that the regular solution has α=m/ 2. (Without loss of generality we can take m≥0since m→−mis a symmetry of the equation.) We get an equationthat is numericallymore tractable if we factor out this behavior.Accordingly we set S=( 1−x 2)m/ 2y (17.4.3 ) We thenfind from (17.4.1) that ysatisfies the equation (1−x2)d2y dx2−2(m+1 )xdy dx+(µ−c2x2)y=0 ( 17.4.4 )