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 first problem deal with the right time I should stop computations since when the dog catches up the jogger we have a singularity in the differential 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 differential 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 find 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 find 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 first the distance closed and then opened up for a while and finally 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 roundoff errors causing q*40/40 to slightly exceed q, in which case you would get an error message.