Linear Systems Solution Methods and Examples

Chapter 2. Solution of Linear Systems AX=B

Example Material Selection Manufacturing a product containing Manganese, silicon, and copper Required amounts [15; 22; 39] lbs/ton manufactured 3 supplies exist with different concentrations of each ingredient Problem: How much to buy from each supplier?

Example - Continued Xj= amount to buy from supplier j Ci= amount of ingredient i required per ton of product manufactured aij= amount of ingredient i in each ton from supplier j Supplier 1 (lb/ton) 1 2 3 Supplier 2 (lb/ton) 3 4 4 Supplier 3 (lb/ton) 2 3 7 Manganese Silicon Copper

Example Cont. Relation among composition from suppliers, the amount needed from each supplier, and the composition of the final product 3 = j = 3 , 2 , 1 = a X C for i ij j i 1 + + = 3 2 15 X X X 1 2 3 + + = 2 4 3 22 X X X 1 2 3 + + = 3 4 7 39 X X X 1 2 3

Example - Truss Truss structure consisting of rigid members interconnected by joints Example B g1 f1 f4 f3 f7 f3 f1 f8 D f4 f6 C A f5 f2 f5 f2 g2

Example Cont. Given values of the load forces, g1and g2, find the values of the forces in the beams, f1, , f8 + = cos 0 f f f 0 cos 1 0 0 0 1 0 0 0 f 1 2 6 1 + = sin 0 f f sin 0 0 0 0 1 0 0 g f 1 7 2 + + = cos cos 0 f f g cos 0 0 cos 0 0 0 0 f 1 f 4 1 3 1 = 0 sin sin + 0 f f sin 0 1 sin 0 0 0 0 0 f 1 3 4 f 4 = = 0 0 f 0 1 0 1 0 0 0 0 g f 2 5 5 = 0 f g 0 1 0 0 0 0 0 f 3 2 6 2 = 0 cos 0 f f 0 0 0 cos 1 0 0 0 0 f 4 5 7 + = sin 0 f f 0 0 0 sin 0 0 1 0 f 4 8 8

Systems of Linear Equations Consider the linear system Ax = b where A is an (n n) matrix, x is the vector of (n) unknown solution values, and b is a column vector of constants a a a x b 11 12 1 1 1 n a a a x b 21 22 2 2 2 n = a a a x b 1 2 n n nn n n + + + = a x a x a x b 11 1 12 2 1 1 n n x + + + = a x a x a b 21 1 22 2 2 2 n n + + + = a x a x a x b 1 1 2 2 n n nn n n

Solution of a Linear System A formal way to obtain a solution using matrix algebra is to multiply each side of the equation by the inverse of A to yield 1 1 = A Ax A b 1 = x A b

Example 2 unknowns in 2 equations 14 + = 3 2 7 x x = 1 2 12 + 4 1 x x 1 2 10 8 Solution 6 3 7 = + y x x 2 1 4 2 2 2 = + 41 x 1 x 2 0 -2 -4 -4 -3 -2 -1 0 1 2 x

Upper-Triangular Linear Systems Def. 2.1. An N N matrix A=[aij] is called upper triangular provided that the elements satisfy aij=0 whenever i>j. A is called lower triangular provided that aij=0 whenever i<j. Thm. 2.1(Back Substitution). Suppose that AX=B is an upper- triangular system with the form as below. If akk 0 for k=1,2, ,N, then there exists a unique solution for the system. + + 3 13 2 12 1 11 + + = a x a x a x a x b 1 1 N N + + + = a x a x a x b 22 2 23 3 2 2 N N + + = a x a x b 33 3 3 3 N N = a x b NN N N

Back Substitution to solve the linear system First, + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N b = N x a + + + = a x a x a x b N 22 2 23 3 2 2 N N N + + = a x a x b 33 3 3 3 N N = a x + a x b , 1 , 1 1 1 1 N N N N N N N = a x b NN N N b - a x , , 1 1 = N N N N x Algorithm: 1 N a , 1 1 N N Finally, = x b /u n n nn n n = = i j = = , ,..., 1 1 , 2 x (b u x )/u i n ( ) b a x 1 1 k k + + + ( ) i i ij j ii b 12 a x a x a x 1 2 13 3 1 2 n n k = = x + 1 1 11 a 11 a

Upper-Triangular Linear Systems The condition that akk 0 is essential because the back-substitu- tion algorithm involves division by akk. If this requirement is not fulfilled, either no solution exists or infinitely many solution exist. Thm. 2.2. If the N N matrix A=[aij] is either upper or lower triangular, then N = = det( ) A 11 22 a a a a NN ii = 1 i

Gaussian Elimination Two linear systems of dimension N N are said to be equivalent provided that their solution sets are the same. When elementary transformations are applied to a given system, the solution sets do not change. A scheme for solving a general system AX=B of N equations and N unknowns is to construct an equivalent upper-triangular system UX=Y that can be solved by the back-substitution algorithm.

Gaussian Elimination A simple example: The back-substitution algorithm is now used to find the unknowns + = 3 2 x 7 x x = 1 2 + 4 1 . x 1 2 The variable x1is eliminated from the second equation by substracting from it 4/3 times the first equation. We arrive at the equivalent system: 2 3 2 1 + x x x = 5 2 and 7 2 5 3 = = 1. x 1 = 7 5 25 = . x 2 3 3

Gaussian Elimination It is efficient to store all the coefficients of the linear system AX=B in an array of dimension N (N+1). The coefficients of B are stored in column N+1 of the array(i.e., ak,N+1=bk). Each row contains all the coefficients necessary to represent an equation in the linear system. The augmented matrix is denoted [A|B] and the linear system is represented as follows: a x a x a + + + + + + = x a x b 11 1 12 2 13 3 1 1 n n x + + = a x a x a x a b 21 1 22 2 23 3 2 2 n n + + + + = a x a x a x a x b 31 1 32 2 33 3 3 3 n n + + + + = . a x a x a x a x b 1 1 2 2 3 3 n n n nn n n The system can be solved by performing row operations on the augmented matrix [A|B]. The variables xk are placeholders for the coefficients and can be omitted until the end of the calculation.

Gaussian Elimination Thm. 2.4(Elementary Row Operations). The following operations applied to the augmented matrix yield an equivalent linear system. Interchanges: The order of two rows can be changed. Scaling: Multiplying a row by a nonzero constant. Replacement: The row can be replaced by the sum of that row and a nonzero multiple of any other row; that is : rowr=rowr-mrp rowp. It is common to use the last operation by replacing a row with the difference of that row and a multiple of another row.

Gauss Elimination Forward Elimination - multiple of one equation is subtracted from another to eliminate unknowns Back Substitution last equation yields unknown, substituting back into the other equations yields the rest

Gaussian Elimination Process + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N x a 21/a + + + + = a x a x a x a b 11 21 1 22 2 23 3 2 2 N N + + + + = a x a x a x a x b 31 1 32 2 33 3 3 3 N N + + + + = a x a x a x a x b 1 1 2 2 3 3 N N N NN N N + + + + + + + + = = = a x a x a x a x a x + a x b a a + + 11 1 0 a x a x a x + 12 2 13 3 1 1 N N x b x b 22 2 23 3 2 2 N N 31 1 32 2 33 3 3 3 N N + x a + a + + = a x a x x b 1 1 2 2 3 3 N N N NN N N

Gaussian Elimination Process-Cont. 31/a a + + + + = a x a x a x a x b 11 11 1 12 2 13 3 1 1 N N x + + + + = a x a x a x a b aN 1/a 21 1 22 2 23 3 2 2 N N 11 + + + + = a x a x a x a x b 31 1 32 2 33 3 3 3 N N + + + + = a x a x a x a x b 1 1 2 2 3 3 N N N NN N N + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N 22 23 2 2 + + + = a x a x a x b 2 3 N N 32 33 3 3 + + + = a x a x a x b 2 3 N N N N NN N + + + = a x a x a x b 2 2 3 3 N

Gaussian Elimination Process-Cont. + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N 22 / ja a 22 23 2 2 + + + = a x a x a x b 2 2 3 N N 32 33 3 3 + + + = a x a x a x b 2 3 N N N N NN N + + + = a x a x a x b 2 2 3 3 N + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N 22 23 2 2 + + + = a x a x a x b 2 3 N N 33 + + = a x a x b 3 3 3 N N = ... NN ... N a x b N

Gaussian Elimination Process-Cont. + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N 22 23 2 2 + + + = a x a x a x b 2 3 N N 33 + + = a x a x b 3 3 3 N N = ... NN ... N a x b N ... N b = x N ... NN a Finally, back-substitution can be used to solve the upper-triangular systems above.

Example Forward Elimination Back Substitution 78 = = 3 x 3 26 28 10 28 10 3 2 + + = 2 4 16 x x x + + = 2 4 + 16 x x x 1 2 3 1 2 x 3 1 + = x - = = = + + 3 2 10 x x = x 5 14 - x x 1 2 3 2 x 3 2 3 2 + + = 3 3 16 x x x + + = - 1 2 3 3 3 16 x x 1 2 3 + 16 ( 4 ) x x = 2 2 3 x + + = + + = 2 4 16 2 4 16 - x x x x x x 1 2 3 1 2 3 1 1 = 10 28 x x = 5 14 - x x 2 3 16 2 4 3 2 1 2 3 2 = = 26 78 x 3 5 + = 8 x x 2 3 2 =

Gauss Elimination Gauss Elimination Input A, b, n Forward Elimination Back Substitution Output Solution Vector, x stop

Gauss Elimination Forward Elimination Loop over all rows in matrix, except last k = 1, 2, ..., n-1 + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N 22 23 2 2 + + + = a x a x a x b 2 3 N N Loop over all rows below the diagonal position (k,k) 33 + + = i = k+1, ..., n a x a x b 3 3 3 N N Compute multiplier for row (i) and column (k) = ... NN ... N a x b mik = aik/akk N Loop over all columns to the right of the diagonal position (k,k) j = k+1, ..., n Update coefficient for row (i) and column (j) aij = aij - mik*akj Update right hand side for row (i) and column (j) bi = bi - mik*bk

Gauss Elimination Back Substitution Compute last unknown xn = bn/ann + + + + = a x a x a x a x b 11 1 12 2 13 3 1 1 N N 22 23 2 2 + + + = a x a x a x b Loop backwards over all rows except last i = n-1, n-2, ..., 1 2 3 N N 33 + + = a x a x b 3 3 3 N N sum = 0 = ... NN ... N a x b N Loop over all columns to the right of the current row j = i+1, ..., n sum = sum + ai,j*xj Compute (i)th unknown xi = (bi - sum)/ai,i

Computational Complexity Count the operations of the first N columns of [A|B]. The outer loop of step p+1 requires N-p=N-(p+1)+1 divisions to compute the multipliers mrp. Inside the loops, but for the first N columns only, a total of (N-p)(N-p) multiplications and the same number of subtractions are required to compute the new row elements ??? carried out for p=1,2, ,N-1. Thus the triangular factorization portion of A=LU requires (?+1). This process is 3 1 N = N N multiplications and divisions ) 1 + = ( )( N p N p 3 1 p and ? 1 (? ?)(? ?) =2?3 3?2+ ? subtractions. 6 ?=1

Example Division by zero: May occur in the forward elimination steps. Consider the set of equations below. Use five significant figures with chopping = 10 7 7 x 099 . 2 x x 10 7 0 7 x 1 x 2 1 . 3 + x + = = 3 x . 2 099 + 6 . 3 901 x 3 6 901 x 1 2 = 3 2 5 1 5 6 x 5 5 6 x 3 1 2 3 At the end of Forward Elimination 10 7 0 7 x 1 = 0 . 0 001 6 . 6 001 x 2 0 0 15005 15004 x 3

Example-Cont. Back Substitution 15004 10 7 0 7 x = = 99993 . 0 x 1 3 15005 = 0 . 0 001 6 . 6 001 x 2 0 0 15005 15004 x . 6 001 6 x 3 = = 5 . 1 3 x 2 . 0 001 + 7 7 0 x x = = 3500 . 0 2 3 x 1 10

Example-Cont. Compare the calculated values with the exact solution . 0 35 x 0 x 1 1 X X = = 99993 5 . 1 x = = 1 1 x calculated 2 exact 2 . 0 x x 3 3 Improvements: Increase the number of significant digits Gaussian Elimination with Pivoting Decreases round off error Avoids division by zero Does not avoid division by zero Reduces round off error

Pivoting Strategy (?)= 0, row p cannot be used to eliminate the elements in column If ??? pbelow the main diagonal. It is necessary to find row k, where ??? 0and k>p, and then interchange row p and row k so that a nonzero pivot element is obtained. If such a nonzero pivot element cannot be found the coefficient matrix is singular so that the system does not have a unique solution. (?) (?) 0, do not switch (?) 0and The trivial pivoting strategy is as follows: If ??? rows. If ??? switch rows k and p. This will result in a nonzero pivot element. (?)= 0, locate the first row below p in which ???

Pitfalls Two Potential Pitfalls - Division by zero: May occur in the forward elimination steps. Consider the set of equations: = x 10 7 7 x x 2 + 3 + = 6 . 2 x 099 + 3 . 3 901 x x 1 2 3 = 5 5 6 x x 1 2 3 - Round-off error: Prone to round-off errors.

Pitfalls: Example Consider the system of equations: Use five significant figures with chopping 099 . 2 10 7 0 7 x 1 3 6 . 3 901 x = 2 5 1 5 6 x 3 At the end of Forward Elimination 7 10 7 0 x 1 . 6 001 0 . 0 001 6 = x 2 15004 0 0 15005 x 3

Pitfalls: Example Back Substitution 15004 10 7 0 7 x = = 99993 . 0 x 1 3 15005 = 0 . 0 001 6 . 6 001 x 2 0 0 15005 15004 x 3 . 6 001 6 x = = 5 . 1 3 x 2 . 0 001 + 7 7 0 x x = = . 0 2 3 3500 x 1 10

Pitfalls: Example Compare the calculated values with the exact solution 0 x 1 X = = 1 1 x exact 2 x 3 . 0 35 x 1 X = = 99993 5 . 1 x calculated 2 . 0 x 3

Improvements Increase the number of significant digits Decreases round off error Does not avoid division by zero Gaussian Elimination with Partial Pivoting Avoids division by zero Reduces round off error

Partial Pivoting Gaussian Elimination with partial pivoting applies row switching to normal Gaussian Elimination. How? At the beginning of the kthstep of forward elimination, find the maximum of a a ,......... , , 1 + ......., a kk k k nk , n k p If the maximum of the values is |apk| In the pthrow, then switch rows p and k.

Partial Pivoting What does it Mean? Gaussian Elimination with Partial Pivoting ensures that each step of Forward Elimination is performed with the pivoting element |akk| having the largest absolute value.

Partial Pivoting: Example Consider the system of equations = 10 x 7 x 7 1 2 + + = 3 x . 2 099 + x 3 x . 3 901 1 2 = 3 5 x x 5 x 6 1 2 3 In matrix form 7 10 7 0 x 1 . 3 901 3 . 2 099 6 x = 2 6 x 5 1 5 3 Solve using Gaussian Elimination with Partial Pivoting using five significant digits with chopping

Partial Pivoting: Example Forward Elimination: Step 1 Examining the values of the first column |10|, |-3|, and |5| or 10, 3, and 5 The largest absolute value is 10, which means, to follow the rules of Partial Pivoting, we switch row1 with row1. Performing Forward Elimination 10 7 0 7 10 7 0 7 x x 1 1 = = 0 . 0 001 6 . 6 001 x 3 . 2 099 6 . 3 901 x 2 2 0 5 . 2 5 5 . 2 x 5 1 5 6 x 3 3

Partial Pivoting: Example Forward Elimination: Step 2 Examining the values of the first column |-0.001| and |2.5| or 0.0001 and 2.5 The largest absolute value is 2.5, so row 2 is switched with row 3 Performing the row swap 5 . 2 10 7 0 7 x 10 7 0 7 x 1 1 = = 0 . 0 001 6 . 6 001 x 0 5 5 . 2 x 2 2 . 6 0 5 . 2 5 5 . 2 x 0 . 0 001 6 001 x 3 3

Partial Pivoting: Example Forward Elimination: Step 2 Performing the Forward Elimination results in: 5 . 2 10 7 0 7 x 1 = 0 5 5 . 2 x 2 . 6 0 0 . 6 002 002 x 3

Partial Pivoting: Example Back Substitution Solving the equations through back substitution . 6 002 5 . 2 = = 10 7 0 7 x 1 x 1 3 . 6 002 = 0 5 5 . 2 x 2 . 6 0 0 . 6 002 002 x 5 . 2 5 x 3 = = 2 1 x 2 5 . 2 + 7 7 0 x x = = 2 3 0 x 1 10

Partial Pivoting: Example Compare the calculated and exact solution The fact that they are equal is coincidence, but it does illustrate the advantage of Partial Pivoting 0 0 x x 1 1 X X = = 1 = = 1 1 x 1 x calculated 2 exact 2 x x 3 3

Triangular Factorization It is easy to solve a problem of the form Lx=b or Ux=b using forward/back elimination. Suppose we are given a matrix A and we are able to find a lower triangular matrix L that has 1 s along the main diagonal and an upper-triangular matrix U with nonzero diagonal elements such that: A = LU Then to solve Ax=b we can solve two simple problems: Ly = b Ux = y (sanity check LUx = L(Ux) = Ly =b)

Triangular Factorization Def. 2.3. The nonsingular matrix A has triangular factorization if it can be expressed as the product of a lower-triangular matrix L and an upper-triangular matrix U: A=LU. In matrix form, this is written as 1 0 1 0 0 1 0 0 0 1 a a a a a a a a a a a a a a a a u u u u u u u u u u 11 12 13 14 11 0 0 0 12 13 14 m m m 21 22 23 24 21 22 0 0 23 24 = m m 31 32 33 34 31 32 33 0 34 m 41 42 43 44 41 42 43 44 The condition that A is nonsingular implies that ukk 0 for all k.

LU Factorization Can all nonsingular matrices A be factored as A=LU? Example: 0 1 1 1 a 0 1 b d c = = A 0 0 The necessary and sufficient condition for the unique LU factorization of an N N matrix A is that the leading principal minors of each order of A are nonzero.

Example of LU Factorization 2 1 3 2 Suppose we are given: = A 2 0 3 1 0 Then we can write A = LU where: = = U L , 0.5 0.5 1 Let s check that: + + 1 0 2 0 3 1*2 0.5*2 1*0 0*0 + 1*3 0*0.5 0.5*3 1*0.5 + 2 1 3 2 = = = LU 0.5 1 0.5

LU Factorization (in words) Construct a sequence of multiplier matrices (M1, M2, M3, M4, ., MN-1) such that: (MN-1 . M4M3M2M1)A is an upper triangle matrix denoted as U. 1 A A A A A A 0 0 0 21 1 0 0 A suitable candidate for M1is the matrix ( 4 by 4 case): 11 = 1 M 1 0 1 0 3 11 1 0 0 1 4 11

LU Factorization-cont. A suitable candidate for M3is the matrix ( 4 by 4 case): A suitable candidate for M2is the matrix ( 4 by 4 case): 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 ( ( ( ( ) ) ) ) 1 M A = 3 M 0 1 0 23 ( ( ) ) = 2 M 2 1 M M A 1 M A 0 0 1 43 22 2 1 M M A 1 M A 33 0 0 1 24 1 M A 22

LU Factorization-cont. So in this case: U=(M3M2M1)A is upper triangle 1 0 0 1 0 0 0 0 Each of the M matrices is lower triangle ( ( ( ( ) ) ) ) 1 M A ( ) 0 1 0 32 1 = 2 M 1 M A The inverse of M2looks like: 22 1 M A In the pthrow, 0 0 1 42 1 M A Using the fact that the inverse of the each M matrix just requires us to negate the below diagonal terms 22 A = (M3M2M1)-1U = (M1)-1 (M2)-1 (M3)-1U Define L = (M1)-1 (M2)-1 (M3)-1and we are done since the product of lower matrices is a lower matrix