clear
echo on
%Solution to Problem 29
%
%In this problem we are to calculate both
%the maximum height reached by a cannon ball
%fired straight up, and the time necessary 
%for it to hit the ground again.  We will employ
%a simple explicit Euler method routine first,
%and then we will use a higher order method for
%comparison.
pause
%To run this problem, we first define a vector
%containing the position and velocity.  We let
%y(1) be the height of the ball and y(2) be the
%velocity.  We also write the attached program
%cannon.m which returns the derivatives of these
%quantities as a function of y and t.  Actually,
%the derivatives -aren't- an explicit function of
%time, but packaged ode solvers require time in
%the function call for generality.
pause
%First we set the initial conditions:
y0=[0,2e4];
%Then we do the integral:
dt=0.01;
t=0;
y=y0;
while y(2)>0,
 ylast=y;
 y=y+cannon(t,y)*dt;
 t=t+dt;
end
pause
%The maximum height is achieved when the upward velocity
%reaches zero.  This time will be slightly less than the
%time after the last iteration, and the maximum height
%will be slightly greater.  We can get a better estimate
%of the time and height by interpolation:
delta=-dt*y(2)/(y(2)-ylast(2));
tmax=t+delta
hmax=y(1)-delta*y(2)+0.5*delta^2*(y(2)-ylast(2))/dt
%where the last little bit is an estimate of the second
%derivative of the height with time.
pause
%We were also asked to calculate the time for the cannon
%ball to hit the ground again.  This is also easy to calculate:
t=0;
y=y0;
while y(1)>=0,
 ylast=y;
 y=y+cannon(t,y)*dt;
 t=t+dt;
end
pause
%In this case we terminate the integration when the ball's
%position becomes negative.  The precise time can be obtained
%with interpolation again:
thit=t-dt*y(1)/(y(1)-ylast(1))
%which is much more than twice the time to get to the maximum
%height.  The difference is due to the air drag.
pause
%We can improve our estimates by using a higher order integration
%scheme.  Let's use the 2-stage Runge-Kutta method.  Thus:
t=0;
y=y0;
while y(2)>0,
 ylast=y;
 k1=cannon(t,y)*dt;
 k2=cannon(t+dt,y+k1)*dt;
 y=y+(k1+k2)/2;
 t=t+dt;
end
delta=-dt*y(2)/(y(2)-ylast(2));
tmax=t+delta
hmax=y(1)-delta*y(2)+0.5*delta^2*(y(2)-ylast(2))/dt
pause
%And for the total time:
t=0;
y=y0;
while y(1)>=0,
 ylast=y;
 k1=cannon(t,y)*dt;
 k2=cannon(t+dt,y+k1)*dt;
 y=y+(k1+k2)/2;
 t=t+dt;
end
thit=t-dt*y(1)/(y(1)-ylast(1))
%which yields values close to those obtained before.
pause
%Note that if we had used even higher order integration
%schemes, the accuracy of our final result would have
%been no better, because it would have been limited by
%the interpolation scheme at the end rather than the
%integration procedure itself!  If we hadn't interpolated
%the final time, for example, the error in this quantity
%would have been O(dt), the same as the order of the
%Euler method.  This is an important point to remember -
%the accuracy of any numerical method is ultimately determined
%by the LEAST accurate step in the process!
echo off



________________________


function f=cannon(t,y)
%This function returns the derivatives of the
%height and velocity of the cannon ball.
%We define the parameter values:
rho=0.001;
r=10;
m=3000;
g=980;
%Thus we have the derivatives:
f(1)=y(2);
f(2)=-g - sign(y(2))*0.5*pi*r^2*(y(2))^2*rho/m;