gaussian elimination

G1BINM Introduction to Numerical Methods 6–1 6 Gaussian Elimination 6.1 Simultaneous linear equations Consider the system of linear equations a11x1+a12x2+· · · +a1nxn=b1 a21x1+a22x2+· · · +a2nxn=b2 ............ an1x1+an2x2+· · · +annxn=bn  (6.1) which is to be solved for x1, . . . , x n. Notice that the number of equations, n, is equal to the number of unknowns, so it is reasonable to expect that there should usually be a unique solution. Such systems arise in many different applications, and it is important to have an efficient method of solution, especially when nis large. The system (6.1) can be written more succinctly, using the summation convention, as aijxj=bi, or alternatively, using matrix notation, as Ax=b, (6.2) where Ais an n×nmatrix with entries {aij},x= (x1, . . . , x n)Tandb= (b1, . . . , b n)T. IfAis invertible, we can formally multiply (6.2) through by A−1to obtain x=A−1b. Example 6.1 The system −x1+ x2+ 2 x3= 1 , 3x1− x2+ x3= 1 , −x1+ 3 x2+ 4 x3= 1 ,(6.3) may be written as Ax=b, where A= −1 1 2 3−1 1 −1 3 4 ,x= x1 x2 x3 ,b= 1 1 1 . The inverse of this Awas found in §2. The solution may thus be determined via x=A−1b=1 10 −7 2 3 −13−2 7 8 2 −2  1 1 1 =1 10 −2 −8 8 , which gives x1=−1/5,x2=−4/5,x3= 4/5. 6–2 School of Mathematical Sciences University of Nottingham The problem with this approach is that, when nis large, it is extremely time- consuming to calculate A−1using determinants. In this section we discuss the method of Gaussian elimination, which provides a much more efficient algorithm for solving systems like (6.1). 6.2 Doing it by hand In practice, one would go about solving a system like (6.3) by eliminating the variables one at a time until just one remains. Then the other variables would be determined by back-substitution. Gaussian elimination is a formal procedure for doing this, which we illustrate with an example. Example 6.2 Consider again the system (6.3). We eliminate the variables one at a time as follows. 1. Eliminate x1from the second and third equations by subtracting suitable multiples of the first equation ( −3and1respectively). This results in the new system −x1+ x2+ 2 x3= 1 , 2x2+ 7 x3= 4 , 2x2+ 2 x3= 0 .(6.4) 2. Subtract a suitable multiple (here 1) of the second equation from the third to eliminate x2: −x1+ x2+ 2 x3= 1 , 2x2+ 7 x3= 4 , −5x3=−4.(6.5) 3. Now we can solve by back-substitution . The third equation gives x3= 4/5; then the second gives 2x2= 4−7×4/5⇒x2=−4/5. Finally, from the first equation, −x1= 1−(−4/5)−2×4/5⇒x1=−1/5. The elimination process just described is performed by applying so-called elementary operations to the equations. These comprise 1. swapping two equations; 2. multiplying an equation by a nonzero constant; 3. adding a multiple of one equation to another equation. The crucial point is that, when any combination of these operations is performed upon a system of equations, the result is an equivalent system, in the sense that it has the same set of solutions. Thus the systems (6.3), (6.4) and (6.5) are all equivalent. The idea of Gaussian elimination is to use elementary operations to reduce the original system to an equivalent one which is in “triangular” form, and can then readily be solved by back-substitution. G1BINM Introduction to Numerical Methods 6–3 6.3 The augmented matrix; elementary row operations In practice, when carrying out this procedure on a general matrix equation of the form (6.2), it is useful to define the so-called augmented matrix /tildewideA, which consists of Awith the right-hand side btacked on as an extra column: /tildewideA= (Ab). Thus each n-dimensional system of the form (6.1) corresponds to an n×(n+1) augmented matrix, and to each elementary operation, which may be performed upon a system of equations, corresponds an elementary row operation , which may be performed upon the augmented matrix /tildewideA: 1. swapping two rows; 2. multiplying a row by a nonzero constant; 3. adding a multiple of one row to another row. In Gaussian elimination we perform a combination of these operations such as to reduce the augmented matrix to a triangular form, known as echelon form . Example 6.3 Consider the system 3x1−4x2+ 5 x3=−1, −3x1+ 2 x2+ x3= 1 , 6x1+ 8 x2− x3= 35 , whose corresponding augmented matrix is  3−4 5 −1 −3 2 1 1 6 8 −1 35 . The first operation is to use the first equation to eliminate x1from the other two. This corre- sponds to subtracting appropriate multiples of row one from rows two and three so as to eliminate their first column. The boxed number 3is called the pivot : it must be divided into the first entry of each subsequent row to determine the appropriate multiple of row one to be subtracted. Here −3/3 =−1and6/3 = 2 , so we subtract (−1)×(row one) from (row two) and 2×(row one) from (row three), resulting in 3−4 5 −1 0−2 6 0 0 16 −11 37 . Next we use row two to eliminate the second element from row three. Here the pivot is −2; 16/(−2) =−8so we subtract (−8)×(row two) from (row three) two obtain  3−4 5 −1 0−2 6 0 0 0 37 37 . This is now in echelon form which completes the Gaussian elimination. It remains to perform back substitution [exercise ]which gives the solution x3= 1,x2= 3,x1= 2. 6–4 School of Mathematical Sciences University of Nottingham Let us reiterate the process just described: 1. use multiples of row 1 to eliminate column 1 from rows 2, 3, . . . ; 2. use multiples of row 2 to eliminate column 2 from rows 3, . . . ; 3. . . . and so on. Note once more the role of the pivot in determining the right multiple to be subtracted at each stage. In particular, the pivot must be nonzero for the process to work. 6.4 Pivoting Example 6.4 The system x1+ x2+ x3=−1, 2x1+ 2 x2+ 5 x3=−8, 4x1+ 6 x2+ 8 x3=−14, corresponds to the augmented matrix  11 1 −1 2 2 5 −8 4 6 8 −14 . First use the pivot 1to clear out column one:  1 1 1 −1 003−6 0 2 4 −10 . Now we appear to be in trouble, since 0cannot be used as a pivot: it does not allow us to clear out the row below. The remedy to this difficulty is simply to switch rows 2and3:  1 1 1 −1 024−10 0 0 3 −6 . Now the pivot is nonzero, and indeed the matrix is in echelon form and ready for back-substitution. The process of swapping two rows as just described is called pivoting . This suggests the following modification of the procedure outlined above: 1. check that the first pivot is nonzero; (a) if it iszero, find a row beneath whose first entry is nonzero and swap it with row 1; 2. use multiples of row 1 to eliminate column 1 from rows 2, 3, . . . ; 3. check that the second pivot is nonzero; G1BINM Introduction to Numerical Methods 6–5 (a) if it iszero, find a row beneath whose second entry is nonzero and swap it with row 2; 4. use multiples of row 2 to eliminate column 2 from rows 3, . . . ; 5. . . . and so on. 6.5 When pivoting fails Notice that the procedure outlined above fails if we encounter a zero entry in the pivot position, but there is nothing nonzero with which to swap it. Example 6.5 The system x1+ x2+ x3=−1, 2x1+ 2 x2+ 5 x3=−8, 4x1+ 4 x2+ 8 x3=−14,(6.6) corresponds to the augmented matrix  11 1 −1 2 2 5 −8 4 4 8 −14 . As before, begin by using the first pivot, 1, to eliminate the first column:  1 1 1 −1 003−6 0 0 4 −10 . (6.7) Now the entry in the pivot position is zero, and the situation cannot be remedied by pivoting. In fact the system (6.6) has no solution ; the final two equations in (6.7) read 3 x3=−6 and 4 x3=−10, and are therefore inconsistent. In general, if no nonzero pivot can be found, this corresponds to the original matrix Abeing singular, i.e.|A|= 0. As an aside, we should point out that Abeing singular does not necessarily imply thatAx=bhas no solution. Example 6.6 If we adjust the right-hand side of (6.6) by changing the −14to−12, the aug- mented matrix  1 1 1 −1 2 2 5 −8 4 4 8 −12  is reduced by Gaussian elimination to  1 1 1 −1 0 0 3 −6 0 0 4 −8 . Now the final two equations both give x3=−2and the first implies x1+x2= 1, i.e. there is an infinity of solutions. 6–6 School of Mathematical Sciences University of Nottingham This situation, though, is nongeneric: any other value apart from −12 in this last entry would have given no solution. Therefore, when performing Gaussian elimination in practice (using real arithmetic with finite precision) it is safest simply to stop the calculation if Aappears to be singular; the system almost certainly has no solution. 6.6 Implementation of Gaussian elimination Now we show how one could write a computer program to carry out Gaussian elimina- tion. We illustrate the algorithm by means of a flowchart in figure 1. It is straightforward to translate this into ( e.g.)Maple . The main loop variable igoes from 1 to n−1; at each stage the first step is to check whether the pivot entry aiiis zero. If it is, we pass to the pivoting subroutine, otherwise the program proceeds directly to the elimination stage. In the pivoting subroutine, j looks at each row in turn below the ithrow until it finds one with a nonzero pivot entry. The program then swaps rows iandjand returns to the elimination stage. If jreaches nand no nonzero pivot entry has been located, then the matrix is singular. The purpose of the elimination stage is to use the ithrow to eliminate the ithcolumn. Thus kcycles through all the rows below the ithone; for each, rcalculates the multiple of row ithat must be subtracted from row k. When kreaches n, elimination of the ithcolumn is completed, and so ican be incremented. When ireaches n, Gaussian elimination is finished, the matrix is in echelon form, and back-substitution may proceed. 6.7 Operations count When evaluating any algorithm, the following criteria should be borne in mind: 1. ease of programming; 2. effect of rounding errors; 3. storage space required; 4. computing time required. Gaussian elimination is straightforward to program, and rounding errors may be con- trolled by the method of partial pivoting , described below. It requires the n×(n+ 1) augmented matrix to be stored, which may be inefficient if nis very large, but most of the coefficients of Aare zero. Such systems are called sparse , and some routines that may be more efficient for large, sparse systems are discussd in §6. Here we address the question of computing time. An estimate of the time it will take to perform a calculation may be obtained by counting the number of arithmetic operations involved. It is standard to count each multiplication and each division as a single operation, but to ignore additions and subtractions (which require significantly less processor time). G1BINM Introduction to Numerical Methods 6–7 START ............................................................................................................................................ .. ...................................................................................................................................................................................................................................................................................... .. . . ...........................................................................................................................................Doi= 1 to n−1 ...................................................................................................... ...................................................................................................... ...................................................................................................... ...................................................................................................... Isaii= 0? ............................................................................................................................................ . . . ....................................................................................................................................................................................................................................................................................... . ............................................................................................................................................Dok=i+ 1 to n Setr=aki/aii Subtract r×(row i) from (row k) ....................................................................................... ........................................................................................................................................................................... . ......................................................................................Next k ....................................................................................... ............................................................................................................................................................................ ......................................................................................Next i ELIMINATION FINISHED............................................................................................................................................ . . . ..................................................................................................................................................................................................................................................................................... .. . ............................................................................................................................................Doj=i+ 1 to n Swap row i with row j....................................................................................... ............................................................................................................................................................................ ......................................................................................Next j MATRIX IS SINGULAR...................................................................................................... ...................................................................................................... ...................................................................................................... ...................................................................................................... Isaji= 0?................... ................... ................... ................... ................... ................... ................... ................... ................... ................... ................... ................... ................... ................... ................... ...................................... ................... ................... ...................................... ................... ................... ...................................... ................... ................... ...................................... ...................................... ................... ................... ................... k < n i < n −1k=n i=n−1Yes No j=nj < n Yes NoPivoting Elimination Figure 1: Flow chart for Gaussian elimination 6–8 School of Mathematical Sciences University of Nottingham The first step in Gaussian elimination is to use the top row to eliminate the first column from each subsequent row. For each, we have to carry out 1 division, to calculate the multiple of row 1 to be substracted ( i.e. the ratio rin figure 1), followed by n multiplication-subtractions. This makes n+ 1 operations for each of the n−1 rows below the top one, leading to a total of ( n−1)(n+ 1) = n2−1 operations needed to eliminate the first column. Next we continue to the second row. Effectively we must perform the same operation as that just described, but with an ( n−1)×nmatrix instead of n×(n+ 1). Thus the number of operations needed to eliminate column 2 is ( n−1)2−1. Therefore, to complete the elimination through all nrows, the total number of operations required is1 Nn=n/summationdisplay k=1/parenleftbig k2−1/parenrightbig =n(n+ 1)(2 n+ 1) 6−n=n(n−1)(2n+ 5) 6. Of particular interest is the behaviour of Nnwhen nis large: Nn=n3 3+n2 2−5n 6∼n3 3 since n3/greatermuchn2ornfor large n. Thus an adequate estimate of the number of operations required for Gaussian elimination of an n-dimensional system is n3/3. This should be compared with the n! operations needed to find the determinant of the same system. It should be clear that Gaussian elimination is vastly preferable to solving using determinants if nis at all large. For example, a 100 ×100 system requires fewer than 106operations, which a 1000MHz machine could perform in less than a millisecond. 6.8 Back-substitution Now suppose the elimination process is complete and the augmented matrix is in echelon form. Thus the reduced system takes the form a11x1+a12x2+· · · + a1(n−1)xn−1 + a1nxn = b1, a22x2+· · · + a2(n−1)xn−1 + a2nxn = b2, ............ a(n−1)(n−1)xn−1+a(n−1)nxn=bn−1, annxn = bn. The next stage is back-substitution, of which the first step is to solve the final equation for xn: xn=−bn/ann. 1Here we use the identity n/summationdisplay k=1k2=n(n+ 1)(2 n+ 1) 6. G1BINM Introduction to Numerical Methods 6–9 Then with xnknown, the penultimate equation gives xn−1=/parenleftbig bn−1−a(n−1)nxn/parenrightbig /a(n−1)(n−1). By continuing in a similar fashion we may generate in turn xn−2, . . . , x 2, x1; the algo- rithm for back substitution may thus be written in the following form. •START •Dok=nto 1 step ( −1) •Setxk=1 akk bk−n/summationdisplay j=k+1akjxj  •End do •END Now let us count the number of operations required for back substitution. To solve forxnrequires just one division, while one multplication-subtraction followed by one division are needed to solve for xn−1. In general, then, at the kthstep, koperations are needed, giving rise to a total of n/summationdisplay k=1k=n(n+ 1) 2∼n2 2 when nis large. The main point is that this is negligible compared to the number of operations used in the elimination stage. 6.9 Partial pivoting We have seen that it may be necessary to swap two rows if a diagonal entry, which we intend to use as a pivot, turns out to be zero. However, when implementing the procedure on a computer with finite-precision arithmetic, testing whether a quantity is equal to zero is unsafe. In any case, inaccuracy may arise if the pivot is too small, so it is prudent at each stage to choose the row with the largest possible entry in the pivot position, and use that as the pivot. Example 6.7 Consider the augmented matrix /parenleftbigg 0.0001 1 1 1 1 2/parenrightbigg , (6.8) in which the pivot entry, although nonzero, is small. The effect of this is that we have to use a large multiple, 10000, in eliminating the first column, resulting in /parenleftbigg 0.0001 1 1 0 −9999 −9998/parenrightbigg . 6–10 School of Mathematical Sciences University of Nottingham Now back-substitution gives x2= 9998 /9999 , and if only three decimal places are stored, this is approximated by x2= 1.000. Then we find 0.0001x1= 0.000, that is x1= 0.000, which is the wrong answer. However if, on the grounds that the pivot entry is too small, we first swap rows: /parenleftbigg 1 1 2 0.0001 1 1/parenrightbigg , then Gaussian elimination leads to /parenleftbigg1 1 2 0 0.9999 0 .9998/parenrightbigg . Now back-substitution results in (to three decimal places) x2= 1.000,x1= 1.000, which is correct. Errors, like those just illustrated, which arise from the fact that a computer can only store real numbers with a finite precision are called rounding errors . Actually, the argument for selecting the largest possible pivot entry is still slightly unsound, since we can multiply a whole row through by any factor we like to make the pivot entry bigger, without alleviating the problem. For example, multiplication of the first row of (6.8) by 10000 results in /parenleftbigg 110000 10000 1 1 2/parenrightbigg . Now the pivot entry in the first row looks fine, but the same loss of accuracy occurs if Gaussian elimination is attempted without pivoting [ exercise ]. What one should really look for, when deciding whether or not to swap rows, is the size of the pivot entry compared with the other entries in the same row. In the method ofpartial pivoting , we choose at each stage the row which maximises (in absolute value) the ratio between the pivot entry and the largest entry in the same row (apart from the last one, i.e.the right-hand side). The row which maximises this so-called pivot ratio is then switched into the pivot position. Example 6.8 [Kreyszig] Consider the augmented matrix  3−4 5 −1 −3 2 1 1 6 8 −1 35 . In row one the largest entry (apart from the last) is 5, giving a pivot ratio 3/5. The corresponding ratios for rows two and three are 3/3and6/8, of which that for row two is the biggest. We therefore swap row two into the pivot position,  −3 2 1 1 3−4 5 −1 6 8 −1 35 , and then eliminate column one:  −3 2 1 1 0−2 6 0 0 12 1 37 . G1BINM Introduction to Numerical Methods 6–11 Now the 6is the largest entry in row two, giving a pivot ratio of 2/6, compared with 12/1for row three. So we switch row three into the pivot position:  −3 2 1 1 0 12 1 37 0−2 6 0 , and eliminate to give −3 2 1 1 0 12 1 37 0 0 37 /6 37 /6 , which is now ready for back-substitution. Note that, as well as avoiding unnecessary rounding errors, this procedure automatically swaps zero pivots where necessary and possible ( i.e.it replaces the pivoting stage in figure 1). Partial pivoting may be written as an algorithm as follows. •START •Doi= 1 to n−1 •Settemp1= 0.0 •Setm= 0 •Doj=iton •Settemp2= max k=i+1,...,n|ajk| •If|aji|>temp1×temp2then ·Settemp1=|aji|/temp2 ·Setm=j •End if •End do •Ifm/negationslash=ithen •Swap row iwith row m •End if •ELIMINATE COLUMN i •End do •ELIMINATION FINISHED This fits inside the main loop of figure 1, at the point where we are about to use rowito eliminate column i. We start by initialising a temporary real variable temp1 and an integer m; these will be used to store the maximum pivot ratio found so far, and the row in which it occurs. Now jcycles through all the rows below and including the ithone, and in each such row, temp2calculates the largest entry. Then the pivot ratio 6–12 School of Mathematical Sciences University of Nottingham for that row is given by |aji|/temp2; we check to see whether this is greater than temp1. If so, row jholds the largest pivot ratio found so far, so temp1and mare replaced accordingly. When jhas run through every row, mtells us which row has the largest pivot ratio; if mis not equal to i, we swap rows iandm. Now we have the best possible row in the pivoting position, and are ready to continue with elimination. Finally, we mention that there is an even more complicated procedure, called total pivoting , which involves swapping columns as well as rows to maximise the pivot ratio at each stage. This is even safer, so far as rounding errors are concerned. However, it is also considerably more expensive in terms of computer time, and since partial pivoting usually gives adequate accuracy, is rarely used.