c19-5 (779620)
Текст из файла
19.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:864Chapter 19.Partial Differential Equationswhere L represents some elliptic operator and ρ is the source term. Rewrite theequation as a diffusion equation,∂u= Lu − ρ∂t(19.5.2)∂u∂2u ∂2u=+ 2 −ρ∂t∂x2∂y(19.5.3)If we use FTCS differencing (cf.
equation 19.2.4), we getnun+1j,l = uj,l +∆t nuj+1,l + unj−1,l + unj,l+1 + unj,l−1 − 4unj,l − ρj,l ∆t (19.5.4)2∆Recall from (19.2.6) that FTCS differencing is stable in one spatial dimension only if∆t/∆2 ≤ 12 . In two dimensions this becomes ∆t/∆2 ≤ 14 . Suppose we try to takethe largest possible timestep, and set ∆t = ∆2 /4. Then equation (19.5.4) becomesun+1j,l = ∆21 nuj+1,l + unj−1,l + unj,l+1 + unj,l−1 −ρj,l44(19.5.5)Thus the algorithm consists of using the average of u at its four nearest-neighborpoints on the grid (plus the contribution from the source). This procedure is theniterated until convergence.This method is in fact a classical method with origins dating back to thelast century, called Jacobi’s method (not to be confused with the Jacobi methodfor eigenvalues).
The method is not practical because it converges too slowly.However, it is the basis for understanding the modern methods, which are alwayscompared with it.Another classical method is the Gauss-Seidel method, which turns out to beimportant in multigrid methods (§19.6). Here we make use of updated values of u onthe right-hand side of (19.5.5) as soon as they become available. In other words, theaveraging is done “in place” instead of being “copied” from an earlier timestep to alater one. If we are proceeding along the rows, incrementing j for fixed l, we haveun+1j,l = ∆21 nnn+1uj+1,l + un+1ρj,l+u+uj,l+1j−1,lj,l−1 −44(19.5.6)This method is also slowly converging and only of theoretical interest when used byitself, but some analysis of it will be instructive.Let us look at the Jacobi and Gauss-Seidel methods in terms of the matrixsplitting concept.
We change notation and call u “x,” to conform to standard matrixnotation. To solveA·x=b(19.5.7)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).An initial distribution u relaxes to an equilibrium solution as t → ∞. Thisequilibrium has all time derivatives vanishing. Therefore it is the solution of theoriginal elliptic problem (19.5.1). We see that all the machinery of §19.2, on diffusiveinitial value equations, can be brought to bear on the solution of boundary valueproblems by relaxation methods.Let us apply this idea to our model problem (19.0.3). The diffusion equation is19.5 Relaxation Methods for Boundary Value Problems865we can consider splitting A asA =L+D+U(19.5.8)D · x(r) = −(L + U) · x(r−1) + b(19.5.9)For our model problem (19.5.5), D is simply the identity matrix.
The Jacobi methodconverges for matrices A that are “diagonally dominant” in a sense that can bemade mathematically precise. For matrices arising from finite differencing, thiscondition is usually met.What is the rate of convergence of the Jacobi method? A detailed analysis isbeyond our scope, but here is some of the flavor: The matrix −D−1 · (L + U) isthe iteration matrix which, apart from an additive term, maps one set of x’s into thenext. The iteration matrix has eigenvalues, each one of which reflects the factor bywhich the amplitude of a particular eigenmode of undesired residual is suppressedduring one iteration. Evidently those factors had better all have modulus < 1 forthe relaxation to work at all! The rate of convergence of the method is set by therate for the slowest-decaying eigenmode, i.e., the factor with largest modulus.
Themodulus of this largest factor, therefore lying between 0 and 1, is called the spectralradius of the relaxation operator, denoted ρs .The number of iterations r required to reduce the overall error by a factor10−p is thus estimated byr≈p ln 10(− ln ρs )(19.5.10)In general, the spectral radius ρs goes asymptotically to the value 1 as the gridsize J is increased, so that more iterations are required. For any given equation,grid geometry, and boundary condition, the spectral radius can, in principle, becomputed analytically.
For example, for equation (19.5.5) on a J × J grid withDirichlet boundary conditions on all four sides, the asymptotic formula for largeJ turns out to beρs ' 1 −π22J 2(19.5.11)The number of iterations r required to reduce the error by a factor of 10−p is thusr'12pJ 2 ln 10' pJ 2π22(19.5.12)In other words, the number of iterations is proportional to the number of mesh points,J 2 . Since 100 × 100 and larger problems are common, it is clear that the Jacobimethod is only of theoretical interest.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).where D is the diagonal part of A, L is the lower triangle of A with zeros on thediagonal, and U is the upper triangle of A with zeros on the diagonal.In the Jacobi method we write for the rth step of iteration866Chapter 19.Partial Differential EquationsThe Gauss-Seidel method, equation (19.5.6), corresponds to the matrix decomposition(L + D) · x(r) = −U · x(r−1) + b(19.5.13)ρs ' 1 −r'π2J21pJ 2 ln 10' pJ 2π24(19.5.14)(19.5.15)The factor of two improvement in the number of iterations over the Jacobi methodstill leaves the method impractical.Successive Overrelaxation (SOR)We get a better algorithm — one that was the standard algorithm until the 1970s— if we make an overcorrection to the value of x(r) at the rth stage of Gauss-Seideliteration, thus anticipating future corrections.
Solve (19.5.13) for x(r) , add andsubtract x(r−1) on the right-hand side, and hence write the Gauss-Seidel method asx(r) = x(r−1) − (L + D)−1 · [(L + D + U) · x(r−1) − b](19.5.16)The term in square brackets is just the residual vector ξ (r−1) , sox(r) = x(r−1) − (L + D)−1 · ξ (r−1)(19.5.17)Now overcorrect, definingx(r) = x(r−1) − ω(L + D)−1 · ξ (r−1)(19.5.18)Here ω is called the overrelaxation parameter, and the method is called successiveoverrelaxation (SOR).The following theorems can be proved [1-3] :• The method is convergent only for 0 < ω < 2. If 0 < ω < 1, we speakof underrelaxation.• Under certain mathematical restrictions generally satisfied by matricesarising from finite differencing, only overrelaxation (1 < ω < 2 ) can givefaster convergence than the Gauss-Seidel method.• If ρJacobi is the spectral radius of the Jacobi iteration (so that the squareof it is the spectral radius of the Gauss-Seidel iteration), then the optimalchoice for ω is given byω=1+p21 − ρ2Jacobi(19.5.19)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 fact that L is on the left-hand side of the equation follows from the updatingin place, as you can easily check if you write out (19.5.13) in components.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
