Heath - Scientific Computing (523150), страница 17

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 17)

Thus, theelimination matrix used for a given column vector a is of the form  1 ··· 0−m10 ··· 0a10.... . ...   ..   ..  ... . . . .... .  .   . .. 0 · · · 1 −m  0··· 0  ak−1   0 k−1  10 · · · 0   ak  =  ak  ,0 ··· 0   0 · · · 0 −mk+1 1 · · · 0   ak+1   0 . ..... . ...   .  ... . .... ..

  ..   .. ..00 · · · 0 −mn0 ··· 1anwhere mi = ai /ak , i = 1, . . . , n. This process requires about n3 /2 multiplications and asimilar number of additions, which is 50 percent more expensive than standard Gaussianelimination.During the elimination phase, the same row operations are also applied to the righthand-side vector (or vectors) of a system of linear equations. Once the elimination phasehas been completed and the matrix is in diagonal form, then the components of the solutionto the linear system can be computed simply by dividing each entry of the transformed righthand side by the corresponding diagonal entry of the matrix. This computation requires atotal of only n divisions, which is significantly cheaper than solving a triangular system, butnot enough to make up for the more costly elimination phase.

Gauss-Jordan eliminationalso has the numerical disadvantage that the multipliers can exceed 1 in magnitude even ifpivoting is used.Despite its higher overall cost, Gauss-Jordan elimination may be preferred in somesituations because of the extreme simplicity of its final solution phase. For example, it is52CHAPTER 2. SYSTEMS OF LINEAR EQUATIONSoccasionally advocated for implementation on parallel computers because it has a uniformworkload throughout the factorization phase, and then all of the solution components canbe computed simultaneously rather than one at a time as in ordinary back-substitution.Gauss-Jordan elimination is also sometimes used to compute the inverse of a matrixexplicitly, if desired. If the right-hand-side matrix is initialized to be the identity matrixI and the given matrix A is reduced to the identity matrix by Gauss-Jordan elimination,then the transformed right-hand-side matrix will be the inverse of A.

For computing theinverse, Gauss-Jordan elimination has about the same operation count as explicit inversionby Gaussian elimination followed by forward- and back-substitution.Example 2.8 Gauss-Jordan Elimination. We illustrate Gauss-Jordan elimination byusing it to compute the inverse of the matrix of Example 2.6. For simplicity, we omitpivoting.

We begin with the matrix A, augmented by the identity matrix I as right-handside, and repeatedly apply elimination matrices to annihilate off-diagonal entries of A untilwe reach diagonal form, then scale by the remaining diagonal entries to produce the identitymatrix on the left, and hence the inverse matrix on the right. 1 0 01 0 024 −2 1 0 02 4 −2 −2 1 0   41 −2 1 0  ,9 −3 0 1 0  =  0 10 151 0 1−2 −37 0 0 11 0 1 1 −4 02 4 −21 0 02 0 −69 −4 0 01 0  0 11 −2 1 0  =  0 11 −21 0 ,0 −1 10 151 0 10 043 −1 132731 02 0 −69 −4 02 0 0− 11222251  0 1 −1   0 11 −21 0  =  0 1 0 − 11,444 −40 010 043 −1 10 0 43 −11 12732 0 0− 112 0 022251  0 1 0   0 1 0 − 11=44 −410 0 40 0 43 −112731 0 0− 1127 −113444151  0 1 0 − 115 −1  ., so A−1 =  −1144 −443110 0 1−43 −11442.2.8Solving Modified ProblemsIn many practical situations linear systems do not occur in isolation but as part of a sequenceof related problems that change in some systematic way.

For example, one may need tosolve a sequence of linear systems Ax = b having the same matrix A but different righthand sides b. After having solved the initial system by Gaussian elimination, then the Land U factors already computed can be used to solve the additional systems by forwardand back-substitution. The factorization phase need not be repeated in solving subsequentlinear systems unless the matrix changes. This procedure represents a substantial savings2.2.

SOLVING LINEAR SYSTEMS53in work, since additional triangular solutions cost only O(n2 ) work, in contrast to the O(n3 )cost of a factorization.In fact, in some important special cases a new factorization can be avoided even whenthe matrix does change. One such case that arises frequently is the addition of a matrix thatis an outer product uv T of two nonzero vectors u and v. This is called a rank-one changebecause the outer product matrix uv T has rank one (i.e., only one linearly independentrow or column), and any rank-one matrix can be expressed as such an outer product of twovectors.

For example, if a single entry of the matrix A changes (say the (j, k) entry changesfrom ajk to ãjk ), then the new matrix is A − αej eTk , where ej and ek are the correspondingcolumns of the identity matrix and α = ajk − ãjk .The Sherman-Morrison formula gives the inverse of a matrix resulting from a rank-onechange to a matrix whose inverse is already known:(A − uv T )−1 = A−1 + A−1 u(1 − v T A−1 u)−1 v T A−1 ,where u and v are n-vectors.

Evaluation of this formula requires only O(n2 ) work (formatrix-vector multiplications) rather than the O(n3 ) work normally required for inversion.To solve a linear system (A−uv T )x = b with the new matrix, we could use the foregoingformula to obtainx = (A − uv T )−1 b = A−1 b + A−1 u(1 − v T A−1 u)−1 v T A−1 b.We would prefer to avoid explicit inversion altogether, however. If we have an LU factorization for A, then the following steps can easily be computed to obtain the solution to themodified system:1. Solve Az = u for z, so that z = A−1 u.2.

Solve Ay = b for y, so that y = A−1 b.3. Compute x = y + [(v T y)/(1 − v T z)]z.Note that the first step is independent of b and hence need not be repeated if there aremultiple right-hand-side vectors b. Again, this procedure requires only triangular solutionsand inner products, so it requires only O(n2 ) work and no explicit inverses.The Woodbury formula, in which u and v become n × k matrices U and V , generalizesthe Sherman-Morrison formula to a rank-k change in the matrix:(A − U V T )−1 = A−1 + A−1 U (I − V T A−1 U )−1 V T A−1 .Using similar techniques, it is possible to update the factorization rather than the inverse orthe solution. Caution must be exercised in using these updating formulas, however, becausein general there is no guarantee of numerical stability through successive updates as thematrix changes.Example 2.9 Rank-One Updating of Solutions.

To illustrate the use of the ShermanMorrison formula, we solve the linear system   24 −2x12 49 −3x2 =8,−2 −17x31054CHAPTER 2. SYSTEMS OF LINEAR EQUATIONSwhich is a rank-one modification of the system in Example 2.6 (only the (3, 2) entry haschanged). We take A to be the matrix of Example 2.6, so we can use the LU factorizationalready computed. One way to choose the update vectors is 00u =  0  and v =  1  ,−20so that the matrix of the new system is A − uv T , and the right-hand-side vector b has notchanged. We can use the previously computed LU factorization of A to solve Az = u toobtain z = [ − 23 12 − 12 ]T , and we had already solved Ay = b to obtain y = [ −1 2 2 ]T .The final step is then to compute the updated solution 3 −7−1−22vT y 1  =  4.z =  2 +x=y+211 − vT z1−22−102We have thus computed the solution to the new system without refactoring the modifiedmatrix.2.32.3.1Norms and Condition NumbersVector NormsTo measure errors and sensitivity in solving linear systems, we need some notion of the “size”of vectors and matrices.

The scalar concept of magnitude, modulus, or absolute value canbe generalized to the concept of norms for vectors and matrices. Although a more generaldefinition is possible, all of the vector norms we will use are instances of p-norms, which foran integer p > 0 and a vector x of dimension n are defined bykxkp =nX|xi |p!1/p.i=1Important special cases are as follows:• 1-norm:kxk1 =nX|xi |,i=1sometimes called the Manhattan norm because in the plane it corresponds to the distancebetween two points as measured in “city blocks.”• 2-norm:!1/2nX,kxk2 =|xi |2i=1which corresponds to the usual notion of distance in Euclidean space.2.3.

NORMS AND CONDITION NUMBERS55• ∞-norm:kxk∞ = max |xi |,1≤i≤nwhich can be viewed as a limiting case as p → ∞.All of these norms give the same results qualitatively, but in certain circumstances a particular norm may be easiest to work with analytically or computationally.

Характеристики

Список файлов книги