clear
echo on
%Solution to problem 7
%
%In this problem we use PLU decomposition
%to solve the problem Ax=b.  We start by
%generating the 7x7 Hilbert matrix:
n=7;
a=zeros(n,n);
%test
for i=1:n
 for j=1:n
  a(i,j)=1/(i+j-1);
  echo off
 end
end
echo on
pause
%This yields the matrix:
a
pause
%Now we do the plu decomposition using the 'lu'
%command:
[l,u,pt]=lu(a);
p=pt';
%Remember that the way matlab sets up the plu
%decomposition yields a permutation matrix which
%is the transpose of the one discussed in class.
%We can look at the matrices:
p
pause
l
pause
u
pause
%And we see that they are of the correct form.
%Now let's solve the problem.  First we set up
%the vector b:
b= [7.00000000000000;5.28214285714286;4.34206349206349;3.71309523809524;...
    3.25382395382395;2.90061327561328;2.61919746919747]
pause
%Now for z.  Remember that p'*p=eye, thus:
z=p'*b
pause
%Next we use forward substitution.  We begin by initializing y:
y=zeros(size(z));
y(1)=z(1)/l(1,1);
for i=2:n
 y(i)=(z(i)-l(i,:)*y)/l(i,i);
 echo off
end
echo on
pause
%We can check this code by multiplying things back out:
z-l*y
%which results in zero (at least to the machine
%precision) as it should.  Note that we did not limit
%the array multiplication of l and y to the [1:i-1]
%elements of the ith row of l and column vector y as
%would be expected from the forward substitution
%formula.  We can get away with this because before
%the ith step is completed, the [i:n] elements of y
%are all zero.  If you don't initialize y with the
%zeros command, this doesn't work!
pause
%We finish by doing the back substitution part:
x=zeros(size(y));
x(n)=y(n)/u(n,n);
for i=(n-1):-1:1
 x(i)=(y(i)-u(i,:)*x)/u(i,i);
 echo off
end
echo on
pause
%and we check the result by multiplying it back out:
y-u*x
%which again is zero.
pause
%Finally, let's compare our result against the solution
%using the operator '\':
x
a\b
%which matches
echo off