The goal is to find the largest ellipse (with given ratio of axes), centered at a given point and with a given orientatiion, that fits inside a specified non-convex oval.
equation of oval
In[1]:= oval[{x_, y_}] = ((x - 1)^2 + y^2) ((x + 1)^2 + y^2) - (21/20)^4;
derive equation of an ellipse with axes "a" and "b", centered at { xc, yc }
with major axis making angle \[Theta] with x-axis.
In[2]:= Thread[{xel, yel} =
DiagonalMatrix[{a, b}^-1].RotationMatrix[-\[Theta]].{x - xc, y - yc}];
In[3]:= eleq[{{a_, b_}, xc_, yc_, \[Theta]_}, {x_, y_}] = xel^2 + yel^2 - 1;
symbolically find largest ellipse with axes "a" and "a/2",
centered at { 1, (1/5) }
oriented at \[Pi]/3
RegionWithin takes about 1 minute.
In[4]:= AbsoluteTiming[
RegionWithin[ImplicitRegion[oval[{x, y}] <= 0, {x, y}],
ImplicitRegion[eleq[{{a, a/2}, 1, 1/5, \[Pi]/3}, {x, y}] <= 0, {x, y}],
GenerateConditions -> True] // N]
Out[4]= {63.4787, 0. < a <= 0.315686 || -0.315686 <= a < 0.}
Calculating it numerically with a Lagrange multiplier
and NSolve takes about 1 second.
The desired answer is the one with the smallest value of a.
In[5]:= AbsoluteTiming[
sln = NSolve[{a >= 0, oval[{x, y}] == 0,
eleq[{{a, a/2}, 1, 1/5, \[Pi]/3}, {x, y}] == 0,
Sequence @@
D[oval[{x,
y}] == \[Lambda] eleq[{{a, a/2}, 1, 1/5, \[Pi]/3}, {x, y}], {{x,
y}}]}, {a, x, y, \[Lambda]}, Reals]]
Out[5]= {0.808361, {{a -> 0.315686, \[Lambda] -> 0.869176, x -> 1.15549,
y -> 0.474695}, {a -> 4.34436, \[Lambda] -> 7.03937, x -> -1.41308,
y -> 0.191823}, {a -> 0.817698, \[Lambda] -> 1.46331, x -> 0.654269,
y -> -0.531984}, {a -> 1.14366, \[Lambda] -> 1.77874, x -> 1.34316,
y -> -0.315728}}}
eliminating the Lagrange multiplier before solving speeds up the calculation.
In[6]:= AbsoluteTiming[
sln = NSolve[{a >= 0, oval[{x, y}] == 0,
eleq[{{a, a/2}, 1, 1/5, \[Pi]/3}, {x, y}] == 0,
Eliminate[
D[oval[{x,
y}] == \[Lambda] eleq[{{a, a/2}, 1, 1/5, \[Pi]/3}, {x, y}], {{x,
y}}], \[Lambda]]}, {a, x, y}, Reals]]
Out[6]= {0.0641214, {{x -> 1.15549, y -> 0.474695, a -> 0.315686}, {x -> -1.41308,
y -> 0.191823, a -> 4.34436}, {x -> 1.34316, y -> -0.315728,
a -> 1.14366}, {x -> 0.654269, y -> -0.531984, a -> 0.817698}}}
Plotting all the results show that the curves are tangent at the intersection point.