# Optimal triangle (minimal area) surrounding a disk

Using RegionWithin it's quite easy to evaluate a parametric triangle with maximal area inside the disk(constant):

dreieck = Triangle[{{x1, y1}, {x2, y2}, {x3, y3}}] (*parametric triangle*)
sol = NMaximize[{Area[dreieck], RegionWithin[Disk[], dreieck]},Flatten[{{x1, y1}, {x2, y2}, {x3, y3}}]]
Graphics[{FaceForm[Lighter[Gray]], Disk[], FaceForm[Darker[Green]],dreieck /. sol}]


Trying to solve a similar problem triangle with minimal area surrounding the disk makes difficulties

NMinimize[{Area[dreieck], RegionWithin[dreieck, Disk[] ]}, Flatten[{{x1, y1}, {x2, y2}, {x3, y3}}]]


RegionWithin[dreieck, Disk[] ] can be evaluated separately (huge expression) but NMinimize starts evaluation without result( EvaluationMonitor shows no progress)

Any idea why Mathematica doesn't solve? Thanks!

• LeafCount[RegionWithin[dreieck, Disk[], GenerateConditions -> False]] results in 377569 Mar 16 at 17:51
• @user64494 Thank you. That means I should first try to reduce the condition... Mar 16 at 18:00
• That reduction is doubtful. Here is a small part of the output -x2^3 y1 + 3 x2^2 x3 y1 - 3 x2 x3^2 y1 + x3^3 y1 + x1 x2^2 y2 - 2 x1 x2 x3 y2 - x2^2 x3 y2 + x1 x3^2 y2 + 2 x2 x3^2 y2 - x3^3 y2 - x2 y1 y2^2 + x3 y1 y2^2 + x1 y2^3 - x3 y2^3 - x1 x2^2 y3 + x2^3 y3 + 2 x1 x2 x3 y3 - 2 x2^2 x3 y3 - x1 x3^2 y3 + x2 x3^2 y3 + 2 x2 y1 y2 y3 - 2 x3 y1 y2 y3 - 3 x1 y2^2 y3 + x2 y2^2 y3 + 2 x3 y2^2 y3 - x2 y1 y3^2 + x3 y1 y3^2 + 3 x1 y2 y3^2 - 2 x2 y2 y3^2 - x3 y2 y3^2 - x1 y3^3 + x2 y3^3 >= 0. Mar 16 at 18:08
• Its LeafCount[%] performs 205. Mar 16 at 18:13
• @user64494 That looks much better,thanks! How did you calculte this simplification? I tried Reduce[RegionWithin[dreieck, Disk[]]] but Mathematica doesn't solve it in finite time. Mar 16 at 18:29

The reason might be that Mathematica fails to find an initial point that satisfies these complex conditions.

a = RegionWithin[dreieck, Disk[]];
FindInstance[a, {x1, y1, x2, y2, x3, y3}, Reals]


Also doesn't terminate. Mathematica doesn't seem to have an option for giving initial points. The more interesting question might be, what can you do solve this problem instead.

You simplify the search for feasible points a lot by choosing different variables. Note that any minimal triangle will touch the circle at 3 points and the lines will be tangents of the circle.

Each point on the unit circle is defined via one real number $$b, c, d$$ for the sake of rotational symmetry we can assume b == 0.

$$x \cos(b) - y \sin(b) + 1 = 0$$ defines a set of $$(x,y)$$ which is the tangent at the point $$(\sin(b),\cos(b))$$.

We can calculate the intersection of two tangents by finding a point $$(x,y)$$ that satisfies both equations.

Solve[{x Cos[b] - y Sin[b] + 1 == 0,  x Cos[c] - y Sin[c] + 1 == 0}, {x, y}]


Has the solution

{x -> -((Sin[b] - Sin[d])/(Cos[c] Sin[b] - Cos[b] Sin[d])),
y -> -((-Cot[b] + Cos[c] Csc[b])/(-Cos[c] + Cot[b] Sin[d]))}


Using those points we can now construct a triangle and calculate the area as you did above. However we also need a constraint that the center of the triangle (average of edge points) must lay in the circle. Futhermore only want to consider triangles bigger than the circle itself.

This results in this query:

NMinimize[{Area[
Triangle[{{-((Sin[b] - Sin[c])/(
Cos[c] Sin[b] - Cos[b] Sin[c])), -((-Cot[c] +
Cos[b] Csc[c])/(-Cos[b] + Cot[c] Sin[b]))}, {-((
Sin[b] - Sin[d])/(
Cos[d] Sin[b] - Cos[b] Sin[d])), -((-Cot[d] +
Cos[b] Csc[d])/(-Cos[b] + Cot[d] Sin[b]))}, {-((
Sin[c] - Sin[d])/(
Cos[d] Sin[c] - Cos[c] Sin[d])), -((-Cot[d] +
Cos[c] Csc[d])/(-Cos[c] + Cot[d] Sin[c]))}}]], ({-((
Sin[b] - Sin[c])/(
Cos[c] Sin[b] - Cos[b] Sin[c])), -((-Cot[c] +
Cos[b] Csc[c])/(-Cos[b] + Cot[c] Sin[b]))} + {-((
Sin[b] - Sin[d])/(
Cos[d] Sin[b] - Cos[b] Sin[d])), -((-Cot[d] +
Cos[b] Csc[d])/(-Cos[b] + Cot[d] Sin[b]))} + {-((
Sin[c] - Sin[d])/(
Cos[d] Sin[c] - Cos[c] Sin[d])), -((-Cot[d] +
Cos[c] Csc[d])/(-Cos[c] + Cot[d] Sin[c]))})/3 \[Element]
Disk[], Area[
Triangle[{{-((Sin[b] - Sin[c])/(
Cos[c] Sin[b] - Cos[b] Sin[c])), -((-Cot[c] +
Cos[b] Csc[c])/(-Cos[b] + Cot[c] Sin[b]))}, {-((
Sin[b] - Sin[d])/(
Cos[d] Sin[b] - Cos[b] Sin[d])), -((-Cot[d] +
Cos[b] Csc[d])/(-Cos[b] + Cot[d] Sin[b]))}, {-((
Sin[c] - Sin[d])/(
Cos[d] Sin[c] - Cos[c] Sin[d])), -((-Cot[d] +
Cos[c] Csc[d])/(-Cos[c] + Cot[d] Sin[c]))}}]] >=
Area[Disk[]], b == 0}, {b, c, d}]


which gives us

{5.19615, {b -> 0, c -> 2.0944, d -> -2.0944}}


which when plotted using RegionPlot looks like:

### An exact solution using Herons Formula

Wikipedia gives an exact formula for the length of the incircle based on the side length a,b,c. Herons formula gives the area of the triangle using the same information.

s := (a + b + c)/2; Heron = Sqrt[s*(s - a) (s - b) (s - c)];

Minimize[{Heron, 1 == Sqrt[(s - a) (s - b) (s - c)/s], a > 0, b > 0, c > 0}, {a, b, c}] Minimizes the area of a triangle with a circle of radius 1 inside. If the resulting expression is Fully Simplified the objective turns out to be 3 Sqrt[3] (matching our numerical calculation before) while all lengths are 2 Sqrt[2]

• Thank you for your answer. Trying NMinimize with the option "InitialPoints" -> { 0, 2, -2, -5/4, 2, -5/4}} also does not lead to success Mar 16 at 15:00
• I followed through with the strategy i gave. Mar 16 at 16:45
• Thank you for your solution. My purpose is to find a more general solution using RegionWithin applied to a parametric Region (for example a Polygon) and a constant region (for example a Polygon) . Mar 16 at 17:58
• I understand. I'm afraid Mathematica's generic implementation of RegionWithin results in to complicated terms for it to be solved by an generic algorithm. I would look into is solving the problem for a fixed outer shape, calculating the sensitives of the solution with regard to those parameters and then re optimize until you see no improvement. While this won't guarantee global optima, you will might get a solutions using NLP and maybe an active set approach. It is possible that more specialized algorithms exist for some shapes i'm not aware off. Maybe see Gilbert–Johnson–Keerthi distance Mar 16 at 20:56
• After some digging (for 2d) i would recommend this approach. You can transform math.stackexchange.com/questions/3191066/… math.stackexchange.com/a/863702/870285 an arbitrary quadrilateral into an unit square. Thus checking that all points are inside an quadrilateral is O(n) (compared to convex hull which has O(n*log(n)) operations. There are formulas to calculate the area of an quadrilateral. I hope this turns out to be feasible. Mar 16 at 22:07