Code covered by the BSD License

# Finite Difference Method to solve Heat Diffusion Equation in Two Dimensions.

### Sathyanarayan Rao (view profile)

Heat diffusion equation of the form Ut=a(Uxx+Uyy) is solved numerically. All units are arbitrary.

heat_2d.m
```%------------------------------------------------------------------------
%--- Heat Equation in two dimensions-------------------------------------
%--- Solves Ut=alpha*(Uxx+Uyy)-------------------------------------------
%------------------------------------------------------------------------

clc;
close all;
clear all;

%--dimensions...........................................................

N = 51;
DX=0.1; % step size
DY=0.1;
Nx=5;
Ny=5;

X=0:DX:Nx;
Y=0:DY:Ny;

alpha=5; % arbitrary thermal diffusivity

%--boundary conditions----------------------------------------------------

U(1:N,1:N) = 0 ;

U(1,1:N) = 100;
U(N,1:N) = 0;
U(1:N,1) = 0;
U(1:N,N) = 0;

%--initial condition------------------------------------------------------

U(23:29,23:29)=1000; % a heated patch at the center

Umax=max(max(U));

%-------------------------------------------------------------------------

DT = DX^2/(2*alpha); % time step

M=2000; % maximum number of allowed iteration

%---finite difference scheme----------------------------------------------

fram=0;
Ncount=0;
loop=1;
while loop==1;
ERR=0;
U_old = U;
for i = 2:N-1
for j = 2:N-1
Residue=(DT*((U_old(i+1,j)-2*U_old(i,j)+U_old(i-1,j))/DX^2 ...
+ (U_old(i,j+1)-2*U_old(i,j)+U_old(i,j-1))/DY^2) ...
+ U_old(i,j))-U(i,j);
ERR=ERR+abs(Residue);
U(i,j)=U(i,j)+Residue;

end
end

if(ERR>=0.01*Umax)  % allowed error limit is 1% of maximum temperature
Ncount=Ncount+1;
if (mod(Ncount,50)==0) % displays movie frame every 50 time steps
fram=fram+1;
surf(U);
axis([1 N 1 N ])
h=gca;
get(h,'FontSize')
set(h,'FontSize',12)
colorbar('location','eastoutside','fontsize',12);
xlabel('X','fontSize',12);
ylabel('Y','fontSize',12);
title('Heat Diffusion','fontsize',12);
fh = figure(1);
set(fh, 'color', 'white');
F=getframe;
end

%--if solution do not converge in 2000 time steps------------------------

if(Ncount>M)
loop=0;
disp(['solution do not reach steady state in ',num2str(M),...
'time steps'])
end

%--if solution converges within 2000 time steps..........................

else
loop=0;
disp(['solution reaches steady state in ',num2str(Ncount) ,'time steps'])
end
end
%------------------------------------------------------------------------
%--display a movie of heat diffusion------------------------------------

movie(F,fram,1)

%------END---------------------------------------------------------------

```