246 Math - ODE for Scientists and Engineers - Section 0101 Instructor: Prof. Dolzmann Fall Term 2003 How to use the event locator option It seems that the syntax for the event locator option is slightly dierent from the one outlined in CHLOS in recent releases of MATLAB. Here is one way to implement the example on pp 108-110. There are three steps involved. a) Denition of the event locator function as an m-le. We dene a function events(t,y) (you can choose any name) which has to return three output arguments, the value of the function at (t, y), information whether the computation is to be continued after the event (in isterminal), and a ag that determines whether the slope of the event function at the time of the event is relevant (in direction). Here is the le for the event y = 1.5 which we write in standard from y 1.5 = 0. % % The definition of the event y-1.5=0. % % isterminal=1 indicates that the evaluation stops % direction=0 indicates that all zeros are computed % (y increasing or decreasing) % function [value,isterminal,direction] = events(t,y) value = y-1.5; isterminal=1; direction=0; b) Denition of f (t, y) as an m-le. If you use the event locator property in ode45, you need to dene the right-hand side f of your ODE y = f (t, y) as an m-le: % % the right hand side of the ODE % function u = f(t,y) u=t./y; c) Setting the parameter and calling the solver. Finally we need to tell the solver that we want to test for an event and how the event function is called. Note that the event function must have the standard form discussed above. Since it is an m-le, we have to use the notation @events in order to get the function handle which is needed as the input argument. The same is true for the function f in the call to ode45. Here is the MATLAB code: % % turn events on and define the events function % options = odeset(events,@events); [t, y, tev, yev, ie] = ode45(@f,[0,2],1,options) plot(t,y)
1
2
Outline for the paratrooper problem You can solve the second part of the paratrooper problem with the events option in ode45. Here is some outline of the code MATLAB you might want to use. % % % % % % % Separate m-file for the event function Definition of the event. Note that y is a vector with two components. Make sure that your condition implies that the height of the paratrooper is equal to zero. function [value,isterminal,direction] = hitground(t,y) value = ...; isterminal=1; direction=0; % % % % % % % % Separate m-file for the right-hand side f(t,y) while the parachute is open. The forces acting on the paratrooper when the parachute is open. Make sure you use the right constant for air resistance. Recall that you need to return a vector. function f = paratrooper(t,y) f = [... ; ...]; % % % % % % % % % % % % Here are the commands to solve the problem free fall, Newtons law F=ma; y(t) F = height at time t [sec] = force due to gravitational acceleration (-mg=-32m) and drag force due to air resistance +\gamma (y)^2
recall that the mass of the paratrooper is given by 195/32 = 6.0938
close all; figure; hold on % % for simplicity we use an inline function in the first part. You can % also write a separate m-file if you wish to do so % force = inline([... ; ... ], ... ); % % call ode45 with the correct parameters and initial conditions %
3
[xp,yp]=ode45( ... ); % % Check: your s...