f17-5

774 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.5 Automated Allocation of Mesh Points In relaxation problems, you have to choose values for the independent variable at the mesh points. This is called allocating the grid or mesh. The usual procedure is to pick a plausible set of values and, if it works, to be content. If it doesn’t work, increasing thenumber of points usually cures the problem. If we know ahead of time where our solutions will be rapidly varying, we can put more gridpointsthereandlesselsewhere. Alternatively,wecansolvetheproblemfirstonauniformmesh and then examine the solution to see where we should add more points. We then repeatthe solution with the improved grid. The object of the exercise is to allocate points in such a way as to represent the solution accurately. It is also possible to automate the allocation of mesh points, so that it is done “dynamically” during the relaxation process. This powerful technique not only improvesthe accuracy of the relaxation method, but also (as we will see in the next section) allowsinternal singularities to be handled in quite a neat way. Here we learn how to accomplishthe automatic allocation. We want to focus attention on the independent variable x, and consider two alternative reparametrizations of it. The first, we term q; this is just the coordinate corresponding to the mesh points themselves, so that q=1atk=1,q=2atk=2, and so on. Between any two meshpointswehave ∆q=1. Inthechangeofindependent variableintheODEsfrom xtoq, dy dx=g (17.5.1 ) becomes dy dq=gdx dq(17.5.2 ) In terms of q, equation (17.5.2) as an FDE might be written yk−yk−1−1 2/bracketleftBigg/parenleftBigg gdx dq/parenrightBigg k+/parenleftBigg gdx dq/parenrightBigg k−1/bracketrightBigg =0 ( 17.5.3 ) or some related version. Note that dx/dqshould accompany g. The transformation between xandqdepends only on the Jacobian dx/dq. Its reciprocal dq/dxis proportional to the density of mesh points. Now, given the function y(x), or its approximation at the current stage of relaxation, we are supposed to have some idea of how we want to specify the density of mesh points.For example, we might want dq/dxto be larger where yis changing rapidly, or near to the boundaries, or both. In fact, we can probably make up a formula for what we would like dq/dxtobeproportionalto. Theproblemisthatwedonotknowtheproportionalityconstant. That is, the formula that we might invent would not have the correct integral over the wholerangeof xsoastomake qvaryfrom 1toM,accordingtoitsdefinition. Tosolvethisproblem we introduce a second reparametrization Q(q), where Qis a new independent variable. The relation between Qandqis taken to be linear, so that a mesh spacing formula for dQ/dx differs only in its unknown proportionality constant. A linear relation implies d 2Q dq2=0 ( 17.5.4 ) or, expressed in the usual manner as coupled first-order equations, dQ(x) dq=ψdψ dq=0 ( 17.5.5 ) where ψis a new intermediate variable. We add these two equations to the set of ODEs being solved. Completing the prescription, we add a third ODE that is just our desired mesh-density function, namely φ(x)=dQ dx=dQ dqdq dx(17.5.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.