862 Chapter 19.Partial Differential Equations The finite-difference form of equation(19.4.28)can be written as a set of vector equations -1+T.u+uj+1=g影△2 (19.4.29) Here the index j comes from differencing in the z-direction,while the y-differencing (denoted by the index I previously)has been left in vector form.The matrix T has the form T=B-21 (19.4.30) where the 21 comes from the r-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) ICAL uj-2+Tu-1+u=g)-1△2 u-1+Tu吗+j+1=g影△2 RECIPES (19.4.31) 马+T.“+1+马+2=87+1△2 9 Matrix-multiplying the middle equation by-T and then adding the three equations, 岂子 we get u-2+T0.马+u+2=g△2 (19.4.32) This is an equation of the same form as(19.4.29),with IENTIFIC 61 T1)=21-T2 (19.4.33) g=△2(g-1-T·g+g+1) After one level of CR.we have reduced the number of equations by a factor of two.Since the resulting equations are of the same form as the original equation,we 热。 10.621 can repeat the process.Taking the number of mesh points to be a power of 2 for Numerica simplicity,we finally end up with a single equation for the central line of variables: 43106 T.2=△2g界--u山 (19.4.34) Here we have moved uo and u to the right-hand side because they are known boundary values.Equation (19.4.34)can be solved for uJ/2 by the standard tridiagonal algorithm.The two equations at level f-1 involve uJ/4 and u3J/4.The equation for uJ/involves uo and uJ/2,both of which are known,and hence can be solved by the usual tridiagonal routine.A similar result holds true at every stage, so we end up solvingJ-1 tridiagonal systems. In practice,equations(19.4.33)should be rewritten to avoid numerical instabil- ity.For these and other practical details,refer to [2].862 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). The finite-difference form of equation (19.4.28) can be written as a set of vector equations uj−1 + T · uj + uj+1 = gj∆2 (19.4.29) 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 T has the form 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∆2 uj−1 + T · uj + uj+1 = gj∆2 uj + T · uj+1 + uj+2 = gj+1∆2 (19.4.31) Matrix-multiplying the middle equation by −T and then adding the three equations, we get uj−2 + T(1) · uj + uj+2 = g (1) j ∆2 (19.4.32) This is an equation of the same form as (19.4.29), with T(1) = 21 − T2 g (1) j = ∆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 of two. Since the resulting equations are of the same form as the original equation, we can repeat the process. Taking the number of mesh points to be a power of 2 for simplicity, we finally end up with a single equation for the central line of variables: T(f) · uJ/2 = ∆2g (f) J/2 − u0 − uJ (19.4.34) Here we have moved u0 and uJ to the right-hand side because they are known boundary values. Equation (19.4.34) can be solved for u J/2 by the standard tridiagonal algorithm. The two equations at level f − 1 involve u J/4 and u3J/4. The equation for uJ/4 involves u0 and uJ/2, both of which are known, and hence can be solved 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 instabil￾ity. For these and other practical details, refer to [2]