c2-4 (779462), страница 2
Текст из файла (страница 2)
The solution vector x overwritesb[1..n]. The other input arrays are not modified, and can be left in place for successive callswith different right-hand sides.{unsigned long i,k,l;int mm;float dum;mm=m1+m2+1;l=m1;for (k=1;k<=n;k++) {Forward substitution, unscrambling the permuted rowsi=indx[k];as we go.if (i != k) SWAP(b[k],b[i])if (l < n) l++;for (i=k+1;i<=l;i++) b[i] -= al[k][i-k]*b[k];}l=1;for (i=n;i>=1;i--) {Backsubstitution.dum=b[i];for (k=2;k<=l;k++) dum -= a[i][k]*b[k+i-1];b[i]=dum/a[i][1];if (l < mm) l++;}}The routines bandec and banbks are based on the Handbook routines bandet1 andbansol1 in [1].CITED REFERENCES AND FURTHER READING:Keller, H.B.
1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA:Blaisdell), p. 74.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),Example 5.4.3, p. 166.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:McGraw-Hill), §9.11.Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol.
II of Handbook for Automatic Computation (New York: Springer-Verlag), Chapter I/6. [1]Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §4.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).Some pivoting is possible within the storage limitations of bandec, and the aboveroutine does take advantage of the opportunity. In general, when TINY is returned as adiagonal element of U , then the original matrix (perhaps as modified by roundoff error)is in fact singular.
In this regard, bandec is somewhat more robust than tridag above,which can fail algorithmically even for nonsingular matrices; bandec is thus also useful (withm1 = m2 = 1) for some ill-behaved tridiagonal systems.Once the matrix A has been decomposed, any number of right-hand sides can be solved inturn by repeated calls to banbks, the backsubstitution routine whose analog in §2.3 is lubksb.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 + δbx.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
