c19-4 (779619)
Текст из файла
19.4 Fourier and Cyclic Reduction Methods857For example, a combined advective-diffusion equation, such as∂u∂2 u∂u= −v+D 2∂t∂x∂x(19.3.21)un+1/m = U1 (un , ∆t/m)un+2/m = U2 (un+1/m , ∆t/m)(19.3.22)···un+1 = Um (un+(m−1)/m , ∆t/m)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 schemefor the operator L. In fact, as a rule of thumb, it is often sufficient to have stable Ui ’sonly for the operator pieces having the highest number of spatial derivatives — theother 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 toboundary 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 forBoundary Value ProblemsAs discussed in §19.0, most boundary value problems (elliptic equations, forexample) reduce to solving large sparse linear systems of the formA·u=b(19.4.1)either once, for boundary value equations that are linear, or iteratively, for boundaryvalue equations that are nonlinear.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).might profitably use an explicit scheme for the advective term combined with aCrank-Nicholson or other implicit scheme for the diffusion term.The alternating-direction implicit (ADI) method, equation (19.3.16), is anexample 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 thatincludes algebraically all the pieces of the total operator L, but which is desirablystable only for the L1 piece; likewise U2 , .
. . Um . Then a method of getting fromun to un+1 is858Chapter 19.Partial Differential EquationsFourier Transform MethodThe discrete inverse Fourier transform in both x and y isujl =J−1 L−11 XXubmn e−2πijm/J e−2πiln/LJL m=0 n=0(19.4.2)This can be computed using the FFT independently in each dimension, or else all atonce via the routine fourn of §12.4 or the routine rlft3 of §12.5.
Similarly,ρjl =J−1 L−11 XXρbmn e−2πijm/J e−2πiln/LJL m=0 n=0(19.4.3)If we substitute expressions (19.4.2) and (19.4.3) in our model problem (19.0.6),we findubmn e2πim/J + e−2πim/J + e2πin/L + e−2πin/L − 4 = ρbmn ∆2orubmn =ρbmn ∆22πn2πm+ cos−22 cosJL(19.4.4)(19.4.5)Thus the strategy for solving equation (19.0.6) by FFT techniques is:• Compute ρbmn as the Fourier transformρbmn =J−1X L−1Xρjl e2πimj/J e2πinl/Lj=0 l=0• Compute ubmn from equation (19.4.5).(19.4.6)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Two important techniques lead to “rapid” solution of equation (19.4.1) whenthe sparse matrix is of certain frequently occurring forms. The Fourier transformmethod is directly applicable when the equations have coefficients that are constantin space. The cyclic reduction method is somewhat more general; its applicabilityis related to the question of whether the equations are separable (in the sense of“separation of variables”).
Both methods require the boundaries to coincide withthe coordinate lines. Finally, for some problems, there is a powerful combinationof these two methods called FACR (Fourier Analysis and Cyclic Reduction). Wenow consider each method in turn, using equation (19.0.3), with finite-differencerepresentation (19.0.6), as a model example.
Generally speaking, the methods in thissection are faster, when they apply, than the simpler relaxation methods discussedin §19.5; but they are not necessarily faster than the more complicated multigridmethods discussed in §19.6.19.4 Fourier and Cyclic Reduction Methods859• Compute ujl by the inverse Fourier transform (19.4.2).The above procedure is valid for periodic boundary conditions. In other words,the solution satisfiesujl = uj+J,l = uj,l+L(19.4.7)ujl =J−1 L−1πln2 2 XXπjmsinubmn sinJ L m=1 n=1JL(19.4.8)This satisfies the boundary conditions that u = 0 at j = 0, J and at l = 0, L.
If wesubstitute this expansion and the analogous one for ρjl into equation (19.0.6), wefind that the solution procedure parallels that for periodic boundary conditions:• Compute ρbmn by the sine transformρbmn =J−1X L−1Xρjl sinj=1 l=1πlnπjmsinJL(19.4.9)(A fast sine transform algorithm was given in §12.3.)• Compute ubmn from the expression analogous to (19.4.5),ubmn =∆2 ρbmnπnπm+ cos−22 cosJL(19.4.10)• Compute ujl by the inverse sine transform (19.4.8).If we have inhomogeneous boundary conditions, for example u = 0 on allboundaries except u = f(y) on the boundary x = J∆, we have to add to the abovesolution a solution uH of the homogeneous equation∂2u ∂2u+ 2 =0∂x2∂y(19.4.11)that satisfies the required boundary conditions.
In the continuum case, this wouldbe an expression of the formuH =XnAn sinhnπynπxsinJ∆L∆(19.4.12)where An would be found by requiring that u = f(y) at x = J∆. In the discretecase, we haveuHjl =L−1πnl2 XπnjsinAn sinhL n=1JL(19.4.13)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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:860Chapter 19.Partial Differential EquationsIf f(y = l∆) ≡ fl , then we get An from the inverse formulaAn =L−1X1πnlfl sinsinh πnL(19.4.14)l=1u = ujl + uHjl(19.4.15)By adding appropriate terms of the form (19.4.12), we can handle inhomogeneousterms on any boundary surface.A much simpler procedure for handling inhomogeneous terms is to note thatwhenever boundary terms appear on the left-hand side of (19.0.6), they can be takenover to the right-hand side since they are known.
The effective source term istherefore ρjl plus a contribution from the boundary terms. To implement this ideaformally, write the solution asu = u0 + uB(19.4.16)where u0 = 0 on the boundary, while uB vanishes everywhere except on theboundary. There it takes on the given boundary value. In the above example, theonly nonzero values of uB would beuBJ,l = fl(19.4.17)The model equation (19.0.3) becomes∇2 u0 = −∇2 uB + ρ(19.4.18)or, in finite-difference form,u0j+1,l + u0j−1,l + u0j,l+1 + u0j,l−1 − 4u0j,l =BBBB2− (uBj+1,l + uj−1,l + uj,l+1 + uj,l−1 − 4uj,l ) + ∆ ρj,l(19.4.19)All the uB terms in equation (19.4.19) vanish except when the equation is evaluatedat j = J − 1, whereu0J,l + u0J−2,l + u0J−1,l+1 + u0J−1,l−1 − 4u0J−1,l = −fl + ∆2 ρJ−1,l(19.4.20)Thus the problem is now equivalent to the case of zero boundary conditions, exceptthat 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 cosineexpansion (12.3.17):ujl =LJπln2 2 X00 X00πjmcosubmn cosJ L m=0 n=0JL(19.4.22)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The complete solution to the problem is19.4 Fourier and Cyclic Reduction Methods861Here the double prime notation means that the terms for m = 0 and m = J shouldbe multiplied by 12 , and similarly for n = 0 and n = L.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
