19.2 Diffusive Initial Value Problems 847 Woodward,P.,and Colella,P.1984,Journal of Computational Physics,vol.54,pp.115-173.[8] Rizzi,A..and Engquist,B.1987,Journal of Computational Physics,vol.72,pp.1-69.[9] 19.2 Diffusive Initial Value Problems Recall the model parabolic equation,the diffusion equation in one space dimension, Ou Ou .com or D 8t (19.2.1) Ox where D is the diffusion coefficient.Actually,this equation is a flux-conservative equation of the form considered in the previous section,with F=-DOu RECIPES (19.2.2) x 线 2 the flux in the z-direction.We will assume D>0.otherwise equation (19.2.1)has physically unstable solutions:A small disturbance evolves to become more and more Press. concentrated instead of dispersing.(Don't make the mistake of trying to find a stable differencing scheme for a problem whose underlying PDEs are themselves unstable!) Even though(19.2.1)is of the form already considered,it is useful to consider it as a model in its own right.The particular form of flux(19.2.2),and its direct OF SCIENTIFIC( generalizations,occur quite frequently in practice.Moreover,we have already seen that numerical viscosity and artificial viscosity can introduce diffusive pieces like 6 the right-hand side of(19.2.1)in many other situations. Consider first the case when D is a constant.Then the equation Ou 02u COMPUTING (ISBN 1992 di=D9 x2 (19.2.3) can be differenced in the obvious way: Numerical 10521 1- 43108 △t D u+1-2u+u3-1 (19.2.4) 天金 (△x)P (outside Recipes This is the FTCS scheme again,except that it is a second derivative that has been Software. differenced on the right-hand side.But this makes a world of difference!The FTCS North scheme was unstable for the hyperbolic equation;however,a quick calculation shows that the amplification factor for equation(19.2.4)is ξ=1- 4D△t k△ (④x2si2 2 (19.2.5) The requirement<1 leads to the stability criterion 2D△t (④x)Ps1 (19.2.6)
19.2 Diffusive Initial Value Problems 847 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). Woodward, P., and Colella, P. 1984, Journal of Computational Physics, vol. 54, pp. 115–173. [8] Rizzi, A., and Engquist, B. 1987, Journal of Computational Physics, vol. 72, pp. 1–69. [9] 19.2 Diffusive Initial Value Problems Recall the model parabolic equation, the diffusion equation in one space dimension, ∂u ∂t = ∂ ∂x D ∂u ∂x (19.2.1) where D is the diffusion coefficient. Actually, this equation is a flux-conservative equation of the form considered in the previous section, with F = −D ∂u ∂x (19.2.2) the flux in the x-direction. We will assume D ≥ 0, otherwise equation (19.2.1) has physically unstable solutions: A small disturbance evolves to become more and more concentrated instead of dispersing. (Don’t make the mistake of trying to find a stable differencing scheme for a problem whose underlying PDEs are themselves unstable!) Even though (19.2.1) is of the form already considered, it is useful to consider it as a model in its own right. The particular form of flux (19.2.2), and its direct generalizations, occur quite frequently in practice. Moreover, we have already seen that numerical viscosity and artificial viscosity can introduce diffusive pieces like the right-hand side of (19.2.1) in many other situations. Consider first the case when D is a constant. Then the equation ∂u ∂t = D ∂2u ∂x2 (19.2.3) can be differenced in the obvious way: un+1 j − un j ∆t = D un j+1 − 2un j + un j−1 (∆x)2 (19.2.4) This is the FTCS scheme again, except that it is a second derivative that has been differenced on the right-hand side. But this makes a world of difference! The FTCS scheme was unstable for the hyperbolic equation; however, a quick calculation shows that the amplification factor for equation (19.2.4) is ξ = 1 − 4D∆t (∆x)2 sin2 k∆x 2 (19.2.5) The requirement |ξ| ≤ 1 leads to the stability criterion 2D∆t (∆x)2 ≤ 1 (19.2.6)
848 Chapter 19.Partial Differential Equations The physical interpretation of the restriction (19.2.6)is that the maximum allowed timestep is,up to a numerical factor,the diffusion time across a cell of width△x. More generally,the diffusion time r across a spatial scale of size A is of order λ2 (19.2.7) Usually we are interested in modeling accurately the evolution of features with spatial scales A>Ax.If we are limited to timesteps satisfying (19.2.6),we will need to evolve through of order 2/(Ax)2 steps before things start to happen on the scale of interest.This number of steps is usually prohibitive.We must therefore find a stable way of taking timesteps comparable to,or perhaps-for accuracy- somewhat smaller than,the time scale of (19.2.7). This goal poses an immediate "philosophical"question.Obviously the large 运州 timesteps that we propose to take are going to be woefully inaccurate for the small scales that we have decided not to be interested in.We want those scales to do something stable,"innocuous,"and perhaps not too physically unreasonable.We RECIPESI want to build this innocuous behavior into our differencing scheme.What should f2a 2 it be? There are two different answers,each of which has its pros and cons.The first answer is to seek a differencing scheme that drives small-scale features to their equilibrium forms,e.g.,satisfying equation (19.2.3)with the left-hand side set to zero.This answer generally makes the best physical sense;but,as we will see,it leads 9 to a differencing scheme("fully implicit")that is only first-order accurate in time for 32.9么○0 the scales that we are interested in.The second answer is to let small-scale features OF SCIENTIFIC( maintain their initial amplitudes,so that the evolution of the larger-scale features of interest takes place superposed with a kind of "frozen in"(though fluctuating) background of small-scale stuff.This answer gives a differencing scheme("Crank- Nicolson")that is second-order accurate in time.Toward the end of an evolution calculation,however,one might want to switch over to some steps of the other kind, to drive the small-scale stuff into equilibrium.Let us now see where these distinct differencing schemes come from: Consider the following differencing of(19.2.3), Numerica 10.621 +1- [}-2+1+u出1 431 (19.2.8) (△x)2 (outside Recipes This is exactly like the FTCS scheme (19.2.4),except that the spatial derivatives on the right-hand side are evaluated at timestep n+1.Schemes with this character are North called fully implicit or backward time,by contrast with FTCS(which is called fully explicit).To solve equation(19.2.8)one has to solve a set of simultaneous linear equationsat each timestep fortheFortunately,this isasimple problem because the system is tridiagonal:Just group the terms in equation(19.2.8)appropriately: -aug++(1+2a)叫写+1-au+}=u吗， j=1,2J-1 (19.2.9) where D△t (△z2 (19.2.10)
848 Chapter 19. Partial Differential Equations 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). The physical interpretation of the restriction (19.2.6) is that the maximum allowed timestep is, up to a numerical factor, the diffusion time across a cell of width ∆x. More generally, the diffusion time τ across a spatial scale of size λ is of order τ ∼ λ2 D (19.2.7) Usually we are interested in modeling accurately the evolution of features with spatial scales λ ∆x. If we are limited to timesteps satisfying (19.2.6), we will need to evolve through of order λ2/(∆x)2 steps before things start to happen on the scale of interest. This number of steps is usually prohibitive. We must therefore find a stable way of taking timesteps comparable to, or perhaps — for accuracy — somewhat smaller than, the time scale of (19.2.7). This goal poses an immediate “philosophical” question. Obviously the large timesteps that we propose to take are going to be woefully inaccurate for the small scales that we have decided not to be interested in. We want those scales to do something stable, “innocuous,” and perhaps not too physically unreasonable. We want to build this innocuous behavior into our differencing scheme. What should it be? There are two different answers, each of which has its pros and cons. The first answer is to seek a differencing scheme that drives small-scale features to their equilibrium forms, e.g., satisfying equation (19.2.3) with the left-hand side set to zero. This answer generally makes the best physical sense; but, as we will see, it leads to a differencing scheme (“fully implicit”) that is only first-order accurate in time for the scales that we are interested in. The second answer is to let small-scale features maintain their initial amplitudes, so that the evolution of the larger-scale features of interest takes place superposed with a kind of “frozen in” (though fluctuating) background of small-scale stuff. This answer gives a differencing scheme (“CrankNicolson”) that is second-order accurate in time. Toward the end of an evolution calculation, however, one might want to switch over to some steps of the other kind, to drive the small-scale stuff into equilibrium. Let us now see where these distinct differencing schemes come from: Consider the following differencing of (19.2.3), un+1 j − un j ∆t = D un+1 j+1 − 2un+1 j + un+1 j−1 (∆x)2 (19.2.8) This is exactly like the FTCS scheme (19.2.4), except that the spatial derivatives on the right-hand side are evaluated at timestep n + 1. Schemes with this character are called fully implicit or backward time, by contrast with FTCS (which is called fully explicit). To solve equation (19.2.8) one has to solve a set of simultaneous linear equations at each timestep for the un+1 j . Fortunately, this is a simple problem because the system is tridiagonal: Just group the terms in equation (19.2.8) appropriately: −αun+1 j−1 + (1 + 2α)un+1 j − αun+1 j+1 = un j , j = 1, 2...J − 1 (19.2.9) where α ≡ D∆t (∆x)2 (19.2.10)
19.2 Diffusive Initial Value Problems 849 Supplemented by Dirichlet or Neumann boundary conditions at j=0 andj=J, equation(19.2.9)is clearly a tridiagonal system,which can easily be solved at each timestep by the method of $2.4. What is the behavior of(19.2.8)for very large timesteps?The answer is seen most clearly in (19.2.9),in the limit aoo(At-oo).Dividing by a,we see that the difference equations are just the finite-difference form of the equilibrium equation 02u =0 0x2 (19.2.11) What about stability?The amplification factor for equation(19.2.8)is 1 E- (19.2.12) 1 +4a sin2 Cam Clearly 1 for any stepsize At.The scheme is unconditionally stable.The details RECIPES I of the small-scale evolution from the initial conditions are obviously inaccurate for 令 large At.But,as advertised,the correct equilibrium solution is obtained.This is the characteristic feature of implicit methods. Here,on the other hand,is how one gets to the second of our above philosophical Press. answers,combining the stability of an implicit method with the accuracy of a method that is second-order in both space and time.Simply form the average of the explicit Progra and implicit FTCS schemes: OF SCIENTIFIC u吗+1-吗 D (u1-2u+1+ug+1)+(u+1-2u+u吗-1) △t △x)2 (19.2.13) Here both the left-and right-hand sides are centered at timestepn+,so the method (ISBN is second-order accurate in time as claimed.The amplification factor is Numerical 105211 1-2a sin2 2 Recipes 43108 k△ (19.2.14) 1+2a sin 2 (outside North Software. so the method is stable for any size At.This scheme is called the Crank-Nicolson scheme,and is our recommended method for any simple diffusion problem(perhaps supplemented by a few fully implicit steps at the end).(See Figure 19.2.1.) Now turn to some generalizations of the simple diffusion equation(19.2.3). Suppose first that the diffusion coefficient D is not constant,say D=D().We can adopt either of two strategies.First,we can make an analytic change of variable dx D(x) (19.2.15)
19.2 Diffusive Initial Value Problems 849 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). Supplemented by Dirichlet or Neumann boundary conditions at j = 0 and j = J, equation (19.2.9) is clearly a tridiagonal system, which can easily be solved at each timestep by the method of §2.4. What is the behavior of (19.2.8) for very large timesteps? The answer is seen most clearly in (19.2.9), in the limit α → ∞ (∆t → ∞). Dividing by α, we see that the difference equations are just the finite-difference form of the equilibrium equation ∂2u ∂x2 =0 (19.2.11) What about stability? The amplification factor for equation (19.2.8) is ξ = 1 1+4α sin2 k∆x 2 (19.2.12) Clearly |ξ| < 1 for any stepsize ∆t. The scheme is unconditionally stable. The details of the small-scale evolution from the initial conditions are obviously inaccurate for large ∆t. But, as advertised, the correct equilibrium solution is obtained. This is the characteristic feature of implicit methods. Here, on the other hand, is how one gets to the second of our above philosophical answers, combining the stability of an implicit method with the accuracy of a method that is second-order in both space and time. Simply form the average of the explicit and implicit FTCS schemes: un+1 j − un j ∆t = D 2 (un+1 j+1 − 2un+1 j + un+1 j−1 )+(un j+1 − 2un j + un j−1) (∆x)2 (19.2.13) Here both the left- and right-hand sides are centered at timestep n+ 1 2 , so the method is second-order accurate in time as claimed. The amplification factor is ξ = 1 − 2α sin2 k∆x 2 1+2α sin2 k∆x 2 (19.2.14) so the method is stable for any size ∆t. This scheme is called the Crank-Nicolson scheme, and is our recommended method for any simple diffusion problem (perhaps supplemented by a few fully implicit steps at the end). (See Figure 19.2.1.) Now turn to some generalizations of the simple diffusion equation (19.2.3). Suppose first that the diffusion coefficient D is not constant, say D = D(x). We can adopt either of two strategies. First, we can make an analytic change of variable y = dx D(x) (19.2.15)
850 Chapter 19.Partial Differential Equations 0 t or n FTCS x orj (a) .com or call granted fori 11-800-872 (including this one) to any server computer, 7423(North America tusers to make one paper 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN C: THE (b) Fully Implicit (c) Crank-Nicolson ART Figure 19.2.1. Three differencing schemes for diffusive problems (shown as in Figure 19.1.2).(a) 是 Forward Time Center Space is first-order accurate,but stable only for sufficiently small timesteps.(b)Fully Programs Implicit is stable for arbitrarily large timesteps,but is still only first-order accurate.(c)Crank-Nicolson is second-order accurate,and is usually stable for large timesteps. Then du D()Ox du (19.2.16 to dir 0t1 Ox becomes OF SCIENTIFIC COMPUTING(ISBN Ou 102u (19.2.17) 192 0t= D(y)ay2 and we evaluate D at the appropriateyHeuristically,the stability criterion(19.2.6) in an explicit scheme becomes FuurrgPgoglrion Numerical Recipes 10621 431085 Au)2 △t≤min (19.2.18) 2D1 (outside Note that constant spacing Ay in y does not imply constant spacing in z. North Software. An alternative method that does not require analytically tractable forms for Ame D is simply to difference equation(19.2.16)as it stands,centering everything appropriately.Thus the FTCS method becomes visit website machine 41-_D+12u+1-)-D,-2(-u- (19.2.19) △t (△x)2 where Dj+1/2三D(x5+1/2) (19.2.20)
850 Chapter 19. Partial Differential Equations 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). t or n x or j FTCS (a) (b) Fully Implicit (c) Crank-Nicolson Figure 19.2.1. Three differencing schemes for diffusive problems (shown as in Figure 19.1.2). (a) Forward Time Center Space is first-order accurate, but stable only for sufficiently small timesteps. (b) Fully Implicit is stable for arbitrarily large timesteps, but is still only first-order accurate. (c) Crank-Nicolson is second-order accurate, and is usually stable for large timesteps. Then ∂u ∂t = ∂ ∂xD(x) ∂u ∂x (19.2.16) becomes ∂u ∂t = 1 D(y) ∂2u ∂y2 (19.2.17) and we evaluate D at the appropriate yj . Heuristically, the stability criterion (19.2.6) in an explicit scheme becomes ∆t ≤ min j (∆y)2 2D−1 j (19.2.18) Note that constant spacing ∆y in y does not imply constant spacing in x. An alternative method that does not require analytically tractable forms for D is simply to difference equation (19.2.16) as it stands, centering everything appropriately. Thus the FTCS method becomes un+1 j − un j ∆t = Dj+1/2(un j+1 − un j ) − Dj−1/2(un j − un j−1) (∆x)2 (19.2.19) where Dj+1/2 ≡ D(xj+1/2) (19.2.20)
19.2 Diffusive Initial Value Problems 851 and the heuristic stability criterion is △t≤min (△x)21 (19.2.21) 2D3+1/2 The Crank-Nicolson method can be generalized similarly. The second complication one can consider is a nonlinear diffusion problem, for example where D=D(u).Explicit schemes can be generalized in the obvious way.For example,in equation (19.2.19)write 81 D+Ve=D+)+D别 (19.2.22) Implicit schemes are not as easy.The replacement(19.2.22)with n-n+1 leaves ICAL us with a nasty set of coupled nonlinear equations to solve at each timestep.Often 3 there is an easier way:If the form of D(u)allows us to integrate RECIPES dz=D(u)du (19.2.23) analytically for z(u),then the right-hand side of (19.2.1)becomes 2z/0x2,which we difference implicitly as 1-2+1+出 (△x)2 (19.2.24) 9超号 Now linearize each term on the right-hand side of equation(19.2.24),for example 梦*=2g+）=2吗)+g1-9) 0z (19.2.25) =z(u)+(u+1-u)D(u) This reduces the problem to tridiagonal form again and in practice usually retains 2鱼 Numerica 10621 the stability advantages of fully implicit differencing. 43126 (outside Recipes Schrodinger Equation Sometimes the physical problem being solved imposes constraints on the North differencing scheme that we have not yet taken into account.For example,consider the time-dependent Schrodinger equation of quantum mechanics.This is basically a parabolic equation for the evolution of a complex quantity For the scattering of a wavepacket by a one-dimensional potential V(),the equation has the form 0y-8b dxz+V()b (19.2.26) (Here we have chosen units so that Planck's constant h=1 and the particle mass m =1/2.)One is given the initial wavepacket,(t =0),together with boundary
19.2 Diffusive Initial Value Problems 851 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). and the heuristic stability criterion is ∆t ≤ min j (∆x)2 2Dj+1/2 (19.2.21) The Crank-Nicolson method can be generalized similarly. The second complication one can consider is a nonlinear diffusion problem, for example where D = D(u). Explicit schemes can be generalized in the obvious way. For example, in equation (19.2.19) write Dj+1/2 = 1 2 D(un j+1) + D(un j ) (19.2.22) Implicit schemes are not as easy. The replacement (19.2.22) with n → n + 1 leaves us with a nasty set of coupled nonlinear equations to solve at each timestep. Often there is an easier way: If the form of D(u) allows us to integrate dz = D(u)du (19.2.23) analytically for z(u), then the right-hand side of (19.2.1) becomes ∂ 2z/∂x2, which we difference implicitly as zn+1 j+1 − 2zn+1 j + zn+1 j−1 (∆x)2 (19.2.24) Now linearize each term on the right-hand side of equation (19.2.24), for example zn+1 j ≡ z(un+1 j ) = z(un j )+(un+1 j − un j ) ∂z ∂u j,n = z(un j )+(un+1 j − un j )D(un j ) (19.2.25) This reduces the problem to tridiagonal form again and in practice usually retains the stability advantages of fully implicit differencing. Schrodinger Equation ¨ Sometimes the physical problem being solved imposes constraints on the differencing scheme that we have not yet taken into account. For example, consider the time-dependent Schrodinger equation of quantum mechanics. This is basically a ¨ parabolic equation for the evolution of a complex quantity ψ. For the scattering of a wavepacket by a one-dimensional potential V (x), the equation has the form i ∂ψ ∂t = −∂2ψ ∂x2 + V (x)ψ (19.2.26) (Here we have chosen units so that Planck’s constant ¯h = 1 and the particle mass m = 1/2.) One is given the initial wavepacket, ψ(x, t = 0), together with boundary