I am about to register PeriodicSystems.jl, a new package dedicated to handle linear control systems with periodic time-varying coefficient matrices. For stability analysis, the computation of characteristic exponents of the state matrix A(t) is a possible approach. I already implemented a function `pseig`

, which computes the characteristic exponents as the logarithms of the eigenvalues of the monodromy matrix (i.e., the state-transition matrix over one period). The following example illustrates the associated computational difficulties and the way to overcome them.

Consider the following periodic matrix of period T = 2pi

`A(t) = [ 0 1; -10*cos(t) -24-10*sin(t)]`

which has characterisitic exponents equal to 0 and -24 (you have to believe me).

However by computing the monodromy matrix and its eigenvalues, one of the resulting characteristic exponents has no **one** exact digit using Float64 based computations with high relative tolerance set for solving the underlying differential equations :

```
julia> using PeriodicSystems
julia> A = PeriodicFunctionMatrix(t -> [0 1; -10*cos(t) -24-10*sin(t)],2pi);
julia> @time ev = psceig(A; reltol = 1.e-10)
0.004462 seconds (21.75 k allocations: 924.010 KiB, 79.31% compilation time)
2-element Vector{Float64}:
-Inf
-1.4052750767605842e-13
```

However, full accuracy can be achieved by determining the mondromy matrix as a product of say 500 matrices and computing the eigenvalues using the periodic Schur decomposition:

```
julia> @time ev = psceig(A,500; reltol = 1.e-10)
0.029436 seconds (329.04 k allocations: 21.171 MiB)
2-element Vector{Float64}:
1.8376538159566882e-15
-23.99999999997412
```

Note that the evaluation of the 500 factors can be done in parallel, in which case a substantial speedup of computations can be achieved.

I recently included the ApproxFun.jl package as an alternative way to define periodic matrices and performed the same computation using operators based computations. I must say, I was pleasantly surprized to obtain almost the same accuracy as before:

```
s = Fourier(0..2π);
a = Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],s);
d = domain(a);
D = Derivative(d);
ND = [D 0I; 0I D];
A = a-ND;
N = 50; n = (2*N+1)*2;
@time ev = eigvals(Matrix(A[1:n,1:n]));
0.044203 seconds (88.73 k allocations: 4.403 MiB)
ev[sortperm(imag(ev),by=abs)][1:4]
4-element Vector{ComplexF64}:
-24.000000000000135 + 0.0im
-23.999999585663453 + 0.0im
-6.042828006006208e-14 + 0.0im
0.0 + 0.0im
```

I guess that the operator A is a Fredholm operator (i.e., sum of a diagonal ND with a Toeplitz operator, but I am not sure) and I only displayed the eigenvalues in the central strip (i.e., in this case, those with null imaginary parts) using a truncated 2*N+1 blocks representation which is typical for Toeplitz operators (not sure if this is the best choice).

My questions are:

- What is the explanation for obtaining the characteristic exponents with double multiplicities?
- There exists a nice compact way to build the diagonal operator ND for an arbitrary order A(t)?