Solving Hill’s equation numerically

In this notebook we will compute fundamental matrices for Hill’s equation, a two-dimensional, nonautonomous linear systems of differential equations.

Hill’s equation

Hill’s equation is a second-order differential equations modelling a linear oscillator with time dependent forcing.  The second-oder differential equation q’’[t]+ w[t]^2 q[t] = 0 is equivalent to a following first order system with a time dependent 2x2 matrix A(t).  Here we let w(t) =  a+b*sin[t] and set up the right-hand side of the ODE.

w[t_] := a + b*Sin[t]
A[t_] := {{0,1},{-w[t]^2,0}}
f[t_] := A[t].{q[t],v[t]}

Next we set up the ODE and initial conditions as a list of equations (using a double == sign).

ode = {q'[t]==f[t][[1]], v'[t]==f[t][[2]]};
ic[{q0_,v0_}] := {q[0]==q0,v[0]==v0}
ode~Join~ic[{1,0}]

Here is how to use Mathematica’s command NDSolve to find two linearly independent solutions.  We also define a fundatmental matrix phi whose columns contain the two solutions.  The plot shows the position coordinates q1[t] and q2[t] for the two solutions.

x1={q[t],v[t]}/.NDSolve[(ode/.{a->1.3,b->0.5})~Join~ic[{1,0}],{q[t],v[t]},{t,0,16*Pi}][[1]];
x2={q[t],v[t]}/.NDSolve[(ode/.{a->1.3,b->0.5})~Join~ic[{0,1}],{q[t],v[t]},{t,0,16*Pi}][[1]];
phi:= {x1,x2}//Transpose;

Plot[{x1[[1]],x2[[1]]},{t,0,16*Pi}]

To get a plot with both positions and velocities we work in the phase plane.  First define a function to make a plot then apply it to one of the solutions.  Note that the curve crosses itself many times.  This does not violate the uniqueness theorem since the system is nonautonomous and the times of the crossings are different.

PhasePlot[{x_,y_},t0_]:= ParametricPlot[{x,y}/.t->s,{s,0,t0},
Axes->True,PlotRange->{{-2.5,2.5},{-2.5,2.5}}]

PhasePlot[x1,16*Pi]

Using the Manipulate command you can make an animation.

Manipulate[PhasePlot[x2,t0],{t0,0.001,16*Pi}]

Here is another plot showing the two columns of the fundamental matrix as vectors.

VWPlot[v_,w_,t0_]:= Show[Graphics[{Arrow[{{0,0},v/.t->t0}],Arrow[{{0,0},w/.t->t0}]}],
Axes->True,PlotRange->{{-2,2},{-2,2}}]

Manipulate[VWPlot[x1,x2,t0],{t0,0,16*Pi}]

In spite of the fact that the matrix changes in a complicated way, its determinant, the Wronsian, is constant.

Plot[Det[phi],{t,0,16*Pi}]

Changing the values of a and b can lead to a different kind of behavior.  Here the solutions grow and decay, but still the determinant is constant.

x1={q[t],v[t]}/.NDSolve[(ode/.{a->1.0,b->0.5})~Join~ic[{1,0}],{q[t],v[t]},{t,0,16*Pi}][[1]];
x2={q[t],v[t]}/.NDSolve[(ode/.{a->1.0,b->0.5})~Join~ic[{0,1}],{q[t],v[t]},{t,0,16*Pi}][[1]];
phi:= {x1,x2}//Transpose;

Plot[{x1[[1]],x2[[1]]},{t,0,16*Pi}]
Manipulate[VWPlot[x1,x2,t0],{t0,0,16*Pi}]
Manipulate[PhasePlot[x2,t0],{t0,0.001,16*Pi}]
Plot[Det[phi],{t,0,16*Pi}]

The Poincare Map

Hill’s equation illustrates another aspect of nonautonomous linear systems.  The forcing function w(t) and the corresponding matrix A(t) are periodic functions of time with period 2*Pi.  In this case the time advance 
matrices Psi( t, t0) are also periodic:  Psi( t + 2*Pi , t0 + 2*Pi) = Psi( t, t_0).  In particular the maps which advance time by 2*Pi satisfy:
Psi( 2*Pi, 0 ) = Psi(4*Pi, 2*pi) = Psi(6*Pi, 4*Pi) =  ....
Call this the Poincare map F = Psi( 2*Pi , 0 ).  We can advance time by multiples of 2*Pi by iterating this linear map, that is, composing the linear map with itself or equalently, taking powers of the matrix of F.

Here are the Poincare maps for the two choices of parameters studied above.

x1={q[t],v[t]}/.NDSolve[(ode/.{a->1.3,b->0.5})~Join~ic[{1,0}],{q[t],v[t]},{t,0,2*Pi}][[1]];
x2={q[t],v[t]}/.NDSolve[(ode/.{a->1.3,b->0.5})~Join~ic[{0,1}],{q[t],v[t]},{t,0,2*Pi}][[1]];

phi:= {x1,x2}//Transpose;
F1= phi/.t->2*Pi;
F1//MatrixForm

x1={q[t],v[t]}/.NDSolve[(ode/.{a->1.,b->0.5})~Join~ic[{1,0}],{q[t],v[t]},{t,0,2*Pi}][[1]];
x2={q[t],v[t]}/.NDSolve[(ode/.{a->1.,b->0.5})~Join~ic[{0,1}],{q[t],v[t]},{t,0,2*Pi}][[1]];

phi:= {x1,x2}//Transpose;
F2 =phi/.t->2*Pi;
F2//MatrixForm

The eigenvalues of the Poincare map determine the long term behavior of the solutions.  In the first case, the eigenvalues are complex numbers on the unit circle and the corresponding solutions remain bounded.  In the second case, they are real numbers lambda_1>1 and lambda_2= 1/lambda_1 and we get exponential growth.

Eigenvalues[F1]
Eigenvalues[F2]