clear
echo on
%In this example we demonstrate the error
%induced by an ill-conditioned matrix.  We
%use as our example the hilbert matrix h
%given below:
n=5;
h=zeros(n,n);
%This just initializes the size of h

for i=1:n
 for j=1:n
  h(i,j)=1/(i+j-1);
 end;
end;
pause
%Now let's display this matrix:
h
pause
%As can be seen, the matrix looks perfectly
%ordinary - it does not reveal the nastiness
%lurking within.  For this example, we set
%up the solution first!  We shall figure out
%the value of b such that the solution is 
%x(i)=i.  So:
pause

x=[1:n]'
%That's the fast way to do it.  Now for b:
pause

b=h*x
%which looks like a perfectly innocent rhs.
pause

%Now let's solve the system using LU decomposition:
format long
x=h\b
format short
pause
%Note the errors in the solution!  This is in spite
%of double precision.  We can see this more clearly
%if we add a bit of error to the vector b:

bp=b+rand(n,1)*10^(-n);
%which adds a random error on the order of 10^-n to
%the vector b.  Now let's look at the solution:
pause
xp=h\bp
%which is wildly different from the original solution!
%Finally, we can compare the norms (we use the manhattan
%norm here):
pause
deltax=norm(xp-x,1)/norm(x,1)
deltab=norm(bp-b,1)/norm(b,1)
%and finally the ratio:
ratio=deltax/deltab
pause
%This ratio is pretty big!  We can compare this to
%the condition number of the hilbert matrix:
cond(h)
%which is of comparable magnitude.
pause
%We can do the same thing for errors in the matrix
%itself, except this time things are even worse.
%Let's add a little noise to the matrix h:
e=rand(n,n)/10^n;
hp=h+e
%Note that the error doesn't even show up when you
%print out the matrix to five decimal places!
pause
%We may calculate the new value of the solution
%vector:
xpp=hp\b
%Which is even further off from the original solution!
%The moral of the story is to watch out for ill-
%conditioned matrices!
pause
%A final note: the command cond(A) in Matlab
%calculates the condition number of the matrix A
%using the 2-norm.  The condition number using the
%1-norm (what was discussed in class) can be obtained
%using condest(A).  For more information about these
%commands, type 'help cond' and 'help condest' at the
%matlab prompt.

echo off