# Band-gap structure of Mathieu equation

**URL:** <https://discourse.julialang.org/t/band-gap-structure-of-mathieu-equation/128360>\
**Category:** General Usage\
**Tags:** package, differentialequation\
**Created:** [April 24, 2025, 10:25am UTC](https://discourse.julialang.org/t/band-gap-structure-of-mathieu-equation/128360 "2025-04-24T10:25:41Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![acubed](https://avatars.discourse-cdn.com/v4/letter/a/b19c9b/32.png) [@acubed](https://discourse.julialang.org/u/acubed)\
**Post date:** [April 24, 2025, 10:25am UTC](https://discourse.julialang.org/t/band-gap-structure-of-mathieu-equation/128360/1 "2025-04-24T10:25:41Z")

</div>

I want to plot band-gap structure (also called zone structure) for the Mathieu equation, u''(t)+(B+A\cos t)u(t)=0. According to the Floquet theory, this structure can be obtained by the computation of the monodromy matrix M trace: if |\text{tr}\,M|\>2, one has a gap (solutions are unstable, they are unbounded) and if |\text{tr}\,M|\<2 one has a band (solutions are stable, they are bounded).

To compute the monodromy matrix, we need to find two solutions with linearly independent initials conditions, say (u\_1(0),u\_1'(0))=(1,0) and (u\_2(0),u\_2'(0))=(0,1). Then, we find the matrix M that consists of columns (u\_1(2\pi),u\_1'(2\pi)) and (u\_2(2\pi),u\_2'(2\pi)). And finally we compute |\text{tr}M|. One can realize such scheme using 1st order ODE or 2nd order ODE. I have tried both.

I define the Mathieu equation as

```julia
function Mathieu_Equation_2ndOrder!(ddu, du, u, p, t)
    omega, A, B = p
    V = [[-(B+A*cos(omega*t)), 0];; [0, -(B+A*cos(omega*t))]]
    ddu .= V*u
    nothing
end

```

The matrix `V` is so-called evolution matrix. Then, I find the monodromy matrix and compute its trace as

```julia
function compute_monodromy!(omega, A, B)
    tspan = (0., 2.0*pi/omega)
    p = (omega, A, B)
    
    prob = SecondOrderODEProblem(Mathieu_Equation!, [0,1], [1,0], tspan, p)
    sol = solve(prob, DPRKN6())
    
    monodromy_matrix = hcat(last(sol.u)[1:2], last(sol.u)[3:4])
    trace_monodromy = tr(monodromy_matrix)
    if abs(trace_monodromy)>2
        return 1
    else
        return 0
    end
end

```

I set `omega=1.0`. Then, I vary `A` and `B` by defining 1D arrays `A_values` and `B_values` and compute the trace of monodromy matrix at each point of my grid,

```julia
@time Threads.@threads for i=1:B_size  
    for j=1:A_size
        band_gap[i,j] = compute_monodromy!(omega, A_values[j], B_values[i])
    end
end

```

So, I have `band_gap=1` if it is the gap and `band_gap=0` if is is the band. But I have a surprise:

 ![Julia_Mathieu_Question](https://global.discourse-cdn.com/julialang/original/3X/e/8/e8c4f09df681d932d3f780274e18168ce374bb03.png)

It is know that band-gap structure should be symmetric, but something goes wrong when B\>0. I do not understand what happens. Any manipulations with symplectic integrators, etc produce the same picture. Cf. with the correct band-gap strcuture:

 ![Stability-chart-of-the-Mathieu-equation-White-dark-regions-correspond-to-unstable-stable](https://global.discourse-cdn.com/julialang/original/3X/f/7/f75afd913e88687bcc535d80d2f6c1b2ba0a8bcc.png)