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

4 次查看（过去 30 天）
Eshna Sasha
Eshna Sasha 2021-4-5
评论： Eshna Sasha 2021-4-6
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 个评论
Eshna Sasha
Eshna Sasha 2021-4-5
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);

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by