0
$\begingroup$

I have this function and I want to see where it is zero. $$\frac{1}{16} \left(\sinh (\pi x) \left(64 \left(x^2-4\right) \cosh \left(\frac{2 \pi x}{3}\right) \cos (y)+\left(x^2+4\right)^2+256 x \sinh \left(\frac{2 \pi x}{3}\right) \sin (y)\right)+\left(x^2-12\right)^2 \sinh \left(\frac{7 \pi x}{3}\right)-2 \left(x^2+4\right)^2 \sinh \left(\frac{5 \pi x}{3}\right)\right)+2 \left(x^2-4\right) \sinh \left(\frac{\pi x}{3}\right)$$ I use ContourPlot

f[x_, y_] := 2 (-4 + x^2) Sinh[(π x)/3] + 1/16 (((4 + x^2)^2 + 64 (-4 + x^2) Cos[y] Cosh[(2 π x)/3] + 256 x Sin[y] Sinh[(2 π x)/3]) Sinh[π x] - 2 (4 + x^2)^2 Sinh[(5 π x)/3] + (-12 + x^2)^2 Sinh[( 7 π x)/3]); ContourPlot[ f[x, y] == 0, {x, 3.465728, 3.465729}, {y, 1.046786, 1.046795}, PlotPoints -> 500]

and I obtain this plot

Now, my question is that can I trust this plot and conclude that the curves do not cross?

Or, I should increase the precision of the plot? And if so, how can I ask Mathematica to give higher precision for the axis in ContourPlot?

Improve this question
$\endgroup$

4 Answers 4

Reset to default
1
$\begingroup$

f[x, y] == 0 can be separated in lhs[y] == rhs[x]. This shows, rhs gets immaginary in a small range of x. So curves do not intersect.

f[x_, y_] = 2 (-4 + x^2) Sinh[(\[Pi] x)/3] + 1/16 (((4 + x^2)^2 + 64 (-4 + x^2) Cos[y] Cosh[(2 \[Pi] x)/3] + 256 x Sin[y] Sinh[(2 \[Pi] x)/3]) Sinh[\[Pi] x] - 2 (4 + x^2)^2 Sinh[(5 \[Pi] x)/ 3] + (-12 + x^2)^2 Sinh[(7 \[Pi] x)/3]); eq = Equal @@@ (First@ Solve[f[x, y] == 0 && Sin[y]^2 + Cos[y]^2 == 1, Cos[y], Sin[y]]) // First (* Cos[y] == (-65536 Cosh[(2 \[Pi] x)/3] Sinh[(\[Pi] x)/ 3] Sinh[\[Pi] x] + 32768 x^2 Cosh[(2 \[Pi] x)/3] Sinh[(\[Pi] x)/3] Sinh[\[Pi] x] - 4096 x^4 Cosh[(2 \[Pi] x)/3] Sinh[(\[Pi] x)/3] Sinh[\[Pi] x] + 8192 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 + 2048 x^2 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 - 512 x^4 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 - 128 x^6 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 - 16384 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] - 4096 x^2 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] + 1024 x^4 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] + 256 x^6 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] + 73728 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] - 30720 x^2 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] + 3584 x^4 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] - 128 x^6 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/ 3] - \[Sqrt]((65536 Cosh[(2 \[Pi] x)/3] Sinh[(\[Pi] x)/ 3] Sinh[\[Pi] x] - 32768 x^2 Cosh[(2 \[Pi] x)/3] Sinh[(\[Pi] x)/ 3] Sinh[\[Pi] x] + 4096 x^4 Cosh[(2 \[Pi] x)/3] Sinh[(\[Pi] x)/ 3] Sinh[\[Pi] x] - 8192 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 - 2048 x^2 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 + 512 x^4 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 + 128 x^6 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x]^2 + 16384 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/ 3] + 4096 x^2 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[( 5 \[Pi] x)/3] - 1024 x^4 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/ 3] - 256 x^6 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[( 5 \[Pi] x)/3] - 73728 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/ 3] + 30720 x^2 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[( 7 \[Pi] x)/3] - 3584 x^4 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/ 3] + 128 x^6 Cosh[(2 \[Pi] x)/3] Sinh[\[Pi] x] Sinh[( 7 \[Pi] x)/3])^2 - 4 (65536 Cosh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2 - 32768 x^2 Cosh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2 + 4096 x^4 Cosh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2 + 65536 x^2 Sinh[(2 \[Pi] x)/ 3]^2 Sinh[\[Pi] x]^2) (16384 Sinh[(\[Pi] x)/3]^2 - 8192 x^2 Sinh[(\[Pi] x)/3]^2 + 1024 x^4 Sinh[(\[Pi] x)/3]^2 - 4096 Sinh[(\[Pi] x)/3] Sinh[\[Pi] x] - 1024 x^2 Sinh[(\[Pi] x)/3] Sinh[\[Pi] x] + 256 x^4 Sinh[(\[Pi] x)/3] Sinh[\[Pi] x] + 64 x^6 Sinh[(\[Pi] x)/3] Sinh[\[Pi] x] + 256 Sinh[\[Pi] x]^2 + 256 x^2 Sinh[\[Pi] x]^2 + 96 x^4 Sinh[\[Pi] x]^2 + 16 x^6 Sinh[\[Pi] x]^2 + x^8 Sinh[\[Pi] x]^2 - 65536 x^2 Sinh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2 + 8192 Sinh[(\[Pi] x)/3] Sinh[(5 \[Pi] x)/3] + 2048 x^2 Sinh[(\[Pi] x)/3] Sinh[(5 \[Pi] x)/3] - 512 x^4 Sinh[(\[Pi] x)/3] Sinh[(5 \[Pi] x)/3] - 128 x^6 Sinh[(\[Pi] x)/3] Sinh[(5 \[Pi] x)/3] - 1024 Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] - 1024 x^2 Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] - 384 x^4 Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] - 64 x^6 Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] - 4 x^8 Sinh[\[Pi] x] Sinh[(5 \[Pi] x)/3] + 1024 Sinh[(5 \[Pi] x)/3]^2 + 1024 x^2 Sinh[(5 \[Pi] x)/3]^2 + 384 x^4 Sinh[(5 \[Pi] x)/3]^2 + 64 x^6 Sinh[(5 \[Pi] x)/3]^2 + 4 x^8 Sinh[(5 \[Pi] x)/3]^2 - 36864 Sinh[(\[Pi] x)/3] Sinh[(7 \[Pi] x)/3] + 15360 x^2 Sinh[(\[Pi] x)/3] Sinh[(7 \[Pi] x)/3] - 1792 x^4 Sinh[(\[Pi] x)/3] Sinh[(7 \[Pi] x)/3] + 64 x^6 Sinh[(\[Pi] x)/3] Sinh[(7 \[Pi] x)/3] + 4608 Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] + 1536 x^2 Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] - 64 x^4 Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] - 32 x^6 Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] + 2 x^8 Sinh[\[Pi] x] Sinh[(7 \[Pi] x)/3] - 9216 Sinh[(5 \[Pi] x)/3] Sinh[(7 \[Pi] x)/3] - 3072 x^2 Sinh[(5 \[Pi] x)/3] Sinh[(7 \[Pi] x)/3] + 128 x^4 Sinh[(5 \[Pi] x)/3] Sinh[(7 \[Pi] x)/3] + 64 x^6 Sinh[(5 \[Pi] x)/3] Sinh[(7 \[Pi] x)/3] - 4 x^8 Sinh[(5 \[Pi] x)/3] Sinh[(7 \[Pi] x)/3] + 20736 Sinh[(7 \[Pi] x)/3]^2 - 6912 x^2 Sinh[(7 \[Pi] x)/3]^2 + 864 x^4 Sinh[(7 \[Pi] x)/3]^2 - 48 x^6 Sinh[(7 \[Pi] x)/3]^2 + x^8 Sinh[(7 \[Pi] x)/3]^2)))/(2 (65536 Cosh[(2 \[Pi] x)/ 3]^2 Sinh[\[Pi] x]^2 - 32768 x^2 Cosh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2 + 4096 x^4 Cosh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2 + 65536 x^2 Sinh[(2 \[Pi] x)/3]^2 Sinh[\[Pi] x]^2)) *) Plot[eq[[2]], {x, 3.465728, 3.465729}, PlotPoints -> 100, WorkingPrecision -> 30] Plot[Im[eq[[2]]], {x, 3.465728, 3.465729}, PlotPoints -> 200, WorkingPrecision -> 60, PlotRange -> All]

enter image description here

eq /. x -> Rationalize[3.46572839, 0] // N[#, 20] & (* Cos[y] == 0.50035241364172394089 - 9.794261144042*10^-8 I *)

Edit

f does not reach zero as it should, in the imaginary range.

Plot[f[Rationalize[3.46572839, 0], y], {y, 1030/1000, 1060/1000}, PlotRange -> 10, PlotPoints -> 200] Maximize[{f[Rationalize[3.46572839, 0], y], 1030/1000 < y < 1060/1000}, y] // N[#, 10] & (* {-7.787746817*10^-6, {y -> 1.046790571}} *) ContourPlot[ f[x, y] == 0, {x, 34657283/10000000, 6931457/2000000}, {y, 1046786/1000000, 1046795/1000000}, PlotPoints -> 200, MaxRecursion -> 5, WorkingPrecision -> 30, FrameLabel -> Automatic]

enter image description here

Improve this answer
$\endgroup$
0
$\begingroup$

Not a complete answer but one direction to go in:

cp = ContourPlot[f[x, y] == 0, {x, 3.465728, 3.465729}, {y, 1.046786, 1.046795}, PlotPoints -> 500] {l1, l2} = Cases[Normal[cp], _Line, Infinity]; {{x1, y1}} = MinimalBy[First@l1, RegionDistance[l2]]; {{x2, y2}} = MinimalBy[First@l2, RegionDistance[l1]]; Show[ cp, ListPlot[{{x1, y1}, {x2, y2},}, PlotStyle -> {PointSize[Large], Red}] ]

Output

sol1 = r /. FindRoot[f[r, y], {r, x1}]; sol2 = r /. FindRoot[f[r, y], {r, x2}]; sol1 - sol2

3.11929*10^-8

At least FindRoot also finds that there are two distinct solutions at that y. However, FindRoot also warns about its ability to find roots with the desired accuracy, so it is still not conclusive. Also, perhaps it should be looking for the root by also adjusting y.

The idea here is that we can extract values from the plot and then try to verify our conclusions about the plot using other functions.

Improve this answer
$\endgroup$
0
$\begingroup$ 
Clear["Global`*"] f[x_, y_] := 2 (-4 + x^2) Sinh[(π x)/3] + 1/16 (((4 + x^2)^2 + 64 (-4 + x^2) Cos[y] Cosh[(2 π x)/3] + 256 x Sin[y] Sinh[(2 π x)/3]) Sinh[π x] - 2 (4 + x^2)^2 Sinh[(5 π x)/ 3] + (-12 + x^2)^2 Sinh[(7 π x)/3]); ContourPlot[f[x, y] == 0, {x, 3464/1000, 3468/1000}, {y, 103/100, 106/100}, PlotPoints -> 100, MaxRecursion -> 5, WorkingPrecision -> 15, FrameLabel -> Automatic]

enter image description here

({min, arg} = NMinimize[{f[x, y]^2, 3464/1000 < x < 3468/1000, 1030/1000 < y < 1060/1000}, {x, y}, WorkingPrecision -> 30]) /. r_Real :> N[r] (* {0, {x -> 3.46572, y -> 1.04691}} *)

EDIT: Verifying the solution,

f[x, y] /. arg (* 0.*10^-20 *)

EDIT 2: Using Plot3D

Plot3D[f[x, y], {x, 3464/1000, 3468/1000}, {y, 103/100, 106/100}, MeshFunctions -> {#3 &}, Mesh -> {{0}}, MeshStyle -> Directive[Red, Thick], PlotPoints -> 100, MaxRecursion -> 5, WorkingPrecision -> 15, ClippingStyle -> None]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Sorry @BobHanlon, your /. r_Real :> N[r] makes result wrong. Do /. r_Real :> N[r, 20] or just leave it, to get only 10^-11 for f[x, y] /. arg $\endgroup$ Commented May 27, 2022 at 6:34
  • $\begingroup$ @Akku14 - If you evaluate Precision[arg] you will see that it has a precision of 30 digits. The N only affects the display (similar to a display wrapper) and is isolated from the definitions of {min, arg} by the parentheses. Subsequently, the evaluation of f[x, y] /. arg only produces 19+ digits of precision due to the computational complexity of f[x, y] $\endgroup$ Commented May 27, 2022 at 14:07
0
$\begingroup$

The following shows that f[x, y] is negative (not zero) along the vertical line x == 3.4657284... through the saddle point between the two branches, and therefore the line separates the two branches where f[x, y] == 0 in the OP's graph:

yAssum = 1.046786`32 < y < 1.046795`32; FindRoot[D[f[x, y], {{x, y}}], {{x, 3.465728}, {y, 1.04679}}, WorkingPrecision -> 32] Simplify[ Reduce[f[x, y] < 0 && yAssum /. First@%, y], yAssum] 
(* Saddle point: {x -> 3.4657284034593587205275273903929, y -> 1.0467905677870295818695381998660} *)

Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.

(* True *)

Another way to show x == 3.465... separates the branches of the curve (different, exact calculation; but same idea, negative maximum along the line):

xAssum = Rationalize[3.465 < x < 3.466, 0]; yAssum = Rationalize[1.046786 < y < 1.046795, 0]; xCP = x /. First@Solve[D[f[x, y], {{x, y}}] == 0 && yAssum && xAssum, {x, y}]; Maximize[{f[xCP, y], yAssum}, y] // N[#, 32] & (* {-0.000037312986399805657836881071055665, {y -> 1.0467905677870295818695381998661}} *)
Improve this answer
$\endgroup$