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
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
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,
@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:
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:
