%%%%% ADI Method for 2D heat equation %%%%%
clc
clear; close all;
uinitial=@(t,x,y) sin(pi*x)*sin(pi*y);
a = 0; b=1; c=0; d=1;
dt = 0.001; tfinal = 0.06;
h = 0.05;
n=1/h;
m=n;
h1 = h*h;
x=a:h:b; y=c:h:d;
% ----- Initial condition -----
t = 0;
for i=1:m+1,
for j=1:m+1,
u1(i,j) = uinitial(t,x(i),y(j));
end
end
% ----- Start Loop -----
k_t = fix(tfinal/dt); % Round toward zero
for k=1:k_t
t1 = t + dt; t2 = t + dt/2;
% ----- x-direction -----
for i=1:m+1 % Boundary condition.
u2(i,1) = uinitial(t2,x(i),y(1));
u2(i,n+1) = uinitial(t2,x(i),y(n+1));
u2(1,i) = uinitial(t2,x(1),y(i));
u2(m+1,i) = uinitial(t2,x(m+1),y(i));
end
for j = 2:n % internal mesh point for fixed y(j)
A = sparse(m-1,m-1); b=zeros(m-1,1);
for i=2:m
b(i-1) = (u1(i,j-1) -2*u1(i,j) + u1(i,j+1))/h1 + 2*u1(i,j)/dt;
if i == 2
b(i-1) = b(i-1) + uinitial(t2,x(i-1),y(j))/h1;
A(i-1,i) = -1/h1;
else
if i==m
b(i-1) = b(i-1) + uinitial(t2,x(i+1),y(j))/h1;
A(i-1,i-2) = -1/h1;
else
A(i-1,i) = -1/h1;
A(i-1,i-2) = -1/h1;
end
end
A(i-1,i-1) = 2/dt + 2/h1;
end
ut = A\b; % Solve the diagonal matrix.
for i=1:m-1,
u2(i+1,j) = ut(i);
end
end % Finish x
% ----- loop in y -direction -----
for i=1:m+1, % Boundary condition
u1(i,1) = uinitial(t1,x(i),y(1));
u1(i,n+1) = uinitial(t1,x(i),y(m+1));
u1(1,i) = uinitial(t1,x(1),y(i));
u1(m+1,i) = uinitial(t1,x(m+1),y(i));
end
for i = 2:m % internal mesh point for fixed x(i)
A = sparse(m-1,m-1); b=zeros(m-1,1);
for j=2:n,
b(j-1) = (u2(i-1,j) -2*u2(i,j) + u2(i+1,j))/h1 + 2*u2(i,j)/dt;
if j == 2
b(j-1) = b(j-1) + uinitial(t1,x(i),y(j-1))/h1;
A(j-1,j) = -1/h1;
else
if j==n
b(j-1) = b(j-1) + uinitial(t1,x(i),y(j+1))/h1;
A(j-1,j-2) = -1/h1;
else
A(j-1,j) = -1/h1;
A(j-1,j-2) = -1/h1;
end
end
A(j-1,j-1) = 2/dt + 2/h1; % Solve the system
end
ut = A\b;
for j=1:n-1,
u1(i,j+1) = ut(j);
end
end % Finish y
t = t + dt;
% ----- finish ADI for one time step -----
end
% ----- Comparison -----
uexact=@(t,x,y) exp(-2*pi*pi*t)*sin(pi*x)*sin(pi*y);
for i=1:m+1,
for j=1:m+1,
ue(i,j) = uexact(tfinal,x(i),y(j));
end
end
e = max(max(abs(u1-ue))) % Calculate error term
mesh(u1); % Plot the numerical solution
figure(2); mesh(u1-ue) % Plot error
