clear
echo on
%Solution to Problem 3
%
%Save the file data.dat to the appropriate directory
%before running this program!
%
%In this problem we must determine the
%initial velocity, exit time, and deceleration
%due to air resistance of a cannonball exiting from
%a cannon.  We want to use linear regression,
%but we must rework the model so that it is
%linear in the modelling parameters.  The original
%model is given by:
%
%  b =  U*(t - t0) - 0.5 * a * (t-t0)^2
%
%which is non-linear in t0, the exit time.
pause
%We can linearize this by multiplying the formula
%out and collecting the terms:
%
%  b = [-U*t0 - 0.5*a*t0^2] + t*[U + a*t0] + t^2*[-0.5*a]
%
%If we let the terms in brackets be the new modelling
%parameters x(1), x(2), and x(3), then the model is
%clearly linear in these parameters.  We may now begin:
pause
%First we read in the data file.  We have stored it
%under the file name 'data.dat'.
load 'data.dat'
t=data(:,1)
b=data(:,2)
pause
%Next we construct the matrix of modelling functions:
n=length(t);
a=[ones(n,1),t,t.^2]
pause
%And we get the solution vector x:
x=a\b
pause
%Now we have to determine the relationship between
%the new modelling parameters and the quantities
%we really want.  Looking at the above equations we
%get:
acc=-2*x(3)
pause
%The next two are harder - a pair of non-linear
%equations.  Rearranging things a bit yields:
%
%  x(3)*t0^2 + x(2)*t0 + x(1) = 0
%
%Thus we can solve for t0:
t0=(-x(2)+(x(2)^2-4*x(1)*x(3))^.5)/2/x(3)
%Where we have chosen the positive root as the
%negtive root lacks physical significance.
pause
%The initial velocity is thus given by:
U=x(2)-acc*t0
%which finishes the problem.
pause
%We can also plot this up:
figure(1)
plot(t,b,'o',t,a*x)
xlabel('time (sec)')
ylabel('height (m)')
echo off