Am I implementing the correct code for the ∂2T/∂x2 + ∂2T/∂y2 = 0 ? The slider on the left of the Matlab command window keep vibrating vigorously for over 10 minutes

I am given

and is tasked to get the solution for comparison to the steady state solution ∂T/∂t = ∂2T/∂x2 + ∂2T/∂y2 . The initial condition can be taken to be T=0.0 at t= 0 for the whole domain

Here's my code:

L=1;
nx=20;
ny=20;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error = 9e9;
tol = 1e-4;
T_L = 0;
T_T = 1;
T_R = 0;
T_B = 0;
T=ones(nx,ny);
 
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx, 2: ny-1) = T_B;
 
T(1,1) = (T_L + T_T)/2;
T(nx,ny) = (T_B+ T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
 
[X,Y] = meshgrid(x,y);
Told= T;
for i_s = 2
tic
if i_s == 2
gs_iterr = 1;
while(error> tol)
       for i = 2:nx-1;
            for j = 2:ny-1;
              T(i,j) = 0.25*((T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1)))
       end
end
error = max(max(abs(Told-T)));
gs_iterr = gs_iterr+1;
time_taken = toc;
end
end
 
figure(1)
[M,N] = contour(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprint('gauss sidel no. of iterations = %d, time taken = %f' , gs_iterr, time_taken);
title(title_text)
xlabel('X');
ylabel('Y');
fprintf('no. of iterations done are %d \ n', gs_iterr);
 
end

[Eshna Sasha]

Here's another code that I copied somewhere. Does it answer to the question? I need interpretation. Which 2 codes should I follow? Sorry, I'm a newbie.

Lx=1; Ly=1;
Nx=10; Ny=10;
nx=Nx+1 ; ny=Ny+1;
dx=Lx/Nx; dy=Ly/Ny;
x=(0:Nx)*dx; y=(0:Ny)*dy;
boundary_index= [ 1:nx, 1:nx:1+(ny-1)*nx,
                               1+(ny-1)*nx:nx*ny, nx:nx:nx*ny ];

diagonals= [4*ones(nx*ny,1), -ones(nx*ny,4)];
A=spdiags(diagonals, [0 -1 1 -nx nx], nx*ny, nx*ny);
 I=speye(nx*ny);
A(boundary_index,:)= I(boundary_index,:);
b=zeros(nx,ny);
b(:,1)=0;
b(1,:)=0;
b(:,ny)=1;
 b(nx,:)=0;
b= reshape(b,nx*ny,1);
Phi= A\b;
Phi=reshape(Phi,nx,ny);
[X,Y]= meshgrid(x,y);
v=[0.8 0.6 0.4 0.2 0.1 0.05 0.01];
contour(X,Y,Phi',v,'ShowText','on');
axis equal;
set(gca,'YTick', [0 0.2 0.4 0.6 0.8 1]);
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1]);
xlabel('$x$','Interpreter', 'latex', 'FontSize', 14);
ylabel('$y$','Interpreter','latex', 'FontSize', 14);
title('Solution of 1b' ,'Interpreter','latex','FontSize',16);

Answer 1

[Alan Stevens]

The following modifications to your code work:

L=1;
nx=20;
ny=20;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error = 1;
tol = 1e-4;
T_L = 0;
T_T = 1;
T_R = 0;
T_B = 0;
T=ones(nx,ny);
 
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx, 2: ny-1) = T_B;
 
T(1,1) = (T_L + T_T)/2;
T(nx,ny) = (T_B+ T_R)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
 
[X,Y] = meshgrid(x,y);
Told= T;
maxits = 500;    %%%%%%%%%%%%%%%%% Safer to include maximum allowed iterations
gs_iterr = 1;
while(error> tol)&&(gs_iterr<maxits)
       for i = 2:nx-1
            for j = 2:ny-1
              T(i,j) = 0.25*((T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1)));
            end
       end    
error = max(max(abs(Told-T)));
Told = T; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Update Told
gs_iterr = gs_iterr+1;
end

 
figure(1)
[M,N] = contour(X,Y,T);
clabel(M,N);
colormap(jet);
colorbar
title_text = sprintf('gauss sidel no. of iterations = %d', gs_iterr);
title(title_text)
xlabel('X');
ylabel('Y');
fprintf('no. of iterations done are %d', gs_iterr);

Note that because of the way the rows are numbered in Matlab , it looks like the high temperature is at the bottom.

[Eshna Sasha]

Thank you so much. It worked.