# 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. May 11, 2015 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)]. May 11, 2015 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. May 11, 2015 at 5:02

I will start with the most general Lissajous figure, which has a parametric equation as follows:

$$x(t) = f_a(at+\phi_a) \\ y(t) = f_b(bt+\phi_b) \\ t \in [0,2\pi)$$

Where:

$$f_a, f_b \in \{\sin, \cos\} \\ a, b \in \mathbb{N} \\ \phi_a, \phi_b \in [0, 2\pi)$$

Note that this curve is periodic, so we can "shift" the range of $$t$$ without changing the curve. Consider one such shift, given by $$t \to t'-\frac{\phi_b}{b}$$. This gives the new parameterization:

$$x(t') = f_a\left(a\left(t'-\frac{\phi_b}{b}\right)+\phi_a\right) =f_a\left(at'+\left(\phi_a-\frac{a}{b}\phi_b\right)\right) \\ y(t') = f_b\left(b\left(t'-\frac{\phi_b}{b}\right)+\phi_b\right) = f_b(bt') \\ t \in [0,2\pi)$$

This shows that we can eliminate one of the phase parameters and still describe the entire set of Lissajous figures. At this point we can also notice that $$\sin$$ and $$\cos$$ are related by a phase shift of $$\pi/2$$, so we can absorb the choice of function into the phase parameter, giving us the simplified equation:

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

Finally, take $$d=\gcd(a,b)$$. This means there exist integers $$a'$$ and $$b'$$ such that $$a=a'd$$ and $$b=b'd$$, so that:

$$x(t) = \cos(a'(dt)+\phi) \\ y(t) = \cos(b'(dt)) \\$$

This is the same curve as the one described by $$a'$$ and $$b'$$, except covered $$d$$ times. Therefore we can describe the entire set of Lissajous figures if we restrict ourselves to cases where $$a$$ and $$b$$ are coprime, i.e. $$\gcd(a,b)=1$$.

Now consider a self-crossing of the Lissajous figure. Such a point occurs when two distinct values of $$t$$ give the same $$x,y$$ coordinates:

$$x(t_1)=x(t_2) \to \cos(at_1+\phi) = \cos(at_2+\phi) \\ y(t_1)=y(t_2) \to \cos(bt_1) = \cos(bt_2)$$

The cosine function is periodic, so one way for two cosines to be equal is for their arguments to differ by a factor of $$2\pi$$:

$$\alpha=\beta+2\pi n\to \cos\alpha=\cos\beta$$

However, the cosine function is also even, so it is also possible for one argument to be the negative of the other (plus a factor of $$2\pi$$):

$$\alpha=-\beta+2\pi n\to \cos\alpha=\cos\beta$$

Note that at all of the self-crossings, one branch of the curve is moving in an "upwards" diagonal direction and the other is moving in a "downwards" diagonal direction. This means that one of $$x$$ or $$y$$ must obey the "negated" equality, so we have one of:

$$at_1+\phi = at_2+\phi+2\pi n \\ bt_1 = -(bt_2)+2\pi m \\ \textrm{or} \\ at_1+\phi = -(at_2+\phi) + 2\pi n\\ bt_1 = bt_2 + 2\pi m$$

These two systems of equations are linear in $$t_1$$ and $$t_2$$ so they each have a unique solution:

$$t_1 = \left(\frac{n}{a}+\frac{m}{b}\right)\pi - \frac{\phi}{a} \\ t_2 = \left(\frac{n}{a}-\frac{m}{b}\right)\pi - \frac{\phi}{a} \\ \textrm{or} \\ t_1 = \left(\frac{n}{a}+\frac{m}{b}\right)\pi \\ t_2 = \left(-\frac{n}{a}+\frac{m}{b}\right)\pi$$

This gives the following sets of intersection points:

$$(x,y) = \left(\cos\left(n\pi + \frac{a}{b}m\pi\right), \cos\left(\frac{b}{a}n\pi + m\pi - \frac{b}{a}\phi)\right)\right) \\ 0\le n\le a-1,\quad1 \le m \le b-1\\ \textrm{and} \\ (x,y) = \left(\cos\left(n\pi + \frac{a}{b}m\pi + \phi\right), \cos\left(\frac{b}{a}n\pi + m\pi\right)\right) \\ 1\le n\le a-1,\quad0 \le m \le b-1$$

Plotting these points in Mathematica:

Manipulate[
ParametricPlot[{Cos[a t + ϕ], Cos[b t]}, {t, 0, 2 Pi},
PlotLabel -> {a, b}, AspectRatio -> Automatic, Axes -> False,
Epilog -> {Red, PointSize[Large],
Table[Point[{Cos[(n + a/b m)Pi], Cos[(b/a n + m)Pi - b/a ϕ]}],
{n, a}, {m, b - 1}],
Table[Point[{Cos[(n + a/b m)Pi + ϕ], Cos[(b/a n + m)Pi]}],
{n, a - 1}, {m, b}]
}
],
{{a, 5}, 2, 20, 1},
{{b, 4}, Select[Range[a], CoprimeQ[a, #] &]},
{{ϕ, Pi/10}, 0, 2 Pi}
] • Neat! and really fast as well! (+1) May 11, 2015 at 3:52
• @Marco Much faster if you use 2.. May 11, 2015 at 3:53
• @2012rcampion You code comes first , but it's a pity that it becomes slow when n is large. May 12, 2015 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). May 12, 2015 at 15:25
• N.B. These intersection points are related to what are now termed Padua points. Jul 19, 2016 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. May 11, 2015 at 3:56
• It's a pity that NSolve takes a very long time. May 11, 2015 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! :) May 11, 2015 at 12:31
• @Michael E2 I have posted my answer below,purely based on Solve. May 11, 2015 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.) May 11, 2015 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. May 11, 2015 at 17:24
• @Michael, thank you for the vote.
– kglr
May 11, 2015 at 17:39
• Glad to see my idea works. :) May 11, 2015 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, 2015 at 21:58
• @WateSoyan, it came up as one of the search results from ??**Intersections*.
– kglr
May 12, 2015 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.) May 11, 2015 at 13:52
• @Michael E2 Yeah,I forgot to check it. May 11, 2015 at 14:46