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

3 views (last 30 days)
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
  1 Comment
Eshna Sasha
Eshna Sasha on 5 Apr 2021
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);

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 6 Apr 2021
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.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!