ME5701程序代写、代做Matlab设计编程

日期：2024-11-03 09:53

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