Assignment for Part 2 in ME5701
——Linear stability analysis of Mathieu equation——
Due time: 23:59:59, Nov. 15th, 2024
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
– (I) From Eq. (12), by setting t = 0, we can have Q(T) = EQ(0). One can assume Q(0) = I, the
identity matrix (why?2
).
– (II) Time-integrate the original equation d
dt
Q = A(t)Q from the initial condition Q(0) = I to get Q(T).
The results will be our monodromy matrix E, according to (I). How to time-integrate an equation?
See the explanation below.
– (III) Calculate the eigenvalue µ of E and do the transformation 2
1
π
log µ, which should be equal to the
eigenvalues obtained in task (b). You have to verify this in the report.
In Matlab, you can use ode45 to time-integrate an ODE. Do “help ode45” in Matlab to understand the
syntax. You will see that ode45 can be executed using
[T OUT, Y OUT] = ode45(ODEF UN, T SP AN, Y 0, OP T IONS).
ODEF UN is to implement Eqs. (1,2). The TSPAN denotes the time span and should be [0, 2π]. The Y0
should be the columns of the identity matrix  1 0
0 1 and so you have to execute ode45 twice with Y0 being

1
0

and  0
1

, respectively. Then you assemble the two resultant columns to form the monodromy matrix.
Follow the above numerical recipe to get the final eigenvalues.
Note that the eigenvalue λ in task (b) is connected to the eigenvalue µ in task (c) by λ = 2
1
π
log µ. This also
means that the stability criterion, when rephrased using µ, would be that the system is unstable if at least
one eigenvalue µ has a magnitude larger than 1. Conversely, the system is stable if the magnitudes of all
the eigenvalues are less than 1. In Matlab, you can use “abs” to get the magnitude of a complex number.
• Task (d) Using either method explained above, do a parametric study of δ, . You can for example generate
a graph showing the stability or instability of the system on the δ −  plane with δ ∈ [0, 2] and  ∈ [0, 1].
Discuss what you find and explain the results.
A sample result can be found in the Appendix.
2You have to provide an explanation in the report.
3
• Task (e) To verify your result, we can also look at the time evolution of the solution u directly. Pick two
arbitrary sets of parameters with δ ∈ [1, 2] and  ∈ [0.5, 1], one being stable and the other being unstable.
In each case, plot u as a function of t and explain what you find.
A sample result can be found in the Appendix.
4
Appendix for suggestions and hints:
• You can use any programming language you are familiar with. Please use double-precision arithmetic in
the computation.
• A late submission will result in a penalty. The complete submission includes a report and all the code
scripts. The code should be executable and generate the graphs to be presented in the report once executed
(this will facilitate my checking of your code).
• This is a group assignment with a maximum of 2 students per group, with only one submission per group.
You can also do it individually.
• This CA makes up 20% of the final mark. (The remaining 30% of the final mark for Part 2 will be the 2
structured questions in the quiz. Similarly, Part 1 will make up the other 50%.)
• In task (a), the complete citation of Deconinck & Nathan (2006) is Deconinck B. & Nathan Kutz, J. 2006
Computing spectra of linear operators using the Floquet–Fourier–Hill method. Journal of Computational
Physics 219 (1), 296–321.
• A sample result for the eigenspectrum in task (b) is shown below. Note that the real part of the eigenvalue
λ is plotted as the x axis and imaginary part the y axis.
Real part of λ
• A multi-dimensional derivation of the Floquet theory in Task (c) is provided below, which may facilitate
your understanding. The equation is
dq
dt = A(t)q, where A(t + T) = A(t)
The linear equation admits a solution of matrix Q(t), that is d
dt
Q = A(t)Q. The columns qk(t) are solutions
to the above equation, i.e., d
dt
qk = A(t)qk(t) or qk(T) = e
R
T
0 A(t)dtqk(0).3 Consider yk(t) = qk(t + T), we
have
dyk
dt =
dqk(t + T)
dt =
dqk(t + T)
d(t + T)
= A(t + T)qk(t + T) = A(t)yk(t)
which is the same equation as the original one. This means yk and qk are a linear combination of each
other; or, we can have
yk(t) = X
1≤j≤n
ek,jqj (t) or Y(t) = EQ(t) or Q(t + T) = EQ(t) (12)
3Note that a time-ordering operator is omitted here since it is implicitly understood that the equation will propagate in the positive
temporal direction.
5
Imaginary part of λ
where Y(t) is the matrix accommodating all the yk(t) vectors. In the Floquet theory, E is called monodromy
matrix which represents the effect of the operator A over one period (i.e. the linearized Poincar´e map).
Thus, a state Q(t) left-multiplied by E will propagate the former to a later state at t + T, that is Q(t + T).
Besides, you can easily understand the matrix E here is the coefficient c in the scalar version of the derivation.
• In task (c), ode45 in Matlab also provides an option for setting up the tolerance criterion for convergence;
for this, you can consider options = odeset(
0 AbsT ol0 , 1e − 7,
0 RelT ol0 , 1e − 8).
• In task (d), a sample result for the stability/instability in the δ −  plane with δ ∈ [0, 1] and  ∈ [0, 0.5]
is provided below. The plot is generated using the “contourf” command in Matlab. The yellow region
indicates stability of the Mathieu equation (i.e., either all λ < 0 or all |µ| < 1). The blue region denotes the
instability. Such a graph is called neutral curve, which means that the interface between the two colours
indicates the neutral condition that λ = 0 or |µ| = 1.
Stable
Unstable
• In task (e), a sample result for an unstable solution u at  = δ = 0.5 is shown below. The black curve
is generated using [t, x] = ode45(0 mathieu0 , [0 6π], [0.2392 0.343], options) in Matlab, where the ODEFUN
is called mathieu in my code, [0 6π] is the TSPAN for the black curve, [0.2392 0.343] is some random
initial condition and options have been explained earlier. The dashed vertical lines indicate 2π, 4π, 6π,
corresponding to 1,2,3 periods. To generate the green and the red curves, you change 6π to 4π and 2π
respectively. You can clearly see that the solution u is amplifying over the periods, so the equation at these
parameters is unstable.
The instability result is consistent with the neutral curve obtained above.
6

请加QQ：99515681  邮箱：99515681@qq.com   WX：codinghelp