This assignment will guide you to study the stability properties of the Mathieu equation. Please read through

carefully the problem description below and understand the derivations provided.

The Mathieu equation is a type of differential equation that is significant in various fields of applied mathe matics and physics. For example, it appears in the analysis of quantum systems with periodic potentials, such as

electrons in a crystal lattice (solid-state physics). The equation is also used in the study of the stability of orbits

in celestial mechanics, particularly in systems with periodic gravitational forces.

Mathematically, the Mathieu equation is a differential equation with periodic coefficients, given by

d

2u

dt2

+ (δ + 2 cost)u = 0

where δ and  are real-valued constant parameters, t is the time and u is the unknown. Note that the coefficient

cost is 2π-periodic. Therefore, understanding and solving the Mathieu equation are essential for predicting and

controlling the behaviour of systems subject to periodic influences. The aim of this assignment is to determine the

stability of the solution u = 0 as a function of t. (You can easily verify that u = 0 is a fixed point of the equation.)

In the following sections, we will employ the Floquet analysis to investigate the stability of the Mathieu equation.

To complete this CA, you will need to develop numerical codes to complete the specified tasks and

summarize your results in a comprehensive report.

• Task (a) To begin, since the Mathieu equation contains a 2nd-order temporal derivative term, we introduce

v =

du

dt to lower the order of time derivative. Consequently, we have the following equations

du

dt = v (1)

dv

dt + (δ + 2 cost)u = 0 (2)

Rewrite the above equations in a matrix form. You have to complete the following matrix

. (3)

The method to be introduced below is called the Floquet–Fourier–Hill method. You can read Deconinck &

Kutz (2006) for more information. Following this method, the solution form for u and v can be expanded

−1 is the imaginary unit, n ∈ (−∞, ∞) is an integer in the Fourier expansion, an, bn are the

expansion coefficients and λ, which will soon become clear that it is the eigenvalue in our stability problem,

is complex-valued. (Recall that the stability of the system is determined by the real parts of its eigenvalues

λ. Specifically, the system is unstable if at least one eigenvalue has a positive real part. Conversely, the

system is stable if the real parts of all the eigenvalues are negative.)

By substituting the solution forms of u(t), v(t) into Eqs. (1,2) and noting that cost =

arrive at an infinite system of equations whose general form is shown below

λan + inan = bn (4)

λbn + inbn = −(δan + an−1 + an+1) (5)

1You have to provide an explanation in the report.

where, again, n can be any integer from −∞ to ∞. Here, the derivation of Eq. (4) is explained:.

By comparing the LHS and RHS of the last equation above, for each term e

(λ+in)t

, their corresponding

coefficients should be the same; otherwise, the equation won’t hold for arbitrary t! This leads to the

underlined part of the equation, which is (λ + in)an = bn, or Eq. (4). Now, your task is to derive Eq. (5)

using the same reasoning. Note that you will need to rename the index n during the derivation process.

• Task (b) Next, we will write the above system of equations in a matrix form. As usual, we only retain the

λ-related term on the LHS and move all the other terms to the RHS, resulting in from Eqs. (4,5)

λan = −inan + bn, (6)

λbn = −inbn − (δan + an−1 + an+1). (7)

Define a vector q =

bn−1

bn

bn

, where, from top to bottom, we stack a−∞ to a∞, followed by b−∞ to b∞. In

this example, 6 representative terms are explicitly shown.

Now we can write Eqs. (6,7) in a matrix form

(8)

or λq = Mq, where I have denoted the big matrix as M. You have to show the complete matrix in the report.

Apparently, this is an eigenvalue problem for M with λ being the eigenvalue. By solving for λ, we can reveal

the stability of the Mathieu equation around u = 0. To numerically solve the problem, we have to truncate

the value of n and we choose n ∈ [−20, 20], which means that we have a−20, a−19..., a−1, a0, a1, ..., a19, a20

and likewise for bn.

Write a numerical code to solve this eigenvalue problem for δ =  = 1 and plot the eigenspectrum. Eigen spectrum means the set of all the eigenvalues. A sample eigenspectrum is shown in the appendix. In Matlab,

you can use “eig” to solve an eigenvalue problem and obtain the whole spectrum by plotting the real part

of the eigenvalue as the x axis and the imaginary part as the y axis. Do “help eig” in Matlab to get more

information.

Try n ∈ [−40, 40] with the same δ,  to see if the results are converged or not. Discuss your result (what’s

the stability of the equation? what do you observe? what do you find interesting? etc.).

• Task (c) Next, we will use another method to arrive at the same conclusion.

We start with a simple case for the illustration. Consider a scalar function q which is governed by

dq

dt = A(t)q (9)

2

where A(t) = A(t + T) is T-periodic. We can solve this equation using the separation of variables, that is,

from dq

q = A(t)dt, we can get in general q(t) = q(0)e

R

t

0 A(t

0 )dt0 . You can verify this result by substitution.

Remember that this is the only solution. From this solution, by assigning t = T, we can also deduce

q(T) = q(0)e

R

T

0 A(t)dt

. (10)

Then, we realise that y(t) = q(t + T) should also be a solution to the original equation because

dy

dt =

dq(t + T)

dt =

dq(t + T)

d(t + T)

= A(t + T)q(t + T) = A(t)y(t)

which is the same equation as the original one. But y(t), q(t) are not necessarily the same and they may

differ by a constant multiplier, which means y(t) = cq(t) = q(t + T), where c is a constant. From this

equation, by considering t = 0, we have q(T) = cq(0). Comparing this equation with Eq. (10), we have

c = e

R 0

T A(t)dt

. (11)

On the other hand, the stability of Eq. (9) refers to the growth of q over a period T, which is q

q

(

(0)

T)

. And

this ratio is exactly c! This means that we can evaluate the value of c to investigate the stability of the

original equation. To do so, we need to time-integrate the matrix A from 0 to T because we have the time

integration R T

0 A(t)dt in Eq. (11).

In our problem, the matrix A is given in Eq. (3) and q =



u

v



is a vector, not a scalar. Thus, the derivation

above should be presented in a multi-dimensional version. This can be found in the Appendix. Read it and

realise that we need to calculate the monodromy matrix E to get its eigenvalues to reveal the stability of

the equation (the monodromy matrix E is equivalent to c above). The numerical recipe is summarised as

follows