Question

how can I change the Gauss-Seidel method to SOR method code in Matlab?

The question has shows that In implementing SOR method in MATLAB, one should not calculate Tw and cw by formulas Tw = (D -wL)^(-1)[(1-w)D+wU)] and Cw = w(D-wL)^(-1)b , where w stands for omega and using MATLAB's built-in inv function, since this function requires O(n^3) flops and therefore the
whole matter loses its point.

I have tried for many times but I can't get the correct answers. Please help Thanks :)

Execute code:

function x = GaussSeidel(A,b,x0,tol,kmax,output)
% This function returns an approximate solution of a
% linear system A*x=b, obtained by the Gauss-Seidel method.
%
% x0 is an initial approximation,
% tol is tolerance,
% kmax is the maximum number of iterations,
% output is a parameter which regulates displays:
% 0 - no display,
% 1 - some display,
% 2 - detailed display.
%
% Input: A, n by n matrix,
% b, n by 1 column,
% x0, initial approximation,
% tol, tolerance,
% kmax, maximum number of iterations,
% output, a parameter which regulates displays.
%
% Output: x, a solution.

[m,n]=size(A);
if m~=n,
error('A is not square');
end
m=length(b);
if m~=n,
error('Dimensions of A and B do not agree');
end
% Write A as D-L-U
[D,L,U] = DLU_decomposition(A);
% Calculate the matrix T and the vector c
T = L+U;
c=b;
for i=1:n,
T(i,:) = T(i,:)/D(i,i); % T = D^(-1)*(L+U)
c(i) = b(i)/D(i,i); % c = D^(-1)*b
end
% Calculate x
x=x0; % initial approximation
if output >= 2
disp(' k TOL x1 x2');
end
for k=1:kmax,
xprev=x;
for i=1:n
x(i)=c(i);
for j=1:n
if j~=i
x(i)=x(i)+T(i,j)*x(j);
end
end
end
TOL=norm(x-xprev,inf);
if output >= 2
out=[k, TOL, x(1), x(2)];
disp(out);
end
if TOL<tol,
if output >= 1
disp('The Gauss-Seidel method has converged.');
end
break
end
end
if output >= 1
if TOL>=tol
disp('The Gauss-Seidel method did not converge.');
end
s=sprintf('The norm of residual vector A*x-b is %e.',norm(A*x-b));
disp(s);
end
end

Answer 1

%%Matlab code for Gauss Siedel Method
clear all
close all

%A_matrix is the Coefficient Marix A,b_matrix is the Result Matrix b,
A_matrix=[3 -1 0 0 0 1/2; -1 3 -1 0 1/2 0 ;0 -1 3 -1 0 0; 0 0 -1 3 -1 0; 0 1/2 0 -1 3 -1; 1/2 0 0 0 -1 3];
b_matrix=[5/2;3/2;1;1;3/2;5/2];
conv_es=10^-12;
w=1.1;
%exact solution
x_exact=[1;1;1;1;1;1];
x0=[0 0 0 0 0 0];
%displaying the matrix
fprintf('The A matrix is \n')
disp(A_matrix)

fprintf('The b matrix is \n')
disp(b_matrix)

%result for Gauss Siedel method
[x_g,error_gauss]=Gauss_method(A_matrix,b_matrix,x0,conv_es);

%result for Gauss Jordan method
[x_sor,error_sor]=SOR(A_matrix,b_matrix,x0,w,conv_es);

%Printing the result
fprintf('\n\tThe result using Gauss Siedel Method after 7 iterations is\n');
disp(x_g)


%Printing the result
fprintf('\n\tThe result using SOR Method for w=1.1 after 7 iterations is\n');
disp(x_sor')

fprintf('Printing the error for each method\n')
fprintf('\tItr \t Gauss \t SOR\n')
%Printing the error

for i=1:length(error_sor)
     fprintf('\t%d\t%f\t%f\n',i,error_gauss(i),error_sor(i))
end

%Function for Gauss Siedel method
function [x,error]=Gauss_method(A,b,x0,conv)
%Jacobi method
ea=1;it=0;xx=ones(length(b),1);
while ea>=conv
    it=it+1;
    for i=1:length(b)
        s=0;
        for j=1:length(b)
            if i~=j
                s=s+A(i,j)*x0(j);
            end
          
        end
         x0(i)=(b(i)-s)/A(i,i);
    end
     if it==8
        break
     end
    if it==1
        ea=10;
    else
        ea = norm(x1 - x0', 'inf');
    end

    error(it) = abs(norm(xx - x0'));
    x1=x0;
end
x=x0';
end


%%Matlab function for Gauss Siedel SOR Method
function [x,error]=SOR(A,b,x0,lambda,conv)
%inputs are A=matrix A; b=result vector, x0=initial guess
%maxit=maximum iteration;conv=error convergence
%outputs are x=solution matrix, it=maximum iteration; lambda=w for sor
%Gauss Siedel method
ea=10;k=0; xx=ones(length(b),1);

while ea>=conv
    k=k+1;
    for i=1:length(b)
        s=0;
        for j=1:length(b)
            if i~=j
                s=s+A(i,j)*x0(j);
            end         
        end
        x11=((lambda*(b(i)-s))/A(i,i))+(1-lambda)*x0(i);
        x0(i)=x11;   
    end
    if k==1
        ea=10;
    else
        ea = norm(x1 - x0', 'inf');
    end
    if k==8
        break
    end
  
    x1=x0';
    error(k)=norm(xx - x1, 'inf');
end
x=x0;
end

%%%%%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%%%%%%%%