Crank Nicolson Method

The code needs debugging

clc
clear

MYU=1;
A=1;
N=100;
M=100;
LX=1;
LY=1;
DX=LX/M;
DY=LY/N;

%-----------------------------INITILIZATION--MATRIX------------------------
t=1;
for i=1:M;
    for j=1:N;
        AN(i,j,t)=(MYU*A/DY);
        AS(i,j,t)=(MYU*A/DY);
        AW(i,j,t)=(MYU*A/DX);
        AE(i,j,t)=(MYU*A/DX);
        SP(i,j,t)=100;
        AP(i,j,t)=AN(i,j,t)+AS(i,j,t)+AE(i,j,t)+AW(i,j,t)-SP(i,j,t);
    end
end

for i=1:M;
    d(i)=1;
end

for i=1:N-2;
    dd(i)=MYU*A/DY;
end

for i=1:N-2;
    ddd(i)=MYU*A/DY;
end

for i=1:M-1;
    dddd(i)=MYU*A/DX;
end

for i=1:M-1;
    ddddd(i)=MYU*A/DX;
end

A1=diag(dd,2)
A2=diag(ddd,-2)
A3=diag(dddd,1)
A4=diag(ddddd,-1)
A5=diag(d,0)
B(1)=12;
for i=2:M-1;
    B(i)=0;
end
B(M)=10
 A=A1+A2+A3+A4+A5
 AB=A^(-1);
 U=B*AB

% t=1;
% for i=3:M-2;
%     for j=3:N-2;
%         U(i,j,t)=0;
%     end
% end
%boundary conditions
U(:,1,t)=100;
U(:,2,t)=100;
U(1,:,t)=100;
U(2,:,t)=100;

U(M-1,:,t)=100;
U(M-2,:,t)=100;
U(:,N-1,t)=100;
U(:,N-2,t)=100;

image(U)
% t=1;
% a=-1;
% b=-1;
% for i=1:M;
%     a=a+1;
%     for j=1:N;
%         b=b+1;
%         x(i,j,t)=a*DX;
%         y(i,j,t)=b*DY;
%     end
% end
%
%
% t=1;
% for i=1:M;
%     for j=1:N;
%         plot(x(i,j,t),y(i,j,t),'*r');
%         hold on
%     end
% end
% grid on
% hold off
%
% figure(2)
% axis([ 0  1  0  1])
%  image(U)

%-----------------------------CALCULATION--MATRIX--------------------------
% t=1;
% for i=1:M;
%     for j=1:N;
%         AN(i,j,t)=(MYU*A/DY);
%         AS(i,j,t)=(MYU*A/DY);
%         AW(i,j,t)=(MYU*A/DX);
%         AE(i,j,t)=(MYU*A/DX);
%         SP(i,j,t)=1;
%         AP(i,j,t)=AN(i,j,t)+AS(i,j,t)+AE(i,j,t)+AW(i,j,t)-SP(i,j,t);
%     end
% end
%
% t=1;
% for i=3:M-2;
%     for j=3:N-2;
%         U(i,j,t)=(AN(i,j,t)*U(i,j+1,t)+AS(i,j,t)*U(i,j-1,t)+AE(i,j,t)*U(i+1,j,t)+AW(i,j,t)*U(i-1,j,t))/(AP(i,j,t));
%     end
% end
%    
% U(M-1,:,t)=0;
% U(M-2,:,t)=0;
% U(:,N-1,t)=0;
% U(:,N-2,t)=0;
% image(U(:,:,1))
% U
%
% %-----------------------------TIME--STEPPING------------------------------
% %--------------------------CRANK-NICKOLSON--------------------------------
%
% L=3;
% THETA=0.5;
% ROH=1;
% DT=0.1;
% B=10;
% for t=1:L
% for i=1:M;
%     for j=1:N;
%         AN(i,j,t)=(MYU*A/DY);
%         AS(i,j,t)=(MYU*A/DY);
%         AW(i,j,t)=(MYU*A/DX);
%         AE(i,j,t)=(MYU*A/DX);
%         SP(i,j,t)=1;
%         AP0(i,j,t)=ROH*(DX/DT);
%     end
% end
% end
%
% for t=1:L-1;
% for i=3:M-2;
%     for j=3:N-2;
%         AP(i,j,t+1)=THETA*(AW(i,j,t+1)+AE(i,j,t+1))+AP0(i,j,t)-0.5*SP(i,j,t);
%     end
% end
% end
%
% for t=1:L-1;
% for i=3:M-2;
%     for j=3:N-2;
%         U1(i,j,t)=AW(i,j,t)*(THETA*(U(i-1,j,t+1)+(1-THETA)*(U(i-1,j,t));
%         U2(i,j,t)=AE(i,j,t)*(THETA*(U(i+1,j,t+1)+(1-THETA)*(U(i+1,j,t));
%         U3(i,j,t)=(AP0(i,j,t)-(1-THETA)*AW(i,j,t)-(1-THETA)*AE(i,j,t))*U(i,j,t);
%         U(i,j,t+1)=(U1(i,j,t)+U2(i,j,t)+U3(i,j,t)+B)/AP(i,j,t+1);
%     end
% end
%end

---

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