Help finding error in gauss seidel method

I have an error which is that for some set of matrices my solutions for the equations instead of converge to a number they diverge... Everything seems rights thou, I would be very grateful if you could help me. I have put some sample equations, the one which is with no comments works, the others dont...

Thank you very much

%Inputs: %- A - N*N Matrix - Inputed by the user %- B - N*1 Matrix - Inputed by the user %- eps - Level of precision of the solution of the system of equation % %Outputs: %- X - N*1 Matrix - Column vector which is the answer of the system of eqns % %

function [x]=gaussseidel(A,B)

eps=1e-10; %This is our convergence factor, when the error between our current and our previous solutions is less than this value, we have found the solution

%A=[2,-5;-3,7]; %B=[3;22];

%A=[1,1,-1;2,-3,4;-3,-1,-2]; %B=[-2;23;-15];

%A=[2,-5;-3,7]; %B=[31;-45];

A=[26,-22,-5;-20,17,4;5,-4,-1]; B=[-10;0;3];

[n,m]=size(A); [k,l]=size(B);

%Check for errors in the input of the matrices

if n~=k error('Dimensions of matrices A and B must agree') elseif n~=m error('Matrix A must be square') end

%Initialisations

x=zeros(m,1); %We make a guess for the solution, all x's are zero z=zeros(m,1); errorx=1; %We also want to know what is the error between our current solution and the previous one to see if its converging iteration=0; %Lastly we want to keep track of the number of iterations for two reasons: Give that data to the user and check the convergance of the system (if after a nº of iterations the error its bigger than a value then it doesnt converge)

   while abs(errorx)>eps

       iteration=iteration+1;

       p=x;% we save the values we have for x in a new variable so we can then compare the new values to the old ones to see convergence

       for i=1:m

           j = 1:m;
           j(i) = [];% we eliminate the value of the column associated with x for its iteration, i.e. if we want to calculate x1, we dont want to take into account the value assocciated with x1 in the calculation
           z=x;
           z(i) = [];% we surpress also the actual value of x when computing its iteration for each value of i, i.e. to calculate x(1) we dont want to involve x(1) in its own estimation

           x(i)=(B(i)-sum(A(i,j)*z))/A(i,i);

       end

     errorx=max(abs(x(1)-p(1)));

     if iteration==100 %We want to check at a given iteration if our solution is going to converge so we pick that point to be 15
         if errorx>10 %If the error at this point is greater than that
             error('The solution for these set of equations diverges')
         end
     end

   end
  iteration

end

Answer 1

[Alberto]

I think it has to fail int the 'while' loop: you dont have a prior value for the variable 'errorx', so it will probably will say that variable doesnt exists, and program will stop.

Iteration was not initialized. There may be other errors.

Some theorical stuff: there are conditions for the convergence, and the matrix you use doesn't verifies. Try this one:

 A=[4,1,-1;2,7,4;-3,-1,20]

This one is diagonal dominant, so the method should converge.

[Jacob]

So what other errors should i be checking for when my function encounters a matrix? What are the differences between the matrices that work and does that dont?

[Alberto]

I did some fast corrections, but at least it runs (properly i hope). Check it and compare.

%Inputs: %- A - N*N Matrix - Inputed by the user %- B - N*1 Matrix - Inputed by the user %- eps - Level of precision of the solution of the system of equation % %Outputs: %- X - N*1 Matrix - Column vector which is the answer of the system of eqns % %

function [x]=gaussseidel(A,B)

eps=1e-10; %This is our convergence factor, when the error between our current and our previous solutions is less than this value, we have found the solution errorx=1; iteration=0; %A=[2,-5;-3,7]; B=[3;22];

%A=[1,1,-1;2,-3,4;-3,-1,-2]; B=[-2;23;-15];

%A=[2,-5;-3,7]; %B=[31;-45];

% A=[26,-22,-5;-20,17,4;5,-4,-1]; B=[-10;0;3];

[n,m]=size(A); [k,l]=size(B);

%Check for errors in the input of the matrices

if n~=k error('Dimensions of matrices A and B must agree') elseif n~=m error('Matrix A must be square') end

%Initialisations

x=zeros(m,1); %We make a guess for the solution, all x's are zero z=zeros(m,1); errorx=1; %We also want to know what is the error between our current solution and the previous one to see if its converging iteration=0; %Lastly we want to keep track of the number of iterations for two reasons: Give that data to the user and check the convergance of the system (if after a nº of iterations the error its bigger than a value then it doesnt converge)

   while abs(errorx)>eps

       iteration=iteration+1

       p=x;% we save the values we have for x in a new variable so we can then compare the new values to the old ones to see convergence

       for i=1:m

           j = 1:m
           j(i) = [];% we eliminate the value of the column associated with x for its iteration, i.e. if we want to calculate x1, we dont want to take into account the value assocciated with x1 in the calculation
           z=x
           z(i) = []% we surpress also the actual value of x when computing its iteration for each value of i, i.e. to calculate x(1) we dont want to involve x(1) in its own estimation
            z
           x(i)=(B(i)-sum(A(i,j)*z))/A(i,i);

       end

     errorx=max(abs(x(1)-p(1)))

     if iteration==5 %We want to check at a given iteration if our solution is going to converge so we pick that point to be 15
         if errorx>10 %If the error at this point is greater than that
             error('The solution for these set of equations diverges')
         end
     end

   end
  iteration

end