Transcript Solving Tridiagonal Matrix

Solving Tridiagonal System
( By LU decomposition and backward substitution )
Course: CS 4MN3
Student: Imam Abdukerim
Professor: Dr. Sanzheng Qiao
Date:
February 27, 2009
Tridiagonal System
 d1 u1
l
 1 d 2 u2

l 2 d 3 u3

 


ln  2

0' s




d n 1
ln 1
u = [u1, u2, …, un-1]
d = [d1, d2, …, dn-1, dn]
l = [l1, l2, …, ln-1]
0' s   x1 
 b1 
 x 
b 
 2 
 2 
  x3 
 b3 

  
  
  
bn 1 
un 1   xn 1 


 
d n   xn 
 bn 
LU decomposition of
Tridiagonal Matrix
1
 l1
 1 1

l12




0' s
 d1 u1
l
 1 d2

l2




0' s
0' s  d11
 
 
 
1

 
 
 
l1n  2
1
 
l1n 1 1   0' s
0' s 

u2


d 3 u3

  

ln  2 d n 1 u n 1 

ln 1 d n 
u11
d12
u12
d13
u13


d1n 1
0' s 





u1n 1 

d1n 
T homasAlgorithm:
l1n 1  ln 1 / d1n 1
d1n  d n  (ln 1 / d1n 1 ) * un 1
LU decomposition of
Tridiagonal Matrix (Matlab Code)
function [u1,d1,l1] = decomt(u,d,l)
n = length(d); % get the length of the diagonal
% initialize the output vectors
u1=u;
d1=d;
l1=l;
%Perform LU decomposition
d1(1)=d(1);
for i=2:n
l1(i-1) = l(i-1)/d1(i-1); % Update the lover triangle vector
d1(i) = d(i) - (l(i-1)/d1(i-1))*u(i-1); %Update the diaginal
end
Complexity: O(n)
LU decomposition of
Tridiagonal Matrix (Testing)
u = [23 43 22 3 77];
d = [3 34 55 18 13 56];
l = [1,2,3,4,5];
>>[u,d,l] = decomt(u,d,l)
u = 23 43 22 3 77
d = 3.0000 26.3333 51.7342 16.7242 12.2825 24.6545
l = 0.3333 0.0759 0.0580 0.2392 0.4071
U=
3.0000 23.0000
0
0
0
0
0 26.3333 43.0000
0
0
0
0
0 51.7342 22.0000
0
0
0
0
0 16.7242 3.0000
0
0
0
0
0 12.2825 77.0000
0
0
0
0
0 24.6545
(Construct L and U using decomposed vectors)
(Construct triangular matrix by multiplying L
and U)
>>L*U
L = diag(l,-1)+eye(6);
U = diag(d)+diag(u,+1);
ans =
L=
1.0000
0
0
0
0
0
0.3333 1.0000
0
0
0
0
0 0.0759 1.0000
0
0
0
0
0 0.0580 1.0000
0
0
0
0
0 0.2392 1.0000
0
0
0
0
0 0.4071 1.0000
3.0000 23.0000
0
0
0
0
1.0000 34.0000 43.0000
0
0
0
0 2.0000 55.0000 22.0000
0
0
0
0 3.0000 18.0000 3.0000
0
0
0
0 4.0000 13.0000 77.0000
0
0
0
0 5.0000 56.0000
Solving Tridiagonal System
TriangularSystem : LUx  b
Step1 : Ly  b
1
 l1
 1 1

l12 1

 


l1n  2

0' s
Step2 : Ux  y
1
l1n 1
0' s   y1   b1 
  y  b 
  2   2 
  y3   b3 

 

 
   
  yn 1  bn 1 
 
  
1   yn   bn 
d11







 0' s
u11
d12
u12
d13 u13


d1n 1
0' s   x1   y1 
 x  y 
  2   2 
  x3   y3 




 
   
u1n 1   xn 1   yn 1 
 
 

d1n   xn   yn 
Solving Tridiagonal System
(Matlab Code - Ditailed)
function [x] = solvet(u,d,l,b)
n = length(d);
x = (1:n);
y = (1:n);
%Solve Tridiagonal system LUx=b;
% Step 1 : Solve Ly=b for y
y(1) = b(1);
for i=2:n
y(i) = b(i) - l(i-1)*y(i-1);
end
% Step 2 : Solve Ux=y for x
x(n) = y(n)/d(n);
for i=(n-1):-1:1
x(i) = (y(i)-u(i)*x(i+1))/d(i);
End
Solving Tridiagonal System
(Testing)
(Using solvet() function)
(Using Matlab standard function)
u = 23
M=
3.0000 23.0000
0
0
0
0
1.0000 34.0000 43.0000
0
0
0
0 2.0000 55.0000 22.0000
0
0
0
0 3.0000 18.0000 3.0000
0
0
0
0 4.0000 13.0000 77.0000
0
0
0
0 5.0000 56.0000
43
22
3
77
d = 3.0000 26.3333 51.7342 16.7242 12.2825 24.6545
l = 0.3333
b = 45
0.0759
43
33
67
0.0580
89
0.2392
0.4071
76
>> solvet(u,d,l,b)
b = 45
y = 45.0000 28.0000 30.8734 65.2097 73.4036 46.1186
>> M\b1
x = -11.9303
ans =
-11.9303
3.5127
-1.5000
4.9307
-5.7506
1.8706
3.5127 -1.5000
ans = -11.9303
4.9307 -5.7506
3.5127 -1.5000
1.8706
4.9307 -5.7506
1.8706
43
33
67
89
76
Solving Tridiagonal System
(Matlab Code – overwrite b)
function [b] = solvet(u,d,l,b)
n = length(d);
%Solve Tridiagonal system LUx=b;
% Step 1 : Solve Ly=b for y (overwrite b)
for i=2:n
b(i) = b(i) - l(i-1)*b(i-1);
end
% Step 2 : Solve Ux=y for x (overwrite b)
b(n) = b(n)/d(n);
for i=(n-1):-1:1
b(i) = (b(i)-u(i)*b(i+1))/d(i);
End
Question(s)?