CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   TDMA code for MATLAB (http://www.cfd-online.com/Forums/main/97367-tdma-code-matlab.html)

amin144 February 15, 2012 14:34

TDMA code for MATLAB
 
TDMA is a quick Aligorithm for solving AX=b when A is Tridiagonal matrix:

Code:

% Written by Amin ShariatKhah 2012 - Shahrood University Of technology
%Feel Free to use it
%This code solve the AX=b (When A is Tridiagonal )

function X=TDMAsolver(A,b)

m=length(b);                % m is the number of rows
X=zeros(m,1);
A(1,2)= A(1,2)  ./ A(1,1);    % Division by zero risk.
b(1)=  b(1)    ./ A(1,1);    % Division by zero would imply a singular matrix

for i=2:m-1
    temp=  A(i,i) - A(i,i-1) .* A(i-1,i);
    A(i,i+1)=  A(i,i+1)  ./ temp;
    b(i)= ( b(i) - A(i,i-1) .* b(i-1) )  ./ temp;
end

i=m;
X(m)=(b(i) - A(i,i-1) .* b(i-1))  ./ (A(i,i) - A(i,i-1) .* A(i-1,i));

for i=m-1:-1:1
X(i)=  -A(i,i+1) .* X(i+1) + b(i);
end
end


there is another version of this code in wikipedia:

Code:

function x = TDMAsolver(a,b,c,d)
%a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
n = length(b); % n is the number of rows
 
% Modify the first-row coefficients
c(1) = c(1) / b(1);    % Division by zero risk.
d(1) = d(1) / b(1);    % Division by zero would imply a singular matrix.
 
for i = 2:n-1
    temp = b(i) - a(i) * c(i-1);
    c(i) = c(i) / temp;
    d(i) = (d(i) - a(i) * d(i-1)) / temp;
end
 
d(n) = (d(n) - a(n) * d(n-1))/( b(n) - a(n) * c(n-1));
 
% Now back substitute.
x(n) = d(n);
for i = n-1:-1:1
    x(i) = d(i) - c(i) * x(i + 1);
end
end


you can find some explanation about TDMA and find this code in C and Fortran90 in :
HTML Code:

http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

hityangsir March 27, 2012 11:33

thank you very much,that is what I need
 
thank you very much,that is what I need

amin144 April 6, 2012 15:27

You're welcome


All times are GMT -4. The time now is 05:49.