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

Eli A. Vazquez, Joaquin Collado

Automatic Control Department, CINVESTAV-IPN, Mexico City, Mexico.

DOI: 10.4236/am.2018.911086                

Abstract

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.

Keywords

Monodromy Operator, Periodic Delay Differential Equations, Walsh Functions, Block Pulse Functions, Finite Rank Approximation

1. 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 semidis- cretization [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 $x={\mathcal{U}}_k{x}_0$, where the directly obtained matrix ${\mathcal{U}}_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 ${\mathcal{U}}_k$ without any need of solving any equation. Convergence of ${\mathcal{U}}_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.

2. Linear Periodic Delay Differential Equations

Consider the linear PDDE:

$\stackrel{\dot{}}{x}\left(t\right)=A\left(t\right)x\left(t\right)+B\left(t\right)x\left(t-\tau \right),$ (1)

where $x\left(t\right)\in {ℝ}^n$, $t\mathrm{,}\tau \in {ℝ}_{+}$, $0<\tau \le \omega$, $A\left(t\right)\mathrm{,}B\left(t\right)$ are $n\times n$ matrices of ω-periodic functions continuous on $\left[\mathrm{0,1}\right]$. Denote as $\mathcal{C}$ the space of conti- nuous functions from $\left[-\tau \mathrm{,0}\right]$ into ${ℝ}^n$, this space is a Banach space with the norm $\Vert \phi \Vert =\underset{-\tau \le \theta \le 0}{max}\Vert \phi \left(\theta \right)\Vert$ [11] .

A solution of (1) with initial condition $\phi \in \mathcal{C}$ at $t_0$ is understood as a mapping $\xi \mathrm{<mo>\; :\; </mo>}\left[-\tau \mathrm{,0}\right]\times ℝ\times C\to {ℝ}^n$, such that $\xi \left(\theta ,t_0,\phi \right)=\phi \left(\theta \right)$ for $-\tau \le \theta \le 0$ and that $\xi \left(\theta \mathrm{,}t\mathrm{,}\phi \right)=x\left(t+\theta \right)$ for $-\tau \le \theta \le 0$ and $t+\theta \ge t_0$, where $x\left(t\right)$ sa- tisfies (1) for $t\ge t_0$ with $x\left(\theta \right)=\phi \left(\theta \right)$ for $\theta \in \left[-\tau \mathrm{,0}\right]$. At time t a solution will be an element of space C. Explicitly we will have:

$\xi \left(\theta \mathrm{,}t\mathrm{,}\xi \left(\theta \mathrm{,}{t}_0\mathrm{,}\phi \right)\right)=T\left(t_0\mathrm{,}{t}_0+t\right)\xi \left(\theta \mathrm{,}{t}_0\mathrm{,}\phi \right)\mathrm{,}$ (2)

with $\theta \in \left[-\tau \mathrm{,0}\right]$. The operator T is called the solution mapping and is analogous to the state transition matrix of the undelayed case.

2.1. Monodromy Operator

Taking into account the periodic nature of (1) it is relevant to consider (2) when $t=t_0+\omega$, in this case we will denote the solution mapping as $U\left(t_0\right)\triangleq T\left(t_0\mathrm{,}{t}_0+\omega \right)$.

Next some properties consequence of $U(\cdot )$ being completely continuous are enlisted, [12] :

• The spectrum $\sigma \left(U\left(t_0\right)\right)$ is a countable compact set in the complex plane.

• $0\in \sigma \left(U\left(t_0\right)\right)$, and if $\lambda \in \sigma \left(U\left(t_0\right)\right)$ then $\left(U\phi \right)\left(t_0\right)=\lambda \phi$ has a nonzero solution $\phi \in C$.

• If the cardinality of the set $\sigma \left(U\left(t_0\right)\right)$ is infinite then the only limit point of $\sigma \left(U\left(t_0\right)\right)$ is 0.

• The cardinality of $\sigma \left(U\left(t_0\right)\right)$ outside any disc $R_r$ of radius r, $R_r=\left\{z\in ℂ,\left|z\right|\le r\right\}$ is finite.

• $\sigma \left(U\left(t_0\right)\right)=\sigma \left(U\left(t_1\right)\right)$ $\forall t_0\mathrm{,}{t}_1\in ℝ$.

$U\left(0\right)$ henceforth denoted simply as U will be called the monodromy opera- tor of (1) and the elements $\lambda \in \sigma \left(U\right)$, $\lambda \ne 0$ will be called characteristic multipliers.

2.2. 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):

Theorem 1. If the characteristic multipliers are situated inside the unit circle $\left\{\lambda ,\Vert \lambda \Vert <1\right\}$, 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 $\tau >\omega$ all the above statements are valid for $U^r$, where $r\omega \ge \tau$ and $r\in \mathbb{Z}$ [12] .

Remark 2. Note that Theorem 1 extends the properties of an undelayed case [14] [15] .

3. Walsh Functions

3.1. Definition

Walsh functions are a set of piecewise constant complete orthonormal functions introduced in [16] and defined on the interval $\left[\mathrm{0,1}\right)$, 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 l-th Walsh function may be defined as [17] :

$w_{\mathcal{l}}\left(t\right)={\displaystyle \underset{i=0}{\overset{k}{\prod }}}{\left[\text{sign}\left(\sin 2^{i-1}\text{π}t\right)\right]}^{{\epsilon }_i},$ (3)

where the ${\epsilon }_i$ are the coefficients of the unique binary expansion of $\mathcal{l}$, namely

$\mathcal{l}={\displaystyle \underset{i=0}{\overset{k}{\sum }}}{\epsilon }_i{2}^i\mathrm{,}\text{with}{\epsilon }_i\in \left\{\mathrm{0,1}\right\}$ (4)

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\in L_{\left[\mathrm{0,1}\right]}^2$ will correspond the Walsh approximation:

$f\left(t\right)~{\displaystyle \underset{i=0}{\overset{n}{\sum }}}{h}_i{w}_i\left(t\right),$ (5)

where $h_i={\displaystyle {\int }_0^1}f\left(t\right)w_i\left(t\right)\text{d}t.$ The right side of (5) will converge in norm to f for any $f\in L_{\left[\mathrm{0,1}\right]}^2$ and if $n=2^k$, it will converge uniformly for any continuous function that has, at most, jump discontinuities at dyadic rationals (numbers of the form $\frac{a}{2^b},a\in \mathbb{N},b\in \mathbb{N}$ ) [16] . When $n=\infty$, (5) will be called a Walsh expan- sion.

3.2. Properties of Walsh Functions

We will define a Walsh matrix $W\left[k\right]$ as the $2^k\times 2^k$ matrix consisting on the discretization values of the Walsh functions over the interval $\left[\mathrm{0,1}\right]$.

For example we will have for $k=3$ :

$W\left[4\right]=\left[\begin{array}{rrrr}\hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill -1& \hfill -1\\ \hfill 1& \hfill -1& \hfill 1& \hfill -1\\ \hfill 1& \hfill -1& \hfill -1& \hfill 1\end{array}\right].$

We can also define the Delayed Walsh Matrix $W\left[-m\mathrm{,}k\right]$ as a Walsh matrix of order $2^k$ shifted m columns to the right and with the first m columns being zeros. In the same manner we define a Forwarded Walsh Matrix $W\left[+m\mathrm{,}k\right]$, but shifted to the left:

Figure 1. Walsh functions in dyadic ordering.

$W\left[-2,4\right]=\left[\begin{array}{rrrr}\hfill 0& \hfill 0& \hfill 1& \hfill 1\\ \hfill 0& \hfill 0& \hfill 1& \hfill 1\\ \hfill 0& \hfill 0& \hfill 1& \hfill -1\\ \hfill 0& \hfill 0& \hfill 1& \hfill -1\end{array}\right];W\left[+2,4\right]=\left[\begin{array}{rrrr}\hfill 1& \hfill 1& \hfill 0& \hfill 0\\ \hfill -1& \hfill -1& \hfill 0& \hfill 0\\ \hfill 1& \hfill -1& \hfill 0& \hfill 0\\ \hfill -1& \hfill 1& \hfill 0& \hfill 0\end{array}\right]$ (6)

We define the vector consisting on the first $2^k$ Walsh functions as:

${\bar{w}}_k\left(t\right)={\left[w_0\left(t\right)w_1\left(t\right)\cdots w_{2^k-1}\left(t\right)\right]}^{\text{T}}.$ (7)

Since the order of the approximation should be clear from context we write simply $\bar{w}\left(t\right)$. We define as well the dyadic sum $\oplus$ between tho nonnegative integers as:

$q\oplus r={\displaystyle \underset{i=0}{\overset{k}{\sum }}}\left|q_i-r_i\right|2^i,$ (8)

where $q_i$ and $r_i$ are respectively the coefficients of the binary expansions of q and r as in (4). To a vector $\sigma ={\left[{\sigma }_0{\sigma }_2\cdots {\sigma }_{2^k-1}\right]}^{\text{T}}$ we associate a $2^k\times 2^k$ sym- metric matrix $\Lambda \left(\sigma \right)$ whose i, jth entry will be given by ${\sigma }_{i-1\oplus j-1}$. If for example $\sigma ={\left[{\sigma }_0{\sigma }_1{\sigma }_2{\sigma }_3\right]}^{\text{T}}$, then:

$\Lambda \left(\sigma \right)=\left[\begin{array}{cccc}{\sigma }_0& {\sigma }_1& {\sigma }_2& {\sigma }_3\\ {\sigma }_1& {\sigma }_0& {\sigma }_3& {\sigma }_2\\ {\sigma }_2& {\sigma }_3& {\sigma }_0& {\sigma }_1\\ {\sigma }_3& {\sigma }_2& {\sigma }_1& {\sigma }_0\end{array}\right].$

We make use of the following Lemma:

Lemma 2. [19] Let $\sigma \in {ℝ}^{2^k}$, then:

$\bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right)\sigma =\Lambda \left(\sigma \right)\bar{w}\left(t\right)\mathrm{.}$ (9)

The integral of a Walsh vector can be approximated by a Walsh expansion of the form:

${\displaystyle {\int }_0^t}\bar{w}\left(s\right)\text{d}s=P\bar{w}\left(t\right)\mathrm{,}$ (10)

where P is given by [20] :

(11)

The approximated integral converges to the exact integral uniformly [21] .

To a function matrix $A\left(t\right)$ we associate the Walsh approximation $\hat{A}\left(t\right)$ :

$\hat{A}\left(t\right)=\alpha \Omega \left(t\right),$ (12)

where

$\alpha =\left[\begin{array}{cccc}{\bar{\alpha }}_{11}^{\text{T}}& {\bar{\alpha }}_{12}^{\text{T}}& \cdots & {\bar{\alpha }}_{1n}^{\text{T}}\\ {\bar{\alpha }}_{21}^{\text{T}}& {\bar{\alpha }}_{22}^{\text{T}}& \cdots & {\bar{\alpha }}_{2n}^{\text{T}}\\ ⋮& ⋮& \ddots & ⋮\\ {\bar{\alpha }}_{n1}^{\text{T}}& {\bar{\alpha }}_{n2}^{\text{T}}& \cdots & {\bar{\alpha }}_{nn}^{\text{T}}\end{array}\right],\Omega \left(t\right)=\left[\begin{array}{cccc}\bar{w}\left(t\right)& 0& \cdots & 0\\ 0& \bar{w}\left(t\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \bar{w}\left(t\right)\end{array}\right],$ (13)

and each ${\bar{\alpha }}_{i\mathrm{,}j}$ is the vector of coefficients of the Walsh approximation of each entry of $A\left(t\right)$.

Let $A\left(t\right)$ be a $n\times n$ function matrix with Walsh approximation $\hat{A}\left(t\right)$ and $f\left(t\right)$ a function vector with Walsh approximation $\hat{f}\left(t\right)$, whose elements are of the form $f_i\left(t\right)=h_i^{\text{T}}\bar{w}\left(t\right)$, $i=1,\cdots ,n$ where $h_i\in {ℝ}^{2^k}$ is the vector of Walsh coefficients of the approximation of each element of f. Define $H={\left[\begin{array}{ccc}{h}_1^{\text{T}}& \cdots & h_n^{\text{T}}\end{array}\right]}^{\text{T}}$, then from (9) we have:

$\begin{array}{c}\hat{A}\left(t\right)\hat{f}\left(t\right)=\alpha \Omega \left(t\right)\left[\begin{array}{c}{\bar{w}}^{\text{T}}\left(t\right)h_1\\ {\bar{w}}^{\text{T}}\left(t\right)h_2\\ ⋮\\ {\bar{w}}^{\text{T}}\left(t\right)h_n\end{array}\right]=\alpha \left[\begin{array}{c}\bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right)h_1\\ \bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right)h_2\\ ⋮\\ \bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right)h_n\end{array}\right]\\ =\left[\begin{array}{c}{\displaystyle \underset{r=1}{\overset{n}{\sum }}}{\left(h_r^{\text{T}}\bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right){\bar{\alpha }}_{1,r}\right)}^{\text{T}}\\ {\displaystyle \underset{r=1}{\overset{n}{\sum }}}{\left(h_r^{\text{T}}\bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right){\bar{\alpha }}_{2,r}\right)}^{\text{T}}\\ ⋮\\ {\displaystyle \underset{r=1}{\overset{n}{\sum }}}{\left(h_r^{\text{T}}\bar{w}\left(t\right){\bar{w}}^{\text{T}}\left(t\right){\bar{\alpha }}_{n,r}\right)}^{\text{T}}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \underset{r=1}{\overset{n}{\sum }}}{\bar{w}}^{\text{T}}\left(t\right)\Lambda \left({\bar{\alpha }}_{1,r}\right)h_r\\ {\displaystyle \underset{r=1}{\overset{n}{\sum }}}{\bar{w}}^{\text{T}}\left(t\right)\Lambda \left({\bar{\alpha }}_{2,r}\right)h_r\\ ⋮\\ {\displaystyle \underset{r=1}{\overset{n}{\sum }}}{\bar{w}}^{\text{T}}\left(t\right)\Lambda \left({\bar{\alpha }}_{n,r}\right)h_r\end{array}\right]\\ ={\Omega }^{\text{T}}\left(t\right)\bar{\Lambda }\left(\alpha \right)H,\end{array}$ (14)

where:

$\bar{\Lambda }\left(\alpha \right)=\left[\begin{array}{cccc}\Lambda \left({\bar{\alpha }}_{11}\right)& \Lambda \left({\bar{\alpha }}_{12}\right)& \cdots & \Lambda \left({\bar{\alpha }}_{1n}\right)\\ \Lambda \left({\bar{\alpha }}_{21}\right)& \Lambda \left({\bar{\alpha }}_{22}\right)& \cdots & \Lambda \left({\bar{\alpha }}_{2n}\right)\\ ⋮& ⋮& \ddots & ⋮\\ \Lambda \left({\bar{\alpha }}_{n1}\right)& \Lambda \left({\bar{\alpha }}_{n2}\right)& \cdots & \Lambda \left({\bar{\alpha }}_{nn}\right)\end{array}\right].$ (15)

3.3. Approximation Error

Regarding the error of Walsh approximation we have:

Lemma 3. [17] If f satisfies the Lipschitz condition, then:

$E_k={\Vert {\displaystyle \underset{i=0}{\overset{p=2^k}{\sum }}}{c}_i{w}_i\left(t\right)-f\left(t\right)\Vert }_{\infty }\le C{2}^{-k},$ (16)

for some constant C.

3.4. Block Pulse Functions

The set of Block pulse functions of order p is defined on the interval $\left[\mathrm{0,1}\right)$ as the set $\left\{{\psi }_0\left(t\right)\mathrm{,}\cdots \mathrm{,}{\psi }_{m-1}\left(t\right)\right\}$, where [22] :

${\psi }_i\left(t\right)=\{\begin{array}{l}\mathrm{1,}t\in \left[\frac{i}{p}\mathrm{,}\frac{i+1}{p}\right)\\ \mathrm{0,}\text{otherwise}\end{array}i=\mathrm{1,}\cdots \mathrm{,}p-1.$ (17)

Block Pulse Functions, shown on Figure 2, are orthogonal, easily normalizable [23] and when $m\to \infty$ they form a complete set [22] . We restrict ourselves to

Figure 2. Block pulse functions.

the case $p=2^k$ 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 $P_{\mathrm{<mo>\; *\; </mo>}}\triangleq W\left[k\right]PW\left[k\right]=\frac{1}{2}{I}_{2^k}+Q+Q^2+\cdots +Q^{2^k-1}$, where Q is a nilpotent matrix

with ones above the diagonal and zeros everywhere else:

$P_{*}=\left[\begin{array}{cccc}\frac{1}{2}& 1& \cdots & 1\\ 0& \frac{1}{2}& \cdots & 1\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{1}{2}\end{array}\right].$ (18)

Lemma 5. [25] Let $\hat{f}\left(t\right)={\sigma }^{\text{T}}\bar{w}\left(t\right)$ be the Walsh approximation of the function $f\left(t\right)$, with $\sigma \in {ℝ}^{2^k}$ being the vector of coefficients of the Walsh approximation. Then $\Lambda \left(\sigma \right)$ is similar to a diagonal matrix:

$\begin{array}{c}{\Lambda }_{\mathrm{<mo>\; *\; </mo>}}\left(\sigma \right)=W^{-1}\left[k\right]\Lambda \left(\sigma \right)W\left[k\right]\\ =\left[\begin{array}{cccc}\hat{f}\left(\frac{0}{2^k}\right)& 0& \cdots & 0\\ 0& \hat{f}\left(\frac{1}{2^k}\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \hat{f}\left(\frac{2^k-1}{2^k}\right)\end{array}\right]\mathrm{.}\end{array}$ (19)

4. Approximation of the Monodromy Operator

Without loss of generality we assume $\omega =1$. We approximate the monodromy operator by projecting (1) on a finite dimensional subspace of $L^2\left[\mathrm{0,1}\right)$ 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\tau =\omega$ for some $m\in \mathbb{N}$.

Integrating (1) from 0 to t we will have:

$x\left(t\right)-x\left(0\right)={\displaystyle {\int }_0^t}A\left(s\right)x\left(s\right)\text{d}s+{\displaystyle {\int }_0^t}B\left(s\right)x\left(s-\tau \right)\text{d}s,$ (20)

with $x\left(\theta \right)=\phi \left(\theta \right)$ for $\theta \in \left[-\tau \mathrm{,0}\right]$. A solution of (20) will correspond to a solution of (1) [26] . Let ${\pi }_k$ denote the projection mapping that takes the Walsh expansion $f\left(t\right)={\displaystyle {\sum }_{\mathcal{l}=0}^{\infty }}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$ to the Walsh approximation

$\hat{f}\left(t\right)={\displaystyle {\sum }_{\mathcal{l}=0}^{2^k}}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$. Along with (20) we introduce its projection onto the space

$M=span{\left\{w_0,w_1,\cdots ,w_{2^k-1}\right\}}^n$ :

$\hat{x}\left(t\right)-x\left(0\right)={\pi }_k{\displaystyle {\int }_0^t}\hat{A}\left(s\right)\hat{x}\left(s\right)\text{d}s+{\pi }_q{\displaystyle {\int }_0^t}\hat{B}\left(s\right)\hat{x}\left(s-\tau \right)\text{d}s,$ (21)

where $\hat{x}\left(t\right)\in M$ and $\phi \left(0\right)\in M$ since it is a constant. $\hat{A}\left(t\right)$ and $\hat{B}\left(t\right)$ correspond to ${\pi }_{k}A\left(t\right)$ and ${\pi }_{k}B\left(t\right)$, this is, the approximations of $A\left(t\right)$ and $B\left(t\right)$, respectively, as in (12). $\hat{x}\left(x\right)\left(\theta \right)$ is still not defined for $\theta \in \left[-\tau \mathrm{,0}\right)$, we 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 $\left[\frac{i}{2^k}\mathrm{,}\frac{i+1}{2^k}\right)$, $i=0,\cdots ,2^k-1$, depends only on the value of the function at that same interval.

Proof. It follows immediately from Theorem 2.1.3 in [17] .

Thus we define the projection ${\pi }_k\phi \left(\theta \right)={\displaystyle {\sum }_{\mathcal{l}=0}^{2^k-1}}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(\theta +\omega \right)\triangleq \hat{\phi }\left(\theta \right)$ as the Walsh approximation of order $2^k$ of an integrable function ${\phi }^{\mathrm{<mo>\; *\; </mo>}}\left(\theta \right)$ defined on $\left[-\omega \mathrm{,0}\right)$, that is equal to $\phi \left(\theta \right)$ in $\left[-\tau \mathrm{,0}\right)$ and 0 everywhere else, and we make $\hat{x}\left(\theta \right)=\hat{\phi }\left(\theta \right)$ for $\theta \in \left[-\tau \mathrm{,0}\right]$. Since Walsh functions were not defined at $t=1$

we simply set $\hat{\phi }\left(0\right)=\hat{\phi }\left(-\frac{1}{2^k}\right)$.

We split the second integral in (20) as:

$\begin{array}{l}{\displaystyle {\int }_0^t}B\left(s\right)x\left(s-\tau \right)\text{d}s\\ ={\displaystyle {\int }_0^{\tau }}B\left(s\right)x\left(s-\tau \right)\text{d}s+{\displaystyle {\int }_{\tau }^t}B\left(s\right)x\left(s-\tau \right)\text{d}s\\ ={\displaystyle {\int }_0^{\tau }}B\left(s\right)\phi \left(s-\tau \right)\text{d}s+{\displaystyle {\int }_{\tau }^t}B\left(s\right)x\left(s-\tau \right)\text{d}s\\ ={\displaystyle {\int }_0^t}\left(1-S\left(s-\tau \right)\right)B\left(s\right)\phi \left(s-\tau \right)\text{d}s+{\displaystyle {\int }_0^t}S\left(s-\tau \right)B\left(s\right)x\left(s-\tau \right)\text{d}s,\end{array}$ (22)

where $S\left(t\right)$ is the unit step function.

We make use of the following Lemma:

Lemma 7. [17] Any function $f\left(t\right)$ constant on the intervals of the form $\left[\frac{i}{2^k}\mathrm{,}\frac{i+1}{2^k}\right)$, $0\le i\le 2^k-1$ can be represented in the form:

$f\left(t\right)={\displaystyle \underset{\mathcal{l}=0}{\overset{2^k-1}{\sum }}}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right),$ (23)

i.e. the non zero coefficients of the Walsh expansion of $f\left(t\right)$ have indices no greater than $2^k-1$. Moreover, this representation is unique.

From linearity of ${\pi }_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:

$\begin{array}{l}{\pi }_k{\displaystyle {\int }_0^t}\hat{B}\left(s\right)\hat{x}\left(s-\tau \right)\text{d}s\\ ={\pi }_k{\displaystyle {\int }_0^t}\left(1-S\left(s-\tau \right)\right)\hat{B}\left(s\right)\hat{\phi }\left(s-\tau \right)\text{d}s+{\pi }_k{\displaystyle {\int }_0^t}S\left(s-\tau \right)\hat{B}\left(s\right)\hat{x}\left(s-\tau \right)\text{d}s,\end{array}$ (24)

Since $\hat{x}\left(t\right)\in M$ we have that $\hat{x}\left(t\right)={\left[h_1^{\text{T}}\bar{w}\left(t\right),\cdots ,h_n\bar{w}\left(t\right)\right]}^{\text{T}}$, we also have that $\hat{\phi }\left(\theta \right)={\left[c_1^{\text{T}}\bar{w}\left(\theta +\omega \right),\cdots ,c_n\bar{w}\left(\theta +\omega \right)\right]}^{\text{T}}$ with each $h_i$ and $c_i$ being vectors of Walsh coefficients. Defining $H={\left[h_1^{\text{T}},\cdots ,h_n^{\text{T}}\right]}^{\text{T}}$, $\Phi ={\left[c_1^{\text{T}}\mathrm{,}\cdots \mathrm{,}{c}_n^{\text{T}}\right]}^{\text{T}}$, and recalling (7) and (12), we can write (21) as:

$\begin{array}{l}{\Omega }^{\text{T}}\left(t\right)H-\hat{x}\left(0\right)\\ ={\pi }_q{\displaystyle {\int }_0^t}\alpha \Omega \left(s\right){\Omega }^{\text{T}}\left(s\right)H\text{d}s+{\pi }_k{\displaystyle {\int }_0^t}\beta \Omega \left(s\right)S\left(s-\tau \right){\Omega }^{\text{T}}\left(s-\tau \right)H\text{d}s\\ +{\pi }_k{\displaystyle {\int }_0^t}\beta \Omega \left(s\right)\left(1-S\left(s-\tau \right)\right){\Omega }^{\text{T}}\left(s-\tau +\omega \right)\Phi \text{d}s,\end{array}$ (25)

From Lemma 7. we have that both, $S\left(t-\tau \right){\Omega }^{\text{T}}\left(t-\tau \right)$ and $\left(1-S\left(t-\tau \right)\right){\Omega }^{\text{T}}\left(s-\tau +\omega \right)$, can be uniquely represented in terms of $\Omega \left(t\right)$ :

$S\left(t-\tau \right){\Omega }^{\text{T}}\left(t-\tau \right)={\Omega }^{\text{T}}\left(t\right){\bar{W}}_D,$

$\left(1-S\left(t-\tau \right)\right){\Omega }^{\text{T}}\left(s-\tau +\omega \right)={\Omega }^{\text{T}}\left(t\right){\bar{W}}_F,$ (26)

${\bar{W}}_F$ and ${\bar{W}}_D$ are $2^{k}n\times 2^{k}n$ block diagonal matrices, whose diagonal entries are given respectively by $W_F^{\text{T}}$ and $W_D^{\text{T}}$, which satisfy $W_F\bar{w}\left(t\right)=\left(1-S\left(t-\tau \right)\right)w\left(s-\tau +\omega \right)$ and $W_D\bar{w}\left(t\right)=S\left(t-\tau \right)\bar{w}\left(t-\tau \right)$ When ap- proximating by Walsh functions, matrix $W_D$ is given by $W_D=\frac{1}{2^k}\left(W\left[-m,k\right]W\left[k\right]\right)$ and it is called the Walsh Shift Operator [10] . $W_F$ is given in an analogous way as $W_F=\frac{1}{2^k}\left(W\left[+\left(1-m\right),k\right]W\left[k\right]\right)$. The term $\hat{x}\left(0\right)$ will also have an unique representation:

$\hat{x}\left(0\right)={\Omega }^{\text{T}}\left(t\right){\bar{W}}_C\Phi \mathrm{,}$ (27)

since $\hat{x}\left(0\right)=\hat{\phi }\left(0\right)$. ${\bar{W}}_C$ will be again a $2^{k}n\times 2^{k}n$ block diagonal matrix, whose diagonal entries are given by $W_C^{\text{T}}$. In the case of Walsh functions, matrix $W_C$ evaluates $\varphi \left(t\right)$ at time $2^k-1$ and assigns its value to the coefficient of the constant Walsh function $w_0\left(t\right)$, hence $W_C^{\text{T}}={\left[\bar{w}\left(2^k-1\right),0,\cdots ,0\right]}^{\text{T}}$.

From (14), (27) and (26) we have that we can write (25) as:

$\begin{array}{l}{\Omega }^{\text{T}}\left(t\right)H-{\Omega }^{\text{T}}\left(t\right){\bar{W}}_C\Phi \\ ={\pi }_k{\displaystyle {\int }_0^t}{\Omega }^{\text{T}}\left(s\right)\text{d}s\bar{\Lambda }\left(\alpha \right)H+{\pi }_k{\displaystyle {\int }_0^t}{\Omega }^{\text{T}}\left(s\right)\text{d}s\bar{\Lambda }\left(\beta \right){\bar{W}}_F\Phi \\ +{\pi }_k{\displaystyle {\int }_0^t}{\Omega }^{\text{T}}\left(s\right)\text{d}s\bar{\Lambda }\left(\beta \right){\bar{W}}_{D}H.\end{array}$ (28)

From (10) we have that:

${\pi }_k{\displaystyle {\int }_0^t}{\Omega }^{\text{T}}\left(s\right)\text{d}s={\Omega }^{\text{T}}\left(t\right)\bar{P},$ (29)

where $\bar{P}$ is a $2^{k}n\times 2^{k}n$ matrix given by:

$\bar{P}=\left[\begin{array}{cccc}{P}^{\text{T}}& 0& \cdots & 0\\ 0& P^{\text{T}}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & P^{\text{T}}\end{array}\right],$ (30)

from which we arrive at:

$\begin{array}{l}{\Omega }^{\text{T}}\left(t\right)H-{\Omega }^{\text{T}}\left(t\right){\bar{W}}_C\Phi \\ ={\Omega }^{\text{T}}\left(t\right)\bar{P}\bar{\Lambda }\left(\alpha \right)H{\Omega }^{\text{T}}\left(t\right)+\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_F\Phi +{\Omega }^{\text{T}}\left(t\right)\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_{D}H.\end{array}$ (31)

Since $\Omega \left(t\right)$ is nonsingular we have:

$H={\left[I-\bar{P}\bar{\Lambda }\left(\alpha \right)-\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_D\right]}^{-1}\left[{\bar{W}}_C+\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_F\right]\Phi .$ (32)

From (32) we can define a mapping ${\mathcal{T}}_k\mathrm{<mo>\; :\; </mo>}M\to M$ :

${\mathcal{T}}_k={\left[I-\bar{P}\bar{\Lambda }\left(\alpha \right)-\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_D\right]}^{-1}\left[{\bar{W}}_C+\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_F\right]$ (33)

However by construction we can see that the domain of ${\mathcal{T}}_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 $t\in \left[\mathrm{0,}\omega -\tau \right)$, we denote this subspace as $M^{\prime }$. Likewise we can extend the domain of definition of the solution map of (1) from the space ${\mathcal{C}}_{\left[-\tau \mathrm{,0}\right]}$, to the subspace of $L_{\left[-\omega \mathrm{,0}\right]}^2$ consisting of al functions that are continuous on $\left[-\tau \mathrm{,0}\right]$ that are equal to 0 for $t\in \left[-\omega \mathrm{,}\tau \right)$, this space, denoted ${\mathcal{C}}^{\prime }$ is isomorphic to the space ${\mathcal{C}}_{\left[-\tau \mathrm{,0}\right]}$ and its projection on the space M corresponds to $M^{\prime }$. We can say now that ${\mathcal{T}}_k$ is an approximation of the solution mapping T of (1), we obtain an approximated solution $\hat{x}\left(t\right)={\Omega }^{\text{T}}\left(t\right)H$ 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=\omega$, this is, we must know the solution of (1) for $t\in \left[\omega -\tau \mathrm{,}\omega \right]$ corresponding to an initial condition. Since $M^{\prime }$ is the projection in the space M of ${\mathcal{C}}^{\prime }$ and since ${\mathcal{C}}^{\prime }$ is isomorphic to ${\mathcal{C}}_{\left[\omega -\tau \mathrm{,}\omega \right]}$, the approximation of the state at $t=\omega$ will be given by the projection of the approximated solution into $M^{\prime }$. Hence, by Lemma 7 we will have:

$H^{\prime }={\bar{W}}_{P}H.$ (34)

where ${\bar{W}}_p$ is a $2^{k}n\times 2^{k}n$ is a block diagonal matrix which projects the approximated solution ${\Omega }^{\text{T}}\left(t\right)H$ into the subspace $M^{\prime }$, with diagonal entries $W_P^{\text{T}}$. In the case of Walsh functions, we will have

$W_P^{\text{T}}=\frac{1}{2^k}{W}_{S\left(t-\left(1-m\right)\right)}\left[k\right]W\left[k\right]$

where $W_{S\left(t-\left(1-m\right)\right)}\left[k\right]$ 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:

${\mathcal{U}}_k={\bar{W}}_P{\left[I-\bar{P}\bar{\Lambda }\left(\alpha \right)-\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_D\right]}^{-1}\left[{\bar{W}}_C+\bar{P}\bar{\Lambda }\left(\beta \right){\bar{W}}_F\right]\mathrm{.}$ (35)

Approximation by Block Pulse Functions

Equation (35) for the approximation of ${\mathcal{U}}_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:

${\bar{W}}_P=diag\left\{W_P^{\text{T}}\mathrm{,}\cdots \mathrm{,}{W}_P^{\text{T}}\right\}\mathrm{,}$ (36)

$W_P^{\text{T}}=\left[\begin{array}{cc}{0}_{2^k-m\times 2^k-m}& 0_{2^k-m\times m}\\ 0_{m\times 2^k-m}& I_{m\times m}\end{array}\right],$

${\bar{W}}_C=diag\left\{W_C^{\text{T}}\mathrm{,}\cdots \mathrm{,}{W}_C^{\text{T}}\right\}\mathrm{,}$ (37)

$W_C^{\text{T}}=\left[\begin{array}{ccccc}0& 0& \cdots & 0& 1\\ 0& 0& \cdots & 0& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 1\end{array}\right],$

$\bar{P}=diag\left\{P^{\text{T}}\mathrm{,}\cdots \mathrm{,}{P}^{\text{T}}\right\}$ (38)

$P=\frac{1}{2^k}\left[\begin{array}{cccc}\frac{1}{2}& 1& \cdots & 1\\ 0& \frac{1}{2}& \cdots & 1\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{1}{2}\end{array}\right].$

If we define:

$\begin{array}{cc}{Q}_{-}=\left[\begin{array}{ccccc}0& \cdots & 0& 0& 0\\ 1& \cdots & 0& 0& 0\\ ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& \cdots & 1& 0& 0\\ 0& \cdots & 0& 1& 0\end{array}\right],& Q_{+}=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 0\\ 0& 0& 0& \cdots & 0\end{array}\right],\end{array}$ (39)

we will have:

${\bar{W}}_D=diag\left\{W_D^{\text{T}},\cdots ,W_D^{\text{T}}\right\}$ (40)

$W_D^{\text{T}}=Q_{-}^m,$

and

${\bar{W}}_F=diag\left\{W_F^{\text{T}},\cdots ,W_F^{\text{T}}\right\}$ (41)

$W_F^{\text{T}}=Q_{+}^{2^k-m}.$

Finally, if in (15) we have $A\left(t\right)=\left[\begin{array}{cccc}{a}_{11}\left(t\right)& a_{12}\left(t\right)& \cdots & a_{1n}\left(t\right)\\ a_{21}\left(t\right)& a_{22}\left(t\right)& \cdots & a_{2n}\left(t\right)\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}\left(t\right)& a_{n2}\left(t\right)& \cdots & a_{nn}\left(t\right)\end{array}\right]$, then $\bar{\Lambda }\left(\alpha \right)$ will be given as in (15) with:

$\Lambda \left({\bar{\alpha }}_{ij}\right)=\left[\begin{array}{cccc}{\hat{a}}_{ij}\left(\frac{0}{2^k}\right)& 0& \cdots & 0\\ 0& {\hat{a}}_{ij}\left(\frac{1}{2^k}\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\hat{a}}_{ij}\left(\frac{2^k-1}{2^k}\right)\end{array}\right],$ (42)

where

${\hat{a}}_{ij}\left(\frac{i}{2^k}\right)={\displaystyle {\int }_{\frac{i}{2^k}}^{\frac{i+1}{2^k}}}{a}_{ij}\left(t\right)\text{d}t,i=0,1,\cdots ,2^k-1,$ (43)

with the same for $\bar{\Lambda }\left(\beta \right)$ and $B\left(t\right)$.

5. Convergence of ${\mathcal{U}}_k$

We denote as ${\mathcal{C}}_{\Delta }$ the space of functions from $\left[\mathrm{0,1}\right]$ to ${ℝ}^n$, whose limit from the left exist at every point in $\left(\mathrm{0,1}\right]$ and whose limit from the right exists at every point in $\left[\mathrm{0,1}\right)$ and that are also continuous at every point that is not a dyadic rational, this is, ${\mathcal{C}}_{\Delta }$ is the space of continuous functions on $\left[\mathrm{0,1}\right]$ that may have jump discontinuities at dyadic rationals.

Proposition 8. The space ${\mathcal{C}}_{\Delta }$ is a Banach space with the sup norm.

Proof. Let $\left\{x_n\right\}$ be a Cauchy sequence in ${\mathcal{C}}_{\Delta }$, then $\forall \epsilon >0$, $\exists N$ such that $n\mathrm{,}m\ge N$ implies:

$\underset{t\in \left[0,1\right]}{\sup }\Vert x_n\left(t\right)-x_m\left(t\right)\Vert <\epsilon .$

Therefore for every $t_0\in \left[\mathrm{0,1}\right]$ we have $\Vert x_n\left(t_0\right)-x_m\left(t_0\right)\Vert <\epsilon$, for $n\mathrm{,}m\ge N$. Thus $\left\{x_n\left(t_0\right)\right\}$ is a Cauchy sequence in ${ℝ}^n$, therefore convergent, we denote such limit as $x^{\mathrm{<mo>\; *\; </mo>}}\left(t_0\right)$, then we have that for every $\epsilon$, there exists N, such that for $n\ge N$ :

$\Vert x_n\left(x_0\right)-x^{\mathrm{<mo>\; *\; </mo>}}\left(x_0\right)\Vert \le \epsilon \mathrm{,}\forall x\in \left[\mathrm{0,1}\right]\mathrm{,}$

this is $\left\{x_n\right\}$ converges to $x^{\mathrm{<mo>\; *\; </mo>}}$ uniformly.

Now let $t_0\in \left[\mathrm{0,1}\right)$, then for every $x_n$ the limit from the right at $t_0$ exists, this is, $\forall x_n\in \left\{x_n\right\}$, $\exists !L_n^{+}$ such that for every $\epsilon >0$, there exist ${\delta }_n>0$ such that:

$0<t-t_0<\delta \Rightarrow \Vert x_n\left(t\right)-L_n^{+}\Vert <\epsilon .$ (44)

The sequence $\left\{L_n^{+}\right\}$ is defined on ${ℝ}^n$, therefore is convergent if and only if it is Cauchy. We assume by contradiction that $\left\{L_n^{+}\right\}$ is not Cauchy, then there exist ${\epsilon }_1$ such that for every $N_1$, there exist $n_1\mathrm{,}{m}_1\ge N_1$ such that:

$\Vert L_{m_1}-L_{n_1}\Vert \ge {\epsilon }_1\mathrm{.}$ (45)

Since $\left\{x_n\right\}$ is Cauchy we have that there exists $N_2$ such that $n\mathrm{,}m\ge N_2$ implies:

$\underset{t\in \left[0,1\right]}{\sup }\Vert x_n\left(t\right)-x_n\left(t\right)\Vert <\frac{{\epsilon }_1}{3}.$

Let $n_1\mathrm{,}{m}_1\ge N_2$ such that (45) holds. Let ${\delta }_1$ such that for $0<t-t_0<{\delta }_1$ we have:

$\Vert x_{n_1}\left(t\right)-L_{n_1}\Vert <\frac{{\epsilon }_1}{3},$

similarly let ${\delta }_2$ such that $0<t-t_0<{\delta }_2$ we have:

$\Vert x_{m_1}\left(t\right)-L_{m_1}\Vert <\frac{{\epsilon }_1}{3}.$

Take $\delta =\min \left\{{\delta }_1,{\delta }_2\right\}$, then for $0<t-t_0<\delta$ we have:

$\Vert x_{n_1}\left(t\right)-x_{m_1}\left(t\right)\Vert <\frac{{\epsilon }_1}{3}$

$\Vert \left(x_{n_1}\left(t\right)-L_{n_1}^{+}+x_{m_1}\left(t\right)-L_{m_1}^{+}\right)-\left(L_{m_1}^{+}-L_{n_1}^{+}\right)\Vert <\frac{{\epsilon }_1}{3}$

$\left|\Vert x_{n_1}\left(t\right)-L_{n_1}^{+}+x_{m_1}\left(t\right)-L_{m_1}^{+}\Vert -\Vert L_{m_1}^{+}-L_{n_1}^{+}\Vert \right|<\frac{{\epsilon }_1}{3},$

but $\Vert L_{m_1}^{+}-L_{n_1}^{+}\Vert \ge {\epsilon }_1$ and $\Vert x_{n_1}\left(t\right)-L_{n_1}^{+}+x_{m_1}\left(t\right)-L_{m_1}^{+}\Vert <\frac{2{\epsilon }_1}{3}$, therefore:

${\epsilon }_1\le \Vert L_{m_1}^{+}-L_{n_1}^{+}\Vert -\Vert x_{n_1}\left(t\right)-L_{n_1}^{+}+x_{m_1}\left(t\right)-L_{m_1}^{+}\Vert <\frac{{\epsilon }_1}{3}$

${\epsilon }_1\le \Vert L_{m_1}^{+}-L_{n_1}^{+}\Vert <\frac{{\epsilon }_1}{3}+\Vert x_{n_1}\left(t\right)-L_{n_1}^{+}+x_{m_1}\left(t\right)-L_{m_1}^{+}\Vert$

${\epsilon }_1\le \Vert L_{m_1}^{+}-L_{n_1}^{+}\Vert <\frac{{\epsilon }_1}{3}+\frac{{\epsilon }_1}{3}+\frac{{\epsilon }_1}{3}={\epsilon }_1,$

which is a contradiction, therefore the sequence $\left\{L_n^{+}\right\}$ converges to a limit $L^{+}$.

Now let $\epsilon >0$, let $N\in \mathbb{N}$ and $\delta$ such that for $0<t+t_0<\delta$ :

$\Vert x_N\left(t\right)-x^{*}\left(t\right)\Vert <\frac{\epsilon }{3}$

$\Vert x_N\left(t\right)-L_N^{+}\Vert <\frac{\epsilon }{3}$

$\Vert L_N^{+}-L^{+}\Vert <\frac{\epsilon }{3},$

then:

$\Vert x^{*}\left(t\right)-L^{+}\Vert \le \Vert x_N\left(t\right)-x^{*}\left(t\right)\Vert +\Vert x_N\left(t\right)-L_N^{+}\Vert +\Vert L_N^{+}-L^{+}\Vert <\frac{\epsilon }{3}+\frac{\epsilon }{3}+\frac{\epsilon }{3}=\epsilon .$

Therefore $L^{+}$ is a right limit of $x^{\mathrm{<mo>\; *\; </mo>}}$ at $t=t_0$, hence the right limits exist. Repeating the process for the left limits we conclude that the left limits of $x^{\mathrm{<mo>\; *\; </mo>}}$ exist where required.

Now let $t_0\in \left[\mathrm{0,1}\right]$, $t_0$ not a dyadic rational, then each $x_n$ is continuous at $t_0$, therefore $\underset{n\to \infty }{\lim }{x}_n=x^{*}$ is continuous at $t_0$. Then we have $x^{\mathrm{<mo>\; *\; </mo>}}\in {\mathcal{C}}_{\Delta }$, which concludes the proof.

We make use of the following:

Lemma 9. [20] Let the dyadic intervals ${\Delta }_k\left(\mathcal{l}\right)=\left[\frac{\mathcal{l}-1}{2^k},\frac{\mathcal{l}}{2^k}\right)$, $\mathcal{l}\in {\Gamma }_k$ with

${\Gamma }_k=\left\{1,\cdots ,2^k\right\}$ and let ${\chi }_{{\Delta }_k\left(\mathcal{l}\right)}\left(t\right)$ be the characteristic function of the interval ${\Delta }_k\left(\mathcal{l}\right)$, then:

${\pi }_k{\displaystyle {\int }_0^t}{\chi }_{{\Delta }_k\left(\mathcal{l}\right)}\left(s\right)\text{d}s=\{\begin{array}{ll}0,\hfill & 0\le t<\frac{\mathcal{l}-1}{2^k}\hfill \\ \frac{1}{2^{k+1}},\hfill & \frac{\mathcal{l}-1}{2^k}\le t<\frac{\mathcal{l}}{2^k}\hfill \\ \frac{1}{2^k},\hfill & \frac{\mathcal{l}}{2^k}\le t<\frac{\mathcal{l}-1}{1}\hfill \end{array}$ (46)

Now we state:

Theorem 10. Let the Banach space ${\mathcal{C}}_{\Delta }$, the approximated monodromy operator ${\mathcal{U}}_k$ converges in the strong operator sense to the monodromy operator U.

Proof. Let ${\pi }_k$ denote the projection mapping that takes the Walsh expansion $f\left(t\right)={\displaystyle {\sum }_{\mathcal{l}=0}^{\infty }}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$ to the Walsh approximation $\hat{f}\left(t\right)={\displaystyle {\sum }_{\mathcal{l}=0}^{2^k}}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$. Without risk of confusion we will denote as ${\mathcal{U}}_k\phi \left(\theta \right)$ the operator ${\mathcal{U}}_k{\pi }_k\phi \left(\theta \right)$ and as ${\mathcal{T}}_k\phi \left(\theta \right)$ to the operator ${\mathcal{T}}_k{\pi }_k\phi \left(\theta \right)$.

Denote as $x\left(t\right)=T\phi \left(\theta \right)$ the solution of (20) corresponding to an initial condition $\phi \left(\theta \right)$ for $\theta \in \left[-\tau \mathrm{,0}\right]$. Likewise denote as $\hat{x}\left(t\right)={\mathcal{T}}_k\hat{\phi }\left(\theta \right)$ the solu- tion of (21). Clearly $x\left(t\right)=U\phi \left(\theta \right)$ and $\hat{x}\left(t\right)={\mathcal{U}}_k\phi \left(\theta \right)$ for $t\in \left[\omega -\tau \mathrm{,}\omega \right]$

Let $\hat{A}\left(t\right)={\pi }_{k}A\left(t\right)$, $\hat{B}\left(t\right)={\pi }_{k}B\left(t\right)$ and $\hat{\phi }\left(t\right)={\pi }_k\phi \left(t\right)$. Then:

$\Vert \hat{x}\left(t\right)-x\left(t\right)\Vert \le \Vert {\pi }_{k}x\left(t\right)-x\left(t\right)\Vert +\Vert \hat{x}\left(t\right)-{\pi }_{k}x\left(t\right)\Vert \mathrm{.}$ (47)

Since the solution $x\left(t\right)$ is continuous we have that the first term on the right converges to 0 as $k\to \infty$. We now observe that $\hat{x}\left(t\right)$ and $x\left(t\right)$ satisfy:

$\begin{array}{c}\hat{x}\left(t\right)=\phi \left(0\right)+{\pi }_k{\displaystyle {\int }_0^t}\hat{A}\left(s\right)\hat{x}\left(s\right)\text{d}s+{\pi }_k{\displaystyle {\int }_0^t}S\left(t-\tau \right)\hat{B}\left(s\right)\hat{x}\left(s-\tau \right)\text{d}s\\ +{\pi }_k{\displaystyle {\int }_0^t}\left(1-S\left(t-\tau \right)\right)\hat{B}\left(s\right)\hat{\phi }\left(s-\tau \right)\text{d}s\end{array}$ (48)

and

$\begin{array}{c}x\left(t\right)=\phi \left(0\right)+{\displaystyle {\int }_0^t}A\left(s\right)x\left(s\right)\text{d}s+{\displaystyle {\int }_0^t}S\left(t-\tau \right)B\left(s\right)x\left(s-\tau \right)\text{d}s\\ +{\displaystyle {\int }_0^t}\left(1-S\left(t-\tau \right)\right)B\left(s\right)\phi \left(s-\tau \right)\text{d}s,\end{array}$ (49)

respectively. Set $\mathcal{A}\left(t\right)=A\left(t\right)+S\left(t-\tau \right)B\left(t\right)$, $\mathcal{B}\left(t\right)=\left(1-S\left(t-\tau \right)\right)B\left(t\right)$, the same way set $\hat{\mathcal{A}}\left(t\right)$ and $\hat{\mathcal{B}}\left(t\right)$ as their respective approximations. Define $y_k\left(t\right)=\hat{x}\left(t\right)-{\pi }_{k}x\left(t\right)$. We have:

$\begin{array}{c}{y}_k\left(t\right)={\pi }_k{\displaystyle {\int }_0^t}\hat{\mathcal{A}}\left(s\right)\hat{x}\left(s\right)\text{d}s+{\pi }_k{\displaystyle {\int }_0^t}\hat{\mathcal{A}}\left(s\right){\pi }_{k}x\left(s\right)\text{d}s-{\pi }_k{\displaystyle {\int }_0^t}\mathcal{A}\left(s\right)x\left(s\right)\text{d}s\\ +{\pi }_k\left[{\displaystyle {\int }_0^t}\hat{\mathcal{B}}\left(s\right)\hat{\phi }\left(s-\tau \right)\text{d}s-\mathcal{B}\left(s\right)\phi \left(s-\tau \right)\text{d}s\right]\mathrm{,}\end{array}$ (50)

$\begin{array}{c}{y}_k\left(t\right)={\pi }_k{\displaystyle {\int }_0^t}\hat{\mathcal{A}}\left(s\right)\hat{x}\left(s\right)\text{d}s+{\pi }_k{\displaystyle {\int }_0^t}\left[\hat{\mathcal{A}}\left(s\right){\pi }_{k}x\left(s\right)-\mathcal{A}\left(s\right)x\left(s\right)\right]\text{d}s\\ +{\pi }_k{\displaystyle {\int }_0^t}\left[\hat{\mathcal{B}}\left(s\right)\hat{\phi }\left(s-\tau \right)-\mathcal{B}\left(s\right)\phi \left(s-\tau \right)\right]\text{d}s.\end{array}$ (51)

We have that the right side of (51) is piecewise constant, so if we define the dyadic intervals ${\Delta }_k\left(\mathcal{l}\right)=\left[\frac{\mathcal{l}-1}{2^k},\frac{\mathcal{l}}{2^k}\right)$, for $\mathcal{l}\in \Gamma$ with $\Gamma =\left\{1,\cdots ,2^k\right\}$, we can write:

${\pi }_k{\displaystyle {\int }_0^t}\hat{\mathcal{A}}\left(s\right)y_k\left(s\right)\text{d}s={\pi }_k{\displaystyle {\int }_0^t}{\displaystyle \underset{q=1}{\overset{2^k}{\sum }}}\hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)y_k\left(\frac{q-1}{2^k}\right){\chi }_{{\Delta }_k\left(q\right)}\left(s\right)\text{d}s\mathrm{,}$ (52)

where ${\chi }_{{\Delta }_k\left(q\right)}\left(t\right)$ is the characteristic function of the interval ${\Delta }_k\left(q\right)$. From Lemma 9, we will have for $t\in {\Delta }_k\left(\mathcal{l}\right)$ :

${\pi }_k{\displaystyle {\int }_0^t}\hat{\mathcal{A}}\left(s\right)y_k\left(s\right)=\frac{1}{2^k}{\displaystyle \underset{q=1}{\overset{\mathcal{l}-1}{\sum }}}\hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)y_k\left(\frac{q-1}{2^k}\right)+\frac{1}{2^{k+1}}\hat{\mathcal{A}}\left(\frac{\mathcal{l}-1}{2^k}\right)y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\mathrm{.}$ (53)

Since $y_k\left(t\right)$ is constant on the dyadic intervals then $y_k\left(t\right)=y_k\left(\frac{\mathcal{l}-1}{2^k}\right)$ for $t\in {\Delta }_k\left(\mathcal{l}\right)$. Therefore:

$\begin{array}{c}\Vert y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \le \frac{1}{2^k}\Vert {\displaystyle \underset{q=1}{\overset{\mathcal{l}-1}{\sum }}}\hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)y_k\left(\frac{q-1}{2^k}\right)\Vert +\frac{1}{2^{k+1}}\Vert \hat{\mathcal{A}}\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \Vert y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \\ +\underset{0\le t\le 1}{\max }\Vert {\pi }_k{\displaystyle {\int }_0^t}\left[\hat{\mathcal{B}}\left(s\right)\hat{\phi }\left(s-\tau \right)-\mathcal{B}\left(s\right)\phi \left(s-\tau \right)\right]\text{d}s\Vert \\ +\underset{0\le t\le 1}{\max }\Vert {\pi }_k{\displaystyle {\int }_0^t}\left[\hat{\mathcal{A}}\left(s\right){\pi }_k{x}_1\left(s\right)-\mathcal{A}\left(t\right)x_1\left(t\right)\right]\text{d}s\Vert .\end{array}$ (54)

Define $C_k^{\mathrm{<mo>\; *\; </mo>}}=1-\frac{1}{2^{k+1}}\underset{q\in \Gamma }{\max }\Vert \hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)\Vert$, and assume $C_k^{*}>0$. This restraint is not significant since $A\left(t\right)$ and $B\left(t\right)$ are bounded, therefore $\hat{A}\left(t\right)$ and $\hat{B}\left(t\right)$ are bounded, and $\frac{1}{2^{k+1}}\to 0$ as $k\to \infty$.

We have:

$\Vert {\displaystyle \underset{q=1}{\overset{\mathcal{l}-1}{\sum }}}\hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)y_k\left(\frac{q-1}{2^k}\right)\Vert \le \underset{q\in \Gamma }{\max }\Vert \hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)\Vert \Vert {\displaystyle \underset{q=1}{\overset{\mathcal{l}-1}{\sum }}}{y}_k\left(\frac{q-1}{2^k}\right)\Vert .$ (55)

From all of the above and defining $D_k^{\mathrm{<mo>\; *\; </mo>}}\triangleq \frac{1}{2^k}\frac{A_k^{\mathrm{<mo>\; *\; </mo>}}}{C_k^{\mathrm{<mo>\; *\; </mo>}}}$, with $A_k^{*}=\underset{q\in \Gamma }{\max }\Vert \hat{\mathcal{A}}\left(\frac{q-1}{2^k}\right)\Vert$ and denoting:

$\begin{array}{l}{B}_k^{\mathrm{<mo>\; *\; </mo>}}=\frac{1}{C_k^{\mathrm{<mo>\; *\; </mo>}}}[\underset{0\le t\le 1}{\max }\Vert {\pi }_k{\displaystyle {\int }_0^t}\left[\hat{\mathcal{A}}\left(s\right){\pi }_{k}x\left(s\right)-\mathcal{A}\left(t\right)x\left(t\right)\right]\text{d}s\Vert \\ +\underset{0\le t\le 1}{\max }\Vert {\pi }_k{\displaystyle {\int }_0^t}\left[\hat{\mathcal{B}}\left(s\right)\hat{\phi }\left(s-\tau \right)-\mathcal{B}\left(s\right)\phi \left(s-\tau \right)\right]\text{d}s\Vert ]\mathrm{,}\end{array}$

we arrive at:

$\Vert y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \le B_k^{*}+D_k^{*}{\displaystyle \underset{q=1}{\overset{\mathcal{l}-1}{\sum }}}\Vert y_k\left(\frac{q-1}{2^k}\right)\Vert .$ (56)

Which gives:

$\Vert y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \le B_k^{\mathrm{<mo>\; *\; </mo>}}{\left(1+D_k^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\mathcal{l}-1}\mathrm{.}$ (57)

Indeed, for $\mathcal{l}=1$, $\Vert y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \le B_k^{\mathrm{<mo>\; *\; </mo>}}$, and assuming $B_k^{*}+D_k^{*}{\displaystyle \underset{q=1}{\overset{\mathcal{l}-2}{\sum }}}\Vert y_k\left(\frac{q-1}{2^k}\right)\Vert \le B_k^{*}{\left(1+D_k^{*}\right)}^{\mathcal{l}-2}$ we have:

$\begin{array}{c}\Vert y_k\left(\frac{\mathcal{l}-1}{2^k}\right)\Vert \le B_k^{*}+D_k^{*}{\displaystyle \underset{q=1}{\overset{\mathcal{l}-2}{\sum }}}\Vert y_k\left(\frac{q-1}{2^k}\right)\Vert +D_k^{*}\Vert y_k\left(\frac{\mathcal{l}-2}{2^k}\right)\Vert \\ \le B_k^{*}{\left(1+D_k^{*}\right)}^{\mathcal{l}-2}+D_k^{*}{B}_k^{*}{\left(1+D_k^{*}\right)}^{\mathcal{l}-2}\\ \le B_k^{*}{\left(1+D_k^{*}\right)}^{\mathcal{l}-1}.\end{array}$ (58)

For $t\in {\Delta }_k\left(\mathcal{l}\right)$ we have $\frac{\mathcal{l}-1}{2^k}\le t$, then $\mathcal{l}-1\le 2^{k}t$, therefore:

$\Vert y_k\left(t\right)\Vert \le B_k^{\mathrm{<mo>\; *\; </mo>}}{\left(1+D_k^{\mathrm{<mo>\; *\; </mo>}}\right)}^{2^{k}t}\mathrm{.}$ (59)

Taking into account the definition of $D_k^{\mathrm{<mo>\; *\; </mo>}}$ we have:

${\left(1+D_k^{*}\right)}^{2^{k}t}\le {\left(1+\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}\right)}^{2^{k}t}\le {\left(1+\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}\right)}^{2^{k}t-\frac{A_k^{*}}{2}t}{\left(1+\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}\right)}^{\frac{A_k^{*}}{2}t}.$ (60)

Since $t\in \left[\mathrm{0,1}\right]$ and we assumed $1-\frac{1}{2^{k+1}}{A}_k^{*}>0$, then:

${\left(1+\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}\right)}^{2^{k}t-\frac{A_k^{*}}{2}t}\le {\left(1+\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}+\frac{1}{2!}{\left(\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}\right)}^2+\cdots \right)}^{2^{k}t-\frac{A_k^{*}}{2}t}\le {\text{e}}^{A_k^{*}t}\le {\text{e}}^{A_k^{*}}$ (61)

and

${\left(1+\frac{A_k^{*}}{2^k-\frac{A_k^{*}}{2}}\right)}^{\frac{A_k^{*}}{2}t}={\left(1+D_k^{*}\right)}^{\frac{A_k^{*}}{2}t}\le {\left(1+D_k^{*}\right)}^{\frac{A_k^{*}}{2}}.$ (62)

From (61), (62) and (59) we arrive at:

$\Vert y_k\left(t\right)\Vert \le B_k^{\mathrm{<mo>\; *\; </mo>}}{\text{e}}^{A_k^{\mathrm{<mo>\; *\; </mo>}}}{\left(1+D_k^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\frac{A_k^{\mathrm{<mo>\; *\; </mo>}}}{2}}\mathrm{,}\forall t\in \left[\mathrm{0,1}\right]\mathrm{.}$ (63)

The expresion ${\text{e}}^{A_k^{\mathrm{<mo>\; *\; </mo>}}}{\left(1+D_k^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\frac{A_k^{\mathrm{<mo>\; *\; </mo>}}}{2}}$ is bounded, and clearly $x\in {\mathcal{C}}_{\Delta }$, and since $\phi \left(t-\omega \right)\in {\mathcal{C}}_{\Delta }$, we have $B_k^{\mathrm{<mo>\; *\; </mo>}}\to 0$ as $k\to \infty$, therefore we have $\Vert y_k\Vert \to 0$ as $k\to \infty$, this is:

$\Vert {\mathcal{T}}_k\phi -T\phi \Vert \underset{k\to \infty }{\to }\mathrm{0,}\forall \phi \in {\mathcal{C}}_{\Delta }\mathrm{,}$ (64)

from where it follows immediately that:

$\Vert {\mathcal{U}}_k\phi -U\phi \Vert \underset{k\to \infty }{\to }\mathrm{0,}\forall \phi \in {\mathcal{C}}_{\Delta }\mathrm{,}$ (65)

which concludes the proof.

Corollary 1. If $\phi$ is Lipschitz, then the error of the approximation ${\mathcal{U}}_k$, satisfies $\Vert {\mathcal{U}}_k\phi -U\phi \Vert \in O\left(\frac{1}{2^k}\right)$.

Proof. From (59), we have that since $\phi$ is Lipschitz and $x\left(t\right)$ 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:

Lemma 11. [21] Let ${\displaystyle {\sum }_{\mathcal{l}=0}^{\infty }}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$ be the Walsh expansion of a function $f\left(t\right)$, if ${\displaystyle {\sum }_{\mathcal{l}=0}^{\infty }}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$ converges to zero everywhere, then $c_i=0$, $i\ge 0$.

Now we are ready to prove:

Theorem 12. Let $x\left(t\right)=T\phi \left(\theta \right)$, the solution of (20) corresponding to an initial condition $\phi \left(\theta \right)$ for $\theta \in \left[-\tau \mathrm{,0}\right]$. Likewise let $\hat{x}\left(t\right)={\mathcal{T}}_k\hat{\phi }\left(\theta \right)$, the solution of (21), then $\hat{x}\left(t\right)={\pi }_{k}x\left(t\right)$.

Proof. We have $\hat{x}\left(t\right)={\displaystyle {\sum }_{\mathcal{l}=0}^{2^k}}{c}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$ and ${\pi }_{k}x\left(t\right)={\displaystyle {\sum }_{\mathcal{l}=0}^{2^k}}{d}_{\mathcal{l}}{w}_{\mathcal{l}}\left(t\right)$, for some coefficients $c_{\mathcal{l}}$ and $d_{\mathcal{l}}$, then :

$\Vert {\pi }_{k}x\left(t\right)-\hat{x}\left(t\right)\Vert \le \Vert x\left(t\right)-\hat{x}\left(t\right)\Vert +\Vert {\pi }_{k}x\left(t\right)-x\left(t\right)\Vert \mathrm{.}$ (66)

Since the solution $x\left(t\right)$ is continuous and from (64) we have that the two terms on the right converge to zero, therefore ${\displaystyle {\sum }_{\mathcal{l}=0}^{2^k}}\left(c_{\mathcal{l}}-d_{\mathcal{l}}\right)w\mathcal{l}\left(t\right)$ converges to zero uniformly, hence everywhere, and from Lemma 11. we have $c_{\mathcal{l}}=d_{\mathcal{l}}$, for $\mathcal{l}=0,\cdots ,2^k$, this is $\hat{x}\left(t\right)={\pi }_{k}x(t)$

Corollary 2. ${\mathcal{U}}_k\phi \left(t\right)={\pi }_{k}U\phi \left(t\right)$,

Lemma 13. Let $X\subset {\mathcal{C}}_{\Delta }$ compact, then for very $\epsilon >0$, there exists $k\in \mathbb{N}$, such that:

$\Vert {\pi }_{k}x-x\Vert <\epsilon ,\forall x\in X.$ (67)

Proof. Let $\epsilon >0$, since X is compact, there exist finite open balls of radius $\frac{\epsilon }{3}$ centered around finite $x_1\mathrm{,}\cdots \mathrm{,}{x}_N\in X$ that cover X, then $\forall x\in X$, $\Vert x-x_i\Vert <\frac{\epsilon }{3}$, for some $i\in \left\{\mathrm{1,}\cdots ,N\right\}$. Let k such that for any $i\in \left\{\mathrm{1,}\cdots ,N\right\}$, $\Vert {\pi }_k{x}_i-x_i\Vert <\frac{\epsilon }{3}$. Let $x\in X$, then for some $i\in \left\{\mathrm{1,}\cdots ,N\right\}$ :

$\begin{array}{c}\Vert {\pi }_{k}x-x\Vert \le \Vert {\pi }_{k}x-{\pi }_k{x}_i\Vert +\Vert x-x_i\Vert +\Vert {\pi }_k{x}_i-x_i\Vert \\ \le \Vert {\pi }_k\Vert \Vert x-x_i\Vert +\Vert x-x_i\Vert +\Vert {\pi }_k{x}_i-x_i\Vert <\epsilon .\end{array}$

Since there are finitely many $x_i$, the desired result follows immediately.

The solution of (1) will be given by:

$x\left(t\right)=X\left(t,0\right)\phi \left(0\right)+{\displaystyle {\int }_{-\tau }^0}X\left(t,s+\tau \right)B\left(s+\tau \right)\phi \left(s\right)\text{d}s,$ (68)

where X is a solution matrix such that $X\left(0,0\right)=I$ and $X\left(t,0\right)=0$ for $t<0$ [13] . If $\phi \in {\mathcal{C}}_{\Delta }$ then the solution $x\left(t\right)$ will be continuous and the solution matrix X will be the same that in $\mathcal{C}$. Furthermore we will have a bounded operator:

$\Vert x\left(t\right)\Vert =\Vert T\left(\phi \left(\theta \right)\right)\Vert \le \Vert T\Vert \Vert \phi \left(\theta \right)\Vert \mathrm{,}$ (69)

and like in $\mathcal{C}$, we will have that T maps arbitrary bounded sequences into equicontinuous sequences, then by the Arzelá-Ascoli Theorem we will have that T is compact in ${\mathcal{C}}_{\Delta }$.

We now have:

Theorem 14 The approximated monodromy operator ${\mathcal{U}}_k$ converges uni- formly to the monodromy operator U.

Proof. Let $\bar{B}$ denote the closed unit ball on ${\mathcal{C}}_{\delta }$. Since T is compact, then the image $T\bar{B}$ is also compact, and by Lemma 13 we will have that for any $\epsilon >0$ there exists k such that:

$\Vert {\mathcal{T}}_k-T\Vert =\underset{\Vert x\Vert =1}{\sup }\Vert {\pi }_{k}Tx-Tx\Vert <\epsilon .$ (70)

The fact that $\Vert {\mathcal{U}}_k-U\Vert \to 0$ as $k\to \infty$ follows immediately.

With this we can establish convergence of the spectrum:

Theorem 15. The spectrum of the approximated monodromy operator ${\mathcal{U}}_k$ converges to the spectrum of the monodromy operator U. More precisely, for any open set $\mathcal{V}$, such that $\sigma \left(\mathcal{U}\right)\subset \mathcal{V}$, then there exist K such that $\sigma \left({\mathcal{U}}_k\right)\subset \mathcal{V}$, for any $k\ge K$. Furthermore, let ${\lambda }_0\in \sigma \left(U\right)$ and let $\Gamma$ be a small circle centered at ${\lambda }_0$ such that any other eigenvalue of U is outside $\Gamma$, then the sum of the multiplicities of the eigenvalues of $\mathcal{U}$ within $\Gamma$ will be equal to the multiplicity of ${\lambda }_0$.

Proof. Let $\Gamma$ be a small circle with center at ${\lambda }_0\in \sigma \left(U\right)$, such that any other eigenvalue of U is outside this circle. Then the spectral projection of ${\lambda }_0$ is given by:

$P=\frac{1}{2\text{π}i}{\displaystyle {\oint }_{\Gamma }}{\left(\lambda I-U\right)}^{-1}\text{d}\lambda \mathrm{.}$ (71)

The spectral projection of the spectrum of ${\mathcal{U}}_k$ inside $\Gamma$ will be given by:

$P_k=\frac{1}{2\text{π}i}{\displaystyle {\oint }_{\Gamma }}{\left(\lambda I-{\mathcal{U}}_k\right)}^{-1}\text{d}\lambda \mathrm{.}$ (72)

Since $\Vert U-{\mathcal{U}}_k\Vert \to 0$, then $P_k$ converges uniformly to P, therefore, there must exist k such that part of the spectrum of ${\mathcal{U}}_k$ is in $\Gamma$.

Since P and $P_k$ are finite dimensional operators and since $P_k$ converges uniformly to P, we will have that for sufficiently large k:

$dim\left(P\right)=dim\left(P_{k}v\right)\mathrm{,}$ (73)

this is, the algebraic multiplicity of ${\lambda }_0$ will be the same as the algebraic multiplicity of the spectrum of ${\mathcal{U}}_k$ contained in $\Gamma$.

6. 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:

$\ddot{x}+\left(\alpha +\beta cost\right)x=\left(\gamma cost\right)x\left(t-\tau \right)\mathrm{.}$ (74)

The operator ${\mathcal{T}}_k$ in (33) provides a natural way to approximate the solution of a periodic delay differential equation. Consider Equation (74), with $\alpha =5.35$,

$\beta =2.5$, $\gamma =0.5$, and $\tau =\frac{\text{2π}}{128}$. Figure 3 and Figure 4, show the comparison

between the solution obtained by simulation and the approximation of the solution obtained from ${\mathcal{T}}_k$, for $k=7$ and $k=10$ respectively, for $t\in \left[\mathrm{0,2}\text{π}\right]$.

Figure 5 and Figure 6, show the magnitude of the error for the appro- ximation of the solution obtained with the approximated operator ${\mathcal{T}}_k$, for $k=7$ and $k=10$, respectively.

Figure 3. Solution of (74) for $t\in \left[\mathrm{0,2}\text{π}\right]$. The solid line corresponds to the simulated solution and the dashed line corresponds to the approximated solution for $k=7$.

Figure 4. Solution of (74) for $t\in \left[\mathrm{0,2}\text{π}\right]$. The solid line corresponds to the simulated solution and the dashed line corresponds to the approximated solution for $k=10$.

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

Figure 6. Magnitude of error of the approximate solution corresponding to $k=10$.

7. 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:

$\ddot{x}\left(t\right)+\left(\alpha +\beta \cos \left(t\right)\right)x\left(t\right)=\gamma \cos \left(t\right)x\left(t-\tau \right).$ (75)

Stability charts are used to to determine the stability of periodic differential equations with respect to some parameters. Figure 7 shows the stability diagram

of (75) for the plane $\alpha \beta$ with $\tau =\frac{\text{2π}}{32}$ and $\gamma =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:

$\ddot{x}\left(t\right)+\left[\alpha +\beta \cos \left(2\text{π}t\right)\right]x\left(t\right)=\delta x\left(t-\tau \right).$ (76)

This equation has been studied in [8] . Figure 8 shows the stability diagram for the parametric plane $\alpha \beta$ with $\delta =0.15$ and $\tau =\frac{\text{2π}}{16}$.

8. 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

Figure 7. Stability diagram for the parametric plane αβ of Equation (75). Asymptotically stable zones (grey) and unstable zones (white) are shown.

Figure 8. Stability diagram for the parametric plane αβ of Equation (76). 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 $A\left(t\right)$ and $B\left(t\right)$ 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.

Fund

This work was supported by CONACyT, Mexico, CVU 418725.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Kolmanovskii, V. and Myshkis, A. (2013) Introduction to the Theory and Applications of Functional Differential Equations. Springer Science & Business Media, Vol. 463.
[2]	Cui, H.-Y., Han, Z.-J. and Xu, G.-Q. (2016) Stabilization for Schrdinger Equation with a Time Delay in the Boundary Input. Applicable Analysis, 95, 963-977. https://doi.org/10.1080/00036811.2015.1047830
[3]	Srividhya, J. and Gopinathan, M. (2006) A Simple Time Delay Model for Eukaryotic Cell Cycle. Journal of Theoretical Biology, 241, 617-627. https://doi.org/10.1016/j.jtbi.2005.12.020
[4]	Datko, R., Lagnese, J. and Polis, M. (1986) An Example on the Effect of Time Delays in Boundary Feedback Stabilization of Wave Equations. SIAM Journal on Control and Optimization, 24, 152-156. https://doi.org/10.1137/0324007
[5]	Datko, R. and You, Y. (1991) Some Second-Order Vibrating Systems Cannot Tolerate Small Time Delays in Their Damping. Journal of Optimization Theory and Applications, 70, 521-537. https://doi.org/10.1007/BF00941300
[6]	Gilsinn, D.E., Potra, F.A., et al. (2006) Integral Operators and Delay Differential Equations. Journal of Integral Equations and Applications, 18, 297-336. https://doi.org/10.1216/jiea/1181075393
[7]	Breda, D., Maset, S. and Vermiglio, R. (2014) Pseudospectral Methods for Stability Analysis of Delayed Dynamical Systems. International Journal of Dynamics and Control, 2, 143-153. https://doi.org/10.1007/s40435-013-0041-x
[8]	Insperger, T. and Stepan, G. (2002) Semi-Discretization Method for Delayed Systems. International Journal for Numerical Methods in Engineering, 55, 503-518. https://doi.org/10.1002/nme.505
[9]	Vazquez, E.A. and Collado, J. (2016) Monodromy Operator Approximation of Periodic Delay Differential Equations by Walsh Functions. 13th International Conference on Electrical Engineering, Computing Science and Automatic Control (CCE), IEEE.
[10]	Chen, W.-L. and Shih, Y.-P. (1978) Shift Walsh Matrix and Delay-Differential Equations. IEEE Transactions on Automatic Control, 23, 1023-1028. https://doi.org/10.1109/TAC.1978.1101888
[11]	Kreyszig, E. (1978) Introductory Functional Analysis with Applications, Ser. Wiley Classics Library. John Wiley & Sons.
[12]	Stokes, A. (1962) A Floquet Theory for Functional Differential Equation. Proceedings of the National Academy of Sciences of the United States of America, 48, 1330. https://doi.org/10.1073/pnas.48.8.1330
[13]	Halanay, A. (1966) Differential Equations: Stability, Oscillations, Time Lags. Academic Press, 23.
[14]	Yakubovich, V.A. and Starzhinski, V.M. (1975) Linear Differential Equations with Periodic Coefficients. Vol. 1, John Wiley, New York.
[15]	Magnus, W. and Winkler, S. (2004) Hill’s Equation, Ser. Dover Books on Mathematics Series. Dover Publications.
[16]	Walsh, J.L. (1923) A Closed Set of Normal Orthogonal Functions. American Journal of Mathematics, 45, 5-24. https://doi.org/10.2307/2387224
[17]	Golubov, B., Efimov, A. and Skvortsov, V. (2012) Walsh Series and Transforms: Theory and Applications. Springer Science & Business Media, Berlin, Vol. 64.
[18]	Ahmed, N. and Rao, K.R. (2012) Orthogonal Transforms for Digital Signal Proces-sing. Springer Science & Business Media, Berlin.
[19]	Karanam, V., Frick, P. and Mohler, R. (1978) Bilinear System Identification by Walsh Functions. IEEE Transactions on Automatic Control, 23, 709-713. https://doi.org/10.1109/TAC.1978.1101806
[20]	Gát, G. and Toledo, R. (2015) Estimating the Error of the Numerical Solution of Linear Differential Equations with Constant Coefficients via Walsh Polynomials. Acta Mathematica Academiae Paedagogicae Nyíregyháziensis, 31, 309-330.
[21]	Fine, N.J. (1949) On the Walsh Functions. Transactions of the American Mathematical Society, 65, 372-414. https://doi.org/10.1090/S0002-9947-1949-0032833-2
[22]	Kwong, C. and Chen, C. (1981) The Convergence Properties of Block-Pulse Series. International Journal of Systems Science, 12, 745-751. https://doi.org/10.1080/00207728108963780
[23]	Rao, G. and Srinivasan, T. (1978) Remarks on “Author’s Reply”; to “Comments on Design of Piecewise Constant Gains for Optimal Control via Walsh Functions”. IEEE Transactions on Automatic Control, 23, 762-763. https://doi.org/10.1109/TAC.1978.1101782
[24]	Cheng, C., Tsay, Y. and Wu, T. (1977) Walsh Operational Matrices for Fractional Calculus and Their Application to Distributed Systems. Journal of the Franklin Institute, 303, 267-284. https://doi.org/10.1016/0016-0032(77)90029-1
[25]	Gulamhusein, M. (1973) Simple Matrix-Theory Proof of the Discrete Dyadic Convolution Theorem. Electronics Letters, 9, 238-239. https://doi.org/10.1049/el:19730172
[26]	Hale, J.K. and Lunel, S.M.V. (2013) Introduction to Functional Differential Equations. Springer Science & Business Media, Berlin, Vol. 99.
[27]	Deb, A., Sarkar, G. and Sen, S.K. (1994) Block Pulse Functions, the Most Fundamental of All Piecewise Constant Basis Functions. International Journal of Systems Science, 25, 351-363. https://doi.org/10.1080/00207729408928964
[28]	Norkin, S. (1973) Introduction to the Theory and Application of Differential Equations with Deviating Arguments. Academic Press, Cambridge, Vol. 105. https://doi.org/10.1016/S0076-5392(08)62969-0
[29]	Franco, C. and Collado, J. (2017) A Novel Discriminant Approximation of Periodic Differential Equations.