I want to make a program which can fill a 2D space with "pursuit polygons".
The following picture will help you understand better what I mean.
enter image description here You can also look up "pursuit curves" or mice problem or watch this gif

What I've tried so far:

First I tried to produce these polygons by rotation of polygons
the square for example

data = {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}};
Graphics[{Table[{Scale[Rotate[Line[data], 3*x Degree], {x, x}]}, {x, 
0, 20}]}]

enter image description here

then I decided to use spirals
I realised that instead of using n spirals for every n-polygon,
I can use only one(!) and produce the effect I want with "fine tuning"
here are the "pursuit polygons" that I made


a = 41.9;
b = 100;
W = Table[{t^2*Cos[a*t], t^2*Sin[a*t]}, {t, 0, b}];
AppendTo[W, W[[-4]]];

enter image description here enter image description here


S = Table[{t^2*Cos[42.4*t], t^2*Sin[42.4*t]}, {t, 0, 160}];
AppendTo[S, S[[-4]]];

enter image description here


P = Table[{t^2*Cos[45.251*t], t^2*Sin[45.251*t]}, {t, 0, 220}];
AppendTo[P, P[[-5]]];

enter image description here

Then I tried to put all these together by rotating and scaling...
(but I'm not satisfied with the result)
notice that for every polygon I use also the "anti-clockwise" version of the polygon, which produces interesting results

T = Table[{t^2*Cos[41.9*t], t^2*Sin[41.9*t]}, {t, 0, 100}];
AppendTo[T, T[[-4]]]; (*Triangle*)
S = Table[{t^2*Cos[42.4*t], t^2*Sin[42.4*t]}, {t, 0, 160}];
AppendTo[S, S[[-4]]];   (*square*)
S1 = Table[{t^2*Cos[-42.4*t], t^2*Sin[-42.4*t]}, {t, 0, 160}];
AppendTo[S1, S1[[-4]]];    (*Anti-clockwise Square*)
P = Table[{t^2*Cos[45.251*t], t^2*Sin[45.251*t]}, {t, 0, 220}];
AppendTo[P, P[[-5]]];      (*Pentagon*)
P1 = Table[{t^2*Cos[-45.251*t], t^2*Sin[-45.251*t]}, {t, 0, 220}];
AppendTo[P1, P1[[-5]]];     (*Anti-clockwise Pentagon*)
Rotate[Scale[Line[(T)], {2.09, 2.09}], -67 Degree], {29500, 3800}],
Rotate[Scale[Line[S + 0], {1, 1}], -30 Degree], 
Rotate[Scale[Line[(S1)], {.98, .98}], -87.5 Degree], {41000, 
Rotate[Scale[Line[P1], {.605, .665}], -20.5 Degree], {76500, 
6500}], Scale[
  Line[(T)], {1.2, 1.7}], -108 Degree], {50900, -12500}], {1.91, 
Rotate[Scale[Line[(P)], {.525, .665}], -64 Degree], {3000, 36000}],
Rotate[Scale[Line[(P)], {.61, .61}], -87 Degree], {59000, 58000}], 
Rotate[Scale[Line[(T)], {2.1, 1.8}], 12 Degree], {82500, 39500}], 
Rotate[Scale[Line[(T)], {2.3, 1.9}], -133 Degree], {32000, 

enter image description here

Can you find a way to divide and fill any given space with pursuit polygons?
The result would look better if this could work with ANY convex polygon and not only the regular polygons that I used...

  • $\begingroup$ Are there any other names for that? Because I think this was asked before. $\endgroup$
    – Kuba
    May 2, 2017 at 10:05
  • $\begingroup$ I provided all the names I could think of. There are not many images outthere.Even the first picture is something that I made a while ago. $\endgroup$
    – ZaMoC
    May 2, 2017 at 10:14
  • 2
    $\begingroup$ @Kuba Lucas, E. (1877). Problème des trois chiens. Nouvelle Correspondance Mathématique, 3, 175–176. $\endgroup$
    – Michael E2
    May 2, 2017 at 11:54
  • 3
    $\begingroup$ A related question. $\endgroup$ May 11, 2017 at 13:15
  • $\begingroup$ And a very related blog-post of mine. And if I'm not completely off here, the word "spacefilling" is incorrect in this context. $\endgroup$
    – halirutan
    Jun 11, 2017 at 15:01

2 Answers 2


The reason that these are called "pursuit polygons" is because they are formed from a dynamical system in which different agents pursue each other. Example:


In this image, one agent starts in each corner of the triangle. The agent starting in the lower right corner pursues the agent starting in the top corner, the agent in the top corner pursues the agent in the lower left corner, and the agent in the lower left corner pursues the agent in the lower right corner.

Drawing lines in between the agents yields a series of triangles which shrink as the agents get closer to each other, and also rotate:

Mathematica graphics

The corresponding differential equations are: $$\begin{align*} \dot{\mathbf{x}}_i &= \mathbf{x}_{i+1} - \mathbf{x}_i,\ i\in\{1,\dots,N-1\}\\ \dot{\mathbf{x}}_i &= \mathbf{x}_1 - \mathbf{x}_N,\ i=N, \end{align*}$$ where $N$ is the number of agents. The agents start in corners for the visualizations we are making, but in general they don't have to.

This code solves the differential equations using Euler integration and plots the lines between the agents:

integrate[pts_, h_] := NestList[Nest[step[#, h] &, #, 10] &, pts, 100]
step[pts_, h_] := pts - h (pts - RotateRight[pts])

wrapAround[pts_] := Append[pts, First[pts]]

PursuitPolygon[pts_, h_: 0.01] := Graphics[{
   Line[wrapAround /@ integrate[pts, h]]

For regular polygons, this results in logarithmic spirals, as you have pointed out:

Grid@Partition[Table[PursuitPolygon[N@CirclePoints[n]], {n, 2, 10}], 2]

Mathematica graphics

But this approach is more general and works for any convex polygons, as you required. The following is an example of how we can fill in a triangulated region with spirals using this system:

reg = DiscretizeRegion[Rectangle[], MaxCellMeasure -> 0.1]

Mathematica graphics

Show[PursuitPolygon @@@ MeshPrimitives[reg, 2]]

Mathematica graphics

Wikipedia's entry on the mice problem suggests that the mice move at unit speed, i.e. like this: $$\begin{align*} \dot{\mathbf{x}}_i &= \frac{\mathbf{x}_{i+1} - \mathbf{x}_i}{\|\mathbf{x}_{i+1} - \mathbf{x}_i\|},\ i\in\{1,\ldots,N-1\}\\ \dot{\mathbf{x}}_i &= \frac{\mathbf{x}_1 - \mathbf{x}_N}{\|\mathbf{x}_1 - \mathbf{x}_N\|},\ i=N, \end{align*}$$ This turns out to be difficult to implement because the denominator grows large as the points approach each other. I spoke to halirutan and MichaelE2 about this. Here is MichaelE2's solution, which I implemented for triangles but it could be implemented for any polygon.

PursuitPolygon[{pt1_, pt2_, pt3_}] := Module[{},

  {sol1x, sol1y, sol2x, sol2y, sol3x, sol3y} = NDSolveValue[{
     {x1'[t], y1'[t]} == Normalize[{x2[t], y2[t]} - {x1[t], y1[t]}],
     {x2'[t], y2'[t]} == Normalize[{x3[t], y3[t]} - {x2[t], y2[t]}],
     {x3'[t], y3'[t]} == Normalize[{x1[t], y1[t]} - {x3[t], y3[t]}],
     {x1[0], y1[0]} == pt1,
     {x2[0], y2[0]} == pt2,
     {x3[0], y3[0]} == pt3,
     WhenEvent[Norm[{x2[t], y2[t]} - {x1[t], y1[t]}] < 1*^-8, "StopIntegration"]},
    {x1, y1, x2, y2, x3, y3},
    {t, 0, 10}

  {tmin, tmax} = MinMax@sol1x["Grid"];

  Graphics@Table[Line[{{sol1x[t], sol1y[t]}, {sol2x[t], sol2y[t]}, {sol3x[t], sol3y[t]}, {sol1x[t], sol1y[t]}}], {t, tmin, tmax, (tmax - tmin)/15}]

Show[PursuitPolygon @@@ MeshPrimitives[reg, 2]]

Mathematica graphics

The solution is qualitatively different when the agents move at unit speed because they no longer meet at the centroid. Note that the difference in the picture is also that the density of the lines is different, this density can be adjusted in both solutions but changing some of the fixed numbers I put in. 15 in (tmax - tmin)/15 controls the number of lines in the last solution.


Here's an extension of @C.E.'s code to non-triangular polygons - polygons taken from a Voronoi mesh of roughly equidistributed (https://mathematica.stackexchange.com/a/141215/3056) 31 points over a square. Code is awful, I didn't really have time to think.

PursuitPolygon[pts_] := Module[{vars, sols},
  vars = Table[Unique[], Length@pts, 2]; 
  sols = NDSolveValue[{Sequence @@ (((#'[t] & /@ #1) == 
           Normalize[Subtract @@ Map[#[t] &, {#2, #1}, {2}]]) & @@@ 
        Partition[Join[vars, {First@vars}], 2, 1]), 
     Sequence @@ MapThread[(#[0] & /@ #1) == #2 &, {vars, pts}], 
     WhenEvent[Norm[x] < 1*^-8, "StopIntegration"] /. 
      x -> Subtract @@ Map[#[t] &, Take[vars, 2], {2}]}, 
    Flatten@vars, {t, 0, 10}];
  {tmin, tmax} = MinMax@First[sols]["Grid"];
     Map[#[t] &, Join[Partition[sols, 2], {Take[sols, 2]}], {2}]], {t,
      tmin, tmax, (tmax - tmin)/15}]]

reg = 
 With[{reg = Rectangle[]}, 
  With[{points = 31, samples = 4000, iterations = 100}, 
    Nest[With[{randoms = Join[#, RandomPoint[reg, samples]]}, 
        Mean@randoms[[#]] & /@ 
         Values@PositionIndex@Nearest[#, randoms]]] &, 
     RandomPoint[reg, points], iterations]] // 
   VoronoiMesh[#, {{0, 1}, {0, 1}}] &];

Show[PursuitPolygon /@ MeshPrimitives[reg, 2][[All, 1]]]

enter image description here


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.