I've made a program that is just some points on a plane. Red points are cops and green points are robbers. Each point has a vector indicating where it is going to be in the next step of $t$. Right now $t = 0$ but I want to see what happens when $t=1,2,3,\ldots$

Here is a picture:

enter image description here

The points are in this list. rob[i] and cop[i] make random pairs of x-y cordinates:

r = Table[rob[i], {i, 1, n}],
c = Table[cop[i], {i, 1, m}]

The next position are given by the formulae:

Table[r[[i]] + Normalize@Sum[(r[[i]] - c[[j]])/Norm[r[[i]] - c[[j]]]^2, {j, 1,Length[c]}],{i,1,Length[r]}]


Table[c[[i]] + Normalize@ Sum[(r[[i]] - c[[j]])/Norm[r[[i]] - c[[j]]]^2,{j, 1,Length[r]}],{i,1,Length[c]}]

for robbers and cops, respectively. I want to do something like add a parameter $t$ and when that parameter changes the next position of the points should be updated by the above formulae.

Thanks for helping.


Using the answer of asterix314 I came accross this:

enter image description here

  • $\begingroup$ Please check the formula for the next position of cops. Should the formula be symmetric for robbers and cops? $\endgroup$
    – asterix314
    Commented Jan 15, 2014 at 2:15
  • $\begingroup$ All formulas are actually valid, also what do you mean by symetric? $\endgroup$
    – Arcotick
    Commented Jan 15, 2014 at 7:40
  • $\begingroup$ Wow, beautiful. So the cops are really chasing after the robbers! But if you look at the definition of r[t] and c[t] in my code, they are symmetric meaning they are defined the same way but with r and c switched. So the robbers should chase after the cops also. Your formula is not strictly symmetric (opposite sign). Also in the second one for c, the iterator inside Sum should be j instead of i? $\endgroup$
    – asterix314
    Commented Jan 16, 2014 at 4:59
  • $\begingroup$ The robber in the upper-right corner first takes a U-turn and then quickly turns right to run away from the 2 cops. How exciting! $\endgroup$
    – asterix314
    Commented Jan 16, 2014 at 5:03
  • $\begingroup$ @asterix314 Yeah you're absolutely right, the iterator inside sum should be j instead of i, so if it's a cop the correct next[] would be: Normalize@Sum[c[t][[i]]-r[t]/...], if it's a robber, Normalize@Sum[r[t][[j]]-c[t]/...]. Just change the r[0] and c[0] position, in all parts of the formula. i=length[c[t],j=length[c[t]]. I guess this explains it all. $\endgroup$
    – Arcotick
    Commented Jan 16, 2014 at 11:29

2 Answers 2


What if we made time continuous? The discrete formulae transform to a system of differential equations:

$$ \mathbf{r_i}'(t) = \text{Normalize@}\sum_j{\frac{\mathbf{r_i}(t)-\mathbf{c_j}(t)}{{|\mathbf{r_i}(t)-\mathbf{c_j}(t)|}^2}} $$

$$ \mathbf{c_i}'(t) = -\text{Normalize@}\sum_j{\frac{\mathbf{c_i}(t)-\mathbf{r_j}(t)}{{|\mathbf{c_i}(t)-\mathbf{r_j}(t)|}^2}} $$

Some observations:

  • The direction of a robber is determined by averaging the direction of all robber-cop pairs, weighted by the inverse of distance (dominated by the nearest cop).
  • Robbers run away from the cops, and cops run towards the robbers, thus the opposite signs of velocity.
  • The magnitude of velocity is always 1 for both robbers and cops (Normalized). This means robbers will never be caught.

NDSolve can solve for complex numbers which we'll use to represent the planar points. To keep the robbers from running too far away, 4 additional guarding cops (who remain invisible) are placed on the sides of a 2x2 square.

solveChase[cops_Integer, robbers_Integer, duration_Integer, opts: OptionsPattern[]] :=
    Module[{b, v, c, r, s, t},
        (* guards on the 4 sides to keep robbers fromm running away. *)
        b[z_] := {Re[z] + I, Re[z] - I, 1 + Im[z] I, -1 + Im[z] I};
        (* direction of velocity *)
        v[z_, g_List] := Normalize@Sum[(z - g[[i]])/Abs[z - g[[i]]]^2, {i, Length[g]}];
        s = NDSolve[Join[
            Table[r[i]'[t] == v[r[i][t], b[r[i][t]]~Join~Array[c[#][t] &, cops]], {i, robbers}],
            Table[c[i]'[t] == -v[c[i][t], Array[r[#][t] &, robbers]], {i, cops}],
            (* initial positions *)    
            Thread[Array[r[#][0] &, robbers] == RandomComplex[{-1 - I, 1 + I}, robbers]],
            Thread[Array[c[#][0] &, cops] == RandomComplex[{-1 - I, 1 + I}, cops]]], 
            Array[c, cops]~Join~Array[r, robbers],
            {t, 0, duration}, 
            FilterRules[{opts}, Options[NDSolve]]];
        {Array[c, cops], Array[r, robbers]} /. s // First]

Let's try

cops = 8;
robbers = 10;
duration = 5;
Clear[r, c]
Evaluate@{Array[c, cops], Array[r, robbers]} =
    solveChase[cops, robbers, duration, PrecisionGoal -> 3];
    {Re[#], Im[#]} & /@ Array[r[#][t] &, robbers],
    {Re[#], Im[#]} & /@ Array[c[#][t] &, cops]},
    Frame -> True, AspectRatio -> 1, Axes -> False, 
    PlotStyle -> {{Blue, PointSize[Medium]}, {Red, PointSize[Medium]}}, 
    PlotRange -> {{-1, 1}, {-1, 1}}],
    {t, 0, duration}]

enter image description here

Numerical Instability

Because the equations are inherently non-linear, the numerical solutions are unstable. NDSolve may sometimes complain of "Maxium number of 10000 steps reached ...", which can be avoided by lowering the PrecisionGoal or raising MaxSteps.

Here we solve from the same initial positions, but use different PrecisionGoals:

Clear[cLow, cHigh, rLow, rHigh]
Evaluate@{Array[cLow, cops], Array[rLow, robbers]} = 
    solveChase[cops, robbers, duration, PrecisionGoal -> 3];
Evaluate@{Array[cHigh, cops], Array[rHigh, robbers]} = 
    solveChase[cops, robbers, duration, PrecisionGoal -> 4, MaxSteps -> 100000];

And we pick the trajectories of robber and cop number 4. In both cases, the trajectories of low and high precisions initially coincide, but inevitably at some point they will split and go their separate ways. As a result these trajectories are only accurate for a very small duration of time.

    {Re[rLow[4][t]], Im[rLow[4][t]]},
    {Re[rHigh[4][t]], Im[rHigh[4][t]]},
    {Re[cLow[4][t]], Im[cLow[4][t]]},
    {Re[cHigh[4][t]], Im[cHigh[4][t]]}},
    {t, 0, duration}, 
    PlotStyle -> {
        {Thickness[0.02], LightGray}, {Dashed, Black},
        {Thickness[0.02], LightBlue}, {Dashed, Red}}, 
    AspectRatio -> 1,
    PlotLegends -> {"robber low", "robber high", "cop low", "cop high"}]

enter image description here

  • $\begingroup$ Very good and complete answer. $\endgroup$
    – Arcotick
    Commented Jan 19, 2014 at 22:59
  • $\begingroup$ Actually, is there a better method to put walls there? Adding 1 cop each side is not natural, and they are running from walls. $\endgroup$
    – Arcotick
    Commented Jan 31, 2014 at 12:31
  • $\begingroup$ Putting guards on walls is very easy to implement :) For a more "natural" wall, perhaps you can play with the WhenEvent conditions of NDSolve? $\endgroup$
    – asterix314
    Commented Feb 5, 2014 at 5:03

Let's make r[t] and c[t] depend on time step t. Suppose we have 5 robbers and 3 cops, and their initial positions are:

r[0] = RandomReal[10, {5, 2}]
c[0] = RandomReal[10, {3, 2}]

Given the current position x of a robber and a list y of that of the cops, define the robber's next step (the arrow) to be:

step[x_, y_List] := Normalize@Sum[(x - y[[i]])/Norm[x - y[[i]]]^2, {i, Length[y]}]

The next position of the robbers and the cops are:

r[t_] := r[t] = r[t - 1] + (step[#, c[t - 1]] & /@ r[t - 1])
c[t_] := c[t] = c[t - 1] - (step[#, r[t - 1]] & /@ c[t - 1])

Note the opposite signs of the step, as robbers are running away from the cops and the cops are running towards the robbers. Memoization is used here to avoid recalculation. If you need to change the initial values then you have to Clear[r, c]. Now we can ask e.g.:

(* {7.63022, 4.67586} *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.