7
$\begingroup$

Let's consider a 2D function we can evaluate numerically, but which is opaque to symbolic manipulation (this is just an example):

ClearAll[rd]
rd[x_?NumericQ, y_?NumericQ] := RegionDistance[Point[{{-1, 0}, {1, 0}}], {x, y}]

I want to find where this function has a value of, say, 2:

ContourPlot[rd[x, y] == 2, {x, -3, 3}, {y, -3, 3}]

enter image description here

FindRoot[rd[x, Sqrt[3]] == 2, {x, #}] & /@ {-3., 0.5, 2.47}

{{x -> -2.}, {x -> -5.47251*10^-17}, {x -> 2.}}

So far so good, I bet I can find lots of solutions to this problem numerically.

My problem is: I want to get only those (and preferably all - in case of my functions, a very limited set) solutions to rd[x, y] == 2 where the function is not differentiable. In the case of this example, points where circle arcs meet at an angle would be the sought answer.

If I could operate symbolically on the function, I would compare Limits of Grad of rd from different directions (x+, x-, y+, y-), construct a constraint that all of these can't be the same, and use Solve to find the symbolic result ({0, Sqrt[3]} and {0, -Sqrt[3]}). I am at loss on how to do this for a black-box numerical function. Any suggestions?

$\endgroup$
4
$\begingroup$

I expect you are after something more general, but for this case we can analyze the lines generated by contourplot: (I dont have RegionDistance but this should be the same )

rd[x_?NumericQ, y_?NumericQ] := 
    Min[EuclideanDistance[{x, y}, #] & /@ {{-1, 0}, {1, 0}}];
points = List @@ 
    First@Cases[
       Normal@First@
          Cases[ ContourPlot[
               rd[x, y] == 2, {x, -3, 3}, {y, -3, 3}] , _GraphicsComplex, 
               Infinity], _Line, Infinity] // First;
corners = 
  Select[ Partition[points, 3, 1] ,
    (#[[2]] - #[[1]]).(#[[3]] - #[[2]])/
      ( Norm[#[[2]]-#[[1]] ] Norm[#[[3]]-#[[2]] ] ) < .9 & ];
Graphics[{Line@points, {Red, PointSize[.02], 
     Point[corners[[All, 2]]]}}]

enter image description here

$\endgroup$
2
$\begingroup$

The curve in the Question, or more precisely the upper half of it, also can be obtained using InverseFunction

Plot[InverseFunction[rd, 2, 2][x, 2.], {x, -3, 3}, PlotRange -> {0, 3}, AxesLabel -> {x, y}]

enter image description here

Unfortunately, attempting to find the discontinuity in the slope by direct computation is both slow and noisy.

Plot[D[InverseFunction[rd, 2, 2][x, 2.], x] /. x -> z, {z, -1, 1}]

enter image description here

However, finding the discontinuity in the slope from a Table is reasonably effective.

rdl = Table[InverseFunction[rd, 2, 2][x, 2.], {x, -3, 3, .1}];
rdl2 = MapIndexed[(i = First@#2; If[i == 1 || i == Length[rdl], 0, 
   (rdl[[i - 1]] - 2 rdl[[i]] + rdl[[i + 1]])/(rdl[[i - 1]] + 2 rdl[[i]] + rdl[[i + 1]])]) &, rdl];
pos = Position[rdl2, Max[rdl2]];
Join[Extract[Range[-3, 3, .1], pos], Extract[rdl, pos]]
(* {0., 1.7320508075688774} *)

The best expression to be calculated in rdl2 depends on the character of the discontinuity sought.

$\endgroup$
2
$\begingroup$

This is a quick and dirty solution. You might want to optimize this and make sure it always works in case you use it seriously. But here's my idea. First get the points from your contour plot.

points = ContourPlot[rd[x, y] == 2, {x, -3, 3}, {y, -3, 3}][[1, 1]];

Now look at the change in the direction from one point to the next:

dist = Table[Total[(Normalize[points[[i]] - points[[i - 1]]] -Normalize[points[[i + 1]] - points[[i]]])^2], {i, 2,  Length[points] - 1}];

You are looking for the points with a pretty dramatic change, i.e.

points[[Flatten[Position[dist, #] & /@ Select[dist, # > 0.9 &]] + 1]]

Hope this helps!

$\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.