278 10.1 LU Decomposition and Matrix Inversion ‘This chapter deals with a class of elimination methods called LU decomposition tech- niques. The primary appeal of LU decomposition is that the time-consuming elimination step can be formulated so that it involves only operations on the matrix of coefficients, [A], Thus, it is well suited for those situations where many right-hand-side vectors (B) ‘must be evaluated for a single value of [A], Although there are a variety of ways in which this is done, we will focus on showing how the Gauss elimination method can be imple- mented as an LU decomposition. One motive for introducing LU decomposition is that il provides an efficient means to compute the matrix inverse. The inverse has a number of valuable applications in engineering practice. It also provides a means for evaluating system condition, LU DECOMPOSITION As described in Chap. 9, Gauss elimination is designed to solve systems of linear alge braic equations, [A]{X) = (8) (0.1) Although it certainly represents a sound way to solve such systems, it becomes ineficient ‘when solving equations with the same coefficients [A], but with different right-hand-side constants (the b's). Recall that Gauss elimination involves two steps: forward elimination and back- substitution (Fig. 9.3). Of these, the forward-elimination step comprises the bulk of the computational effort (recall Table 9.1), This is particularly true for large systems of equations. LU decomposition methods separate the time-consuming elimination of the matix [4] from the manipulations of the right-hand side (B). Thus, once [A] has been “decom posed,” muliple right-hand-side vectors can be evaluated in an ellicient manner. Interestingly, Gauss elimination itself can be expressed as an LU decomposition Before showing how this can be done, let us frst provide a mathematical overview of the decomposition strategy.10.1_1U DECOMPOSITION 279 10.1.1 Overview of LU Decomposition Just as was the case with Gauss elimination, LU decomposition requires pivoting to avoid division by zero. However, to simplify the following description, we will defer the issue of pivoting until after the fundamental approach is elaborated, In addition, the following ‘explanation is limited to a set of three simultaneous equations, The results can be directly extended to n-dimensional systems. Equation (10.1) can be rearranged to give [A](X} ~ (B} = 0 02) Suppose that Eq, (10.2) could be expressed as an upper triangular system Bead o Recognize that this is similar to the manipulation that occurs in the first step of Gauss climination, That is, elimination is used (o reduce the system to upper triangular form, Equation (10.3) can also be expressed in matrix notation and rearranged to give [UX — {(D} = 0 uo) Now, assume that there is a lower diagonal matrix with 1°s on the diagonal, 1 00 W1=]iy 1 0 «0. Is that has the property that when Eg, (10.4) is premultiplied by it, Eq. (10.2) is the result. That is, [LLUKX) (DY) = (AEX) — (8) 0.6) IC this equation holds, it follows from the rules for matrix multiplication that [10] = [4] aon and [KD = (8) os) A two-step strategy (see Fig. 10.1) for obtaining solutions can be based on Eqs. (10.4), (20.7), and 0.8) 1. LU decomposition step. [A] is factored or “decomposed” into lower [Z] and upper [U] triangular matrices. 2. Substitution step. (L] and |U} ate used to determine a solution (X} for a right-hand- side (B). This step itself consists of two steps. Fist, Eg. (10.8) is used to generate an intermediate vector {D} by forward substitution. Then, the result is substituted into Eq, (10.4), which can be solved by back substitution for (X} ‘Now, let us show how Gauss elimination can be implemented in this way.280 FIGURE 10.1 The sleps in IU decomposition. LU DECOMPOSITION AND MATRIX INVERSION, 10.1.2 LU Decomposition Version of Gauss Elimination Although it might appear at face value to be unrelated to LU decomposition, Gauss ‘climination can be used to decompose [A] into [Z] and [U]. This can be easily seen for IU], which is a direct product of the forward elimination, Recall that the forward- climination step is intended to reduce the original coefficient matrix [A] to the form ay an ay m-[9 ay “| 09) 0 0 ay ‘which is in the desired upper triangular format. ‘Though it might not be as apparent, the matrix [L] is also produced during the step. This can be readily illustrated for a thtee-equation system, S228 ‘The first step in Gauss elimination is to multiply row 1 by the factor [recall Eq. (9.13)] Juan and subtract the result from the second row to eliminate ay. Similarly, row 1 is multiplied by -m PeayEXAMPLE 10.1 10.1_1U DECOMPOSITION 281 and the result subtracted from the thitd row to eliminate as,. The final step is to multiply the modified second row by fa and subtract the result from the third row to eliminate ay Now suppose that we merely perform all these manipulations on the matrix [A], Clearly, if we do not want to change the equation, we also have to do the same to the right-hand side {B). But there is absolutely no reason that we have o perform the ma- nipolations simultaneously. Thus, we could save the f*s and manipulate (8) later. ‘Where do we store the factors fi, fu, and fiz? Recall thatthe whole idea behind the aa 3] eo Sa fia as. [4] >[L][U] (10.11), m/s aby “| 0 0 as and 1 0 0) we [i ro fa fo 1S ‘The following example confirms that [A] = [LJ[U1. LU Decomposition with Gauss Elimination Problem Slolemen!. Derive an LU decomposition based on the Gauss elimination formed in Example 9.5 Solulion, In Example 9.5, we solved the matrix 3-01 02 ol 7-03 03-02 10 la) = After forward elimination, the following upper tiangular matrix was obtained: 3-01 -02 [ul =| 0 7.00333 | -0.293333 o 0 10.0120282 LU DECOMPOSITION AND MATRIX INVERSION, ‘The factors employed to obtain the upper triangular matrix can be assembled into a lower triangular matrix. The elements a,, and ay, were eliminated by using the factors ou 03 fa = 52 = 003333333 fy, = 5 = 0.100000 and the element az was eliminated by using the factor 0.19 fu = SS = -0.0271300 7700333 Thus, the lower triangular matrix is 1 ° o [L] = | 0.0333333 1 ° 0.100000 -0.0271300. 1 Consequently, the LU decomposition is 1 ° Oya -04 -02 4] = (L][U] = | 0.033333, 1 0|/0 7.00333 — 0.293333 0.100000 -0.0271300 1],0 0 10.0120 ‘This result can be verified by performing the multiplication of [J[U] to give 3 -o1 -02 ) [Z11U] = | 0.099999 7 =~ -03 03 © 02 9.99996 | where the minor discrepancies are due to roundoff The following is pseudocode for a subroutine to implement the decomposition phase: ep 00 Ep 0 Notice that this algorithm is “naive” in the sense that pivoting is not included. This feature will be added later when we develop the full algorithm for LU decomposition. ‘After the matrix is decomposed, a solution ean be generated for a particular right- hhand-side vector (B}. This is done in two steps. First, a forward-substitution step is executed by solving Eq, (10.8) for (D). It is important to recognize that this merely10.1_1U DECOMPOSITION 283 amounts to performing the elimination manipulations on (B}. Thus, at the end of this step, the right-hand side will be in the same state that it would have been had we per- formed forward manipulation on [A] and (B) simultaneously. The forward-substitution step can be represented concisely as - Sah fori 23.40 aoa) The second step then merely amounts to implementing back substitution, as in Eq. (10.4). Again, it is important to recognize that this is identical to the back-substitution phase of conventional Gauss elimination. Thus, in a fashion similar to Eqs. (9.16) and (0.17), the back-substitution step can be represented concisely as By = dy dyn (00.13) a Daisy x Sood 0.14) EXAMPLENIO12) The Substitution Steps Problem Siaiement Complete the problem initiated in Example 10.1 by generating the final solution with forward and back substitution, Solution, As stated above, the intent of forward substitution is to impose the elimination ‘manipulations, that we had formerly applied to [A], on the righ-hand-side vector {B) Recall that the system being solved in Example 9.5 was [3-01 02) (x: 7.85 01 7-03) x; -193 Los -02 10 jks, na and that the forward-climination phase of conventional Gauss climination resulted in 3-01 -02 Jf 785 0 7.00333 -0.293333 |{ x2 = 4 -19.5617 102.) 0 0 10.0120 J lay, 70,0843, The forward-substitution phase is implemented by applying Eq. (10.7) to our problem, 1 0 0} (ar 7.85 0.033333 1 0] {dsp = 9-193 Lo.100000 -0.0271300 1} lay, nA or multiplying out the left-hand side, 4 = 785 00333333, + dd} = 193 O.1d; ~ 0.027134, + dy = 714284 LU DECOMPOSITION AND MATRIX INVERSION, We can solve the first equation for ds, dy = 7.85 which can be substituted into the second equation to solve for d; = ~19.3 ~ 0.0333333(7.85) = -19.5617 Both dy and d; can be substituted into the third equation to give dy = 14 ~ 0.1(7,85) + 0.02713( 19.5617) = 70,0843, Thus, 785 (D} =) -19.5617 70.0843, which is identical to the right-hand side of Eq. (10.2.1) This result can then be substituted into Eg, (10.4), [UK{X) = {D}. to give 3-01 -02 | fm 7.85 0 7.00333 -0,293333 | mp = { -19.5617 0 0 10.0120 | (ss 70.0843 which can be solved by back substitution (see Example 9.5 for details) for the final solution, 3 (x)=) -25 7.00003 ‘The following is pseudocode for a subroutine to implement both substitution phases: SUB Substitute (@, n, 8, x) “forward substitution x) = (0) — sumifay, E00 00 END Substitute10.1_1U DECOMPOSITION 285 The LU decomposition algorithm requires the same total multiply/divide flops as for Gauss elimination. The only difference is that a litle less effort is expended in the de- ‘composition phase since the operations are not applied to the right-hand side. Thus, the number of multply/divide flops involved in the decomposition phase can be calculated + O(n) (0.15) Conversely, the substitution phase takes a litle mote effort. Thus, the number of flops for forward and back substitution is n*, The total effort is therefore identical to Gauss elimination y73¢ 3 O(n?) (00.16) 10.1.3 LU Decomposition Algorithm An algorithm to implement an LU decomposition expression of Gauss elimination is listed in Fig. 10.2. Four features of this algorithm bear mention ‘+ The factors generated during the elimination phase are stored in the lower part of the matrix. This can be done because these are converted to zeros anyway and are unnecessary for the final solution, This storage saves space. + This algorithm keeps track of pivoting by using an order vector o. This greatly speeds up the algorithm because only the order vector (as opposed to the whole row) is pivoted. + The equations are not scaled, but scaled values of the elements are used to determine whether pivoting is to be implemented. + The diagonal term is monitored during the pivoting phase to detect near-zero occurrences in order to flag singular systems. If it passes back a value of er = —1, a singular matrix has been detected and the computation should be terminated. A parameter fol is set by the user to a small number in order (© detect near-zero 10.1.4 Crout Decomposition Notice that for the LU decomposition implementation of Gauss elimination, the [Z] matrix has 1's on the diagonal. This is formally referred to as a Doolittle decomposition, or fac- torization, An alternative approach involves a [U] matsix with 1's on the diagonal, This is called Crout decomposition. Although there are some differences between the approaches (Atkinson, 1978; Ralston and Rabinowitz, 1978), their performance is comparable. The Crout decomposition approach generates [U] and [L] by sweeping through the matrix by columns and rows, as depicted in Fig. 10.3. It can be implemented by the following concise series of formulas: asa fori 12.00 (0.17) forj = 23,00 (0.18)286 LU DECOMPOSITION AND MATRIX INVERSION, oreo UB Pivot (a, 2, s,m, K) Meee or et 19 = ABS(aacus.ad Sacur) it seit, 2 9 Bo DoFOR fi kL atemy = ABS(09i1)4/ 52) E10 Ledecom IF day > ig THEN SUP Decenpase (2, 9, tol, 0, 5, er) 219 = amy DorOR = 1,0 pai as eo IF $y = ABSC2,) evo 00 o0FOR J = 2, 0 tomy = 05 TF ABS(2,,;)>5) THEN sy = ABS(a,,3) o ew 00 diemy a0 0 a0 Pivot DOFOR k= 1, 0 — 1 SUB Substitute (2, 0, 0, b,x) Gis Pival(a, 0, s,m, &) took TF ABS (aaa, afSuns) < tOT THEN ne PRINT aoces,e/ Sore) mong ant * mu sein = si = ay o wo we ec ad rai factor = aural B= Dwailtvnn Sona 5 factor aire Fonte tot o eno 09 ven ah sa END 00 8 oo 5 88 fioy) < tl TBH few nD subst tate © PRINT agus alSo FIGURE 10.2 Pseudocode for an IU decomposition algorithm. @ For j y= ay— Shaws fori + dan 10.19) rioure 10.3 Bylay ehh ‘ ‘A schematic depicting the vt evaluations involved in Crout an Die {U decomposition we fork ftp am 020)10.2 THE MATRIX INVERSE 287 FIGURE 10.4 Peeudacade for Crout's (U decomposiion lgorihm DOOR k= 1,.0-1 10.2 to Sat asap ‘Aside from the fact that it consists of a few concise loops, the foregoing approach also hhas the benefit that storage space can be economized. There is no need to store the 1's on the diagonal of [U] or the 0's for [Z] or [U] because they are givens in the method, Con- sequently, the values of [U] ean be stored in the zero space of [L]. Further, close examina tion of the foregoing derivation makes it clear that alter each element of [A] is employed once, it is never used again, Therefore, as each element of [L] and [U] is computed, it can be substituted for the corresponding element (as designated by its subscripts) of [A] Pseudocode to accomplish this is presented in Fig. 10.4. Notice that Eq. (10.17) is not included in the pseudocode because the first column of [] is already stored in [A]. Otherwise, the algorithm directly follows from Bqs. (10.18) through (10.21). THE MATRIX INVERSE In our discussion of matrix operations (Sec. PT3.2.2). we introduced the notion that if a matrix (A] is square, there is another maui, {A]"', called the inverse of [A], for which (Eq, @T3.3)] (aay [Ay'1A] =288 LU DECOMPOSITION AND MATRIX INVERSION, EXAMPLE 10.3 Now we will focus on how the inverse can be computed numerically. Then we will explore how it can be used for enginecring analysis. 10.2.1 Calculating the Inverse The inverse can be computed in a column-by-column fashion by generating sotutions with unit vectors as the right-hand-side constants. For example, if the right-hand-side constant has a 1 in the first position and zeros elsewhere, 1 {8) = 0 0 the resulting solution will be the first column of the matrix inverse, Similarly, if a unit veetor with a 1 at the second row is used 0 (by = 41 0 the result will be the second column of the matrix inverse, ‘The best way to implement such a calculation is with the LU decomposition algorithm. described at the beginning of this chapter. Recall that one of the great strengths of LU decomposition is that it provides a very efficient means to evaluate multiple right- hand-side vectors, Thus, itis ideal for evaluating the multiple unit vectors needed to compute the inverse. Matrix Inversion Problem Statement. Employ LU decomposition to determine the matrix inverse for the system from Example 10.2. 3-01 02 lj=)01 7-03 03-02 10 Recall that the decomposition resulted in the following lower and upper triangular matrices: 3-01 -02 1 0 0 [ul =| 7.00333 —0.293333) [£] = | 0.033333 1 0 o 0 10.0120 0.100000 —0.0271300 1) Solution, The fist column of the matrix inverse can be determined by performing the forward-substitution solution procedure with a unit vector (with 1 in the first row) as the right-hand-side vector. Thus, Bg, (10.8), the lower-siangular system, ean be set up a3 1 0 of (a) 2 0.033333 1 ol fap = fo o.100000 00271300 1} la) bo10.2 THE MATRIX INVERSE 289 and solved with forward substitution for (D}7 = |1_ 0.03333 -0.1008 |. This vector ccan then be used as the right-hand side of Bq. (10.3), -02 J fx 1 =0,293333 |x p = 4 -0.03333, 100120 | (x, =0.1009 which can be solved by back substitution for {X}¥ = [0.33249 -0.00518 0.01008 ‘which is the frst column of the matrix, 0.33249 0 0] [art = | -000518 0 0 001008 0 0) ‘To determine the second column, Eq, (10.8) is formulated as 1 0 071 (4 0 010333333 1 olay st 100000 -o0271300 1)la) bo ‘This ean be solved for (D), and the results are used with Eq. (10.3) t0 determine {497 = [0.004944 0.142903 0.00271 J, which is the second column of the matrix, 0.33249 0.004944 0 [Al" = | -0.00518 0.142903 0 =0.01008 0.00271 0, Finally, the forward- and back-substitution procedures can be implemented with {B}"=|0 0 1) to solve for {X}" = [0.006798 0.004183 0.09988 |, which is the final column of the matrix, 0.33249 0.004944 0.006798 [ar = | -0.00518 0.142903 0.004183 0.01008 0,00271 0.09988 ‘The validity of this result can be checked by verifying that [AJIAI-" = (1) Pseudocode to generate the matrix inverse is shown in Fig. 10.5, Notice how the decomposition subroutine from Fig. 10.2 is called to perform the decomposition and then generates the inverse by repeatedly calling the substitution algorithm with unit vectors, ‘The effort required for this algorithm is simply computed as. noe’) (10.22) = decomposition +» X substitutions where from Sec. 10.1.2, the decomposition is defined by Eq. (10.15) and the effort in- volved with every right-hand-side evaluation involves n® multiply/divide flops.290 LU DECOMPOSITION AND MATRIX INVERSION, CALL Decampese (a, n, tol, 0, 5. er) IF er = 0 THEW boroR i = 1, 1 DoFOR § = 1, 1 IF t= 3 THEW biy= 4 FISe b(5) = 0 eno iF 00 00 CALL Substitute f@, 0, n, b,x) BoFOR J = 3, 1 ai. 1 = x) ND 00 N06 00 Output ai, if desired fuse PR conditioned systen" 0 IF FIGURE 10.5 Driver program that uses some ofthe subprograms fiom Fig. 10.2 lo genercle a mat’x inverse 10.2.2 Stimulus-Response Computations AAs discussed in Sec. PT3.1.2, many of the linear systems of equations confronted in engi- neering practice are derived from conservation laws, The mathematical expression of these laws is some form of balance equation to ensure that a particular property—mass, force, heat, momentum, or other—is conserved. For a force balance on a structure, the properties might be horizontal or vertical components of the forces acting on each node ofthe structure (see Sec. 122). For a mass balance, the properties might be the mass in each reactor of a chemical process (see Sec, 12.1). Other fields of engineering would yield similar examples ‘A single balance equation can be written for each part of the system, resulting in a set of equations defining the behavior of the property forthe entite system. These equa- tions are interrelated, or coupled, in that each equation may include one or more of the variables from the other equations. For many cases, these systems are linear and, there- fore, of the exact form dealt with in this chapter: (A(X) = (B) 023) Now, for balance equations, the terms of Eq. (10.23) have a definite physical interpreta tion. For example, the elements of (X) are the levels of the property being balanced for each part of the system. Ina fore balance ofa structure, they represent the horizontal and vertical forces in each member. For the mass balance, they are the mass of chemical in each reactor. In either case, they represent the system's state or esponse, which we are trying to determine. The right-hand-side vector {B) contains those elements of the balance that are in- dependent of behavior of the system—ihat is, they are constants. As such, they often represent the external forces or stimuli that drive the system.10.3 ERROR ANALYSIS AND SYSTEM CONDITION, 291 10.3 Finally, the matrix of coefficients [A] usually contains the parameters that express hhow the parts of the system interact or are coupled. Consequently, Eq. (10.23) might be reexpressed as [Interactions] {response} = { stimuli) ‘Thus, Bg, (10.23) can be seen as an expression of the fundamental mathematical model that we formulated previously as a single equation in Chap. | [recall Eq. (1.1)]. We ean now see that Eq. (10.23) represents a version that is designed for coupled systems involv. ing several dependent variables (X). ‘As we know from this chapter and Chap. 9, there arc a variety of ways to solve Eq, (10.23). However, using the matrix inverse yields a particularly interesting result. ‘The formal solution can be expressed as (x) = (AB) or (recalling our definition of matrix multiplication from Box PT3.2) wi = aiihy + aitbs + ats w= agi, + aise + be + a5itby ays ails + af ‘Thus, we find that the inverted matrix itself, aside from providing a solution, has ex- tuemely useful properties. That is, each of its elements represents the response of a single part of the system (o a unit stimulus of any other part of the system, Notice that these formulations are linear and, therefore, superposition and propor- lionality hold. Superposition means that if a system is subject o several different stimuli (the 6's), the responses can be computed individually and the results summed to obtain a total response. Proportionality means that multiplying the stimuli by a quantity results, in the response to those stimuli being multiplied by the same quantity. Thus, the cocf- ficient aj is a proportionality constant that gives the value of x; due to a unit level of +b, This result is independent of the effects of b, and b, on x;, which are rellected in the coefficients aj;! and aj;', respectively. Therefore, we can draw the general conclusion that the element aj* of the inverted matrix represents the value of x due to a unit quan- tity of by. Using the example of the structure, clement aj'of the matrix inverse would represent the force in member i due to a unit extemal force at node j. Even for small systems, such behavior of individual stimulus-response interactions would not be intui- lively obvious, As such, the matrix inverse provides a powerful technique for understand- ing the interrelationships of component parts of complicated systems, This power will be demonstrated in Secs. 12.1 and 12.2. ERROR ANALYSIS AND SYSTEM CONDITION Aside from its engineering applications, the inverse also provides a means to discern whether systems are ill-conditioned, Three methods are available for this purpose: 1. Scale the matrix of coefficients [A] so that the largest element in each row is 1. Invert te scaled matrix and if there are elements of [AJ that are several orders of magnitude sreater than one, it is likely thatthe system is il-conditioned (See Box 10.1).292 LU DECOMPOSITION AND MATRIX INVERSION, ox 10.1 of IlkConditioning ‘One method for assessing a system’s condition is to scale [A] so ‘that the largest element in each row is 1 and then compute [A] clement of [A] are several orders of magnitude greater than the clements ofthe original scaled matrix, iis likely thatthe system is illeconditoned. Insight into thi approach can be gained by recelling that away to check whether an approximate solution {X} is acceptable is to substitute it into the original equations and see whether the origi- sal ight-hand-side constants tesult, This is equivalent to 3) = A) (B10.1.1) where (R) isthe residual etween the right-hand-ide constants and the values computed with the solution (). IF (R) #8 smal, we sight conse tht the (1) values ate adequate. However, sappose thas (X) ithe ex «oy olution that yields a zero residual, as in (B) - A) «B10.12) Interpreting the Elements of the Matrix Inverse as a Measure Sobuatisg Ba (B10.12) om (B101.1 ils (R= arf ~ (a) ‘Multiplying both sides of this equation by [AI gives (x) = 0) = aR) “This result indicates why checking a solution by substitution can bbe misleading, For cases where elements of (Al are lange, @ sonal discrepancy inthe ight-hand-sde residual (R} could cor- respond to «lage exror (X) ~ (X) inthe calculate value ofthe unknowns. In other words, a small residual doesnot guarantee an sccurate solution, However, we can conclude that if the lazgeat clement of [AJ is onthe order of magnitude of unity. the systems can be considered to he well-sonditioned. Conversely, if (4|"" includes elements much larger than wnity, we conclude thatthe system is ill-conéitioned 2, Multiply the inverse by the original coefficient matrix and assess whether the result is close to the identity matrix. If not, it indicates ill-conditioning, 3. Invert the inverted matrix and assess whether the result is sulliciently close to the original coefficient matrix. If not, it again indicates that the system is ill-conditioned, Although these methods can indicate ill-conditioning, it would be preferable to ob- tain a single number (such as the condition number from Sec. 4.2.3) that could serve as, an indicator of the problem, Attempts to formulate such a matrix condition number are based on the mathematical concept of the norm, 10.3.1 Vector and Matrix Norms A norm is a real-valued function that provides a measure of the size or “length” of multicomponent mathematical entities such as vectors and maltices (see Box 10.2) A simple example is a vector in three-dimensional Euclidean space (Fig. 10.6) that can be represented as Fl=|la be where a, b, and c are the distances along the x, y, and z axes, respectively. The length of this vector—that is, the distance from the coordinate (0, 0, 0) to (a, b, e)—can he simply computed as Fle= Ve th +e where the nomenclature ||F||, indicates that this length is referred to as the Euclidean norm of [F]10.3 ERROR ANALYSIS AND SYSTEM CONDITION, 293 Box 10.2 Matrix Norms ‘As developed inthis section, Buclidean norms ean be employed to quantify the size of a vector, Por vectors, there are alleratives called p norms that ean be represented generally by IX = ( = r)" We can also see thatthe Euclidean norm and the 2 norm, |X] are identical for vectors. ‘Other important examples are Which represents the norm asthe sum of the absolute values ofthe elements. Another is the maximum-magnitude or uniform-vector UX = pay Dsl Which defines the norm as the element with the largest absolute value Using a similar approach, norms ean be developed for matrices For example, Wal = pox Se ‘That is, «summation of the absolute values of the coeficiens is performed for each column, and the largest ofthese summations is {aken a8 the norm. This is called the column sum norm. ‘A similar determination can be made forthe sows, resulting in a tuniformematrix oF row-sum norm, lial ass Diol i shouldbe noted tha, in contrast to vectors, tbe 2 norm and the Euclidean norm foramatix are not the same. Whereas the Euclidean norm {ll can be easily determined by Eq. (10.24), the matrix 2 norm |All is calculated as IAs = on) Where fs isthe largest eigenvalue of (AVT [A]. In Chap. 27, we will learn mare about eigenvalues. For the time being, the impar- tant point is that the Al , or speetral, norm isthe minimum norm and therefore, provides the tightest measure of size (Oxtegs 1972), FIGURE 10.6 Graphical deoiction of o vector [F]=|a Belin Evcldeon294 LU DECOMPOSITION AND MATRIX INVERSION, Similarly, for an n-dimensional vector |X| = |x; x2 “* %,, 8 Buclidean norm would be computed as iXlle ‘The concept can be extended further to a matrix [A], as in Alle = 024) which is given a special name—the Frobenius norm. However, as with the other vector norms, it provides a single value to quantify the “size” of [A] It should be noted that there are alternatives (o the Euclidean and Frobenius norms (see Box 10.2). For example, a uniform vector norm is defined as Xe = That is, the element with the largest absolute value is taken as the measure of the veetor's size. Similarly, a uniform matrix norm of row-sum norm is defined as IAL = pax > lo] o2s) In this case, the sum of the absolute value of the elements is computed for each row, and the largest of these is taken as the norm, Although there are theoretical benefits for using certain of the norms, the choice is sometimes influenced by practical considerations. For example, the uniform-row norm is widely used because of the ease with which it can be calculated and the fact that it usu ally provides an adequate measure of matrix size. 10.3.2 Matrix Condition Number Now that we have introduced the concept of the norm, we can use it to define Cond [A] = \1All- JA“) (10.26) where Cond [A] is called the matrix condition number. Note that for a matsix [A], this umber will be greater than or equal (0 1. It can be shown (Ralston and Rabinowitz, 1978: Gerald and Wheatley, 2004) that HAxll og, <> = Cond [4] x (Al) That is, the relative error of the norm of the computed solution can be as large as the relative error of the norm of the coefficients of [A] multiplied by the condition number. For example, if the coefficients of [A] are known to rigit precision (that is, rounding errors are on the order of 10) and Cond [A] = 10%, the solution [X] may be valid to only 1 ~ ¢ digits (sounding errors ~10EXAMPLE 10.4 10.3 ERROR ANALYSIS AND SYSTEM CONDITION, 295 ‘Matrix Condition Evaluation Problem Stolomen|. The Hilbert matrix, which is notoriously ill-conditioned, can be represented generally as 1 12 1/3 In 1213 1/4 Vat) I/m Af(n +1) A/a +2) 1/Qn- 1) } Use the row-sum norm to estimate the matrix condition number for the 3 X 3 Hilbert matrix, 1 4/2 4/3 [l= 1/2 1/3 1/4 1/3 1/4 as Soluiion, First, the matrix can be normalized so that the maximum element in each row is 1 142 1/3 lay |i 2/3 4/2 1 3/4 3/5 Summing each of the rows gives 1,833, 2.1667, and 2.35. Thus, the third row has the largest sum and the row-sum norm is The inverse of the scaled matrix can be computed as 9-18 10 fart =|-36 96-60 30-90 60, Note that the elements of this matrix are larger than the original matrix. This is also reflected in ils row-sum norm, which is computed as JA7tla. = [=36] + [96] + [60] = 192 ‘Thus, the condition number can be calculated as Cond [A] = 2.35192) = 451.2 ‘The fact that the condition number is considerably greater than unity suggests that the system is ill-conditioned. The extent of the ill-conditioning can be quantified by calculating ¢ = log 451.2 = 2.65. Computers using IEEE floating-point representation296 LU DECOMPOSITION AND MATRIX INVERSION, have approximately 1 = log 2” = 7.2 significant base-10 digits (recall Sec. 3.4.1). Therefore, the solution could exhibit rounding errors of up to 10°*7? = 3 x 107 Note that such estimates almost always overpredict the actual ertor. However, they are useful in alerting you to the possiblity that round-off errors may be significant Practically speaking, the problem with implementing Eg. (10.26) is the computa- tional price required to obtain | A~'|. Rice (1983) outlines some possible strategies to mitigate this problem. Further, he suggests an altemative way to assess system condi- tion: run the same solution on two different compilers. Because the resulting codes will likely implement the arithmetic differently, the effect of ill-conditioning should be evi- dent from such an experiment. Finally, it should be mentioned that software packages such as MATLAB sofiware and Mathcad have the capability o conveniently compute matrix condition. We will review these capabilities when we review such packages at the end of Chap. 11 10.3.3 Iterative Refinement In some eases, round-olf errors can be reduced by the following procedure. Suppose that we are solving the following set of equations: by at days = by oan aux + an + am aux; + angry + anny = By For conciseness, we will limit the following discussion to this small (3 X 3) system, However, the approach is generally applicable to larger sets of linear equations. Suppose an approximate solution vector is given by {X}" = ¥ % %: |: This solution can be substituted into Eq. (10.27) to give anh + ak + ak = by ak, + danke + agyhy = by (20.28) anh + asaks + ayyky = Now, suppose that the exact solution {X) is expressed as a function of the approximate solution and a vector of correction factors (AX), where ay = + Ay pt Ax 10.29) Mah + Ay If these results are substituted into Eq, (10.27), the following system results: ay (% + Ax) + anle + Ay) + a(t + Ax) = ayy(¥, + AX) + ay + Bag) + ays(dy + Axs) ~ By (1030) ay (iy + Ax) + anGy + Avy) + an(iy + Avy)PROBLEMS 297 Now, Eq. (10.28) can be subtracted from Eg, (10.30) to yield aya + andn + adn = by — by = By ayn, + andes + asx, = by ~ By = By (03) ayyAay + ann + ayydxy = by — by = By This system itself is a sot of simultaneous linear equations that can be solved (0 obtain the cortection factors. The factors can then be applied to improve the solution, as specified by Bq, (10.29) Its relatively straightforward to integrate an iterative refinement procedure into com- puter programs for elimination methods, I is especially effective forthe LU decomposition approaches described earlier, which are designed to evaluate different right-hand-side vec- tors efficiently. Note that to be effective for correcting ill-conditioned systems, the £°s in Eq, (10.31) must be expressed in double precision, PROBLEMS 10.1 Use the rues of matsix muliplication to prove that Eqs. (10.7) and (10.8) follow from Eq, (10.6) 10.2 (a) Use naive Gauss elimination to decompose the following system according to the description in See. 10.1.2. 10x + 2m — ay = 27 3x, = 6m + 2x5 = 615 xy tay + $n 215 “Then, mulply the resulting 2] and (U) maces to determine that {A} ss produced. (by Use LU decomposition to solve the system. ‘Show al te step in the computation, () Also solve the system for fn allemative rghtshand-ide vector: (B)* = [12 18 ~6) 103 (a) Solve the following system of equations by LV decomposition without pivoting By t4y— y= Atn tS tas dy =u + by =7 (b) Determine the matrix inverse, Check your esulls by verifying ‘hat [AIlA)~" = [1 104 Solve the following system of equations using LU decompo- sition with partial pivoting: 10,5 Determine the total flops as a function of the number of equations n for the (a) decomposition, (b) forward-substitution, and (e) back-substitution phases of the LU decomposition version ff Gauss elimination, 10.6 Use LU decomposition o determine the matin inverse forthe following system, Do not use a pivoting strategy, and check your results by verifying that [AILA}"® = (J 10x; + 2a — ay = 27 3x, — 6x + 25 = 61S ay tag + Sy = 215 10.7 Pesform Crout decomposition on ‘Then, multiply the resulting (L] and [U7] matrices to determine that [A] is produced. 10.8 The following system of equations is designed to determine concentrations (he c's in g/m’) in a eres of coupled reactors as a function of the amount of mass npato each reactor (the igh-hand sides in g/439), Seq — Sez — 6 = 3800 Bey + 18; — 6c, = 1200 ey — e+ Wey = 2350 (2) Determine the matrix inverse, ()) Use the inverse to determine the solution (6) Determine how mich the rate of mass input to reactor § raust be increased to induce a 10 gi’ rise inthe concentration of eactr 1. (@) How much will he concentration in reactor 3 be reduced ifthe rate of mass input to reactors I and 2 is reduced by 500 and 250 piday, respectively?298 10.9 Solve the following set of equations with LU decomposition: Bx -2e bay = 10 Day Fxg — 4 = 48 may — Bag} Sey = —26 10.10 (@) Determine the LU decoraposition without piveting by ‘hand forthe following matrix and check your results by validating ‘hat [LIU] = (A. s 21 372 2309 (b) Eraploy the result of (a) to compute the determinant. (6) Repeat (a) and (b) using MATLAB. 10.11 Use the following LU decomposition to (a) compute the de- terminant and (b) solve [A]{x) = {5} with {b}" = |-10.44 —26) 1 32 01 [a] = Eytw) = foss67 13333 — 4.6667 0.3333 —03636 1 3.6364 10.12 Determine A Al ad Al for s 2 -10 [aj=)-9 103 1-16 Seale the matsix by making the maximum element in each row equal to one. 10.13 Determine the Frobenius and the row-sum norms for the systems in Probs, 10.3 and 10.4, Scale the matrices by making the ‘maximum element in each row equal t one, 10.14 A matix [A] is defined as ous 025s ay =| 0018625 0825028 0.00463 002777 0.16667 1 0.001983 0.015625 0.125 1 Using the coluran-sum norm, compute the condition number and how many suspect digits would be generated by this matrix. 10.15 (@) Determine the condition aumber for the following system using the row-sum norm, Do not normalize the system, 1 4 9 16 25 49 16 25 36 2 16 25 36 49 16 25 36 49 6s 25 36 49 64 81 LU DECOMPOSITION AND MATRIX INVERSION, How many digits of precision will be lost due to ill-conditioning? (b) Repeat (a), but scale the matix by making the maximum ele- sent in each row equal to one 40.16 Determine the condition number based on the sow-sum nonin forthe normalized § X 5 Hilbert matrix. How many signili- cant digits of precision willbe lost due lo il-conditioning? 410.17 Besides the Hilbert matrix, there are ather matrices that are inherently ill-conditioned. One such case is the Vandermonde ‘matrix, which has the following forms Ee] (@) Determine the condition number based on the row-sum norm forthe case where xy = 4, = 2, and xy = 7. (b) Use MATLAB or Mathcad software to compute the spectral and Frobenius condition numbers 40.18 Develop a user-friendly program for LU decomposition based an the pseudocode from Fig. 10.2 10.19 Develop a user-friendly program for LU decomposition, cluding the capability to evaluate the matrix inverse, Base the pro- gram on Figs. 10.2 and 10.5. 110.20 Use iterative refinement techniques to improve 1 = 2, 3, and.x, = 8, which are approximate solutions of 2x + Sy tn = —5 2 Sx 4 2y tay at Qu tas 110.21 Consider vectors Ana 3] bak Beavis j-4k, Bi tej + 2k Vecioy A, is perpendicular to Bas well as to C . tis also known ‘hat BC = 2, Use any method studied inthis chapter to solve for ‘he thiee unknowns, a, bande 110.22 Consider the following vector: Anais bjs Bana + jak Cai asp 4k where Ais an unknown vector. IE Axi) + AX C)= 60461 + Ob-D] + Che + DE tase any method leamed in this chapter to salve for the three un knowns, a,b, and ¢PROBLEMS 299 10.28 Let the function be defined om the interval [0,2] as follows: arth Osxs1 oetd 1Sx52 Determine the constants a,b ¢, and dso that the function fsatistes the following: + AO) = $2) = 1 + Fis continuous on the entire interval. satbed, Jey = { Derive and solve a system of linear algebraic equations witha ma tix form ideatieal t By. 10.1), 1024 (@) Create 43 X 3 Hilbert matsx. Ths will be your matsx (A) Muliply the matrix by the column veewr (x) ~ (1 1,1] The solution of [A]{x) will be another column vector (b). Using any numerical package and Gauss elimination find the solution {oTAltx) = {busing the Hilber matrix andthe vector (8} tht you calculated. Compare the result o your known (x) vector Use sufficient precision in displaying results to allow you to detect imprecision. (b) Repeat part (a) using a7 7 Hilbert matrix (6) Repeat part (a) using a 10 X 10 Lilbert matin 10.25 Polynomial interpolation consists of determining the unique (= Disorder polynomial tha fits m data points. Such polynomi- als have the general form, POD =P" + po + Pease t Pe (P1025) where the p's are constant coefficients, A straightforward way for computing the coefficients is to generate n linear algebraic equations that we can solve simultaneously for the coefiients. Suppose that we want todetermine the coefiiens ofthe fourth-order polynomial Sl3) = pix’ + pax? + psx? + pax + ps that passes through the following five points: (200, 0.746, (250, 0.675), (300, 0.616), (40, (0,525), and (500, 0457), Each ofthese pairs canbe substituted into Eq, (P10.25) to yield a system of five equations wit five unknowns (the p's) Use this approach to solve forthe coeticients. In adition, determine and interpret the condition nurber