# Find convenient numeric solution to intersection contour?

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!

• You can extract the points from the plot and refine them using FindRoot. See, for instance, here. This is actually a pretty common method; there should be many examples of this on this very site. – march Oct 3 '16 at 18:13
• @march Thank you, this helped me alot! I adopted the solution in this link to my case, with a few tweaks. Since I need a lot of precision, I'll post my tweaks below in case if others might be interested. – Kagaratsch Oct 3 '16 at 18:40
• – george2079 Oct 3 '16 at 19:43
• Just how do you want to use the "numerical intersection"? I don't see the purpose in having a few scattered points to very high precision. (1) For plotting, increasing precision beyond a few digits does not improve the image. (2) Given a point close the curve, should one keep x or y fixed in trying to improve its position, or should one move parallel to the gradient or some other direction? Your example fixes x, but that seems arbitrary and inappropriate for all possible curves. – Michael E2 Oct 4 '16 at 15:24
• @MichaelE2 my intent was to find the slope of the curve close to the spot where it bends back. Having a few numeric points in this subregion allows to do a fit. – Kagaratsch Oct 6 '16 at 14:24

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.

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.

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]


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}]}}
*)