ME5701 Linear stability analysis of Mathieu equation

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 throughcarefully 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 mathematics and physics. For example, it appears in the analysis of quantum systems with periodic potentials, such aselectrons in a crystal lattice (solid-state physics). Theequation is also used in the study of the stability of orbitsn celestial mechanics, particularly in systems with periodic gravitational forces.Mathematically, the Mathieu equation is a differential equation with periodic coefficients, given byd 2u dt2+ (δ + 2 cost)u = 0where δ and  are real-valued constant parameters, t is the time and u is the unknown. Note that the coefficientcost is 2π-periodic. Therefore, understanding and solving the Mathieu equation are essential for predicting andcontrolling the behaviour of systems subject to periodic influences. The aim of this assignment is to determine thestability 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.

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 expandedwhere i = √ −1 is the imaginary unit, n ∈ (−∞, ∞) is an integer in the Fourier expansion, an, bn are theexpansion coefficients and λ, which will soon become clear that it is the eigenvalue in our stability problemis complex-valued. (Recall that the stability of the system is determined by the real parts of its eigenvaluesλ. Specifically, the system is 代写ME5701 Linear stability analysis of Mathieu equationunstable if at least one eigenvalue has a positive real part. Conversely, thesystem 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 = e arrive at an infinite system of equations whose general form is shown below1where, 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 correspondingcoefficients should be the same; otherwise, the equation won’t hold for arbitrary t! This leads to theunderlined 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 thindex 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 th

(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 revealthe 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, a20and likewise for bnWrite a numerical code to solve this eigenvalue problem for δ =  = 1 and plot the eigenspectrum. Eigenspectrum 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

informationry n ∈ [−40, 40] with the same δ,  to see if the results are converged or not. Discuss your result (what’sthe 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 bydq dt = A(t)q

9)

2where 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 Rt 0A(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 deduceq(T) = q(0)e R

= 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 maydiffer by a constant multiplier, which means y(t) = cq(t) = q(t + T), where c is a constant. From thisequation, by considering t = 0, we have q(T) = cq(0). Comparing this equation with Eq. (10), we havec = 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) . Andthis ratio is exactly c! This means that we can evaluate the value of c to investigate the stability of theoriginal equation. To do so, we need to time-integrate the matrix A from 0 to T because we have the timeintegration RT 0A(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 derivationabove should be presented in a multi-dimensional version. This can be found in the Appendix. Read it anrealise that we need to calculate the monodromy matrix E to get its eigenvalues to reveal the stability ofthe 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, theidentity 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 theeigenvalues 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 thesyntax. 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 Y0should be the columns of the identity matrix  1 0 0 1 and soyouhave 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 alsomeans that the stability criterion, when rephrased using µ, would be that the system is unstable if at leastone eigenvalue µ has a magnitude larger than 1. Conversely, the system is stable if the magnitudes of allthe 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 generatea 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 sampleresultcan be found in the Appendix.2You have to provide an explanation in the report.

• Task (e) To verify your result, we can also look at the time evolution of the solution u directly. Pick twoarbitrary 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.4Appendixorsuggestions and hints:You can use any programming language you are familiar with. Please use double-precision arithmetic inthe computation.

• A late submission will result in a penalty. The complete submission includes a report and all the codescripts. 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 2tructured 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. 2006Computing 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 facilitateyour understanding. The equation isi.e., d dt qk = A(t)qk(t) or qk(T) = e = 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 eachother; or, we can have(t) = X1≤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 wilpropagate in the positivetemporal direction.5Imaginary part of λwhere Y(t) is the matrix accommodating all the yk(t) vectors. In the Floquet theory, E is called monodromymatrix 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 he 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 regionstability of the Mathieu equation (i.e., either all λ < 0 or all |µ| < 1). The blue region denotes theinstability. Such a graph is called neutral curve, which means that the interface between the two coloursindicates the neutral condition that λ = 0 or |µ| = 1.StableUnstable

• In task (e), a sample result for an unstable solution u at  = δ = 0.5 is shown below. The black curveis generated using [t, x] = ode45(0 mathieu0 , [0 6π], [0.2392 0.343], options) in Matlab, where the ODEFUNis called mathieu in my code, [0 6π] is the TSPAN for the black curve, [0.2392 0.343] is some randominitial condition and options have been explained earlier. Thedashed 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 theseparameters is unstable.The instability result is consistent with the neutral curve obtained above.6