clear
echo on
%Solution to problem 6
%This script uses the attached function
%glyvisc.m and the data table (in ascii
%format) gly_table.dat to give an inter-
%polated viscosity as a function of 
%concentration for aqueous glycerin solutions
%The viscosity is evaluated for the 
%concentrations of 30%, 50% and 95%, and
%then we play around with the function a bit.
%First, we do the assigned part.
%
%At 30% concentration the viscosity is:
viscosity=glyvisc(.3)
%At 50% concentration the viscosity is:
viscosity=glyvisc(.5)
%and at 95% concentration the viscosity is:
viscosity=glyvisc(.95)
pause

%OK, now let's see how good this interpolation
%really is.  Let's plot the interpolated concen-
%trations over the range of [0:.002:1]:
conc=[0:.002:1];
plot(conc,glyvisc(conc))
ylabel('viscosity')
xlabel('wt% glycerin')
pause
%As can be seen, the viscosity increases -very-
%rapidly with concentration at high concentrations.
%We can see this better with a semilog plot:
semilogy(conc,glyvisc(conc))
ylabel('viscosity')
xlabel('wt% glycerin')
pause
%We can also put in the data points from the table:
load 'gly_table.dat';
dataconc=gly_table(:,1);
datavisc=gly_table(:,2)*0.01002;
hold on
semilogy(dataconc,datavisc,'o')
pause
%Note that the interpolated fit goes through all
%of the data points, but that the dependence on
%concentration isn't all that smooth.  This is
%probably due to measurement errors in the original
%data rather than any errors in the interpolation.
hold off
echo off



____________________________



function y=glyvisc(c)
%This function calculates the viscosity of an
%aqueous glycerin solution by parabolic inter-
%polation of a table.  We pick the three elements
%of the table which are closest to the input 
%concentration c.  We adapt the program so that
%we can take an array of values of c as input
%arguments.
%
%First we put in the arrays of concentrations
%and viscosities from the table.  It is easier
%to store these as an ascii file under the name
%'gly_table.dat'.  Thus the data will come in using
%the array name gly_table.

load 'gly_table.dat';
conc=gly_table(:,1);
visc=log(gly_table(:,2));
ntable=length(conc);

%This last bit puts the data into the two arrays
%conc and visc, and then determines the number of
%points that were put in.  It is written in a more
%general way, so that if the user wishes more data
%can be added to the viscosity table without having
%to rewrite the program.  Note that we are actually
%operating on the log of the viscosity rather than
%the viscosity itself.  This is because of the strong
%concentration dependence of the viscosity.

%We can also sort the data, to be sure that the 
%concentrations are all in order.  You can figure
%out how to use this command by invoking help sort:

[conc,order]=sort(conc);
visc=visc(order);
%This insures that all of the concentrations are
%in consecutive order in concentration, even if the
%concentrations in the table are all mixed up - such
%as might arise if new data points are added to the
%bottom of the file.

%OK, now for the meat of the program:
npt=length(c);
y=zeros(size(c));

for kk=1:npt
 ct=c(kk);
 
 %We find the closest data point in the array conc:
 [z,index]=min(abs(conc-ct));
 
 %Now we set up the array of three values we will use
 %in fitting a parabola to the data:
 ifit=[index-1:index+1]';
 
 %We also have to check to make sure that all of these
 %values are actually in the table:
 if index == 1
  ifit=ifit+1;
 end;
 if index==ntable
  ifit=ifit-1;
 end;
 
 %Now we can set up the matrix a:
 a=[ones(3,1),conc(ifit),conc(ifit).^2];
 %Note that this only works if ifit is a column array
 %rather than a row array, e.g. a (3,1) vector.
 
 %Finally we set up the right hand side:
 b=visc(ifit);
 
 %and then we get the coefficients:
 coeff=a\b;
 
 %and then calculate the viscosity at the desired
 %input concentration:
 y(kk)=coeff(1)+coeff(2)*ct+coeff(3)*ct^2;

%and then we loop back to the next concentration:
end;

%Finally, we convert from the log of the viscosity
%to the actual viscosity and multiply by the viscosity
%of water:
y=exp(y)*.01002;


_________________________________


Save the next bit under the name gly_table.dat


_________________________________


        .0      1.0
        .05     1.125
        .10     1.288
        .20     1.734
        .32     2.632
        .40     3.646
        .44     4.434
        .48     5.402
        .52     6.653
        .56     8.332
        .60     10.66
        .64     13.63
        .68     18.42
        .72     27.57
        .80     59.78
        .84     84.17
        .88     147.2
        .92     383.7
        .96     778.9
        1.00    1759.6