# Finding solutions at non-differentiable points of a numerical function

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}] FindRoot[rd[x, Sqrt] == 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} and {0, -Sqrt}). I am at loss on how to do this for a black-box numerical function. Any suggestions?

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] ,
(#[] - #[]).(#[] - #[])/
( Norm[#[]-#[] ] Norm[#[]-#[] ] ) < .9 & ];
Graphics[{Line@points, {Red, PointSize[.02],
Point[corners[[All, 2]]]}}] 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}] 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}] 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.

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!