clear
echo on
%Solution to Problem 2
%
%In this problem you were asked to calculate the
%error in the dimensionless drag on a sphere
%given that there was some known error in each
%of the quantities that went into the calculation.
%The problem is to determine the error in f(x)
%where x is a vector of parameters, each with
%known independent error.  The error in the
%calculated quantity may be determined by taking
%the partial derivative with respect to each of
%the parameters in x.  We can do this analytically
%or numerically.  We will do it numerically first.
pause
x=[.679,1.22,.051,2.40,1.215,980];
%These are the parameter values, with x(1) the velocity,
%x(2) the viscosity, x(3) the radius, x(4) the sphere
%density, x(5) the fluid density, and x(6) gravity.
pause
sigma=[.012,.03,.004,.05,.005,0]/2;
%These are the corresponding standard deviations.  Note
%that we divided by 2 because the error was specified as
%the 95% confidence level.  In the absence of further
%information we assume that all of these measurements
%are independent.  The matrix of covariance of these
%parameters is thus just a diagonal matrix.
pause
%Now we calculate the dimensionless drag (or velocity, if
%you want to think of it that way):
drag=x(1)/((x(4)-x(5))*x(6)*x(3)^2/x(2))
%which isn't all that close to the theoretical value
%of 2/9 = 0.2222...
pause
%The error can be obtained similarly.  We first take the
%derivatives.  This can be done analytically or numerically.
%For this problem it is probably easier to do them analytically,
%but we will do it numerically just to show you how it is done.
ep=x/10^8;
%(This works fine provided that none of the x values are zero)
n=length(x);
deriv=zeros(1,n);
for i=1:n
 xt=x;
 xt(i)=xt(i)+ep(i);
 dragt=xt(1)/((xt(4)-xt(5))*xt(6)*xt(3)^2/xt(2));
 deriv(i)=(dragt-drag)/ep(i);
echo off
end
echo on
pause
%Finally we use these derivatives to get the variance in the drag.
%We employ the matrix form given in class.  Note that this would
%work equally well if we had some covariance between the variables,
%only in that case the matrix of covariance would be more complex.
%Thus:

var=deriv*diag(sigma.^2)*deriv'

%The 95% confidence interval is just twice the square root of
%this variance:
[drag-2*sqrt(var),drag+2*sqrt(var)]
%so the interval doesn't quite contain the theoretical result at
%low Reynolds numbers, although it is really close.  Since, by the
%minimum dissipation theorem, you -must- have a velocity lower
%than that calculated for Stokes Law, either the errors were
%underestimated or there was something (systematic error) wrong
%with the experiment.
pause
%We can answer the second question, which measurement to put 
%energy into improving the accuracy of, by looking at the 
%contribution of each element to the overall variance.  We 
%simply break it out term-by-term:

vararray=deriv.^2.0.*sigma.^2.0

pause
%The overall variance is just the sum of these values.  As you
%can see, the error is totally dominated by the third element
%of the array, which corresponds to the radius of the sphere.
%This is usually the case for this sort of problem, and suggests
%That we can best improve our overall error by making the
%measurement of this quantity more accurate.
pause
%As an alternative, we can do the same problem analytically
%using the propagation of errors formula for addition and
%subtraction, multiplication and division.  This gets a little
%messy, but basically you add the variances for the case
%where you subtract the densities, and you add the fractional
%variances where you are multiplying and dividing.  You also
%have 4x the fractional variance when you square the radius!
%This last occurs because a and a have perfect covariance
%when calculating a*a, and this must be accounted for.
%Anyway, the formula is:
pause

varanal=sigma(1)^2/x(1)^2 + sigma(2)^2/x(2)^2 + 4*sigma(3)^2/x(3)^2;
varanal=drag^2*(varanal + (sigma(4)^2 + sigma(5)^2)/(x(4)-x(5))^2)

%which matches the variance calculated before:
var
pause
%Either formula is correct, but if there is covariance between
%any of the variables, then it is -much- better to use the matrix
%formula!  Writing out all of the covariance terms by hand would
%be a real pain!
echo off