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
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 behaviour 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 funcion, 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]
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 eperiodic systems

The general theory of n-dimensional linear, periodic systems is very similar to the theory of one-dimensional linear periodic equations. We considerthe 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} = {\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) 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 ∈ ℝ).
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.2} {\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).


There are two famous typical examples of differential equations with periodic coefficients: Hill's equation
\begin{equation} \label{EqFloquet.3} y'' + q(t)\,y = 0 \end{equation}
and its particular but very important case,the Mathieu equation
\begin{equation} \label{EqFloquet.4} y'' + k \left( 1 - 2m\,\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.5} y'' + k \left( 1 - 2m\,\cos (\omega t) \right) y (t) = 0 . \end{equation}
Eq.\eqref{EqFloquet.3} 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).


Hill's equation

By defining \( y = \dot{x} , \) we may convert Hill's equation into the standard vector form
\begin{equation} \label{EqFloquet.6} \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) & 0 \end{bmatrix} . \end{equation}
Since the trace of the matrix A in Eq. is zero, according to Liouville--Ostrogradski theorem, the Wronskian of two linearly independent solutions is 1.
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 ,i.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 increasedby ∆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 idealisation 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}]}];
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_1 \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_1 \dot{y}^2 \tag{4.2} \]
\[ \mbox{K}_{\theta} = \frac{1}{2}\,m_2 \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_2 \ell^2 \cos^2 \theta \,\dot{\theta}^2 + \frac{1}{2}\,m_2 \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_1 g\, y + \frac{1}{2}\,k\,y^2 \tag{4.4} \]
\[ \Pi_{\theta} = m_2 g \left( y - \ell \cos\theta \right) . \tag{4.5} \]
The viscous friction for support mas 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_2 \left( \dot{\theta}^2 \cos \theta + \ddot{\theta} \sin\theta \right) &= -d\, \dot{y} + F_0 \tag{4.6} \\ m_2 \ell^2 \ddot{\theta} + g \ell m_2 \sin\theta + \ell m_2 \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 eqution, American Journal of Physics, 86, Isuue 4, 257 (2018);
  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
  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