PoiseuilleFlow.m

% Author: Housam Binous

% Poiseuille Flow in a rectangular channel

% National Institute of Applied Sciences and Technology, Tunis, TUNISIA

% Email: [e-mail]

% LU decomposition leads to fill-ins, which make solution by 
% elimination not manageable

A=zeros(900,900);

dx=1/29.

for k=31:870
    if( mod(k,30)~=1 && mod(k,30)~=0)
        A(k,k)=4/dx^2;
    end
end

for k=31:870
    if( mod(k,30)==1 || mod(k,30)==0)
        A(k,k)=1;
    end
end

for k=1:30
        A(k,k)=1;
end

for k=871:900
        A(k,k)=1;
end

for k=31:870
    if( mod(k,30)~=1 && mod(k,30)~=0)
        A(k,k-1)=-1/dx^2;
        A(k,k+1)=-1/dx^2;
        A(k,k+30)=-1/dx^2;
        A(k,k-30)=-1/dx^2;
    end
end

[L,U,P]=lu(A);
NZ_A=nnz(A)
NZ_U=nnz(U)
NZ_L=nnz(L)
NZ_P=nnz(P)
figure(1);
spy(A)
figure(2);
spy(L)
figure(3);
spy(U)
figure(4);
spy(P)
b=5*ones(900,1);
for i=1:30
    b(i)=0;
end
for i=871:900
    b(i)=0;
end
for i=31:870
    if (mod(i,30)==0 || mod(i,30)==1) b(i)=0; end
end

% Numerical solution 

c=A\b;

for i=1:30
    for j=1:30
        d(i,j)=c(i+(j-1)*30);
    end
end
figure(5);
surf(d)

% Analytical solution 

figure(6)
for i=1:30
    x=(i-1)*dx;
    for j=1:30
        y=(j-1)*dx;
        da(i,j)=-5/2*(x - 1)*x;
        for n=1:200
            da(i,j)=da(i,j)+(-4*1^2*5*cosh((n*pi*(1 - 2*y))/(2*1))...
                *sech(1*n*pi/(2*1))*sin(n*pi/2)^2 ... 
                 *sin(n*pi*x/1))/(n^3*pi^3);
        end
    end
end
surf(da)

% Comparing analytical and numerical solutions

figure(7);
plot(d(:,13),'r')
hold on
plot(da(:,13),'b')