c17-3 (779603), страница 2
Текст из файла (страница 2)
We illustratethe steps schematically by figures.Consider the block of equations describing corrections available from the initial boundaryconditions. We have n1 equations for N unknown corrections. We wish to transform the firstblock so that its left-hand n1 × n1 square section becomes unity along the diagonal, and zeroin off-diagonal elements. Figure 17.3.3 shows the original and final form of the first blockof the matrix. In the figure we designate matrix elements that are subject to diagonalizationby “D”, and elements that will be altered by “A”; in the final block, elements that are storedare labeled by “S”. We get from start to finish by selecting in turn n1 “pivot” elements fromamong the first n1 columns, normalizing the pivot row so that the value of the “pivot” elementis unity, and adding appropriate multiples of this row to the remaining rows so that theycontain zeros in the pivot column.
In its final form, the reduced block expresses values for thecorrections to the first n1 variables at mesh point 1 in terms of values for the remaining n2unknown corrections at point 1, i.e., we now know what the first n1 elements are in terms ofthe remaining n2 elements. We store only the final set of n2 nonzero columns from the initialblock, plus the column for the altered right-hand side of the matrix equation.We must emphasize here an important detail of the method. To exploit the reducedstorage allowed by operating on blocks, it is essential that the ordering of columns in thes matrix of derivatives be such that pivot elements can be found among the first n1 rowsof the matrix. This means that the n1 boundary conditions at the first point must containsome dependence on the first j=1,2,...,n1 dependent variables, y[j][1].
If not, then theoriginal square n1 × n1 subsection of the first block will appear to be singular, and the methodwill fail. Alternatively, we would have to allow the search for pivot elements to involve allN columns of the block, and this would require column swapping and far more bookkeeping.The code provides a simple method of reordering the variables, i.e., the columns of the smatrix, so that this can be done easily. End of important detail.76517.3 Relaxation MethodsXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXVVVVVVVVVVVVVVVVVVVVBBBBBBBBBBBBBBBBBBBBFigure 17.3.1. Matrix structure of a set of linear finite-difference equations (FDEs) with boundaryconditions imposed at both endpoints. Here X represents a coefficient of the FDEs, V represents acomponent of the unknown solution vector, and B is a component of the known right-hand side.
Emptyspaces represent zeros. The matrix equation is to be solved by a special form of Gaussian elimination.(See text for details.)11X XX X1 X X11XXX1X11 X1XXXXX1XXX1X11 X1XXXXX1XXX1X11 X1XXXXX1VVVVVVVVVVVVVVVVVVVVBBBBBBBBBBBBBBBBBBBBFigure 17.3.2. Target structure of the Gaussian elimination. Once the matrix of Figure 17.3.1 has beenreduced to this form, the solution follows quickly by backsubstitution.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).XXXXXXXX766Chapter 17.Two Point Boundary Value ProblemsVVVAAA(b ) 1 0 0 S S0 1 0 S S0 0 1 S SVVVSSSFigure 17.3.3. Reduction process for the first (upper left) block of the matrix in Figure 17.3.1. (a)Original form of the block, (b) final form. (See text for explanation.)(a)(b )100ZZZZZ010ZZZZZ001ZZZZZSSSDDDDDSSSDDDDD100000000100000000100000SSS10000SSS01000DDDDD00100DDDDD00010DDDDD00001AAAAASSSSSAAAAAVVVVVVVVSSSAAAAASSSSSVVVVVVVVSSSSSSSSFigure 17.3.4. Reduction process for intermediate blocks of the matrix in Figure 17.3.1.
(a) Originalform, (b) final form. (See text for explanation.)(a)000000000000000100000100000100ZZ00010ZZ00001ZZSSSSSDDSSSSSDDVVVVVVVSSSSSAA(b )0000000000000001000001000001000000010000000100SSSSS10SSSSS01VVVVVVVSSSSSSSFigure 17.3.5. Reduction process for the last (lower right) block of the matrix in Figure 17.3.1. (a)Original form, (b) final form. (See text for explanation.)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).(a) D D D A AD D D A AD D D A A17.3 Relaxation Methods767err =nem XX1|∆Y [j][k]|m × nescalv[j]k=1 j=1(17.3.13)When err≤conv, the method has converged. Note that the user gets to supply an arrayscalv 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 ideato believe that the corrections generated from a first-order Taylor series are accurate.
Thenumber slowc modulates application of corrections. After each iteration we apply only aSample 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).Next consider the block of N equations representing the FDEs that describe the relationbetween the 2N corrections at points 2 and 1. The elements of that block, together with resultsfrom the previous step, are illustrated in Figure 17.3.4. Note that by adding suitable multiplesof rows from the first block we can reduce to zero the first n1 columns of the block (labeledby “Z”), and, to do so, we will need to alter only the columns from n1 + 1 to N and thevector element on the right-hand side. Of the remaining columns we can diagonalize a squaresubsection of N × N elements, labeled by “D” in the figure. In the process we alter the finalset of n2 + 1 columns, denoted “A” in the figure.
The second half of the figure shows theblock 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 correctionsat points 3 and 2, we see that the state of available results and new equations exactly reproducesthe situation described in the previous paragraph. Thus, we can carry out those steps againfor each block in turn through block M . Finally on block M + 1 we encounter the remainingboundary conditions.Figure 17.3.5 shows the final block of n2 FDEs relating the N corrections for variablesat mesh point M , together with the result of reducing the previous block.
Again, we can firstuse the prior results to zero the first n1 columns of the block. Now, when we diagonalizethe remaining square section, we strike gold: We get values for the final n2 correctionsat mesh point M .With the final block reduced, the matrix has the desired form shown previously inFigure 17.3.2, and the matrix is ripe for backsubstitution. Starting with the bottom row andworking up towards the top, at each stage we can simply determine one unknown correctionin terms of known quantities.The function solvde organizes the steps described above. The principal proceduresused in the algorithm are performed by functions called internally by solvde.
The functionred eliminates leading columns of the s matrix using results from prior blocks. pinvsdiagonalizes the square subsection of s and stores unreduced coefficients. bksub carriesout 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 smatrix for each block.Most of the arguments in the call to solvde have already been described, but somerequire discussion.
Array y[j][k] contains the initial guess for the solution, with j labelingthe dependent variables at mesh points k. The problem involves ne FDEs spanning pointsk=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 nbboundary conditions at k=1 involve the dependent variables referenced by the first nb columnsof the s matrix. Thus, columns j of the s matrix can be ordered by the user in difeq to referto derivatives with respect to the dependent variable indexv[j].The function only attempts itmax correction cycles before returning, even if the solutionhas not converged.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
