How to plot the graph for different value of t?

Question

kaps 2021-6-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/854700-how-to-plot-the-graph-for-different-value-of-t

回答： Nipun 2024-5-16

I wrote 2D heat equation backward finite difference code. How do I plot the graph for each t?

clc
% Construction of the meshgrid
Lx = 1;
Ly = 1;
Ny = 200;         % Number of grid points in y direction
Nx = 200;         % Number of grid points in x direction
T = 1;
Nt = 3;
dx = Lx/(Nx+1); % Spacing in x direction
dy = Ly/(Ny+1); % Spacing in y direction
dt = T/(Nt+1); % Spacing for the t 
x=0:dx:Lx;
y=0:dy:Ly;
t=0:dt:T;
alpha = 1/(2*pi*pi);
 
[X,Y] = meshgrid(x,y);%2d array of x,y values
X = X' ;             % X(i,j), Y(i,j) are coordinates of (i,j) point
Y = Y';
sx = (alpha*dt)/dx^2;
sy = (alpha*dt)/dy^2;
%-----------------------------------------------------------
Iint = 2:Nx+1; % Indices of interior point in x direction
Jint = 2:Ny+1; % Indices of interior point in y direction
 
Xint = X(Iint,Jint);
Yint = Y(Iint,Jint);
 
%--------Intial condition
u0 = sin(pi*X).*cos(pi*Y);
%uexact = sin(pi*X).*cos(pi*Y).*exp(-1);
 
%------adjust the rhs to include boundary terms
rhs = rhs(u0,Iint, Jint, Nx, Ny ,sx , sy);
 
%-----------------------------------------------
 
rhs;
%convert the rhs into column vector
F = reshape(rhs, Nx*Ny,1);
 
%------assembly of the tridiagonal matrix(LHS)---------
e=ones(Nx,1);
T=spdiags([-sx*e,(2.*sx+2.*sy+1)*e,-sx*e],-1:1,Nx,Nx);
D=spdiags(-sy*e,0,Nx,Nx);
A=blktridiag(T,D,D,Ny);
%-------------------------------------------------
U=zeros(Nx+2,Ny+2,length(t));
for k=1:length(t)
    U(:,:,k)=u0;
end
 
 
for k=2:length(t)
    
    uvec = A\F;
    U(Iint, Jint,k)= reshape(uvec, Nx,Ny);
    rhs = rhs(U,Iint, Jint, Nx, Ny ,sx , sy);
    F = reshape(rhs, Nx*Ny,1);
    
end
 

4 个评论
显示 2更早的评论隐藏 2更早的评论

dpb 2021-6-13

That won't fix your problems -- you can't have both a function and a variable of the same name--and I hadn't noticed you had already wiped out the rhs array you started with way early on with the first reassignment.

You need to define that array once and then refer to it elsewhere and not reassign to it anywhere, at least with out subscripting to modify, if needed, certain pieces of it.

Then, if you're going to use a function to do that, name it something different than the variable.

kaps 2021-6-13

在 MATLAB Online 中打开

In the actual code, I defined rhs_new function file, then I call it for different U(:,:,k) then I stored this value in rhs variable. Is there any mistake in the following block of code?

rhs = rhs_new(U(:,:,k),Iint, Jint, Nx, Ny ,sx , sy); % calling

%definition 
function rhs= rhs_new(u0,Iint, Jint, Nx, Ny ,sx , sy)
rhs = u0(Iint, Jint);
rhs(:,1) = rhs(:,1) + sy.*u0(Iint,1);
rhs(:,Ny)= rhs(:,Ny) + sy.*u0(Iint,Ny+2);
rhs(1,:) = rhs(1,:) + sx.*u0(1,Jint);
rhs(Nx,:) = rhs(Nx,:) + sx.*u0(Nx+2,Jint);
end

请先登录，再进行评论。

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

Answer 1

Nipun 2024-5-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/854700-how-to-plot-the-graph-for-different-value-of-t#answer_1458656

在 MATLAB Online 中打开

Hi Kaps,

I understand that you are trying to solve the 2D heat equation using the backward finite difference method and want to plot the graph for each time step (t). To achieve this, you can use MATLAB's surf or mesh function within your time-stepping loop to plot the temperature distribution at each time step. I'll guide you on how to integrate the plotting part into your code.

First, ensure you have a function rhs defined somewhere in your code or script that calculates the right-hand side of your equation. This function is crucial for your code to work but is not shown in your snippet.

Then, to plot the graph for each (t), you can add a plotting command right after updating your solution U for each time step. I recommend the following code snippet to achieve the desired result:

clc; clear; % It's good practice to clear the workspace too
% Your existing code setup here...
U=zeros(Nx+2,Ny+2,length(t));
for k=1:length(t)
    U(:,:,k)=u0;
end
figure; % Open a new figure window
for k=2:length(t)
    
    uvec = A\F;
    U(Iint, Jint,k)= reshape(uvec, Nx,Ny);
    
    % Update rhs for the next time step
    rhs = rhs(U(:,:,k),Iint, Jint, Nx, Ny ,sx , sy); % Make sure your rhs function can handle U(:,:,k)
    F = reshape(rhs, Nx*Ny,1);
    
    % Plotting the current state
    surf(X, Y, U(:,:,k)); % You can use mesh(X, Y, U(:,:,k)) if you prefer
    title(['Solution at t = ', num2str(t(k))]);
    xlabel('X');
    ylabel('Y');
    zlabel('Temperature');
    drawnow; % This updates the figure window dynamically
    
    pause(0.1); % Pause for a short while to visually inspect the plot. Adjust as needed.
end