19.4 Fourier and Cyclic Reduction Methods 857 For example,a combined advective-diffusion equation,such as Ou oi=-v 02u +D (19.3.21) Ox 0x2 might profitably use an explicit scheme for the advective term combined with a Crank-Nicolson or other implicit scheme for the diffusion term. The alternating-direction implicit (ADI)method,equation(19.3.16),is an example of operator splitting with a slightly different twist.Let us reinterpret (19.3.19)to have a different meaning:Let i now denote an updating method that 8 includes algebraically all the pieces of the total operator C,but which is desirably stable only for the Ci piece;likewise U2,...Um.Then a method of getting from un to un+1 is B un+1/m=h(u”,△t/m) ICAL un+2/m =U(un+1/m,At/m) (19.3.22) un+1=Um(un+(m-1)/m,At/m) 8 9 The timestep for each fractional step in(19.3.22)is now only 1/m ofthe full timestep, because each partial operation acts with all the terms of the original operator. Equation(19.3.22)is usually,though not always,stable as a differencing scheme 》 9 for the operator C.In fact,as a rule of thumb,it is often sufficient to have stable i's only for the operator pieces having the highest number of spatial derivatives-the other i's can be unstable-to make the overall scheme stable! 6 It is at this point that we turn our attention from initial value problems to boundary value problems.These will occupy us for the remainder of the chapter. CITED REFERENCES AND FURTHER READING: Ames,W.F.1977,Numerical Methods for Partial Differential Equations,2nd ed.(New York: Academic Press). 3四三a Numerica 10.621 4312 (outside Recipes 19.4 Fourier and Cyclic Reduction Methods for North Boundary Value Problems As discussed in $19.0,most boundary value problems(elliptic equations,for example)reduce to solving large sparse linear systems of the form A·u=b (19.4.1) either once,for boundary value equations that are linear,or iteratively,for boundary value equations that are nonlinear
19.4 Fourier and Cyclic Reduction Methods 857 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). For example, a combined advective-diffusion equation, such as ∂u ∂t = −v ∂u ∂x + D ∂2u ∂x2 (19.3.21) might profitably use an explicit scheme for the advective term combined with a Crank-Nicolson or other implicit scheme for the diffusion term. The alternating-direction implicit (ADI) method, equation (19.3.16), is an example of operator splitting with a slightly different twist. Let us reinterpret (19.3.19) to have a different meaning: Let U1 now denote an updating method that includes algebraically all the pieces of the total operator L, but which is desirably stable only for the L1 piece; likewise U2,... Um. Then a method of getting from un to un+1 is un+1/m = U1(un, ∆t/m) un+2/m = U2(un+1/m, ∆t/m) ··· un+1 = Um(un+(m−1)/m, ∆t/m) (19.3.22) The timestep for each fractional step in (19.3.22) is now only 1/m of the full timestep, because each partial operation acts with all the terms of the original operator. Equation (19.3.22) is usually, though not always, stable as a differencing scheme for the operator L. In fact, as a rule of thumb, it is often sufficient to have stable U i’s only for the operator pieces having the highest number of spatial derivatives — the other Ui’s can be unstable — to make the overall scheme stable! It is at this point that we turn our attention from initial value problems to boundary value problems. These will occupy us for the remainder of the chapter. CITED REFERENCES AND FURTHER READING: Ames, W.F. 1977, Numerical Methods for Partial Differential Equations, 2nd ed. (New York: Academic Press). 19.4 Fourier and Cyclic Reduction Methods for Boundary Value Problems As discussed in §19.0, most boundary value problems (elliptic equations, for example) reduce to solving large sparse linear systems of the form A · u = b (19.4.1) either once, for boundary value equations that are linear, or iteratively, for boundary value equations that are nonlinear
858 Chapter 19.Partial Differential Equations Two important techniques lead to "rapid"solution of equation (19.4.1)when the sparse matrix is of certain frequently occurring forms.The Fourier transform method is directly applicable when the equations have coefficients that are constant in space.The cyclic reduction method is somewhat more general;its applicability is related to the question of whether the equations are separable(in the sense of "separation of variables").Both methods require the boundaries to coincide with the coordinate lines.Finally,for some problems,there is a powerful combination of these two methods called FACR (Fourier Analysis and Cyclic Reduction).We now consider each method in turn,using equation (19.0.3),with finite-difference representation(19.0.6),as a model example.Generally speaking,the methods in this section are faster,when they apply,than the simpler relaxation methods discussed in $19.5;but they are not necessarily faster than the more complicated multigrid methods discussed in 819.6. Fourier Transform Method RECIPES I The discrete inverse Fourier transform in both z and y is 1 J-1L-1 ujl Umne-2xijm/Je-2xiln/L (19.4.2) 。 Press. n=0n=0 This can be computed using the FFT independently in each dimension,or else all at once via the routine fourn of $12.4 or the routine rlft3 of $12.5.Similarly, IENTIFIC J-1L-1 Pmne-2xijm/Je-2xiln/L 19.4.3) JL n=0n=0 If we substitute expressions (19.4.2)and (19.4.3)in our model problem (19.0.6), we find Numerica 10621 立nn(e2mim/W+e-2mim/J+e2min/L+e-2min/L-4）=mn△2 (19.4.4) 431 or (outside Recipes pmn△2 (19.4.5) 2πm +cos 20一2 L North Thus the strategy for solving equation(19.0.6)by FFT techniques is: .Compute Pmn as the Fourier transform J-1L-1 Pmn= pil e2ximj/Je2xinl/L (19.4.6) j=0=0 .Compute timn from equation (19.4.5)
858 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). Two important techniques lead to “rapid” solution of equation (19.4.1) when the sparse matrix is of certain frequently occurring forms. The Fourier transform method is directly applicable when the equations have coefficients that are constant in space. The cyclic reduction method is somewhat more general; its applicability is related to the question of whether the equations are separable (in the sense of “separation of variables”). Both methods require the boundaries to coincide with the coordinate lines. Finally, for some problems, there is a powerful combination of these two methods called FACR (Fourier Analysis and Cyclic Reduction). We now consider each method in turn, using equation (19.0.3), with finite-difference representation (19.0.6), as a model example. Generally speaking, the methods in this section are faster, when they apply, than the simpler relaxation methods discussed in §19.5; but they are not necessarily faster than the more complicated multigrid methods discussed in §19.6. Fourier Transform Method The discrete inverse Fourier transform in both x and y is ujl = 1 JL J −1 m=0 L −1 n=0 umne−2πijm/J e−2πiln/L (19.4.2) This can be computed using the FFT independently in each dimension, or else all at once via the routine fourn of §12.4 or the routine rlft3 of §12.5. Similarly, ρjl = 1 JL J −1 m=0 L −1 n=0 ρmne−2πijm/J e−2πiln/L (19.4.3) If we substitute expressions (19.4.2) and (19.4.3) in our model problem (19.0.6), we find umn e2πim/J + e−2πim/J + e2πin/L + e−2πin/L − 4 = ρmn∆2 (19.4.4) or umn = ρmn∆2 2 cos 2πm J + cos 2πn L − 2 (19.4.5) Thus the strategy for solving equation (19.0.6) by FFT techniques is: • Compute ρmn as the Fourier transform ρmn = J −1 j=0 L −1 l=0 ρjl e2πimj/J e2πinl/L (19.4.6) • Compute umn from equation (19.4.5).
19.4 Fourier and Cyclic Reduction Methods 859 .Compute ui by the inverse Fourier transform (19.4.2). The above procedure is valid for periodic boundary conditions.In other words, the solution satisfies ujl uj+Jl=ujl+L (19.4.7) Next consider a Dirichlet boundary condition u=0 on the rectangular boundary. Instead of the expansion(19.4.2),we now need an expansion in sine waves: 22 J-1L-1 uj1了工 ∑立mn sin m sin tln si元 (19.4.8) m=1n=1 This satisfies the boundary conditions that u=0 atj=0,J and at I=0,L.If we g ICAL substitute this expansion and the analogous one for pj into equation(19.0.6),we find that the solution procedure parallels that for periodic boundary conditions: Compute pmn by the sine transform RECIPES 1-1L-1 9 pmn=】 sin rln (19.4.9) J j=1l=1 (A fast sine transform algorithm was given in 812.3.) Compute amn from the expression analogous to (19.4.5), e多口 umn= △2pmn /πm xn (19.4.10) 2(cos+cos元-2 Compute ujt by the inverse sine transform(19.4.8). If we have inhomogeneous boundary conditions,for example u =0 on all boundaries except u=f(y)on the boundary =JA,we have to add to the above solution a solution u of the homogeneous equation Numerica 10621 ￥uu=0 0r2+ 0y2 (19.4.11) that satisfies the required boundary conditions.In the continuum case,this would be an expression of the form ”=∑Asnh2贤血2器 (19.4.12) where An would be found by requiring that u=f(y)at x=JA.In the discrete case,we have -1 (19.4.13)
19.4 Fourier and Cyclic Reduction Methods 859 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). • Compute ujl by the inverse Fourier transform (19.4.2). The above procedure is valid for periodic boundary conditions. In other words, the solution satisfies ujl = uj+J,l = uj,l+L (19.4.7) Next consider a Dirichlet boundary condition u = 0 on the rectangular boundary. Instead of the expansion (19.4.2), we now need an expansion in sine waves: ujl = 2 J 2 L J −1 m=1 L −1 n=1 umn sin πjm J sin πln L (19.4.8) This satisfies the boundary conditions that u = 0 at j = 0, J and at l = 0, L. If we substitute this expansion and the analogous one for ρjl into equation (19.0.6), we find that the solution procedure parallels that for periodic boundary conditions: • Compute ρmn by the sine transform ρmn = J −1 j=1 L −1 l=1 ρjl sin πjm J sin πln L (19.4.9) (A fast sine transform algorithm was given in §12.3.) • Compute umn from the expression analogous to (19.4.5), umn = ∆2ρmn 2 cos πm J + cos πn L − 2 (19.4.10) • Compute ujl by the inverse sine transform (19.4.8). If we have inhomogeneous boundary conditions, for example u = 0 on all boundaries except u = f(y) on the boundary x = J∆, we have to add to the above solution a solution uH of the homogeneous equation ∂2u ∂x2 + ∂2u ∂y2 =0 (19.4.11) that satisfies the required boundary conditions. In the continuum case, this would be an expression of the form uH = n An sinh nπx L∆ sin nπy L∆ (19.4.12) where An would be found by requiring that u = f(y) at x = J∆. In the discrete case, we have uH jl = 2 L L −1 n=1 An sinh πnj L sin πnl L (19.4.13)
860 Chapter 19.Partial Differential Equations If f(y =IA)=fi,then we get An from the inverse formula 1 L-1 πnl An= sinh (nJ/L) (19.4.14) l=1 The complete solution to the problem is u ujl uf (19.4.15) By adding appropriate terms of the form(19.4.12),we can handle inhomogeneous 81 terms on any boundary surface. A much simpler procedure for handling inhomogeneous terms is to note that whenever boundary terms appear on the left-hand side of(19.0.6),they can be taken over to the right-hand side since they are known.The effective source term is ICAL therefore pit plus a contribution from the boundary terms.To implement this idea formally,write the solution as u=u'+uB (19.4.16) 9 where u'=0 on the boundary,while uB vanishes everywhere except on the boundary.There it takes on the given boundary value.In the above example,the only nonzero values of u would be u明=f无 (19.4.17) The model equation (19.0.3)becomes V2=-72uB+p (19.4.18) or,in finite-difference form, 4+1,1+-1l++1+叫-1-4.1= -(明1l+u月l+u+1+u唱-1-4u)+△2p (19.4.19) 容室 Nume 10621 All the uB terms in equation(19.4.19)vanish except when the equation is evaluated at j=J-1,where +-21+-141+4-1,1-1-4d-11=-f+△2p1-1l(（19.4.20) Thus the problem is now equivalent to the case of zero boundary conditions,except that one row of the source term is modified by the replacement △2pJ-1,l→△2pJ-1,l-fh (19.4.21) The case of Neumann boundary conditions Vu =0 is handled by the cosine expansion (12.3.17): L erjm。πln (19.4.22) m=0n=0
860 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). If f(y = l∆) ≡ fl, then we get An from the inverse formula An = 1 sinh (πnJ/L) L −1 l=1 fl sin πnl L (19.4.14) The complete solution to the problem is u = ujl + uH jl (19.4.15) By adding appropriate terms of the form (19.4.12), we can handle inhomogeneous terms on any boundary surface. A much simpler procedure for handling inhomogeneous terms is to note that whenever boundary terms appear on the left-hand side of (19.0.6), they can be taken over to the right-hand side since they are known. The effective source term is therefore ρjl plus a contribution from the boundary terms. To implement this idea formally, write the solution as u = u + uB (19.4.16) where u = 0 on the boundary, while uB vanishes everywhere except on the boundary. There it takes on the given boundary value. In the above example, the only nonzero values of uB would be uB J,l = fl (19.4.17) The model equation (19.0.3) becomes ∇2u = −∇2uB + ρ (19.4.18) or, in finite-difference form, u j+1,l + u j−1,l + u j,l+1 + u j,l−1 − 4u j,l = − (uB j+1,l + uB j−1,l + uB j,l+1 + uB j,l−1 − 4uB j,l)+∆2ρj,l (19.4.19) All the uB terms in equation (19.4.19) vanish except when the equation is evaluated at j = J − 1, where u J,l + u J−2,l + u J−1,l+1 + u J−1,l−1 − 4u J−1,l = −fl + ∆2ρJ−1,l (19.4.20) Thus the problem is now equivalent to the case of zero boundary conditions, except that one row of the source term is modified by the replacement ∆2ρJ−1,l → ∆2ρJ−1,l − fl (19.4.21) The case of Neumann boundary conditions ∇u = 0 is handled by the cosine expansion (12.3.17): ujl = 2 J 2 L J m=0 L n=0 umn cos πjm J cos πln L (19.4.22)
19.4 Fourier and Cyclic Reduction Methods 861 Here the double prime notation means that the terms for m=0 and m =should be multiplied by and similarly for n=0 and n=L.Inhomogeneous terms Vu=g can be again included by adding a suitable solution of the homogeneous equation,or more simply by taking boundary terms over to the right-hand side.For example,the condition Ou =9(y at x=0 (19.4.23) OT becomes 4-u-14=g0 (19.4.24) 2△ where g=g(y =IA).Once again we write the solution in the form (19.4.16), where now Vu'=0 on the boundary.This time VuB takes on the prescribed value on the boundary,but uB vanishes everywhere except just outside the boundary. Thus equation (19.4.24)gives RECIPES uB11=-2△g (19.4.25) 9 All the uB terms in equation (19.4.19)vanish except whenj=0: 三学 .1+∠11+6.+1+6.l-1-4u6.1=2△g+△2p0.l (19.4.26 Thus u'is the solution of a zero-gradient problem,with the source term modified by the replacement IENTIFIC △2p0.1→△2p0.1+2△9m (19.4.27) 6 Sometimes Neumann boundary conditions are handled by using a staggered grid,with the u's defined midway between zone boundaries so that first derivatives are centered on the mesh points.You can solve such problems using similar 9 techniques to those described above if you use the alternative form of the cosine transform,equation (12.3.23). Numerical 10521 uction 43106 Cyclic Reduction (outside Recipes Evidently the FFT method works only when the original PDE has constant coefficients.and boundaries that coincide with the coordinate lines.An alternative North Software. algorithm,which can be used on somewhat more general equations,is called cyclic reduction (CR). We illustrate cyclic reduction on the equation 82u,02u 0m2+2 +b(y) -+c()u=g(x,y) (19.4.28) This form arises very often in practice from the Helmholtz or Poisson equations in polar,cylindrical,or spherical coordinate systems.More general separable equations are treated in [1]
19.4 Fourier and Cyclic Reduction Methods 861 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). Here the double prime notation means that the terms for m = 0 and m = J should be multiplied by 1 2 , and similarly for n = 0 and n = L. Inhomogeneous terms ∇u = g can be again included by adding a suitable solution of the homogeneous equation, or more simply by taking boundary terms over to the right-hand side. For example, the condition ∂u ∂x = g(y) at x =0 (19.4.23) becomes u1,l − u−1,l 2∆ = gl (19.4.24) where gl ≡ g(y = l∆). Once again we write the solution in the form (19.4.16), where now ∇u = 0 on the boundary. This time ∇uB takes on the prescribed value on the boundary, but uB vanishes everywhere except just outside the boundary. Thus equation (19.4.24) gives uB −1,l = −2∆gl (19.4.25) All the uB terms in equation (19.4.19) vanish except when j = 0: u 1,l + u −1,l + u 0,l+1 + u 0,l−1 − 4u 0,l = 2∆gl + ∆2ρ0,l (19.4.26) Thus u is the solution of a zero-gradient problem, with the source term modified by the replacement ∆2ρ0,l → ∆2ρ0,l + 2∆gl (19.4.27) Sometimes Neumann boundary conditions are handled by using a staggered grid, with the u’s defined midway between zone boundaries so that first derivatives are centered on the mesh points. You can solve such problems using similar techniques to those described above if you use the alternative form of the cosine transform, equation (12.3.23). Cyclic Reduction Evidently the FFT method works only when the original PDE has constant coefficients, and boundaries that coincide with the coordinate lines. An alternative algorithm, which can be used on somewhat more general equations, is called cyclic reduction (CR). We illustrate cyclic reduction on the equation ∂2u ∂x2 + ∂2u ∂y2 + b(y) ∂u ∂y + c(y)u = g(x, y) (19.4.28) This form arises very often in practice from the Helmholtz or Poisson equations in polar, cylindrical, or spherical coordinate systems. More general separable equations are treated in [1].