# Can't figure out how to plot this function

I've got a real nasty function that I can't figure out how to plot.

The function is

$\qquad f(k) = \pm\sqrt{1+4\cos^{2}(k/2)+4\cos(k/2)\cos(p)}$

from $k=-\pi$ to $k=\pi$, which doesn't look so bad, except $p$ is all the roots of

$\qquad \sin(pn)+2\cos(k/2)\sin(p(n+1))$

where $0<p<\pi$. There are $n$ or $n-1$ such solutions for $p$, and they depend on $k$, and I need to plot the lines for the different values of $p$ all on the same plot. I get to choose the value of $n$, which is an integer, probably around 30.

The way I've tried to do this so far is to start $k$ off at $-\pi$, slowly increment it towards $\pi$, and at each step, work out the values of $p$, and then plot those specific points for that $k$ and the corresponding $p$'s. By making the steps small enough, I get points close enough together that they look like a solid line. However, this approach is horribly slow and feels wrong. Is there a better way to go about this, preferably one where I end up with actual curves and not just a bunch of points?

Also, the function clearly has mirror symmetries, where the axes of symmetry are the axes of the plot. Is there some way to work out the plot for just the top right quadrant, and then mirror it onto the other three quadrants?

This function describes energy bands in graphene nanoribbons with zigzag edges. See Appendix B here for details.

• Please provide the equations in copy-and-pastable code form with correct Mathematica syntax. That way, potential answerers don't have to re-type in all the equations, and they can get you a solution faster. Please edit your post by clicking the grey edit button below you question and click the grey question mark on the right-side of the editing toolbar for help on how to insert code blocks. – march Nov 3 '15 at 21:08
• ...where did this "nasty function" come from? Some additional context would be nice. – J. M. is away Nov 3 '15 at 21:11
• with "p is all the roots of..." you mean the number of roots? – paw Nov 3 '15 at 21:11
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Nov 3 '15 at 21:13
• @paw. Probably the OP means the actual values of $p$ at the roots, but of course that makes the thing they're calling a function not a function, but we already knew that anyway because of the $\pm$. – march Nov 3 '15 at 21:16

Here is my take on the problem that takes the classic approach of leveraging Mathematica's plotting routines to find the roots, refine the roots using FindRoot, and then feeding these roots to the function f. Here's how it goes.

Define

num = 10;
f[k_, p_] := Sqrt[1 + 4 Cos[k/2]^2 + 4 Cos[k/2] Cos[p]]
rootFunction[k_, p_, n_] := Sin[p n] + 2 Cos[k/2] Sin[p (n + 1)]


We first use ContourPlot to very quickly find the roots of the second function.

ClearAll[pl]
pl[n_] := pl[n] = ContourPlot[rootFunction[k, p, n]
, {k, -π, π}, {p, 0, π}
, RegionFunction -> (0 < #2 < π &)
, Contours -> {0}, ContourShading -> None, PlotPoints -> 50
, Frame -> False]
pl[num] We extract the points from the plot using

roots = Cases[Normal@pl[num], Line[a_] :> a, Infinity];


and refine the roots by feeding these points to FindRoot:

refinedRoots = {#1, p /. FindRoot[rootFunction[#1, p, num], {p, #2}]} & @@@ # & /@ roots;


Note that at this point, it might spit out errors that FindRoot isn't working quite correctly. We can always check to make sure that the refinedRoots are good:

Total[Abs@rootFunction[#1, #2, 20] & @@@ #] & /@ refinedRoots Looks good. Finally, we feed these to our function

pts = {#1, f[#1, #2]} & @@@ # & /@ refinedRoots;


and plot:

ListLinePlot[pts] Since the negative of the function is just flipped over the k axis, I won't include them.

• Thank you, this is great. I've just spoken with someone about this problem, and you see that bottom band there, where it vanishes above a certain value of k? Apparently this is because p goes complex for such values of k. Is there some way to extract the complex root of this rootFunction? – Redingold Nov 4 '15 at 13:50
• @Redingold. I don't know how to fix that problem. Complex root finding is hard. If p goes imaginary, then you can replace the Sin's withSinh's and do the same thing I already did. But if they go full-on complex, the problem is trickier, and off the top of my head, I don't know how to do it. – march Nov 4 '15 at 16:23
• So, turns out that p goes to Pi + qi, and therefore, the function f goes to Sqrt[1+4Cos[k/2]^2 -4Cos[k/2]Cosh[q]], and q is given by the solution to Sinh[q n] - 2Cos[k/2]Sinh[q(n+1)] = 0, but there's a problem. The problem is that the method that you suggested for the first problem doesn't seem to work for this slightly different problem. When I implement the same logic, but with the new functions, I get errors about precision and accuracy for any n > 4. What should I do? – Redingold Nov 9 '15 at 22:09
• After some playing around, it seems that it's the values of q near k = pi that are the problem. If I simply don't search for solutions within a certain distance of k=pi (k = pi - 0.03*n seems to work as a cutoff), then I get something that looks right, because the value of f should go to 0 near k = pi. – Redingold Nov 10 '15 at 12:18
• @Redingold. It seems like you've figured it out, but a point of clarification: are you saying that in the region where the roots go complex, you make the replacement $p = \pi+iq$ and search for the roots on $q$? – march Nov 10 '15 at 17:17

If you compute the xy-points in the 1st quadrant as pairs {x, y}, then you can use the approach I take here, where I plot an entire circle from points generated in the only 1st quadrant.

Module[{q1, q2, upper, lower},
q1 = Table[{Cos[u], Sin[u]}, {u, 0., 90. °, 5. °}];
q2 = Reverse[ReflectionTransform[{1, 0}] /@ q1]; (* reflect about the y-axis *)
upper = Join[q1, q2];
lower = ReflectionTransform[{0, 1}] /@ upper; (* reflect about the x-axis *)
Graphics[{Line[upper], Line[lower]}, Frame -> True] Define

p[k_, n_] := Solve[Sin [pp n] + 2 Cos[k/2] Sin[pp (n + 1)] == 0 && 0 < pp < Pi, pp]


Then, the plot for n = 1 (for instance) is

Plot[{Sqrt[1 + 4 Cos[k/2]^2 + 4 Cos[k/2] Cos[pp /. p[k, 1]]],
-Sqrt[1 + 4 Cos[k/2]^2 + 4 Cos[k/2] Cos[pp /. p[k, 1]]]},
{k, -Pi, Pi}, PlotStyle -> Red, AxesLabel -> {k}] For n = 30, on the other hand, the following is more efficient, running in seconds.

pts = DeleteCases[MapIndexed[{Pi/60 (Last[#2] - 61), #1} &, Transpose[PadRight[
Table[Sqrt[1 + 4 Cos[k/2]^2 + 4 Cos[k/2] Cos[pp /. p[k, 30]]], {k, -Pi, Pi, Pi/60}]]],
{2}], {_, 0}, {2}];
Show[ListLinePlot[pts], ListLinePlot[-pts],
AxesLabel -> {k}, PlotRange -> All, AspectRatio -> 1] In response to the OP's final question, the plot takes advantage of symmetry about the origin. Taking advantage of symmetry about the axes also is straightforward but hardly worthwhile in view of how fast the present computation runs.

• The OP says in his post that n will be about 30. Will this approach be robust enough to handle n that large? – m_goldberg Nov 3 '15 at 22:23
• @m_goldberg Yes. I shall post it in a few minutes. – bbgodfrey Nov 3 '15 at 22:24