# How to plot the graph for different value of t?

4 views (last 30 days)
kaps on 12 Jun 2021
Commented: kaps on 13 Jun 2021
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
kaps on 13 Jun 2021
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