clear
echo on
format compact
% In this problem we examine the behavior of a dilute 
% suspension of particles which is being sheared.  As the 
% suspension is sheared (the motion which is produced 
% when a fluid is confined between two concentric cylinders 
% and one of them is rotated) the particles will tumble over 
% one another.  This will lead to a random walk of the 
% particles which can be characterized by a diffusion 
% coefficient, much like a molecular diffusivity.  We can 
% measure this quantity - the shear-induced self-diffusivity - 
% by examining the random walk of a single tracer particle in 
% a suspension of otherwise identical particles.  Because the 
% random motion is due to the interaction of the tracer with 
% other particles, the diffusivity will be identically zero if the 
% concentration of the other particles is zero.
pause
% In our lab we measured the shear-induced self-diffusivity in 
% the dilute limit.  If the particles are perfect spheres, theory 
% suggests that the leading order term in the diffusivity 
% should go as c^2 where c is the concentration.  Other 
% models suggest that if the spheres are not perfect, the 
% diffusivity should be proportional to c to leading order.  More
% recently, models have suggested that even for smooth spheres the
% presence of bounding walls gives rise to an O(c) diffusivity.
% 
% In this example we examine data to determine which model best
% describes the experiments:
% 
% c       diffusivity     error
% 0.01    2.37e-4         2.2e-5
% 0.025   5.63e-4         6.5e-5
% 0.05    1.42e-3         1.6e-4
% 0.075   2.27e-3         4.0e-4
% 0.10    4.06e-3         7.6e-4
% 0.15    9.96e-3         1.8e-3
% 
% The errors given above are the one standard deviation errors
% in the measured diffusivities calculated from the statistics
% governing the measurement process.
pause
% Using this data, we wish to fit a constitutive equation for the 
% diffusivity of the form:
% 
%         Diffusivity = c * x(1) + c^2 * x(2) + c^3 * x(3)
% 
% using weighted linear regression, and determine the error 
% in each of the fitting coefficients.  In particular, we need to
% determine if the difference between the O(c) coefficient x(1) and 
% zero is statistically significant.

% First we put in the data:
c=[.01,.025,.05,.075,.1,.15]';
diff=[2.37e-4,5.63e-4,1.4153e-3,2.27e-3,4.055e-3,9.965e-3]';
sigma=[2.2e-5,6.491228e-5, 1.617643e-4,3.98208e-4,7.605929e-4,1.775587e-3]';
pause
%Let's plot this data up:
figure(1)
errorbar (c,diff,sigma)
set(gca,'FontSize',18)
xlabel('concentration')
ylabel('diffusivity')
grid on
%As you can see, it certainly goes to zero as the
%concentration goes to zero.
pause
%The asymptotic behavior can be visualized better
%by plotting diff/c rather than just diff:
figure(2)
errorbar(c,diff./c,sigma./c,sigma./c,'o')
set(gca,'FontSize',18)
xlabel('concentration')
ylabel('diffusivity / concentration')
grid on
hold on
%Here there is clearly a non-zero asymptotic value
%of diff/c as c goes to zero.  This is what we are
%trying to capture.
pause
%OK, now for the weighted linear regression.  The
%unweighted linear regression problem is of the 
%form:
%
%   min || A x - b ||
%
%In the weighted regression problem we weight each
%data point by the inverse of its variance when forming
%the sum of the squares of the deviation.  Let's do this:
%
%First we set up the matrix of modeling functions.  We
%have a linear term, a quadratic term, and a cubic term:
a=[c,c.^2,c.^3]
pause
%This is the weighting function:
vardiff=diag(sigma.^2);
weight=inv(vardiff)
pause
%We now define a matrix kw:
kw=inv(a'*weight*a)*a'*weight
%which can be used the conventional way:
x=kw*diff
%which has a lead coefficient different from zero.  
pause
%Now for the error:
varx=kw*vardiff*kw'
xerror=diag(varx).^.5;
%So the coefficients and errors are:
[x,xerror]
%and in particular,
x(1)/xerror(1)
%So the order c coefficient is statistically different from
%zero (about 8 standard deviations).
pause
%Let's finish off with a bit of graphics.  Let's plot up
%the model predictions together with its uncertainty.
cp=[.001:.001:.15]';
ap=[cp,cp.^2,cp.^3];
dp=ap*x;
perror=diag(ap*varx*ap').^.5;
figure(2)
plot(cp,dp./cp,'b')
plot(cp,(dp+perror)./cp,'g')
plot(cp,(dp-perror)./cp,'g')
axis([0 .15 0 .1]);
%where we have plotted the 1-sigma confidence interval
%of the model function, based on the given standard deviations.
hold off
echo off