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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%