Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

A Modified Thomas Algorithm by MATLAB Codes

Modified Thomas Algorithm: For special matrices such as tridiagonal matrix, the Thomas algorithm may be applied. The given matrix in the question is not in tri-diagonal format. So, in the following program, the matrix is made tridiagonal by taking coefficients of the upper and lower triangles to the right side of the equation and then the algorithm is implemented. The initial guesses for the solutions are assumed which is corrected iteratively in the program. This means the initial guesses are substituted to the right hand side variables and then by Thomas algorithm new solutions are obtained which then again are considered as initial values based on the no of iterations and error limit. The process is simply iterative. The matrix in the following is not convergent by Thomas algorithm.


System of Algebraic Equations











A MATLAB Program for Modified Thomas Algorithm

% Modified Thomas Algorithm to solve system of equations iteratively
 
a=[1 -1 4 0 2 9; 0 5 -2 7 8 4; 1 0 5 7 3 -2; 6 -1 2 3 0 8; -4 2 0 5 -5 3; 0 7 -1 5 4 -2]; % Given Matrix
 
b=[19; 2; 13; -7; -9; 2]; % System input
 
a1=[1 -1 0 0 0 0; 0 5 -2 0 0 0; 0 0 5 7 0 0; 0 0 2 3 0 0; 0 0 0 5 -5 3; 0 0 0 0 4 -2]; % Making original 
                                 % matrix into tridiagonal form by taking coefficients to the right side
 
c=1;
 
N=4; % No of iterations
 
e=[0.1; 0.1; 0.1; 0.1; 0.1; 0.1]; % Error tolerance
 
x(1,1)=0; x(2,1)=0; x(3,1)=0; x(4,1)=0; x(5,1)=0; x(6,1)=0; % Assuming initial solutions to continue 
                    % algorithm as the right hand side of equation contains variables with coefficient
 
 
while(c<N)
    
    x0=[x(1,1); x(2,1); x(3,1); x(4,1); x(5,1); x(6,1)];         % Initial solution
 
b1=[19-(4*x(3,1))-(2*x(5,1))-(9*x(6,1));                % Modified input matrix 
    2-(7*x(4,1))-(8*x(5,1))-(4*x(6,1));
    13-x(1,1)-(3*x(5,1))+(2*x(6,1));
    -7-(6*x(1,1))+x(2,1)-(8*x(6,1));
    -9+(4*x(1,1))-(2*x(2,1));
    2-(7*x(2,1))+x(3,1)-(5*x(4,1))];
 
%disp(b1);
 
for i=2:6           % Decomposition
    
a1(i,i-1)= a1(i,i-1)/a1(i-1,i-1);
 
a1(i,i)= a1(i,i)-(a1(i,i-1)*a1(i-1,i));
 
end
 
%a1(k,k-1)= a1(k,k-1)/a1(k-1,k-1);    
 
for i=2:6           % Forward Substitution
    
    b1(i,1) = b1(i,1)-(a1(i,i-1)*b1(i-1,1));
    
end
 
%disp(a1);
 
%disp(b1);
 
x(6,1)=b1(6,1)/a1(6,6);
 
%disp(x(6,1));
 
for i=5:-1:1        % Back Substitution
    
    x(i,1)=(b1(i,1)-(a1(i,i+1)*x(i+1,1)))/a1(i,i);
    
end
 
disp(x);
 
if ((abs(x(1:6,1))-abs(x0(1:6,1))) < e(1:6,1))  % Error checking
    
    break;
    
end
 
disp((abs(x(1:6,1))-abs(x0(1:6,1))));
 
x0=x;
 
disp(x0);
 
c=c+1;                    % loop continuation if condition is not satisfied
 
end                     % Program end

 
Program Outputs:
 

   1.0e+05 *% absolute error

   0.150253955555556
   0.087403955555556
   0.196764888888889
   0.140929777777778
   3.545057444444467
   7.086479888888927



   1.0e+05 * % solutions

  -0.150799955555556
  -0.087759955555556
  -0.197644888888889
   0.141539777777778
  -3.548047444444467
  -7.092449888888927 

No comments:

Post a Comment