Finite Dimensional Approximation of the Monodromy Operator of a Periodic Delay Differential Equation with Piecewise Constant Orthonormal Functions

Using piecewise constant orthonormal functions, an approximation of the monodromy operator of a Linear Periodic Delay Differential Equation (PDDE) is obtained by approximating the integral equation corresponding to the PDDE as a linear operator over the space of initial conditions. This approximation allows us to consider the state space as finite dimensional resulting in a finite matrix approximation whose spectrum converges to the spectrum of the monodromy operator.


Introduction
Linear Periodic Differential Equations (PDDEs) have been of importance for studying problems of vibration, mechanics, astronomy, electric circuits, biology among others in [1] several examples of delay effects on mechanical systems are given, in [2] and [3] effects of the delay in physics and biological processes are considered.Neglecting the fact that interaction between particles does not occur instantaneously sometimes is no longer possible or practical, these finite velocity interactions bring new behaviors that modify significantly the behavior of the system, see for example [4] and [5].In the study of these Delay Differential Equations many problems arise, mainly due to the infinite dimensional nature of the system, in the case of linear PDDEs the stability depends on the spectrum of the monodromy operator, which in the non delayed case corresponds to the monodromy matrix.
Approximation methods of the monodromy operator have been proposed a number of times, [6] and [7] make use of pseudospectral collocation methods to approximate the monodromy operator.The well known method of semidiscretization [8] has also been used to determine the stability of a PDDE.In [9] using a Walsh approximation method from [10] a set of approximated solutions of a PDDE was used to construct an approximation of the monodromy operator by numerical integration.
In this work, the main contribution is an approximation of the monodromy operator of the PDDE by a linear Equation (35) of the form 0 k x x =  , where the directly obtained matrix k  will correspond to the approximated monodromy operator, with no need of approximating solutions or numerical integration.Stability of the PDDE can then be determined by the spectrum of k  without any need of solving any equation.Convergence of k  and its spectrum is stated in Theorems 10, 14 and 15.This approximation is made by projecting the integral equation corresponding to the PDDE in to a a finite dimensional subspace spanned by finitely many piecewise constant functions.The utilized functions must belong to a complete set of piecewise constant orthonormal functions with discontinuities only in dyadic rationals, such as Haar, Walsh or Block Pulse Functions (BPF).The theoretical framework of this paper will be based on Walsh functions since most results are stated for these functions.Once obtained the finite dimensional approximation of the monodromy operator will be stated in terms of BPF to reduce the computational cost of the stability analysis.
The main goal of this paper is to provide a computationally light method, with straightforward implementation, to approximate the monodromy operator of a delay differential equation, in order to facilitate the computation of stability diagrams used to study the behavior of the equation with respect to changes in its parameters.
A solution of (1) with initial condition ϕ ∈  at 0 t is understood as a mapping [ ] ) ( ) .The operator T is called the solution mapping and is analogous to the state transition matrix of the undelayed case.

Monodromy Operator
Taking into account the periodic nature of (1) it is relevant to consider (2) when 0 t t ω = + , in this case we will denote the solution mapping as ( ) ( ) Next some properties consequence of ( ) U ⋅ being completely continuous are enlisted, [12]: , and if ( ) 0 U henceforth denoted simply as U will be called the monodromy operator of (1) and the elements

Stability
It was shown [12] that Floquet theory is valid for PDDEs in a certain sense, and that stability of (1) will be related to the spectrum of U. The following theorem from [13] establishes the conditions for the stability of ( 1

Definition
Walsh functions are a set of piecewise constant complete orthonormal functions introduced in [16] and defined on the interval [ ) 0,1 , although easily translated to any other interval.Formal definition of the Walsh functions may be done in many ways [17], with the use of Rademacher functions the -th Walsh function may be defined as [17]: where the i ε are the coefficients of the unique binary expansion of , namely { } 0 2 , with 0,1 There are as well many ways of ordering the Walsh functions, in this paper we will use the so called Dyadic ordering [18].Figure 1 shows the first eight Walsh functions in this ordering.

To a function
f L ∈ will correspond the Walsh approximation: where ( ) ( ) The right side of (5) will converge in norm to f for any 2 k n = , it will converge uniformly for any continuous function that has, at most, jump discontinuities at dyadic rationals (numbers of the form , , 2 b a a b ∈ ∈  ) [16].When n = ∞ , (5) will be called a Walsh expan- sion.

Properties of Walsh Functions
We will define a Walsh matrix

[ ]
W k as the 2 2 k k × matrix consisting on the discretization values of the Walsh functions over the interval [ ] 0,1 .For example we will have for We can also define the Delayed Walsh Matrix [ ] [ ] We define the vector consisting on the first 2 k Walsh functions as: Since the order of the approximation should be clear from context we write simply ( ) w t .We define as well the dyadic sum ⊕ between tho nonnegative integers as: where i q and i r are respectively the coefficients of the binary expansions of q and r as in (4).To a vector whose i, jth entry will be given by σ σ σ σ σ = , then: ( ) We make use of the following Lemma: The integral of a Walsh vector can be approximated by a Walsh expansion of the form: where P is given by [20]: (11) The approximated integral converges to the exact integral uniformly [21].
To a function matrix ( ) A t we associate the Walsh approximation ( ) where where is the vector of Walsh coefficients of the approximation of each element of f.Define  , then from (9) we have: where: .

Approximation Error
Regarding the error of Walsh approximation we have: Lemma 3. [17] If f satisfies the Lipschitz condition, then: for some constant C.

Block Pulse Functions
The set of Block pulse functions of order p is defined on the interval [ ) , where [22]: ( )

0, otherwise
Block Pulse Functions, shown on Figure 2, are orthogonal, easily normalizable [23] and when m → ∞ they form a complete set [22].We restrict ourselves to the case 2 k p = for k integer.Walsh functions and Block Pulse Functions are related by a one on one correspondence, meaning that there exists a unique bijective linear transformation that maps the first 2 k Walsh functions on to the set of Block Pulse Functions of order 2 k [24].The existence of this transfor- mation and the completeness of the Block Pulse functions ensure that the properties of Walsh functions are inherited to Block Pulse Functions.In particular we have the following: Lemma 4. [24] Matrix P in ( 10) is similar to the upper triangular matrix , where Q is a nilpotent matrix with ones above the diagonal and zeros everywhere else: σ ∈  being the vector of coefficients of the Walsh approximation.Then ( ) σ Λ is similar to a diagonal matrix:

Approximation of the Monodromy Operator
Without loss of generality we assume 1 ω = .We approximate the monodromy operator by projecting (1) on a finite dimensional subspace of [ ) 2 0,1 L formed by the span of 2 k piecewise constant orthonormal functions.We will assume the orthonormal functions to be Walsh functions, the analysis can be carried to the case of Block Pulse Functions or Haar functions by means of similarity transformations.We also restrict ourselves to the case of commensurable delays, that is mτ ω = for some m ∈  .Integrating (1) from 0 to t we will have: with ( ) ( ) . A solution of (20) will correspond to a solution of (1) [26].Let k π denote the projection mapping that takes the Walsh expansion ( ) ( ) to the Walsh approximation . Along with (20) we introduce its projection onto the space where ( ) ( ) ( ) , this is, the approximations of ( ) ( ) B t , respectively, as in (12).( )( ) cant define yet the projection of the initial condition since its defined on a different domain than the domain of definition of the Walsh functions.For this we have: Proposition 6.The value of the Walsh approximation of order 2 k at an interval 1 , 2 2 τ − and 0 everywhere else, and we make ( ) ( ) . Since Walsh functions were not defined at We split the second integral in (20) as: where ( ) S t is the unit step function.We make use of the following Lemma: Lemma 7. [17] Any function ( ) i.e. the non zero coefficients of the Walsh expansion of ( ) From linearity of k π and from Lemma 7. we have that we can split the second integral in (21) as in (22) with this being consistent with the projection: Since ( )  with each i h and i c being vectors of Walsh coefficients.Defining  , and recalling ( 7) and ( 12), we can write (21) as: From Lemma 7. we have that both, ( ) ( ) given in an analogous way as . The term ( ) will also have an unique representation: since ( ) ( ) From ( 14), ( 27) and ( 26) we have that we can write (25) as: From (10) we have that: ( ) ( ) where matrix given by: T T T 0 0 0 0 , 0 0 from which we arrive at: is nonsingular we have: From (32) we can define a mapping However by construction we can see that the domain of k  is in fact the subspace of M consisting of al functions generated by linear combinations of the first 2 k Walsh functions, that are equal to 0 for  is an approximation of the solution mapping T of (1), we obtain an approximated solution ( ) ( ) which satisfies (21), and thus, we have an approximation of the solution map of (1).
In order to study the monodromy operator of (1), we must study the state at t ω = , this is, we must know the solution of (

( )
T t H Ω into the subspace M ′ , with diagonal entries T P W .In the case of Walsh functions, we will have are given by Walsh matrices of order k that instead have 0s in the first (1 − m)th columns.
From (34) we obtain our approximated monodromy operator, given by:

Approximation by Block Pulse Functions
Equation (35) for the approximation of k  is valid for all sets of orthonormal functions that can be obtained by linear combinations of Walsh functions.For the numerical calculation of the approximation of the monodromy operator, Block Pulse Functions are the most advantageous set since with their simplicity comes a lesser computational load.In particular for the Block Pulse Function approximation we will have: , , , , , P diag P P = (38) If we define: we will have: Finally, if in (15) we have ( ) will be given as in (15) with: ( ) where ( )

Convergence of k 
We denote as ∆  the space of functions from [ ] 0,1 to n  , whose limit from the left exist at every point in ( ] 0,1 and whose limit from the right exists at every point in [ ) 0,1 and that are also continuous at every point that is not a dyadic rational, this is, ∆  is the space of continuous functions on [ ] 0,1 that may have jump discontinuities at dyadic rationals.Proposition 8.The space ∆  is a Banach space with the sup norm.
Proof.Let { } n x be a Cauchy sequence in ∆  , then is a Cauchy sequence in n  , therefore convergent, we denote such limit as ( ) * 0 x t , then we have that for every ε , there exists N, such that for n N ≥ : , then for every n x the limit from the right at 0 t exists, this is, such that for every 0 ε > , there exist The sequence { } n L + is defined on n  , therefore is convergent if and only if it is Cauchy.We assume by contradiction that { } n L + is not Cauchy, then there exist 1 ε such that for every 1 N , there exist 1 1 1 , n m N ≥ such that: ( ) ( ) min , δ δ δ = , then for 0 0 t t δ < − < we have: ( ) ( ) ( ) ( ) ( ) ( ) which is a contradiction, therefore the sequence { } n L + converges to a limit L + .Now let 0 ε > , let N ∈  and δ such that for Therefore L + is a right limit of * x at 0 t t = , hence the right limits exist.
Repeating the process for the left limits we conclude that the left limits of * x exist where required.

Now let
We make use of the following: Lemma 9. [20] Let the dyadic intervals ( ) be the characteristic function of the interval , then: Now we state: Theorem 10.Let the Banach space ∆  , the approximated monodromy operator k  converges in the strong operator sense to the monodromy operator U.
Proof.Let k π denote the projection mapping that takes the Walsh expansion ( ) ( ) to the Walsh approximation .
Without risk of confusion we will denote as ( ) Denote as ( ) ( ) x t Tϕ θ = the solution of ( 20) corresponding to an initial condition ( ) . Likewise denote as ( ) ( ) the solution of (21).Clearly ( ) ( ) Since the solution ( ) x t is continuous we have that the first term on the right converges to 0 as k → ∞ .We now observe that ( ) x t and ( ) x t satisfy: We have that the right side of (51) is piecewise constant, so if we define the dyadic intervals ( ) , we can write: Lemma 9, we will have for k y t is constant on the dyadic intervals then ( ) , and assume * 0 k C > .This restraint is not significant since ( ) A t and ( ) B t are bounded, therefore ( ) ( ) B t are bounded, and We have: From all of the above and defining we arrive at: Which gives: ( ) Indeed, for , and assuming ( ) we have: For ( ) Taking into account the definition of * k D we have: and From ( 61), ( 62) and (59) we arrive at: The expresion as k → ∞ , therefore we have from where it follows immediately that: which concludes the proof.
Corollary 1.If ϕ is Lipschitz, then the error of the approximation k Proof.From (59), we have that since ϕ is Lipschitz and ( ) x t is differentiable, the result follows immediately from Lemma 3.
We now prove that the approximated solution is indeed equal to the Walsh approximation of the exact solution.First we state: .Likewise let ( ) ( ) , the solution of ( 21), then ( ) ( ) Proof.We have ( ) ( ) , for some coefficients c and d , then : Since the solution ( ) x t is continuous and from (64) we have that the two terms on the right converge to zero, therefore ( ) ( ) Proof.Let 0 ε > , since X is compact, there exist finite open balls of radius Since there are finitely many i x , the desired result follows immediately.
Theorem 14 The approximated monodromy operator k  converges uniformly to the monodromy operator U.
Proof.Let B denote the closed unit ball on δ  .Since T is compact, then the image TB is also compact, and by Lemma 13 we will have that for any 0 ε > there exists k such that: The fact that With this we can establish convergence of the spectrum: λ such that any other eigenvalue of U is outside Γ , then the sum of the multiplicities of the eigenvalues of  within Γ will be equal to the multiplicity of 0 λ .Proof.Let Γ be a small circle with center at ( ) , such that any other eigenvalue of U is outside this circle.Then the spectral projection of 0 λ is given by: ( ) The spectral projection of the spectrum of k  inside Γ will be given by: ( ) this is, the algebraic multiplicity of 0 λ will be the same as the algebraic multiplicity of the spectrum of k  contained in Γ .

Approximation of the Solution of a Delayed Mathieu Equation
If one considers a generalization of a Mathieu Differential Equation to a Functional Differential Equation, we will find that we can consider different cases depending on where the parametric excitation is placed [9].We consider the scalar equation: The operator k  in (33) provides a natural way to approximate the solution of a periodic delay differential equation.Consider Equation (74), with

Stability Chart of the Delayed Mathieu Equation
We will use the approximated monodromy operator to determine the stability chart of a particular case of delayed Mathieu equation: This equation has been studied in [8]. Figure 8 shows the stability diagram for the parametric plane αβ with 0.15 δ = and 2π 16 τ = .

Conclusion
The use of Walsh functions provides the finite dimensional approximation (35) of the monodromy operator.This approximation and the analysis leading to it are virtually the same for any piecewise constant orthonormal basis which can be formed by linear combinations of Walsh functions such as Block Pulse Functions.The use of Block Pulse Functions provides a computationally  Asymptotically stable zones (grey) and unstable zones (white) are shown.Asymptotically stable zones (grey) and unstable zones (white) are shown.
inexpensive method that is useful when obtaining stability diagrams.The rate of decay of the error will be maintained regardless of the orthonormal set used [27].The use of Block Pulse Functions provides a method with a light computational load, due to the simplicity of the functions and the sparse structure of the involved matrices.Implementation of the algorithm is straightforward, it is only necessary to compute matrices (42) corresponding to the matrices ( ) ( ) B t in (1), and to substitute the remaining matrices in Section 4.1.Furthermore, the approximated monodromy operator might be used to provide insight in the nature of the PDDE, specially with a second order equation if the solution space is confined to a two dimensional vector space [28], since a similar approximation with the use of Block Pulse Functions has been used to analytically prove properties of Periodic Ordinary Differential Equations.[29].However a downside of the proposed method lies in the rate of convergence, which for certain cases is slower than the convergence of Fourier functions, and is certainly slower than the rate of convergence of approximations with, for example, Chebyshev polynomials.

,
A t B t are n n × matrices of ω-periodic functions continuous on [ ] 0,1 .Denote as  the space of conti- nuous functions from [ ]

): Theorem 1 .
If the characteristic multipliers are situated inside the unit circle { } , 1 λ λ < , then the zero solution of the system is uniformly asymptotically stable.If the multipliers of the system are inside the closed unit circle, and if multipliers situated on the unit circumference correspond to simple elementary divisors then the zero solution is uniformly stable.Remark 1.If τ ω > all the above statements are valid for r U , where rω τ ≥ and r ∈  [12].Remark 2. Note that Theorem 1 extends the properties of an undelayed case [14] [15].
only on the value of the function at that same interval.Proof.It follows immediately from Theorem 2.1.3in[17].Thus we define the projection approximation of order 2 k of an integrable function ( )

F
denote this subspace as M ′ .Likewise we can extend the domain of definition of the solution map of (1) from the space [ ] denoted ′  is isomorphic to the space [ ] on the space M corresponds to M ′ .We can say now that k , hence everywhere, and from Lemma 11. we have c d =

Theorem 15 .
The spectrum of the approximated monodromy operator k  converges to the spectrum of the monodromy operator U.More precisely, for any open set  , such that

.
Figure 3 and Figure 4, show the comparison between the solution obtained by simulation and the approximation of the solution obtained from k 

Figure 5 and
Figure 5 and Figure 6, show the magnitude of the error for the approximation of the solution obtained with the approximated operator k  , for 7 k = and 10 k = , respectively.

Figure 4 .
Figure 4. Solution of (74) for t ∈ .The solid line corresponds to the simulated solution and the dashed line corresponds to the approximated solution for 10 k = .

Figure 5 .
Figure 5. Magnitude of error of the approximated solution corresponding to 7 k = .

Figure 6 .
Figure 6.Magnitude of error of the approximate solution corresponding to 10 k = .
used to to determine the stability of periodic differential equations with respect to some parameters.Figure7shows the stability diagram of (75) for the plane αβ with 2π 32 τ = and 1.5 γ = .For this equation stability zones are disconnected and there is no symmetry with respect to the horizontal axis, contrary to the case of the undelayed Mathieu equation.Now consider the equation:
Cauchy we have that there exists 2 N such that