Gauss Seidel Matlab

clc 
clear 

% Constants Input 
L=1; 
N=32; 
alpha=32*(10)^(-6); 
dt=0.1;
dx=L/N;
r=(0.5*alpha*dt)/(dx)^2;
umax=10; 

% velocity Profile p=0;
 for pp=-N/2:N/2; 
p=p+1;
 x(p)=pp*dx; 
u(p)=umax*(1-(x(p))^2); 
end 


% checking paramters 
% plot(u) 
% x' 
% u'
pp=0; 
j=1;
for p=-N/2:N/2; pp=pp+1; if (p==-N/2) pp; dd(pp,j)=(1-2*r)*u(pp)+r*u(pp+1); 
elseif ((p>-N/2) & (p<N/2) ) 
pp;
dd(pp,j)=(1-2*r)*u(pp)+r*(u(pp+1)+u(pp-1)); 
elseif (p==N/2) 
pp; 
dd(pp,j)=(1-2*r)*u(pp)+r*u(pp-1); 
end 
end 

% checking parameters 
% plot(dd)
 % dd 
% p' 
% pp' 
% Assigning values to the diagonal Matrix

j=1;
for i=1:N+1;
a(i)=-r; 
b(i)=1-2*r;
c(i)=-r; 
d(i)=dd(i,j);
end 

% checking paramters 
% a'
 % b' 
% c' 
% d' 
% plot(d) 
% Calculating the first velocity output after Initialization

[u]= TDMAsolver(a,b,c,d)'; 


% Assigning the calculated velocity output into a 2D matrix
 j=1;
 for i=1:N+1;
 U(i,j)=u(i); 
 end 
 % Checking the output 
 plot(U(1:N+1,j)) 
 hold on
 grid on
 title('Velocity Profiles'); 
 xlabel('Xdirection');
 ylabel('Velocity'); 
 % Time Stepping Section
 j=0; 
 for time=1:200;
 j=j+1;
 for i=1:N+1; 
 dd(i,j+1)=U(i,j);
 d(i)=dd(i,j+1); 
 end
 [u]= TDMAsolver(a,b,c,d)';
 for i=1:N+1; 
 U(i,j+1)=u(i); 
 end 
 plot(U(1:N+1,j+1)) 
 hold on
 end

TDMA Function

function x = TDMAsolver(a,b,c,d) 
%a, b, c, and d are the column vectors for the compressed tridiagonal 
matrix 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 
 id = 1 / (b(i) - c(i-1) * a(i));
 % Division by zero risk. 
 c(i) = c(i)* id; 
 % Last value calculated is redundant. 
 d(i) = (d(i) - d(i-1) * a(i)) * id; 
end
 % Now back substitute. 
x(n) = d(n); 
for i = n-1:-1:1 
 x(i) = d(i) - c(i) * x(i + 1); 
end 
end

---

Unless otherwise noted, all content on this site is @Copyright by Ahmed Al Makky 2012-2015 - http://cfd2012.com