# Issue with ContourPlot

I am trying to Contour plot the solutions "$$x$$" of an equation versus $$k$$ for some fixed $$T$$ values. However, the contour plot is giving incorrect plot perhaps due to some accuracy issues. e.g., at small $$T$$ ($$T\leq 1$$) I should get three solutions, but I am not getting three on the contour plot.? The code:

eq[x_, k_, T_] := -Sin[3*k + x]/Sin[2*k + x] + z + 2*Cos[k] + T^2 + (A*T^2*Sin[k]^2)/(Sin[2*k + x]^2 + B*T^4*Sin[k]^2) ==  0 /. {A -> 1/2, B -> 0.00001, z -> -2.37}
ContourPlot[Evaluate[eq[x, k, 0.1]], {x, 0, Pi}, {k, 0, Pi}]

• You say three, but I get six: NSolve[{eq[x, k, 1/10] /. x -> 1.5, 0 < k < Pi}, k, WorkingPrecision -> 16] – Michael E2 Jul 25 '19 at 0:21
• Close, but distinct: wp = 100; k0 = k /. NSolve[{eq[x, k, 1/10] == 0 /. x -> x0, 0 < k < Pi}, k, WorkingPrecision -> wp]; ListLinePlot[ Table[ eq[1, k, 1/10], {k, k0[[2]] - 1*^-7, k0[[2]] + 1*^-7, 1*^-9}], GridLines -> {k0, None}, DataRange -> {k0[[2]] - 1*^-7, k0[[2]] + 1*^-7}] --- ListLinePlot[ Table[ eq[1, k, 1/10], {k, k0[[3]] - 1*^-7, k0[[3]] + 1*^-7, 1*^-9}], GridLines -> {k0, None}, DataRange -> {k0[[3]] - 1*^-7, k0[[3]] + 1*^-7}] – Michael E2 Jul 25 '19 at 1:09
• @MichaelE2 Thanks for your elaborations. Actually x are the sought solutions (a phase to be found between $0\leq x\leq \pi$). You fixed it to x=1.5, When I replace xwith k I get the three solutions for each fixed k: cf. NSolve[{eq[x, k, 1/10] /. k -> 1.5, 0 < x < Pi}, x, WorkingPrecision -> 16] – AtoZ Jul 25 '19 at 7:00
• Oops, my bad. Excuse me. I was paying more attention to the ContourPlot, and thinking you were solving for k as a function of x (and I was getting six lines), than to what you wrote. I get only three roots $x$ for $k(x) = k_0$, as well. The problem with ContourPlot is that two of solutions are too close to be resolved by the typical sampling density of ContourPlot, which finds curves by sign changes. The plot should still look like the plot in my answer, but, as you can see, two solutions overlap each other at normal image resolution. – Michael E2 Jul 25 '19 at 15:15
• I updated my answer. I'm not sure entirely what you would like. I think I explained why ContourPlot failed. Theoretically, ContourPlot should start to resolve the roots with PlotPoints -> 3000 ($\Delta x \approx 0.0005-0.005$, interval length $\pi \approx 3$), maybe considerably less, but that's a prohibitive computation for my laptop. And you would need to make the image roughly 3000 pixels wide, or zoom in. Please look at my answer, and if there's something you'd like me to show you, I'll try. – Michael E2 Jul 25 '19 at 16:04

It appears that the solution of x -> 2.94159 that you mention in a comment in @Hugh 's answer is incorrect (or maybe better said as "inappropriate"). (@Hugh 's answer essentially says it all.)

That particular solution is really at a discontinuity. If you plug in $$\pi-2/10$$ for $$x$$ (i.e., essentially what x -> 2.94159 means), you'll get

Consider a slight change in your equation (rationalizing the constants and removing the ==0):

eq[x_, k_, T_] := -Sin[3*k + x]/Sin[2*k + x] + z + 2*Cos[k] +
T^2 + (A*T^2*Sin[k]^2)/(Sin[2*k + x]^2 + B*T^4*Sin[k]^2) /.
{A -> 1/2, B -> 10^-5, z -> -237/100}


Then we have

eq[x, 1/10, 1]


If $$\pi-2/10$$ is substituted for $$x$$, then $$\csc \left(x+\frac{1}{5}\right)$$ results in ComplexInfinity:

Csc[1/5 + x] /. x -> π - 2/10


So, in short, there are only two contour lines of zero, not three, when $$T=0.1$$.

Oooops! I guess there are more than two contours. See @MichaelE2 's answer.

In general we have for eq[x,k,T]

The discontinuity occurs when Csc[2 k + x] is ComplexInfinity or whenever $$2k+x=\pi$$ or $$2k+x=2\pi$$ no matter what the value of $$T$$ happens to be for the ranges of interest of $$x$$ and $$k$$. So you could always include those (dotted) lines of discontinuity to your contour plots:

Here is the code to produce the animated contour plot:

t = {"0.025", "0.050", "0.075", "0.100", "0.125", "0.150", "0.175", "0.200", "0.225",
"0.250", "0.275", "0.300", "0.325", "0.350", "0.375", "0.400", "0.425", "0.450",
"0.475", "0.500", "0.525", "0.550", "0.575", "0.600", "0.625", "0.650", "0.675",
"0.700", "0.725", "0.750", "0.775", "0.800", "0.825", "0.850", "0.875", "0.900",
"0.925", "0.950", "0.975", "1.000"};
eq[x_, k_, T_] := -Sin[3*k + x]/Sin[2*k + x] + z + 2*Cos[k] +
T^2 + (A*T^2*Sin[k]^2)/(Sin[2*k + x]^2 + B*T^4*Sin[k]^2) /. {A -> 1/2, B -> 10^-5, z -> -237/100}
g = Table[Show[ContourPlot[Evaluate[eq[x, k, T/40] == 0], {x, 0, π}, {k, 0, π},
PlotPoints -> 100, AspectRatio -> 1,
FrameLabel -> (Style[#, Bold, 18] &) /@ {"x", "k"},
PlotLabel -> Style["T = " <> t[[T]], Bold, 24]],
Plot[{π/2 - x/2, π - x/2}, {x, 0, π},
PlotRange -> {{0, π}, {0, π}}, AspectRatio -> 1,
PlotStyle -> {{Gray, Dotted}}]], {T, 40}];
Export["contours.gif", Flatten[{g, Table[g[[i]], {i, Length[g], 2, -1}]}],
"DisplayDurations" -> 0.25, AnimationRepetitions -> Infinity]

• Thanks. How did you create the above plot, with the dotted lines of discontinuity? May I see? – AtoZ Jul 25 '19 at 8:12
• Thank you very much ! – AtoZ Jul 26 '19 at 20:25

Not an answer but if you do

Plot3D[Evaluate[eq[x, k, 0.1]], {x, 0, Pi}, {k, 0, Pi},
PlotRange -> All]


Division by the sin function has zeros which dominate the plot. What are you hoping for?

• Thanks. When employing the numerical solution by sol[T_?NumericQ, k_ /; 0 <= k <= Pi] := x /. NSolve[{eq[x, k, T], 0 <= x <= Pi}, x, Reals] sol[1, 0.1] for example gives solution for $T=1$, $k=0.1$, which is {2.64172, 2.94159, 2.98454} i.e., three solutions. But on the contour plot, I see only two? – AtoZ Jul 24 '19 at 17:20

Update: I misread the OP and counted the six distinct function for $$k(x)$$ as solutions. The OP was trying to solve equations of the form $$k(x) = k_0$$, for which I get only three solutions for $$0 < k_0 < \pi$$, just as expected. ContourPlot has a rather weak ability to resolve roots/curves that are close together. It detects solutions via sign changes, and to detect a pair of roots it has to sample the function between them; otherwise, there would be no sign change and ContourPlot will assume there are no solutions there.

I proposed the NDSolve method below as an alternative to ContourPlot. One can use the power of NSolve to resolve the roots, and NDSolve can trace roots throughout the domain. If preferred, one could switch the dependent and independent variables and solve for x[k] as a function of k. There is a discontinuity at k = Pi/2 and one would have to solve twice, once for each of the domains 0 < k < Pi/2 and Pi/2 < k < Pi (or one could integrate over {x, 0, 2 Pi} and use Mod[x, Pi] to construct the solutions over the lower domain.

For instance, one can zoom in with PlotRange and show the three roots around $$k_0 = 3.1355$$:

I find six distinct solutions functions $$k(x)$$ for T = 1/10 at high working precision:

eq[x_, k_, T_] := -Sin[3*k + x]/Sin[2*k + x] + z + 2*Cos[k] +
T^2 + (A*T^2*Sin[k]^2)/(Sin[2*k + x]^2 + B*T^4*Sin[k]^2) /.
{A -> 1/2, B -> 1/100000, z -> -237/100};
x0 = 1;
wp = 100;   (* working precision *)
k0 = k /. NSolve[{eq[x, k, 1/10] == 0 /. x -> x0, 0 < k < Pi},  (* six roots *)
k, WorkingPrecision -> wp];
Quiet[
ndsol = First@NDSolve[{
D[eq[x, k[x], 1/10], x] == 0,
k[x0] == #},
k, {x, 0, Pi},
PrecisionGoal -> 20, WorkingPrecision -> wp] & /@ k0,
{Power::infy, Infinity::indet}];
k["Domain"] /. ndsol // N
(*  messages Power::infy, Infinity::indet, NDSolve::ndsz,... omitted  *)
(*  Domains of solutions:
{{{0., 3.14159}}, {{0., 3.14159}}, {{0., 3.14159}},
{{8.10181*10^-99, 3.14159}}, {{4.85285*10^-55, 3.14159}},
{{7.61263*10^-99, 3.14159}}}
*)

ListLinePlot[k /. ndsol, PlotLegends -> Thread[k[x0] == N@k0], ImageSize -> Large]


Two pair are close to each other but distinct. For instance:

ListLinePlot[
Table[
eq[1, k, 1/10], {k, k0[[2]] - 1*^-7, k0[[2]] + 1*^-7, 1*^-9}],
GridLines -> {k0, None},
DataRange -> {k0[[2]] - 1*^-7, k0[[2]] + 1*^-7}]

ListLinePlot[
Table[
eq[1, k, 1/10], {k, k0[[3]] - 1*^-7, k0[[3]] + 1*^-7, 1*^-9}],
GridLines -> {k0, None},
DataRange -> {k0[[3]] - 1*^-7, k0[[3]] + 1*^-7}]


• My last numerical analysis class was 45 years ago so please bear with me. I'm thinking that a zero doesn't occur at k = k0[[2]]. k0[[2]] is essentially $\frac{1}{2} (\pi -1)$ (a point of discontinuity). This is what the following plot tells me: Show[Plot[eq[1, k, 1/10], {k, 1, 1.2}, PlotPoints -> 1000, PlotRange -> {{1.06, 1.08}, {-100, 1000}}, MaxRecursion -> 5, WorkingPrecision -> 100], ListPlot[{{{k0[[2]], 0}, {k0[[2]], 1000}}, {{k0[[3]], 0}, {k0[[3]], 1000}}}, Joined -> True, PlotStyle -> {Red, Green}, PlotLegends -> {"k0[[2]]", "k0[[3]]"}]]. – JimB Jul 25 '19 at 3:36
• @JimB The picture is not perfectly clear (to me), and the appearance of the roots k0[[2]], k0[[3]] may even both be spurious. But the roots are definitely there for the exact parameters I used. They are also stable under perturbations of the parameters {A -> 1/2 (1 + eA), B -> 1/100000 (1 + eB), z -> -237/100 (1 + ez), T -> 1/10 (1 + eT)} for eA, eB, ez, eT up to ±10^-3 (k0[[2]] is about $8 \times 10^{-8}$ above $(\pi-1)/2$). I think what is happening is that as T -> 0, both x0[[2]] and x0[[3]] approach $(\pi-1)/2$. But I suspect they are distinct roots for T > 0. – Michael E2 Jul 25 '19 at 4:41
• Thanks. Mathematics is not my strong point. I'll look at it and what you wrote in the comment. – JimB Jul 25 '19 at 4:44