c2-5 (779463)
Текст из файла
552.5 Iterative Improvement of a Solution to Linear EquationsAδxx+bδxδbA−1Figure 2.5.1. Iterative improvement of the solution to A · x = b. The first guess x + δx is multiplied byA to produce b + δb. The known vector b is subtracted, giving δb. The linear set with this right-handside is inverted, giving δx. This is subtracted from the first guess giving an improved solution x.2.5 Iterative Improvement of a Solution toLinear EquationsObviously it is not easy to obtain greater precision for the solution of a linearset than the precision of your computer’s floating-point word.
Unfortunately, forlarge sets of linear equations, it is not always easy to obtain precision equal to, oreven comparable to, the computer’s limit. In direct methods of solution, roundofferrors accumulate, and they are magnified to the extent that your matrix is closeto singular.
You can easily lose two or three significant figures for matrices which(you thought) were far from singular.If this happens to you, there is a neat trick to restore the full machine precision,called iterative improvement of the solution. The theory is very straightforward (seeFigure 2.5.1): Suppose that a vector x is the exact solution of the linear setA·x=b(2.5.1)You don’t, however, know x.
You only know some slightly wrong solution x + δx,where δx is the unknown error. When multiplied by the matrix A, your slightly wrongsolution gives a product slightly discrepant from the desired right-hand side b, namelyA · (x + δx) = b + δb(2.5.2)Subtracting (2.5.1) from (2.5.2) givesA · δx = δb(2.5.3)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).b + δbx56Chapter 2.Solution of Linear Algebraic EquationsBut (2.5.2) can also be solved, trivially, for δb. Substituting this into (2.5.3) givesA · δx = A · (x + δx) − b(2.5.4)#include "nrutil.h"void mprove(float **a, float **alud, int n, int indx[], float b[], float x[])Improves a solution vector x[1..n] of the linear set of equations A · X = B.
The matrixa[1..n][1..n], and the vectors b[1..n] and x[1..n] are input, as is the dimension n.Also input is alud[1..n][1..n], the LU decomposition of a as returned by ludcmp, andthe vector indx[1..n] also returned by that routine. On output, only x[1..n] is modified,to an improved set of values.{void lubksb(float **a, int n, int *indx, float b[]);int j,i;double sdp;float *r;r=vector(1,n);for (i=1;i<=n;i++) {Calculate the right-hand side, accumulatingsdp = -b[i];the residual in double precision.for (j=1;j<=n;j++) sdp += a[i][j]*x[j];r[i]=sdp;}lubksb(alud,n,indx,r);Solve for the error term,for (i=1;i<=n;i++) x[i] -= r[i];and subtract it from the old solution.free_vector(r,1,n);}You should note that the routine ludcmp in §2.3 destroys the input matrix as itLU decomposes it. Since iterative improvement requires both the original matrixand its LU decomposition, you will need to copy A before calling ludcmp.
Likewiselubksb destroys b in obtaining x, so make a copy of b also. If you don’t mindthis extra storage, iterative improvement is highly recommended: It is a processof order only N 2 operations (multiply vector by matrix, and backsubstitute — seediscussion following equation 2.3.7); it never hurts; and it can really give you yourmoney’s worth if it saves an otherwise ruined solution on which you have alreadyspent of order N 3 operations.You can call mprove several times in succession if you want. Unless you arestarting quite far from the true solution, one call is generally enough; but a secondcall to verify convergence can be reassuring.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).In this equation, the whole right-hand side is known, since x + δx is the wrongsolution that you want to improve. It is essential to calculate the right-hand sidein double precision, since there will be a lot of cancellation in the subtraction of b.Then, we need only solve (2.5.4) for the error δx, then subtract this from the wrongsolution to get an improved solution.An important extra benefit occurs if we obtained the original solution by LUdecomposition.
In this case we already have the LU decomposed form of A, and allwe need do to solve (2.5.4) is compute the right-hand side and backsubstitute!The code to do all this is concise and straightforward:572.5 Iterative Improvement of a Solution to Linear EquationsMore on Iterative ImprovementR ≡ 1 − B0 · A(2.5.5)which is supposed to be “small” (we will be more precise below).
Note that thereforeB0 · A = 1 − R(2.5.6)Next consider the following formal manipulation:−1−1· B−1· B0A−1 = A−1 · (B−10 · B0 ) = (A0 ) · B0 = (B0 · A)= (1 − R)−1 · B0 = (1 + R + R2 + R3 + · · ·) · B0(2.5.7)We can define the nth partial sum of the last expression byBn ≡ (1 + R + · · · + Rn ) · B0(2.5.8)−1so that B∞ → A , if the limit exists.It now is straightforward to verify that equation (2.5.8) satisfies some interestingrecurrence relations. As regards solving A · x = b, where x and b are vectors, definex n ≡ Bn · b(2.5.9)xn+1 = xn + B0 · (b − A · xn )(2.5.10)Then it is easy to show thatThis is immediately recognizable as equation (2.5.4), with −δx = xn+1 − xn , and with B0taking the role of A−1 . We see, therefore, that equation (2.5.4) does not require that the LUdecompositon of A be exact, but only that the implied residual R be small.
In rough terms, ifthe residual is smaller than the square root of your computer’s roundoff error, then after oneapplication of equation (2.5.10) (that is, going from x0 ≡ B0 · b to x1 ) the first neglected term,of order R2 , will be smaller than the roundoff error. Equation (2.5.10), like equation (2.5.4),moreover, can be applied more than once, since it uses only B0 , and not any of the higher B’s.A much more surprising recurrence which follows from equation (2.5.8) is one thatmore than doubles the order n at each stage:B2n+1 = 2Bn − Bn · A · Bnn = 0, 1, 3, 7, . .
.(2.5.11)Repeated application of equation (2.5.11), from a suitable starting matrix B0 , convergesquadratically to the unknown inverse matrix A−1 (see §9.4 for the definition of “quadratically”). Equation (2.5.11) goes by various names, including Schultz’s Method and Hotelling’sMethod; see Pan and Reif [1] for references. In fact, equation (2.5.11) is simply the iterativeNewton-Raphson method of root-finding (§9.4) applied to matrix inversion.Before you get too excited about equation (2.5.11), however, you should notice that itinvolves two full matrix multiplications at each iteration.
Each matrix multiplication involvesN 3 adds and multiplies. But we already saw in §§2.1–2.3 that direct inversion of A requiresonly N 3 adds and N 3 multiplies in toto. Equation (2.5.11) is therefore practical only whenspecial circumstances allow it to be evaluated much more rapidly than is the case for generalmatrices. We will meet such circumstances later, in §13.10.In the spirit of delayed gratification, let us nevertheless pursue the two related issues:When does the series in equation (2.5.7) converge; and what is a suitable initial guess B0 (if,for example, an initial LU decomposition is not feasible)?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).It is illuminating (and will be useful later in the book) to give a somewhat more solidanalytical foundation for equation (2.5.4), and also to give some additional results. Implicit inthe previous discussion was the notion that the solution vector x + δx has an error term; butwe neglected the fact that the LU decomposition of A is itself not exact.A different analytical approach starts with some matrix B0 that is assumed to be anapproximate inverse of the matrix A, so that B0 · A is approximately the identity matrix 1.Define the residual matrix R of B0 as58Chapter 2.Solution of Linear Algebraic EquationsWe can define the norm of a matrix as the largest amplification of length that it isable to induce on a vector,kRk ≡ maxv6=0|R · v||v|(2.5.12)kRk < 1(2.5.13)Pan and Reif [1] point out that a suitable initial guess for B0 is any sufficiently small constant times the matrix transpose of A, that is,B0 = ATR = 1 − AT · Aor(2.5.14)To see why this is so involves concepts from Chapter 11; we give here only the briefestsketch: AT · A is a symmetric, positive definite matrix, so it has real, positive eigenvalues.In its diagonal representation, R takes the formR = diag(1 − λ1 , 1 − λ2 , .
. . , 1 − λN )(2.5.15)where all the λi ’s are positive. Evidently any satisfying 0 < < 2/(maxi λi ) will givekRk < 1. It is not difficult to show that the optimal choice for , giving the most rapidconvergence for equation (2.5.11), is = 2/(max λi + min λi )i(2.5.16)iRarely does one know the eigenvalues of AT · A in equation (2.5.16). Pan and Reifderive several interesting bounds, which are computable directly from A. The followingchoices guarantee the convergence of Bn as n → ∞,XXX≤1a2jkor≤1max|aij | × max|aij |(2.5.17)j,kijjiThe latter expression is truly a remarkable formula, which Pan and Reif derive by noting thatthe vector norm in equation (2.5.12) need not be the usual L2 norm, but can instead be eitherthe L∞ (max) norm, or the L1 (absolute value) norm.
See their work for details.Another approach, with which we have had some success, is to estimate the largesteigenvalue statistically, by calculating si ≡ |A · vi |2 for several unit vector vi ’s with randomlychosen directions in N -space. The largest eigenvalue λ can then be bounded by the maximumof 2 max si and 2N Var(si )/µ(si ), where Var and µ denote the sample variance and mean,respectively.CITED REFERENCES AND FURTHER READING:Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), §2.3.4, p.
55.Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), p. 74.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§5.5.6, p. 183.Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Englewood Cliffs, NJ: Prentice-Hall), Chapter 13.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
