868 Chapter 19.Partial Differential Equations Consider a general second-order elliptic equation in z and y,finite differenced on a square as for our model equation.Corresponding to each row of the matrix A is an equation of the form aj.luj+i,l+bj.luj-1,1+cj.luj+1+dj,luj.1-1+ej.luj.l fj.l (19.5.25) For our model equation,we had a =b=c=d =1,e =-4.The quantity f is proportional to the source term.The iterative procedure is defined by solving (19.5.25)for uj,. 三 *l= (G1-a西+11-b西-1u-69+1-d4-) (19.5.26) eil Then u is a weighted average ugg"wru"jt+(1-w)ug (19.5.27) We calculate it as follows:The residual at any stage is ⊙】 RECIPES 5,1=a,luj+1,1+b5.14j-1,1+C3.1u,+1+4.1u.l-1+ej1u.1-f.1(19.5.28) 9 and the SOR algorithm (19.5.18)or (19.5.27)is 5"=H-w赳 (19.5.29) ei.l This formulation is very easy to program,and the norm of the residual vector j. IENTIFIC can be used as a criterion for terminating the iteration. Another practical point concerns the order in which mesh points are processed. The obvious strategy is simply to proceed in order down the rows (or columns). Alternatively,suppose we divide the mesh into"odd"and"even"meshes,like the red and black squares of a checkerboard.Then equation(19.5.26)shows that the odd points depend only on the even mesh values and vice versa.Accordingly, we can carry out one half-sweep updating the odd points,say,and then another half-sweep updating the even points with the new odd values.For the version of Numerica 10.621 SOR implemented below,we shall adopt odd-even ordering. 431 The last practical point is that in practice the asymptotic rate of convergence in SOR is not attained until of order J iterations.The error often grows by a factor of 20 before convergence sets in.A trivial modification to SOR resolves this (outside Recipes problem.It is based on the observation that,while w is the optimum asymptotic relaxation parameter,it is not necessarily a good initial choice.In SOR with North Chebyshev acceleration,one uses odd-even ordering and changes w at each half- sweep according to the following prescription: w(0)=1 w1/2)=1/(1-Pjacobi/2） wn+1/2)=1/1-p7 acobi(n)/4), (19.5.30) n=1/2,1,,0∞ w(oo）一Woptimal868 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 machine￾readable 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). Consider a general second-order elliptic equation in x and y, finite differenced on a square as for our model equation. Corresponding to each row of the matrix A is an equation of the form aj,luj+1,l + bj,luj−1,l + cj,luj,l+1 + dj,luj,l−1 + ej,luj,l = fj,l (19.5.25) For our model equation, we had a = b = c = d = 1, e = −4. The quantity f is proportional to the source term. The iterative procedure is defined by solving (19.5.25) for uj,l: u*j,l = 1 ej,l (fj,l − aj,luj+1,l − bj,luj−1,l − cj,luj,l+1 − dj,luj,l−1) (19.5.26) Then unew j,l is a weighted average unew j,l = ωu*j,l + (1 − ω)uold j,l (19.5.27) We calculate it as follows: The residual at any stage is ξj,l = aj,luj+1,l + bj,luj−1,l + cj,luj,l+1 + dj,luj,l−1 + ej,luj,l − fj,l (19.5.28) and the SOR algorithm (19.5.18) or (19.5.27) is unew j,l = uold j,l − ω ξj,l ej,l (19.5.29) This formulation is very easy to program, and the norm of the residual vector ξ j,l can be used as a criterion for terminating the iteration. Another practical point concerns the order in which mesh points are processed. The obvious strategy is simply to proceed in order down the rows (or columns). Alternatively, suppose we divide the mesh into “odd” and “even” meshes, like the red and black squares of a checkerboard. Then equation (19.5.26) shows that the odd points depend only on the even mesh values and vice versa. Accordingly, we can carry out one half-sweep updating the odd points, say, and then another half-sweep updating the even points with the new odd values. For the version of SOR implemented below, we shall adopt odd-even ordering. The last practical point is that in practice the asymptotic rate of convergence in SOR is not attained until of order J iterations. The error often grows by a factor of 20 before convergence sets in. A trivial modification to SOR resolves this problem. It is based on the observation that, while ω is the optimum asymptotic relaxation parameter, it is not necessarily a good initial choice. In SOR with Chebyshev acceleration, one uses odd-even ordering and changes ω at each half￾sweep according to the following prescription: ω(0) = 1 ω(1/2) = 1/(1 − ρ2 Jacobi/2) ω(n+1/2) = 1/(1 − ρ2 Jacobiω(n) /4), n = 1/2, 1, ..., ∞ ω(∞) → ωoptimal (19.5.30)