I'm working with the Mathieu equation as a part of my research and as part of the analysis I'd like to include a stability chart such as Fig. 1 in this paper. Unfortunately I'm stuck with Mathematica and the only other questions on here regarding this topic are extremely specific and don't really help me. Despite Mathematica having Mathieu functions built in, plotting these does not result in anything close to what I'm looking for. Any help would be appreciated!

  • 2
    $\begingroup$ Hi @jj7510! Can you show what you've tried and where you got stuck? $\endgroup$
    – Chris K
    Commented Jan 15, 2019 at 14:36
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) 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. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Chris K
    Commented Jan 15, 2019 at 14:37
  • 1
    $\begingroup$ What is the formula being plotted in Fig. 1? $\endgroup$
    – Michael E2
    Commented Jan 15, 2019 at 14:56
  • $\begingroup$ Granted, I called it a diagram, but this question may be a duplicate. $\endgroup$ Commented Jan 16, 2019 at 2:32
  • $\begingroup$ Hey @ChrisK, it's less getting stuck and more not knowing where to begin. Mathematica has solution functions for the Mathieu equation but I can't figure out how one might go about getting a diagram in the Mathieu parameters from these. $\endgroup$
    – jj7510
    Commented Jan 18, 2019 at 13:11

1 Answer 1


It looks to me like the authors are plotting where the equations are stable (periodic) and where they are unstable (grow exponentially). The red in Fig. 1 indicates where the solutions are unstable with $\epsilon$ vs. $\delta$. In Mathematica's notation for MathieuC and MathieuS, those would be the terms q and a, respectively. I'm not sure I understand the difference between the two Mathieu functions, but if you plot them for a and q in the white region, you get something periodic:

a = 7.5 (* delta *)
q = 2.5 (* epsilon *)
Plot[{Re[MathieuC[a, q, z], MathieuS[a, q, z]},
    {z, 0, 20Pi},

MathieuC and S plots.

You can see that these plots seem to oscillate fairly nicely. In order to take a stab at replicating the plot of Fig. 1, I basically used a table to calculate the value of MathieuC at z = 1000 for a given $a$ and $q$. If the value was greater than 1000 at that z-point (edit: I realized it might be confusing because I chose 1000 to be both the z-value on the x-axis at which the function is evaluated as well as the cutoff value on the y-axis, but they are different), I assumed it was increasing exponentially and set it to "1". If it was less than that, I assumed it was oscillating nicely and set the value to "0".

table = ParallelTable[
   {a, q, If[Re[MathieuC[a, q, 1000]] >= 1000, 1, 0]}, {a, 0, 15, 
    0.1}, {q, 0, 15, 0.1}];
ListDensityPlot[Flatten[table, 1], ImageSize -> 400,
 ColorFunction -> (If[# < 0.5, White, Red] &)]

My attempt at recreating Fig. 1.

I realize that the code may not be the most elegant, but I'm hoping it will work as a starting point for your own explorations, or else a jumping-off point for someone who knows more about these sorts of functions.

Edit 2: I'm including a plot with more points just for the fun of it. For this plot I had $a$ ranging from -2.5 to 16 and $q$ ranging from 0 to 18 to more closely match the figure in the paper. I decreased the step size to 0.025 for both, and I now take the Abs[] of the value as I found that some parameters yielded large negative values (~ $-10^{562}$).

amin = -2.5;
amax = 16;
qmin = 0;
qmax = 18,
stepsize = 0.025;
table2 = ParallelTable[
   {a, q, If[Abs[Re[MathieuC[a, q, 1000.]]] >= 1000, 1, 0]}, {a, amin,
     amax, stepsize}, {q, qmin, qmax, stepsize}];

Improved Mathieu Density Plot.

  • $\begingroup$ Thank you so much! Btw, MathieuC and MathieuS are the even/odd solutions to the Mathieu eqn, similar to the role filled by $\cos x$ and $\sin x$ for the equation $x''+x=0$. Everything you said until after the first plot is what I've been able to figure out, but it was the stuff you did after that serves as the perfect starting point for me and answers my question. Thanks again! $\endgroup$
    – jj7510
    Commented Jan 18, 2019 at 13:20
  • $\begingroup$ @jj7510 If this answers your question, you can accept it by clicking on the checkmark sign. $\endgroup$
    – Chris K
    Commented Jan 18, 2019 at 13:55
  • $\begingroup$ @jj7510 No problem, glad I could help! Ah, thanks for the explanation! $\endgroup$
    – MassDefect
    Commented Jan 18, 2019 at 22:46
  • $\begingroup$ I am having a problem while plotting a stability chart. I am using your code. I want to see the intersection of the stable area with the ordinate axis. For this, I choose the interval q from 40 to 42, for example. As a result, the chart is missing and an empty chart window is displayed. Could you suggest what needs to be changed in the program to display the graph? $\endgroup$ Commented Jul 28, 2020 at 17:58
  • $\begingroup$ @AndrewMarphichev Looking back at it, a density plot probably isn't the best method to plot this due to interpolation. I would recommend switching to ArrayPlot or MatrixPlot and using a very small step size in order to see the white points. You could do something like ArrayPlot[ Reverse[table2[[All, All, 3]]\[Transpose]], DataRange -> {{amin, amax}, {qmin, qmax}}, FrameTicks -> All ] $\endgroup$
    – MassDefect
    Commented Jul 28, 2020 at 19:09

Your Answer

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

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