Use the TDMA method to solve the given equation12 partialdif

Use the TDMA method to solve the given equation:(12%) partialdifferential^2T/partialdifferentialx^2 + partialdifferential^2T/partialdifferentialy^2 = 1 L_x = 10m, Ly = 10m Boundary conditions: T_x = 0 = 0, T_y = 0 = 0, partialdifferential T/partialdifferential x_x = L_x = 0, partialdifferential T/partialdifferential y_y = L_y = 0

Solution

The question is solved using MATLAB. Standard TDMA solver function is used

clear all
clc
%%

% as mentioned in the hint, n order to solve for T at x and y, we proceed
% in half steps i.e. first advance k+1/2 for x direction and then from
% k+1/2 to k in y direction. In each step, we formulate a tridiagonal
% matrix in x and y.

Lx=10;
Ly=10;
dx=0.1;
dy=0.1;
nx=Lx/dx; % number of points in x
ny=Ly/dy; % number of points in y

T=zeros(ny,nx); % initialize T matrix rows represent y axis and coloumn x axis

T(1,:)=0 ;     % first row of T is x=0 and boundary condition T=0 at x=0
T(:,1)=0; % T=0 at y=0
T(:,nx)=T(:,nx-1); % gradient in x-direction is 0
T(ny,:)=T(ny-1,:); % gradient in y-direction is 0

% Procedure
% We can follow two approaches to solve this via TDMA.
% Methos 1
% Here we olve along all x nodes for each constant j.Once we have solved for T
% along x direction, we change the direction of solution and then solve for
% all T in y direction.

% Method 2
% In this we solve for T in x along one row, then change the direction and
% solve for a coloumn, then again swap and so on. This is alternating
% technique.

% In this we will be using Method 1

% Also we use TDMA algorithm (and fucntion) which solves
%A_i*X_(i-1) + B_i*X_i + C_i*X_(i+1) = D_i (where A_1 = 0, C_n = 0)
% A,B,C,D are input vectors. X is the solution, also a vector.

for jj=2:ny-1
A(1:nx-3)=1/dx^2; % we are going to solve for x from i=2 to nx-1.
% so A will be for i=2:nx-2 or size of A is 1:nx-3
C(1:nx-3)=1/dx^2;
  
B(1:nx-2)=(-2/dx^2)-(2/dy^2);
D(1:nx-2)=1-((T(jj-1,2:nx-1)-T(jj+1,2:nx-1))/dy^2);
A=[0 A];
C=[C 0];
T(jj,2:nx-1)=TDMAsolver(A,B,C,D);
A=[];B=[];C=[];D=[];
end

% Boundary conditions
T(1,:)=0 ;     % first row of T is x=0 and boundary condition T=0 at x=0
T(:,1)=0; % T=0 at y=0
T(:,nx)=T(:,nx-1); % gradient in x-direction is 0
T(ny,:)=T(ny-1,:); % gradient in y-direction is 0

% solving along Y (step 2)
for ii=2:nx-1
A(1:ny-3)=1/dy^2; % we are going to solve for x from i=2 to nx-1.
% so A will be for i=2:nx-2 or size of A is 1:nx-3
C(1:ny-3)=1/dy^2;
  
B(1:ny-2)=(-2/dx^2)-(2/dy^2);
D(1:ny-2)=1-((T(2:ny-1,ii-1)-T(2:ny-1,ii+1))/dy^2);
A=[0 A];
C=[C 0];
T(jj,2:nx-1)=TDMAsolver(A,B,C,D);
A=[];B=[];C=[];D=[];
end

% Boundary conditions
T(1,:)=0 ;     % first row of T is x=0 and boundary condition T=0 at x=0
T(:,1)=0; % T=0 at y=0
T(:,nx)=T(:,nx-1); % gradient in x-direction is 0
T(ny,:)=T(ny-1,:); % gradient in y-direction is 0

TDMA solver:

function x = TDMAsolver(a,b,c,d)
%a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
n = length(b); % n is the number of rows

% Modify the first-row coefficients
c(1) = c(1) / b(1); % Division by zero risk.
d(1) = d(1) / b(1); % Division by zero would imply a singular matrix.

for i = 2:n-1
temp = b(i) - a(i) * c(i-1);
c(i) = c(i) / temp;
d(i) = (d(i) - a(i) * d(i-1))/temp;
end

d(n) = (d(n) - a(n) * d(n-1))/( b(n) - a(n) * c(n-1));

% Now back substitute.
x(n) = d(n);
for i = n-1:-1:1
x(i) = d(i) - c(i) * x(i + 1);
end
end