Transcript Gauss Elimination Procedure

PART 3.
Linear Algebraic Equations
(선형대수 방정식)
Se-Woong Chung, Chungbuk National University
1
Lecture Overview
 What is linear algebraic equations?
•
Determine the values x1,x2,…,xn that satisfy a set of equations
a11 x1  a12 x2  ............  a1n xn  b1
a21 x1  a22 x2  ............  a2 n xn  b2
.
.
an1 x1  an 2 x2  ............  ann xn  bn
 Engineering Applications (공학적 응용)
• Conservation laws for mass, momentum, and energy
• A large set of equations that must be solved simultaneously
• Related numerical techniques
: Regression Analysis, Spline interpolation
Se-Woong Chung, Chungbuk National University
2
 Applications to Environmental Engineering
(a) lumped variable system (CSTRs)
(b) distributed variable system that involves a continuum
미지수?
방정식?
Se-Woong Chung, Chungbuk National University
3
Mathematical Background
 A matrix
Se-Woong Chung, Chungbuk National University
4
 Special Types of Square Matrices
• 대칭행렬 (symmetric matrix) p.172
ai , j  a j ,i
• 대각행렬 (diagonal matrix)
• 단위행렬 (identity matrix)
1

 1



A 

1 


1


a11

a22
A  



• 상부삼각행렬 (upper triangular matrix)
• 하부삼각행렬 (lower triangular matrix)
• 띠행렬 (banded matrix)
• 삼중대각행렬 (tridiagonal matrix)
 a11 a12
a
a
A   21 22

a32


Se-Woong Chung, Chungbuk National University
a33
a11 a12

a22
A  



 a11
a
a
A   21 22
a31 a32

a41 a42
a23
a33
a43





a44 
a13
a23
a33
a33
a43
a14 
a24 
a34 

a44 





a44 



a34 

a44 
5
 Matrix Operating Rules
• Multiplying Two Matrices
Se-Woong Chung, Chungbuk National University
6
 Matrix Operating Rules
• 덧셈의 교환법칙
A  B  B  A
• 뺄셈의 결합법칙
E  F   F   E
• 곱셈의 교환법칙
(AB)C   A(BC )
• 곱셈의 분배법칙
A(B  C )  AB  AC 
• 역행렬 (inverse)
AA1  A1 A  I 
• 전치행렬(transpose)
• 행렬의 trace
 AT
n
tr[ A]   ai ,i
i 1
Se-Woong Chung, Chungbuk National University
7
PART 3 Lecture Contents
• Chap. 8 : Gauss Elimination (가우스 소거법)
• Chap. 9: LU Decomposition (LU 분해법)
• Chap.10: Special Matrices and Gauss Seidel(반복법)
• Chap.11: Engineering Applications
Se-Woong Chung, Chungbuk National University
8
PART 3. Linear Algebraic Eq.
CHAPTER 8
Gauss Elimination Method
Se-Woong Chung, Chungbuk National University
9
1. Solving Small Numbers of Equations
• The Graphical method
(도식적 방법)
plot x1 versus x2
Ex 8.1
• Cramer’s rule
see text
Ex 8.3
• Elimination of unknowns
(미지수의 소거)
Ex 8.4
Se-Woong Chung, Chungbuk National University
10
 singular & ill-conditioned (특이 행렬시스템 )
(a) No solution (parallel) : 행렬식 (Determinant) = 0
(b) Infinite solutions(coincident) : 행렬식 (Determinant) = 0
(a) and (b) are called “singular”
Ex 8.2
(c) ill-conditioned system (close): 행렬식 (Determinant)  0
Se-Woong Chung, Chungbuk National University
11
 Cramer’s rule (Ex. 8.3)
AX  B
x1 
a11 a12
D  a21 a22
a13
a23
a31 a32
a33
 a11 a12
a
 21 a22
a31 a32
a13   x1  b1 
   

a23   x2   b2 
a33   x3  b3 
식 8.4 참조
예제 8.4 미지수 소거법 = Cramer’s rule
b1
a12
a13
a11
b1
a13
a11
a12
b1
b2
a22
a23
a12
b2
a23
a12
a22
b2
b3
a32
a33
a13
b3
a33
a13
a23
b3
D
x2 
D
x3 
Se-Woong Chung, Chungbuk National University
D
12
2. Naive Gauss Elimination
 General form
a11 x1  a12 x2  ............  a1n xn  b1
a21 x1  a22 x2  ............  a2 n xn  b2
.
.
an1 x1  an 2 x2  ............  ann xn  bn
 Solution procedure
1. Forward Elimination
(미지수의 전진소거 )
2. Back Substitution
(후진대입)
Se-Woong Chung, Chungbuk National University
13
 Gauss Elimination Procedure (1/3)
• Forward Elimination
a11x1  a12 x2  ............  a1n xn  b1 ...... (1)
Pivot element
Pivot row
a21x1  a22 x2  ............  a2 n xn  b2 ...... (2)
.
.
an1 x1  an 2 x2  ............  ann xn  bn ...... (n)
•
Step1. 두번째 방정식 에서 n번째 방정식(식 2~ n)까지 x1 소거
(2)식 - a21/a11*(1)식  (2)’
(3)식 - a31/a11*(1) 식  (3)’
(n)식 - an1/a11*(1) 식  (n)’
a22’
a2n’
b2’
a21
a21
a21
(a22 
a12 ) x2  .......  (a2 n 
a1n ) xn  b2 
b1 ...... (2)'
a11
a11
a11
Se-Woong Chung, Chungbuk National University
14
 Gauss Elimination Procedure (2/3)
a11x1  a12 x2  ............  a1n xn  b1 ...... (1)
a'22 x2  ............  a'2 n xn  b'2 ...... (2)'
Pivot element
Pivot row
a'32 x2  ............  a'3n xn  b'3 ...... (3)'
.
.
a'n 2 x2  ............  a'nn xn  b'n ...... (n)'
•
Step 2. (3)’식부터 (n)’식까지 x2 소거
(3)’식 - a’32/a’22*(2)’식  (3)’’
(n)’식 - a’32/a’22*(n)’식  (n)’’
n번째 방정식으로부터 xn-1을 소거하기까지
“동일한 작업을 반복”
Se-Woong Chung, Chungbuk National University
15
 Gauss Elimination Procedure (3/3)
• Final form
a11 x1  a12 x2  a13 x3  .....  a1n xn  b1
a '22 x2  a '23 x3  .....  a'2 n xn  b'2
a ' '33 x3  .......  a' '3n xn  b' '3
.
(4)
(5)
(6)
.
( n 1)
ann
xn  bn( n 1) (7)
• Back Substitution
bnn 1
xn  n 1
ann
( i 1)
i
b
xi 

Number of modifications
n
( i 1)
a
 ij x j
j i 1
( i 1)
ii
a
for i  n  1, n  2,...,1
Se-Woong Chung, Chungbuk National University
16
 Gauss Elimination Algorithm
Pivot row k=1 to n-1
Row I= k+1 to n
Forward
Elimination
column j= k+1 to n
(2)’ = (2) – factor * (1)
bnn 1
xn  n 1
ann
Backward
Substitution
Ex.8.5 (수계산, 프로그램)
b
xi 
Se-Woong Chung, Chungbuk National University
(i1)
i

n
a iji

j i
a
 1
(i1)
ii
( 1)
xj
for i  n  1,...,
1
17
 Determinant Evaluation Using Gauss Elimination
Gauss 소거법의 전진소거 후 얻은 삼각행렬로부터
행렬식은 대각선 요소들의 곱으로 주어짐
a11
a12
a13
D 0
0
a22
0
a23  a11a22a33
a33
General form
'
''
( n 1)
D  a11a22
a33
     ann
(1) p
P = 행이 pivot화된 횟수
Se-Woong Chung, Chungbuk National University
18
3. Pitfalls of Gauss Elimination Method
 Gauss 소거법의 문제점
1) 8.3.1 Division by Zero (0으로 나누는 경우)
In case a coefficient of x is 0 or close to 0 (a11 = 0)
2) 8.3.2 Round-Off Errors (반올림 오차)
Use of limited number of significant figures causes error
development
예제 8.9
3) 8.3.3 Ill-Conditioned Systems (불량조건 시스템)
Small changes in coefficients result in large changes in the
solution
예제 8.6
Effect of scale on the determinant 예제 8.7 예제 8.8
Se-Woong Chung, Chungbuk National University
19
4. Techniques for Improving Solutions
1) Use of more significant figures for ill-conditioning
: computation price increase (계산량 및 메모리 증가)
2) Partial Pivoting
 Before each row is normalized, determine the largest available
coefficient in the column below the pivot element, and then
switch the rows so that the largest element is the pivot element
 Ex. 8.9
0.0003x1  3.0000 x2  2.0001
Switch the
pivot
element
1.0000 x1  1.0000 x2  1.0000
1.0000 x1  1.0000 x2  1.0000
0.0003x1  3.0000 x2  2.0001
Se-Woong Chung, Chungbuk National University
20
Pseudocode to Implement Partial Pivoting
초기 Pivot행 가정
초기 Pivot요소 가정
최대값으로 Pivot 행 교체
새로운 Pivot 행과 초기 가정의 위치 교체
Se-Woong Chung, Chungbuk National University
21
각행에서 제일
큰 계수 si 저장
Gauss Elimination
with Partial
Pivoting
전진소거
정규화와
피벗행 검사
피벗행 교체
부분피벗화 호출
대각선행 검사
(0에 가까우면 특이행렬 발생)
후진대입
Se-Woong Chung, Chungbuk National University
22
Ex. 8.11 Use of Gauss Elimination for a team of Parachutists
Se-Woong Chung, Chungbuk National University
23
Ex. 8.11: Development of linear algebraic equations
Parachutist
Mass
(kg)
Drag Coeff.
(kg/s)
1
70
10
2
60
14
3
40
17
m1 g  T  c1v
 m1a
m2 g  T  c2 v  R  m2 a
m3 g
 c3v  R  m3 a
70 1 0   a  636
60  1 1  T   518

   
40 0  1 R  307
Se-Woong Chung, Chungbuk National University
24
5. Gauss-Jordan
 A variation of Gauss
Elimination method
 The elimination step results
in an identity matrix
 Back substitution is not
necessary
 Gauss소거법보다 50% 이상
연산횟수 증가
Ex.9.12
Se-Woong Chung, Chungbuk National University
25
PART 3. Linear Algebraic Eq.
CHAPTER 10
LU Decomposition
Matrix Inversion
Se-Woong Chung, Chungbuk National University
26
1. LU Decomposition
AX  B  0
• Make upper diagonal matrix using Gauss Elimination method
u11 u12 u13   x1   d1 
0 u
  x   d 
u
 2
22
23   2 

 d 
 0
0 u33  
x
 3  3
U X  D  0
Transformed from B
• Assume a lower diagonal matrix with 1’s on the diagonal
1 0
L  l21 1
l31 l32
Elimination
factor=a21/a11
0
0
1
LU X  D  AX  B
LU  A
Ex. 10.1
LD B
Se-Woong Chung, Chungbuk National University
27
 The Steps in LU Decomposition
1) [A] is decomposed into lower [L] and upper [U] triangular
matrix
2) [L] {D} = {B} is used to generate {D} by forward substitution
3) [U] {X} = {D} is solved to determine {X} by back substitution
Se-Woong Chung, Chungbuk National University
28
 LU Decomposition Algorithm
Forward
sustitution
Save [L]
Backward
substitution
Se-Woong Chung, Chungbuk National University
29
 The Matrix Inverse using LU Decomposition
AA1  A1 A  I 
• 단위벡터를 우변 1열벡터로 놓고
• 전진대입과 후진대입으로
역행렬의 첫번째 열 계산
• 순차적으로 다음 열 계산
Ex. 10.3
Se-Woong Chung, Chungbuk National University
30
PART 3. Linear Algebraic Eq.
CHAPTER 11
Special Matrices and
Gauss Seidel
Se-Woong Chung, Chungbuk National University
31
1. Special Matrices
 Banded Matrix(띠행렬)
BW (Band width): 대역폭
HBW(Half-band width): 반-대역폭
BW = 2HBW + 1
i  j  HBW 이면, aij  0
행과 열의 차이가 반대역폭 보다 큰
행렬 요소(entry)의 값은 0
띠 밖의 “0”의 요소들에 대해서는
연산을 수행하지 않는 것이 효율적
Se-Woong Chung, Chungbuk National University
32
 Tridiagonal Matrix (삼중대각 시스템)
(BW = 3)
 f1
e
 2









g1
f2
e3
g2
f3

g3







en 1

f n 1
en
  x1   r1 
 x   r 
 2   2 
  x3   r3 
   
     
     
   
     
   

g n 1  xn 1  rn 1 

f n   xn   rn 
We have been changed our notation for the coefficients from a’s and
b’s to e’s, f’s, g’s to avoid storing large number of useless zeros
Se-Woong Chung, Chungbuk National University
33
 Thomas Algorithm for Tridiagonal System
Decomposition
Forward Substitution
LD B
Backward Substitution
U X  D  0
Ex. 11.1 : Practice by hand
Se-Woong Chung, Chungbuk National University
34
Thomas Algorithm : Excel-VBA
Sub Thomas(e, f, g, r, n, x)
Call Decomp(e, f, g, n)
Call Substitute(e, f, g, r, n, x)
End Sub
Sub Decomp(e, f, g, n)
Dim k As Integer
For k = 2 To n
e(k) = e(k) / f(k - 1)
f(k) = f(k) - e(k) * g(k - 1)
Next k
End Sub
Se-Woong Chung, Chungbuk National University
35
Thomas Algorithm : Excel-VBA
Sub Substitute(e, f, g, r, n, x)
Dim k As Integer
For k = 2 To n
r(k) = r(k) - e(k) * r(k - 1)
Next k
x(n) = r(n) / f(n)
For k = n - 1 To 1 Step -1
x(k) = (r(k) - g(k) * x(k + 1)) / f(k)
Next k
End Sub
Se-Woong Chung, Chungbuk National University
36
2. Gauss-Seidel
• Most commonly used iterative method
• Guess a value and then use a systematic method to
obtain a refined estimate of the root
AX   B
 a11 a12
a
 21 a22
a31 a32
b1  a12 x2  a13 x3
x1 
a11
a13   x1  b1 
   

a23   x2   b2 
 b 
a33  
x
3
   3
b2  a21 x1  a23 x3
x2 
a22
b3  a31 x1  a32 x2
x3 
a33
Se-Woong Chung, Chungbuk National University
37
Gauss-Seidel Iterative Procedure
1) Start the solution process by choosing guesses
for x1, x2, x3 (ex, 0)
2) Calculate a new x1 and use it (new x1) and x3
to calculate a new x2
3) Repeat the same process for x3 calculation
4) Until our solutions xi values converges
 a ,i
xij  xij 1

100 %   s
j
xi
Ex. 10.3 : By hand
Se-Woong Chung, Chungbuk National University
38
Graphical Depiction of Gauss-Seidel Procedure
(a) Gauss-Seidel method
(b) Jacobi iterative method
Se-Woong Chung, Chungbuk National University
39
Gauss-Seidel 반복법의 수렴조건
•
각각의 행에서 대각선 상의 요소 절대값이 비대각선 요소들
의 절대값 합보다 커야 수렴이 보장된다
n
aii   aij
j 1
j i
a11 x1  a12 x2  ............  a1n xn  b1
a21 x1  a22 x2  ............  a2 n xn  b2
.
.
an1 x1  an 2 x2  ............  ann xn  bn
Se-Woong Chung, Chungbuk National University
40
Relaxation method (이완법)
•
A slight modification of Gauss-Seidel method
to enhance convergence
•
새로운 값은 이전과 현재의 반복계산에서 얻은 값
의 가중평균값으로 사용
xinew  xinew  (1   ) xiold
•
0 <  < 1 : Underrelaxation (하이완법)
수렴이 잘 안되는 시스템을 수렴하도록 하거나, 진동을 감소
시켜 수렴을 빠르게 할 때 사용
•
1 <  < 2 : Overrelaxation (상이완법)
현재값에 더 큰 가중치를 주는 것이며, 수렴하고 있는 시스템
의 수렴성을 가속화하기 위해 사용
Se-Woong Chung, Chungbuk National University
41
Gauss-Seidel with Relaxation Method
Sub Gseid(a, b, n, x, imax, es, lambda)
Dim i As Integer, j As Integer, iter As Integer, sentinel As Integer
Dim dummy As Single, sum As Single, ea As Single, old As Single
For i = 1 To n
dummy = a(i, i)
‘각 행별 대각행렬요소 저장
For j = 1 To n
a(i, j) = a(i, j) / dummy
‘각요소를 대각행렬로 나눔
Next j
b(i) = b(i) / dummy
‘b 행렬요소를 대각행렬로 나눔
Next i
For i = 1 To n
sum = b(i)
For j = 1 To n
If i <> j Then sum = sum - a(i, j) * x(j) ‘새로운 x값 계산
Next j
x(i) = sum
Next i
iter = 1
Se-Woong Chung, Chungbuk National University
42
Gauss-Seidel with Relaxation Method
Do
sentinel = 1
‘x값들의 오차점검 ->반복횟수 줄임
For i = 1 To n
old = x(i)
sum = b(i)
For j = 1 To n
If i <> j Then sum = sum - a(i, j) * x(j)
Next j
x(i) = lambda * sum + (1# - lambda) * old
‘이완법 적용
If sentinel = 1 And x(i) <> 0 Then ‘x(i) 허용오차는 하나만 계산
ea = Abs((x(i) - old) / x(i)) * 100
‘오차 계산
If ea > es Then sentinel = 0 ‘허용오차 초과시 sentinel값 “0” 그 후 오차계산 안 함
End If
Next i
iter = iter + 1
‘xi 모두 허용오차 만족시 exit
If sentinel = 1 Or iter >= imax Then Exit Do
Loop
Se-Woong Chung, Chungbuk National University
43
3. Using Excel for Linear Algebraic Eqs.
1 1
2

1 2
 33
1
4
1  x
3   1  1.83333 
1   x   2.16667
2 2
3   x3  2.35000
5
Se-Woong Chung, Chungbuk National University
Ex. 11.4
44
PART 3. Linear Algebraic Eq.
CHAPTER 12
Engineering Applications
Se-Woong Chung, Chungbuk National University
45
1. Steady-state Analysis of a System of Reactors
Accumulation = Inputs - Outputs
dM d (C  V )

 M in  M out  Cin  Qin  Cout  Qout
dt
dt
Se-Woong Chung, Chungbuk National University
46
 A Completely Mixed CSTR
dM
 M in  M out  0
dt
At steady-state
d (C  V )
 Cin  Qin  Cout  Qout  0
dt
C1  Q1  C2  Q2  C3  Q3
Easy to calculate by hand
C3 = 18.6 mg/m3
Se-Woong Chung, Chungbuk National University
47
 Five Reactors Linked by Pipes
Steady state problem
Se-Woong Chung, Chungbuk National University
48
 Develop Mass Balance Equations
5(10)  C3  Q31  C1  Q12  C1  Q15
Reactor 1.
Similar equations can be developed for other reactors
6C1  C3  50
 3C1  3C2  0
 C2  9C3  160
 C2  8C3  11C4  2C5  0
 3C1  C2  4C5  0
11.51
11.51


{c}  19.06


17
.
00


11.51
Practice : Gauss Elimination, EXCEL Matrix, Gauss-Seidel
Se-Woong Chung, Chungbuk National University
49
Homework #4
8장: 8.8, 8.14
10장: 10.11
11장: 11.6, 11.10
(1) Develop mass balance equations for each problems
(2) Use Gauss Elimination method, Gauss-Seidel, Excel
package to solve the problems
(3) Hand out your report including (a) problem statement, (b)
your mass balance equations, (c) description of numerical
methods that you used, and (d) solutions
Se-Woong Chung, Chungbuk National University
50