## 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. 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

#### 1 comment:

1. solution is diverging