Gaussian Elimination

Download Report

Transcript Gaussian Elimination

Gaussian Elimination
Civil Engineering Majors
Author(s): Autar Kaw
http://numericalmethods.eng.usf.edu
Transforming Numerical Methods Education for STEM
Undergraduates
Naïve Gauss Elimination
http://numericalmethods.eng.usf.edu
Naïve Gaussian Elimination
A method to solve simultaneous linear
equations of the form [A][X]=[C]
Two steps
1. Forward Elimination
2. Back Substitution
Forward Elimination
The goal of forward elimination is to transform the
coefficient matrix into an upper triangular matrix
 25 5 1  x1  106.8 
 64 8 1  x   177.2 

  2 

144 12 1  x3  279.2
5
1   x1   106.8 
25
 0  4.8  1.56  x     96.21

  2 

0
0.7   x3   0.735 
 0
Forward Elimination
A set of n equations and n unknowns
a11x1  a12 x2  a13 x3  ...  a1n xn  b1
a21x1  a22 x2  a23 x3  ...  a2n xn  b2
.
.
.
.
.
.
an1x1  an2 x2  an3 x3  ...  ann xn  bn
(n-1) steps of forward elimination
Forward Elimination
Step 1
For Equation 2, divide Equation 1 by a11 and
multiply by a21 .
 a21 
 a (a11 x1  a12 x2  a13 x3  ...  a1n xn  b1 )
 11 
a21
a21
a21
a21 x1 
a12 x2  ... 
a1n xn 
b1
a11
a11
a11
Forward Elimination
Subtract the result from Equation 2.
a21x1  a22 x2  a23 x3  ...  a2n xn  b2
a21
a21
a21
− a21 x1  a a12 x2  ...  a a1n xn  a b1
11
11
11
_________________________________________________


a21 
a21 
a21
 a22 
a12  x2  ...   a2 n 
a1n  xn  b2 
b1
a11 
a11 
a11


or
'
a22
x2  ...  a2' n xn  b2'
Forward Elimination
Repeat this procedure for the remaining
equations to reduce the set of equations as
a11x1  a12 x2  a13 x3  ...  a1n xn  b1
'
'
a22
x2  a23
x3  ...  a2' n xn  b2'
'
'
a32
x2  a33
x3  ...  a3' n xn  b3'
.
.
.
.
.
.
.
.
.
'
an' 2 x2  an' 3 x3  ...  ann
xn  bn'
End of Step 1
Forward Elimination
Step 2
Repeat the same procedure for the 3rd term of
Equation 3.
a11x1  a12 x2  a13 x3  ...  a1n xn  b1
'
'
a22
x2  a23
x3  ...  a2' n xn  b2'
"
a33
x3  ...  a3"n xn  b3"
.
.
.
.
.
.
"
an" 3 x3  ...  ann
xn  bn"
End of Step 2
Forward Elimination
At the end of (n-1) Forward Elimination steps, the
system of equations will look like
a11 x1  a12 x2  a13 x3  ...  a1n xn  b1
'
'
a22
x2  a23
x3  ...  a2' n xn  b2'
a x  ...  a x  b
"
33 3
"
3n n
.
.
.
"
3
.
.
.
n 1 
 n 1
ann
xn  bn
End of Step (n-1)
Matrix Form at End of Forward
Elimination
a11
0

0



 0
a12
'
22
a1n   x1   b1 
'
' 




 a 2 n x2
b2
  

 a"3n   x3    b3" 
  


      
(n 1 )
  xn  bn(n-1 ) 
0 ann
a13 
'
23
"
33
a
0
a
a


0
0
Back Substitution
Solve each equation starting from the last equation
5
1   x1   106.8 
25
 0  4.8  1.56  x     96.21

  2 

0
0.7   x3   0.735 
 0
Example of a system of 3 equations
Back Substitution Starting Eqns
a11 x1  a12 x2  a13 x3  ...  a1n xn  b1
'
'
a22
x2  a23
x3  ...  a2' n xn  b2'
"
a33
x3  ...  an" xn  b3"
.
.
.
 n 1
.
.
.
n 1 
ann xn  bn
Back Substitution
Start with the last equation because it has only one unknown
( n 1)
n
( n 1)
nn
b
xn 
a
Back Substitution
( n 1)
n
( n 1)
nn
b
xn 
a
xi 
bii 1  ai,ii11 xi 1  ai,ii12 xi 2  ... ai,in1 xn
i 1
aii
i 1
xi 
bi
n
i 1
  aij x j
j i 1
i 1
ii
a
for i  n  1,...,1
fori  n  1,...,1
THE END
http://numericalmethods.eng.usf.edu
Naïve Gauss Elimination
Example
http://numericalmethods.eng.usf.edu
Example: Cylinder Stresses
To find the maximum stresses in a compound cylinder, the following
four simultaneous linear equations need to be solved.
4.2857107

7
4
.
2857

10


 6.5

0

 9.2307105
0
 5.4619105
 0.15384
 4.2857107
6.5
0
4.2857107
  c1   7.887103 
  

5.4619105  c2  
0


0.15384  c3   0.007 

 
5 
 3.605710  c4  
0

0
Example: Cylinder Stresses
In the compound cylinder, the inner cylinder has an internal radius of
a = 5” and an outer radius c = 6.5”, while the outer cylinder has an
internal radius of c = 6.5” and outer radius of b = 8”. Given
E = 30×106 psi, ν = 0.3, and that the hoop stress in outer cylinder is
given by
E
 
1  2

 1 
c3 1    c4  r 2





find the stress on the inside radius of the outer cylinder.
Find the values of c1, c2, c3 and c4 using Naïve Gauss Elimination.
Example: Cylinder Stresses
Forward Elimination: Step 1
 4.2857107 
Row2  
 ( Row1) 
7
 4.285710 
Yields
4.2857107

0

  6.5

0

 9.2307105
0
3.7688105
 0.15384
 4.2857107
6.5
0
4.2857107
  c1   7.887103 
  

5.4619105  c2   7.887103 


 0.007 


c
0.15384
3

 
5 
c
 3.605710   4  
0

0
Example: Cylinder Stresses
Forward Elimination: Step 1
 6.5


Row3  
 ( Row1) 
7
 4.285710 
Yields
4.2857107

0


0

0

 9.2307105
0
3.7688105
 0.29384
 4.2857107
6.5
0
4.2857107
  c1    7.887103 
  

5.4619105  c 2   7.887103 


5.8038103 


c3
0.15384
  

 3.6057105  c 4  
0

0
Example: Cylinder Stresses
Forward Elimination: Step 1
0


Row4  
 ( Row1) 
7
 4.285710 
Yields
4.2857107

0


0

0

 9.2307105
0
3.7688105
 0.29384
 4.2857107
6.5
0
4.2857107
  c1    7.887103 
  

5.4619105  c 2   7.887103 


c3  5.8038103 
0.15384
  

 3.6057105  c 4  
0

0
Example: Cylinder Stresses
Forward Elimination: Step 2
  0.29384 
Row3  
 ( Row2) 
5
 3.768810 
Yields
4.2857107

0


0

0

 9.2307105
0
3.7688105
 4.2857107
0
 26.914
0
4.2857107
  c1    7.887103 
  

5.4619105  c 2   7.887103 


c3  1.195310 2 
0.57968
  

 3.6057105  c 4  
0

0
Example: Cylinder Stresses
Forward Elimination: Step 2
0


Row4  
 ( Row2) 
5
 3.768810 
Yields
4.2857107

0


0

0

 9.2307105
0
3.7688105
0
 4.2857107
 26.914
0
4.2857107
  c1    7.887103 
  

5.4619105  c 2   7.887103 

0.57968  c3  1.195310 2 

 
5 
c
 3.605710   4  
0

0
Example: Cylinder Stresses
Forward Elimination: Step 3
 4.2857107 
Row4  
  ( Row3) 
  26.914 
Yields
4.2857107

0


0

0

 9.2307105
0
3.7688105
0
 4.2857107
 26.914
0
0
  c1    7.887103 
  

5.4619105  c 2   7.887103 


1.195310 2 


c
0.57968
3
 
5  
4 
5.62510  c 4   1.903410 
0
This is now ready for Back Substitution.
Example: Cylinder Stresses
Back Substitution: Solve for c4 using the fourth equation
5.625105 c 4  1.9034104
1.9034104
c4 
5.625105
2
 3.383710
Example: Cylinder Stresses
Back Substitution: Solve for c3 using the third equation
 26.914c3  0.57968c4  1.1953102
1.1953102  0.57968c4
c3 
 26.914
1.1953102  0.57968 3.3837102

 26.914
 2.8469104
Example: Cylinder Stresses
Back Substitution: Solve for c2 using the second equation


3.7688105 c2   4.2857107 c3  5.4619105 c4  7.887103



 
7.887103   4.2857 107 c3  5.4619 105 c4
c2 
3.7688105
 
 
7.887103   4.2857 107  2.8469 10 4  5.4619 105  3.3837 10 2

3.7688 105
 4.2615 103

Example: Cylinder Stresses
Back Substitution: Solve for c1 using the first equation
4.2857107 c1  9.2307105 c2  0c3  0c4  7.887103



 
 7.887103   9.2307105 c2  0 c3  0 c4
c1 
4.2857 107
 7.887103   9.2307105  4.2615 103

4.2857107
 9.2244105

Example: Cylinder Stresses
 c1   9.2244105 
c  
3 
 2    4.261510 
c3   2.846910 4 
  
2 
c4   3.383710 
The solution vector is
The stress on the inside radius of the outer cylinder is then
given by
E 
 1   
 
c 1    c

1  
2
30106

1  0.32
3
4
 2 
 r 

4
 2  1  0.3  
2.846910 1  0.3  3.383710  6.52 



 30683psi
THE END
http://numericalmethods.eng.usf.edu
Naïve Gauss Elimination
Pitfalls
http://numericalmethods.eng.usf.edu
Pitfall#1. Division by zero
10x2  7 x3  3
6 x1  2 x2  3x3  11
5 x1  x2  5 x3  9
0 10  7  x1   3 
6 2
3   x2   11

   
5  1 5   x3   9 
Is division by zero an issue here?
12x1  10x2  7 x3  15
6 x1  5x2  3x3  14
5x1  x2  5x3  9
12 10  7  x1  15
6 5
3   x2   14

   
 5  1 5   x3   9 
Is division by zero an issue here?
YES
12x1  10x2  7 x3  15
6 x1  5x2  3x3  14
24x1  x2  5x3  28
12 10  7  x1  15
6 5
3   x2   14

   
24  1 5   x3  28
12 10  7  x1   15 
0
0
6.5  x2   6.5

   
12  21 19   x3   2
Division by zero is a possibility at any step
of forward elimination
Pitfall#2. Large Round-off Errors
15
10  x1   45 
 20
 3  2.249 7   x   1.751

  2 

1
3   x3   9 
 5
Exact Solution
 x1  1
 x   1
 2  
 x3  1
Pitfall#2. Large Round-off Errors
15
10  x1   45 
 20
 3  2.249 7   x   1.751

  2 

1
3   x3   9 
 5
Solve it on a computer using
6 significant digits with chopping
 x1   0.9625 
 x    1.05 
 2 

 x3  0.999995
Pitfall#2. Large Round-off Errors
15
10  x1   45 
 20
 3  2.249 7   x   1.751

  2 

1
3   x3   9 
 5
Solve it on a computer using
5 significant digits with chopping
 x1   0.625 
 x    1 .5 
 2 

 x3  0.99995
Is there a way to reduce the round off error?
Avoiding Pitfalls
Increase the number of significant digits
• Decreases round-off error
• Does not avoid division by zero
Avoiding Pitfalls
Gaussian Elimination with Partial Pivoting
• Avoids division by zero
• Reduces round off error
THE END
http://numericalmethods.eng.usf.edu
Gauss Elimination with
Partial Pivoting
http://numericalmethods.eng.usf.edu
Pitfalls of Naïve Gauss Elimination
• Possible division by zero
• Large round-off errors
Avoiding Pitfalls
Increase the number of significant digits
•
Decreases round-off error
• Does not avoid division by zero
Avoiding Pitfalls
Gaussian Elimination with Partial Pivoting
• Avoids division by zero
• Reduces round off error
What is Different About Partial
Pivoting?
At the beginning of the kth step of forward elimination,
find the maximum of
akk , ak 1,k ,................, ank
If the maximum of the values is a pk
in the p
th
row, k  p  n, then switch rows p and k.
Matrix Form at Beginning of 2nd
Step of Forward Elimination
a11 a12
 0 a'
22

'
0
a

32




 0 a'n 2
a13

'
23
'
33
a
a




an' 3 an' 4
a1n   x1   b1 
' 
' 


a 2 n x2
b2
   
a3' n   x3    b3' 
   
      
'
  xn  bn' 
ann
Example (2nd step of FE)
6 14
0  7

0 4

0 9
0  17
6   x1   5 
2   x2   6
   
1 11  x3    8 
   
6
8   x4   9 
11 43  x5   3 
5.1 3.7
6
1
12
23
12
Which two rows would you switch?
Example (2nd step of FE)
6 14
0  17

0 4

0
9

0  7
5.1 3.7 6   x1   5 





12 11 43 x2
3
   
12 1 11  x3    8 
   
23 6
8   x4   9 
6
1
2   x5   6
Switched Rows
Gaussian Elimination
with Partial Pivoting
A method to solve simultaneous linear
equations of the form [A][X]=[C]
Two steps
1. Forward Elimination
2. Back Substitution
Forward Elimination
Same as naïve Gauss elimination method
except that we switch rows before each
of the (n-1) steps of forward elimination.
Example: Matrix Form at Beginning
of 2nd Step of Forward Elimination
a11 a12
 0 a'
22

'
 0 a32




 0 a'n 2
a13

'
23
'
33
a
a




a
'
n3
a
'
n4
a1n   x1   b1 
' 
' 


a 2 n x2
b2
   
'
'
a3n   x3    b3 
   
      
'
'




ann   xn  bn 
Matrix Form at End of Forward
Elimination
a11 a12
 0 a'
22

0
0




 0
0
a1n   x1   b1 
'
' 




 a 2 n x2
b2
  

"
"
 a3n   x3    b3 
  


      
(n 1 )
(n-1 )




0 ann   xn  bn 
a13 
'
23
"
33
a
a

0
Back Substitution Starting Eqns
a11 x1  a12 x2  a13 x3  ...  a1n xn  b1
'
'
a22
x2  a23
x3  ...  a2' n xn  b2'
"
a33
x3  ...  an" xn  b3"
.
.
.
 n 1
.
.
.
n 1 
ann xn  bn
Back Substitution
( n 1)
n
( n 1)
nn
b
xn 
a
i 1
xi 
bi
n
i 1
  aij x j
j i 1
i 1
ii
a
for i  n  1,...,1
THE END
http://numericalmethods.eng.usf.edu
Gauss Elimination with
Partial Pivoting
Example
http://numericalmethods.eng.usf.edu
Example 2
Solve the following set of equations
by Gaussian elimination with partial
pivoting
 25 5 1  a1  106.8 
 64 8 1 a   177.2 

  2 

144 12 1  a 3  279.2
Example 2 Cont.
 25 5 1  a1  106.8 
 25 5 1  106.8 
 64 8 1 a   177.2   

64
8
1

177
.
2

  2 



144 12 1  a3  279.2
144 12 1  279.2
1. Forward Elimination
2. Back Substitution
Forward Elimination
Number of Steps of Forward
Elimination
Number of steps of forward elimination is
(n1)=(31)=2
Forward Elimination: Step 1
• Examine absolute values of first column, first row
and below.
25, 64, 144
• Largest absolute value is 144 and exists in row 3.
• Switch row 1 and row 3.
 25 5 1  106.8 
144 12 1  279.2
 64 8 1  177.2    64 8 1  177.2 




144 12 1  279.2
 25 5 1  106.8 
Forward Elimination: Step 1 (cont.)
144 12 1  279.2
 64 8 1  177.2 


 25 5 1  106.8 
144
Divide Equation 1 by 144 and
64
 0.4444 .
multiply it by 64,
144
12 1  279.2 0.4444 63.99 5.333 0.4444  124.1
.
Subtract the result from
Equation 2
64
 63.99
0
1  177.2
8
5.333 0.4444 
124.1
2.667 0.5556  53.10
1
 279.2
Substitute new equation for 144 12
 0 2.667 0.5556  53.10
Equation 2


 25
5
1
 106.8 
Forward Elimination: Step 1 (cont.)
1
 279.2 Divide Equation 1 by 144 and
144 12
 0 2.667 0.5556  53.10
25

 multiply it by 25,
 0.1736 .
144
 25
5
1
 106.8 
144
12 1  279.2 0.1736 25.00 2.083 0.1736  48.47
.
Subtract the result from
Equation 3
Substitute new equation for
Equation 3
25
 25
0
5
1  106.8
2.917
0.8264  58.33
2.083 0.1736  48.47
1
 279.2
144 12
 0 2.667 0.5556  53.10


 0 2.917 0.8264  58.33
Forward Elimination: Step 2
• Examine absolute values of second column, second row
and below.
2.667, 2.917
• Largest absolute value is 2.917 and exists in row 3.
• Switch row 2 and row 3.
1
 279.2
1
 279.2
144 12
144 12
 0 2.667 0.5556  53.10   0 2.917 0.8264  58.33




 0 2.917 0.8264  58.33
 0 2.667 0.5556  53.10
Forward Elimination: Step 2 (cont.)
1
 279.2
144 12
 0 2.917 0.8264  58.33


 0 2.667 0.5556  53.10
0
Divide Equation 2 by 2.917 and
multiply it by 2.667,
2.667
 0.9143 .
2.917
2.917 0.8264  58.33 0.9143 0 2.667 0.7556  53.33
.
Subtract the result from
Equation 3
Substitute new equation for
Equation 3
0
 0
0
2.667 0.5556 
2.667 0.7556 
0
53.10
53.33
 0.2   0.23
1
 279.2 
144 12
 0 2.917 0.8264  58.33 


 0
0
 0.2   0.23
Back Substitution
Back Substitution
1
 279.2 
1 
144 12
144 12
 0 2.917 0.8264  58.33    0 2.917 0.8264




 0
 0
0
 0.2   0.23
0
 0.2 
Solving for a3
 0.2a3  0.23
 0.23
a3 
 0.2
 1.15
 a1   279.2 
a    58.33 
 2 

 a3   0.23
Back Substitution (cont.)
1 
144 12
 0 2.917 0.8264


0
 0.2 
 0
 a1   279.2 
a    58.33 
 2 

 a 3    0.23
Solving for a2
2.917a2  0.8264a3  58.33
58.33  0.8264a3
a2 
2.917
58.33  0.82641.15

2.917
 19.67
Back Substitution (cont.)
1 
144 12
 0 2.917 0.8264


0
 0.2 
 0
 a1   279.2 
a    58.33 
 2 

 a 3    0.23
Solving for a1
144a1  12a2  a3  279.2
279.2  12a2  a3
a1 
144
279.2  1219.67  1.15

144
 0.2917
Gaussian Elimination with Partial
Pivoting Solution
 25 5 1  a1  106.8 
 64 8 1 a   177.2 

  2 

144 12 1  a 3  279.2
 a1  0.2917
a    19.67 
 2 

 a3   1.15 
Gauss Elimination with
Partial Pivoting
Another Example
http://numericalmethods.eng.usf.edu
Partial Pivoting: Example
Consider the system of equations
10x1  7 x2  7
 3x1  2.099x2  6 x3  3.901
5x1  x2  5x3  6
In matrix form
 7 0  x1 
 7 
 10
3.901
 3 2.099 6  x 


  2 = 
 6 
 5
 1 5  x 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
 7 0  x1   7 
 10
 3 2.099 6  x   3.901

 2  

 5
 1 5  x3   6 

7
0  x1   7 
10
 0  0.001 6  x   6.001

 2  

 0
2.5
5  x3   2.5 
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
7
0  x1   7 
10
 0  0.001 6  x   6.001

 2  

 0
2.5
5  x3   2.5 

7
0  x1   7 
10
0
  x    2.5 
2
.
5
5

 2  

 0  0.001 6  x3  6.001
Partial Pivoting: Example
Forward Elimination: Step 2
Performing the Forward Elimination results in:
0   x1   7 
10  7
 0 2. 5
  x    2.5 
5

 2  

 0
0 6.002  x3  6.002
Partial Pivoting: Example
Back Substitution
Solving the equations through back substitution
0   x1   7 
10  7
 0 2. 5
  x    2.5 
5

 2  

 0
0 6.002  x3  6.002
6.002
x3 
1
6.002
2.5  5 x3
x2 
 1
2.5
7  7 x 2  0 x3
x1 
0
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
 x1   0 
X  calculated   x2    1
 x3   1 
X  exact
 x1   0 
  x 2    1
 x3   1 
THE END
http://numericalmethods.eng.usf.edu
Determinant of a Square Matrix
Using Naïve Gauss Elimination
Example
http://numericalmethods.eng.usf.edu
Theorem of Determinants
If a multiple of one row of [A]nxn is added or
subtracted to another row of [A]nxn to result in
[B]nxn then det(A)=det(B)
Theorem of Determinants
The determinant of an upper triangular matrix
[A]nxn is given by
det A  a11  a22  ...  aii  ...  ann
n
  a ii
i 1
Forward Elimination of a
Square Matrix
Using forward elimination to transform [A]nxn to an
upper triangular matrix, [U]nxn.
Ann  U  nn
det  A  det U 
Example
Using naïve Gaussian elimination find the
determinant of the following square
matrix.
 25 5 1
 64 8 1


144 12 1
Forward Elimination
Forward Elimination: Step 1
 25 5 1
 64 8 1


144 12 1
25
Divide Equation 1 by 25 and
64
 2.56 .
multiply it by 64,
25
5 1 2.56  64 12.8 2.56
Subtract the result from
Equation 2
64
8
 64 12.8
0  4.8
Substitute new equation for
Equation 2
5
1 
 25
 0  4.8  1.56


144 12
1 
.
1
2.56
 1.56
Forward Elimination: Step 1 (cont.)
5
1 
 25
 0  4.8  1.56 Divide Equation 1 by 25 and

 multiply it by 144, 144  5.76 .
144 12
1 
25
25 5 1 5.76  144 28.8 5.76
.
Subtract the result from
Equation 3
Substitute new equation for
Equation 3
144 12
 144 28.8
0  16.8
1
5.76
 4.76
5
1 
25
 0  4.8  1.56


 0  16.8  4.76
Forward Elimination: Step 2
5
1 
25
 0  4.8  1.56


 0  16.8  4.76
0
Divide Equation 2 by −4.8
and multiply it by −16.8,
 16 .8
 3 .5 .
 4 .8
 4.8  1.56  3.5  0  16.8  5.46
.
Subtract the result from
Equation 3
Substitute new equation for
Equation 3
0
 0
0
 16.8  4.76
 16.8  5.46
0
0.7
5
1 
25
 0  4.8  1.56


 0
0
0.7 
Finding the Determinant
After forward elimination
5
1 
 25 5 1
25
 64 8 1   0  4.8  1.56




0
0.7 
144 12 1
 0
.
det A   u11  u22  u33
 25  4.8 0.7
 84.00
Summary
-Forward Elimination
-Back Substitution
-Pitfalls
-Improvements
-Partial Pivoting
-Determinant of a Matrix
Additional Resources
For all resources on this topic such as digital audiovisual
lectures, primers, textbook chapters, multiple-choice tests,
worksheets in MATLAB, MATHEMATICA, MathCad and
MAPLE, blogs, related physical problems, please visit
http://numericalmethods.eng.usf.edu/topics/gaussian_elimi
nation.html
THE END
http://numericalmethods.eng.usf.edu