# 2D Heat equation - MatLab implementation (FD in space, Expl. Euler in time)

I'm trying to solve the heat equation in 2D in $$\Omega=[0,1] \times [0,1]$$, with homogeneous Dirichlet boundary conditions, and initial condition $$u(x,y,0)=\sin(2 \pi x y)$$ i.e.

$$\begin{cases} u_t=u_{xx}+u_{yy} \\ u(x,y,0)=\sin(2 \pi x y), \quad (x,y) \in \Omega \\ u_{|d\Omega}=0 \end{cases}$$

with finite differences in space and explicit euler in time.

I'd like to check if my implementation, i.e. updating and boundary conditions imposition is correct. From the following plot it seems "good", but I've never done it in 2D before, so every comment on the code is really appreciated. I discretize the unit square with $$Nx=N_y=N$$ points, and after initializing the numerical paramters $$\text{dx},\text{ dt}$$ (checking that stability condition is fulfilled), I write the finite difference code where my solution at time $$t_n$$ is stored in a $$N \times N$$ matrix.

In particular, I solve in the "internal points" (i.e from $$2$$ ti $$N-1$$), and then add the boundary conditions outside of the loops.

In order to impose an homogeneous Dirichlet condition as $$u(0,y,t)=0, y \in [0,1]$$, I proceed in this way :

$$u(0,y,t)=0, y \in [0,1]$$ means $$u(x_1=0,y_j,t_n)=0$$ for every time $$t_n$$ and for every $$j \in {1, \ldots, N}$$ (here $$x_1$$ means just the first grid point in $$x$$). But this means that I have to fill the first line with 0's !

I do the same argument for $$u(x,0,t)=0, x \in [0,1]$$ and get that

\begin{align} u(x_i,y_1=0,t_n)=0 \end{align} for every time $$t_n$$ and every $$i \in \{1,\ldots,N \}$$. Hence the first column of $$u$$ has to be filled with 0's.

EDIT

I have to do the argument before also for $$u(1,y,t)=0, y \in [0,1]$$ and $$u(x,1,t)=0, x \in [0,1]$$, and so I set to zeros also the last row and the last column!

clear all
close all

N=100; %steps in x and y
Nt=1000; %steps in time
x=linspace(0,1,N);
y=x;
dx=x(2)-x(1);
dy=y(2)-y(1);

[x2d,y2d] = meshgrid(0:dx:1, 0:dy:1);  % create grid

u0=@(x,y) sin(2*pi*x.*y);
u=zeros(N,N);
dt= min([dx,dy])^2/4;

c1=dt/(dx^2);
c2=dt/(dy^2);

for i=1:N
for j=1:N
u(i,j)=u0(x(i),y(j)); #build initial condition matrix
end
end

t=0;
for n=1:Nt

%update
for j=2:N-1
for i=2:N-1
uNew(i,j)=u(i,j)+c1*(u(i,j+1)-2*u(i,j)+u(i, j-1)) + c2*(u(i+1,j)-2*u(i,j)+u(i-1,j));
end
end

%boundary conditions
for i=1:N
uNew(i,1)=0;
uNew(i,N)=0;
end
for j=1:N
uNew(1,j)=0;
uNew(N,j)=0;
end
u=uNew;
t=t+dt;

pause(0.01)
surf(x,y,uNew)
xlabel('x')
ylabel('y')
title(sprintf('Heat equation at time t=%0.5f',t))

end

• OK... What's your question now? – Alone Programmer Oct 18 at 0:02
• Is the implementation right? – VoB Oct 18 at 0:19
• @VoB: Do you have an analytic solution against which you can compare? – Richard Oct 18 at 0:24
• Treatment of boundary conditions in the revised post looks right. – Maxim Umansky Oct 18 at 16:45
• We're not going to read over your code to debug it in detail. At first glance, it looks fine. But if you want to check if your implementation is correct, find an exact solution and compare to it. Given your boundary conditions, you should find that for an infinite time that the heat becomes 0, so that could also help check your implementation. Those would be the typical ways to try to validate something like this. – EMP Oct 18 at 19:07