f19-0

Sample 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).Chapter 19. Partial Differential Equations 19.0 Introduction The numerical treatment of partial differential equations is, by itself, a vast subject. Partial differential equations are at the heart of many, if not most, computer analyses or simulations of continuous physical systems, such as fluids,electromagnetic fields, the human body, and so on. The intent of this chapter is to givethebriefestpossibleusefulintroduction. Ideally,therewouldbeanentiresecond volumeof Numerical Recipes dealingwithpartialdifferentialequationsalone. (The references [1-4]provide, of course, available alternatives.) In most mathematics books, partial differentialequations (PDEs) are classified into the three categories, hyperbolic, parabolic, andelliptic, on the basis of their characteristics , or curves of information propagation. The prototypical example of a hyperbolic equation is the one-dimensional waveequation ∂2u ∂t2=v2∂2u ∂x2(19.0.1 ) where v=constant is thevelocityofwavepropagation. Theprototypicalparabolic equation is the diffusion equation ∂u ∂t=∂ ∂x/parenleftbigg D∂u ∂x/parenrightbigg (19.0.2 ) where Dis the diffusion coefficient. The prototypical elliptic equation is the Poissonequation ∂2u ∂x2+∂2u ∂y2=ρ(x, y)( 19.0.3 ) where the source term ρis given. If the source term is equal to zero, the equation isLaplace’s equation . Fromacomputationalpointofview,theclassificationintothesethreecanonical types is not very meaningful — or at least not as important as some other essentialdistinctions. Equations (19.0.1) and (19.0.2) both define initial value orCauchy problems: If information on u(perhaps including time derivative information) is 818 19.0Introduction 819Sample 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).. ... ... ... ... ... ... .. boundary conditions initial values (a) boundary values (b) Figure 19.0.1. Initial value problem (a) and boundary value problem (b) are contrasted. In (a) initial valuesaregivenonone “timeslice, ”anditisdesiredtoadvancethesolutionintime,computingsuccessive rows of open dots in the direction shown by the arrows. Boundary conditions at the left and right edges of each row ( ⊗) must also be supplied, but only one row at a time. Only one, or a few, previous rows need be maintained in memory. In (b), boundary values are speci fied around the edge of a grid, and an iterative process is employed to find the values of all the internal points (open circles). All grid points must be maintained in memory. given at some initial time t0for all x, then the equations describe how u(x, t) propagates itself forward in time. In other words, equations (19.0.1) and (19.0.2) describe time evolution. The goal of a numerical code should be to track that time evolution with some desired accuracy. Bycontrast,equation(19.0.3)directsusto findasingle “static”function u(x, y) which satis fies the equationwithin some (x, y)regionof interest, and which —one must also specify —has some desired behavior on the boundary of the region of interest. These problems are called boundary value problems . In general it is not 820 Chapter19. PartialDifferentialEquationsSample 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).possible stably to just “integrate in from the boundary ”in the same sense that an initial value problem can be “integrated forward in time. ”Therefore, the goal of a numericalcodeis somehowtoconvergeonthecorrectsolutioneverywhereat once. This, then, is the most important classi fication from a computational point of view: Is the problem at hand an initial value (time evolution) problem? or is it aboundary value (static solution) problem? Figure 19.0.1 emphasizes the distinction. Noticethatwhiletheitalicizedterminologyis standard,theterminology in parentheses is a much better description of the dichotomy from a computational perspective. The subclassi fication of initial value problems into parabolic and hyperbolicis much less important because (i) many actual problems are of a mixedtype, and (ii) as we will see, most hyperbolic problems get parabolic pieces mixed into them by the time one is discussing practical computationalschemes. InitialValue Problems An initial valueproblemis de finedby answers to the followingquestions: •What are the dependentvariables to be propagatedforward in time? •What is the evolution equation for each variable? Usually the evolution equations will all be coupled, with more than one dependent variable appearing on the right-hand side of each equation. •Whatisthehighesttimederivativethatoccursineachvariable ’sevolution equation? If possible, this time derivative should be put alone on the equation’s left-hand side. Not only the value of a variable, but also the value of all its time derivatives —up to the highest one —must be specified to define the evolution. •Whatspecialequations(boundaryconditions)governtheevolutionintime of points on the boundary of the spatial region of interest? Examples: Dirichletconditions specifythevaluesoftheboundarypointsasafunction oftime;Neumannconditions specifythevaluesofthenormalgradientson the boundary; outgoing-waveboundaryconditions arejust whattheysay. Sections 19.1 –19.3 of this chapter deal with initial value problems of several differentforms. We make no pretence of completeness, but rather hope to conveyacertain amount of generalizable information through a few carefully chosen model examples. These examples will illustrate an important point: One ’s principal computational concern must be the stabilityof the algorithm. Many reasonable- lookingalgorithmsforinitialvalueproblemsjustdon ’twork—theyarenumerically unstable. Boundary Value Problems The questions that de fine a boundary value problem are: •What are the variables? •What equationsare satis fied in the interior of the regionof interest? •What equations are satis fied by points on the boundary of the region of interest? (HereDirichletandNeumannconditionsarepossiblechoicesfor ellipticsecond-orderequations,butmorecomplicatedboundaryconditions can also be encountered.) 19.0Introduction 821Sample 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).In contrast to initial value problems, stability is relatively easy to achieve for boundary value problems. Thus, the efficiency of the algorithms, both in computationalload and storage requirements,becomes the principal concern. Because all the conditions on a boundary value problem must be satis fied “simultaneously, ”these problems usually boil down, at least conceptually, to the solutionoflargenumbersofsimultaneousalgebraicequations. Whensuchequations arenonlinear,theyareusuallysolvedbylinearizationanditeration;sowithoutmuch loss of generality we can view the problem as being the solution of special, large linear sets of equations. As an example, one which we will refer to in §§19.4–19.6 as our “model problem,”let us consider the solution of equation (19.0.3) by the finite-difference method. We representthe function u(x, y)by its valuesat the discreteset ofpoints xj=x0+j∆,j =0,1, ..., J yl=y0+l∆,l =0,1, ..., L(19.0.4 ) where ∆is thegrid spacing . From now on, we will write uj,lforu(xj,y l), and ρj,lforρ(xj,y l). For (19.0.3) we substitute a finite-difference representation (see Figure 19.0.2), uj+1 ,l−2uj,l+uj−1,l ∆2+uj,l+1−2uj,l+uj,l−1 ∆2=ρj,l (19.0.5 ) or equivalently uj+1 ,l+uj−1,l+uj,l+1+uj,l−1−4uj,l=∆2ρj,l (19.0.6 ) To write this system of linear equations in matrix form we need to make a vector out of u. Let us number the two dimensions of grid points in a single one-dimensional sequence by de fining i≡j(L+1 )+ lfor j=0,1, ..., J, l =0,1, ..., L (19.0.7 ) In other words, iincreases most rapidly along the columns representing yvalues. Equation (19.0.6) now becomes ui+L+1+ui−(L+1)+ui+1+ui−1−4ui=∆2ρi (19.0.8 ) This equation holds only at the interior points j=1,2, ..., J −1;l=1,2, ..., L−1. The points where j=0 j=J l=0 l=L[i.e.,i=0, ..., L ] [i.e.,i=J(L+1 ) , ..., J (L+1 )+ L] [i.e.,i=0,L+1, ..., J (L+1 ) ] [i.e.,i=L, L +1+ L, ..., J (L+1 )+ L](19.0.9 ) 822 Chapter19. PartialDifferentialEquationsSample 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).yL ∆ y1 y0x0 xJ x1...∆ A B Figure 19.0.2. Finite-difference representation of a second-order elliptic equation on a two-dimensional grid. Thesecond derivatives atthe point Aareevaluated usingthe points towhich Aisshownconnected. The second derivatives at point Bare evaluated using the connected points and also using “right-hand side”boundary information, shown schematically as ⊗. are boundary points where either uor its derivative has been speci fied. If we pull all this“known”information over to the right-hand side of equation (19.0.8), then the equation takes the form A·u=b (19.0.10 ) whereAhas the form shown in Figure 19.0.3. The matrix Ais called“tridiagonal with fringes. ”A general linear second-order elliptic equation a(x, y)∂2u ∂x2+b(x, y)∂u ∂x+c(x, y)∂2u ∂y2+d(x, y)∂u ∂y +e(x, y)∂2u ∂x∂y+f(x, y)u=g(x, y)(19.0.11 ) will lead to a matrix of similar structure except that the nonzero entries will not be constants. As a rough classi fication, there are three different approaches to the solution of equation (19.0.10), not all applicable in all cases: relaxation methods, “rapid” methods (e.g., Fourier methods), and direct matrix methods. 19.0Introduction 823Sample 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).−4 11 −4 •1 ••• • 1• −4 11 −4 J + 1 blocksincreasing iincreasing j J + 1 blocks1 1 • • • •−4 11 −4 •1 ••• ••• ••• • • •• ••• ••• ••• ••• • • •• ••• ••• • 1• −4 11 −4 −4 11 −4 •1 ••• • 1• −4 11 −4• • • • • • 1 • • • • 1• • • • 1 11 1 • • • 1 • • • • • •1 1 • • • • • • • • 1 1each block(L + 1) × (L + 1) Figure 19.0.3. Matrix structure derived froma second-order elliptic equation (here equation 19.0.6). All elements not shown are zero. The matrix has diagonal blocks that are themselves tridiagonal, and sub- and super-diagonal blocks that are diagonal. This form is called “tridiagonal with fringes. ”A matrix this sparse would never be stored in its full form as shown here. Relaxation methods make immediate use of the structure of the sparse matrix A. The matrix is split into two parts A=E−F (19.0.12 ) whereEis easily invertibleand Fis the remainder. Then (19.0.10)becomes E·u=F·u+b (19.0.13 ) The relaxation method involves choosing an initial guess u(0)and then solving successively for iterates u(r)from E·u(r)=F·u(r−1)+b (19.0.14 ) SinceEis chosen to be easily invertible, each iteration is fast. We will discuss relaxation methods in some detail in §19.5 and §19.6. 824 Chapter19. PartialDifferentialEquationsSample 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).So-called rapid methods [5]apply for only a rather special class of equations: those with constant coef ficients, or, more generally, those that are separable in the chosencoordinates. Inaddition,theboundariesmustcoincidewithcoordinatelines. This special class of equations is met quite often in practice. We defer detailed discussion to §19.4. Note, however,that the multigridrelaxationmethodsdiscussed in§19.6 can be faster than “rapid”methods. Matrix methods attempt to solve the equation A·x=b (19.0.15 ) directly. The degree to which this is practical depends very strongly on the exact structureofthematrix Afortheproblemathand,soourdiscussioncangonofarther than a few remarks and references at this point. Sparseness of the matrix mustbe the guiding force. Otherwise the matrix problem is prohibitively large. For example, the simplest problem on a 100×100 spatial gridwouldinvolve10000unknown uj,l’s, implyinga 10000 ×10000matrix A, containing 108elements! As we discussed at the end of §2.7, ifAis symmetric and positive de finite (as it usually is in elliptic problems), the conjugate-gradient algorithm can be used. In practice, rounding error often spoils the effectiveness of the conjugate gradientalgorithmforsolving finite-differenceequations. However,itisusefulwhen incorporatedin methodsthat first rewritethe equationsso that Ais transformedto a matrixA/primethat is close to the identity matrix. The quadratic surface de fined by the equations then has almost spherical contours, and the conjugate gradient algorithm works very well. In §2.7, in the routine linbcg, an analogous preconditioner was exploited for non-positivede finite problems with the more general biconjugate gradient method. For the positive de finite case that arises in PDEs, an example of a successful implementation is the incomplete Cholesky conjugate gradient method (ICCG)(see[6-8]). Anothermethodthatreliesonatransformationapproachisthe stronglyimplicit procedure of Stone [9]. A program called SIPSOL that implements this routine has been published [10]. A third class of matrix methods is the Analyze-Factorize-Operateapproach as described in §2.7. Generally speaking, when you have the storage available to implement these methods—not nearly as much as the 108above, but usually much more than is requiredbyrelaxationmethods —thenyoushouldconsiderdoingso. Onlymultigrid relaxation methods ( §19.6) are competitive with the best matrix methods. For grids larger than, say, 300×300, however, it is generally found that only relaxation methods, or “rapid”methods when they are applicable, are possible. There Is Moreto Life than FiniteDifferencing Besidesfinite differencing, there are other methods for solving PDEs. Most importantare finiteelement,MonteCarlo,spectral,andvariationalmethods. Unfor- tunately, we shall barely be able to do justice to finite differencing in this chapter, and so shall not be able to discuss these other methods in this book. Finite element methods [11-12]areoftenpreferredbypractitionersinsolidmechanicsandstructural 19.1Flux-ConservativeInitialValueProblems 825Sample 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).engineering; these methods allow considerable freedom in putting computational elementswhereyouwantthem,importantwhendealingwithhighlyirregulargeome-tries. Spectral methods [13-15]are preferredfor veryregulargeometriesand smooth functions;theyconvergemorerapidlythan finite-differencemethods(cf. §19.4),but they do not work well for problems with discontinuities. CITED REFERENCES AND FURTHER READING: Ames, W.F. 1977, Numerical Methods for Partial Differential Equations , 2nd ed. (New York: Academic Press). [1] Richtmyer, R.D., andMorton, K.W.1967, DifferenceMethods for InitialValue Problems ,2nded. (New York: Wiley-Interscience). [2] Roache, P.J. 1976, Computational Fluid Dynamics (Albuquerque: Hermosa). [3] Mitchell,A.R., andGriffiths, D.F.1980, TheFiniteDifferenceMethod inPartialDifferentialEqua- tions(New York: Wiley) [includes discussion of finite element methods]. [4] Dorr, F.W. 1970, SIAM Review , vol. 12, pp. 248–263. [5] Meijerink, J.A., and van der Vorst, H.A. 1977, Mathematics of Computation , vol. 31, pp. 148– 162. [6] van der Vorst, H.A. 1981, Journalof Computational Physics , vol. 44, pp. 1–19[review of sparse iterative methods]. [7] Kershaw, D.S. 1970, Journal of Computational Physics , vol. 26, pp. 43–65. [8] Stone, H.J. 1968, SIAM Journal on Numerical Analysis , vol. 5, pp. 530–558. [9] Jesshope, C.R. 1979, Computer Physics Communications , vol. 17, pp. 383–391. [10] Strang, G., and Fix, G. 1973, An Analysis of the Finite Element Method (Englewood Cliffs, NJ: Prentice-Hall). [11] Burnett, D.S. 1987, Finite Element Analysis: From Concepts to Applications (Reading, MA: Addison-Wesley). [12] Gottlieb, D. andOrszag, S.A. 1977, NumericalAnalysis of Spectral Methods: Theoryand Appli- cations(Philadelphia: S.I.A.M.). [13] Canuto, C., Hussaini, M.Y., Quarteroni, A., and Zang, T.A. 1988, Spectral Methods in Fluid Dynamics (New York: Springer-Verlag). [14] Boyd, J.P. 1989, Chebyshev and FourierSpectral Methods (New York: Springer-Verlag). [15] 19.1 Flux-Conservative InitialValue Problems Alargeclassofinitialvalue(time-evolution)PDEsinonespacedimensioncan be cast into the form of a flux-conservative equation , ∂u ∂t=−∂F(u) ∂x(19.1.1 ) whereuandFare vectors, and where (in some cases) Fmay depend not only on u but also on spatial derivativesof u. The vector Fis called the conserved flux. For example, the prototypical hyperbolic equation, the one-dimensional wave equation with constant velocity of propagation v ∂2u ∂t2=v2∂2u ∂x2(19.1.2 )