6
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ 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. $\endgroup$
    – march
    Nov 3, 2015 at 21:08
  • $\begingroup$ ...where did this "nasty function" come from? Some additional context would be nice. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Nov 3, 2015 at 21:11
  • $\begingroup$ with "p is all the roots of..." you mean the number of roots? $\endgroup$
    – paw
    Nov 3, 2015 at 21:11
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Nov 3, 2015 at 21:13
  • $\begingroup$ @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$. $\endgroup$
    – march
    Nov 3, 2015 at 21:16

3 Answers 3

Reset to default
6
$\begingroup$

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]

enter image description here

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

enter image description here

Looks good. Finally, we feed these to our function

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

and plot:

ListLinePlot[pts]

enter image description here

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

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – Redingold
    Nov 4, 2015 at 13:50
  • $\begingroup$ @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. $\endgroup$
    – march
    Nov 4, 2015 at 16:23
  • $\begingroup$ 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? $\endgroup$
    – Redingold
    Nov 9, 2015 at 22:09
  • $\begingroup$ 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. $\endgroup$
    – Redingold
    Nov 10, 2015 at 12:18
  • $\begingroup$ @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$? $\endgroup$
    – march
    Nov 10, 2015 at 17:17
4
$\begingroup$

This answers your 2nd question -- about using symmetry to reduce the amount of computation needed.

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]

plot

Improve this answer
$\endgroup$
4
$\begingroup$

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

enter image description here

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]

enter image description here

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ The OP says in his post that n will be about 30. Will this approach be robust enough to handle n that large? $\endgroup$
    – m_goldberg
    Nov 3, 2015 at 22:23
  • $\begingroup$ @m_goldberg Yes. I shall post it in a few minutes. $\endgroup$
    – bbgodfrey
    Nov 3, 2015 at 22:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.