How to plot the graph for different value of t?

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
 

[dpb]

Which/what style of graph do you want to plot; I've not tried to read the code to know for certain what it returns...

But, I'll note that in the last for...end loop that

 rhs = rhs(U,Iint, Jint, Nx, Ny ,sx , sy);

overwrites the rhs array the first pass through so nothing else will change after that point...

[kaps]

I want to plot 3D plot of the output. I observed this mistake. I replace this as

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

function rhs= rhs(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

[dpb]

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]

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]

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

This approach will allow you to visually inspect how the solution evolves at each time step directly in MATLAB.

Hope this helps.

Regards,

Nipun