f17-6

17.6HandlingInternalBoundaryConditionsorSingularPoints 775Sample 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 φ(x)is chosen by us. Written in terms of the mesh variable q, this equation is dx dq=ψ φ(x)(17.5.7 ) Notice that φ(x)should be chosen to be positive definite, so that the density of mesh points is everywhere positive. Otherwise (17.5.7) can have a zero in its denominator. To use automated mesh spacing, you add the three ODEs (17.5.5) and (17.5.7) to your set of equations, i.e., to the array y(j,k).N o w xbecomes a dependent variable! Qandψ also become new dependent variables. Normally, evaluating φrequires littleextra work since itwillbecomposed frompiecesofthe g’sthatexistanyway. Theautomatedprocedure allows one to investigate quickly how the numerical results might be affected by various strategiesformeshspacing. (Aspecialcaseoccurs ifthedesired meshspacing function Qcanbefound analytically, i.e., dQ/dxis directly integrable. Then, you need to add only two equations, those in 17.5.5, and two new variables x, ψ.) As an example of a typical strategy for implementing this scheme, consider a system with one dependent variable y(x). We could set dQ=dx ∆+|dlny| δ(17.5.8 ) or φ(x)=dQ dx=1 ∆+/vextendsingle/vextendsingle/vextendsingle/vextendsingledy/dx yδ/vextendsingle/vextendsingle/vextendsingle/vextendsingle(17.5.9 ) where ∆andδare constants that we choose. The first term would give a uniform spacing inxif it alone were present. The second term forces more grid points to be used where yis changing rapidly. The constants act to make every logarithmic change in yof an amount δ about as “attractive” to a grid point as a change in xof amount ∆. You adjust the constants according to taste. Other strategies are possible, such as alogarithmic spacing in x, replacing dxin the first term with dlnx. CITED REFERENCES AND FURTHER READING: Eggleton, P. P. 1971, Monthly Notices ofthe RoyalAstronomical Society , vol. 151,pp. 351–364. Kippenhan, R., Weigert, A., and Hofmeister, E. 1968, in Methods in Computational Physics , vol. 7 (New York: Academic Press), pp. 129ff. 17.6 Handling Internal Boundary Conditions or Singular Points Singularitiescanoccurintheinteriorsoftwopointboundary valueproblems. Typically, there is a point xsat which a derivative must be evaluated by an expression of the form S(xs)=N(xs,y) D(xs,y)(17.6.1 ) where the denominator D(xs,y)=0. In physical problems with finite answers, singular points usually come with their own cure: Where D→0, there the physical solution ymust be such as to make N→0simultaneously, insuch away thatthe ratiotakes on ameaningful value. This constraint on the solution yis often called a regularity condition . The condition thatD(xs,y)satisfy some special constraint at xsis entirely analogous to an extra boundary condition, an algebraic relation among the dependent variables that must hold at a point. We discussed a related situation earlier, in §17.2, when we described the “fitting point method” to handle the task of integrating equations with singular behavior at the boundaries.In those problems you are unable to integrate from one side of the domain to the other. 776 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).1 X XXXX XXXXXXX XXXXXX XXXXXX XXXXXX XXXXXX XXXXX XXXXX XXXXX XXXXX XXXXV VVVVVB BBBBB1 XXXXXX XXXXXXX XXXXXXX XXXXXXX XXXXX XXXXX XXXXX XXXXX XXXXV VVVVVV V VVVVVV VVVVVV X XXXXX XXXX111X 1X XXXX X 1X XXXX 1X XXXX11X X 1X X 1X X 111(b)B BBBBBB special block special block(a) BBBBBBBBBBBBB Figure 17.6.1. FDE matrix structure with an internal boundary condition. The internal condition introduces a special block. (a) Original form, compare with Figure 17.3.1; (b) final form, compare with Figure 17.3.2. However, the ODEs do have well-behaved derivatives and solutions in the neighborhood of the singularity, so it is readily possible to integrate away from the point. Both the relaxationmethod and the method of “shooting”to afitting point handle such problems easily. Also, in those problems the presence of singular behavior served to isolate some special boundaryvalues that had to be satis fied to solve the equations. The difference here is that we are concerned with singularities arising at intermediate points, where the location of the singular point depends on the solution, so is not known a priori. Consequently, we face a circular task: The singularity prevents us from finding a numerical solution, but we need a numerical solution to find its location. Such singularities are also associated with selecting a special value for some variable which allows the solution to satisfy the regularity condition at the singular point. Thus, internal singularities take onaspects of being internal boundary conditions. One way of handling internal singularities is to treat the problem as a free boundary problem, as discussed at the end of §17.0. Suppose, as a simple example, we consider the equation dy dx=N(x, y) D(x, y)(17.6.2 ) where NandDare required to pass through zero at some unknown point xs. We add the equation z≡xs−x1dz dx=0 ( 17.6.3 ) 17.6HandlingInternalBoundaryConditionsorSingularPoints 777Sample 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 xsis the unknown location of the singularity, and change the independent variable totby setting x−x1=tz, 0≤t≤1( 17.6.4 ) The boundary conditions at t=1become N(x, y)=0 ,D (x, y)=0 ( 17.6.5 ) Useofanadaptivemeshasdiscussed intheprevioussectionisanotherwaytoovercome the difficulties of an internal singularity. For the problem (17.6.2), we add the mesh spacing equations dQ dq=ψ (17.6.6 ) dψ dq=0 ( 17.6.7 ) with a simple mesh spacing function that maps xuniformly into q, where qruns from 1to M, the number of mesh points: Q(x)=x−x1,dQ dx=1 ( 17.6.8 ) Having added three first-order differential equations, we must also add their corresponding boundary conditions. If there were no singularity, these could simply be atq=1: x=x1,Q =0 ( 17.6.9 ) atq=M: x=x2 (17.6.10 ) and a total of Nvalues yispecified at q=1. In this case the problem is essentially an initial value problem with all boundary conditions speci fied at x1and the mesh spacing function is super fluous. However, in the actual case at hand we impose the conditions atq=1: x=x1,Q =0 ( 17.6.11 ) atq=M:N(x, y)=0 ,D (x, y)=0 ( 17.6.12 ) andN−1values yiatq=1. The“missing”yiis to be adjusted, in other words, so as to make the solution go through the singular point in a regular (zero-over-zero) rather thanirregular( finite-over-zero)manner. Notice alsothatthese boundary conditions do notdirectly impose a value for x 2, which becomes an adjustable parameter that the code varies in an attempt to match the regularity condition. In this example the singularity occurred at a boundary, and the complication arose because the location of the boundary was unknown. In other problems we might wish tocontinue the integration beyond the internal singularity. For the example given above, wecouldsimplyintegratetheODEstothesingularpoint,thenasaseparateproblemrecommencethe integration from the singular point on as far we care to go. However, in other cases thesingularity occurs internally, but does not completely determine the problem: There are stillsome more boundary conditions to be satis fied further along in the mesh. Such cases present nodifficultyinprinciple,butdorequiresomeadaptationoftherelaxationcodegivenin §17.3. In effect all you need to do is to add a “special”block of equations at the mesh point where the internal boundary conditions occur, and do the proper bookkeeping. Figure 17.6.1 illustrates a concrete example where the overall problem contains 5 equationswith2boundary conditionsatthe firstpoint,one “internal”boundarycondition,and twofinal boundary conditions. The figure shows the structure of the overall matrix equations along the diagonal in the vicinity of the special block. In the middle of the domain, blockstypically involve 5 equations (rows) in 10 unknowns (columns). For each block prior to thespecial block, the initial boundary conditions provided enough information to zero the first two columns of the blocks. The five FDEs eliminate five more columns, and the final three columns need tobe stored forthe backsubstitution step (as described in §17.3). Tohandle the extra condition we break the normal cycle and add a special block with only one equation: 778 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 internal boundary condition. This effectively reduces the required storage of unreduced coefficientsbyonecolumnfortherestofthegrid,andallowsustoreducetozerothe firstthree columns of subsequent blocks. The subroutines red, pinvs, bksub can readily handle these cases with minor recoding, but each problem makes for a special case, and you willhave to make the modi fications as required. CITED REFERENCES AND FURTHER READING: London, R.A., and Flannery, B.P. 1982, Astrophysical Journal , vol. 258, pp. 260–269.