Chapter: 3c System of Linear Equations Dr. Asaf Varol

1

Pivoting

Some disadvantages of Gaussian elimination are as follows: Since each result follows and depends on the previous step, for large systems the errors introduced due to round off (or chop off) errors lead to loss of significant figures and hence to in accurate results. The error committed in any one step propagates till the final step and it is amplified. This is especially true for ill-conditioned systems. Of course if any of the diagonal elements is zero, the method will not work unless the system is rearranged so to avoid zero elements being on the diagonal. The practice of interchanging rows with each other so that the diagonal elements are the dominant elements is called

partial pivoting

. The goal here is to put the largest possible coefficient along the diagonal by manipulating the order of the rows. It is also possible to change the order of variables, i.e. instead of letting the unknown vector {X} T = (x, y, z) we may let it be {X} T = ( y, z, x) When this is done in addition to partial pivoting, this practice is called

full pivoting

. In this case only the meaning of each variable changes but the system of equations remain the same.

2

Example E3.4.2

Consider the following set of equations 0.0003 x + 3.0000 y = 2.0001

1.0000 x + 1.0000 y = 1.0000

The exact solution to which is x = 1/3 and y = 2/3 Solving this system by Gaussian elimination with a three significant figure mantissa yields x = -3.33, y = 0.667

with four significant figure mantissa yields x = 0.0000 and y = 0.6667

3

Example E3.4.2 (continued)

Partial pivoting

(switch the rows so that the diagonal elements are largest) 1.0000 x + 1.0000 y = 1.0000

0.0003 x + 3.0000 y = 2.0001

The solution of this set of equations using Gaussian elimination gives y = 0.667 and x = 0.333 with three significant figure arithmetic y = 0.6667 and x = 0.3333 with four significant figure arithmetic 4

Example E3.4.3

Problem:

Apply full pivoting to the following system to achieve a well conditioned matrix.

3x + 5y - 5z = 3 2x - 4y - z = -3 6x - 5y + z = 2     3 2 6 5  4  5   1 5 1      

x y

 

z

 3   2 3   5

Example E3.4.3 (continued)

Solution:

First switch the first column with the second column, then switch the first column with the third column to obtain -5z + 3x + 5y = 3 - z + 2x - 4y = -3 z + 6x - 5y = 2 Then switch second row with the third row -5z + 3x + 5y = 3 z + 6x - 5y = 2 -z + 2x - 4y = -3 Yielding finally a well conditioned system given by       1 5 1 3 6 2   5 5 4       z x y       3 2 3   6

Gauss – Jordan Elimination

This method is very similar to Gaussian elimination method. The only difference is in that the elimination procedure is extended to the upper diagonal elements so that a backward substitution is no longer necessary. The elimination process begins with the augmented matrix, and continued until the original matrix turns into an identity matrix, of course with necessary modifications to the right hand side.

In short our goal is to start with the general augmented system and arrive at the right side after appropriate algebraic manipulations.

 a    a a aa 21 31 Start a 12 a 22 a 32 a 13 a 23 a 33       c c c 1 2 3   ===> The solution can be written at once as x 1  c 1 *

;

x 2  c * 2

;

x 3  c * 3  1    0 0 arrive 0 1 0 0 0 1        c c * c * 2 * 3 1   7

Pseudo Code for Gauss – Jordan Elimination

do for k = 1 to n do for j= k+1 to n+1 end do a kj = a kj /a kk do for i = 1 to n ;

i is not equal to k

do for j = k+1 to n+1 a ij = a ij - (a ik )(a kj ) end do end do end do cc-----

The solution vector is saved in a(i,n+1), i=1, to n

! Important note ! a(i,n+1) represents the ! the right hand side 8

Example E3.4.4

Problem:

Solve the problem of Equations (3.4.3) using Gauss-Jordan elimination 4x 1 + x 2 + x 3 = 6 -x 1 - 5x 2 + 6x 3 = 0 2x 1 - 4x 2 + x 3 = -1

Solution:

To do that we start with the augmented the coefficient matrix

a

     4  1 2   1 5 4 1 6 1 6 0  1     (i) First multiply the first row by 1/4, then multiply it by -1 and subtract the result from the second row and replace the result of the subtraction with the second row. Similarly multiply the first row by 2 and subtract the result from the third row and replace the third row with the result of the subtraction. These operations lead to   a (1)      1 .

0 0 0 0 .

25  4 .

75  4 .

50 0 0 .

.

25 6 .

25 50 1 1  .

.

5 5 4     9

Example E3.4.4 (continued)

(ii) Multiply the second row by -1/4.75, then multiply the second row by 0.25 and subtract the result from the first row and replace the result with the first row. Similarly multiply the second row by -4.5 and subtract the result from the third row and replace the result with the third row to obtain   a      1 0 0 0 1 0  0 .

5 5789  1 .

3158 .

4211   1 .

0 5789 5 .

.

3158 4211     (iii) Multiply the third row by -1/5.4211, then multiply the third row by 0.5789 and subtract the result from the first row, and then replace the first row by the result. Similarly multiply the third row by -1.3158 and subtract the result from the second row, and then replace the second row by the result to finally arrive at   a      1 0 0 0 1 0 0 0 1     Hence we found the expected solution x 1 =1.0, x 2 =1.0, and x 3 = 1.0. Note, however, that the solutions of 1.00 were just the solutions. Of course, each matrix will have its own solutions, and 10 not always equal to one!

Finding the Inverse using Gauss-Jordan Elimination

The inverse, [B] of a matrix, [A], is defined such that [A][B] = [B][A] = [I] where [I] is the identity matrix. In other words, to find the inverse of a matrix we need to find the elements of the matrix [B] such that when [A] is multiplied by [B] the result should equal the identity matrix. If we recall what matrix multiplication is, this amounts to the same thing as saying the first column of [B] multiplied by [A] must be equal to the first column of [I] ; the second column of [B] multiplied by [A] must be equal to the second column of [I] and so on.

For summary, using Gauss-Jordan elimination     4 1 2 1 5 6 | 4 1 1 | 1 0 0 | 0 1 0 0 0   1   We obtain     1 0 0 0 0 1 0 0 1 | 0.185

| 0.1262

0.04854

0.01942

| 0.1359

0.17476

0.1068

0.2427

0.1845

    11

LU Decomposition

Lower and Upper (LU) triangular decomposition technique is one of the most widely used techniques used for solution of linear systems of equations due to its generality and efficiency. This method calls for first decomposing a given matrix into a product of a lower and an upper triangular matrices such that [A] = [L][U] Given that [A]{X} = {R} [U]{X} = {Q} for an unknown vector {Q} The additional unknown vector {Q} can be found from [L]{Q} = {R} The solution procedure (i.e. the algorithm) can now be summarized as (i) Determine [L] and [U] for a given system of equations [A]{X} = {R} (ii) (iii) Calculate {Q} from [L]{Q} = {R} by forward substitution Calculate {X} from [U]{X} = {Q} by backward substitution 12

Example E3.4.5

[A] =     2  1 3   3 5 4 1  1 2     ; {R} =   12  8 16   [L] =     2  1 3 To find the solution we start with [L]{Q} = {R}; that is 0 0  0 4    ; 2 q 1 -q 1 3q 1 + 0 + q 2 + 7q /2 + 0 2 + 0 /2+ 4q 3 = 12 = -8 = 16 We compute {Q} by forward substitution q 1 = 12/2 = 6 Hence {Q} T = {6, -4, 3} q 2 = (-8 + q 1 )/(1/2) = -4 Then we solve [U]{X} = {Q} by backward substitution [U] =  1    0 0  / 1 0 q 3 = (16 - 3q 1 - 7q 2 /2)/4 = 3  /  1 1    x 1 - (5/2) x 2 0 + x 2 0 + 0 +(1/2) x 3 - x 3 + x 3 = 6 = -4 = 3 x 3 = 3 x 2 = -4 + x 3 = -1 Hence the final solution is x 1 = 6 + (5/2) x 2 -(1/2) x 3 = 2 {X} T = { 2, -1, 3 } 13

Crout Decomposition

We illustrate Crout-decomposition in detail for a 3x3 general matrix  l    l l 11 21 31 [L] [U] = [A] 0 l 22 l 32 0 0 l 33      1    0 0 u 1 0 12 u 13 u 23 1    

=

    a a a 11 21 31 a 12 a 22 a 32 a 13 a a 23 33     Multiplying [L] by [U] and equating the corresponding elements of the product with that of the matrix [A] we get the following equations l 11 = a 11 ; l 21 = a 21 ; l 31 = a 31 ;

in index notation: l

i1 = a i1 , I=1,2, ... , n (j=1, first column)

(The first column of [L] is equal to the first column of [A]

equation l 11 u 12 = a 12 l 11 u 13 = a 13 ====> ====>

in index notation

solution u u 12 13

u 1j = a 1j /l 11

= a = a 12 13 /l /l 11 11

j=2, 3, ... , n (i=1, first row)

l 21 u 12 + l 22 = a 22 l 31 u 12 + l 32 = a 32 ====> ====>

in index notation

l l 22 32 = a = a 22 32 - l - l 21 31 u u 12 12

l i2 = a i2 - l i1 u 12 ; i=2,3, ... , n (j=2, second column)

14

Pseudo Code for LU-Decomposition

cc--- Simple coding

do for i = 1 to n l i1 = a i1 end do do for j = 2 to n u 1j = a 1j /l 11 end do do for j = 2 to n-1 l jj = a jj j  1 - l u k  1 kj c--- Forward substitution q 11 = r 1 / l 11 do for i = 2 to n q i = r i i  1 - ( l q j  1 j )/l ii end do c--- Back substitution xn = qn do for i = n-1 to 1 x i = q i - n  1 u x ij j end do do for i = j+1 to n l ij = a ij j  1 - l u k  1 kj end do u ji = a ji j  1 - ( l u k 1 ki )/l jj end do l nn = a nn n k  1 - l u kn 15

Cholesky Decomposition for Symmetric Matrices

For symmetric matrices the LU-decomposition can be further simplified to take advantage of the symmetric property that [A] = [A] T ; or a ij = a ji (3.4.9) Then, it follows that That is [A] = [L][U] = [A] T = ([L][U]) T = [U] T [L] T [L] = [U] T ; [U] = [L] T (3.4.10) or in index notation u ij = l ji (3.4.11) It should be noted that the diagonal elements of [U] are no longer, necessarily, equal to one.

Comment:

This algorithm should be programmed using complex arithmetic

. 16

Tridiagonal Matrix Algorithm (TDMA),

also known as

Thomas Algorithm

[A] =     

b a

2 0 0 1

c

1

b

2

a

3 0 0

c

2

b

3

a

4 0 0

c

3

b

4      If we apply LU-decomposition to this matrix taking into account of its special structure and let [L] =      e f 0 0 1 2 0 f 2 e 3 0 f e 0 0 3 4 0 0 0 f 4      [U] =  1     0 0 0 g 1 1 0 0 0 g 2 1 0 0 0 g 3 1      By multiplying [L] and [U] and equating the corresponding elements of the product to that of the original matrix [A] it can be shown that f 1 = b 1 g 1 = c 1 /f 1 e 2 = a 2 and for k=2,3, ... , n-1 e k = a k f k = b k - e k g k-1 g k = c k / f k further e f n n = a = b n n - e n g n-1 17

TDMA Algorithm

c----- Decomposition c 1 = c 1 /b 1 do for k=2 to n-1 b k = b k - a k c k-1 c k = c k /b k end do b n = b n - a n c n-1 c----- (note here r k ; k=1,2,3, ..., n is the right hand side of the equations) c----- Forward substitution r 1 = r 1 /b 1 do for k = 2, n r k = (r k - a k r k-1 )/b k end do c----- Back substitution x n = r n do for k = (n-1) to 1 x k = (r k - c k x k+1 ) end do 18

Example E3.4.6

Problem:

Determine the LU-decomposition for the following matrix using TDMA [A] =  2    1 0 2 4 1 0  4 2   

Solution:

n=3 c 1 = c 1 /b 1 = 2/2 = 1 k = 2 b 2 = b 2 - a 2 c 1 = 4 - (1)(1) = 3 c 2 = c 2 /b 2 = 4/3 k = n = 3 b 3 = b 3 - a 3 c 2 = 2 - (1)(4/3) = 2/3 The other elements remain the same; hence [L] =  2    1 0 0 3 1 0 0     ; [U] =  1    0 0 1 1 0 0 1     It can be verified that the product [L][U] is indeed equal to [A].

19

Iterative Methods for Solving Linear Systems

Iterative methods

are those where an initial guess is made for the solution vector followed by a correction sequentially until a certain bound for error tolerance is reached. The concept of iterative computations was introduced in Chapter 2 for nonlinear equations.

The idea is the same here except that we deal with linear systems of equations. In this regard the fixed point iteration method presented in Chapter 2 for two (nonlinear) equations and two unknowns will be repeated here, because that is precisely what Jacobi iteration is about.

20

Jacobi Iteration

Two linear equations, in general, can be written as a 11 x 1 + a 12 x 2 = c 1 a 21 x 1 + a 22 x 2 = c 2 We now rearrange these equations such that each unknown is written as a function of the others, i.e.

x 1 = ( c 1 - a 12 x 2 )/ a 11 = f(x 1 ,x 2 ) x 2 = ( c 2 - a 21 x 1 )/ a 22 = g(x 1 ,x 2 ) which are in the same form as those used for the fixed point iteration (Chapter 2). The functions f(x 1 ,x 2 ) and g(x 1 ,x 2 ) are introduced for generality. It is implicitly assumed that the diagonal elements of the coefficient matrix are not zero. In case there is a zero diagonal element this should be avoided by changing the order of the equations, i.e. by

pivoting

. Jacobi iteration calls for starting with a guess for the unknowns, x 1 and x 2 , and then finding new values using these equations iteratively. If certain conditions are satisfied this iteration procedure converges to the right answer as it did in case of fixed point iteration method.

21

Example E3.5.1

Problem:

Solve the following system of equations using Jacobi iteration 10 x 1 + 2x 2 + 3x 3 = 23 2x 1 - 10x 2 + 3x 3 = -9 - x 1 - x 2 + 5x 3 = 12

Solution:

Rearranging x 1 = ( 23 - 2x 2 - 3x 3 )/10 x 2 = ( -9 - 2x 1 - 3x 3 )/(-10) x 3 = (12 + x 1 + x 2 )/5 Starting with a somewhat arbitrary guess, x 1 = 0 , x 2 = 0 , and x 3 = 0. and iterating yield in the values shown in the table given below ITER X 1 X 2 X 3 Error Norm, E = i n   1 x i new  x i old 0 0 0 0 -- 1 2.300000 0.900000 2.400000 5.600000

2 1.400000 2.080000 3.040000 2.720000

3 0.972000 2.092000 3.096000 4.960001E-01 4 0.952800 2.023200 3.012800 1.712000E-01 5 0.991520 1.994400 2.995200 8.512014E-02 6 1.002560 1.996864 2.997184 1.548803E-02 7 1.001472 1.999667 2.999885 6.592035E-03 8 1.000101 2.000260 3.000228 2.306700E-03 9 0.9998797 2.000089 3.000072 5.483031E-04 10 0.9999606 1.999998 2.999994 2.506971E-04 22

Convergence Criteria for Jacobi Iteration

 f  x 1   f  x 2 =  a 12 /a 11  < 1  g  x 1   g  x 2 =  a 21 /a 22  < 1 In general a sufficient (but not necessary) condition for convergence of Jacobi iteration is ( j n   1, a ij ) / a ii  1 In other words, the sum of the absolute value of the off diagonal elements on each row must be less than the absolute value of the diagonal element of the coefficient matrix.

These matrices are called

strictly diagonally dominant

matrices. The reader can show easily that the matrix of example E3.5.1 does satisfy the criteria given.

23

Gauss - Seidel Iteration

This method is essentially the same as the Jacobi iteration method except that the new values of the variables (i.e. the most recently updated value) is used in subsequent calculations without waiting for completion of one iteration for all variables. In cases the convergence criteria is satisfied, Gauss-Seidel iteration takes fewer iteration to achieve the same error bound compared to Jacobi iteration. To illustrate how this procedure works we present the first two iterations of the Gauss-Seidel method applied to the example problem: x 1 = ( 23 - 2x 2 - 3x 3 )/10 x 2 = ( -9 - 2x 1 - 3x 3 )/(-10) x 3 = (12 + x 1 + x 2 )/5 ITER=0, x1=0, x2=0, x3=0 (initial guess) ITER=1, x 1 = ( 23 - 0 -0 )/10 = 2.3

24

Gauss - Seidel Iteration

(II) Now this new value for x 1 is used in the second equation as opposed to using x 1 = 0 (the previous value of x 1 ) in Jacobi iteration, that is x 2 = [-9 - (2)(2.3) - 0 )/(-10) = 1.36

Similarly in the third equation the most recent values of x 1 and x 2 are used. Hence, x 3 = (12 + 2.3 + 1.36 )/5 = 3.132

And so an for the rest of the iterations ITER=2 x 1 = [23 - (2)(1.36) - (3)(3.132) ] /10 x 2 = [ -9 - (2)(0.8884) - (3)(3.132) ] /(-10) x 3 = (12 + 0.8884 + 2.10728) /5 = 0.8884

= 2.10728

= 2.998256

25

Example E3.5.2

Problem:

Solve the following system of equations (from Example E3.5.1) using Gauss Seidel iteration.

10 x 1 + 2x 2 + 3x 3 = 23 2x 1 - 10x 2 + 3x 3 = -9 - x 1 - x 2 + 5x 3 = 12

Table E3.5.1 Convergence trend of Gauss-Seidel method for Example E3.5.1

ITER X 1 X 2 X 3 Error Norm, E 0 0 0 0 -- 1 2.300000E+00 1.360000E+00 3.132000E+00 6.792000E+00 2 1.088400E+00 2.057280E+00 3.029136E+00 2.011744E+00 3 9.798032E-01 2.004701E+00 2.996901E+00 1.934104E-01 4 9.999894E-01 1.999068E+00 2.999811E+00 2.872980E-02 5 1.000243E+00 1.999992E+00 3.000047E+00 1.412988E-03 6 9.999875E-01 2.000012E+00 3.000000E+00 3.223419E-04 7 9.999977E-01 1.999999E+00 3.000000E+00 2.276897E-05 8 1.000000E+00 2.000000E+00 3.000000E+00 3.457069E-06 9 1.000000E+00 2.000000E+00 3.000000E+00 3.576279E-07 10 1.000000E+00 2.000000E+00 3.000000E+00 0.000000E+00 26

Convergence Criteria for Gauss - Seidel Iteration

Convergence criteria for Gauss-Seidel method are more relaxed than the equation for Jacobi Iteration convergence. Here a sufficient condition is 

j n

  1 ,

a ij



a ii

     1 1

for all equations except one

.

for at least one equation

.

This condition is also known as the Scarborough criterion. We also introduce the following theorem concerning the convergence of Gauss-Seidel iteration method.

Theorem: Let [A] be a real symmetric matrix with positive diagonal elements. Then the Gauss-Seidel method solving [A]{x} = {c} will converge for any choice of initial guess if and only if [A] is a positive definite matrix.

Definition: If for all nonzero complex vectors {x}, {x} T [A]{x} > 0 for a real symmetric matrix [A] then [A] is said to be positive definite matrix. (Note also that [A] is positive definite if and only if all of its eigenvalues are real and positive.

27

Relaxation Concept

The iterative methods such as Jacobi and Gauss-Seidel iteration are successively prediction correction methods in that an initial guess is corrected to compute a new approximate solution, then the new solution is corrected and so on. The Gauss-Seidel method, for example can be formulated as

x i



k

 1  

x i

 

x i

where k denotes the number of iterations and  x i (k) is the correction to the k th solution. Which is obtained from 

x i

 

x i e



x i

 To interpret and understand this equation better we let

x i e



x i new

,

x i



x i old

and write it as x i  new    x i new   1    x i old 28

Example E3.5.4

Problem:

Show that the Gauss-Seidel iteration method converges in 50 iterations when applied to the following system without relaxation (i.e.  = 1) x 1 + 2x 2 = -5 x 1 – 3x 2 = 20 Find an optimum relaxation factor that will minimize the number of iterations to achieve a solution with the error bound.

|| E i || < 0.5 x 10 -5 where || E i || = 2 

i

 1

x i new



x i old

Solution:

A few iterations are shown here with  = 0.8

Let x 1 old = 0 and x 2 old = 0 (initial guess) x 1 new = (-5 - 2x 2 old ) = -5 After a few iterations we reach Under relax x 2 new = (0.8)(-4.853) + (0.2)(-6.4) = -5.162

x 2 old = -5.162 (switch old and new) 29

Example E3.5.4 (continued)

Table E3.5.4b Number of Iterations versus Relaxation Factor for Example E3.5.4

 0.20

0.40

0.60

0.70

0.80

0.90

1.00

1.05

1.10

Number of iterations 62 30 17 14 11 15 50 84 >1000 It is seen that minimum number of iterations is obtained with a relaxation factor of  = 0.80.

30

Case Study - The Problem of Falling Objects

T1 T2 m1 m2 m3 g - gravity

-m

1

a = c

1

v

2

- m

1

g - T

1

-m

2

a = c

2

v

2

- m

2

g - T

1

+ T

2

-m

3

a = c

3

v

2

- m

3

g + T

2 31

Case Study - The Problem of Falling Objects

These three equations can be used to solve for any three of the unknowns parameters: m1, m2, m3, c1, c2, c3, v, a, T1, and T2. As an example here we set up a problem by specifying the masses, the drag coefficients, and the acceleration as m 1 = 10 kg, m 2 = 6 kg, m 3 = 14 kg, c 1 = 1.0 kg/m, c 2 = 1.5 kg/m, c 3 = 2.0 kg/m, and a = 0.5g and find the remaining three unknowns; v, T 1 and T 2 . The acceleration of gravity g=9.81 m/sec 2 .

Substituting these values in the given set of equations and rearranging yield v 2 - T 1 = 49.05

1.5 v 2 + T 1 - T 2 = 29.43

2.0 v 2 + T 2 = 68.67

v 2 = 32.7 (v = 5.72 m/s), T 1 = -16.35 N, T 2 = 3.27 N (N=newton=kg.m/sec 2 ) 32

References

1.

2.

3.

4.

5.

6.

Celik, Ismail, B., “Introductory Numerical Methods for Engineering Applications”, Ararat Books & Publishing, LCC., Morgantown, 2001 Fausett, Laurene, V. “Numerical Methods, Algorithms and Applications”, Prentice Hall, 2003 by Pearson Education, Inc., Upper Saddle River, NJ 07458 Rao, Singiresu, S., “Applied Numerical Methods for Engineers and Scientists, 2002 Prentice Hall, Upper Saddle River, NJ 07458 Mathews, John, H.; Fink, Kurtis, D., “Numerical Methods Using MATLAB” Fourth Edition, 2004 Prentice Hall, Upper Saddle River, NJ 07458 Varol, A., “Sayisal Analiz (Numerical Analysis), in Turkish, Course notes, Firat University, 2001 http://mathonweb.com/help/backgd3e.htm

33