11
$\begingroup$

I have a cubic equation as below, which I am plotting:

Plot[(x + 1) (x - 1) (x - 2), {x, -2, 3}]

I like Mathematica to help me locate the position/equation of a circle which is on the lower part of this curve as shown, which would fall somewhere in between {x,-1,1}, which is tangent to the cubic at the 2 given points shown in red arrows. I like to find those points and the circle equation.

enter image description here

I tried generalizing the circle equation of (x-h)^2+(y-k)^2=r^2 and setting that in an equation to intersect the cube to be able to see the intersection points as such: (expanded the above circle equation and set a few conditions which I know are true)

eqn = FullSimplify[{x^3 - 2 x^2 - x + 2 == x^2 + h^2 - 2 x*h + y^2 + k^2 - 2 y*k - r^2, 1 > x > -1,0 <= r <= 1, h > 0}]

Table[FindInstance[eqn, {x, y, h, k, r}, Reals, 5] ] // N

But I cannot get a set of results from it...

However I had a bit more luck just Solving it straight as:

NSolve[x^3 - 2 x^2 - x + 2 == x^2 + h^2 - 2 x*h + y^2 + k^2 - 2 y*k - r^2]

which Produces:

{{h -> -1.90052, k -> -0.21556, r -> -0.918928, x -> -0.7017, 
  y -> 0.66687}, 

{h -> -0.314487, k -> 1.05695, r -> -0.172924, 
  x -> 0.621468, y -> 1.05907}, 

{h -> 1.63094, k -> 2.61779, 
  r -> 0.742119, x -> 2.24446, y -> 1.54014}}

First 2 appear to be valid intersect points for x and y But r is negative, which makes no sense. 3rd one has r>0 , but x is 2.2 which does not visually fall in the range that should be on the chart.

How can I get this to work and produce the results?

$\endgroup$
2
  • 1
    $\begingroup$ 1) start by adding AspectRatio -> Automatic to your Plot. You will see that the “circle” you drew would really have to be an ellipse. 2) wouldn’t there be an infinite number of such circles, with decreasing radii, until the two points of tangency “fuse” into one? Am I missing a constraint that would make it possible to identify a single circle? $\endgroup$
    – MarcoB
    Dec 27, 2020 at 7:00
  • $\begingroup$ Yes there definitely may be infinite. I am still trying to study answers below! $\endgroup$
    – Steve237
    Dec 27, 2020 at 16:07

2 Answers 2

13
$\begingroup$

First define a parametric curve and regions above and below the curve:

ClearAll[curve, disk, region1,  region2, radius]

curve[x_] := {x, (x + 1) (x - 1) (x - 2)}

region1 = ImplicitRegion[(x + 1) (x - 1) (x - 2) >= y, {{x, -20, 30}, {y, -80, 80}}];

region2 = ImplicitRegion[(x + 1) (x - 1) (x - 2) <= y, {{x, -20, 30}, {y, -80, 80}}];

and a disk with radius r tangent to the curve at curve[t]:

disk[dir_ : -1][t_, r_] := Module[{rr = Rationalize[r, 0],tt = Rationalize[t, 0]}, 
  Disk[curve[tt] + dir rr Cross@Normalize[curve'[tt]], rr]]

For given input region reg and parameter value t find the maximal radius r such that the disk with radius r tangent to the curve at curve[t] stays within reg:

radius[reg_, dir_ : -1][t_?NumericQ] :=
    NMaxValue[{r, RegionWithin[reg, disk[dir][t, r]]}, r]

Examples:

radius[region1][-1]
0.992403  
radius[region1][3/2]
 2.29173 
radius[region2, 1][1]
 0.555255  
pp = ParametricPlot[curve[x], {x, -2, 3}, ImageSize -> Medium, 
   PlotRange -> {{-5, 5}, {-6, 5}}];

frames = Table[Show[pp, 
    Graphics @ {PointSize[Large], Red, Point[curve[t]], FaceForm[], 
      EdgeForm[Red], disk[][t, radius[region1][t]]}], 
  {t, -3/2, 3/2, 1/100}];

Export["diskoncurve.gif", frames]

enter image description here

ParametricPlot[curve[x], {x, -2, 3}, AspectRatio -> Automatic, 
  Epilog -> {FaceForm[], 
    EdgeForm[Red], disk[][1 + 5/10, radius[region1][1 + 5/10]], 
    EdgeForm[Green], disk[][1, radius[region1][1]], 
    EdgeForm[Orange], disk[][1/5, radius[region1][1/5]], 
    EdgeForm[{Dashed, Red}], disk[1][1 + 5/10, radius[region2, 1][1 + 5/10]], 
    EdgeForm[{Dashed, Green}], disk[1][1, radius[region2, 1][1]], 
    EdgeForm[{Dashed, Orange}], disk[1][1/5, radius[region2, 1][1/5]]}]

enter image description here

$\endgroup$
6
  • 1
    $\begingroup$ Cool idea to use regions here! $\endgroup$
    – mgamer
    Dec 27, 2020 at 12:29
  • $\begingroup$ Thanks @kglr, like WOW! Even animations. So If I am interested to zoom in on the circle which gets split in half by the X axis, I assume I have to overlay that with my k being zero in: (x-h)^2+(y-k)^2=r^2, right? any quick way to set that? In your Parametric equation, where do you control that "K" value? .... $\endgroup$
    – Steve237
    Dec 27, 2020 at 16:12
  • 1
    $\begingroup$ @Steve237, if you want to add a constraint that the y-coordinate of the center of the disk is above a given threshold, say y0, Iyou can change the constraint in NMaximize to (curve[t] + dir r Cross @ Normalize[curve'[t]) >=y0 &&RegionWithin[reg, disk[dir][t, r]] (I think). $\endgroup$
    – kglr
    Dec 27, 2020 at 16:26
  • $\begingroup$ Thanks @kglr Lastly: I think you switched the circle/disk to parametric plot to make it easier to maneuver the equations and solve it, right? Would it have been possible to solve the same equation in the original cartesian format which I had given at start? Would the equation not be valid still? Why Mathematica did not solve that to get the same answer as yours? $\endgroup$
    – Steve237
    Dec 27, 2020 at 16:39
  • $\begingroup$ @Steve237, I used a parametric curve to represent the cubic, because finding the normal direction becomes simple (just use Cross curve'[x]) , and the centers of circles tangent to the curve at x lie on the line normal to the curve at x. When we pick a maximal radius subject to the constraint that the resulting disk lies within region with a smooth boundary we get tangency at a second point guaranteed (without having to compute the second tangency) . $\endgroup$
    – kglr
    Dec 27, 2020 at 16:54
6
$\begingroup$

Edit

f[x_] = (x + 1) (x - 1) (x - 2);
r[x_] = {x, f[x]};
eq1 = r[x1] + t1*Cross@Normalize[r'[x1]];
eq2 = r[x2] + t2*Cross@Normalize[r'[x2]];
plot = ParametricPlot[r[x], {x, -2, 3}, PlotStyle -> Blue];
xmin = x /. Solve[f'[x] == 0, x] // Last;
ymin = f[xmin];
xmax = x /. Solve[f'[x] == 0, x] // First;
ymax = f[xmax];
circles1 = 
  Graphics[Table[{Orange, Circle[eq1, Abs@t1]} /. 
     NMinimize[{0, t1^2 == t2^2, -2 < x1 < xmin < x2 < 3, eq1 == eq2, 
        eq1[[2]] == eq2[[2]] == d}, {x1, x2, t1, t2}][[2]], {d, 
     ymin + 1.5, ymin + 5, .2}]];
circles2 = 
  Graphics[Table[{Green, Circle[eq1, Abs@t1]} /. 
     NMinimize[{0, t1^2 == t2^2, -2 < x1 < xmax < x2 < 3, eq1 == eq2, 
        eq1[[2]] == eq2[[2]] == d}, {x1, x2, t1, t2}][[2]], {d, 
     ymax - 1, ymax - 5, -.2}]];
Show[plot, circles1, circles2, AspectRatio -> Automatic]

enter image description here

Original

Using the method from Find the equidistance curve between two curves

f[x_] = (x + 1) (x - 1) (x - 2);
r[x_] = {x, f[x]};
eq1 = r[x1] + t1*Cross@Normalize[r'[x1]];
eq2 = r[x2] + t2*Cross@Normalize[r'[x2]];
plot = ParametricPlot[r[x], {x, -2, 3}];
Manipulate[
 sol = NMinimize[{0, eq1 == eq2, t1^2 == t2^2, x2 - x1 == c}, {x1, x2,
     t1, t2}, Method -> Automatic];
 disk = Graphics[{Cyan, Opacity[0.2], EdgeForm[Red], 
    Disk[eq1, Abs@t1] /. Last[sol]}];
 Show[plot, disk, AspectRatio -> Automatic], {c, -2.5, 2.5, .01}, 
 ControlPlacement -> Bottom]

enter image description here

$\endgroup$

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.