Strange Behavior of Matlab

The following exercise function solve linear system Ax=b with J.O.R. method:

function [x]=milspojor(A,b,x0,om,tol,itmax)

%Input-Check/Inizializzazione

[n,m]=size(A);
if n~=m
   error('Error:Non-Squared Matrix!');
end
if min(abs(eig(A)))==0
   error('Error:Singular Matrix!');
end
if isrow(b)==1
   b=b';
end
if isrow(x0)==1
   x0=x0';
end

it=0;r=b-A*x0;err0=norm(r,2); 
err=err0;x=x0;

%Calcolo

while err > tol && it < itmax
    it = it + 1;xit=zeros(n,1);
    for i=1:n
        xit(i)=om*...
               (b(i)-(A(i,1:i-1)*x(1:i-1)+A(i,i+1:n)*x(i+1:n)))/A(i,i)+...
               (1-om)*x(i);       
    end
    x=xit;r=b-A*x;err=norm(r,2)/err0;
end

if it  < itmax
   fprintf('\nConvergence-Reached @at n=% d',it);
end
if it >= itmax
   error('Error:Convergence-Not-Reached!');
end

end

Using as an Example: A=[-3 3 -6;-4 7 -8;5 7 -9] b=[-6 -5 3] (so exact solution is [1 1 1]) tol=eps, itmax=500;om=1 (Jacobi method) and x0=rand(3,1)=[0.8147 0.9058 0.1270]

the function works and the output is : Convergence-Reached @at n= 181 x =[1.0000 1.0000 1.0000]

BUT If I change the order of the scalar division (/A(i,i)) in the for loop in this way:

xit(i)=(om/A(i,i))*...
               (b(i)-(A(i,1:i-1)*x(1:i-1)+A(i,i+1:n)*x(i+1:n)))+...
               (1-om)*x(i);

(This is a fully licensed and apparently insignificant operation!!)

I get an unexpected result: the method don't converge and i get the Error Message described in the last if cicle.

I tought that this is caused by approximation error but I try to execute a single cicle of the for loop and in each way i get the same result.

So help me.... why this strange behavior?

Answer 1

[OCDER]

Nothing's wrong with the equation, but it's an approximation error caused when you're dividing/multiplying numbers with many decimal numbers that cannot be represented by the computer. I ran both versions of the equation and had the code stop when the two equations started to return different xit (it's not final solution). EX:

xit_1 =

   2.150133809189486
  -0.102098536507900
   0.422401406963743

xit_2 =

   2.150133809189486
  -0.102098536507900
   0.422401406963744

Doing computation with floating point numbers can be tricky, and so you have to account for rounding error (it's an issue in other programming languages too - not just MATLAB). More info about this from MathWorks here: https://www.mathworks.com/help/matlab/matlab_prog/floating-point-numbers.html#f2-98690

To fix your code, increase the tolerance tol. Ex:

while err > tol + 3*eps && it < itmax
   ....
end

Answer 2

[Leonardo Gherlinzoni]

Ok, so it's an error due to floating point operations like I expected. I don't find it immediately because I only try to compare the two version of the function with x0 without considering the propagation of the rounding error while the loop goes on. Anyway I make some other test and (like you suggest) the problems rise up when tolerance is shut down over 2*eps ( very restricting condition), before this limits (3eps 4eps exc..) the two method converge and differ for some iteration.

Thanks to mach for the answer!