2D Heat equation Crank Nicolson method

Question

Ouvrir dans MATLAB Online

0 votes

Hi everyone

I'm trying to code te 2D heat equation using the crank nicolson method on

with test solution

and Dirichlet boundary conditions.

Here is a reference for this method: https://www.iist.ac.in/sites/default/files/people/parabolic.pdf

It envolves solving a system of equations

for which i'm pretty sure programmed it correctly including the boundary conditions.

The solution its wrong. Can someone help me?

Thank you in advance :)

%define variables
xmin = 0;
xmax = 1;
ymin= 0;
ymax= 1;
tmin= 0;
tmax= 1;
tsteps=100; %nº of time steps
N=10;%nº of space steps (1D)
%discretise the domain
dx = (xmax-xmin)/(N-1); %space step
x = xmin : dx : xmax ; %x grid
dy = (ymax-ymin)/(N-1); %space step
y = ymin : dy : ymax ; %y grid
dt = (tmax-tmin)/(tsteps-1); %time step
t = tmin : dt : tmax ; %time grid
%matrix3s to calculate implicit part
A=zeros(N^2,N^2);
B=zeros(N^2,N^2);
c=(dt/(dx^2));
 
for i=1:N^2
    if i<=N||i>=N^2-N
         A(i,i)=1;
         B(i,i)=1;
    elseif i/N - floor(i/N) == 0 || (i-1)/N - floor((i-1)/N)==0 
         A(i,i)=1;
         B(i,i)=1;
        
    elseif (i-2)/N - floor((i-2)/N) == 0 
         A(i,i)=1-2*c;
         A(i,i+1)=c/2;
         B(i,i)=1+2*c;
         B(i,i+1)=-c/2;
         A(i,i+N)=c/2;
         B(i,i+N)=-c/2;
         A(i+N,i)=c/2;
         B(i+N,i)=-c/2;
    elseif (i+1)/N - floor((i+1)/N) == 0
        A(i,i)=1-2*c;
        A(i,i-1)=c/2;
        B(i,i)=1+2*c;
        B(i,i-1)=-c/2;
        A(i,i+N)=c/2;
        B(i,i+N)=-c/2;
        A(i+N,i)=c/2;
        B(i+N,i)=-c/2;
    else
         A(i,i)=1-2*c;
         A(i,i+1)=c/2;
         A(i,i-1)=c/2;
         B(i,i)=1+2*c;
         B(i,i+1)=-c/2;
         B(i,i-1)=-c/2;
         A(i,i+N)=c/2;
         B(i,i+N)=-c/2;
         A(i+N,i)=c/2;
         B(i+N,i)=-c/2;
    end    
end    
for i=N^2-2*N : N^2-N
 
   A(i,i+N)=0;
   B(i,i+N)=0; 
   A(i+N,i)=0;      
  B(i+N,i)=0;
    
    
    
end      
%matrix to store the 2D initial condition
XY=zeros(N,N);
for j=2:N-1
    for i=2:N-1
    XY(j,i)=u(x(j),y(i),0);
    end
end
Unrecognized function or variable 'u'.
step1=reshape(XY,[],1);
mesh(XY)
hold off
pause(1)
F=zeros(N,N); % F vector
for k=2:tsteps-1
    
   
    for j=2:N-1
        for i=2:N-1
            F(j,i)=f(x(j),y(i),t(k)+dt/2);  %update F vecotr
        end
    end
       
    step2= B\A*(step1) + B\reshape(F,[],1)*dt; %solve system of equations
    step1=step2;
 
 
   mesh(reshape(step2,[N N]))
   title(num2str(k/tsteps))
   hold off
   
   
   pause(0.01)
   
  
end  