Stability analysis for Hill's Equation

The closed cell contains a few general functions which take a formula for a differential equation (a vectorfield), convert it to an appropriate form forMathematicaand then use it to compute solutions using the built - in function NDSolve.The details of how these functions work are not so important.

MakeTFunction[x_] := x[t]
MakeTDerivs[x_] := x'[t]
MakeTFunctionSubst[x_] := x->x[t]

MakeOdeFromVf[vf_,varnames_,params_]:=Module[{varderivs,varsub,vfwithparams,i},
varderivs = Map[MakeTDerivs,varnames];
varsub= Map[MakeTFunctionSubst,varnames];
vfwithparams = vf[t,varnames,params]/.varsub;
Table[varderivs[[i]]==vfwithparams[[i]],{i,1,Length[varnames]}]]

OrbitM[vf_,varnames_,x0_,time_,params_,ndsolveopts___] :=  Module[{varfuncs,odewithparams,ic,i},
odewithparams = MakeOdeFromVf[vf,varnames,params];
varfuncs = Map[MakeTFunction,varnames];
ic = Table[varnames[[i]][0]==x0[[i]],{i,1,Length[varnames]}];
varfuncs/.NDSolve[odewithparams~Join~ic,varfuncs,{t,0,time},MaxSteps->50000][[1]]]

PoincareMap[vf_,varnames_,x0_,time_,params_,ndsolveopts___] :=  OrbitM[vf,varnames,x0,time, params,ndsolveopts]/.t->time

PoincareMapDerivative[varvf_,varnames_,x0_,time_,params_,ndsolveopts___] :=  Module[{z,point,phimatrix = {{0,0},{0,0}}},
z = PoincareMap[varvf,varnames,x0~Join~{1,0},time,params,ndsolveopts];
point = {z[[1]],z[[2]]};
phimatrix[[1,1]]=z[[3]];phimatrix[[2,1]]=z[[4]];
z = PoincareMap[varvf,varnames,x0~Join~{0,1},time,params,ndsolveopts];
phimatrix[[1,2]]=z[[3]];phimatrix[[2,2]]=z[[4]];
{point,phimatrix}]

Vectorfield for Hill’s equation and some functions for calculating the monodromy matrix and its trace for given parameter values.

Hillvf[t_,{q_,v_},{alpha_,lambda_}]={v,-(lambda + alpha*Sin[t])*q};

HillMonodromyMatrix[alpha_,lambda_]:=  Module[{z,phimatrix = {{0,0},{0,0}}},
z = PoincareMap[Hillvf,{q,v},{1,0},2*Pi,{alpha,lambda}];
phimatrix[[1,1]]=z[[1]];phimatrix[[2,1]]=z[[2]];
z = PoincareMap[Hillvf,{q,v},{0,1},2*Pi,{alpha,lambda}];
phimatrix[[1,2]]=z[[1]];phimatrix[[2,2]]=z[[2]];
phimatrix]

HillMonodromyMatrixTrace[alpha_,lambda_]:=Module[{z,tau},
z = PoincareMap[Hillvf,{q,v},{1,0},2*Pi,{alpha,lambda},WorkingPrecision->30,Method->"ImplicitRungeKutta"];
tau=z[[1]];
z = PoincareMap[Hillvf,{q,v},{0,1},2*Pi,{alpha,lambda},WorkingPrecision->30,Method->"ImplicitRungeKutta"];
tau = tau+z[[2]];
tau]

Here is the monodromy matrix for parameters alpha = 0.5 and lambda = 1.1.  Trace less than two means the origin is an elliptic fixed point.

HillMonodromyMatrix[0.5,1.1]
HillMonodromyMatrixTrace[0.5,1.1]

Here are some graphs of the trace as a function of lambda for fixed alpha = 0 and alpha = 0.5.   For the latter the trace goes slightly into the hyperbolic region near the resonances.

Plot[{2,-2,HillMonodromyMatrixTrace[0,lambda]},{lambda,-0.1,5}]

GraphicsRow[{Plot[{2,-2,HillMonodromyMatrixTrace[0.5,lambda]},{lambda,-0.1,5}],
Plot[{-2,HillMonodromyMatrixTrace[0.5,lambda]},{lambda,2.24,2.29}]}]

The Arnold tongues or resonance horns in the lambda-alpha plane.  We get a small region of instablility near each resonant forcing frequency, λ = (k/2)^2.  For larger lambda, they are too thin to show up clearly in this kind of plot.

ContourPlot[HillMonodromyMatrixTrace[alpha,lambda],{lambda,-0.2,1.2},{alpha,-0.3,0.3},Contours->{2,-2},PlotPoints->150,ContourShading->{GrayLevel[0.25],White,GrayLevel[0.75]},ContourStyle->{Thick,Black},AxesLabel->{"λ","α"},AspectRatio->1/2]