# Plotting Stability of damped mathieu equation

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


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


• MaxRecursion->3 ? Jul 11 '20 at 20:24

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 *)
ClearAll[mat]
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]


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},
{e1,e2} = eig[a,q] ;
If[And[Chop[Re[e1]] <= 0,Chop[Re[e2]] <= 0],1.0,I]
]
gamma = 0.0 ;
RegionPlot[Quiet[test[a,q]]\[Element]Reals,{a,-10,10},{q,-5,5}]
(* 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]


Plot[
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}}],
{t,0,10*Pi},
PlotStyle -> {Blue,Red}
]


• 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}$. Jul 12 '20 at 4:49
• 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}]  Jul 12 '20 at 4:55
• @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.
– I.M.
Jul 12 '20 at 6:46