We consider the following problem: a jogger is running along his favorite trail on the plane in order to get his daily exercise. Suddenly, he is being attacked by a dog. The dog is running with constant speed towards the jogger. Compute the orbit of the dog.
The orbit of the dog has the property that the velocity vector of the dog points at every time to its goal, the jogger. We assume that the jogger is running on some trail and that his motion is described by the two functions and .
Let us assume that for the dog is at the point , and that at time his position will be . The following equations hold:
function [zs,isterminal,direction] = dog(t,z,flag); % global w % w = speed of the dog X= jogger(t); h= X-z; nh= norm(h); if nargin < 3 | isempty(flag) % normal output zs= (w/nh)*h; else switch(flag) case 'events' % at norm(h)=0 there is a singularity zs= nh-1e-3; % zero crossing at pos_dog=pos_jogger isterminal= 1; % this is a stopping event direction= 0; % don't care if decrease or increase otherwise error(['Unknown flag: ' flag]); end end
The main program main2.m defines the initial conditions and calls ode23 for the integration. We have to provide an upper bound of the time for the integration.
% main2.m global w y0 = [60;70]; % initial conditions, starting point of the dog w = 10; % w speed of the dog options= odeset('RelTol',1e-5,'Events','on'); [t,Y] = ode23('dog',[0,20],y0,options); clf; hold on; axis([-10,100,-10,70]); plot(Y(:,1),Y(:,2)); J=[]; for h= 1: length(t), w = jogger(t(h)); J = [J; w']; end; plot(J(:,1), J(:,2),':');The integration will stop either if the upper bound for the time is reached or if the dog catches up with the jogger. For the latter case we set the flag Events of the ODE options to 'on'. This tells the integrator to check for zero crossings of the function dog called with flag = 'events'. After the call to ode23 the variable Y contains a table with the values of the two functions and . We plot the orbit of the dog simply by the statement plot(Y(:,1),Y(:,2)). In order to show also the orbit of the jogger we have to compute it again using the vector and the function jogger.
Let us now compute a few examples. First we let the jogger run along the -axis:
function s = jogger(t); s = [8*t; 0];In the above main program we chose the speed of the dog as , and since here we have the jogger is slower. As we can see the dog is catching the poor jogger.
If we wish to indicate the position of the jogger's troubles, (perhaps to build a small memorial), we can make use of the following file cross.m
function cross(Cx,Cy,v) % draws at position Cx,Cy a cross of height 2.5v % and width 2*v Kx = [Cx Cx Cx Cx-v Cx+v]; Ky = [Cy Cy+2.5*v Cy+1.5*v Cy+1.5*v Cy+1.5*v]; plot(Kx,Ky); plot(Cx,Cy,'o');The cross in the plot was generated by appending the statements
p = max(size(Y)); cross(Y(p,1),Y(p,2),2)to the main program.
The next example shows the situation where the jogger turns around and tries to run back home:
function s = jogger1(t); % if t<6, s = [8*t; 0]; else s = [8*(12-t) ;0]; endHowever, using the same main program as before the dog catches up with the jogger at time .
Let us now consider a faster jogger running on an ellipse
function s = jogger2(t); s = [ 10+20*cos(t) 20 + 15*sin(t)];If the dog also runs fast (), he manages to reach the jogger at time .
We finally consider an old, slow dog (). He tries to catch a jogger running on a elliptic track. However, instead of waiting for the jogger somewhere on the ellipse, he runs (too slow) after his target, and we can see a steady state developing where the dog is running on a closed orbit inside the ellipse.
Peter Arbenz 2008-09-24