function tDevour = timeToDeath_2( angle ) % t_devour = timeToDeath( angle ) % angle should be in degrees % Solves part (a) for problem 1, HW3, ACM11 % to find the best (or worse) angle to run at, you can do the following: % [ang,val,flg] = fminbnd(@timeToDeath,-180,180) % (which is the fastest way to die), or % [ang,val,flg] = fminbnd(@(t)-timeToDeath(t),-180,180) % (the best way to run) % % This version uses the events function for odeset % == The human starts at (0,0) and the raptors start at... == a = 10*sqrt(3)/3; R1 = [-10,-a]; R2 = [10,-a]; R3 = [0,2*a]; % == Setup the human's path == theta = deg2rad(angle); c = [cos(theta),sin(theta)]; c = c/norm(c); H_SPEED = 6; h = @(t) H_SPEED*t*c; % person runs is a straight line % == Setup the raptors' equations of motion == R_SPEED = 25; R_SPEED2 = 20; % This is for the injured raptor odefun = @(t,y) R_SPEED * (h(t)-[y(1),y(2)])' ./ norm(h(t)-[y(1),y(2)]); odefun2 = @(t,y) R_SPEED2* (h(t)-[y(1),y(2)])' ./ norm(h(t)-[y(1),y(2)]); tol = .1; % The human is caught when the raptors are this close % When the raptor has caught the human, we'll stop modelling the true % motion, because dividing by norm( x - y ) will result in huge numbers % An alternative is to use the 'Events' option -- see "odeset". % This is what we did in version 1: % cutoff = @(t,y) norm( h(t)-[y(1),y(2)]) > tol ; % odefun = @(t,y) odefun(t,y) .* cutoff(t,y); % odefun2 = @(t,y) odefun2(t,y) .* cutoff(t,y); % == Run the simulation == opts = odeset('abstol',1e-3); % set the accuracy of the solver % We need to tell ode45 the start and end time. The start time can be 0, % wlog. The end time can be anything greater than t_devour -- but we don't % know what d_devour is yet! So, pick something conservative, like 10 % seconds. Once you run this a few times, you can pick a more accurate % estimate, like 1 second, or even .5 seconds. Here, we use 1 second, so % our time vector is [0,1]. % Add our events (see the file caught.m) % Here, we pass the trajector of the human to caught.m ev = @(t,y) caught(t,y,h); opts = odeset(opts,'Events', ev ); [T1,Y1] = ode45(odefun ,[0,1],R1,opts); % raptor #1 [T2,Y2] = ode45(odefun ,[0,1],R2,opts); % raptor #2 [T3,Y3] = ode45(odefun2,[0,1],R3,opts); % raptor #3 (injured) % == Calculate t_devour using the results == n = sqrt(sum( (h(T1)-Y1).^2 ,2)); f = find( n < tol ); t1 = T1(f(1)); n = sqrt(sum( (h(T2)-Y2).^2 ,2)); f = find( n < tol ); t2 = T2(f(1)); n = sqrt(sum( (h(T3)-Y3).^2 ,2)); f = find( n < tol ); t3 = T3(f(1)); tDevour = min( [t1,t2,t3] ); if nargout < 1 fprintf('Time until raptors eat you: %.4f seconds\n',tDevour ); end