How to plot the graph for different value of t?

4 views (last 30 days)
kaps
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
  4 Comments
kaps
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

Sign in to comment.

Answers (0)

Categories

Find more on Language Fundamentals in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!