clear
echo on
%Solution to problem 10
%
%In this problem we use QR factorization
%to develop a correlation for the viscosity
%of liquid sodium.  This script uses the
%attached ascii table sodium_table.dat which
%contains the viscosity data.  First we
%read in the data:
load 'sodium_table.dat';
t=sodium_table(:,1);
visc=sodium_table(:,2);
pause
%We can plot up the viscosity:
figure(1)
plot(t,visc,'o')
xlabel('temperature')
ylabel('viscosity')
pause
%Notice that the data is almost hyperbolic in
%character.  This makes it difficult to fit with
%a polynomial.  Thus (as was suggested) we look
%at the inverse of the viscosity:
plot(t,1./visc,'o')
xlabel('temperature')
ylabel('1/viscosity')
%Note that this looks much closer to polynomial
%form!
%
pause
%OK, now for part a:
%First we set up the matrix.  We want to fit the
%data to a parabola, so:
m=length(t);
a=[ones(m,1),t,t.^2]
pause
%The normal equations are just A'Ax=A'b, so:
format short e            %This gives us the answer in
x=(a'*a)\(a'*(1./visc))   %scientific notation!
pause
%We can compare this to the data points:
hold on
plot(t,a*x)
%Notice that the residual is fairly small.
pause
%We can also compare the viscosity directly:
figure(2)
plot(t,visc,'o')
xlabel('temperature')
ylabel('viscosity')
hold on
plot(t,1./(a*x))
hold off
%which is also a pretty good fit.
%
%Now for part b:
%First we do the qr factorization:
[q,r]=qr(a)
pause
%Next we set up the new right hand side:
b=q'*(1./visc)
pause
%and finally we use back substitution:
n=3;
xqr(n)=b(n)/r(n,n);
for i=(n-1):-1:1
 temp=b(i);
 for j=(i+1):n
  temp=temp-r(i,j)*x(j);
 end
 xqr(i)=temp/r(i,i);
end
%which gives the same answer:
xqr
x'
echo off



______________________________________

save this data under the name sodium_table.dat



100         .705
200         .450
300         .345
400         .284
500         .234
600         .210
700         .186
800         .165
900         .150