# Hightlight all the self-intersections of a Lissajous figure

This graph–also known as a Lissajous figure–contains so many self-intersections.How can I highlight them?

ParametricPlot[{Sin[100 t], Sin[99 t]}, {t, 0, 2 π},
PlotRange -> All]

• Too tired to tackle this now, but a hint: for general $n$ (this is case $n=100$) the number of nodes is $2n^2-4n+1$. They are distributed on a $2n-1$ by $2n-3$ grid, occupying a diagonal pattern (like the white squares on a chessboard). The key is going to be to find the positions of the lines of the grid. – 2012rcampion May 11 '15 at 3:27
• Looks like they have x-coordinates Cos[Pi Range[2 (n - 1) - 1]/(2 (n - 1))] and y-coordinates Cos[Pi Range[2 (n - 1) + 1]/(2 (n - 1) + 2)]. – 2012rcampion May 11 '15 at 3:34
• A quick idea: extract the Line[] objects from a plot of the curve (Cases[] is useful here), split any polylines present into simple lines of the form Line[{pt1, pt2}], and then use a line intersection algorithm on the lines produced. Polishing with FindRoot[] is optional. – J. M. will be back soon May 11 '15 at 5:02

Manipulate[
ParametricPlot[({Sin[n t1], Sin[(n - 1) t1]}), {t1, 0, 2 Pi},
Epilog -> {Red, PointSize[Large],
Table[If[OddQ[i + j],
Point[{Cos[Pi i/(2 (n - 1))], Cos[Pi j/(2 (n))]}]], {i,
2 n - 3}, {j, 2 n - 1}]}], {{n, 5}, 2, 20, 1}]


## General Case I

We can generalize to the Lissajous curve specified by the two non-negative integers $a$ and $b$:

$$x = \sin at \\ y = \sin bt \\ t \in [0,2\pi)$$

Without loss of generality, I will assume $b<a$ and $b\nmid a$.

We can start by making a table of small cases:

Column[Row /@
Table[ParametricPlot[({Sin[a t], Sin[b t]}), {t, 0, 2 Pi},
Epilog -> {}, PlotLabel -> {a, b}, Axes -> False,
ImageSize -> Tiny], {a, 10}, {b,
Select[Range[a - 1], CoprimeQ[#, a] &]}], Alignment -> Center]


When both $a$ and $b$ are odd, we get a degenerate curve that traces itself twice. I'll handle those cases later.

It looks like each self-intersection occurs on a horizontal and vertical line shared with several other solutions. We can make a table mapping $a$ and $b$ to the number of horizontal and vertical grid lines (ignoring $b=1$ as a special case for now):

$$3,2\to 3,5\\ 4,3\to 5,7\\ 5,2\to 3,9\\ 5,4\to 7,9\\ 6,5\to 9,11$$

It's fairly evident that the number of grid lines is merely:

$$2b-1,\,2a-1$$

The spacing of the grid lines looks mathematically like it might be more difficult. However, the spacing looks familiar to me: like the spacing of points in an airfoil .dat file:

Graphics[Point@
Rest[Import[


I remember from AE311 (incompressible flow) that this spacing follows the transformation:

$$x\mapsto \frac{c}{2}\left(1-\cos(\theta)\right)$$

with the points evenly spaced in $\theta$. Could it really be that simple?

Manipulate[
ParametricPlot[({Sin[a t], Sin[b t]}), {t, 0, 2 Pi},
GridLines -> {Cos[Pi Range[2 b - 1]/(2 b)],
Cos[Pi Range[2 a - 1]/(2 a)]}, PlotLabel -> {a, b},
Axes -> False], {{a, 5}, 2, 20, 1}, {{b, 4},
Select[Range[a - 1], CoprimeQ[#, a] &]}]


Heck yeah it is: lucky guess!

Note that only every other grid node contains an intersection; they form a sort of checkerboard pattern. This accounts for the seeming fewer number of grid lines when $b=1$: only every other line is occupied, so there are twice (plus one) as many grid lines as intersections.

We can also take a look at the odd-odd special cases:

We can see that they follow a double-size checkerboard pattern, with adjacent intersections two diagonals apart.

With all this in mind, we can now extend the code from the original example:

Manipulate[
ParametricPlot[
{Sin[a t], Sin[b t]}, {t, 0, 2 Pi},
GridLines -> {Cos[Pi Range[2 b - 1]/(2 b)], Cos[Pi Range[2 a - 1]/(2 a)]},
PlotLabel -> {a, b}, Axes -> False,
Epilog -> {Red, PointSize[Large],
Table[If[
If[OddQ[a] && OddQ[b],
EvenQ[i] && Divisible[i + j + a + b + 2, 4],
OddQ[i + j]],
Point[Cos[Pi/2 {i/b, j/a}]]
],
{i, 2 b - 1}, {j, 2 a - 1}]}
],
{{a, 5}, 2, 20, 1}, {{b, 4}, Select[Range[a - 1], CoprimeQ[#, a] &]}
]


## General Case II

We can follow a similar procedure for phased Lissajous curves. Without loss of generality, we can apply a phase $\phi$ to the $x$-coordinate:

$$x = \sin(at +\phi) \\ y = \sin(bt) \\ t \in [0,2\pi)$$

If we apply phases $\phi_a$ and $\phi_b$ to the $x$ and $y$-coordinates, respectively, this is equivalent to a curve with $\phi=\phi_a-\frac a b \phi_b$ and $t'=t+\frac{\phi_b}b$.

First we'll take a look at what's going on:

Manipulate[
ParametricPlot[{Sin[a t + ϕ], Sin[b t]}, {t, 0, 2 Pi},
PlotLabel -> {a, b}, Axes -> False], {{a, 5}, 2, 20, 1}, {{b, 4},
Select[Range[a - 1], CoprimeQ[#, a] &]}, {ϕ, 0, 2 Pi}]


I like to visualize this as the projection of a pattern on the surface of a vertical cylinder, rotating about its axis:

A little bit of work transforms the original solution to follow the intersections of the cylinder pattern:

Note that we're missing half of the intersections now! The missing intersections are where lines from the 'front' half of the cylinder overlap the 'back' half. We can get those via a similar process, treating the pattern as a projection from the surface of a horizontal cylinder. In the image above, we essentially want to reflect the 'missing' intersections across the diagonal:

This gives us our final result:

Manipulate[
With[{gcd = GCD[a, b]},
With[{a = a/gcd, b = b/gcd},
ParametricPlot[
{Sin[a t + ϕ], Sin[b t]}, {t, 0, 2 Pi},
PlotLabel -> {a, b}, Axes -> False,
Epilog -> {
PointSize[Large],
Red,
Table[
If[EvenQ[i + j],
Point[{Sin[2 Pi (i + a)/(2 b) + ϕ], Cos[Pi j/a]}]
],
{i, 2 b}, {j, a - 1}
],
Orange,
Table[
If[EvenQ[i + j],
Point[{Cos[Pi i/b], Sin[2 Pi (j + b)/(2 a) - b/a ϕ]}]
],
{i, b - 1}, {j, 2 a}
]
}
]
]
],
{{a, 6}, 1, 20, 1}, {{b, 13}, 1, 20, 1}, {{ϕ, Pi/10}, 0, 2 Pi}
]


(Note that for some values of ϕ, you will see repeated intersections or intersections at the edge of the curve. This happens when the curve becomes degenerate and overlaps itself.)

• Neat! and really fast as well! (+1) – MarcoB May 11 '15 at 3:52
• @Marco Much faster if you use 2.. – 2012rcampion May 11 '15 at 3:53
• @2012rcampion You code comes first , but it's a pity that it becomes slow when n is large. – WateSoyan May 12 '15 at 11:02
• @Wate Actually the answer to that one is quite simple: the number of intersections for $a$ and $b$ both odd is $(a-1)(b-1)/2$ and the number of intersections in the other cases is $2ab-(a+b)$. The number of regions is simply the number of intersections plus one (every time you close a curve by intersection, you split a previous region into two, and you start with one region). – 2012rcampion May 12 '15 at 15:25
• N.B. These intersection points are related to what are now termed Padua points. – J. M. will be back soon Jul 19 '16 at 14:08

One way (whew, there are a lot of intersections! -- here's a shorter version):

sol = NSolve[{Sin[10 t], Sin[9 t]} == ({Sin[10 t], Sin[9 t]} /. t -> s) &&
0 <= t < s < 2 Pi, {t, s}];

ParametricPlot[{Sin[10 t], Sin[9 t]}, {t, 0, 2 π},
Epilog -> {Red, PointSize[Large], Point[{Sin[10 t], Sin[9 t]} /. sol]}]


({Sin[100 t], Sin[99 t]} will take a lot longer.)

General solution via Mathematica that is not too slow

Solve returns solutions in the form ConditionalExpression, and the condition can be used to generate the values {m, n} for each point of intersection (via Solve inside Block).

gensols = Cases[
Solve[(-a t + a s == 2 Pi m || a t + a s == Pi + 2 Pi m) &&
(b t - b s == 2 Pi n || b t + b s == Pi + 2 Pi n) &&
a > b > 0 && {a, b, m, n} ∈ Integers &&
0 <= t < s < 2 Pi, {s, t}],
HoldPattern[t -> t0_] :> {t -> t0}, 2];

Block[{a = 100, b = 99},
pts = Flatten[
Hold[{Sin[a t], Sin[b t]} /. Solve[Last[t], {m, n}]] /. gensols // ReleaseHold,
1]
] // Length // AbsoluteTiming

(*  {1.55879, 19601}  *)

ParametricPlot[{Sin[100 t], Sin[99 t]}, {t, 0, 2 Pi},
PlotStyle -> {Black, Thickness[0.0015]}, PlotPoints -> 3000,
Epilog -> {GraphicsComplex[N@pts, {Red, PointSize[0.003], Point[Range@Length@pts]}]}]


• 19,601 intersections if my math is correct. – 2012rcampion May 11 '15 at 3:56
• It's a pity that NSolve takes a very long time. – WateSoyan May 11 '15 at 10:45
• @WateSoyan Well, 20K points are a lot to track down and NSolve might be converting the trig. eqns. to polynomial ones -- ouch. The difference between solving (this answer) and having been solved (rcampion2012's) is to be expected. You must be making a pretty big poster to show all those points! :) – Michael E2 May 11 '15 at 12:31
• @Michael E2 I have posted my answer below,purely based on Solve. – WateSoyan May 11 '15 at 13:43
• @WateSoyan To me it's a pity that I preferred to figure out how to get Mathematica to solve the problem, than to solve it mathematically myself, like rcampion2012. (But then again, this site is about how to use M.) – Michael E2 May 11 '15 at 15:20
GraphicsMeshMeshInit[];

eps = 1/1000000;
pp = ParametricPlot[{Sin[10 t], Sin[9 t]}, {t, eps, 2 π}];
intersections = GraphicsMeshFindIntersections[pp];

Show[pp, Epilog -> {Red, PointSize[Large], Point@intersections}]


GraphicsMeshFindIntersections[ParametricPlot[{Sin[100 t], Sin[99 t]},
{t, eps, 2 π}]] // Length // Timing


{0.078125, 20330}

Row[Show[plt = ParametricPlot[{Sin[# t], Sin[(# - 1) t]}, {t, 0, 2 π}],
ImageSize -> 300, Epilog -> {Red, PointSize[.2/#],
Point@GraphicsMeshFindIntersections[plt]}] & /@ {5, 10, 20, 50}]


• Now I remember that function! +1. – Michael E2 May 11 '15 at 17:24
• @Michael, thank you for the vote. – kglr May 11 '15 at 17:39
• Glad to see my idea works. :) – J. M. will be back soon May 11 '15 at 21:44
• @J.M, sorry I was distracted by Cases[...] in your comment above and failed to click the link you provided. Yes, it does work; and directly on the graphics input without the extra need to extract the lines. – kglr May 11 '15 at 21:58
• @WateSoyan, it came up as one of the search results from ??**Intersections*. – kglr May 12 '15 at 11:11

I find a workaround which can find all that exact self-intersections:

    sol = Solve[(100 (t1 - t2) == 2 k1 \[Or]
100 (t1 + t2) == (2 k1 + 1)) && (99 (t1 - t2) == 2 k2 \[Or]
99 (t1 + t2) == (2 k2 + 1)), {t1, t2}];
Flatten[({t1, t2} /. # /. Solve[(0 <= t1 < t2 < 2) /. #, {k1, k2}, Integers]) & /@ sol, 1]

• You have slightly more than twice as many solutions as rcampion2012. I think you want Flatten[({t1, t2} /. # /. Solve[(0 <= t1 < t2 < 2) /. #, {k1, k2}, Integers]) & /@ sol, 1]` to get each solution exactly once. (I was working along similar lines, but I was trying to get the solution even faster.) – Michael E2 May 11 '15 at 13:52
• @Michael E2 Yeah,I forgot to check it. – WateSoyan May 11 '15 at 14:46