Monday, April 8, 2019

geometry - Minimum eccentricity of ellipses around another ellipse



Six circles can surround another circle of equal size, with each circle touching both the central circle and its two neighbouring outer circles.



For sufficiently eccentric ellipses, it is possible to do something similar, but with four ellipses surrounding another ellipse of the same size.




What is the minimum eccentricity of ellipses for which this is possible?


Answer



I'll concentrate on the symmetric case, as depicted here:



Illustration of the symmetric case





First, here are some numeric results:




\begin{align*}
a &= 1 \\
b &\approx 0.384369194474690789828391313191545078531 \\
c = \sqrt{1-\frac{b^2}{a^2}}
&\approx 0.923179463776614417385720356966316449484 \\
\varphi
&\approx 0.662140513907384715377580828031180874720 \,\text{rad} \\
\varphi
&\approx 37.93785689151653274240856839763865301015° \\

x &\approx 0.823320543612498211991170347260337013360 \\
y &\approx 0.685480094624740433428838597906083940535 \\
\alpha = \frac{b^2}{a^2} = 1-c^2
&\approx 0.147739677661122668796455680683973472705 \\
\beta = \tan\varphi
&\approx 0.779540453991525714639344299875163938705
\end{align*}



$a$ and $b$ are the length of the major and minor semiaxis. $c$ is the eccentricity, and $\pm\varphi$ is the angle of rotation between the central and one of the surrounding ellipses. $(\pm x,\pm y)$ is the center of one of the surrounding ellipses. $\alpha$ and $\beta$ are two parameters discussed below.




The above is the optimum I found using 150 iterations with 215 bits of numeric precision along the way. I later on verified that they are indeed accurate, using the exact algebraic computation described further below.



How to obtain them



Here is what I did: I defined all my geometry in terms of two parameters, $\alpha\approx 0.14778$ and $\beta\approx 0.77656$. The first defines the conic, as $\alpha x^2+y^2=1$. So contrary to my illustration and the figures above, I'm using $1$ as the length of the minor semiaxis and $\sqrt{\frac1\alpha}$ for the major internally, and the figures were scaled afterwards. I then interpret everything projectively, and use the second parameter to define $(1,\beta,0)^T$ and $(-\beta,1,0)^T$, two points at infinity which represent bundles of parallel lines, the two bundles orthogonal to one another. From these points I construct tangents to the central ellipse, and their point of intersection. These are the blue lines and red dots in the figure above. I then compute the transformation which maps these tangent lines to the coordinate lines, and apply it to the central ellipse to obtain the surrounding one which will certainly touch the coordinate axes.



Next I theoretically consider the two points of intersection which this ellipse might have with the central one. But more importantly, I consider the line joining these two points. The benefit here is on the one hand that this line is in fact easier to compute than the points of intersection, and that on the other hand the line will be real even if the ellipses no longer intersect in any real points.



Now the task can be formulated in terms of the position of this line wrt. the central ellipse: it must touch the central ellipse, and any change of $\beta$ in either direction must cause it to intersect the central ellipse, so it must be optimal in a certain sense. The adjoint of the matrix of the central ellipse will represent a quadratic form, and plugging the line into that form will result in a single number. According to the above conditions, this number must be both zero and optimal. In my choice of signs it must be maximal, but sign changes will occur easily.




Now comes the numeric refinement via bisection. Each loop adjusts first $\beta$ then $\alpha$. $\beta$ is adjusted towards optimality, while $\alpha$ is adjusted towards a root, where the result is zero. So in other words, I choose the angle such that the line is as close to the center of the central ellipse, then choose the eccentricity so that for that angle it still only touches the central ellipse. The adjustment is a pretty stupid bisection, which works since I always stay close to the actual solution and won't have to worry about other optima. Speed of convergence might be improved, but I'd rather keep the simplicity of my approach and wait a bit longer.



I've uploaded the sage code for my computations here. If anyone were to find symbolic expression for the minimal eccentricity, then I guess one could adjust my code to verify that.





As achille hui pointed out in a comment below, a post by Jim Ferry already contains an exact solution, i.e. one in terms of algebraic numbers. I've now been able to reproduce that result. The key idea is to work less in terms of construction steps, and instead more in terms of a characterization of desired properties. This allows one to avoid using square roots, so one can express everything (except angles) as polynomials.



Formulating the touching conditions




I started with two ellipses, one with horizontal semimajor axis of length $\sqrt{\frac1\alpha}$ and vertical semiminor axis of $1$, the other rotated by $\arctan\beta$ and translated by $(x,y)$. I formulated both as $3\times 3$ matrices in the polynomial ring $\mathbb Q[\alpha,\beta,x,y]$. The rotation by $\arctan\beta$ would usually involve a square root, but that multiplication is applied to the matrix of the conic from both sides (i.e. the conic is conjugated by it), and if done right (i.e. the square root placed in the homogenization term, in the bottom right corner of the transformation matrix), the two square roots will multiply so there is no square root left in the result. Here are my two ellipses. (In case you wonder, $E_2$ is the rotated but not translated step in between.)



$$E_1=\begin{pmatrix}
\alpha & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & -1
\end{pmatrix}
\\
E_3=
{\scriptsize\begin{pmatrix}

\beta^{2} + \alpha & \alpha \beta - \beta & - \beta^{2} x - \alpha \beta y - \alpha x + \beta y \\
\alpha \beta - \beta & \alpha \beta^{2} + 1 & - \alpha \beta^{2} y - \alpha \beta x + \beta x - y \\
- \beta^{2} x - \alpha \beta y - \alpha x + \beta y & - \alpha \beta^{2} y - \alpha \beta x + \beta x - y & \alpha \beta^{2} y^{2} + \beta^{2} x^{2} + 2 \alpha \beta x y + \alpha x^{2} - 2 \beta x y - \beta^{2} + y^{2} - 1
\end{pmatrix}}
$$



(Note that I guess I might have rotated my ellipse in the opposite direction, but that's irrelevant for now.) The fact that a given line touches a given ellipse can be computed via the dual quadric, which can be easily (i.e. without divisions or square roots) computed as the adjoint matrix. So the rotated ellipse will touch the axes of the coordinate system iff



$$
(1,0,0)\cdot\operatorname{adj}(E_3)\cdot\begin{pmatrix}1\\0\\0\end{pmatrix}=0

\qquad
(0,1,0)\cdot\operatorname{adj}(E_3)\cdot\begin{pmatrix}0\\1\\0\end{pmatrix}=0
$$



Next comes the condition that the two ellipses touch. When intersecting two conics, one approach starts by looking for roots of $\det(\mu E_1+E_3)=0$. These values of $\mu$ result in a degenerate conic $\mu E_1+E_3$ consisting of a pair of lines and passing through the points of intersection. If two points of intersection coincide, i.e. if you have a point of higher multiplicity, then this third degree polynomial in $\mu$ must have a zero of higher multiplicity as well since some of the degenerate conics will coincide. So to find out the touching condition, we look at the discriminant of $\det(\mu E_1+E_3)$ with respect to $\mu$. This is the third condition.



Computing the solution



Now one can use resultants to eliminate variables, namely to eliminate $x$ and $y$, thereby combining the three conditions into a single polynomial in $a$ and $b$. For this step, Mr. Lewis used a method he developed himself and called “Dixon-EDF” (for Early Detection of Factors). But although his computer algebra system Fermat comes with a file implementing that method, I haven't yet worked out how to adjust the various knobs to get it to solve a different set of equations, and the user interface kept annoying me. So I sticked with sage, and the naive computation took quite some time, about 10 minutes for the computation up to the final result but mostly for this step here. In the end, the resulting polynomial had 2069 terms. I factorized it and cheated a bit by not considering all potential factors, but only the one which had the smallest absolute value for the approximate numeric values obtained above.




We are looking for the minimal value of $\alpha$ satisfying all requirements. This optimality constraint corresponds to a root of the condition we just found which is a double root with respect to rotation. After all, we want an intersecting situation on both sides of the solution $\beta$. So I computed the discriminant of this polynomial with respect to $\beta$, again factorized and used the numerical approximation to find the correct factor, which is this:



$$151632\alpha^{12} - 556632\alpha^{11} - 4029183\alpha^{10} + 5710568\alpha^9 + 2456300\alpha^8 - 8614032\alpha^7 - 40073338\alpha^6 - 8614032\alpha^5 + 2456300\alpha^4 + 5710568\alpha^3 - 4029183\alpha^2 - 556632\alpha + 151632$$



This is essentially the same polynomial Mr. Ferry used in his approach as well, except that since my $\alpha$ relates to the squared ratio of axes, my exponents are half of his. He describes how one can use the symmetry of the polynomial (he called it a “palindrome”) to use a degree 6 polynomial from which everything else can be computed. But for me this was enough, so I obtained real roots, expressed as algebraic numbers in sage, and chose the one closest to the approximate value.



Plugging that value into the polynomial before, the relevant factor of the big polynomial in $\alpha$ and $\beta$, I obtained a value for $\beta$. Its minimal polynomial is particularly nice and easy:



$$16\beta^{12} - 100\beta^{10} + 615\beta^8 - 2468\beta^6 + 3186\beta^4 - 1728\beta^2 + 351$$




(Here you notice that $\beta$ is defined only up to sign, which is the reason why you don't have to worry too much about which direction is which when rotating the ellipse.)
Having found $\alpha$ and $\beta$, I could compute all the rest, and in fact verify that the digits given above were indeed accurate. I will update my uploaded code to include the code and output for this computation.


No comments:

Post a Comment

analysis - Injection, making bijection

I have injection $f \colon A \rightarrow B$ and I want to get bijection. Can I just resting codomain to $f(A)$? I know that every function i...