3
$\begingroup$

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

enter image description here

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!

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – march Oct 3 '16 at 18:13
  • $\begingroup$ @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. $\endgroup$ – Kagaratsch Oct 3 '16 at 18:40
  • $\begingroup$ related . mathematica.stackexchange.com/q/75352/2079 $\endgroup$ – george2079 Oct 3 '16 at 19:43
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Oct 4 '16 at 15:24
  • $\begingroup$ @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. $\endgroup$ – Kagaratsch Oct 6 '16 at 14:24
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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]

enter image description here

or random-ish points:

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

enter image description here

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

$\endgroup$
2
$\begingroup$

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]

enter image description here

$\endgroup$
1
$\begingroup$

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]

Mathematica graphics

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