function canonballbig
% Canon ball flight. Zero drag, constant gravity ballistics.
% Like the homework due Sept 14, 2006.
% OK, this one is easy with pencil and paper.
% Its just computer practice.
% Fancy stuff here: we end the integration when something interesting happens,
% namely when the projectile hits the ground.
% To do this we need to
% 1) use odeset to set options to include events
% and also name the file where the 'event' is defined.
% In this case the 'event' is going to be y passing through zero
% from above. Canon ball hits ground. Stop calculating.
% 2) use a fancier calling of ode45 that returns not
% just the value of t and z but also what happenned at the 'event'.
% tevent = the time the event occured
% zevent = the value of z (the state) when the event occured
% ievent = which event occurred. In this example we only have
% one event, so this is dumb.
% 3) Define an 'event' function. Here we call it fstopevent.
% has to return 3 things.
% value = the thing which, when it passes through
% zero, should make the integration stop.
% In this case its the height y.
% isterminal = set to 1 if you want the integration to stop,
% and we do in this example.
% direction = set to -1 if you only want to end the integration
% if y is decreasing when it passes through zero.
% Setting to zero will stop the integration at
% any time y=0, which could be right at the start.
% In another file is the same program, functionally, just without
% all the comments and stuff done to make it all more useable and
% readable.
%Set parameters used in the ODE solution
% Another fancy thing: we are going to pass the parameters
% to the ODE instead of defining them inside the ODE function.
% This adds parameters to the call of ode45 and to the definition
% of the right hand side. These are the parameters
g = 10; c=2;
%Set initial conditions
v0= 100; % initial speed
theta = 45*pi/180; % launch angle
vx0 = v0*cos(theta);
vy0 = v0*sin(theta);
x0=0;
y0=0;
zzero=[x0 vx0 y0 vy0]; % pack initial conditions into z0
% Note the huge integration time below. Why?
% Because we are using 'events' to stop the integration when
% the canan ball hits the ground. Thus we need to leave enough
% time to make sure that happens. By using only 2 numbers for
% tspan we are letting the computer decide what times to give
% us output.
tspan=[0 1000]; %specify start and stop integration time
% Setup options for solving the ODE
% The new thing here is 'events'. In the command that follows
% we are saying that the 'even' is in a function called fstopevent
% that is further down in this same text file.
options=odeset('abstol',1e-12,'reltol',1e-12,'events',@fstopevent);
% This next line is the core of the whole file.
% Note that the parameters g and c are passed to
% the right-hand-side function f.
% SOLVE THE ODE
[t,z,tevent, zevent, ievent] = ode45(@f,tspan,zzero,options,g,c);
% In this case, with a terminating 'event', z(end) = zevent.
%compare solution with the analytic solution
format long % show lots of digits in printout
analyticsolution = 2*v0^2*sin(theta)*cos(theta)/g;% anaylic range
computersolution = zevent(1); % also = z(end,1)
numericalerror= computersolution - analyticsolution % This gets printed, no semicolon
%Unpack the solution so you can read it easily
x=z(:,1);
y=z(:,3);
plot(x,y)
title('Andy''s Canon ball flight 9/14/06: no air friction')
xlabel('x'); ylabel('y')
axis('equal') % make x and y axis have same scale so shape looks right
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function defines the right hand side of the ODE
function zdot = f(t,z,g,c);
% Note, after the t and z is the same list of parameters
% used in the calling ode 45: g and c.
% We don't use c in this example, but you might want
% to for your homework.
%Unpack z into variables one can understand
x = z(1);
vx = z(2);
y = z(3);
vy = z(4);
% calculate derivatives, this is the differenctial equation
xdot = vx;
vxdot = 0;
ydot = vy;
vydot = -g;
% pack up the derivatives into zdot.
zdot = [xdot vxdot ydot vydot]';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function defines the 'events'
% In this case the 'event' is to stop integrating
% when the projectile hits the ground.
function [value,isterminal,direction] = fstopevent(t,z,g,c)
value = z(3) ; % look for zero crossings of y=z(3)
isterminal = 1; % 1 = stop, 0 = evaluate and move on
direction = -1; % -1 means decreasing, 0 = who cares, 1 = increasing
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%