Consider the n-dimensional linear system of ordinary differential equations: $$x'(t)=A(t)x(t)$$ where in addition, the matrix $A(t)$ is periodic with period $T$, i.e. $A(t+T)=A(t)$ for all real values of $t$. Let $Φ(t)$ denote the principal fundamental matrix solution such that the columns of $Φ(t)$ are linearly independent, and $Φ(0)=I$. Then, Floquet's theorem decomposes the principal fundamental matrix as the product of a periodic matrix $P(t)$ with period $T$ and an exponential matrix $e^{tB}$.



In [1]:

    
using Plots,ApproxFun



In [3]:

    
function Base.expm{T<:Number}(A::Matrix{T},f::Fun)
    D,V = eig(A)
    fC = D == real(D) ? fill(f,n) : fill(f+0im,n)
    for i = 1:length(D)
        fC[i] = Fun(t->exp(D[i]*fC[i](t)),domain(f))
    end
    V*diagm(fC)*inv(V)
end;

Consider the system of coupled harmonic oscillators with periodic parameterized excitation: $$x'' + (1+a \cos(2t))x = x-y$$ $$y'' + (1+a \cos(2t))y = y-x$$ We calculate the Principal Fundamental Matrix $Φ$.



In [4]:

    
T = π;a=0.15
t = Fun(identity,[0,T])
d=domain(t)
D=Derivative(d)
B=ldirichlet(d)
n=4
A=[ diagm(fill(ldirichlet(d),n));
    D             -I  0I            0I;
   (2+a*cos(2t))   D  -I            0I;
   0I             0I   D            -I;
   -I             0I  (2+a*cos(2t))  D]

Φ = A\eye(n);

It can be further expressed as $Φ(t) = P(t) e^{tB}$ where $P(t)$ is a periodic matrix, and $B$ is constant. Using the values of $Φ(T)$:



In [5]:

    
ΦT = Φ(T)









    Out[5]:





4x4 Array{Float64,2}:
 -0.170879  -0.148885  -0.836059   0.265569
  0.732284  -0.170879  -0.612945  -0.836059
 -0.836059   0.265569  -0.170879  -0.148885
 -0.612945  -0.836059   0.732284  -0.170879

we calculate the matrix $B$ and the Floquet exponents and multipliers.



In [6]:

    
B = logm(ΦT)/T
D,V = eig(B)

println("Floquet Exponents = \n")
[println(D[i]) for i=1:n]
println("\nFloquet Multipliers = \n")
[println(exp(D[i]*T)) for i=1:n];









    



WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix logarithm will be returned.






    



Floquet Exponents = 

-1.6653345369377348e-16 - 0.2683546905319053im
-2.6075062248276382e-17 + 0.2683546905319057im
0.03747531973228921 + 0.9999999999999994im
-0.03747531973228907 + 1.0000000000000002im

Floquet Multipliers = 

0.6651802570066602 - 0.746682814646589im
0.6651802570066595 + 0.7466828146465901im
-1.1249427991479368 + 2.136065595020538e-15im
-0.888934086921956 - 6.806690441197217e-16im

We can build the matrix exponential $e^{-tB}$ and use it to find the periodic matrix $P(t)$. Periodicity is numerically confirmed by Fun instantiation with low number of coefficients.



In [7]:

    
PI = Φ*expm(B,-t)
P = Fun(t->PI(t),PeriodicInterval(d))
P|>ncoefficients









    Out[7]:





208

With the matrix B and the periodic matrix P(t), we can construct the solution with any initial conditions for as long as we want!



In [9]:

    
t = Fun(identity,[0,10T])
d = domain(t)
x0 = rand(n)

PI = Fun(t->P(t),d)
xsol = real(PI*expm(B,t)*x0)

plot(mat(xsol))









    



[Plots.jl] Initializing backend: gr






    Out[9]: