17
$\begingroup$

I have three function $f(x)=x^x$, $g(x)=\ln(x)^{\ln(x)}$, $h(x)=\ln(x)^2$ and I want to find (numerically) $C$ the center of the circle tangent to the three curves.

I think that one way is to find two equidistance curves. Let's say $E_1$ the equidistance curve between $f(x)$ and $h(x)$, and $E_2$ the equidistance curve between $g(x)$ and $h(x)$. Then the point of intersection between $E_1$ and $E_2$ would be $C$ at the same Euclidean distance from the three curves.

I'm new to Mathematica, what are the methods to do so? Does the equidistance curve need to be traced point by point? If so how can I find such points?

To start I tried the following in order to find just one point in $E_1$ but it doesn't work because the first argument in NMinimize is not a number...

f[x_] := x^x;
h[x_] := Log[x]^2;
NSolve[First@NMinimize[EuclideanDistance[{0.9,y}, {x,f[x]}], x] == First@NMinimize[EuclideanDistance[{0.9,y}, {x,h[x]}], x], y]
$\endgroup$
1
  • 1
    $\begingroup$ Regarding my attempt that doesn't work, to those who face a similar situation, this is the proper way to handle an optimization nested inside another one. $\endgroup$ Dec 6 '20 at 15:37
16
$\begingroup$
f[x_] := x^x
g[x_] := Log[x]^Log[x]
h[x_] := Log[x]^2

plot = Plot[{f[x], g[x], h[x]}, {x, 0, E}, PlotRange -> {0, 1}]

enter image description here

Extract the three lines from plot and use them to construct three RegionDistance functions:

{linef, lineg, lineh} = Cases[plot, _Line, All];

{rdf, rdg, rdh} = RegionDistance /@ {linef, lineg, lineh};

ContourPlot to get the lines where rdf[{x, y}] == rdh[{x, y}] and rdg[{x, y}] == rdh[{x, y}]:

cp = ContourPlot[{ConditionalExpression[rdf[{x, y}] - rdh[{x, y}], 
      y <= f[x]] == 0, rdg[{x, y}] - rdh[{x, y}] == 0}, 
 {x, 0.434, 2.5}, {y, 0, 1}, ContourStyle -> {Red, Green}];

Show[plot, cp]

enter image description here

Find the intersection of the two contour lines:

center = First @ Graphics`Mesh`FindIntersections @ cp
{0.912936, 0.468359}

The distances to the three curves are

Through[{rdf, rdg, rdh} @ center]
 {0.374443, 0.374455, 0.374516}

Display the curves with the circle found, points of tangency and normals:

{ptf, ptg, pth} = RegionNearest[#, center] & /@ {linef, lineg, lineh};

Show[plot, Graphics[{AbsolutePointSize[10], Thick, 
   Red, Circle[intersection, rdf@center],
   Dashed, ColorData[97]@1, 
   InfiniteLine[{ptf, ptf + Cross @ {1, f'[ptf[[1]]]}}],
   ColorData[97]@2, InfiniteLine[{ptg, ptg + Cross @ {1, g'[ptg[[1]]]}}],
   ColorData[97]@3, InfiniteLine[{pth, pth + Cross @ {1, h'[pth[[1]]]}}],
   Black, Point@intersection, Black, Point @ {ptf, ptg, pth}}],
 AspectRatio -> Automatic, ImageSize -> 700]

enter image description here

A slower alternative: extract the two contourlines and find their intersection using RegionIntersection:

edlines = Cases[Normal @ cp, _Line, All];

center = (RegionIntersection @@ edlines)[[1, 1]]
 {0.912936, 0.468359}

An aside: To find the largest circle enclosed within the region (without the constraint that the circle touches all three curves) we can do

pw[x_] := Piecewise[{{f[x], x <= 1}}, g[x]];

plot = Plot[{f[x], g[x], h[x], pw[x], Min[pw[x], h[x]]}, {x, 0, E}, 
  Exclusions -> None, PlotRange -> {0, 1}, 
  Filling -> {4 -> {{3}, {None, Yellow}}}]

enter image description here

Extract the polygon from plot:

poly = First @ Cases[Normal[plot], _Polygon, All];

Use SignedRegionDistance with poly to get an objective function to be used with NMinimize:

srd[{x_, y_}]:= SignedRegionDistance[poly][{x, y}]

Find the coordinates within poly that maximizes the distance to the boundary of poly:

sol = NMinimize[{srd[{x, y}], {x, y} ∈ poly}, {x, y}]
{-0.394924,{x->0.991457,y->0.395402}}

Process to get the radius and center of largest circle within poly:

{radius, center} = {Abs @ #, {x, y} /. #2} & @@ sol
 {0.394924,{0.991457,0.395402}}

Display:

Graphics[{EdgeForm[Gray], Yellow, poly, Red, Circle[center, radius]}]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks, there are so many new things for me to learn from your answer! Also, how would you increase the precision of the solution? $\endgroup$ Dec 6 '20 at 15:27
17
$\begingroup$

Here we manual create the normal line of the three differentiable curves.

f[x_] := x^x;
g[x_] := Log[x]^Log[x]
h[x_] := Log[x]^2
eq1 = {x1, f[x1]} + t*Normalize[{-f'[x1], 1}];
eq2 = {x2, g[x2]} + t*Normalize[{-g'[x2], 1}];
eq3 = {x3, h[x3]} - t*Normalize[{-h'[x3], 1}];
sol = NMinimize[{0, (eq1 - eq2)^2 + (eq2 - eq3)^2 + (eq3 - eq1)^2 == 
    0}, {x1, x2, x3, t}]
disk = Disk[eq1, Abs[t]] /. Last[sol];
plot = Plot[{f[x], g[x], h[x]}, {x, 0, 3}, PlotRange -> {0, 3}, 
   AspectRatio -> Automatic, PlotRange -> All];
Show[plot, Graphics[{EdgeForm[Red], FaceForm[Cyan], disk}]]

{0., {x1 -> 0.733078, x2 -> 1.13371, x3 -> 0.583456, t -> -0.374512}}

enter image description here

$\endgroup$
7
  • $\begingroup$ Very clean and straightforward method! This is how I thought of doing it in the first place, but didn't know the full power of NMinimize to "converge" the three normal lines $\endgroup$ Dec 6 '20 at 15:31
  • $\begingroup$ Can you explain why NMinimize has to be "tricked" that way, can't we create a proper function to be minimized? $\endgroup$ Dec 6 '20 at 17:16
  • 1
    $\begingroup$ @DomenicoModica Yes,I am trying to do this for general three parametric curves. Maybe need more times to do this. The difficult come from how to determine the inside or outside of region. Here we use +t,+t and -t indicate the inside or outside of the normal vector. $\endgroup$
    – cvgmt
    Dec 6 '20 at 23:07
  • $\begingroup$ Oh, I'm interested in a general approach! By the way I've noted that if EuclideanDistance[eq1,eq2]+(*...*)==0 is used instead of (eq1-eq2)^2+(*...*) == 0 the solutions are very much different. The one obtained with EuclideanDistance has increasingly secure decimal digits with higher precision whereas the other can't secure not even the seventh decimal digit with WorkingPrecision 100 $\endgroup$ Dec 6 '20 at 23:27
  • $\begingroup$ @DomenicoModica The definition of EuclideanDistance or Norm or Normalize is use the Complex Function Abs instead of RealAbs,so I always trying to avoid use them. $\endgroup$
    – cvgmt
    Dec 6 '20 at 23:36
3
$\begingroup$

Here is a slightly more general attempt to solve the problem with NMinimize without "tricks". It is based on the idea used in @cvgmt interesting answer.

Because the "sign" of the normal directions isn't known one has to introduce three parameters t1,t2,t3 instead of t, all having the same absolut value:

f[x_] := x^x;
g[x_] := Log[x]^Log[x]
h[x_] := Log[x]^2
eq1 = {x1, f[x1]} + t1*Normalize[{-f'[x1], 1}];
eq2 = {x2, g[x2]} + t2*Normalize[{-g'[x2], 1}];
eq3 = {x3, h[x3]} + t3*Normalize[{-h'[x3], 1}];

sol = NMinimize[{  ((t1^2 - t2^2)^2 + (t2^2 - t3^2)^2 + (t3^2 -t1^2)^2),
eq1 == eq2 , eq3 == eq2 ,  eq1 == eq3 }, {x1, x2,x3, t1 , t2 , t3 } 
, Method ->"DifferentialEvolution" ]  
(*{1.63224*10^-13, {x1 -> 0.732916, x2 -> 1.13392, x3 -> 0.58397, 
t1 -> -0.374526, t2 -> -0.374526, t3 -> 0.374526}}*)

The result agrees with @cvgmt`s, perhaps the approach is a little bit more straighforward!

$\endgroup$
2
  • $\begingroup$ Nice! But it seems that increasing working precision with DifferentialEvolution, messes up the result... $\endgroup$ Dec 6 '20 at 19:44
  • $\begingroup$ @ DomenicoModica Yes, 'til now I didn't understand, why NMinimize- "Automatic"-modus doesn't work. Perhaps there exists a more robust formulation of the minimization problem? $\endgroup$ Dec 7 '20 at 7:10
3
$\begingroup$

For general parametric curves, NMinimize sometimes work. Here we give an example.

We change the sign of {t1,t2,t3}, for example t1 < 0, t2 < 0, t3 > 0.

But there still a result is not correct,that is t1 < 0, t2 < 0, t3 < 0 can not work,in this case, the result should be a large circle which just contain all the three circles.(The minimal circle contain all the circles).

And I have test some Method such as

Method ->Automatic;
Method -> "RandomSearch";
Method -> "DifferentialEvolution"
Method -> "SimulatedAnnealing";
Method -> "NelderMead"
f[s_] = {2 + Cos[s], 2 + Sin[s]};
g[s_] = {3 + Cos[s], 2 + Sin[s]};
h[s_] = {2 + Cos[s], 1 + Sin[s]};
eq1 = f[s1] + t1*Cross@Normalize[f'[s1], Sqrt[#.#] &];
eq2 = g[s2] + t2*Cross@Normalize[g'[s2], Sqrt[#.#] &];
eq3 = h[s3] + t3*Cross@Normalize[h'[s3], Sqrt[#.#] &];
sol = NMinimize[{0, Total[(eq1 - eq2)^2] == 0, 
    Total[(eq2 - eq3)^2] == 0, Total[(eq3 - eq1)^2] == 0, 
    t1^2 == t2^2 == t3^2, 0 <= s1 <= 2 π, 0 <= s2 <= 2 π, 
    0 <= s3 <= 
     2 π, {Sign[t1], Sign[t2], Sign[t3]} == {-1, -1, 1}}, {s1, s2,
     s3, t1, t2, t3}, Method -> Automatic];
disk = Disk[eq1, Abs[t1]] /. Last[sol];
plot = ParametricPlot[{f[s], g[s], h[s]}, {s, 0, 2 π}, 
   PlotRange -> {0, 4}, AspectRatio -> Automatic];
Show[plot, 
 Graphics[{EdgeForm[Red], FaceForm[Cyan], Opacity[0.2], disk}]]

enter image description here

With[{sign = #}, f[s_] = {2 + Cos[s], 2 + Sin[s]};
   g[s_] = {3 + Cos[s], 2 + Sin[s]};
   h[s_] = {2 + Cos[s], 1 + Sin[s]};
   eq1 = f[s1] + t1*Cross@Normalize[f'[s1], Sqrt[#.#] &];
   eq2 = g[s2] + t2*Cross@Normalize[g'[s2], Sqrt[#.#] &];
   eq3 = h[s3] + t3*Cross@Normalize[h'[s3], Sqrt[#.#] &];
   sol = NMinimize[{0, Total[(eq1 - eq2)^2] == 0, 
      Total[(eq2 - eq3)^2] == 0, Total[(eq3 - eq1)^2] == 0, 
      t1^2 == t2^2 == t3^2, 0 <= s1 <= 2 π, 0 <= s2 <= 2 π, 
      0 <= s3 <= 2 π, {Sign[t1], Sign[t2], Sign[t3]} == 
       sign}, {s1, s2, s3, t1, t2, t3}, Method -> Automatic];
   disk = Disk[eq1, Abs[t1]] /. Last[sol];
   plot = 
    ParametricPlot[{f[s], g[s], h[s]}, {s, 0, 2 π}, 
     PlotRange -> {0, 4}, AspectRatio -> Automatic];
   Show[plot, 
    Graphics[{EdgeForm[Red], FaceForm[Cyan], Opacity[0.2], 
      disk}]]] & /@ Tuples[{-1, 1}, 3]

enter image description here

Another example.

enter image description here

enter image description here

$\endgroup$
3
  • $\begingroup$ Nice! I went fancy and I'm working on a program that creates fractals finding the circle in every new area that gets create, see Apollonian gasket. Maybe next we can accorpate the two together! $\endgroup$ Dec 9 '20 at 14:59
  • $\begingroup$ For the case that can't work, it would be a second case of t1 > 0, t2 > 0, t3 > 0, because the engulfing circumference radius would traverse the three and not be outside. $\endgroup$ Dec 9 '20 at 15:36
  • $\begingroup$ @DomenicoModica Thanks for your advice. I will keep on study this question. $\endgroup$
    – cvgmt
    Dec 9 '20 at 15:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.