How to plot 2D line graph to compare the approximate solution with the actual solution?

1 view (last 30 days)
I am trying to solve the following problem.
I wrote the following code to solve this problem. I stored my approximate solution in the matrix U of size Nx+2 times Ny+2. I want to compare my approximated solution with the actual solution. I am not sure how to plot a matrix in to 2D line graph. Can anyone please help me with this step?
%------Construction of meshgrids-------------------
Lx = 1;
Ly = 1;
%Ny=Nx;
Nx = 3; % Number of interior grid points in x direction
Ny = 3; % Number of interior grid points in y direction
dx = Lx/(Nx+1); % Spacing in the x direction
dy = Ly/(Ny+1); % Spacing in the y direction
x=0:dx:Lx;
y=0:dy:Ly;
[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';
% -------------------------------------------------------
Iint = 1:Nx+1; % Indices of interior point in x direction
Jint = 1:Ny+2; % Indices of interior point in y direction
Xint = X(Iint,Jint);
Yint = Y(Iint,Jint);
U = zeros(Nx+2,Ny+2); % Define U to store the answer
%--------------------------------------------------------
uinit = zeros(Nx+1,Ny+2);
u0 = @(x,y) cos(2*pi*x).*cos(2*pi*y);
U(1,:) = u0(x(1),y(Jint));
U(Nx+2,:)=U(1,:);
uinit(1,:)=U(1,:);
uinit(Nx,:)=uinit(1,:);
F1 = reshape(uinit, (Nx+1)*(Ny+2),1);
%---------------------------------------------
%assembly of the tridiagonal matrix(LHS)
sx = 1/(dx^2);
sy = 1/(dy^2);
e=ones(Nx,1);
A = zeros(Nx,Nx+2);
B = zeros(Nx+1,Nx+2);
T=spdiags([sx.*e,(-2*sx)+(-2*sy).*e,sx.*e],-1:1,Nx,Nx);
A(:,2:Nx+1) = T;
B(1:Nx,:)=A;
B(Nx+1,1)=-1;
B(Nx+1,2)=1;
B(Nx+1,Nx+1)=1;
B(Nx+1,Nx+2)=-1;
%T(1,Nx+1)= sx;
%T(Nx+1,1)= sx;
D = zeros(Nx+1,Nx+2);
D1 = zeros(Nx,Nx+2);
D2=spdiags(sy.*e,0,Nx,Nx);
D1(:,2:Nx+1)=D2;
D(1:Nx,:)=D1;
C=blktridiag(B,D,D,Ny+2);
for i=1:Nx
for j=1:Nx+1
C(i,j)=(1/2).*C(i,j);
C((Nx+1)*(Ny+1)+i,(Nx+2)*(Ny+1)+j)=(1/2).*C((Nx+1)*(Ny+1)+i,(Nx+2)*(Ny+1)+j);
end
end
%---------------------------------------------------------------
%----------Compuating the R.H.s term----------------------------
f = @(x,y) (-8*pi*pi).*cos(2*pi*x).*cos(2*pi*y);
%Solve the linear system
rhs = f(Xint,Yint);
rhs_new = zeros(Nx+1,Ny+2);
rhs_new(1:Nx,:)=rhs(2:Nx+1,:);
%for j=1:Ny+2
% for i=2:Nx+3
% rhs(i,j) = (-8*pi*pi)*cos(2*pi*x(i)).*cos(2*pi*y(j));
% end
%end
for i=1:Nx
rhs_new(i,1)=(1/2).*rhs_new(i,1);
rhs_new(i,Ny+2)=(1/2).*rhs_new(i,Ny+2);
end
%convert the rhs into column vector
F = reshape(rhs_new, (Nx+1)*(Ny+2),1);
F2 = F - sx.* F1;
uvec = C\F2;
v= reshape(uvec, Nx+2,Ny+2);
U(2:Nx+1,:)=v(2:Nx+1,:);
%----------Computation of Exact solution--------------------
ue = zeros(Nx+2,Ny+2);
ue(:,:) = cos(2*pi*X) .* cos(2*pi*Y); % Exact Solution
%Error
err = norm(U(:,:) - ue(:,:))/norm(ue(:,:));

Answers (0)

Categories

Find more on Vector Fields in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!