Chapter 2Linear EquationsOne of the problems encountered most frequently in scientific computation is thesolution of systems of simultaneous linear equations. This chapter covers the solu-tion of linear systems by Gaussian elimination and the sensitivity of the solution toerrors in the data and roundoff errors in the computation.2.1 Solving Linear SystemsWith matrix notation, a system of simultaneous linear equations is writtenAx = b.In the most frequent case, when there are as many equations as unknowns, A is agiven square matrix of order n, b is a given column vector of n components, and xis an unknown column vector of n components.Students of linear algebra learn that the solution to Ax = b can be writtenx = A−1b, where A−1is the inverse of A. However, in the vast majority of practicalcomputational problems, it is unnecessary and inadvisable to actually computeA−1. As an extreme but illustrative example, consider a system consisting of justone equation, such as7x = 21.The best way to solve such a system is by division:x =217= 3.Use of the matrix inverse would lead tox = 7−1× 21 = 0.142857 × 21 = 2.99997.The inverse requires more arithmetic—a division and a multiplication instead ofjust a division—and produces a less accurate answer. Similar considerations applyDecember 26, 200512 Chapter 2. Linear Equationsto systems of more than one equation. This is even true in the common situationwhere there are several systems of equations with the same matrix A but differentright-hand sides b. Consequently, we shall concentrate on the direct solution ofsystems of equations rather than the computation of the inverse.2.2 The MATLAB Backslash OperatorTo emphasize the distinction between solving linear equations and computing in-verses, Matlab has introduced nonstandard notation using backward slash andforward slash operators, “\” and “/”.If A is a matrix of any size and shape and B is a matrix with as many rowsas A, then the solution to the system of simultaneous equationsAX = Bis denoted byX = A\B.Think of this as dividing both sides of the equation by the coefficient matrix A.Because matrix multiplication is not commutative and A occurs on the left in theoriginal equation, this is left division.Similarly, the solution to a system with A on the right and B with as manycolumns as A,XA = B,is obtained by right division,X = B/A.This notation applies even if A is not square, so that the number of equations is notthe same as the number of unknowns. However, in this chapter, we limit ourselvesto systems with square coefficient matrices.2.3 A 3-by-3 ExampleTo illustrate the general linear equation solution algorithm, consider an example oforder three:10 −7 0−3 2 65 −1 5x1x2x3=746.This, of course, represents the three simultaneous equations10x1− 7x2= 7,−3x1+ 2x2+ 6x3= 4,5x1− x2+ 5x3= 6.The first step of the solution algorithm uses the first equation to eliminate x1fromthe other equations. This is accomplished by adding 0.3 times the first equation2.3. A 3-by-3 Example 3to the second equation and subtracting 0.5 times the first equation from the thirdequation. The coefficient 10 of x1in the first equation is called the first pivotand the quantities −0.3 and 0.5, obtained by dividing the coefficients of x1in theother equations by the pivot, are called the multipliers. The first step changes theequations to10 −7 00 −0.1 60 2.5 5x1x2x3=76.12.5.The second step might use the second equation to eliminate x2from the thirdequation. However, the second pivot, which is the coefficient of x2in the secondequation, would be −0.1, which is smaller than the other coefficients. Consequently,the last two equations are interchanged. This is called pivoting. It is not actuallynecessary in this example because there are no roundoff errors, but it is crucial ingeneral:10 −7 00 2.5 50 −0.1 6x1x2x3=72.56.1.Now the second pivot is 2.5 and the second equation can be used to eliminate x2fromthe third equation. This is accomplished by adding 0.04 times the second equationto the third equation. (What would the multiplier have been if the equations hadnot been interchanged?)10 −7 00 2.5 50 0 6.2x1x2x3=72.56.2.The last equation is now6.2x3= 6.2.This can be solved to give x3= 1. This value is substituted into the second equation:2.5x2+ (5)(1) = 2.5.Hence x2= −1. Finally, the values of x2and x3are substituted into the firstequation:10x1+ (−7)(−1) = 7.Hence x1= 0. The solution isx =0−11.This solution can be easily checked using the original equations:10 −7 0−3 2 65 −1 50−11=746.4 Chapter 2. Linear EquationsThe entire algorithm can be compactly expressed in matrix notation. For thisexample, letL =1 0 00.5 1 0−0.3 −0.04 1, U =10 −7 00 2.5 50 0 6.2, P =1 0 00 0 10 1 0.The matrix L contains the multipliers used during the elimination, the matrix Uis the final coefficient matrix, and the matrix P describes the pivoting. With thesethree matrices, we haveLU = P A.In other words, the original coefficient matrix can be expressed in terms of productsinvolving matrices with simpler structure.2.4 Permutation and Triangular MatricesA permutation matrix is an identity matrix with the rows and columns interchanged.It has exactly one 1 in each row and column; all the other elements are 0. Forexample,P =0 0 0 11 0 0 00 0 1 00 1 0 0.Multiplying a matrix A on the left by a permutation matrix to give P A permutesthe rows of A. Multiplying on the right, AP , permutes the columns of A.Matlab can also use a permutation vector as a row or column index to re-arrange the rows or columns of a matrix. Continuing with the P above, let p bethe vectorp = [4 1 3 2]Then P*A and A(p,:) are equal. The resulting matrix has the fourth row of A asits first row, the first row of A as its second row, and so on. Similarly, A*P andA(:,p) both produce the same permutation of the columns of A. The P*A notationis closer to traditional mathematics, P A, while the A(p,:) notation is faster anduses less memory.Linear equations involving permutation matrices are trivial to solve. Thesolution toP x = bis simply a rearrangement of the components of b:x = PTb.An upper triangular matrix has all its nonzero elements above or on the maindiagonal. A unit lower triangular