Steven C. Chapra Raymond P. Canale Berger Chair in Computing and Engineering Professor Emeritus of Civil Engineering Tufts University University of Michigan Numerical Methods for Engineers With Software and Programming Applications Fourth Edition Boston Burr Ridge, |L_ Dubuque, IA Madison, WI New York San Francisco St. Louis Bangkok Bogoté Caracas Kuala Lumpur Lisbon Londen Madrid Mexico City Milan Montreal New Delhi Santiago Seoul Singapore Sydney Taipei TorontoJind the’ imple- oftware CHAPTER 9 SEE Gauss Elimination ‘This chapter deals with simultaneous linear algebraic equations that can be represented generally as = aust + aary +++ + any = by 21% + ane + anny = eo» yt gata ++ + Oke Where the a’s are constant coefficients and the b’s are constants The technique described in this chapter is called Gauss elimination because it involves ‘combining equations to eliminate unknowns. Although itis one of the earliest methods for solving simultaneous equations, it remains among the most important algorithms in use today and isthe basis for inear equation solving on many popular software packages. 9.1 SOLVING SMALL NUMBERS OF EQUATIONS Before proceeding to the computer methods, we will describe several methods that are appropriate for solving small (n << 3) sets of simultaneous equations and that do not re- quire a computer. These are the graphical method, Cramer's rule, and the elimination of unknowns. 9.1.1 The Graphical Method ‘A graphical solution is obtainable for two equations by plotting them on Cartesian coordi- nates with one axis corresponding to x, and the other to x2. Because we are dealing with lin- ear systems, each equation is a straight line, This can be easily illustrated for the general equations ona + aja%2 = by nix + am%2 = ba| EXAMPLE 9.1 GAUSS ELIMINATION. Both equations can be solved for x2 a ‘Thus, the equations are now in the form of straight lines; that is, x2 = (slope) x1 + inter- cept. These lines can be graphed on Cartesian coordinates with x2 as the ordinate and x; as the abscissa. The values of x; and 2 at the intersection of the lines represent the solution, ‘The Graphical Method for Two Equations Problem Statement. Use the graphical method to solve 18 wo.) 2 9.12) Solution. Letx; be the abscissa, Solve Eq (E9.1.1) for x a1 +9 which, when plotted on Fig, 9.1, is straight line with an intercept of 9 and a slope of —3/2. FIGURE 9.1 Graphical solution ofa set of wo simultaneous lines rep a1 algebraic equations. The intersection oftheon) a1) 9.1_ SOLVING SMALL NUMBERS OF EQUATIONS Equation (E9.1.2) can also be solved for x2: meine 2 Which is also plotted on Fig. 9.1. The solution is the intersection of the two lines at xy = 4 and x2 = 3. This result can be checked by substituting these values into the original equa- tions to yield 344) +28) = 18 -@) +28) =2 ‘Thus, the results are equivalent to the right-hand sides of the original equations. For three simultaneous equations, each equation would be represented by a plane in 4 three-dimensional coordinate system, The point where the three planes intersect would ‘represent the solution. Beyond three equations, graphical methods break down and, conse- 4Quently, have little practical value for solving simultaneous equations, However, they some- times prove useful in visualizing properties of the solutions. For example, Fig. 9.2 depicts three cases that ean pose problems when solving sets of linear equations. Figure 9.24 shows the case where the two equations represent parallel lines. For such situations, there is no solution because the lines never cross. Figure 9.26 depicts the case where the two lines are coincident, For such situations there is an infinite number of solutions. Both types of systems are said to be singular. In addition, systems that are very close to being singular (Fig, 9.2c) can also cause problems. These systems are said to be ill-conditioned. Graphically, this corresponds to the fact that it is difficult to identify the exact point at ‘which the lines intersect. Ill-conditioned systems will also pose problems when they are encountered during the numerical solution of linear equations. This is because they will be extremely sensitive to round-off error (recall Sec, 4.2.3). FIGURE 9.2 dificult 0 detect visually Grophical depiction of singular and ikcondiioned systems: (al no soluton, (b infinite solutions, ‘rd (ileondltioned system where the slopes re s0 clase thatthe point of intersection IsEXAMPLE 9.2 GAUSS EUMINATION. 9.1.2 Determinants and Cramer's Rule Cramer's rule is another solution technique that is best suited to small numbers of eq tions, Before describing this method, we will briefly introduce the concept ofthe determi nant, which is used to implement Cramer's rule. In addition, the determinant has relevance to the evaluation of the ill-conditioning of a matrix, Determinants. The determinant can be illustrated for a set of three equations: [A}0) = (B) Where [A] isthe coefficient matrix: an an ay a) ‘The determinant D of this system is formed from the coefficients of the equation, as in D=|ax an an 02) Although the determinant D and the coefficient matrix [4] are composed of the same ele- ‘ments, they are completely different mathematical concepts. That is why they are distin- guished visually by using brackets to enclose the matrix and straight lines to enclose the determinant, In contrast to a matrix, the determinant is a single number. For example, the value of the second-order determinant p-|™ | is calculated by D= anon — ana 0 For the third-order case (Eq. (9.2) a single numerical value for the determinant can be ‘computed as D=ay =a 4) where the 2 by 2 determinants are called minors. Determinants Problem Statement. Compute values for the determinants of the systems represented it Figs. 9.1 and 9.2. Solution. For Fig. 9.1 3Q)- 2-1) =8e distin. ‘ose the. aple, the 69) vt can be EXAMPLE 9.3 9.1_ SOLVING SMALL NUMBERS OF EQUATIONS For Fig. 92a: 1/2 12 For Fig. 9.2b: 1/2 -1 For Fig. 9.2: | -72 D= 0.04 In the foregoing example, the singular systems had zero determinants. Additionally, the results suggest that the system that is almost singular (Fig. 9.2c) has a determinant that is close to zero. These ideas will be pursued further in our subsequent discussion of il- conditioning (Sec. 9.3.3), Cramer's Rule. This rule states that each unknown in a system of linear algebraic equa- tions may be expressed as a fraction of two determinants with denominator D and with the ‘numerator obtained from D by replacing the column of coefficients of the unknown in «vestion by the constants by, bay... By For example, x; would be computed as os) Cramer's Rule Problem Statement, Use Cramer's rule to solve 0.31; +0.52e +23 = 0.01 O15x1 4.12 + 1.935 = 0.67 O14; +0.3x2 +0.5; 0.44 Solution. ‘The determinant D can be written as (Eq. (9.2)] os 052 1 D=|os 1 19 01 03 05 The minors are (Ea. (9.3) 119 . 5) - 1.9003) = 007 AL loo fa 0.5) - 190.3) 42=|05 12| =0505) - 190.1) =006 01 0s(GAUSS EUMINATION. 5) 1) *=lo1 03] =0.5(0.3) — 1(0.1) = 0.05 ‘These can be used to evaluate the determinant, as in (Bq. (9.4)} D =0,3(-0.07) — 0.52(0.06) + 1(0.05) = ~0.0022 Applying Eq. (9.5), the solution is [oor 0.52 0.67 44 0s 067 OL -0.44 05) =0.0022 jo 052-001 os 1 067 jos 001 1 19| Joa_o3 | -nose Sioa Sh FFor more than three equations, Cramer's rule becomes impractical because, a8 the ‘number of equations increases, the determinants are time consuming to evaluate by hand (or by computer). Consequently, more efficient alternatives are used. Some of these alter natives are based on the last noncomputer solution technique covered in the next section— the elimination of unknowns. 9.1.3 The Elimination of Unknowns ‘The elimination of unknowns by combining equations is an algebraic approach that can be illustrated fora set of two equations: aux; + any ana + a2 00 “The basic strategy is to multiply the equations by constants so that one of the unknoWs willbe eliminated when the two equations are combined. The result is a single equatie? that can be solved for the remaining unknown. This value can then be substituted in. cither ofthe original equations to compute the other variable. For example, Eq, (9.6) might be multiplied by a2; and Eq. (9.7) by au to give ayant + aiadait2 = bras anyays4y + anayn¥2 = brantEXAMPLE 9.4 9.1_ SOLVING SMALL NUMBERS OF EQUATIONS tions to yield Smaaiit2 ~ ai2aa.%2 = beaiy — bya, Which can be solved for ay = itd = andy Fauation (9.10) can then be substituted into Eg. (9.6), which canbe solved for aby ~ ay. Novice that Eqs. (9.10) and (9.11) follow diretly from Cramer's ul, which states jor 6 bia — arabs Elimination of Unknowns Problem Statement. Use the elimination of uninowns to solve (real Example 91) 3x1 42m =18 an +2 =2 Solution. Using Eqs. (9.11) and (9.10), = 2018) ~2@) © 3)=2-) _3Q)- Cis = 32) — (-1) which is consistent with our graphical solution (Fig. 9.1, ‘The elimination of unknowns can be extended to systems with more than two or three equations. However, the numerous calculations that are required fot lamer systems make the method extremely tedious to implement by hand. However, as described in the noxe section, the technique can be formalized and reaiy programmed for he computer(GAUSS EUMINATION. 9.2 NAIVE GAUSS ELIMINATION In the previous section, the elimination of unknowns was used to solve a pair of simulta- neous equations. The procedure consisted of two steps: 1L. The equations were manipulated to eliminate one of the unknowns from the equations, ‘The result ofthis elimination step was that we had one equation with one unknown, | 2. Consequently, this equation could be solved directly and the result back-substinuted into one of the original equations to solve for the remaining unknown. ‘This basic approach can be extended to large sets of equations by developing a sys: tematic scheme ot algorithm to eliminate unknowns and to back-substitute. Gauss elimi- ration isthe most basic of these schemes. This section includes the systematic techniques for forward elimination and back sub- stitution that comprise Gauss elimination. Although these techniques are ideally suited for implementation on computers, some modifications will be required to obtain a reliable al- ‘gorithm. In particular, the computer program must avoid division by zero. The following ‘method is called “naive” Gauss elimination because it does not avoid this problem. Sub- sequent sections will deal with the additional features required for an effective computer | program, ‘The approach is designed to solve a general set of x equations: by (0.120) be 128) aya tank + ars tot aan any t+ anane + aaats toot dann ays ga %2 + Ops Ho + Manin = Oo @.122) 'As was the case with the solution of two equations, the technique for m equations consists ‘of two phases: elimination of unknowns and solution through back substitution Forward Elimination of Unknowns. The first phase is designed to reduce the set of equations to an upper triangular system (Fig. 9.3). The initial step will be to eliminate the first unknown, 1, from the second through the nth equations. To do this, multiply Eq, (9.122) by ai/ai: to give aan + PB aaa + Now this equation can be subtracted from Bq, (9.120) t give : +( - cn) 2 | digrr + 11+ Oighe = Be where the prime indicates that the elements have been changed from their original valoes “The procedure is then repeated for the remaining equations, For instance, Bg, (9.124) ccan be multiplied by a2;/ai; and the result subtracted from the third equation, Repeatinéq ___9.2_ NAIVE GAUSS EUMINATION 239 ted x FIGURE 9.3 mi The two phases of Gouss elim: rofon: forward elimination and ob. fock subsinsion. The primas for indicate the number of mes at batho coolicients and ing constant have been modified. ab ater the procedure forthe remaining equations results inthe following modified system: | 0) aux) Fane + ast + by 140) me ahgte + digts + % 140) digtr + ages + 6 (146) 129 sist gta + digs +--+ ay = by (14d) For the foregoing steps, Eq, (9.12a) is called the pivot equation and ay is ealled the pivot ne coefficient or element, Note that the process of multiplying the first row by ani/ay. is ethe equivalent to dividing it by az; and multiplying ity az). Sometimes the division operation wily is referred to as normalization. We make this distinction because a zero pivot element can interfere with normalization by causing a division by zero. We will return to this important issue alter we complete our description of naive Gauss elimination, cal ‘Now repeat the above to eliminate the second unknown from Eq, (9.14¢) through (9.144). To do this multiply Eq. (9.14) by a’, /a and subtract the result from Eq. (9.14c) Perform a similar elimination for the remaining equations to yield aust ant baits +++ ante =, algae ars + alues. 3.124) asa bo dy | eating where the double prime indicates that the elements have been modified twice. |‘GAUSS EUMINATION _ ‘The procedure can be continued using the remaining pivot equations. The final ma nipulation in the sequence is to use the (n — 1)th equation to eliminate the x,—1 term from the nth equation. At this point, the system will have been transformed to an upper triangu- Jar system (recall Box PT3.1): auiay tyne + ansty +--+ ante = BL (9.150) b, (9.156) 150) dhaka + agts FF diye = abyts toot alee af My 15a) Pseudocode to implement forward elimination is presented in Fig. 9.4a. Notice that three nested loops provide a concise representation of the process. The outer loop moves down the matrix from one pivot row to the next, The middle loop moves below the pivet row to each of the subsequent rows where elimination is to take place. Finally the inner ost loop progresses across the columns fo eliminate or transform the elements of & pat ticular row. Back Substitution. Equation (9.15d) can now be solved for o-0) bas 16 “This result can be back-substituted into the (n — I}th equation to solve for dure, which is repeated to evaluate the remaining x's, can be represented by the following |. The proce FIGURE 9.4 Pseudecade fo perform (al for ward elimination and (6) Bock substiution @ eno't0 by = by ~ “factor + by exo 00 exo 00 ® Xe = Do faa sun = 0 exo 00 0 00016 owing EXAMPLE 9.5 1 9.2_NAIVE GAUSS EUMINATION 241 formula: fori=n—1n 1 ol) Pseudocode to implement Eqs. (9.16) and (9.17) is presented in Fig. 9.4b, Notice the similarity between this pseudocode and that in Fig. PT34 for matrix multiplication. As ‘with Fig. PT3.4, a temporary variable, sum, is used to accumulate the summation from Eq. (9.17). This results ina somewhat faster execution time than ifthe summation were ac- cumulated in b; More importantly it allows efficient improvement in precision if the vari- able, sum, is declared in double pre Naive Gauss Elimination Problem Statement, Use Gauss elimination to solve 3x; O12 — 020 9.5.1) O.1xy + 7x2 03x) 3 £9.52) 03x, = 0.22 + 10x = 714 (953) Camry six significant figures during the computation, Solution. The first part of the procedure is forward elimination. Multiply Eq. (E9.5.1) by (0.1)/3 and subtract the result from Eq. (E9.5.2) to give 7.00333x) ~ 0.29333315 = -19.5617 ‘Then multiply Eg. (E9.5.1) by (0.3)/3 and subtract it from Eq, (E9.5.3) to eliminate x, After these operations, the set of equations is 3 -O.1x, — 0.2m = 7.85 E95.) 7.00333" ~ 0.29333313 = ~19.5617 (E953) =0,190000x; + 10.0200x3 = 70.6150 955) ‘To complete the forward elimination, x must be removed from Eq, (9.5.6). To accom- plish this, multiply Eq, (E9.5.5) by ~0.190000/7.00333 and subtract the result from Eq, (E9.5.6). This eliminates x from the thi equation and reduces the system to an upper ‘wiangular form, as in 3x, -Olm 02x =7.85 95.) 7.003332 ~ 0.293333x3 = —19,5617 (£9.58) 10.0203 = 70.0843 e959) We can now solve these equations by back substitution. First, Eq, (E9.5.9) can be solved for 70.0843 = 70,0200 7.00003 95.10)GAUSS EUMINATION ‘This result can be back-substituted into Eg. (E9.5.8) 7.00333x2 ~ 0.293333(7.00003) = —19.5617 which can be solved for =19.5617 + 0.293333(7.00003) : 7.003: 2.50000 95.1) Finally, Eqs. (E9.5.10) and (E9.5.11) can be substituted into Eq. (E9.5.4) 3x1 — 0.1(—2.50000) ~ 0.2(7.00003) = 7.85 which can be solved for 785+ 0.1 50000) + 0.2(7.00003) 3 Although there isa slight round-off error in Bq. (E9.5.10), the results are very close tothe exact solution of x; = 3, x2 = ~2.5, and x5 = 7. This can be verified by substituting the results into the original equation set 3(3) — 0.1(—2.5) — 0.2(7,00003) = 7.84999 = 7.85 0.1(3) + 7(-2.5) ~ 0.3(7.00003) = ~19.3000 = -19.3, 0.3(3) ~ 0.2(—2.5) + 10(7.00003) = 71.4003 = 71.4 = 3.00000 9.2.1 Operation Counting ‘The execution time of Gauss elimination depends on the amount of floating-point opett tions (or FLOPS) involved in the algorithm. In general, the time consumed to perfor multiplications and divisions is about the same, and is larger than for additions and su tractions. Before analyzing naive Gauss elimination, we will fist define some quantities thi Ferme Erosw-Sro+Feo — iw Voiete tte tiem Wisis2esttm=™ where O(m") means “terms of order m* and lower” ‘Now letus examine the naive Gauss elimination algorithm in detail. As in Fig. will ist count the multiplication/divsion FLOPS inthe elimination stage. On the fis ee ear 9.40.8 ps9s. 9.2_NAIVE GAUSS EUMINATION 243 {frovah the outer loop, k= 1. Therefor, the limits on the middle loop are from I= 2 ton According to Bq. (9.184), this means thatthe numberof iterations ofthe middle oop will be yal Now for every one of these iterations, there is | division to define factor = terior loop then performs a single multiplication (factor «ay, ton. Finally, there is one addtional multipication ofthe righ ‘Thus, for every iteration of the middle loop, the number of T+ 244 @.19) 444/24, The in- .) for each iteration from j = 2 tshand-side value (factor » by), multiplications is +n (920) ZB total for the first passthrough the outer loop is therefore obtained by multiplying Eq, (9.19) by (9.20) to give fn — 1](1 +n) ‘A similar procedure can be used to estimate the multiply/divide FLOPS for the subse ‘quent iterations of the outer loop. These can be summatized as Outer Loop ‘Midalle Loop k i Flops 1 Qa fe +) 2 30 (fo 2H0) Therefore the total FLOPs for elimination canbe computed as Sten 42-6) = Finn + 2)— ken 42) 4) oa Applying some of the relationships from By. (9.18) yes 3 Thus, the total number of multiply/divide FLOPS is equal to n°/3 plus an additional com- Ponent proportional to terms of order n? and lower. The result is written in this way because asm ges large, the O(n?) terms become negligible, We are therefore justified in concluding that for large n, the effort involved in forward elimination is n°/3. Because only a single loop is used, back substitution is much simpler to evaluate, The ‘number of multiplication FLOPS can be directly taken from Eq. (9.18e), (+ 002) ~ + OM) + {= + oun} =F +00 om) (n+)9.3 _GAUSS EUMINATION, TABLE 9.1 Number of FLOPs for naive Gauss elimination, Back Total Substitution 10 55 100 5050 3 3 000 500500 7 9 “Thus, the total effort in naive Gauss elimination can be represented as Be oats oc) seme. ™ + O(n 7 5 TOW) +S + 00 5 Foe) 023) Forward Back elimination substiution ‘Two useful general conclusions can be drawn from this analy 1. As the system gets larger, the computation time increases greatly. As in Table 9.1, the ‘amount of FLOPS increases nearly 3 orders of magnitude for every order of magnitude ‘nerease in the dimension, Most of the effort is incurred in the elimination step. Thus, efforts to make the method ‘mote efficient should probably focus on this step. ‘Throughout the remainder of this part, we will make operation counts to compare a- temative solution methods, Although we may not go info the detail of the above analysis, the same general approach will be employed. PITFALLS OF ELIMINATION METHODS ‘Whereas there are many systems of equations that can be solved with naive Gauss elimi nation, there are some pitfalls that must be explored before writing a general computer Pro ‘gram to implement the method, Although the following material relates directly to nai ‘Gauss elimination, the information is relevant for other elimination techniques as well 9.3.1 Division by Zero ‘The primary reason thatthe foregoing technique is called “nave” is that daring bot elimination and the back-substitution phases, itis possible that a division by zero ca foe, For example, i we use naive Gauss elimination to solve g 8 Dep +323 4m +6 +70 dn bm tO = the normalization of the first row would involve division by ax = 0. Probler arise when a coefficient is very close to zero, The technique of pivoring has been devel to partially avoid these problems. It will be described in Sec. 9.4.2.ak EXAMPLE 9.6 9.3. PITFALLS OF EUMINATION METHODS 9.3.2 Round-Off Errors -Bven though the solution in Example 9.5 was close co the true answer, there was 8 slight “isoepaney in th result for (Eg. (E9.5.10). This discrepancy, which amounted w oer or of --0.00043 percent, was due to our use of six significant figures during the Computation. IF we had used more significant figures, the error inthe results would be sefunher. If we had used fractions instead of decimals (and consequently avoided ound-off altogether), the answers would have been exact. However, because computers Cary only a Fimited number of significant figures (recall Sec. 3.4.1), round-off errors can ‘occur and must be considered when evaluating the resus “The problem of round-off error can become particulary important when large NUE bers of equations are fo be solved. This is due to the fact that every result s dependent 09 previous results. Consequently, an error inthe carly steps wil tend to propagate —that si ‘will eause errors in subsequent steps. Specifying the system size where round-off err becomes significant is complicated by the fact thatthe type of computer and the properties of the equations are determining factors. A rough rule of thamb is that round-off error may be important when dealing with 100 or more equations. In any event, you should always substitute your answers back into the origina equations to check whether a substantial error has ocurred. However, as dis- assed below, the magnitudes ofthe coefficients themselves can influence whether such an error check ensures a reliable result. 9.3.3 Il-Conditioned Systems “The adequacy ofthe solution depends on the condition ofthe system. In Sec. 9.1.1, a graphy jcal depiction of system condition was developed. As discussed in Sec. 4.2.3, wel: ‘Conditioned systems are those where a small change in one or more of the coefficients results fn a similar small change in the solution, U-conditioned systems are those where small Changes in coefficients result in large changes inthe solution. An altemative interpretation Of illsconcitioning is that a wide range of answers can approximately satisfy the equations. Because round-off errors can induce small changes in the coefficients, these artificial changes can lead to large solution errors for il-conditioned systems, as illustrated in the following example. IIL Conditioned Systems Problem Statement. Solve the following system: x $2x =10 96.1) Lx; +21 = 104 (£9.62) ‘Then, solve it again, but with the coefficient of xy in the second equation modified slightly to 105 Solution. Using Eqs. 9.10) and (9.11), the solution is (10) — 2110.4) 2) — 21.1)GAUSS ELMINATION 1004) ~ 1.100) 1@) =20.) n However, with the slight change of the coefficient a2) from 1.1 to 1.05, the result is changed dramatically to 1@) — 201.05) Notice that the primary reason for the discrepancy between the two results is that the denominator represents the difference of two almost-equal numbers. As illustrated previ ously in See. 3.4.2, such differences are highly sensitive to slight variations in the numbers being manipulated. At this point, you might suggest that substitution of the results into the original equa- tions would alert you to the problem. Unfortunately, for ill-conditioned systems this is often not the case, Substitution of the erroneous values of x Tinto Eqs. (29.6.1) and (E9.6.2) yields 842(1) = 10= 10 1.18) +2) = 108 = 10.4 ‘Therefore, although xj = 8 and x) = | is not the true solution to the original problem, the ‘error check is close enough to possibly mislead you into believing that your solutions are adequate. ‘As was done previously in the section on graphical methods, a visual representative of ill-conditioning can be developed by plotting Eqs. (E9.6.1) and (E9.6.2) (recall Fig. 92). Because the slopes of the ines are almost equal, itis visually difficult to see exactly where they intersect. This visual difficulty is reflected quantitatively in the nebulous results of Example 9.6. We can mathematically characterize tis situation by writing the two equ tions in general form: 29 035) aya + ay Dividing Bq. (9.24) by ap and Eq, (9 that are in the format of straight lines (x ) by aye and rearranging yields alternative versions slope) x1 + intercept}thatthe “d previ- aumbers val equ. sis often 9.6.1) stem, the tions ae atative of Fig. 92) ly whee’ -esuls of wo equa EXAMPLE 9.7 9.3 PITFALLS OF EUMINATION METHODS 207 or, cross-multiplying, ‘which can be also expressed as and — A120) =O 026) Now, recalling that a:ja22 ~ ayaa is the determinant of a two-dimensional system (Eq. (0.3), we arrive at the general conclusion that an il-conditioned system is one with a determinant close to zero. In fat, if the determinant is exactly zero, the two slopes are iden- tical, which connotes either no solution or an infinite number of solutions, as is the case for the singular systems depicted in Fig. 9.2a and b. It is difficult to Specify how close to zero the determinant must be to indicate ill- conditioning. This is complicated by the fact thatthe determinant can be changed by mul- tiplying one or more of the equations by a scale factor without changing the solution. Con sequently, the determinant is a relative value that is influenced by the magnitude of the coefficients, Effect of Scale on the Determinant Problem Statement. Evaluate the determinant of the following systems: (a) From Example 9.1 3x +2 = 18 E97.) =n +2 97.2) (b) From Example 9.6: xy 42x = 10 (e973) Lx +2x = 104 e974) (©) Repeat (b) but with the equations multiplied by 10 Solution. (a) The determinant of Eqs. (E9.7.1) and (E9.7.2), which are well-conditioned, is D -2A-)=8 36 (b) The determinant of Eqs, (E9.7.3) and (E9.7.4), which are ill-conditioned, is D=12)-20.) 02 (©) The results of (a) and (b) seem to bear out the contention that ill-conditioned systems have near-zero determinants. However, suppose that the ill-conditioned system in (0) is multiplied by 10 to give 10x, + 20x, = 100 Lay + 20% = 104 ‘The multiplication of an equation by a constant has no effect on its solution. In ‘addition, i is stil il-conditioned, This can be verified by the fact that multiplying by aEXAMPLE 9.8 PN GAUSS ELIMINATION, Constant has no effect on the graphical solution. However, the determinant is drama. ically affected D = 100) OC) = 26 Not only hts it been raised two orders of magnitude, but itis now over twice as large as the determinant of the well-conditioned system in (a). As illustrated by the previous example, the magnitude of the coefficients interjects 2 scale effect that complicates the relationship between system condition and determinant size. One way to partially circumvent this difficulty is to scale the equations so that the 4 maximum element in any row is equal to 1 Sealing Problem Statement. Scale the systems of equations in Example 9.7 to a maximum vale of | and recompute their determinants Solution, (a) For the well-conditioned system, scaling results in 21 + 0.6672 =6 -05y+ n=l for which the determinant is D=1(1) ~0.667(-0.5) = 1,333 (b) For the ill-conditioned system, scaling gives 0.5m +12 058m +32 for which the determinant is D =05(1) ~ 160.55) = —0.05 (©) For the last case, scaling changes the system to the same form as in (b) and she ‘determinant is also ~0.05. Thus, the scale effect is removed. In a previous section (Sec. 9.1.2), we suggested that the determinant is difficult (0 compute for more than three sinjultaneous equations. Therefore, it might seem that it oes not provide a practical means for evaluating system condition, However, as describe! it Box 9.1, there isa simple algorithm that results from Gauss elimination that can be used'® evaluate the determinant Aside from the approach used in the previous example, there are a variety of other wa to evaluate system condition. For example, there are alternative methods for normalizing tb? clements (se Stark, 1970). In addition, as described inthe next chapter (Sec. 10.3), hem 4 trix inverse and matrix norms ean be employed to evaluate system condition, Final). 9 simple (but time-consuming) testis to modify the coefficients slightly and repeatBox 9.1 In Sec. 9.1.2, we stated that determinant evaluation by expansion of iors was impractical for large sets of equations. Thus, we con cluded that Cramer's rule would be applicable only to small sys tems. However, as mentioned in Sec. 9.3.3, the determinant has value in assessing system condition, It would, therefore, be useful {have a practical method for computing this quan. Fortunately, Gaus elimination provides a simple way to do thi ‘The method is based on the fact that the determinant of a trian. gular matix can be simply computed as the produetof its diagonal 9.3 PITFALLS OF EUMINATION METHODS Determinant Evaluation Using Gauss Elimination Recall thatthe forward-elimination step of Gauss elimination results in an upper triangular system, Because the value of the de- terminant is not changed by the forwatd-eimination process, the determinant can be simply evelusted atthe end of this stp via ce 39.12) D= ayaa, where the superscrips signify the number of times that the ele- ‘ents have been modified by the elimination process, Thus, we can clement: DH anand *+ dye Jan ae ais] 0 am ay 0 0 ayl an) [0 0 ay] ‘The valty ofthis formulation canbe ilustated for a3 by 3 system: where the determinant can be evaluated as (recall Eg (9.4)) ox by evaluating the minors (hats, the 2 by 2 determinants), capitalize on the effort tha has already been expended in reducing the system to tangular form and, in the bargin, come up with a simple estimate of the determinant, There is a slight modification tothe above approach when the program employs partial pivoting (Sec. 9.4.2) For ths ease, the de ‘erminant changes sign every time a row is pivoted. One way to represent this is to modify Eq. B9.1.2 oo Da ay May 013) Where p represents the number of times that rows are pivoted ‘This modification can be incorporated simply into a program, merely keep track ofthe numberof pivots tha take place during the course ofthe computation and then use Eq, (B9.1.3) to evaluate the eterminant +a] 2 D = an an0s5 ~ a7 (0) +419(0) = ayanars solution. If such modifications lead to drastically different results, the system is likely to be ill-conditioned, ‘As you might gather from the foregoing discussion, ill-conditioned systems are prob- Jematic. Fortunately, most linear algebraic equations derived from engineering-problem settings are naturally well-conditioned. In addition, some of the techniques outlined in Sec. 9-4 help wo alleviate the problem, 9.3.4 Singular Systems In the previous section, we learned that one way in which a system of equations can be ill- conditioned is when two or more ofthe equations are nearly identical. Obviously, itis even ‘worse when the two are identical. In such cases, we would lose one degree of freedom, and ‘would be dealing with the impossible case of n ~ 1 equations with m unknowns. Such cases might not be obvious to you, particularly when dealing with large equation sets. Conse- ‘quently, it would be nice to have some way of automatically detecting singularity ‘The answer to this problem is neatly offered by the fact thatthe determinant of a sin- gular system is zero. This idea can, in turn, be connected to Gauss elimination by recog- nizing that aftr the elimination step, the determinant can be evaluated as the product ofthe diagonal elements (recall Box 9.1). Thus, a computer algorithm can test to discern whether 1 zero diagonal element is created during the elimination stage. If one is discovered, the calculation can be immediately terminated and a message displayed alerting the user. We9.4 EXAMPLE 9.9 GAUSS ELIMINATION. will show the details of how this is done when we present a full algorithm for Gauss elin- ination later in this chapter. TECHNIQUES FOR IMPROVING SOLUTIONS. ‘The following techniques can be incorporated into the naive Gauss elimination algorithm to.circumvent some of the pitfalls discussed in the previous section. 9.4.1 Use of More Significant Figures ‘The simplest remedy for ill-conditioning is to use more significant figures in the compute- tion. If your application can be extended to handle Iarger word size, such @ feature will greatly reduce the problem. However, a price must be paid in the form of the computationel and memory overhead connected with using extended precision (recall See. 3.4.1) 9.4.2 Pivoting AAs mentioned at the beginning of Sec. 9.3, obvious problems occur when a pivot element is zero because the normalization step leads to division by zero. Problems may also aise ‘when the pivot element is close to, rather than exactly equal to, zero because if the magni- tude of the pivot element is small compared to the other elements, then round-off errors can be introduced, ‘Therefore, before each row is normalized, itis advantageous to determine the largest available coefficient inthe column below the pivot element. The rows can then be switched 80 that the largest element is the pivot element, This is called partial pivoting. if columns as well as rows are searched for the largest element and then switched, the procedure is called complete pivoting. Complete pivoting is rarely used because switching columns changes the order of the x's and, consequently, adds significant and usually unjustified complexity to the computer program. The following example illustrates the advantages of partial pivoting. Aside from avoiding division by zero, pivoting also minimizes round-off error. As such, it also serves as a partial remedy for ill-conditioning Partial Pivoting Problem Statement. Use Gauss elimination to solve 0.003%; + 3.000%, 1.0000; + 1.0000x3 2,001 1.0000 Note that in this form the first pivot element, a1; = 0.0003, is very close to zero. Then repeat the computation, but partial pivot by reversing the order of the equations. The exact solution is x; = 1/3 and x9 = 2/3. Solution, Multiplying the first equation by 1/(0.0003) yields, x1-+ 10,000% = 6667 ‘which can be used to eliminate 2; from the second equation: 999922 = ~6666puta will tional cement 2 arise naghi largest ritebed >tumns dure is slurs astified ages of and-off 9.4_TECHNIQUES FOR IMPROVING SOLUTIONS 251 ‘which can be solved for ‘This result can be substituted back into the first equation to evaluate x) 2.0001 - 32/3) 0c 99.1) However, due to subtractive cancellation, the result is very sensitive to the number of sig nificant figures carried in the computation: ‘Absolute Value ‘of Percent Relative Error 2 x for xi 0067, 3.23 1099 0.6667 100 esse? 10 | 0.066667 1 0.666667 0.233000 a [Note how the solution for x is highly dependent on the number of significant figures. This is because in Eq. (E9.9.1), we are subtracting two almost-equal numbers. On the other hhand, if the equations are solved in reverse order, the row with the larger pivot element is normalized. The equations are | Elimination and substitution yield x; = 2/3. For different numbers of significant figures, xy ccan be computed from the first equation, asin (@99.2) Significant Figures ‘Thus, a pivot strategy is much more satisfactory.JF (dummy > big) FIGURE 9.5 Pseudocode to implement par ‘cl pi EXAMPLE 9.10 (GAUSS EUMINATION General-purpose computer programs must include a pivot strategy. Figure 9.5 pro- Vides a simple algorithm to implement such a strategy. Notice that the algorithm consists cof two major loops. After storing the current pivot element and its row number as the Vari- ables, big and p, the first loop compares the pivot element with the elements below it to check whether any of these is larger than the pivot element. If so, the new largest element and its row number are stored in big and p. Then, the second loop switches the original pivot row with the one with the largest element so that the latter becomes the new pivot row. This pseudocode can be integrated into a program based on the other elements of Gauss elimination outlined in Fig. 9.4. The best way to do this is to employ a modular ap: proach and write Fig. 9.5 as a subroutine (or procedure) that would be called directly after the beginning ofthe first loop in Fig. 94a Note that the second IF/THEN construct in Fig. 9.5 physically interchanges the rows, For large matrices, this can become quite time consuming. Consequently, most codes do rot actually exchange rows but rather keep track of the pivot rows by storing the appropri ‘ate subscripts in a vector. This vector then provides a basis for specifying the proper row ordering during the forward-elimination and back-substitution operations. Thus, the oper- ations are said to be implemented in place. 9.4.3 Scaling In Sec. 9.3.3, we proposed that scaling had value in standardizing the size of the deter ‘inant. Beyond this application, it has utility in minimizing round-off errors for those ‘cases where some of the equations in a system have much larger coefficients than others. Such situations are frequently encountered in engineering practice when widely different units are used in the development of simultaneous equations, For instance, in elect circuit problems, the unknown voltages can be expressed in units ranging from microvols to kilovolts. Similar examples can arise in all fields of engineering. As long as each ‘equation is consistent, the system will be technically correct and solvable. However, the use of widely differing units can lead to coefficients of widely differing magnitudes. This, Jn tum, ean have an impact on round-off esror as it affects pivoting, as iMlustrated by the following example. Effect of Scaling on Pivoting and Round-Off Problem Statement. (4) Solve the following set of equations using Gauss elimination and a pivoting sated 2xy + 100,000% = 100,000 nt m=? (b) Repeat the solution after scaling the equations so that the maximum coefficient in each (©) Fly, use the scaled cooficins 10 determine wheter pivoing is nese However, actually solve the equations withthe orginal coefficient values For ceases, retain only three significant figures. Note that the correct answers are 1: 1,00002 and x; = 0.99998 or, for three significant figures, xj = x9 = 1.00. %9.4 TECHNIQUES FOR IMPROVING SOLUTIONS Solution. (@) Without scaling, forward elimination is applied to give 100,000 50,000 Although x2 is correct, x; is 100 percent in eror because of round-off. (b) Scaling transforms the original equations to Forward elimination yields mit =2 22 = 1.00 which can be solved for ‘Thus, scaling leads to the correct answer, (©) The scaled coefficients indicate that pivoting is necessary, We therefore pivot but retain the original coefficients to give at m=2 2x, + 100,000x2 = 100,000 Forwatd climination yields at n= 100,000. (00,000 Which can be solved for the correct answer: x; = 12 = 1. Thus, scaling was useful in determining whether pivoting was necessary, but the equations themselves did not require scaling to arrive at a correct result.(GAUSS ELIMINATION. | DIMENSION = (a) 5 = ABS(22) ms=20 4 TF ABS( ay s)>5) THEN 5; = BSC i 0 09 a exo 09 i cate er) i IFer # -1 7H inate (2, k= Ln-1 )-< to) THEN 1 FIGURE 9.6 | Peeudocode to inlement Gaus elimination wih por photng. UB Pivot (2, b pek big = ABS(a4n/5e) OHek+La uminy = ABS( Arya) 9 THEN IF dummy > big = dunay P &N 00 IF pk THEY acy = dummy ew 00