c19-4 (779619), страница 2
Текст из файла (страница 2)
Inhomogeneous terms∇u = g can be again included by adding a suitable solution of the homogeneousequation, or more simply by taking boundary terms over to the right-hand side.For example, the conditionbecomesat x = 0u1,l − u−1,l= gl2∆(19.4.23)(19.4.24)where gl ≡ g(y = l∆). Once again we write the solution in the form (19.4.16),where now ∇u0 = 0 on the boundary. This time ∇uB takes on the prescribed valueon the boundary, but uB vanishes everywhere except just outside the boundary.Thus equation (19.4.24) givesuB−1,l = −2∆gl(19.4.25)All the uB terms in equation (19.4.19) vanish except when j = 0:u01,l + u0−1,l + u00,l+1 + u00,l−1 − 4u00,l = 2∆gl + ∆2 ρ0,l(19.4.26)Thus u0 is the solution of a zero-gradient problem, with the source term modifiedby the replacement∆2 ρ0,l → ∆2 ρ0,l + 2∆gl(19.4.27)Sometimes Neumann boundary conditions are handled by using a staggeredgrid, with the u’s defined midway between zone boundaries so that first derivativesare centered on the mesh points.
You can solve such problems using similartechniques to those described above if you use the alternative form of the cosinetransform, equation (12.3.23).Cyclic ReductionEvidently the FFT method works only when the original PDE has constantcoefficients, and boundaries that coincide with the coordinate lines. An alternativealgorithm, which can be used on somewhat more general equations, is called cyclicreduction (CR).We illustrate cyclic reduction on the equation∂u∂2u ∂2u+ c(y)u = g(x, y)+ 2 + b(y)2∂x∂y∂y(19.4.28)This form arises very often in practice from the Helmholtz or Poisson equations inpolar, cylindrical, or spherical coordinate systems. More general separable equationsare treated in [1].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).∂u= g(y)∂x862Chapter 19.Partial Differential EquationsThe finite-difference form of equation (19.4.28) can be written as a set ofvector equationsuj−1 + T · uj + uj+1 = gj ∆2(19.4.29)T = B − 21(19.4.30)where the 21 comes from the x-differencing and the matrix B from the y-differencing.The matrix B, and hence T, is tridiagonal with variable coefficients.The CR method is derived by writing down three successive equations like(19.4.29):uj−2 + T · uj−1 + uj = gj−1 ∆2uj−1 + T · uj + uj+1 = gj ∆2uj + T · uj+1 + uj+2 = gj+1 ∆(19.4.31)2Matrix-multiplying the middle equation by −T and then adding the three equations,we get(1)uj−2 + T(1) · uj + uj+2 = gj ∆2(19.4.32)This is an equation of the same form as (19.4.29), withT(1) = 21 − T2(1)gj= ∆2 (gj−1 − T · gj + gj+1 )(19.4.33)After one level of CR, we have reduced the number of equations by a factor oftwo.
Since the resulting equations are of the same form as the original equation, wecan repeat the process. Taking the number of mesh points to be a power of 2 forsimplicity, we finally end up with a single equation for the central line of variables:(f)T(f) · uJ/2 = ∆2 gJ/2 − u0 − uJ(19.4.34)Here we have moved u0 and uJ to the right-hand side because they are knownboundary values.
Equation (19.4.34) can be solved for uJ/2 by the standardtridiagonal algorithm. The two equations at level f − 1 involve uJ/4 and u3J/4 . Theequation for uJ/4 involves u0 and uJ/2 , both of which are known, and hence can besolved by the usual tridiagonal routine. A similar result holds true at every stage,so we end up solving J − 1 tridiagonal systems.In practice, equations (19.4.33) should be rewritten to avoid numerical instability. For these and other practical details, refer to [2].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).Here the index j comes from differencing in the x-direction, while the y-differencing(denoted by the index l previously) has been left in vector form. The matrix Thas the form19.5 Relaxation Methods for Boundary Value Problems863FACR Method(r)(r)kbkj + ubkj+2r = ∆2 gjubkj−2r + λk u(19.4.35)(r)Here λk is the eigenvalue of T(r) corresponding to the kth Fourier mode.
For(r)the equation (19.0.3), equation (19.4.5) shows that λk will involve terms likecos(2πk/L) − 2 raised to a power. Solve the tridiagonal systems for ubkj at the levelsrrrrj = 2 , 2 × 2 , 4 × 2 , ..., J − 2 . Fourier synthesize to get the y-values on thesex-lines. Then fill in the intermediate x-lines as in the original CR algorithm.The trick is to choose the number of levels of CR so as to minimize the totalnumber of arithmetic operations.
One can show that for a typical case of a 128×128mesh, the optimal level is r = 2; asymptotically, r → log2 (log2 J).A rough estimate of running times for these algorithms for equation (19.0.3)is as follows: The FFT method (in both x and y) and the CR method are roughlycomparable. FACR with r = 0 (that is, FFT in one dimension and solve thetridiagonal equations by the usual algorithm in the other dimension) gives about afactor of two gain in speed. The optimal FACR with r = 2 gives another factorof two gain in speed.CITED REFERENCES AND FURTHER READING:Swartzrauber, P.N. 1977, SIAM Review, vol.
19, pp. 490–501. [1]Buzbee, B.L, Golub, G.H., and Nielson, C.W. 1970, SIAM Journal on Numerical Analysis, vol. 7,pp. 627–656; see also op. cit. vol. 11, pp. 753–763. [2]Hockney, R.W. 1965, Journal of the Association for Computing Machinery, vol. 12, pp. 95–113.
[3]Hockney, R.W. 1970, in Methods of Computational Physics, vol. 9 (New York: Academic Press),pp. 135–211. [4]Hockney, R.W., and Eastwood, J.W. 1981, Computer Simulation Using Particles (New York:McGraw-Hill), Chapter 6. [5]Temperton, C. 1980, Journal of Computational Physics, vol. 34, pp. 314–329. [6]19.5 Relaxation Methods for Boundary ValueProblemsAs we mentioned in §19.0, relaxation methods involve splitting the sparsematrix that arises from finite differencing and then iterating until a solution is found.There is another way of thinking about relaxation methods that is somewhatmore physical. Suppose we wish to solve the elliptic equationLu = ρ(19.5.1)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 best way to solve equations of the form (19.4.28), including the constantcoefficient problem (19.0.3), is a combination of Fourier analysis and cyclic reduction,the FACR method [3-6].
If at the rth stage of CR we Fourier analyze the equations ofthe form (19.4.32) along y, that is, with respect to the suppressed vector index, wewill have a tridiagonal system in the x-direction for each y-Fourier mode:.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
