Periodic matrix utilities

 pseig(A, K = 1; lifting = false, solver, reltol, abstol, dt) -> cm

Compute the characteristic multipliers of a continuous-time periodic matrix.

For the given square periodic matrix A(t) of period T, the characteristic multipliers cm are the eigenvalues of the monodromy matrix Ψ = Φ(T,0), where Φ(t,τ) is the state transition matrix satisfying the homogeneous linear ODE

dΦ(t,τ)/dt = A(t)Φ(t,τ),  Φ(τ,τ) = I.

If lifting = false, Ψ is computed as a product of K state transition matrices Ψ = Ψ_K*...*Ψ_1 (see monodromy with the associated keyword arguments). The eigenvalues are computed using the periodic Schur decomposition method of [1]. A may be given as a PeriodicFunctionMatrix, a HarmonicArray, a PeriodicSymbolicMatrix or a FourierFunctionMatrix.

If lifting = true, Ψ is (implicitly) expressed as Ψ = inv(N)*M, where M-λN is a regular pencil with N invertible and the eigenvalues of M-λN are the same as those of the matrix product Ψ := Ψ_K*...*Ψ_1. An efficient version of the structure exploiting fast reduction method of [2] is employed, which embeds the determination of transition matrices into the reduction algorithm. This option may occasionally lead to inaccurate results for large values of K.

References

[1] A. Bojanczyk, G. Golub, and P. Van Dooren, The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems and Control Letters, 50:371-381, 2003.

 cm = pseig(A::PeriodicArray; fast = false)

Compute the characteristic multipliers of a discrete-time periodic matrix A(t). A is given as a PeriodicArray. The characteristic multipliers are computed as the eigenvalues of the reverse product of component matrices, without evaluating the product. If fast = false (default) then the eigenvalues are computed using an approach based on the periodic Schur decomposition [1], while if fast = true the structure exploiting reduction [2] of an appropriate lifted pencil is employed. This later option may occasionally lead to inaccurate results for large number of matrices.

References

[1] A. Bojanczyk, G. Golub, and P. Van Dooren, The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems and Control Letters, 50:371-381, 2003.

cm = pseig(A::PeriodicMatrix[, k = 1]; fast = false)

Compute the characteristic multipliers of a discrete-time periodic matrix A(t). A is given as a PeriodicMatrix. The characteristic multipliers are computed as the eigenvalues of the reverse product of component matrices, without evaluating the product. The argument k specifies the starting index of component matrices (default: k = 1). If fast = false (default) then the eigenvalues are computed using an approach based on the periodic Schur decomposition [1], while if fast = true the structure exploiting reduction [2] of an appropriate lifted pencil is employed. This later option may occasionally lead to inaccurate results for large number of matrices.

Note: The first nmin components of cm contains the core characteristic multipliers, where nmin is the minimum row dimensions of component matrices, while the last ncur-nmin components of cm are zero, where ncur is the column dimension of the k-th component matrix.

References

[1] A. Bojanczyk, G. Golub, and P. Van Dooren, The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems and Control Letters, 50:371-381, 2003.

 psceig(A[, K = 1]; lifting = false, solver, reltol, abstol, dt) -> ce

Compute the characteristic exponents of a continuous-time periodic matrix.

For a given square continuous-time periodic matrix A(t) of period T, the characteristic exponents ce are computed as log.(ev)/T, where ev are the characteristic multipliers (i.e., the eigenvalues of the monodromy matrix of A(t)). For available options see pseig(::PeriodicFunctionMatrix). A may be given as a PeriodicFunctionMatrix, a HarmonicArray, a PeriodicSymbolicMatrix or a FourierFunctionMatrix. For a given square discrete-time periodic matrix A(t) of discrete period N, the characteristic exponents ce are computed as ev.^-N.

psceig(A[, k]; fast = false) -> ce

Compute the characteristic exponents of a cyclic matrix product of p matrices.

Compute the characteristic exponents of a discrete-time periodic matrix A(t). The characteristic exponents are computed from the characteristic multipliers as determined by calling the function pseig(::PeriodicMatrix). If fast = false (default) then the characteristic multipliers are computed using an approach based on the periodic Schur decomposition [1], while if fast = true the structure exploiting reduction [2] of an appropriate lifted pencil is employed. This later option may occasionally lead to inaccurate results for large number of matrices.

The argument k specifies the starting index of the component matrices (default: k = 1).

Note: The first nmin components of ce contains the core characteristic exponents, where nmin is the minimum row dimensions of the component matrices, while the last components of ce are zero.

 pseigsm(PeriodicSymbolicMatrix[, K = 1]; lifting = false, solver, reltol, abstol, dt) -> ev

Compute the characteristic multipliers of a continuous-time periodic symbolic matrix.

For the given square periodic matrix A(t) of period T, the characteristic multipliers ev are the eigenvalues of the monodromy matrix Ψ = Φ(T,0), where Φ(t,τ) is the state transition matrix satisfying the homogeneous linear ODE

dΦ(t,τ)/dt = A(t)Φ(t,τ),  Φ(τ,τ) = I.

If lifting = false, Ψ is computed as a product of K state transition matrices Ψ = Ψ_K*...*Ψ_1 (see monodromy with the associated keyword arguments). The eigenvalues are computed using the periodic Schur decomposition method of [1].

If lifting = true, Ψ is (implicitly) expressed as Ψ = inv(N)*M, where M-λN is a regular pencil with N invertible and the eigenvalues of M-λN are the same as those of the matrix product Ψ := Ψ_K*...*Ψ_1. An efficient version of the structure exploiting fast reduction method of [2] is employed, which embeds the determination of transition matrices into the reduction algorithm. This option may occasionally lead to inaccurate results for large values of K. A has the type PeriodicSymbolicMatrix.

References

[1] A. Bojanczyk, G. Golub, and P. Van Dooren, The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems and Control Letters, 50:371-381, 2003.

psceighr(Ahr::HarmonicArray[, N]; P, nperiod, shift, atol) -> ce

Compute the characteristic exponents of a continuous-time periodic matrix in Harmonic Array representation.

For a given square continuous-time periodic function matrix Ahr(t) of period T in a Harmonic Array representation, the characteristic exponents ce are computed as the eigenvalues of a truncated harmonic state operator A(N)-E(N) lying in the fundamental strip -ω/2 < Im(λ) ≤ ω/2, where ω = 2π/T. If Ahr(t) has the harmonic components A_0, A_1, ..., A_p, then for N ≥ p, P = 1 and nperiod = 1, the matrices A(N) and E(N) are built as

       ( A_0  A_{-1} …  A_{-p}        0    )           ( -im*ϕ_{-N}I                                 0        )
       ( A_1   A_0             ⋱           )           (     ⋮       ⋱                                        )
       (  ⋮         ⋱            ⋱         )           (               -im*ϕ_{-1}*I                           )
A(N) = ( A_p             ⋱          A_{-p} ) , E(N) =  (                           -im*ϕ_{0}*I                )
       (        ⋱           ⋱         ⋮    )           (     ⋮                                  ⋱              )
       (  0        A_p      …         A_0  )           (     0                                   -im*ϕ_{N}I   )

with ϕ_{i} := shift+i*ω. If N < p, then a truncated full block Toeplitz matrix A(N) is built using the first N harmonic components. The default value used for N is N = max(10,p-1).

Generally, for given P ≥ 1 and nperiod ≥ 1, the block Toeplitz matrix A(N) (and also E(N)) is constructed with (2N*np+1)×(2N*np+1) blocks, with np = P*nperiod, such that each A_i is preceeded in its column by np-1 zero blocks, each A_{-i} is preceeded in its row by np-1 zero blocks and all diagonal blocks are equal toA_0.

The keyword argument atol (default: atol = 1.e-10) is a tolerance on the magnitude of the trailing components of the associated eigenvectors used to validate their asymptotic (exponential) decay. Only eigenvalues satisfying this check are returned in ce.

References

[1] J. Zhou, T. Hagiwara, and M. Araki. Spectral characteristics and eigenvalues computation of the harmonic state operators in continuous-time periodic systems. Systems and Control Letters, 53:141–155, 2004.

psceigfr(A::FourierFunctionMatrix[, N]; P, atol) -> ce

Compute the characteristic exponents of a continuous-time periodic matrix in Fourier Function Matrix representation.

For a given square continuous-time periodic function matrix A(t) of period T in a Fourier Function Matrix representation, the characteristic exponents ce are computed as the eigenvalues of the state operator A(t)-D*I lying in the fundamental strip -ω/2 < Im(λ) ≤ ω/2, where ω = 2π/T. A finite dimensional truncated matrix of order n*(2*N*P+1) is built to approximate A(t)-D*I, where n is the order of A(t), N is the number of selected harmonic components in the Fourier representation and P is the period multiplication number (default: P = 1). The default value used for N is N = max(10,p-1), where p the number of harmonics terms of A(t) (see FourierFunctionMatrix).

The keyword argument atol (default: atol = 1.e-10) is a tolerance on the magnitude of the trailing components of the associated eigenvectors used to validate their asymptotic (exponential) decay. Only eigenvalues satisfying this check are returned in ce.

 psceigsm(A::PeriodicSymbolicMatrix[, K = 1]; lifting = false, solver, reltol, abstol, dt) -> ce

Compute the characteristic exponents of a periodic symbolic matrix.

For a given square continuous-time periodic function matrix A(t) of period T, the characteristic exponents ce are computed as log.(ev)/T, where ev are the characteristic multipliers (i.e., the eigenvalues of the monodromy matrix of A(t)). For available options see pseig(::PeriodicFunctionMatrix).

 monodromy(A[, K = 1]; solver, reltol, abstol, dt) -> Ψ::PeriodicArray

Compute the monodromy matrix for a linear ODE with periodic time-varying coefficients.

For the given square periodic matrix A(t) of period T and subperiod T′ = T/k, where k is the number of subperiods, the monodromy matrix Ψ = Φ(T′,0) is computed, where Φ(t,τ) is the state transition matrix satisfying the homogeneous linear ODE

dΦ(t,τ)/dt = A(t)Φ(t,τ),  Φ(τ,τ) = I.

A may be given as a PeriodicFunctionMatrix, a HarmonicArray, a PeriodicSymbolicMatrix or a FourierFunctionMatrix.

If K > 1, then Ψ = Φ(T′,0) is determined as a product of K matrices Ψ = Ψ_K*...*Ψ_1, where for Δ := T′/K, Ψ_i = Φ(iΔ,(i-1)Δ) is the state transition matrix on the time interval [(i-1)Δ,iΔ].

The state transition matrices Φ(iΔ,(i-1)Δ) are computed by integrating numerically the above homogeneous linear ODE. The ODE solver to be employed can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-3), absolute accuracy abstol (default: abstol = 1.e-7) and/or the fixed step length dt (default: dt = min(Δ, Δ*K′/100)) (see tvstm). For large values of K, parallel computation of factors can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

 monodromy(A::PeriodicTimeSeriesMatrix) -> Ψ::PeriodicArray

Compute the monodromy matrix for a continuous-time time series matrix.

For the given square periodic matrix A(t) of period T and subperiod T′ = T/k, where k is the number of subperiods, the monodromy matrix Ψ = Φ(T′,0) is computed, where Φ(t,τ) is the state transition matrix satisfying the homogeneous linear ODE

dΦ(t,τ)/dt = A(t)Φ(t,τ),  Φ(τ,τ) = I.

A is a PeriodicTimeSeriesMatrix with K component matices and the resulting monodromy matrix Ψ is stored as a discrete-time periodic array with K component matrices, of period T and k subperiods.

Ψ = Φ(T′,0) is determined as a product of K matrices Ψ = Ψ_K*...*Ψ_1, where for Δ := T′/K, Ψ_i = Φ(iΔ,(i-1)Δ) is the state transition matrix on the time interval [(i-1)Δ,iΔ]. Each state transition matrix is computed as a matrix exponential Φ(iΔ,(i-1)Δ) = exp(A[i]*Δ), where A[i] is the i-th component matrix of the time series representation. For large values of K, parallel computation of factors can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

 tvstm(A, tf, t0; solver, reltol, abstol, dt) -> Φ

Compute the state transition matrix for a linear ODE with periodic time-varying coefficients. For the given square periodic continuous-time matrix A(t), initial time t0 and final time tf, the state transition matrix Φ(tf,t0) is computed by integrating numerically the homogeneous linear ODE

  dΦ(t,t0)/dt = A(t)Φ(t,t0),  Φ(t0,t0) = I

on the time interval [t0,tf]. A may be given as a PeriodicFunctionMatrix, a HarmonicArray, a PeriodicSymbolicMatrix or a FourierFunctionMatrix.

The ODE solver to be employed can be specified using the keyword argument solver (see below), together with the required relative accuracy reltol (default: reltol = 1.e-3), absolute accuracy abstol (default: abstol = 1.e-7) and/or the fixed step length dt (default: dt = tf-t0). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "linear" - use a special solver for linear ODEs (MagnusGL6()) with fixed time step dt;

solver = "symplectic" - use a symplectic Hamiltonian structure preserving solver (IRKGL16());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

 tvmeval(At::PeriodicTimeSeriesMatrix, t; method = "linear") -> A::Vector{Matrix}

Evaluate the time response of a periodic time series matrix.

For the periodic time series matrix At and the vector of time values t, an interpolation/extrapolation based approximation A[i] is evaluated for each time value t[i]. The keyword parameter method specifies the interpolation/extrapolation method to be used for periodic data. The following interpolation methods from the Interpolations.jl package can be selected:

method = "constant" - use periodic B-splines of degree 0;

method = "linear" - use periodic B-splines of degree 1 (periodic linear interpolation) (default);

method = "quadratic" - use periodic B-splines of degree 2 (periodic quadratic interpolation);

method = "cubic" - use periodic B-splines of degree 3 (periodic cubic interpolation).

 tvmeval(Ahr::HarmonicArray, t; ntrunc, exact = true) -> A::Vector{Matrix}

Evaluate the time response of a harmonic array.

For the harmonic array Ahr representing representing a continuous-time time periodic matrix A(t) and the vector of time values t, A[i] = A(t[i]) is computed for each time value t[i]. If exact = true (default) an exact evaluation is computed, while for exact = false, a linear interpolation based approximation is computed (potentially more accurate in intersample points). The keyword argument ntrunc specifies the number of harmonics to be used for evaluation (default: maximum possible number of harmonics).

 tvmeval(A, t) -> Aval::Vector{Matrix}

Evaluate the time response of a periodic matrix.

For the periodic matrix A(t) and the vector of time values t, the vector Aval of time values is computed such that Aval[i] = A(t[i]) for each time value t[i].

 tvmeval(Asym::PeriodicSymbolicMatrix, t) -> A::Vector{Matrix}

Evaluate the time response of a periodic symbolic matrix.

For the periodic symbolic matrix Asym representing a time periodic matrix A(t) and the vector of time values t, A[i] = A(t[i]) is evaluated for each time value t[i].

 hreval(Ahr::HarmonicArray, t; ntrunc, exact = true) -> A::Matrix

Evaluate the harmonic array Ahr representing a continuous-time time periodic matrix A(t) for a time value t. For a real value t, if exact = true (default) an exact evaluation is computed, while for exact = false, a linear interpolation based approximation is computed (potentially more accurate in intersample points). The keyword argument ntrunc specifies the number of harmonics to be used for the evaluation (default: maximum possible number).

 hreval(Ahr::HarmonicArray, t; ntrunc, exact = true) -> A::Matrix

Evaluate the harmonic array Ahr representing a continuous-time time periodic matrix A(t) for a symbolic time value t. A symbolic evaluation of A(t) is performed (see also hr2psm) The keyword argument ntrunc specifies the number of harmonics to be used for the evaluation (default: maximum possible number).

 hrchop(Ahr::HarmonicArray; tol) -> Ahrtrunc::HarmonicArray

Remove the trailing terms of a harmonic representation by deleting those whose norms are below a certain tolerance.

 hrtrunc(Ahr::HarmonicArray, n) -> Ahrtrunc::HarmonicArray

Truncate a harmonic representation by deleting the trailing terms whose indices exceed certain number n of harmonics.

pmaverage(A; rtol = sqrt(eps())) -> Am

Compute for the continuous-time periodic matrix A(t) the corresponding time averaged matrix Am over one period. The Gauss-Konrod quadratue method is employed for numerical integration using a relative accuracy tolerance specified by rtol.

pmaverage(A::PeriodicSymbolicMatrix) -> Am

Compute for the continuous-time periodic matrix A(t) the corresponding time averaged matrix Am over one period.