17.3 Relaxation Methods 767 Next consider the block of N equations representing the FDEs that describe the relation between the 2N corrections at points 2 and 1.The elements of that block,together with results from the previous step,are illustrated in Figure 17.3.4.Note that by adding suitable multiples of rows from the first block we can reduce to zero the first n columns of the block (labeled by "Z"),and,to do so,we will need to alter only the columns from n+1 to N and the vector element on the right-hand side.Of the remaining columns we can diagonalize a square subsection of N x N elements,labeled by "D"in the figure.In the process we alter the final set of n2 +1 columns,denoted "A"in the figure.The second half of the figure shows the block when we finish operating on it,with the stored (n2+1)x N elements labeled by "S." If we operate on the next set ofequations corresponding to the FDEs coupling corrections at points 3 and 2,we see that the state of available results and new equations exactly reproduces the situation described in the previous paragraph.Thus,we can carry out those steps again 8 for each block in turn through block M.Finally on block M+1 we encounter the remaining boundary conditions. Figure 17.3.5 shows the final block of n2 FDEs relating the N corrections for variables at mesh point M,together with the result of reducing the previous block.Again,we can first use the prior results to zero the first n columns of the block.Now,when we diagonalize the remaining square section,we strike gold:We get values for the final n2 corrections at mesh point M. With the final block reduced,the matrix has the desired form shown previously in Figure 17.3.2,and the matrix is ripe for backsubstitution.Starting with the bottom row and 令 working up towards the top,at each stage we can simply determine one unknown correction in terms of known quantities. The function solvde organizes the steps described above.The principal procedures used in the algorithm are performed by functions called internally by solvde.The function Press. ART red eliminates leading columns of the s matrix using results from prior blocks.pinvs diagonalizes the square subsection of s and stores unreduced coefficients.bksub carries out the backsubstitution step.The user of solvde must understand the calling arguments, 9 Programs as described below,and supply a function difeg,called by solvde,that evaluates the s matrix for each block. Most of the arguments in the call to solvde have already been described,but some require discussion.Array y[j][k]contains the initial guess for the solution,with j labeling the dependent variables at mesh points k.The problem involves ne FDEs spanning points k=1,...,m.nb boundary conditions apply at the first point k=1.The array indexv[j] establishes the correspondence between columns of the s matrix,equations(17.3.8),(17.3.10). SCIENTIFIC COMPUTING(ISBN and (17.3.12),and the dependent variables.As described above it is essential that the nb 198918920 boundary conditions at k=1 involve the dependent variables referenced by the first nb columns of the s matrix.Thus,columns j of the s matrix can be ordered by the user in difeg to refer to derivatives with respect to the dependent variable indexv [j]. Numerical Recipes 10-621 The function only attempts itmax correction cycles before returning,even if the solution has not converged.The parameters conv,slowc,scalv relate to convergence.Each -43108 inversion of the matrix produces corrections for ne variables at m mesh points.We want these to become vanishingly small as the iterations proceed,but we must define a measure for the size of corrections. This error"norm"is very problem specific,so the user might wish to (outside rewrite this section of the code as appropriate.In the program below we compute a value for the average correction err by summing the absolute value of all corrections,weighted by North Software. a scale factor appropriate to each type of variable 言 m ne 1 |△Y[][k] err= (17.3.13) m×ne scalv[j] k=1j=1 When err<conv,the method has converged.Note that the user gets to supply an array scalv which measures the typical size of each variable. Obviously,if err is large,we are far from a solution,and perhaps it is a bad idea to believe that the corrections generated from a first-order Taylor series are accurate.The number slowc modulates application of corrections.After each iteration we apply only a17.3 Relaxation Methods 767 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). Next consider the block of N equations representing the FDEs that describe the relation between the 2N corrections at points 2 and 1. The elements of that block, together with results from the previous step, are illustrated in Figure 17.3.4. Note that by adding suitable multiples of rows from the first block we can reduce to zero the first n1 columns of the block (labeled by “Z”), and, to do so, we will need to alter only the columns from n1 + 1 to N and the vector element on the right-hand side. Of the remaining columns we can diagonalize a square subsection of N × N elements, labeled by “D” in the figure. In the process we alter the final set of n2 + 1 columns, denoted “A” in the figure. The second half of the figure shows the block when we finish operating on it, with the stored (n2 + 1) × N elements labeled by “S.” If we operate on the next set of equations corresponding to the FDEs coupling corrections at points 3 and 2, we see that the state of available results and new equations exactly reproduces the situation described in the previous paragraph. Thus, we can carry out those steps again for each block in turn through block M. Finally on block M + 1 we encounter the remaining boundary conditions. Figure 17.3.5 shows the final block of n2 FDEs relating the N corrections for variables at mesh point M, together with the result of reducing the previous block. Again, we can first use the prior results to zero the first n1 columns of the block. Now, when we diagonalize the remaining square section, we strike gold: We get values for the final n2 corrections at mesh point M. With the final block reduced, the matrix has the desired form shown previously in Figure 17.3.2, and the matrix is ripe for backsubstitution. Starting with the bottom row and working up towards the top, at each stage we can simply determine one unknown correction in terms of known quantities. The function solvde organizes the steps described above. The principal procedures used in the algorithm are performed by functions called internally by solvde. The function red eliminates leading columns of the s matrix using results from prior blocks. pinvs diagonalizes the square subsection of s and stores unreduced coefficients. bksub carries out the backsubstitution step. The user of solvde must understand the calling arguments, as described below, and supply a function difeq, called by solvde, that evaluates the s matrix for each block. Most of the arguments in the call to solvde have already been described, but some require discussion. Array y[j][k] contains the initial guess for the solution, with j labeling the dependent variables at mesh points k. The problem involves ne FDEs spanning points k=1,..., m. nb boundary conditions apply at the first point k=1. The array indexv[j] establishes the correspondence between columns of the s matrix, equations (17.3.8), (17.3.10), and (17.3.12), and the dependent variables. As described above it is essential that the nb boundary conditions at k=1 involve the dependent variables referenced by the first nb columns of the s matrix. Thus, columns j of the s matrix can be ordered by the user in difeq to refer to derivatives with respect to the dependent variable indexv[j]. The function only attempts itmax correction cycles before returning, even if the solution has not converged. The parameters conv, slowc, scalv relate to convergence. Each inversion of the matrix produces corrections for ne variables at m mesh points. We want these to become vanishingly small as the iterations proceed, but we must define a measure for the size of corrections. This error “norm” is very problem specific, so the user might wish to rewrite this section of the code as appropriate. In the program below we compute a value for the average correction err by summing the absolute value of all corrections, weighted by a scale factor appropriate to each type of variable: err = 1 m × ne m k=1 ne j=1 |∆Y [j][k]| scalv[j] (17.3.13) When err≤conv, the method has converged. Note that the user gets to supply an array scalv which measures the typical size of each variable. Obviously, if err is large, we are far from a solution, and perhaps it is a bad idea to believe that the corrections generated from a first-order Taylor series are accurate. The number slowc modulates application of corrections. After each iteration we apply only a