%--------------------------------------------------------------------------
% 2D- Laplace Equation by Using Sparse Matrix
% Constraints: 1:Quadratic Mesh
% 2:Drichlet Boundary Condition(Known)
%--------------------------------------------------------------------------
clc; clear all; close all;
tic; % Time start
%---------------INPUTS-----------------------------------------------------
L = 2; % Lenght
H = 1; % Height
N = 21; % the grid size or number of nodes
% N = input('Number of Iteration, (N) =');
Tb1 = 0; % Initialize the boundary conditions
Tb2 = 10;
Tb3 = 0;
Tb4 = 0;
%--------------------------------------------------------------------------
Nx= N-2; %number of unknown points along x direction
Ny= N-2; %number of unknown points in the y direction
dx = L/(N-1); %the grid spacing
dy = H/(N-1);
x = [0:dx:L]; y=[H:-dy:0]; % x & y points in solution domain
[X,Y] = meshgrid(x,y);
Ae = dy/dx;
Aw = dy/dx;
As = dx/dy;
An = dx/dy;
Ap = Ae+An+Aw+As;
%--------------------------------------------------------------------------
%set the matrix A and the Right hand side vector
e = ones(Nx*Ny,1);
A = spdiags([-e*An,-e*Aw,e*Ap,-e*Ae,-e*As],[-Nx,-1,0,1,Nx],Nx*Ny,Nx*Ny);
b = linspace(0,0,Nx*Ny)'; % initialize b matrix
for i = 1:Nx
%border 1
j = 1;
k = j + (Nx)*(i-1);
b(k) = b(k) + Tb1*Aw;
if (k-1>0)
A(k,k-1) = 0;
end
%border 2
j = Ny;
k = j + (Nx)*(i-1);
b(k) = b(k) + Tb3*Ae;
if (k+1<Nx*Ny)
A(k,k+1) = 0;
end
end
for j = 1:Ny
%border 3
i = 1;
k = j + (Nx)*(i-1);
b(k) = b(k) + Tb2*An;
%border 4
i = Nx;
k = j + (Nx)*(i-1);
b(k) = b(k) + Tb4*As;
end
%--------------------------------------------------------------------------
% AFUL = full(A); % Full SParse Matrix
TT = A\b; % Unknown Temperature
toc; % time end
%-------------------------------------------------------------------------
S=0;
for i=1:Nx;
for j=1:1:Ny
TN(i,j) = TT(j+S);
end
S= S+Nx;
end
%--------------------------------------------------------------------------
TNN = zeros(N,N);
TNN(:,1) = Tb3;
TNN(:,N) = Tb4;
TNN(1,:) = Tb2;
TNN(N,:) = Tb4;
for i= 2:Nx+1
for j = 2:Ny+1
TNN(i,j) = TN(i-1,j-1);
end
end
% -------------------Analytical Solution-----------------------------------
R=200; % Number of iteration
theta = zeros(N);
for r = 1:R
tint = (2/pi)*(((-1)^(r+1)+1)/r)*(sin(r*pi*X/L).*(sinh(r*pi*Y/L)/sinh(r*pi*H/L)));
theta = theta + tint;
end
t = Tb1 + (Tb2-Tb1)*theta;
t(1,1:N-1) = Tb2;
t(N,1:N-1) = Tb4;
t(1:N,1) = Tb1;
t(1:N,N) = Tb3;
t = t;
% --------------------PLOTS------------------------------------------------
figure(1);
subplot(2,2,3)
surf(X,Y,TNN); %Numerical Solution
title('Temp. Distribution (Numerical Solution)');
xlabel('x');
ylabel('y');
zlabel('t(x,y)');
%--------------------------------------------------------------------------
% figure(2);
hold on
subplot(2,2,4)
surf(X,Y,t); %Analytical Solution
title('Temp. Distribution (Analytical Solution)');
xlabel('x');
ylabel('y');
zlabel('T(x,y)');
%--------------------------------------------------------------------------
% figure(2)
% isolignes
hold on
subplot(2,2,2)
contour(X,Y,TNN)
% figure(3);
hold on
subplot(2,2,1)
spy(A) % Sparse Matrix Patteren
%--------------------------------------------------------------------------