CODEWITHSUNDEEP
matlabmatrixlinear-algebranumerical-methods
Sofia
Sofia

Reputation: 456

Gauss-Seidel Method exceeds Machine Numbers?

Happy New Year everyone! :)

I'm writing a Gauss-Seidel function in Matlab and I'm encountering some problems.

The iteration must stop when we reach 6 decimal digits of precision. It means that the infinite norm (asked to use it) of x-xprevious must be less than 0.5*10^(-6). Firstly, here's my function: 

function [x] = ex1_3(A,b)

format long 

sizeA=size(A,1);

x=zeros(sizeA,1);  

%Just a check for the conditions of the Gauss-Seidel Method
for i=1:sizeA
    sum=0;
    for j=1:sizeA
        if i~=j 
            sum=sum+A(i,j);
        end
    end
    if A(i,i)<sum
        fprintf('\nGauss-Seidel''s conditions not met!\n');
        return     
    end
end

%Actual Gauss-Seidel Method

max_temp=10^(-6); %Pass first iteration
while max_temp>(0.5*10^(-6))
    xprevious=x;
    for i=1:sizeA
        x(i,1)=b(i,1);
        for j=1:sizeA
            if i~=j
              x(i,1)=x(i,1)-A(i,j)*x(j,1); 
            end
        end
        x(i,1)=x(i,1)/A(i,i);
    end
    x
    %Calculating infinite norm of vector x-xprevious 

    temp=x-xprevious;
    max_temp=temp(1,1);
    for i=2:sizeA
       if abs(temp(i,1))>max_temp
           max_temp=abs(temp(i,1));
       end
    end
end

And now the problems! When I call the function for a 3x3 array, I think it works. However, when I call it for a 10x10 array x becomes Inf (I guess it's out of machine numbers limits). Is there anything I can do to prevent this, except for changing the infinite norm and the 6 decimal digits precision (I must use these two, because my tutor told me so) ?

In the array I use (which was given to me) the entries outside the diagonal are -1 and the ones on the diagonal are 3. b is like this b=[2;1;1;1;1;1;1;1;1;2] (for n=10)

Upvotes: 3

Views: 899

Answers (1)

lennon310
lennon310

Reputation: 12689

Your condition of the Gauss-Seidel Method is not correct:

D=diag(diag(A));
L=-tril(A,-1);U=-triu(A,1);
B=(D-L)\U;
R = max(abs(eig(B)));
if R>=1
    fprintf('\nGauss-Seidel''s conditions not met!\n');
    return     
end

R is called spectral radius of iterative matrix B. It has to be less than 1 that Gauss-Seidel converges. Actually the matrix A in your test case has the R=1.8092, thus Gauss-Seidel method won't converge.

Check this slide from page 18 for more details.

EDIT

According to @LutzL's comment, you may use Gershgorin circle theorem to estimate the eigenvalue rather than calculate them with computational cost.

Upvotes: 1

Related Questions