I am trying to create the stability diagram of the damped Mathieu equation using Mathematica.The Mathieu equation is $$D(y)+(a-2q \cos(2t))y=0$$ where D(y) is the second-order derivative with respect to t. The damped Mathieu equation is $$D(y)+(a-f-2q \cos(2t))y=0$$ where f is the damping constant. The real values of the Mathieu characteristic exponent will give a stable region. Damping increases the stability region. But I am not getting that. The Mathematica codes are the following

RegionPlot[MathieuCharacteristicExponent[a , q]\[Element]Reals,{a,-10,10},{q,-5,5}]

enter image description here

RegionPlot[MathieuCharacteristicExponent[(a-2) , q]\[Element]Reals,{a,-10,10},{q,-5,5}]

enter image description here

  • 1
    $\begingroup$ MaxRecursion->3 ? $\endgroup$
    – yarchik
    Jul 11, 2020 at 20:24

1 Answer 1


Your 'dumping' constant is a shift of $a$ parameter as it can be seen from your pics. To some get actual dumping your DE should be of the form:

$$y'' + g y' + (a-2 q \cos(2 t)) y = 0$$

This equation is linear with periodic coefficients, stability info can be obtained from eigenvalues of the monodromy matrix.

(* t-map *)
map = DSolveValue[{y''[t] +g y'[t] +  (a - 2*q*Cos[2 *t])  y[t] == 0,y[0]==X,y'[0]==Y},{y[t],y'[t]},t] ;
(* 1-map *)
one = map /. t -> Pi ;
(* monodromy matrix *)
mat[g_][a_,q_] = D[one,{{X,Y}}] ;
(* eigen values *)
ClearAll[eig] ;
eig[a_,q_] :=Eigenvalues[(1/2 Pi MatrixLog[mat[gamma][a,q]])] ;

For $g=0$, something similar to a characteristic exponent can be used to determine the stability region in parameter space:

gamma = 0.0 ;
RegionPlot[Quiet[Chop[First[eig[a , q]]/I]]\[Element]Reals,{a,-10,10},{q,-5,5},ImageSize -> Medium]

enter image description here

Another option is to define stable region as a region where real parts of both eigenvalues are less or equal to zero:

ClearAll[test] ;
test[a_,q_] := Block[
{e1,e2} = eig[a,q] ;
 If[And[Chop[Re[e1]] <= 0,Chop[Re[e2]] <= 0],1.0,I]
gamma = 0.0 ;
(* same output *)

For $g < 0$ all points are unstable in the sence $y(t\to\infty)\to\infty$, for some $g > 0$ the region with both negative real parts:

gamma = 0.5 ;
plot = RegionPlot[Quiet[test[a,q]]\[Element]Reals,{a,-10,10},{q,-5,5}] ;
p1 = {4.5,2.0} ;
p2 = {4.5,2.5} ;
Show[plot,Graphics[{PointSize[Large],Blue,Point[p1],Red,Point[p2]}],ImageSize -> Medium]

enter image description here

Evaluate[{First[map]  /. {X -> 1,Y->0, g -> 0.5, a -> 4.5,q -> 2.0},First[map]  /. {X -> 1,Y->0, g -> 0.5, a -> 4.5,q -> 2.5}}],
PlotStyle -> {Blue,Red}

enter image description here

  • 1
    $\begingroup$ The damped Mathieu equation of $y$ can be transformed into Mathieu equation of $x$ if one use the transformation $y=xe^{-g t}$. For that $a$ will become $a-\frac{g^2}{4}$. $\endgroup$
    – Sarthok
    Jul 12, 2020 at 4:49
  • $\begingroup$ I want to understand the stability of damped Mathieu equation from the Floquet theory. The real values of the Floquet exponent or Mathieu Characteristic exponent describe the stability zones. If damping increases the stability, then why the stability zones are decreased with increasing damping in the code RegionPlot[MathieuCharacteristicExponent[a-2 , q]\[Element]Reals,{a,-10,10},{q,-5,5}] $\endgroup$
    – Sarthok
    Jul 12, 2020 at 4:55
  • 1
    $\begingroup$ @sarthok, there is no dumping in MathieuCharacteristicExponent function. If you want analytical results, study powers of monodromy matrix mat, no need to compute eigenvalues, without dumping you'll get a constraint for matrix trace. $\endgroup$
    – I.M.
    Jul 12, 2020 at 6:46

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.