c2-10 (779468)
Текст из файла
98Chapter 2.Solution of Linear Algebraic Equationsx[i]=sum/p[i];}}for (i=1;i<=n;i++) {a[i][i]=1.0/p[i];for (j=i+1;j<=n;j++) {sum=0.0;for (k=i;k<j;k++) sum -= a[j][k]*a[k][i];a[j][i]=sum/p[j];}}CITED REFERENCES AND FURTHER READING:Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag), Chapter I/1.Gill, P.E., Murray, W., and Wright, M.H. 1991, Numerical Linear Algebra and Optimization, vol. 1(Redwood City, CA: Addison-Wesley), §4.9.2.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§5.3.5.Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §4.2.2.10 QR DecompositionThere is another matrix factorization that is sometimes very useful, the so-called QRdecomposition,A = Q·R(2.10.1)Here R is upper triangular, while Q is orthogonal, that is,QT · Q = 1(2.10.2)where QT is the transpose matrix of Q.
Although the decomposition exists for a generalrectangular matrix, we shall restrict our treatment to the case when all the matrices are square,with dimensions N × N .Like the other matrix factorizations we have met (LU , SVD, Cholesky), QR decomposition can be used to solve systems of linear equations. To solveA·x =b(2.10.3)R · x = QT · b(2.10.4)first form Q · b and then solveTby backsubstitution. Since QR decomposition involves about twice as many operations asLU decomposition, it is not used for typical systems of linear equations.
However, we willmeet special cases where QR is the method of choice.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 typical use of choldc and cholsl is in the inversion of covariance matrices describingthe fit of data to a model; see, e.g., §15.6. In this, and many other applications, one often needsL−1 .
The lower triangle of this matrix can be efficiently found from the output of choldc:2.10 QR Decomposition99R = Qn−1 · · · Q1 · A(2.10.5)Since the Householder matrices are orthogonal,Q = (Qn−1 · · · Q1 )−1 = Q1 · · · Qn−1(2.10.6)In most applications we don’t need to form Q explicitly; we instead store it in the factoredform (2.10.6). Pivoting is not usually necessary unless the matrix A is very close to singular.A general QR algorithm for rectangular matrices including pivoting is given in [1]. For squarematrices, an implementation is the following:#include <math.h>#include "nrutil.h"void qrdcmp(float **a, int n, float *c, float *d, int *sing)Constructs the QR decomposition of a[1..n][1..n].
The upper triangular matrix R is returned in the upper triangle of a, except for the diagonal elements of R which are returned ind[1..n]. The orthogonal matrix Q is represented as a product of n − 1 Householder matricesQ1 . . . Qn−1 , where Qj = 1 − uj ⊗ uj /cj . The ith component of uj is zero for i = 1, . . . , j − 1while the nonzero components are returned in a[i][j] for i = j, . .
. , n. sing returns astrue (1) if singularity is encountered during the decomposition, but the decomposition is stillcompleted in this case; otherwise it returns false (0).{int i,j,k;float scale,sigma,sum,tau;*sing=0;for (k=1;k<n;k++) {scale=0.0;for (i=k;i<=n;i++) scale=FMAX(scale,fabs(a[i][k]));if (scale == 0.0) {Singular case.*sing=1;c[k]=d[k]=0.0;} else {Form Qk and Qk · A.for (i=k;i<=n;i++) a[i][k] /= scale;for (sum=0.0,i=k;i<=n;i++) sum += SQR(a[i][k]);sigma=SIGN(sqrt(sum),a[k][k]);a[k][k] += sigma;c[k]=sigma*a[k][k];d[k] = -scale*sigma;for (j=k+1;j<=n;j++) {for (sum=0.0,i=k;i<=n;i++) sum += a[i][k]*a[i][j];tau=sum/c[k];for (i=k;i<=n;i++) a[i][j] -= tau*a[i][k];}}}d[n]=a[n][n];if (d[n] == 0.0) *sing=1;}The next routine, qrsolv, is used to solve linear systems.
In many applications only thepart (2.10.4) of the algorithm is needed, so we separate it off into its own routine rsolv.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).The standard algorithm for the QR decomposition involves successive Householdertransformations (to be discussed later in §11.2). We write a Householder matrix in the form1 − u ⊗ u/c where c = 12 u · u.
An appropriate Householder matrix applied to a given matrixcan zero all elements in a column of the matrix situated below a chosen element. Thus wearrange for the first Householder matrix Q1 to zero all elements in the first column of Abelow the first element. Similarly Q2 zeroes all elements in the second column below thesecond element, and so on up to Qn−1 .
Thus100Chapter 2.Solution of Linear Algebraic Equationsfor (j=1;j<n;j++) {Form QT · b.for (sum=0.0,i=j;i<=n;i++) sum += a[i][j]*b[i];tau=sum/c[j];for (i=j;i<=n;i++) b[i] -= tau*a[i][j];}rsolv(a,n,d,b);Solve R · x = QT · b.}void rsolv(float **a, int n, float d[], float b[])Solves the set of n linear equations R · x = b, where R is an upper triangular matrix stored ina and d. a[1..n][1..n] and d[1..n] are input as the output of the routine qrdcmp andare not modified. b[1..n] is input as the right-hand side vector, and is overwritten with thesolution vector on output.{int i,j;float sum;b[n] /= d[n];for (i=n-1;i>=1;i--) {for (sum=0.0,j=i+1;j<=n;j++) sum += a[i][j]*b[j];b[i]=(b[i]-sum)/d[i];}}See [2] for details on how to use QR decomposition for constructing orthogonal bases,and for solving least-squares problems.
(We prefer to use SVD, §2.6, for these purposes,because of its greater diagnostic capability in pathological cases.)Updating a QR decompositionSome numerical algorithms involve solving a succession of linear systems each of whichdiffers only slightly from its predecessor. Instead of doing O(N 3 ) operations each timeto solve the equations from scratch, one can often update a matrix factorization in O(N 2 )operations and use the new factorization to solve the next set of linear equations.
The LUdecomposition is complicated to update because of pivoting. However, QR turns out to bequite simple for a very common kind of update,A → A+s⊗t(2.10.7)(compare equation 2.7.1). In practice it is more convenient to work with the equivalent formA = Q·R→A0 = Q0 · R0 = Q · (R + u ⊗ v)(2.10.8)One can go back and forth between equations (2.10.7) and (2.10.8) using the fact that Qis orthogonal, givingt = v and eithers = Q · u oru = QT · s(2.10.9)The algorithm [2] has two phases. In the first we apply N − 1 Jacobi rotations (§11.1) toreduce R + u ⊗ v to upper Hessenberg form. Another N − 1 Jacobi rotations transform thisupper Hessenberg matrix to the new upper triangular matrix R0 .
The matrix Q0 is simply theproduct of Q with the 2(N − 1) Jacobi rotations. In applications we usually want QT , andthe algorithm can easily be rearranged to work with this matrix instead of with Q.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).void qrsolv(float **a, int n, float c[], float d[], float b[])Solves the set of n linear equations A · x = b. a[1..n][1..n], c[1..n], and d[1..n] areinput as the output of the routine qrdcmp and are not modified. b[1..n] is input as theright-hand side vector, and is overwritten with the solution vector on output.{void rsolv(float **a, int n, float d[], float b[]);int i,j;float sum,tau;2.10 QR Decomposition101#include <math.h>#include "nrutil.h"for (k=n;k>=1;k--) {Find largest k such that u[k] 6= 0.if (u[k]) break;}if (k < 1) k=1;for (i=k-1;i>=1;i--) {Transform R + u ⊗ v to upper Hessenberg.rotate(r,qt,n,i,u[i],-u[i+1]);if (u[i] == 0.0) u[i]=fabs(u[i+1]);else if (fabs(u[i]) > fabs(u[i+1]))u[i]=fabs(u[i])*sqrt(1.0+SQR(u[i+1]/u[i]));else u[i]=fabs(u[i+1])*sqrt(1.0+SQR(u[i]/u[i+1]));}for (j=1;j<=n;j++) r[1][j] += u[1]*v[j];for (i=1;i<k;i++)Transform upper Hessenberg matrix to upper trirotate(r,qt,n,i,r[i][i],-r[i+1][i]);angular.}#include <math.h>#include "nrutil.h"void rotate(float **r, float **qt, int n, int i, float a, float b)Given matrices r[1..n][1..n] and qt[1..n][1..n], carry out a Jacobi rotation√ on rowsi and i + 1√of each matrix.
a and b are the parameters of the rotation: cos θ = a/ a2 + b2 ,sin θ = b/ a2 + b2 .{int j;float c,fact,s,w,y;if (a == 0.0) {Avoid unnecessary overflow or underflow.c=0.0;s=(b >= 0.0 ? 1.0 : -1.0);} else if (fabs(a) > fabs(b)) {fact=b/a;c=SIGN(1.0/sqrt(1.0+(fact*fact)),a);s=fact*c;} else {fact=a/b;s=SIGN(1.0/sqrt(1.0+(fact*fact)),b);c=fact*s;}for (j=i;j<=n;j++) {Premultiply r by Jacobi rotation.y=r[i][j];w=r[i+1][j];r[i][j]=c*y-s*w;r[i+1][j]=s*y+c*w;}for (j=1;j<=n;j++) {Premultiply qt by Jacobi rotation.y=qt[i][j];w=qt[i+1][j];qt[i][j]=c*y-s*w;qt[i+1][j]=s*y+c*w;}}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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
