Preface

This section provides an introduction to the Floquet--Lyapunov theory.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part I of the course APMA0340
Introduction to Linear Algebra with Mathematica

Glossary

A classical context of the Floquet theory was discovered in 1883 by a French mathematician Gaston Floquet (1847--1920). It was intended to describe the behavior of a set of linear differential equations with a time-periodic coefficients, which was in turn originating from the problem of the stability of periodic orbits in classical mechanics. In the quantum world, where the linearity of the Schrödinger equation is guaranteed from the start, the Floquet theory applies whenever the Hamiltonian governing the system is time-periodic.

 

One dimensional case

The most general homogeneous first order differential equation has the form:
\begin{equation} \label{EqOne.1} \frac{{\text d}x(t)}{{\text d}t} = a(t) \,x(t) , \qquad a(t+T) = a(t) , \end{equation}
where 𝑎(t) is a T-periodic function of time. This separable equation can be solved by direct integration
\begin{equation} \label{EqOne.2} x(t) = x(t_0 )\,\exp \left\{ \int_{t_0}^t a(\tau )\,{\text d}\tau \right\} . \end{equation}

If x(t) is a solution of Eq.\eqref{EqOne.1}, then the function y(t) = x(t+T) is also a solution because

\[ \frac{{\text d}y}{{\text d}t} = \frac{{\text d}}{{\text d}t} \, x(t+T) = \frac{{\text d}x(t+T)}{{\text d}(t+T)} = a(t+T)\, x(t+T) . \]
Hence, x(t) and y(t) satisfy the same equation. They are, however, not necessarily the same function, but because the differential equation \eqref{EqOne.1} is of first order and linear, there is only one linearly independent solution. Therefore, we must have
\[ x(t+T) = c\,x(t) , \]
where c is a constant independent of time. Then for any integer n ∈ ℤ, we have
\[ x(t+ nT) = c^n x(t) . \]
This equation has a particularly simple interpretation. If c > 1, then cn grows exponentially with increasing n, so solution x(t) also increases exponentially. However, if c < 1, then cn → 0, and x(t) → 0 as t → +∞. In the special case c = 1, x(t) = x(t+T), for all t, and the solution is periodic. This occurs only when the mean value of 𝑎(t) is zero. So from the solution formula \eqref{EqOne.2}, we get
\[ c = \exp\left\{ \int_{0}^T a(\tau )\,{\text d}\tau \right\} . \]
Example 1: Let us consider a periodic function
\[ a(t) = a+ b\,\sin \omega t , \]
where 𝑎, b, and ω are real constants. In this case, Eq.\eqref{EqOne.2} becomes
\[ x(t) = x(0) \,\exp\left\{ \int_0^t \left( a + b\,\sin \omega\tau \right) {\text d}\tau \right\} = x(0) \,\exp\left\{ at - \frac{b}{\omega}\,\cos \cos \omega t + \frac{b}{\omega} \right\} . \tag{1.1} \]
If 𝑎 < 0, then c < 1 and x(t) → 0 as t → +∞. If 𝑎 > 0, then c > 1 and the solution x(t) grows without bound being 2π/ω-periodic. When 𝑎 = 0, then the solution
\[ x(t) = x(0) \,e^{b/\omega} \, e^{-(b/\omega )\,\cos \omega t} \]
is 2π/ω-periodic.

      We plot solutions for some values of parameter 𝑎, but fixing b = ω = 1.
Plot[{Callout[Exp[t/3 - Cos[t]], "a=1/3", Above, LabelStyle -> Directive[Italic, Medium, ColorData[100, 1]]], Callout[Exp[-Cos[t]], "a=0"], Callout[Exp[-3*t/2 - Cos[t]], "a=-3/2"]}, {t, 0, 10}, PlotRange -> {0, 3.6}, PlotTheme -> "Web"]
       Three solutions of one-dimensional Floquet equation.            Mathematica code

    ◾
Example 2: Let us consider the initial value problem
\[ \frac{{\text d}x}{{\text d}t} = x(t)\,\cos t\,\sin \omega t , \qquad x(0) = 1, \qquad 0 < \omega \le 1, \tag{2.1} \]
with ω = p/q, where p and q are prime positive integers.

We expect that its solution to be periodic with period T = πq. According to Eq.\eqref{EqOne.2}, its solution is

\[ x(t) = \exp \left\{ \int_0^t \cos\tau\,\sin \omega\tau\,{\text d}\tau \right\} = \exp \left\{ \frac{\omega}{\omega^2 -1} \left( 1 - \cos t\,\cos \omega t \right) - \frac{1}{\omega^2 -1} \,\sin t\,\sin \omega t \right\} . \]
Integrate[Cos[s]*Sin[a*s], {s, 0, t}]
(a - a Cos[t] Cos[a t] - Sin[t] Sin[a t])/(-1 + a^2)
Since the limit
\[ \lim_{\omega \to 1} \left[ \frac{\omega}{\omega^2 -1} \left( 1 - \cos t\,\cos \omega t \right) - \frac{1}{\omega^2 -1} \,\sin t\,\sin \omega t \right] = \frac{1}{2}\,\sin^2 t , \]
f[a_, x_] = (a - a Cos[t] Cos[a t] - Sin[t] Sin[a t])/(-1 + a^2);
Limit[f[a, t], a -> 1]
Sin[t]^2/2
is finite for any t, the function x(t) is bounded for all t ∈ ℝ.
      We plot solutions for some values of parameter ω.
f[a_, x_] = (a - a Cos[t] Cos[a t] - Sin[t] Sin[a t])/(-1 + a^2);
Plot[{Callout[f[1/3, t], "a=1/3", Above, LabelStyle -> Directive[Italic, Medium, ColorData[100, 1]]], Callout[f[1/2, t], "a=1/2"], Callout[f[0.999, t], "a=0.999"]}, {t, 0, 4}, PlotRange -> {0, 0.6}, PlotTheme -> "Web"]
       Three solutions of one-dimensional Floquet equation.            Mathematica code

    ◾

 

Multi-dimensional linear periodic systems

The general theory of n-dimensional linear, periodic systems is very similar to the theory of one-dimensional linear periodic equations. We consider the linear system of differential equations with periodic coefficients

\begin{equation} \label{EqFloquet.1} \frac{{\text d}{\bf x}(t)}{{\text d}t} = {\bf A}(t) \,{\bf x}(t) , \qquad {\bf A}(t+T) = {\bf A}(t) , \end{equation}
where A(t) is a real, non-singular n×n matrix with elements that are T-periodic functions of t.

This linear system of equation \eqref{EqFloquet.1} has a fundamental matrix that we denote by X(t). This means that detX(t) ≠ 0 and it satisfies the matrix differential equation

\[ \frac{{\text d}}{{\text d}t}\, {\bf X}(t) = {\bf A}(t) \,{\bf X}(t) , \qquad {\bf A}(t+T) = {\bf A}(t) . \]
So every column xk(t), k = 1, 2, …, n, of matrix X(t) is a solution of the first order system of equation \eqref{EqFloquet.1}. The vectors yk(t) = xk(t+T), k = 1, 2, …, n, are also solutions of the original vector equation \eqref{EqFloquet.1} because
\[ \frac{{\text d}{\bf y}_k (t)}{{\text d}t} = \frac{{\text d}}{{\text d}t}\, {\bf x}_k (t+T) = \frac{{\text d}{\bf x}_k (t+T)}{{\text d}(t+T)} = {\bf A}(t+T) \,{\bf x}_k (t+T) = {\bf A}(t) \,{\bf y}_k (t) . \]
Hence, every column-vector yk(t) = xk(t+T), k = 1, 2, …, n, must be a linear combination of the columns of the fundamental matrix X(t):
\[ {\bf y}_k (t) = \sum_{1 \le j \le n} e_{j,k} {\bf x}_j (t) \]
for some constants ej,k. These new solutions may be used to form another fundamental matrix (because the fundamental matrix has non zero determinant for all t ∈ ℝ).
\begin{equation} \label{EqFloquet.2} {\bf X}(t+T) = {\bf X}(t) \,{\bf E} \qquad \Longleftrightarrow \qquad {\bf E} = {\bf X}^{-1}(t) \,{\bf X}(t+T) , \end{equation}
where E is the constant matrix with elements ei,j. Since X(t) is a fundamental matrix for Eq.\eqref{EqFloquet.1}, matrix E is invertible; it is usually referred to as the monodromy matrix.

For any fundamental matrix Φ(t) of Eq.\eqref{EqFloquet.1}, the nonsingular matrix
\[ {\bf E} = {\bf \Phi}(t+T)\,{\bf \Phi}^{-1}(t) \]
is called a monodromy operator or matrix. The eigenvalues of a monodromy operator are called characteristic multipliers of the corresponding time-periodic homogeneous system \eqref{EqFloquet.1}. Moreover, all monodromy operators have the same eigenvalues. In particular, there are exactly n characteristic multipliers, counting multiplicities. A complex number μ is called a characteristic exponent (or a Floquet exponent) of the system, if ρ is a characteristic multiplier and eμT = ρ.

If λ is an eigenvalue of E and v is the corresponding eigenvector, E v = λv, then for every fundamental matrix X(t), the solution

\[ {\bf u}(t) = {\bf X}(t)\,{\bf v} \qquad\mbox{has the property }\quad {\bf u}(t+T) = \lambda {\bf u}(t) \qquad\mbox{for all }t. \]
It follows from
\[ {\bf u}(t+T) = {\bf X}(t+T){\bf v} = {\bf X}(t){\bf E}\,{\bf v} = \lambda {\bf u}(t) \qquad\mbox{for all }t. \]
The eigenvalues of the monodromy matrix E are therefore named characteristic multipliers or numbers of the system \eqref{EqFloquet.1}. They are eigenvalues of eTK. The eigenvalues of E are independent of the choice of fundamental matrix so are a property of the system, not a particular fundamental matrix. Indeed, if two different fundamental matrices X(t) and Y(t) give rise to two distinct monodromy matrices E1 and E2, then X(t) = CY(t), for some nonsingular matrix C. We have
\[ {\bf \Phi}_{1}(t+T) = {\bf \Phi}_{2}(t+T){\bf C} = {\bf \Phi}_{2}(t) {\bf E}_2 {\bf C} \qquad \Longrightarrow \qquad {\bf E}_2 = {\bf C}^{-1} {\bf E}_1 {\bf C} . \]
Hence, E1 and E2 are similar matrices, so they have the same eigenvalues.

It is convenient to write the eigenvalues in the form

\[ \lambda_k = e^{T\rho_k} , \]
where ρk is made unique by choosing its imaginary part to satisfy −π < ℑ(Tρk) ≤ π, so we write
\[ \rho_k = \frac{1}{T}\,\ln \left( \lambda_k \right) \qquad \mbox{or} \qquadT\,\rho_k = \ln \left\vert \lambda_k \right\vert + {\bf j}\,\arg \left( \lambda_k \right) , \]
and ln denoted the principal branch of the natural logarithm. Real numbers ρk are called characteristic exponents. While there are exactly n characteristic multipliers for the periodic linear system \eqref{EqFloquet.1}, there are infinitely many Floquet exponents μ + 2jkπ/T for each integer k ∈ ℤ.

If n×n matrix E has n distinct eigenvalues λk, the Eq.\eqref{EqFloquet.1} has n linearly independent solutions that may be written as

\begin{equation} \label{EqFloquet.3} {\bf u}_k (t) = {\bf p}_k (t)\,e^{\rho_k t} , \end{equation}
where pk(t) is a T-periodic function of time. Equation \eqref{EqFloquet.3} shows that the long time behavior of the solution uk(t) is determined solely by the magnitude of ρk: Since arbitrary solution to the Floquet equation \eqref{EqFloquet.1} is a linear combination of vector functions uk(t), its behavior will depend on which of these functions is missing and what remains. If all |ρk| > 1, then all solutions are unbounded. If all |ρk| < 1, then all solutions tend to zero as t → ∞. If all |ρk| = 1, then all solutions are bounded. Using Liouville--Ostrogradski theorem, we derive
\begin{equation} \label{EqFloquet.4} \det {\bf E} = \exp \left\{ \int \mbox{tr}({\bf A}(t))\,{\text d}t \right\} . \end{equation}
Floquet--Lyapunov Theorem: For the equation \( \dot{\bf x} = {\bf A}(t) {\bf x} , \) where the elements of the n×n matrix A(t) are T-periodic, piecewise continuous functions of t, with a finite number of discontinuities on (−∞, ∞) and integrable at each discontinuity, a fundamental matrix may be expressed in the form
\begin{equation} \label{EqFloquet.5} {\bf \Phi}(t) = {\bf P}(t) \, e^{{\bf K}t} , \end{equation}
where P(t) is a T-periodic, n×n matrix, nonsingular for all t and with elements continuous with integrable, piecewise continuous derivatives. Also, K is a constant n×n matrix.
Remark: If A(t) is real, then P(t) and K are real and
\[ {\bf P}(t+T) = {\bf P}(t)\,{\bf R} , \]
for some real matrix R satisfying R² = I, the identity matrix, and which commutes with K so RK = KR. If R = I, then P(t + T) = P(t), but otherwise P(t + 2T) = P(t).

The representation Φ(t) = P(t)etK in Floquet’s theorem is called a Floquet normal form for the fundamental matrix Φ(t).
Floquet Theorem: If Φ(t) is a fundamental matrix solution of the T-periodic system \( \dot{\bf x} = {\bf A}(t) {\bf x} , \) then, for all t ∈ ℝ
\[ {\bf \Phi}(t+T) = {\bf \Phi}(t)\,{\bf \Phi}^{-1}(0)\,{\bf \Phi}(T) . \]
In addition, there is a matrix K (which may be complex) such that
\[ e^{T{\bf K}} = {\bf \Phi}^{-1}(0)\,{\bf \Phi}(T) \]
Since the function t ↦ A(t) is periodic, it is defined for all t ∈ ℝ. Since the system of equations \eqref{EqFloquet.1} is linear, all solutions of the system are defined for t ∈ ℝ.

If Ψ(t) = Φ(t + T), then Ψ(t) is a matrix solution. Indeed, we have

\[ \dot{\bf \Psi}(t) = \dot{\bf \Phi}(t+T) = {\bf A}(t+T)\,{\bf \Phi}(t+T) = {\bf A}(t)\,{\bf \Phi}(t) , \]
as required.

Define

\[ {\bf C} = {\bf \Phi}^{-1}(0)\,{\bf \Phi}(T) = {\bf \Phi}^{-1}(0)\,{\bf \Psi}(0) , \]
and note that C is nonsingular. The matrix function t ↦ Φ(t)C is clearly a matrix solution of the linear system with initial value Φ(0)C = Ψ(0). By the uniqueness of solutions, Φ(t) = Φ(t)C for all t ∈ ℝ. In particular, we have
\begin{align*} {\bf \Phi}(t+T) &= {\bf \Phi}(t)\,{\bf C} = {\bf \Phi}(t)\, {\bf \Phi}^{-1}(0)\, {\bf \Phi}(T) , \\ {\bf \Phi}(t+2T) &= {\bf \Phi}(t+T +T) = {\bf \Phi}(t+T)\,{\bf C} = {\bf \Phi}(t)\,{\bf C}^2 . \end{align*}
There exists a matrix B, possibly complex, such that
\[ e^{T{\bf B}} = {\bf C} . \]
Also, there is a real matrix R such that
\[ e^{2T{\bf R}} = {\bf C}^2 . \]
If P(t) = Φ(t)e-tB and Q(t) = e-tR, then
\begin{align*} {\bf P}(t+T) &= {\bf \Phi}(t+T)\,e^{-T{\bf B}}\,e^{-t{\bf B}} = {\bf \Phi}(t)\,{\bf C}\,e^{-T{\bf B}}\,e^{-t{\bf B}} = {\bf \Phi}(t)\,,e^{-t{\bf B}} = {\bf P}(t) , \\ {\bf Q}(t+2T) &= {\bf \Phi}(t+2T)\,,e^{-2T{\bf R}}\,e^{-t{\bf R} = {\bf \Phi}(t)\,,e^{-t{\bf R} = {\bf Q}(t) . \end{align*}
Thus, we have
\[ {\bf P}(t+T) = {\bf P}(t), \qquad {\bf Q}(t+2T) = {\bf Q}(t) \]
and
\[ {\bf \Phi}(t) = {\bf P}(t)\, e^{t{\bf B}} = {\bf Q}(t)\, e^{t{\bf R}} , \]
as required.

 

There are two famous typical examples of differential equations with periodic coefficients: Hill's equation
\begin{equation} \label{EqFloquet.6} y'' + \lambda y - q(t)\,y = 0 , \end{equation}
and its particular but very important case,the Mathieu equation
\begin{equation} \label{EqFloquet.7} y'' + \left( \lambda - 2q\,\cos (2t) \right) y (t) = 0 . \end{equation}
In Hill's equation, function q(t) is a periodic continuous function. The latter can be slightly generalized as
\begin{equation} \label{EqFloquet.8} y'' + \left( \lambda - 2q\,\cos (\omega t) \right) y (t) = 0 . \end{equation}
Eq.\eqref{EqFloquet.6} is named after the American astronomer and mathematician George William Hill (1838--1914), who introduced it in 1886. Eq.\eqref{EqFloquet.4} was first introduced in 1868 by the French mathematician Émile Léonard Mathieu (1835--1890), who encountered them while studying vibrating elliptical drumheads. This equation is encountered in many different issues in physics, engineering and industry, including the stability of floating ships and railroad trains, the motion of charged particles in electromagnetic Paul traps, the theory of resonant inertial sensors, and many other problems (see Ruby's article).

Example 3: WE consider an example of Lawrence Marcus and Hidehiko Yamabe:
\[ \frac{{\text d}{\bf x}}{{\text d}t} = {\bf A}(t)\,{\bf x} , \qquad \mbox{where} \qquad {\bf A} = \begin{bmatrix} -1 + \frac{3}{2}\,\cos^2 t & 1 -\frac{3}{4}\,\sin 2t \\ 1 -\frac{3}{4}\,\sin 2t & -1 + \frac{3}{2}\,\sin^2 t \end{bmatrix} \tag{3.1} \]
is a π-periodic matrix-function. Since the trace of the matrix A(t) is −½, the Wronskian of any two solutions is proportional to
\[ W(t) = W_0 e^{-t/2} . \]
However, one solution of Eq.(3.1) is unstable
\[ {\bf x} (t) = e^{t/2} \begin{bmatrix} -\cos t \\ \phantom{-}\sin t \end{bmatrix} . \]
    ◾

 

Hill's equation

A famous example where Floquet theory applies to give good stability results is Hill’s equation \eqref{EqFloquet.6}. It was introduced by George William Hill in his study of the motions of the moon. Roughly speaking, the motion of the moon can be viewed as a harmonic oscillator in a periodic gravitational field. However, this model equation arises in many areas of applied mathematics where the stability of periodic motions is an issue. A prime example is the stability analysis of small oscillations of a pendulum whose length varies with time that leads to the Mathieu equation.

Hill's equation \eqref{EqFloquet.6} is written in a standard Sturm--Liouvulle form for self-adjoint linear differential operator \( L\left[ t, \texttt{D} \right] = \texttt{D}^2 - q(t)\,\texttt{I} , \quad \mbox{with} \quad \texttt{D} = {\text d}/{\text d}t , \) which is negative when q > 0. It is convenient to transfer Hill's equation into a system of first order differential equations. By setting \( y = \dot{x} , \) Hill's equation is written in the standard vector form

\begin{equation} \label{EqFloquet.9} \frac{{\text d}{\bf x}}{{\text d}t} = {\bf A}(t)\,{\bf x}(t) , \qquad {\bf x}(t) = \begin{bmatrix} x(t) \\ y(t) \end{bmatrix} , \qquad {\bf A}(t) = \begin{bmatrix} 0 & 1 \\ q(t) - \lambda & 0 \end{bmatrix} . \end{equation}
Since the trace of the matrix A in Eq.\eqref{EqFloquet.9} is zero, according to Liouville--Ostrogradski theorem, the Wronskian of two linearly independent solutions is 1. We immediately deduce that the determinant of any fundamental matrix is a constant. So the determinant of the monodromy matrix E is 1. Its eigenvalues are roots of the characteristic equation
\begin{equation} \label{EqFloquet.10} \lambda^2 - \mbox{tr}({\bf E}) + 1 = 0 \qquad \Longrightarrow \qquad 2\lambda = \mbox{tr}({\bf E}) \pm \sqrt{\mbox{tr}({\bf E})^2 -4} . \end{equation}
So the long-time behavior of the solutions is determined mainly by the single real number, the trace of the monodromy matrix E.

Let us consider two independent solutions of Hill's equation \eqref{EqFloquet.6} that satisfy mutually complement initial conditions

\[ \xi (0) = 1 , \quad \dot{\xi} (0) = 0 \qquad\mbox{and} \qquad \eta (0) = 0, \quad \dot{\eta}(0) = 1. \]
Then the corresponding fundamental matrix \( {\bf \Phi}(t) = \left[ {\bf \xi} , {\bf \eta} (t) \right] \) satisfies the initial condition &Phi(0) = I and
\[ \mbox{tr}({\bf E}) = \xi (T) + \eta (T) . \]
Then we observe five separate cases according the value of the trace of the monodromy matrix E.
  1. The trace of the monodromy matrix E is greater than 2. The eigenvalues are positive, different not equal to +1 and satisfy 0 < λ1 < 1 < λ2. The characteristic exponents are ±ρ, where
    \[ T\,\rho = \ln \lambda_2 > 0 . \]
    Two linearly independent solutions are
    \[ \xi_1 (t) = e^{-\rho t} p_1 (t) , \qquad \xi_2 (t) = e^{-\rho t} p_2 (t) , \]
    where p1(t) and p2(t) are T-periodic functions.
  2. Tr(E) = 2. the eigenvalues of the monodromy matrix are identical and equal to +1, so ρ = 0. Now the behavior of solutions depends whether matrix E is defective or not.
    • The monodromy matrix E is not defective. Then there are two linearly independent T-periodic solutions:
      \[ \xi_1 (t) = p_1 (t) , \qquad \xi_2 (t) = p_2 (t) , \]
      where p1(t) and p2(t) are T-periodic functions.
    • The monodromy matrix E is defective. Then there is only one eigenvector and two solutions are
      \[ \xi_1 (t) = p_1 (t) , \qquad \xi_2 (t) = t\,p_1 (t) + p_2 (t) , \]
      where p1(t) and p2(t) are T-periodic functions. The first solution ξ1(t) is bounded; the amplitude of the second solution increases linearly with t, so it is unbounded.
  3. The absolute value of the trace of E is less than 2. The eigenvalues of the monodromy matrix are complex and may be written in the form \( \lambda = e^{\pm{\bf j}\theta}, \) 0 < θ < π, with the characteristic exponents \( T\,\rho = \pm \mbox{arccos}(\mbox{tr}(E)/2). Now two independent solutions are
    \[ \xi_1 (t) = e^{{\bf j}\rho t} p_1 (t) , \qquad \xi_2 (t) = e^{-{\bf j}\rho t} p_1 (t) p_2 (t) , \]
    where p1(t) and p2(t) are T-periodic functions. In this case, all solutions are bounded for all times.
  4. Tr(E) = −2. The eigenvalues of the monodromy matrix E are identical and equal to −1. Again the behavior of solutions depends upon whether matrix E is defective or not.
    • The monodromy matrix E is not defective. Then there are two linearly independent T-periodic solutions:
      \[ \xi_1 (t+T) = - \xi_1 (t) , \qquad \xi_2 (t+T) = - \xi_2 (t) . \]
    • The monodromy matrix E is defective. Then there is only one eigenvector and two solutions are
      \[ \xi_1 (t) = p_1 (t) , \qquad \xi_2 (t) = t\,p_1 (t) + p_2 (t) , \]
      where p1(t) and p2(t) are 2T-periodic functions. In this case, the first solution is bounded, and another one is unbounded.
  5. Tr(E) < −2. The eigenvalues of the monodromy matrix E are negative, distinct, not equal to −1, and satisfy &lambda2 < −1 < λ1 < 0. The characteristic exponents are
    \[ T\,\rho = \pm \ln \left( -\lambda_2 \right) + {\bf j}\pi \]
    and the two linearly independent solutions are as in part (I) with p1(t) and p2(t) increasing as t → ∞.
This characterization of stability gives the impression that there is a sharp distinction between stable and unstable solutions. However, this is valid only for long-term behavior as t → ∞.
Example 1a: The differential equation
\[ \texttt{D}^2 x(t) + q(t)\, x(t) = 0 , \tag{2.1} \]
where \( \ddot{x} = \frac{{\text d}^2 x}{{\text d} t^2} = \texttt{D}^2 x(t) \) stands for the second derivative with respect to time variable t, and q(t) is a single-valued periodic function.

One of the basic questions related to the equation (2.1) is if it has periodic solutions, either one class of periodic solutions or all periodic solutions, when we say that the solutions are in co-existence. The necessary condition for periodicity of solutions of the differential equations requires all coefficients to be periodic with the same or commensurable period. Since the coefficient q(x) is periodic, it means that the necessary condition is fulfilled.

If y1(x) is one particular integral of the differential equation (2.1), then a second particular solution is found according to Liouville’s formula

\[ y_1 (x) = y_1 (x) \,\int \frac{{\text d}x}{y_1^2 (x)} . \tag{2.2} \]
The Wronskian determinant W(x) is either identically equal to zero or a constant different from zero for each value of the variable x from the observed interval, then, without reducing the general nature, we can take that W(x) = 1 and obtain the following:
\[ q(x) = y''_1 \cdot y'_2 - y'_1 \cdot y''_2 . \]
Upon dividing by y'2² , we get
\[ \frac{y'_1}{y'_2} = \int \frac{q(x)}{(y'_2)^2} \,{\text d}x . \]
    ◾

 

Mathieu equation

When the angular frequency ω is close to \( 2\sqrt{k}/n , \) where n is a positive integer, the solutions of the Mathieu equation \eqref{EqFloquet.5} grow exponentially in time.
Example 2a: We consider a particular mechanical system that leads to the Mathieu equation.
rod = Graphics[{Black, Thickness[0.012], Line[{{0, 0}, {1, -2}}]}];
bob = Graphics[{Thick, Gray, Disk[{1.11, -2.18}, 0.2]}];
r1 = Graphics[{Gray, Rectangle[{1.7, 0.025}, {0.05, -0.025}]}];
r2 = Graphics[{Gray, Rectangle[{-1.7, 0.025}, {-0.7, -0.025}]}];
circle = Graphics[{Thickness[0.01], Gray, Circle[{0, 0.7}, 0.15]}];
disk = Graphics[{Thick, Red, Disk[{0, 0.7}, 0.07]}];
line1 = Graphics[{Dashed, Line[{{0, 1.0}, {0, -2}}]}];
line2 = Graphics[{Black, Thickness[0.012], Line[{{0, 0}, {-0.156, 0.7}}]}];
line3 = Graphics[{Dashed, Line[{{-0.156, 0.7}, {-0.156, 1.1}}]}];
line4 = Graphics[{Dashed, Line[{{1, -2.18}, {-1.5, -2.18}}]}];
line5 = Graphics[{Dashed, Line[{{1, -2.5}, {-1.5, -2.5}}]}];
arrow = Graphics[{Dashed, Arrowheads[{{-0.05}, {0.05}}], Arrow[{{-1.5, 0}, {-1.5, -2.5}}]}];
ar1 = Graphics[{Gray, Arrowheads[0.03], Arrow[{{0.333, 0.75}, {0.34, 0.7}}]}];
ar = Graphics[{Thick, Gray, Circle[{0, 0.7}, 0.34, {Pi/4, 0}]}];
ar2 = Graphics[{Gray, Arrowheads[0.03], Arrow[{{0, 1}, {-0.156, 1}}]}];
t1 = Graphics[{Black, Text[Style["r", 18], {-0.08, 1.1}]}];
t2 = Graphics[{Black, Text[Style["\[ScriptL]", 18, FontFamily -> "Mathematica1"], {-1.3, -1.3}]}];
t3 = Graphics[{Black, Text[Style["\[ScriptL](1 - cos\[Theta])", 18, FontFamily -> "Mathematica1"], {-0.8, -2.35}]}];
t4 = Graphics[{Black, Text[Style["\[Theta]", 18], {0.13, -0.6}]}];
t5 = Graphics[{Black, Text[Style["m", 18], {1.13, -2.18}]}];
Show[rod, circle, line1, line2, line3, line4, line5, r1, r2, bob, disk, ar, ar1, ar2, t1, t2, t3, t4, t5, arrow]

As the simplest instance of parametric resonance, consider a mass m hanging from a string of length ℓ. Let θ be the angle between the string and the vertical. The pendulum’s equation of motion is
\[ \ddot{\theta} + \frac{g}{\ell}\,\sin\theta = 0 , \tag{2.1} \]
where g is the gravitational acceleration. Suppose that the string is threaded through a hole in a solid plate above the pendulum and that one end of the string is attached to the rim of a wheel with radius r ≪ ℓ, turned by a motor, as shown in Figure above. When the motor runs, the length ℓ of the string varies with time, giving an equation of motion of the form of Eq.\eqref{EqFloquet.5} when θ is small.     ◾

 

Parametric Resonance

It is evident that solutions of the Hill equation are not time-translation invariant, and analytic solutions can be expressed only via special functions or obtained by approximation. For some values of parameters, solutions of the Hill or _Mathieu equation may grow unboundely. This is the phenomenon of parametric resonance. When a mechanical system has at least two vibrating components, the vibration of one of the components may influence the other component. This influence effect which might stabilize or destabilize the system is called autoparametric resonance. Beside purely mechanical systems, the same phenomenon is observed in an electrical system with two coupled resonators. The phenomenon of supplying energy to a system whereby system parameters are varied periodically is known as parametric pumping. Parametric resonance resembles self-oscillation in that the growth of the amplitude of small oscillations is exponential in time, as long as there is some initial perturbation away from the unstable equilibrium at y = 0.

A physical example of parametric resonance is a pendulum with a variable length oscillating in time. This is the example with which you are bound to be familiar---a child on a playground swing. WE can make an observation that the amplitude of the motion is increased by the rhythmical bending and straightening of the child's body with the effect that the center of mass is raised as swing passes through its lowest point and lowered when the swing reaches its highest point. Experimental investigations were carried out by Lord Rayleigh (J. Strutt) in 1880’s using vibrating strings and wave propagation apparatus.

Example 3a: A person sitting on a swing and shifting his/her center of mass up and down with the double frequency is doing exactly this. Near the vertical position (θ= 0) its length being decreased by ∆l, while near turning points it is being increased-by ∆l. Working against gravity, the person is injecting the energy ∆E= 2mg∆l(1−cos(φmax) during the period of motion. The greater is the amplitude φmax, the higher is the energy input. This is why the amplitude is increasing exponentially.

An idealization of this motion is obtained by treating the swing and child as a vertical pendulum with shifts in the center of mass taking place instantaneously at the lowest and the highest points.     ◾

Traditionally, the study of parametric resonance has circled around its control and prevention in regards to structural failure. For example, a ship sailing in a sea may encounter an unexpected roll motion. This roll motion can build up to 40∘ or even more and can cause serious damage to cargo and ship. Parametric resonance is distinct due to an instability phenomenon governed by the Mathieu function.

Example 4a: Let us consider a pendulum with a stiff rod of length ℓ and mass m2 suspended on a spring damper system with mass m1 and stiffness k and damping linear factor d. Mass m1 can move only in the vertical y direction and pendulum only in the x-y plane.
rod = Graphics[{Black, Thickness[0.012], Line[{{0, 0}, {2, -1.25}}]}];
bob = Graphics[{Thick, Gray, Disk[{2.15, -1.32}, 0.2]}];
disk = Graphics[{Thick, Red, Disk[{0, 0}, 0.07]}];
box = Graphics[{LightGray, Rectangle[{1.75, 0.25}, {-1.75, -0.25}]}];
t1 = Graphics[{Black, Text[Style[Subscript[m, 1], 18], {2.15, -1.0}]}];
t2 = Graphics[{Black, Text[Style[Subscript[m, 2], 18], {1.15, 0}]}];
tl = Graphics[{Black, Text[Style["\[ScriptL]", 18, FontFamily -> "Mathematica1"], {0.8, -0.7}]}];
t3 = Graphics[{Black, Text[Style["y(t)", 18, FontFamily -> "Mathematica1"], {-0.3, -1.35}]}];
t4 = Graphics[{Black, Text[Style["\[Theta](t)", 18], {-0.28, 0.6}]}]; line1 = Graphics[{Dashed, Line[{{0, 0}, {-1, 1}}]}];
line2 = Graphics[{Dashed, Line[{{0, 0}, {0, 1}}]}];
line3 = Graphics[{Green, Thickness[0.02], Line[{{-2.0, -2.55}, {2.0, -2.55}}]}];
line4 = Graphics[{Black, Thickness[0.012], Line[{{-1.2, 0}, {-1.2, -1}}]}];
box2 = Graphics[{Gray, Rectangle[{-1.45, -1.14}, {-1.0, -1}]}];
ine5 = Graphics[{Black, Thickness[0.012], Line[{{-1.55, -0.95}, {-1.55, -1.2}, {-0.9, -1.2}, {-0.9, -0.95}}]}];
line6 = Graphics[{Black, Thickness[0.012], Line[{{-1.25, -1.2}, {-1.25, -2.5}}]}];
string = Graphics[{Black, Thickness[0.012], Line[{{1.25, -2.5}, {1.25, -2.0}, {1.0, -1.8}, {1.5, -1.6}, {1.0, -1.4}, {1.5, -1.2}, {1.0, -1.0}, {1.5, -0.8}, {1.25, -0.7}, {1.25, -0.25}}]}];
arrow = Graphics[{Dashed, Arrowheads[0.05], Arrow[{{0, -2.5}, {0, -0.25}}]}];
ar2 = Graphics[{Gray, Thickness[0.02], Arrowheads[0.08], Arrow[{{0, -3}, {0, -2.5}}]}];
t5 = Graphics[{Black, Text[Style[Subscript[F, 0], 18], {0.26, -2.8}]}];
Show[rod, line1, line2, line3, line4, line5, line6, box, box2, \ string, bob, disk, t1, t2, tl, t3, t4, arrow, ar2, t5]
Spring-mass-pendulum.
The governing equations are derived using the Lagrange formalism. Considering the frame of reference at the origin shown in Figure above and defining in Cartesian coordinates first the two degrees of freedom vector ry and rθ
\[ {\bf r}_y = \begin{bmatrix} 0 \\ y(t) \end{bmatrix} \qquad \mbox{and} \qquad {\bf r}_{\theta} = \begin{bmatrix} \ell\,\sin\theta \\ y(t) - \ell\,\cos\theta \end{bmatrix} \tag{4.1} \]
The kinetic energy for both masses are
\[ \mbox{K}_{y} = \frac{1}{2}\,m_2 \left[ \left( \frac{\text d}{{\text d}t} \,{\bf r}_{xy} \right)^2 + \left( \frac{\text d}{{\text d}t} \,{\bf r}_{yy} \right)^2 \right] = \frac{1}{2}\,m_2 \dot{y}^2 \tag{4.2} \]
and
\[ \mbox{K}_{\theta} = \frac{1}{2}\,m_1 \left[ \left( \frac{\text d}{{\text d}t} \,{\bf r}_{x\theta} \right)^2 + \left( \frac{\text d}{{\text d}t} \,{\bf r}_{y\theta} \right)^2 \right] = \frac{1}{2}\,m_1 \ell^2 \cos^2 \theta \,\dot{\theta}^2 + \frac{1}{2}\,m_1 \left( \dot{y} + \ell\,\sin\theta\,\dot{\theta} \right)^2 . \tag{4.3} \]
The potential energies are derived from the same vectors
\[ \Pi_y = m_2 g\, y + \frac{1}{2}\,k\,y^2 \tag{4.4} \]
and
\[ \Pi_{\theta} = m_1 g \left( y - \ell \cos\theta \right) . \tag{4.5} \]
The viscous friction for support mass is
\[ F_{dy} = d\,\dot{y} . \]
The excited driving force is given by
\[ F_{0} = \left( m_1 + m_2 \right) A\,\omega^2 \cos (\omega t) . \]
Applying the Lagrange formalism, we derive the governing system of differential equations
\begin{align*} \left( m_1 + m_2 \right) \ddot{y} + g \left( m_1 + m_2 \right) + k\,y + \ell m_1 \left( \dot{\theta}^2 \cos \theta + \ddot{\theta} \sin\theta \right) &= -d\, \dot{y} + F_0 \tag{4.6} \\ m_1 \ell^2 \ddot{\theta} + g \ell m_1 \sin\theta + \ell m_1 \ddot{y} \sin\theta &= -\delta \dot{\theta} . \tag{4.7} \end{align*}
    ◾
  1. Bukov, M., D’Alessio, L., and Polkovnikov, A., Universal high-frequency behav-ior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics, 64(2):139–226, 2015.
  2. Butikov, Eugene I., Analytical expressions for stability regions in the Ince–Strutt diagram of Mathieu equation, American Journal of Physics, 86, Issue 4, 257 (2018); https://doi.org/10.1119/1.5021895
  3. Giulio Casati and Luca Molinari. “quantum chaos” with time-periodic hamiltonians. Prog. Theor. Phys. Suppl., 98:287, 1989
  4. M. Grifoni and P. Hänggi. Driven quantum tunneling. Physics Reports, 304:229–354, 1998.
  5. Hill, G. W. "On the Part of the Motion of Lunar Perigee Which is a Function of the Mean Motions of the Sun and Moon." Acta Math. 8, 1-36, 1886.
  6. Holthaus, Martin, Floquet engineering with quasienergy bands of periodically driven optical lattices, Journal of Physics B: Atomic, Molecular and Optical Physics, 2015, Volume 49, Number 1, 013001
  7. McLachlan, N.W., Theory of Application of Mathieu Functions, Dover, New York, 1964.
  8. Rand, R.H., Lecture Notes in Nonlinear Vibrations (new edition). Published on-line by The Internet-First University Press (Cornell's digital repository), 2012.
  9. Ruby, L., Applications of the Mathieu equation, American Journal of Physics, 1995, Volume 64, Issue 1 https://doi.org/10.1119/1.18290
  10. V.A. Yakubovich, V.M. Starzhinskii, "Linear differential equations with periodic coefficients" , Wiley (1975) (Translated from Russian)

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions