Find convenient numeric solution to intersection contour?

Question

Consider the following contour of intersections between two analytic expressions:

ContourPlot[ Sinh[x^2] (x^2 + Cosh[y^2] y^2) == y^2, {x, -1, 1}, {y, -2, 2}]

The above nicely visualizes the intersection set, and also implies that the implementation of ContourPlot contains a very efficient algorithm that evaluates the intersection region numerically before plotting it. If I try to find numeric solutions to this problem, I fix i.e. x=1/2 and use something like FindRoot to get y->0.290862, or x=1/3 and y->0.117988, always one point at a time. This is a very inconvenient way to approach this problem, since it also tends to jump whenever several intersections exist and numerically break down whenever the slope becomes large.

Surely, there must be a better way to evaluate the intersection region numerically? Is there a function in Mathematica that would return a convenient representation of the numeric intersection region? Or maybe one can write such a function efficiently? Thanks for any suggestion!

Answer 1

As pointed out in the comments by march, one solution is to extract the points from the contour plot and refine them, as described i.e. here.

In my case I need more than machine precision, so here is what I do if I need the points to be 300 digits precise (mostly following the link above):

f[x_,y_]:=Sinh[x^2] (x^2 + Cosh[y^2] y^2) - y^2;
precis=300;
cp = ContourPlot[f[x,y] == 0, {x, -1, 1}, {y, -2, 2}];
cpRef = {#1, y /. FindRoot[ f[x,y] /. x :> #1, {y, #2}, WorkingPrecision -> precis]} & @@@ SetPrecision[cp[[1, 1]], precis+50];

Note how the original points from the plot are set to have a larger precision. We can check that the results indeed are very good:

Sum[f[x, y] /. x :> cpRef[[i, 1]] /. y :> cpRef[[i, 2]], {i, 1, cpRef // Length}]

0.*10^-297

It is also worth noting that one can increase the number of points in the ContourPlot with the PlotPoints-># option, as needed.

Answer 2

you can do some things with ImplicitRegion

r = ImplicitRegion[
  Sinh[x^2] (x^2 + Cosh[y^2] y^2) == y^2 , {{x, -1, 1}, {y, -4, 4}}]
d = DiscretizeRegion[r, MaxCellMeasure -> .01]
Graphics@Point@MeshCoordinates[d]

or random-ish points:

 Graphics@Point[{x, y} /. 
    FindInstance[Element[{x, y}, r], {x, y}, Reals, 100]]

ImplicitRegion tends to be a bit flakey however. Simply changing the domain sometimes causes this to fail. The ContourPlot approach is probably more robust.

Answer 3

NSolve can find all of the roots within a restricted region. Set the WorkingPrecision to whatever is required.

eqn = Sinh[x^2] (x^2 + Cosh[y^2] y^2) == y^2;

You can fix x and solve for y

data1 = Select[
   Table[{x, y} /. NSolve[{eqn, -3 <= y <= 3}, y], {x, -0.71, 0.71, .025}] //
          Flatten // Partition[#, 2] &,
   FreeQ[#, y] &];

and/or you can fix y and solve for x

data2 = Select[
   Table[{x, y} /. NSolve[{eqn, -0.71 <= x <= 0.71}, x], {y, -3, 3, .025}] //
          Flatten // Partition[#, 2] &,
   FreeQ[#, x] &];

ListPlot[{data1, data2},
 ImageSize -> 300,
 Frame -> True,
 Axes -> False,
 PlotRange -> {{-1, 1}, {-3, 3}},
 AspectRatio -> 2]

Answer 4

Solve returns exact solutions, which can be evaluated to arbitrary precision.

Block[{x = 1/2},
 Solve[Sinh[x^2] (x^2 + Cosh[y^2] y^2) == y^2, y, Reals]]
(*
  {{y -> Root[{Sinh[1/4] - 4 #1^2 + 2 E^-#1^2 Sinh[1/4] #1^2 + 
          2 E^#1^2 Sinh[1/4] #1^2 &, -1.42120738768396141533}]},
   {y -> Root[{Sinh[1/4] - 4 #1^2 + 2 E^-#1^2 Sinh[1/4] #1^2 + 
          2 E^#1^2 Sinh[1/4] #1^2 &, -0.29086218728559318745}]},
   {y -> Root[{Sinh[1/4] - 4 #1^2 + 2 E^-#1^2 Sinh[1/4] #1^2 + 
          2 E^#1^2 Sinh[1/4] #1^2 &, 0.29086218728559318745}]},
   {y -> Root[{Sinh[1/4] - 4 #1^2 + 2 E^-#1^2 Sinh[1/4] #1^2 + 
          2 E^#1^2 Sinh[1/4] #1^2 &, 1.42120738768396141533}]}}
*)

To get exact values from Solve, you need an exact value for x, which can achieved with SetPrecision[] or Rationalize[].

Plot[With[{x = SetPrecision[x0, Infinity]}, 
  y /. Solve[Sinh[x^2] (x^2 + Cosh[y^2] y^2) == y^2, y, Reals] /.
   {y -> {Indeterminate, Indeterminate, Indeterminate, Indeterminate}}],
 {x0, -1, 1}, Evaluated -> False, AspectRatio -> 1]

If you need to focus on one branch, you may include inequalities:

Block[{x = 1/2},
 Solve[Sinh[x^2] (x^2 + Cosh[y^2] y^2) == y^2 && y > 874/1000, y, Reals]]
(*
  {{y -> Root[{Sinh[1/4] - 4 #1^2 + 2 E^-#1^2 Sinh[1/4] #1^2 + 
          2 E^#1^2 Sinh[1/4] #1^2 &, 1.42120738768396141533}]}}
*)