Sol LeWitt and his Incomplete Cubes

 The recent video on 3Blue1Brown is something I enjoyed after a long time in the channel. While any post, article or video about Burnside's Lemma is always a delight, the introduction of Sol Lewitt and his mathematical art made the video all the more satisfying to me.

Naturally, I too wanted to explore this idea. I was interested in couple of things - what if we include reflections as well? what about cubes that are unconnected?

While I tried to find the answers mathematically, I realized that I'm both stupid and lazy to figure these myself. I decided to take the easy way and use Mathematica to aid myself, and these results form the crux of this post.

While I understand that from an artist's perpective, Sol might have been interested about (only) connected cubes, we don't have to constrain ourselves with that. With a straightforward bruteforce approach, Mathematica gave me the following.

Count of Sol LeWitt's cubes by Connected Components and Number of Edges

Rows in the above matrix gives the number of connected components and the Columns represents the number of edges. For example, the top left value shows that there is exactly one cube that has $0$ connected components and $0$ edges - (obviously) the empty cube. Similarly, we can see that $23$ cubes have $6$ edges and $2$ connected components.

Some examples of Sol LeWitt's Incomplete cubes

The second row of the matrix (representing the number of cubes with exactly one connected component) was what Sol was interested. We see the two sums to $127$ of which $4$ are planar and $1$ is the complete cube giving the $122$ cubes that Sol found.

From here, I realized that the code I wrote should be easily extendable to identify reflections as well. With the assumption that the Octahedral group (3D equivalent to Dihedral group) is a combination of rotation and point reflection, I modified the same and got the following 'count' matrix.

Count of Sol LeWitt's cubes by Connected Components and Number of Edges with reflections identified

We see that there are a total of $144$ unique cubes after identifying rotations and reflections. Summing the second row, we see that there are $82$ cubes that have exactly one connected component of which we know $4$ are planar and one is the complete cube.

All $82$ connected cubes identifying reflections and rotations

First Version of Mathematica code
Clear["Global`*"];
vertexcoordinates = {1 -> {-1, -1, -1}, 2 -> {1, -1, -1}, 3 -> {1, 1, -1}, 4 -> {-1, 1, -1},
    5 -> {-1, -1, 1}, 6 -> {1, -1, 1}, 7 -> {1, 1, 1}, 8 -> {-1, 1, 1}};
cubegraph = {UndirectedEdge[1, 2], UndirectedEdge[2, 3], UndirectedEdge[1, 4], UndirectedEdge[3, 4],
    UndirectedEdge[1, 5], UndirectedEdge[2, 6], UndirectedEdge[5, 6], UndirectedEdge[3, 7],
    UndirectedEdge[6, 7], UndirectedEdge[4, 8], UndirectedEdge[5, 8], UndirectedEdge[7, 8]};
ReOrderGraph[g0_] := Module[{g = g0, temp},
    temp = Table[UndirectedEdge[Min[First[k], Last[k]], Max[First[k], Last[k]]], {k, g}];
    temp = SortBy[temp, {First, Last}];
    temp
];
cubegraph = ReOrderGraph[cubegraph];
cubegraphsubsets = Subsets[cubegraph];
cubegraphsubsets = Rest[cubegraphsubsets];
RotateGraph[g0_, p0_] := Module[{g = g0, p = p0, res},
    If[Length[g] <= 0, Return[{}];];
    res = Table[UndirectedEdge[p[[First[k]]], p[[Last[k]]]], {k, g}];
    res = ReOrderGraph[res];
    res
];
cuberotations = {
    {1, 2, 3, 4, 5, 6, 7, 8},
    {2, 3, 4, 1, 6, 7, 8, 5},
    {4, 1, 2, 3, 8, 5, 6, 7},
    {5, 6, 2, 1, 8, 7, 3, 4},
    {4, 3, 7, 8, 1, 2, 6, 5},
    {5, 1, 4, 8, 6, 2, 3, 7},
    {2, 6, 7, 3, 1, 5, 8, 4},
    {3, 4, 1, 2, 7, 8, 5, 6},
    {6, 5, 8, 7, 2, 1, 4, 3},
    {8, 7, 6, 5, 4, 3, 2, 1},
    {1, 4, 8, 5, 2, 3, 7, 6},
    {1, 5, 6, 2, 4, 8, 7, 3},
    {6, 2, 1, 5, 7, 3, 4, 8},
    {3, 2, 6, 7, 4, 1, 5, 8},
    {6, 7, 3, 2, 5, 8, 4, 1},
    {8, 4, 3, 7, 5, 1, 2, 6},
    {3, 7, 8, 4, 2, 6, 5, 1},
    {8, 5, 1, 4, 7, 6, 2, 3},
    {7, 3, 2, 6, 8, 4, 1, 5},
    {5, 8, 7, 6, 1, 4, 3, 2},
    {2, 1, 5, 6, 3, 4, 8, 7},
    {4, 8, 5, 1, 3, 7, 6, 2},
    {7, 6, 5, 8, 3, 2, 1, 4},
    {7, 8, 4, 3, 6, 5, 1, 2}
};
isomorphcounts = Association[{}];
While[Length[cubegraphsubsets] > 0,
    k = First[cubegraphsubsets];
    temp = Table[RotateGraph[k, j], {j, cuberotations}];
    temp = DeleteDuplicates[temp];
    isomorphcounts[k] = Length[temp];
    cubegraphsubsets = Complement[cubegraphsubsets, temp];
];
isomorphgraphedgelist = Keys[isomorphcounts];
graphcnt = Table[0, {5}, {13}];
graphcnt[[1, 1]] += 1;
Do[
    graphcnt[[1 + Length[ConnectedComponents[k]], 1 + EdgeCount[k]]] += 1;
, {k, isomorphgraphedgelist}
];
graphcnt // MatrixForm
(* isomorphgraphs = Table[GraphPlot3D[k, VertexCoordinates -> vertexcoordinates, ViewProjection -> "Orthographic", VertexLabels -> "Name", Boxed -> True, PlotRange -> {{-17/16, 17/16}, {-17/16, 17/16}, {-17/16, 17/16}}], {k, isomorphgraphedgelist}]; *)
(* isomorphgraphs = Table[GraphPlot3D[k, VertexCoordinates -> vertexcoordinates, ViewProjection -> "Orthographic", Boxed -> False, PlotRange -> {{-17/16, 17/16}, {-17/16, 17/16}, {-17/16, 17/16}}], {k, isomorphgraphedgelist}]; *)
isomorphgraphs = Table[GraphPlot3D[k, VertexCoordinates -> vertexcoordinates, ViewProjection -> "Orthographic", Boxed -> True, BoxStyle -> Directive[Dashed, White], PlotRange -> {{-17/16, 17/16}, {-17/16, 17/16}, {-17/16, 17/16}}], {k, isomorphgraphedgelist}];
(* GraphicsGrid[Partition[isomorphgraphs, 15]] *)
(* GraphicsGrid[{Take[isomorphgraphs, {200, 202}]}] *)
(* Table[GraphPlot3D[k, VertexCoordinates -> vertexcoordinates, ViewProjection -> "Orthographic", VertexLabels -> "Name", Boxed -> True, PlotRange -> {{-17/16, 17/16}, {-17/16, 17/16}, {-17/16, 17/16}}], {k, Take[isomorphgraphedgelist, 140]}] *)
(* GraphPlot3D[isomorphgraphedgelist[[3]], VertexCoordinates -> vertexcoordinates, ViewProjection -> "Orthographic", VertexLabels -> "Name"] *)
(* GraphicsGrid[Partition[Take[isomorphgraphs, {200, 210}], 4]] *)
Take[isomorphgraphs, {200, 211}]

Second Version:
Clear["Global`*"];
vertexcoordinates = {1 -> {-1, -1, -1}, 2 -> {1, -1, -1}, 3 -> {1, 1, -1}, 4 -> {-1, 1, -1},
    5 -> {-1, -1, 1}, 6 -> {1, -1, 1}, 7 -> {1, 1, 1}, 8 -> {-1, 1, 1}};
cubegraph = {UndirectedEdge[1, 2], UndirectedEdge[2, 3], UndirectedEdge[1, 4], UndirectedEdge[3, 4],
    UndirectedEdge[1, 5], UndirectedEdge[2, 6], UndirectedEdge[5, 6], UndirectedEdge[3, 7],
    UndirectedEdge[6, 7], UndirectedEdge[4, 8], UndirectedEdge[5, 8], UndirectedEdge[7, 8]};
ReOrderGraph[g0_] := Module[{g = g0, temp},
    temp = Table[UndirectedEdge[Min[First[k], Last[k]], Max[First[k], Last[k]]], {k, g}];
    temp = SortBy[temp, {First, Last}];
    temp
];
cubegraph = ReOrderGraph[cubegraph];
cubegraphsubsets = Subsets[cubegraph];
cubegraphsubsets = Rest[cubegraphsubsets];
RotateGraph[g0_, p0_] := Module[{g = g0, p = p0, res},
    If[Length[g] <= 0, Return[{}];];
    res = Table[UndirectedEdge[p[[First[k]]], p[[Last[k]]]], {k, g}];
    res = ReOrderGraph[res];
    res
];
cuberotations = {
    {1, 2, 3, 4, 5, 6, 7, 8},
    {2, 3, 4, 1, 6, 7, 8, 5},
    {4, 1, 2, 3, 8, 5, 6, 7},
    {5, 6, 2, 1, 8, 7, 3, 4},
    {4, 3, 7, 8, 1, 2, 6, 5},
    {5, 1, 4, 8, 6, 2, 3, 7},
    {2, 6, 7, 3, 1, 5, 8, 4},
    {3, 4, 1, 2, 7, 8, 5, 6},
    {6, 5, 8, 7, 2, 1, 4, 3},
    {8, 7, 6, 5, 4, 3, 2, 1},
    {1, 4, 8, 5, 2, 3, 7, 6},
    {1, 5, 6, 2, 4, 8, 7, 3},
    {6, 2, 1, 5, 7, 3, 4, 8},
    {3, 2, 6, 7, 4, 1, 5, 8},
    {6, 7, 3, 2, 5, 8, 4, 1},
    {8, 4, 3, 7, 5, 1, 2, 6},
    {3, 7, 8, 4, 2, 6, 5, 1},
    {8, 5, 1, 4, 7, 6, 2, 3},
    {7, 3, 2, 6, 8, 4, 1, 5},
    {5, 8, 7, 6, 1, 4, 3, 2},
    {2, 1, 5, 6, 3, 4, 8, 7},
    {4, 8, 5, 1, 3, 7, 6, 2},
    {7, 6, 5, 8, 3, 2, 1, 4},
    {7, 8, 4, 3, 6, 5, 1, 2}
};
ReflectVertex[k_] := {7, 8, 5, 6, 3, 4, 1, 2}[[k]];
cubereflections = Table[Map[ReflectVertex, k], {k, cuberotations}];
cuberotations = Join[cuberotations, cubereflections];
isomorphcounts = Association[{}];
While[Length[cubegraphsubsets] > 0,
    k = First[cubegraphsubsets];
    temp = Table[RotateGraph[k, j], {j, cuberotations}];
    temp = DeleteDuplicates[temp];
    isomorphcounts[k] = Length[temp];
    cubegraphsubsets = Complement[cubegraphsubsets, temp];
];
isomorphgraphedgelist = Keys[isomorphcounts];
graphcnt = Table[0, {5}, {13}];
graphcnt[[1, 1]] += 1;
Do[
    graphcnt[[1 + Length[ConnectedComponents[k]], 1 + EdgeCount[k]]] += 1;
, {k, isomorphgraphedgelist}
];
graphcnt // MatrixForm
lewitts = Select[isomorphgraphedgelist, Length[ConnectedComponents[#]] == 1 &];
isomorphgraphs = Table[GraphPlot3D[k, VertexCoordinates -> vertexcoordinates, ViewProjection -> "Orthographic", Boxed -> True, BoxStyle -> Directive[Dashed, White], PlotRange -> {{-17/16, 17/16}, {-17/16, 17/16}, {-17/16, 17/16}}], {k, lewitts}];
isomorphgraphs

Yours Aye

Me

Drain times with Toricelli's law

 It's easy to calculate the drain time of water in a cylindrical tank using Toricelli's law, $v=\sqrt{2gh}$. In fact, using a simple differential equation, we can see that

$\displaystyle \sqrt{\frac{h}{H}}=1-\frac{t}{t_H}$

where $H$ is the initial height and $t_H$ is the time it takes for the tank to drain a water column of height $H$.

I got curious about the case when there is a constant inflow of water. Then, we have the differential equation,

$a\sqrt{2gh_0} \,dt-a\sqrt{2gh}\,dt=A\,dh$

where $a$ is the area of discharge, $A$ the cylindrical area and $h_0$ the measure of constant inflow. Solving this gives,

$\displaystyle \frac{t}{t_0}=\sqrt{H/h_0}-\sqrt{h/h_0}+\ln\left(\frac{1-\sqrt{H/h_0}}{1-\sqrt{h/h_0}}\right)$

where $t_0$ is the time taken to drain a water column of height $h_0$.

Plotting this shows that irrespective of the height of the water column we start with, we always end up with a steady state water column of height $h_0$.

This brings us to the next case. Suppose we have a tank with an initial water column of height $H$ which drains water into a identical tank. We know that the first tank takes time $t_H$ to drain completely and therefore, will be restricting our analysis to this timeframe.

In this case, the differential equation (using the first equation of this post) becomes

$\displaystyle a\sqrt{2gH}\left(1-\frac{t}{t_H}\right)\,dt-a\sqrt{2gh}\,dt=A\,dh$

This can be simplified to

$\displaystyle 1 - \frac{t}{t_H}-\sqrt{\frac{h}{H}}=\frac{t_H}{2H}\frac{dh}{dt}$

We now makes substitutions, $x=1 - t/t_H$ and $y=\sqrt{h/H}$ so that the above equation becomes

$\displaystyle y-x=y\frac{dy}{dx}$

Solving this homogenous differential equation, we get

$\displaystyle\ln(y^2-yx+x^2)+\frac{2}{\sqrt{3}}\tan^{-1}\left(\frac{2}{\sqrt{3}}\frac{y}{x}-\frac{1}{\sqrt{3}}\right)=C$

Let's see some special cases of height of the water column in the second tank when the first tank runs out of water.

If the second tank is empty to start with, we have $h=0,t=0 \implies y=0,x=1$ as the initial conditions. With this we get,

$\displaystyle\ln(y^2-yx+x^2)+\frac{2}{\sqrt{3}}\tan^{-1}\left(\frac{2}{\sqrt{3}}\frac{y}{x}-\frac{1}{\sqrt{3}}\right)=-\frac{\pi}{3\sqrt{3}}$

Putting $x=0$ in the above (to find the water column at $t=t_H$), we see that $\displaystyle h(t_H) =\exp\left(-\frac{4\pi}{3\sqrt{3}}\right)H$

Similarly, if the two tanks start with the same level from the start, we use the initial conditions $h=H,t=0 \implies y=1,x=1$. This time we get,

$\displaystyle\ln(y^2-yx+x^2)+\frac{2}{\sqrt{3}}\tan^{-1}\left(\frac{2}{\sqrt{3}}\frac{y}{x}-\frac{1}{\sqrt{3}}\right)=\frac{\pi}{3\sqrt{3}}$

Now, if we put $x=0$ in the above we get, $\displaystyle h(t_H) =\exp\left(-\frac{2\pi}{3\sqrt{3}}\right)H$

Finally, note that if the initial water level in the second tank is lower than that of the first tank, it's water level rises initially, reaches a maximum and then falls off.

To find the maximum height it reaches, we need $dh=0 \implies dy=0 \implies y=x$.

If we let $y_{\text{max}}$ to denote the maximum $y$ and $y_1$ to denote the value of $y$ at $x=1$, we see that

$\displaystyle\ln(y^2-yx+x^2)+\frac{2}{\sqrt{3}}\tan^{-1}\left(\frac{2}{\sqrt{3}}\frac{y}{x}-\frac{1}{\sqrt{3}}\right)=\ln y_1^2+\frac{\pi}{\sqrt{3}}$

Putting $y=x$ in the above, we have

$\displaystyle \ln y_{\text{max}}^2+\frac{\pi}{3\sqrt{3}}=\ln y_1^2+\frac{\pi}{\sqrt{3}}$

This then shows that

$\displaystyle h_{\text{max}}=\exp\left(\frac{2\pi}{3\sqrt{3}}\right)h_1$

While it was already surprising that both $\pi$ and $e$ made an appearance in this problem, this result that the max. height and final height are in a constant ratio made it all the more satisfying for me.

Until then

Yours Aye

Me

Visual proof an integer sided $120^\circ$ triangle

 

Inspired by Proof Without Words: Integer Right Triangle Hypotenuses Without Pythagoras.

Until then

Yours Aye

Me

Sphere on a freely spinning turntable

Steve Mould somewhat recently published a video of a sphere rolling on a turntable rotating at constant rate and demonstrated that the path of the ball is, quite surprisingly, a circle. This demonstration was based on a 1978 paper titled Stable circular orbits of freely moving balls on rotating discs, Klaus Weltner. The same problem was also solved in Harvard Physics forum as Ball on turntable.

Its even more amazing that the path of a spherical ball rolling on a freely spinning turntable is a conic section as shown in A ball on a freely spinning turntable, Warren Weckesser [1]. Weckesser's request for a geometric derivation was answered about 20 years later in Comment on "A ball on a freely spinning turntable", Luis Rodriguez. A 2011 paper by Hector Munera titled A ball rolling on a freely spinning turntable: Insights from a solution in polar coordinates [2] further discusses some aspects that were not covered previously.

In this post, we aim to show 'the cone and the plane' of the conic section which has not been noted in any of the quoted papers. The key equations we need are (11) and (5) of [2] which we note here for ease.

$\displaystyle\Omega^2=\frac{\delta+\alpha r_0^2}{\delta+\alpha r^2}$ and $\mathbf{r}'=\alpha\Omega\mathbf{k}\times\mathbf{r}+\mathbf{c}$

For our purposes, we slightly alter these expression by differentiating the second and using some algebraic manipulations in the first.

$\displaystyle\Omega^2\left(\frac{\delta}{\alpha}+r^2\right)=\Omega_0^2\left(\frac{\delta}{\alpha}+r_0^2\right)$ and $\mathbf{r}''=\alpha\mathbf{k}\times(\Omega\mathbf{r})'$

Note that $\Omega_0=1$ by definition and therefore will be ignored at will from now on in this post.

To bring the cone out from hiding, we take the problem to the third dimension. To this end, we define

$\mathbf{p}=\mathbf{r}+\sqrt{\delta/\alpha} \text{ } \mathbf{k}$

Then, our two equations become

$\Omega p=\Omega_0p_0$  and $\mathbf{p}''=\alpha p_0 \mathbf{k}\times\widehat{\mathbf{p}}'$

where $p=|\mathbf{p}|$.

Because $\mathbf{p}$ and $\mathbf{r}$ are just a translation away, they have the same shape. Also, because of our definition of $\mathbf{p}$, we know that the curve defined by $\mathbf{p}$ lies on a plane perpendicular to the z-axis.

We also like to highlight that because of the definition of $\mathbf{p}$, the derivatives of $\mathbf{p}$ and $\mathbf{r}$ coincide and because $\mathbf{r} \cdot \mathbf{k}=0$, the derivatives of $\mathbf{p}$ are also perpendicular to $\mathbf{k}$.

Now we define a new vector $\mathbf{m}$ such that

$\mathbf{m}=\mathbf{p}\times\mathbf{p}'+\sqrt{\delta\alpha}p_0\text{ }\widehat{\mathbf{p}}$

Surprisingly, $\mathbf{m}$ is a constant of motion. This can be proved in multiple ways of which we show one here.

$\begin{align}(\mathbf{p}\times\mathbf{p}')' &= \mathbf{p}\times\mathbf{p}'' \\ &= \mathbf{p}\times (\alpha p_0 \mathbf{k} \times \widehat{\mathbf{p}}') \\ &= \alpha p_0 \bigl[ (\mathbf{p}\cdot\widehat{\mathbf{p}}')\mathbf{k}-(\mathbf{p}\cdot\mathbf{k})\widehat{\mathbf{p}}'\bigr] \\ &= -\sqrt{\delta\alpha} p_0\text{ }\widehat{\mathbf{p}}' \end{align}$

and therefore $\mathbf{m}'=(\mathbf{p}\times\mathbf{p}'+\sqrt{\delta\alpha} p_0\text{ }\widehat{\mathbf{p}})'=\mathbf{0}$.

We have used the fact that $\mathbf{q}\cdot\widehat{\mathbf{q}}'=0$ for any vector $\mathbf{q}$ because

$\displaystyle \widehat{\mathbf{q}}'=\frac{(\mathbf{q}\times\mathbf{q}')\times \mathbf{q}}{|\mathbf{q}|^3}$

Now that we know $\mathbf{m}$ is a constant vector, we can also easily see that $\widehat{\mathbf{p}}\cdot \mathbf{m}=\sqrt{\delta\alpha}p_0$, a constant which shows that the unit vectors of $\mathbf{p}$ have the same angle with $\mathbf{m}$ at all times.

Using these facts, we can finally see that the conic section that results from the path of the ball is the intersection of the $z=\sqrt{\delta/\alpha}$ plane and the cone whose vertex is the origin, $\mathbf{m}$ as its axis and $\mathbf{p} _0$ at its generator.

While the constancy of $\mathbf{m}$ could've been derived directly from (1), (4) and (5) of [2], taking the problem to the third dimension brings out the cone in the conic just beautifully.

Update (28-Mar-25)

The expression for $\mathbf{m}$ can rewritten as

$\mathbf{p}\times\mathbf{p}'=\mathbf{m}-\sqrt{\delta\alpha}p_0\text{ }\widehat{\mathbf{p}}$

As $\mathbf{m}$ is constant and $\widehat{\mathbf{p}}$ traces a circle with its axis along  $\mathbf{m}$, we note that the magnitude of the RHS of the above expression is a constant.

Therefore,

$|\mathbf{p}\times\mathbf{p}'|=|\mathbf{p}_0\times\mathbf{p}'_0|$

Since the cross product of two vectors gives the area of the parallelogram with the two vectors as its side, the curved surface of the cone $dA$ swept in time $dt$ is given by,

$\displaystyle \frac{dA}{dt}=\frac{1}{2}|\mathbf{p}\times\mathbf{p}'|=\frac{1}{2}|\mathbf{p}_0\times\mathbf{p}'_0|$

This shows the 'areal velocity' of $p$ is constant, much like that of the Kepler's laws of planetory motion. In particular, the orbital period $T$ of the sphere rolling on a freely spinning turntable is given by,

$\displaystyle T=\frac{A}{\frac{1}{2}|\mathbf{p}_0\times\mathbf{p}'_0|}$

The curved surface $A$ of the cone can be calculated using the formula given in Finkel's Solution Book as referenced in Problem 336, Mar 1914, American Mathematical Monthly.

To the best of my knowledge, this seems to be first time that an explicit expression for the time period of a sphere rolling on a freely spinning turntable!

Hope you enjoyed the discussion. See ya later.

Until then

Yours Aye

Me

An interesting Probability problem

I recently came across an interesting problem concerning the probability of a random sphere lying inside an unit ball in stackexchange.

The solution methodology was so fascinating that inspired me to come up with the following, almost similar, problem: Given a triangle and three points selected uniformly randomly from the triangle, what is the probability that the circle passing through those points lie completely inside the triangle?

We can characterise the three points the same way as in the SE post. Then,

$\displaystyle\mathbb{P}(\text{circle inside triangle})=\frac{1}{A^3}\iiint\limits_{\sqrt{x^2+y^2}+r \in \triangle}\,dx\,dy\,dr\int\limits_0^{2\pi}\int\limits_0^{2\pi}\int\limits_0^{2\pi}\frac{\partial (x_1,y_1,x_2,y_2,x_3,y_3)}{\partial (x,y,r,\phi_1,\phi_2,\phi_3)}\,d\phi_1\,d\phi_2\,d\phi_3$

where $A$ is the area of the given triangle $\triangle$ and the coordinates $(x,y)$ are w.r.t some arbitrary origin.

The integral involving $\phi$'s can be handled the same way as shown in the SE post which reduces the above expression to,

$\begin{align}\displaystyle\mathbb{P}(\text{circle inside triangle}) &= \frac{1}{A^3}\cdot (2\pi)^3\cdot 2 \cdot \frac{3}{2\pi}\iiint\limits_{\sqrt{x^2+y^2}+r \in \triangle}r^3 \,dx\,dy\,dr \\ &= \frac{24\pi^2}{A^3}\iiint\limits_{\sqrt{x^2+y^2}+r \in \triangle}r^3 \,dx\,dy\,dr \\ \end{align}$

If we now choose the incenter of the triangle to be the origin and $R$ to be inradius of the triangle, then the integral above w.r.t $\,dx\,dy$ is just the area of the set of points that are at the distance of $r$ inward from the triangle. But that is just the same triangle scaled such that its inradius is $R-r$. That is,

$\displaystyle \iint\limits_{\sqrt{x^2+y^2}+r \in \triangle}\,dx\,dy=\left(1-\frac{r}{R}\right)^2\cdot A$

Therefore,

$\begin{align}\displaystyle\mathbb{P}(\text{circle inside triangle}) &= \frac{24\pi^2}{A^3}\int\limits_0^R \left(1-\frac{r}{R}\right)^2\cdot A\,dr \\ \end{align}$

Evaluating this simple integral, we finally get,

$\displaystyle \mathbb{P}(\text{circle inside triangle}) = \frac{2}{5}\left(\frac{\pi R^2}{A}\right)^2$

A similar approach can be used to find the probability that a circle formed by choosing two points uniformly randomly inside a triangle completely lies within the triangle. In this case,

$\displaystyle \mathbb{P}(\text{circle inside triangle}) = \frac{2}{3}\frac{\pi R^2}{A}$

It should be easy to see from here that the above formulas, in fact, works for any tangential polygon. If we take it a step further, the same approach can be made to work for any closed convex curve. The only tricky part is the integral w.r.t '$\,dx\,dy$' which would now be the area enclosed by the 'inward' parallel curve of the given curve, which is not terribly difficult but isn't straightforward either.

An equally interesting but a tad more challenging problem is the following: Choose two points uniformly randomly inside a given triangle. If we now construct the square with the line joining the chosen points as the diagonal, what is the probability that the square lies within the triangle?

We'll see this problem in a separate post. If you figure out the answer, please share it with us in the comments section below.

Until then

Yours Aye

Me

Monkey Dog Problem - Grid Multiplication

The following problems was recently posted in the IIM-Cal puzzle group (and eventually in my puzzle group as well): Can we fill an $n \times n$ grid with the first $n^2$ integers such that the product of the numbers in the $k$th row equals that of the $k$th column for $1 \leq k \leq n$?

The fact that this is impossible for $n=2$ is easy to show by inspection. Vamshi Jandhyala in our group posted grids with the required property for $3 \times 3$ and $4 \times 4$ (and eventually for $10 \times 10$) using a Python non-linear solver. However, this approach ended up time-expensive for larger grids.

He also showed such grids are impossible whenever $n=9$ or $n>10$ using a nice argument involving the number of primes between $n^2/2$ and $n^2$.

Not wanting to get into the realm of non-linearity, I log'ed the entire grid so that multiplication becomes addition and the problem becomes linear. Note that this made an Integer Non-Linear problem to a Non-Integer Linear problem. Despite precision errors, this approach worked for $n=3,4$ but proved elusive for larger $n$.

Running out of options, I eventually realised that each grid that satisfies the required property can be seen only for a particular prime (their exponents in all the $n^2$ numbers to be specific) while ignoring others.

The top half of the above snapshot shows a solution to a $3 \times 3$ grid for a particular prime. The bottom half shows the same but in terms of the exponents of the prime.

The prime exponents form makes it clear that the problem can indeed be modelled as an Integer Linear Problem which can be solved extremely efficiently using free solvers in Python. In fact, doing so enables us to solve the problem for all the cases in a matter of seconds.

A Solution for $5 \times 5$ grid

A Solution for $6 \times 6$ grid

A Solution for $7 \times 7$ grid

A Solution for $8 \times 8$ grid

Python code used for solving this problem.

from ortools.linear_solver import pywraplp
import numpy as np


def main():
    # Data
    n = 3
    vals = []
    # primes = [2, 3, 5, 7, 11, 13, 17, 19, 23]
    # primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    # primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61]
    primes = [2, 3, 5, 7]
    primepi = len(primes)
    
    def vals(i0, k):
        i, cnt, p = i0, 0, primes[k]
        while i % p == 0:
            i //= p
            cnt += 1
        return cnt

    # Solver
    # Create the mip solver with the SCIP backend.
    solver = pywraplp.Solver.CreateSolver("SCIP")

    x = {}
    for i in range(n * n):
        for j in range(n * n):
            x[i, j] = solver.IntVar(0, 1, "")
    
    for i in range(n * n):
        solver.Add(solver.Sum([x[i, j] for j in range(n * n)]) == 1)
    
    for j in range(n * n):
        solver.Add(solver.Sum([x[i, j] for i in range(n * n)]) == 1)
    
    y = {}
    for i in range(n * n):
        for k in range(primepi):
            y[i, k] = solver.Sum([x[i, j] * vals(j + 1, k) for j in range(n * n)])

    
    rowsum = {}
    for i in range(n):
        for k in range(primepi):
            rowsum[i, k] = solver.Sum([y[j, k] for j in range(n * i, n * i + n)])
    
    colsum = {}
    for i in range(n):
        for k in range(primepi):
            colsum[i, k] = solver.Sum([y[j, k] for j in range(i, n * n, n)])
    
    for i in range(n):
        for k in range(primepi):
            solver.Add(rowsum[i, k] == colsum[i, k])
    
    solver.Minimize(1)
    status = solver.Solve()

    # Print solution.
    if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
        res = []
        for i in range(n * n):
            for j in range(n * n):
                if x[i, j].solution_value() > 0:
                    res += [j + 1]
        print(res)
    else:
        print("No solution found.")


if __name__ == "__main__":
    main()

That eventually settles the Monkey Dog Problem of Grid multiplication. There is also an interesting question of finding $n$ which has the lowest number of solutions if rotations and reflections are identified. That could be left as an exercise to the reader.

Hope you enjoyed this. See ya later.

Until then

Yours Aye

Me

Five ways of solving an Expected Value problem

 Let's say an Independent Aurelian is investigating butterflies on two types of traits, A and B. It is known that the probability of a butterfly exhibiting trait A is $p$ and that of trait B is $q$ (independent of each other). Note that this means the probability of a butterfly exhibiting both these traits is $pq$.

If she wants to study both these traits, how many butterflies, on average, would she need to catch? She can do this either by catching two butterflies each with exactly on of these traits or by catching a butterfly with both these traits.

The problem is simple and that subject of this post is look at five different methods of this problem. A key fact that will be of great use in all these methods is that for an event with success probability $r$, we would need, on average, $1/r$ trials to get the first success.

Simple equation

Let $N$ be the expected value we are looking for. We can condition this based on the nature of the first butterfly that is caught.

$\displaystyle N=(1-p)(1-q)\cdot(1+N)+(1-p)q\cdot \left(1+\frac{1}{p}\right)+p(1-q)\cdot\left(1+\frac{1}{q}\right)+pq\cdot 1$

Transfer Matrix

It is easy to setup the transfer matrix for this problem with the states being $\{\text{None, Only trait }A\text{, Only trait }B\text{, Both}\}$. Then the transfer matrix $\mathbf{T}$ is,

$\displaystyle \mathbf{T}=\begin{pmatrix}1-(1-p)(1-q)&p(1-q)&(1-p)q&pq \\ 0&1-q&0&q \\ 0&0&1-p&p \\ 0&0&0&1\end{pmatrix}$

It is then well known that the expected value is given by $\mathbf{1}(\mathbf{I}-\mathbf{T})^{-1}$ where $\mathbf{I}$ is the identity matrix and $\mathbf{1}$ is a row matrix of 1's.

Infinite Series

Let's say we've caught $n$ butterflies. The probability that none of them have trait A is $(1-p)^n$ which means that atleast one of them have trait A is $1-(1-p)^n$ and the same applies for trait B as well (with $p$ replaced with $q$).

$\mathbb{P}(\text{atleast of one each of both traits in }n\text{ trials})=(1-(1-p)^n)(1-(1-q)^n)$

With this, we see that $\mathbb{P}(N>n)=1-(1-(1-p)^n)(1-(1-q)^n)$

Then, $\displaystyle \mathbb{E}(N)=\sum_{n=0}^\infty \mathbb{P}(N>n)=\sum_{n=0}^\infty \Bigl (1-(1-(1-p)^n)(1-(1-q)^n) \Bigl )$

Story proof

We can answer this more directly. The expected number of butterflies needed can be broken into parts. First, expected number to get atleast one of the traits and second, the remaining value.

The probability of a butterfly exhibiting one of the traits is $1-(1-p)(1-q)=p(1-q)+(1-p)q+pq$. Therefore, we would need $1/(1-(1-p)(1-q))$ trials to get atleast of the traits. Now,

$\displaystyle \mathbb{P}(\text{Only trait A}|\text{atleast one of the traits})=\frac{p(1-q)}{p(1-q)+(1-p)q+pq}=\frac{p(1-q)}{1-(1-p)(1-q)}$

After having got trait A, we would be looking for trait B which we would get, on average, in $1/q$ trails. This approach applies for the other two cases as well. Therefore,

$\displaystyle\mathbb{E}(N)=\frac{1}{1-(1-p)(1-q)}+\frac{p(1-q)}{1-(1-p)(1-q)}\cdot\frac{1}{q}+\frac{(1-p)q}{1-(1-p)(1-q)}\cdot\frac{1}{p}+\frac{pq}{1-(1-p)(1-q)}\cdot 0$

(Shorter story proof using) Min-Max identity

It is easy to see that $x_1+x_2=\max (x_1,x_2)+\min (x_1,x_2)$ for any two real $x_1,x_2$.

Using this and the linearity of expectations, we have,

$\mathbb{E}(N)=\mathbb{E}(\text{trials to get trait A})+\mathbb{E}(\text{trials to get trait B})-\mathbb{E}(\text{trials to get atleast one of the traits})$

All of these are Geometric random variables with probabilities $p$, $q$ and $p+q-pq$. Therefore,

$\displaystyle \mathbb{E}(N)=\frac{1}{p}+\frac{1}{q}-\frac{1}{p+q-pq}$

Where's the fun in solving a problem if getting an answer is the only objective?

Yours Aye

Me