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
solution is diverging
ReplyDelete