#### 6.43 animation problem (15.5.00)

##### 6.43.1 Aldrovando Azeredo Araujo

I’ve listed the code of a simple program I’ve been dealing with. I am trying to show my students simple models like the one below that treats chasing situations. Here we have a dog that chases a jogger. The ﬁrst problem deal with the right time I should stop computations since when the dog catches up the jogger we have a singularity in the diﬀerential equation.

So I would like to know the time the dog reach the jogger to plot the two trajectories using insequence=true. I want to put all this comands in a procedure. How do I do that?

As anyone can see I had to make a sequence of plots that take a lot of time to be generated. Is it possible to get this animation faster? I can do that in matlab quite fast.


> restart:
> with(plots):
> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> x0:=2: y0:=1:
>
>
> K:=sqrt((X(t)-x(t))^2+(Y(t)-y(t))^2);

2                2
K := sqrt((X(t) - x(t))  + (Y(t) - y(t)) )

> eq1:=diff(x(t),t)=(v/K)*(X(t)-x(t));

d                    v (X(t) - x(t))
eq1 := -- x(t) = -------------------------------------
dt                          2                2
sqrt((X(t) - x(t))  + (Y(t) - y(t)) )

> eq2:=diff(y(t),t)=(v/K)*(Y(t)-y(t));

d                    v (Y(t) - y(t))
eq2 := -- y(t) = -------------------------------------
dt                          2                2
sqrt((X(t) - x(t))  + (Y(t) - y(t)) )

> X:=t->(t+1)*sin(2.5*t);

X := t -> (t + 1) sin(2.5 t)

> Y:=t->(t+1)*cos(2.5*t);

Y := t -> (t + 1) cos(2.5 t)

> v:=4.6*t:
> var:={x(t),y(t)}:
> init:={x(0)=x0,y(0)=y0}:
> sys:={eq1,eq2}:
> Sol:=dsolve(sys union init, var, numeric, output=listprocedure);

Sol := [t = (proc(t)  ...  end), x(t) = (proc(t)  ...  end),

y(t) = (proc(t)  ...  end)]

> xt:=subs(Sol,x(t));

xt := proc(t)  ...  end

> yt:=subs(Sol,y(t)):
> Tdog[0]:=[x0,y0]:Tman[0]:=[X(0),Y(0)]:
> xt(1.98);solve(r*0.05=1.98);

-2.897729591544956

39.60000000

> for i from 1 by 1 to 39 do
Tdog[i]:=Tdog[i-1],[xt(i*0.05),yt(i*0.05)]:
Tman_dog[i]:=plots[display]({plot([Tdog[i]],color=blue),
plot([X(t),Y(t),t=0..i*0.05],color=red)}):
od:
> display([seq(Tman_dog[i],i=1..39)],insequence=true);



##### 6.43.2 Willard, Daniel Dr (17.5.00)

There is a book on diﬀerential gaming that you might wish to read: "Advances in Missile Guidance" purchased from AIAA Publications, 9 Jay Gould Court, Waldorf, MD 20602. It treats such problems and also ones where the dog tries to escape, and gives some Maple programs to go along with it.

##### 6.43.3 John Reinmann (19.5.00)

In answer to your question "how can I ﬁnd the time (T_catchup) when the dog catches up with the jogger?", I propose the following approach:

When you obtain the numerical solutions for xt and yt, you ﬁnd that when the time exceeds the T_catchup, Maple gives the error messages "sqrt of a negative number" or "division by zero".

Thus, you can use this error message to signal when you have overestimated the value of T_catchup. You use Maple’s "traperror" command.

Taking both error messages into account and improving the way I change "del" in the "for" loop, I now recommend the following Maple code to determine q = T_catchup (where, q denotes the time):

>q:=-1:
del:=1:
A:=sqrt of a negative number:
B:=division by zero:
test:=false:
for i from 1 to n do
while test=false do
q:=q+del:
result:=traperror(yt(q)):  # xt(q) could be used instead of yt(q)
test:=evalb((result=A) or (result=B)):
od:
q:=q-del:
del:=del/10:
test:=false:
od:
q;



The result depends on the number of Digits used. The following results are for Digits=12.

for n=4, q=1989331/1000000

for n=5, q=19833319891/10000000000



In either case, if you increase the last digit in the numerator by 1, you get the error message above; so these results are very good.

The test I used to see if the value of T_catchup was good, was to evaluate the expression for the distance between the jogger and the dog, which is

distance = sqrt((X(q)-xt(q))^2+(Y(q)-yt(q))^2)

for n=4, distance = 0.16594...e-5

for n=5, distance = 0.63480...e-7

For n=11, the above code gives q=19893319961/100000000000 = 1.9893319961



When I plotted distance versus time, I got the interesting result that at ﬁrst the distance closed and then opened up for a while and ﬁnally the distance closed again and the dog caught up with the jogger.

For an animation that can be generated quickly, I suggest using "pointplot", which can show the position of the dog and the jogger simultaneously as symbols (say a blue circle for the jogger and a red diamond for the dog).

Actually, this is a lot of fun to observe. Here’s the code:

let q be the time at which the dog catches up with the jogger, and xt(t), yt(t), X(t), Y(t) be the same as in your code.

>display(seq(display({pointplot([X(q*j/40.0001),Y(q*j/40.0001)],color=blue,
symbol=circle),pointplot([xt(q*j/40.0001),yt(q*j/40.0001)],color=red,
symbol=diamond)}),j=0..40),insequence=true,title=red-dog blue-jogger,
scaling=constrained);



I used q*j/40.0001 rather than q*j/40 to avoid the possibility of roundoﬀ errors causing q*40/40 to slightly exceed q, in which case you would get an error message.