How to make a point's position time-dependent given a formula for the next step?

Question

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:

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]}]

and

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.

EDIT:

Using the answer of asterix314 I came accross this:

Answer 1

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];
Manipulate[ListPlot[{
    {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}]

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]
SeedRandom[0]
Evaluate@{Array[cLow, cops], Array[rLow, robbers]} = 
    solveChase[cops, robbers, duration, PrecisionGoal -> 3];
SeedRandom[0]
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.

ParametricPlot[{
    {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"}]

Answer 2

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.:

c[8][[2]]
(* {7.63022, 4.67586} *)