c2-3 (779461)
Текст из файла
432.3 LU Decomposition and Its ApplicationsIsaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §2.1.Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), §2.2.1.Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New York: Wiley).Suppose we are able to write the matrix A as a product of two matrices,L·U=A(2.3.1)where L is lower triangular (has elements only on the diagonal and below) and Uis upper triangular (has elements only on the diagonal and above). For the case ofa 4 × 4 matrix A, for example, equation (2.3.1) would look like this: α11 α21α31α410α22α32α4200α33α430β110 0· 00α440β12β2200β13β23β330β14β24 β34β44a11a= 21a31a41a12a22a32a42a13a23a33a43a14a24 a34a44(2.3.2)We can use a decomposition such as (2.3.1) to solve the linear setA · x = (L · U) · x = L · (U · x) = b(2.3.3)by first solving for the vector y such thatL·y=b(2.3.4)U·x=y(2.3.5)and then solvingWhat is the advantage of breaking up one linear set into two successive ones?The advantage is that the solution of a triangular set of equations is quite trivial, aswe have already seen in §2.2 (equation 2.2.4).
Thus, equation (2.3.4) can be solvedby forward substitution as follows,y1 =b1α11i−1X1 yi =αij yj bi −αiij=1(2.3.6)i = 2, 3, . . . , Nwhile (2.3.5) can then be solved by backsubstitution exactly as in equations (2.2.2)–(2.2.4),yNxN =βNNN(2.3.7)X1 xi =βij xj yi −i = N − 1, N − 2, . . . , 1βiij=i+1Sample 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).2.3 LU Decomposition and Its Applications44Chapter 2.Solution of Linear Algebraic EquationsPerforming the LU DecompositionHow then can we solve for L and U, given A? First, we write out thei, jth component of equation (2.3.1) or (2.3.2). That component always is a sumbeginning withαi1 β1j + · · · = aijThe number of terms in the sum depends, however, on whether i or j is the smallernumber. We have, in fact, the three cases,i<j:i=j:i>j:αi1 β1j + αi2 β2j + · · · + αii βij = aijαi1 β1j + αi2 β2j + · · · + αii βjj = aijαi1 β1j + αi2 β2j + · · · + αij βjj = aij(2.3.8)(2.3.9)(2.3.10)Equations (2.3.8)–(2.3.10) total N 2 equations for the N 2 + N unknown α’s andβ’s (the diagonal being represented twice).
Since the number of unknowns is greaterthan the number of equations, we are invited to specify N of the unknowns arbitrarilyand then try to solve for the others. In fact, as we shall see, it is always possible to takeαii ≡ 1i = 1, . . . , N(2.3.11)A surprising procedure, now, is Crout’s algorithm, which quite trivially solvesthe set of N 2 + N equations (2.3.8)–(2.3.11) for all the α’s and β’s by just arrangingthe equations in a certain order! That order is as follows:• Set αii = 1, i = 1, .
. . , N (equation 2.3.11).• For each j = 1, 2, 3, . . . , N do these two procedures: First, for i =1, 2, . . ., j, use (2.3.8), (2.3.9), and (2.3.11) to solve for βij , namelyβij = aij −i−1Xαik βkj .(2.3.12)k=1(When i = 1 in 2.3.12 the summation term is taken to mean zero.) Second,for i = j + 1, j + 2, . . . , N use (2.3.10) to solve for αij , namely1αij =βjjaij −j−1X!αik βkj.k=1Be sure to do both procedures before going on to the next j.(2.3.13)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).Equations (2.3.6) and (2.3.7) total (for each right-hand side b) N 2 executionsof an inner loop containing one multiply and one add. If we have N right-handsides which are the unit column vectors (which is the case when we are inverting amatrix), then taking into account the leading zeros reduces the total execution countof (2.3.6) from 12 N 3 to 16 N 3 , while (2.3.7) is unchanged at 12 N 3 .Notice that, once we have the LU decomposition of A, we can solve with asmany right-hand sides as we then care to, one at a time. This is a distinct advantageover the methods of §2.1 and §2.2.452.3 LU Decomposition and Its Applicationsacegibddifaghjonalsubdelemeniatsgonaxleleetc.mentsFigure 2.3.1.
Crout’s algorithm for LU decomposition of a matrix. Elements of the original matrix aremodified in the order indicated by lower case letters: a, b, c, etc. Shaded boxes show the previouslymodified elements that are used in modifying two typical elements, each indicated by an “x”.If you work through a few iterations of the above procedure, you will see thatthe α’s and β’s that occur on the right-hand side of equations (2.3.12) and (2.3.13)are already determined by the time they are needed.
You will also see that every aijis used only once and never again. This means that the corresponding αij or βij canbe stored in the location that the a used to occupy: the decomposition is “in place.”[The diagonal unity elements αii (equation 2.3.11) are not stored at all.] In brief,Crout’s method fills in the combined matrix of α’s and β’s,β11 α21α31α41β12β22α32α42β13β23β33α43β14β24 β34β44(2.3.14)by columns from left to right, and within each column from top to bottom (seeFigure 2.3.1).What about pivoting? Pivoting (i.e., selection of a salubrious pivot elementfor the division in equation 2.3.13) is absolutely essential for the stability of Crout’sSample 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).etc.x46Chapter 2.Solution of Linear Algebraic Equations#include <math.h>#include "nrutil.h"#define TINY 1.0e-20;A small number.void ludcmp(float **a, int n, int *indx, float *d)Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwisepermutation of itself.
a and n are input. a is output, arranged as in equation (2.3.14) above;indx[1..n] is an output vector that records the row permutation effected by the partialpivoting; d is output as ±1 depending on whether the number of row interchanges was evenor odd, respectively. This routine is used in combination with lubksb to solve linear equationsor invert a matrix.{int i,imax,j,k;float big,dum,sum,temp;float *vv;vv stores the implicit scaling of each row.vv=vector(1,n);*d=1.0;No row interchanges yet.for (i=1;i<=n;i++) {Loop over rows to get the implicit scaling informabig=0.0;tion.for (j=1;j<=n;j++)if ((temp=fabs(a[i][j])) > big) big=temp;if (big == 0.0) nrerror("Singular matrix in routine ludcmp");No nonzero largest element.vv[i]=1.0/big;Save the scaling.}for (j=1;j<=n;j++) {This is the loop over columns of Crout’s method.for (i=1;i<j;i++) {This is equation (2.3.12) except for i = j.sum=a[i][j];for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];a[i][j]=sum;}big=0.0;Initialize for the search for largest pivot element.for (i=j;i<=n;i++) {This is i = j of equation (2.3.12) and i = j +1 .
. . Nsum=a[i][j];of equation (2.3.13).for (k=1;k<j;k++)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).method. Only partial pivoting (interchange of rows) can be implemented efficiently.However this is enough to make the method stable.
This means, incidentally, thatwe don’t actually decompose the matrix A into LU form, but rather we decomposea rowwise permutation of A. (If we keep track of what that permutation is, thisdecomposition is just as useful as the original one would have been.)Pivoting is slightly subtle in Crout’s algorithm. The key point to notice is thatequation (2.3.12) in the case of i = j (its final application) is exactly the same asequation (2.3.13) except for the division in the latter equation; in both cases theupper limit of the sum is k = j − 1 (= i − 1). This means that we don’t have tocommit ourselves as to whether the diagonal element βjj is the one that happensto fall on the diagonal in the first instance, or whether one of the (undivided) αij ’sbelow it in the column, i = j + 1, . . .
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
