Pseudospectral discretization of delay differential equations in sun-star formulation: Results and conjectures

In this paper we study the pseudospectral approximation of delay differential equations formulated as abstract differential equations in the \begin{document}$ \odot* $\end{document} -space. This formalism also allows us to define rigorously the abstract variation-of-constants formula, where the \begin{document}$ \odot* $\end{document} -shift operator plays a fundamental role. By applying the pseudospectral discretization technique we derive a system of ordinary differential equations, whose dynamics can be efficiently analyzed by existing bifurcation tools. To better understand to what extent the resulting finite-dimensional system "mimics" the dynamics of the original infinite-dimensional one, we study the pseudospectral approximations of the \begin{document}$ \odot* $\end{document} -shift operator and of the \begin{document}$ \odot* $\end{document} -generator in the supremum norm, which is the natural choice for delay differential equations, when the discretization parameter increases. In this context there are still open questions. We collect the most relevant results from the literature and we present some conjectures, supported by various numerical experiments, to illustrate the behavior w.r.t. the discretization parameter and to indicate the direction of ongoing and future research.

1. Introduction. Delay differential equations (DDEs), together with renewal equations (REs) and systems which couple REs and DDEs (REs/DDEs), are recognized as a fundamental tool for modelling phenomena in many fields, including, for instance, population dynamics and control theory. For this reason, in the last years the interest in the study of the dynamics of delay models has been increasing and important challenges, in particular numerical, have been identified. Indeed, delay equations describe infinite-dimensional dynamical systems, and theoretical results should be complemented with efficient numerical methods to approximate solutions of initial value problems [3,4,5,7,17,15,16,19,44,43], boundary value problems [48,49,50], and to investigate the stability of equilibria and periodic solutions [10,11,12,13,14,47,52,53,66]. In applications the attention is focused not only on the approximation of the dynamical properties for some given parameter values, but also on how such properties change when varying some parameters. In particular, it is interesting to identify the critical thresholds, called bifurcation points, and to draw stability charts in two or more parameters.
While the theory of DDEs is well developed, see for instance [6,25,39,65], the numerical methods are still unsatisfactory: some software packages, like DDE-BIFTOOL [27] and Knut [59], are available for the stability and bifurcation analysis of equations with discrete delays, but cannot be applied, for instance, in the case of distributed delays. To address this problem, the authors of [8,56] proposed an alternative method, which consists in reformulating the original equation as a nonlinear abstract differential equation (ADE) in the state space, and then applying the pseudospectral discretization (PSD) approach. Then, the dynamics of the resulting nonlinear systems of ordinary differential equations (ODEs) can be investigated by one of the continuation packages for the bifurcation analysis of ODEs, like MatCont [22]. In [8] the approach has been applied also to REs and REs/DDEs. Moreover, it has been proved that there is a one-to-one correspondence between the equilibria, and the stability is accurately described. Numerous experiments give evidence that the PSD approach converges for periodic solutions, too [8,9], but at the moment the rigorous proof is not yet completed. In this context, the approximation of the so-called shift operator, which translates a function to the right while extending it with a constant value, plays a fundamental role. An important question is how accurately the shift operator is approximated by the resulting discrete operator, w.r.t. the discretization parameter. The action of the shift operator is described by the trivial equation y (t) = 0 considered as a DDE with delay τ > 0, and is related to the partial differential equation (PDE)  [25, pag. 39]. The PSD approach has been widely applied for the numerical approximation of PDEs [18,31], and many results are available in the literature. Nevertheless, the estimates of the approximation of the shift operator regard weighted 2-norms in suitable Sobolev spaces [26,29,30,35,33,34,45,40,58,62,64], whereas, to our knowledge, no results are available for the supremum norm. The supremum norm is largely used when studying DDEs with the Banach space of continuous functions as state space [25,39]. Nevertheless this choice is not unique and, having in mind especially some applications to control theory, L p -state space has been considered in [6,19,65].
2. * -formulation of delay differential equations. In this section we briefly summarize the basic results on the * -reformulation of DDEs, with the aim to introduce the ADE and to emphasize the role of the * -shift operator in this abstract framework, in view of the introduction of the PSD approach in Section 3. For a deeper insight into the theory, we refer to [20,23,25]. Let d be a positive integer, τ > 0, and consider the Banach space A nonlinear autonomous DDE is where G : Y → R d is a continuous nonlinear function, and the history y t ∈ Y is defined as y t (θ) = y(t + θ), θ ∈ [−τ, 0].  (2) is either a continuously differentiable function y : R → R d satisfying (2) everywhere, or a continuous function y : [−τ, t + ) → R d , 0 < t + ≤ ∞, which is differentiable on (0, t + ) and satisfies (2) for 0 < t < t + .
An example of a solution of the first type in Definition 2.1 is given by a periodic solution y of (2), which is differentiable and periodic with period h, i.e. y(t + h) = y(t), t ∈ R, and satisfies (2) everywhere.
The * -formulation of DDEs aims to represent the equation (2) in an abstract framework, where the action of the dynamical system can be described in terms of two different components: the action of "shift", which is linear and exactly the same for every equation, and an action of "extension", which captures the effect of the nonlinear function G. In the * -framework the perturbation theory can be applied, and the variation-of-constants formula can be rigorously defined [20,23,25].
The main idea of the * -formulation is to embed Y into the larger space Y * given by which is a Banach space equipped with the norm The embedding j : Y → Y * is given by j(ψ) = (ψ(0), ψ). As for the "shift" component, it is well known that the generator of translation is differentiation. So we introduce the operator A * 0 : where "Lipschitz continuous" is a short-hand way to indicate an equivalence class that contains a Lipschitz continuous function. The operator A * 0 is the infinitesimal generator in the weak * -sense of the * -shift semigroup {T * 0 (t)} t≥0 , where the * -shift operator is defined by But we do not want to lose sight of Y , so we introduce the space where, as before, "continuous" indicates an equivalence class that contains a continuous function. Since the delay is bounded, the space Y is -reflexive in the sense that j is a bijection between Y and Y , with j −1 (α, ψ) = ψ, see e.g. [20,25]. Finally, by describing the rule for the extension in Y * through the nonlinear operator G : Y → Y * given by we can construct the operator A * : whose action describes both the "shift" and "extension", through A * 0 and G respectively, while its domain does not depend on the particular equation. Notice that G inherits from G the Lipschitz continuity. The operator (5) allows to represent the DDE (2) as an evolutionary system in the Y * space, namely the semilinear ADE in Y * djv(t) dt = A * (jv(t)).
For ψ ∈ Y , we define mild solution of the Cauchy problem the unique solution of the abstract integral equation (AIE) where the integral is a weak * -Riemann integral [20,Theorem 3.1]. From [20, Proposition 2.1] we have that the integral at the right-hand side of (8) defines an element of Y , so we can apply j −1 to both sides of (8) and write the following AIE in where T 0 (t) := j −1 T * 0 (t)j is the shift operator on Y . Following [25, Chapter III], we can conclude that v(t) = y t (·; ψ), where y(·; ψ) is a solution of the initial value problem (2) & (3).
The part of A * in Y is the operator jAj −1 , with A defined by whose action is independent of G, whereas now the information about the particular DDE is incorporated in the domain. For ψ ∈ D(A), the solution jv(t) of (7) is a classical solution, i.e. a continuously differentiable function jv(t) satisfying (6) for t ≥ 0 [20, Theorem 3.6].
The key role of the * -shift operator T * 0 (t) and of the generator A * 0 comes to light clearly in respectively (8) and (6). This will be also useful in the analysis of their discretization in Section 3. Notice that they describe the * -formulation of the trivial equation interpreted as a DDE.
We remark that the nonlinearity in the domain of (9) makes it hard to handle the perturbations of (2), since they amount to perturbing the domain of an infinitesimal generator. The main advantage of the * -formulation is indeed to avoid the dependence of the domains on G, thus treating operators defined on the same domain. So, although the * -framework requires some additional technical effort, the result yields a powerful abstract variation-of-constants formula for a larger class of perturbations than those bounded from Y into Y . The method also applies to other classes of delay equations, namely REs of the form with x t ∈ X := L 1 ((−τ, 0), R d ) and F : X → R d a smooth function, and REs/DDEs 3. * -pseudospectral discretization. In this section we introduce the PSD of the operator A * in (5) with the aim to derive a nonlinear system of ODEs that "mimics" the dynamics of the original nonlinear DDE (2). The basic idea of the PSD approach consists in approximating a function with a polynomial, which is represented as the interpolating polynomial on a suitable set of nodes in the domain of definition. Then "do to the interpolating polynomial what you would do to the function" [2].
In order to start from pointwise defined functions and to keep in touch with the original DDE, we will develop the PSD approach in the subspace represents the discretization of index M of Y , in the sense that every (α, ψ) is discretized by the element (α, ψ(Θ M )) ∈ Y M . The choice of the norm (13) is motivated by the norm of the state space Y * .
As further ingredients we introduce the restriction operator R M : Y → Y M as where and the prolongation operator P M : Y M → D(A * 0 ) ⊂ Y , which associates to (α, Ψ) ∈ Y M the pair (α, ψ M ), where ψ M is the M -degree polynomial interpolating (α, Ψ) at the nodes (11). It is important to underline that the construction of ψ M takes into account the domain condition in (4). By introducing the Lagrange (see e.g. [60, p. 34]). Hereafter, when it will be necessary to emphasize the dependence on (α, Ψ), we also use the notation where P M := j −1 P M . Note that R M P M is the identity operator on Y M , while P M R M j is the interpolation operator on Y relative to the Chebyshev nodes (11). Now we are ready to define both the PSD of A * 0 as the bounded operator and the nonlinear map where • denotes the composition of (nonlinear) operators. Their sum furnishes the discretization A M of A * , i.e.
So the PSD approach leads to the following nonlinear To write (18) explicitly, we need to find a matrix representation of the operator (16). Let us define the elements d ij := j (θ i ), for i, j = 0, . . . , M and the matrix which will play a fundamental role in the definition and analysis of the discretized equation (18).
. . , M , by denoting 1 the vector with all entries equal to 1 (we do not specify the dimension, since it is clear from the context), from (15) we get where ⊗ is the tensor product and I d denotes the identity matrix of dimension d.

Remark 1.
Since we have to handle ψ such that jψ ∈ D(A * 0 ), the Chebyshev points (11) are a good choice as interpolation nodes. Let us define the Lebesgue function  (11), we get the following error bound

and the Lebesgue constant
and, when also ψ is a Lipschitz continuous function we have with C a positive constant independent of M and ψ [51, pag. 269], and ψ − ψ M Y → 0 as M → ∞. Notice that the convergence results mentioned above apply also to ψ =ȳ t , whereȳ is a periodic solution of (2).
Finally, similarly to the abstract framework in Section 2, we complete the path by deriving also the integral equation associated to (18). Given the initial condition the initial value problem (18) ,M t and admits the following representation We underline the close analogy in structure between the integral formulations (22) and (8). We will return to this in Section 5.
In the paper [8] the PSD approach has been applied to the operator (9) obtaining the same ODE (18). So what did we gain by working in the * framework? It has furnished a deep insight into the connection between the infinite-dimensional and the finite-dimensional equations, including not only the differential equations (6) and (18), but also the integral equations (8) and (22). Moreover, the regularity condition of the functions in the domain has been relaxed w.r.t. those in D(A) in (9), and so we can handle Lipschitz continuous functions.
In order to use the ODE (18) to gain insight into the dynamics of the original DDE (2), we need to understand how the dynamical properties of the two equations are related. In the paper [8] it has been proved that the equilibria of (2) and (18) are in one-to-one correspondence, and the stability is the same for M large enough. In order to understand the behavior of more complicated solutions, like periodic solutions, it is fundamental to study the solution operators generated by the nonlinear equations, which are described by the integral equations (8) and (22). The integral formulation is indeed one of the advantages of the * -formalism and it is particularly useful for the convergence analysis of solution operators and, in turn, of periodic solutions. Indeed, thanks to the integral equations, we can relate the convergence of the nonlinear semigroups to the convergence of the trivial linear semigroups T 0,M (t) to T * 0 (t). Moreover (4) and (20) suggest that the matrix (16) "mimics" the infinitesimal generator A * 0 of the semigroup {T * 0 (t)} t≥0 associated to the trivial DDE (10). Therefore, to understand the dynamical behavior of (18), it is crucial to investigate the properties of the matrix A 0,M and of the corresponding  (20) it is clear that the matrix D M plays a major role. We remark that D M can be viewed as the Lagrange representation of the PSD w.r.t. the nodes (12) of the derivative operator D : with domain D(D) the set of equivalence classes containing a Lipschitz continuous function ψ such that ψ(0) = 0 [25, pag. 124]. In the following we collect some of the results about the matrix D M available in the literature, which are mainly obtained in the context of the approximations of PDEs of the form (1). For the sake of presentation hereafter we assume d = 1, and we use · for all the norms, omitting the subscripts, since the space we are working with will be clear from the context. Moreover in the numerical experiments we consider τ = 1.
Indeed the results for general τ > 0 can be easily obtained by a suitable scaling of the time variable.
All the experiments are made with Matlab 2018b, whose machine precision is Finally, hereafter we consider complex-valued functions and complex Banach spaces. For a detailed presentation of the complexification procedure, the interested reader is referred to [25,Sections III.7 and IV.2].
Since the derivative operator (24) is unbounded, we expect D M to diverge as M → ∞. The following proposition confirms this claim and defines the order of divergence w.r.t. M , which is determined by the asymptotic behavior of the Chebyshev polynomials.
Indeed the numerical results in Table 1 give evidence that D M = 2M 2 − 1 for τ = 1. We now study whether the discrete operator preserves the spectral properties of the continuous operator. The spectrum of the operator A * 0 contains only the eigenvalue λ = 0, and every non-null constant function is an eigenfunction. The spectrum of A 0,M contains, in addition to the eigenvalue λ = 0, also the eigenvalues of the matrix D M : the pseudospectral approximation introduces M eigenvalues which are "spurious", i.e., due only to the discretization procedure. Therefore it is important to study the behavior of the spectrum of D M when M → ∞.
The eigenvalues of D M are not explicitly known, but many theoretical and empirical results have been derived by various authors. In fact their spectral properties are relevant in the convergence analysis of the pseudospectral methods, otherwise known as spectral collocation, for first order hyperbolic PDEs.
So the first fundamental results have been established within this context [26,35,33,34,40,58]. Various further properties can be found also in [30,29,31,18,45,54,63,64]. Both the derivative operator in Hilbert space and the matrix D M are non-normal, and we recall the illuminating book [62], where the fundamental properties of non-normal operators in Hilbert spaces and matrices have been investigated through the analysis of their pseudospectra.
The following proposition allows us to analyze the position in the complex plane of the eigenvalues of D M . Concerning the asymptotic behavior of the spectral abscissa α(D M ), which is defined as where σ(D M ) denotes the spectrum of D M , the first result goes back to Dubiner [26], and was revised recently by Wang & Waleffe [64].  [62,63], the computation of the spectrum of pseudospectral differentiation matrices is extremely sensitive to rounding errors. In particular, the authors observed that eigenvalues λ ∈ C with Reλ < log are affected by rounding errors because the precision in the corresponding eigenvectors is lost. This phenomenon is observable in Figure 1 (left): for M = 64, a part of the computed spectrum falls to the left of the line Reλ = log and we expect such results to be affected by rounding errors. We note that the rounding errors affect mostly the eigenvalues which are large in modulus, while the computation of the α(D M ) seems to be more stable, (see Figure  1). The right panel of Figure 1 shows the oscillatory convergence of −α(D M ) log(M ) . We now focus our attention on the resolvent operators and their approximations. For any λ ∈ C \ {0}, the resolvent of A * 0 exists and it is the bounded operator given by for all θ ∈ [−τ, 0], which leads to the following bound of the resolvent norm Note that the upper bound depends on Reλ only, and, though it is finite for every λ ∈ C \ {0}, it increases exponentially when Reλ → −∞. This phenomenon has been already documented w.r.t. L 2 -norm in [62]. When β = 0, we get What can we say about the resolvent operator of the discrete operator A 0,M ? Given λ = 0, Proposition 3 states that there exists M 0 := M 0 (λ) such that λ / ∈ σ(D M ) for all M ≥ M 0 . Therefore for all M ≥ M 0 the resolvent operator of A 0,M is defined by the matrix Let λ ∈ C \ {0} and (β, ϕ) ∈ Y . The function is the solution of the following problem ψ (θ) = λψ(θ) − ϕ(θ), θ ∈ [−τ, 0], We denote by ψ M (θ) := ψ M (θ; λ, β, ϕ), θ ∈ [−τ, 0], the M -degree collocation polynomial for the problem (29) such that We state the following theorem. where and L M −1 is the interpolation operator relative to the nodes Θ M in (12). Moreover, for all M ≥ M 0 , ψ M admits the following representation which is the discrete counterpart of (28).
Proof. The proof is based on the approach of [14, Proposition 5.1]. By using the integral operator (26), we can rewrite both (29) and (30) where r M is defined in (32). We have that e M is a solution of (34) if and only if e M = Vz M with z M ∈ Y solution of the equation By observing that the operator I − λV is invertible with inverse   (36) and (37) we easily get the bound and (31) follows. Finally, by using the Lagrange representation (15) of ψ M and the matrix (19), we can easily obtain from the collocation equation (30) that ψ M is given by (33).
The bound (31) states that the error depends on the interpolation error of both, ψ and ϕ, and so we need to assume more regularity on ϕ to obtain convergence. It is sufficient to consider jϕ ∈ D(A * 0 ). Moreover when ϕ is an analytic function, we get the so-called spectral accuracy, i.e. the error ψ M − ψ decays as O(ρ −M ) for some ρ > 1 [60, Chapter 8]. In Figure 2 we plot the error ψ M − ψ when β = 0 and ϕ = 1. In this case ψ = 1−e λ· λ is analytic and L M −1 1 = 1. Notice the spectral accuracy and that the convergence is slower when λ has larger modulus. Finally for β = 0 and λ = 0, we get that

5.
The solution operators T * 0 (t) and T 0,M (t): Results and conjectures. In this section, we focus our attention on the connection between T * 0 and T 0,M , when M → ∞, which is a crucial question for the convergence analysis. Indeed this is motivated by the integral equations (8) and (22), which can be reformulated respectively as and where z(t) is G(v(t)) in (38) and G(P M v M (t)) in (39). We first introduce the operator H 0 (t) : and we reformulate the action of the * −shift operator as follows As usual α denotes either a scalar or the constant function. Similarly, by using (23) we can rewrite the discrete counterpart as which hightlighs that e D M t plays the role of H 0 (t) for all t ≥ 0. According to these reformulations and without loss of generality, we can assume that α = 0 and focus on the behavior of e D M t w.r.t. M for t ≥ 0. As already pointed out in Section 4, many results about D M and e D M t have been derived in the applications of the spectral collocation methods to hyperbolic PDEs. In fact, within this context, the space stability, i.e. the norm boundedness of e D M t w.r.t. M, plays a key role in the convergence analysis [18,28,40]. PDEs are usually studied in suitable Hilbert spaces and so the results mainly regard weighted 2-norms in Y M , which correspond to the discretization of L 2 norms (see [31,33,36,45,58]). Hyperbolic PDEs are strongly related to DDEs, but the latter ones require to work in infinite-dimensional Banach spaces endowed with supremum norm and, to our knowledge, no bounds are available in this case.
To address the question and have a first insight, we have performed some numerical simulations. The computation of the exponential matrices, especially nonnormal, can be difficult and strongly affected by rounding errors, and the research on this subject is still in progress (see for instance [42] and the references therein). Here we use the built-in Matlab function expm, which is is based on the algorithms in [1,41].
The results in Figure 3 give also evidence that there existst M where e D M t changes its behavior from not convergent to convergent, and that the convergence is attained for t ≥t M . The relevant fact emerging from the simulations is that t M ∈ (τ, 1.5τ ) for sufficiently large M. In truth the upper bound is quite safe.
Since the norm of the interpolation operator for Chebyshev nodes gives Λ M = O(log(M )), the divergence behavior of e D M t w.r.t. M in the first interval [0, τ ] is not totally surprising and it is quite natural to wonder if there is a connection between the two phenomena. Indeed, even if Faber's Theorem states that Chebyshev interpolation points could not ensure the convergence for all functions, quoting Trefethen [61], the "polynomial interpolation in Chebyshev points is a powerful and reliable method for approximation of functions", since "for Lipschitz continuous functions or better, as is easily done in almost any application, Faber's Theorem ceases to be applicable". By extending this remark to e D M t , we conjecture that its divergence w.r.t. M in the transient phase is due to a "non-generic" set of vectors Ψ, which are constructed in a particularly nasty way. The left panel of Figure 4 makes evident that also random data do not capture the "bad" data, whereas the right panel of Figure 4 suggests that e D M t diverges as log(M ).
Having in mind Proposition 3 and the simulations in Figure 3 & Figure 4, we come to the following conjecture: there existτ ∈ (τ, 1.5τ ) and positive constants a, K independent of M such that for M → ∞.
But in the first term of the right-hand side of (39), the vector Ψ derives from the discretization of a function ψ, namely Ψ = ψ(Θ M ), and ψ is at least Lipschitz continuous with ψ(0) = 0. May we expect to attain the space stability or, even better, its convergence w.r.t. M in the transient phase? In case of a positive answer, is the asymptotic convergence rate w.r.t. M faster for smoother functions? In other words, we wonder if e D M t R M ψ retains memory of the interpolants. Indeed we have that where D is defined in (24) and The formula above makes clear the role of interpolation through P 0 M R M , and that the action of every term corresponds to interpolation with null value at θ = 0 and derivation, imposing, to some extent, that the value at θ = 0 of next derivatives is zero.
So let us get another actor on stage to give us some suggestions, namely the operator H 0, From (15), (43) admits the following representation Notice that, contrary to {e D M t } t≥0 , the family {H 0,M (t)} t≥0 does not define a semigroup on Y M . Moreover H 0,M (t) = 0 for t ≥ τ, while e D M t is not. From (44) we easily get the following bound and sup t≥0 C M (t) ≤ (1 + 2 π log(M + 1)). In the right panel of Figure 5, we notice that the behavior of e D M t resembles the behavior of the function C M (t), and that, moreover, sup t≥0 C M (t) "dominates" sup t≥0 e D M t . So, encouraged by this experiment, we also conjecture that H 0,M (t)Ψ can predict, to some extent, the behavior of e D M t Ψ in the transient phase. To support this we first analyze H 0,M (t)Ψ for Ψ = ψ(Θ M ) by using wellknown interpolation results. Second we experimentally explore the behavior of e D M t ψ(Θ M ) by selecting some test functions ψ with different regularity properties.
Suppose that ψ and its derivatives through ψ (k−1) are absolutely continuous on [−τ, 0] and that ψ (k) has bounded variation V (ψ (k) ), for k ≥ 1. Then for any M > k we have the interpolation process does not converge, but the interpolation error is uniformly bounded in supremum norm, i.e., lim M →∞ for a suitable constant c 1 [60,Chapter 9], and therefore we can conclude that . This effect is due to the well-known Gibbs phenomenon for step functions: the interpolants always "overshoot" their target and, as M increases, the overshooting does not get lower, although the region of overshooting gets narrower. This phenomenon prevents the norm-convergence of the interpolants and is responsible for the uniform bound. Although we will not go further into details here, we note that the Gibbs phenomenon still ensures the pointwise convergence of the interpolants of H for all θ ∈ [−τ, 0].
For our simulations we choose some functions ψ with different regularity properties and we estimate the behavior of the error e D M t Ψ − R M H 0 (t)Ψ with Ψ = ψ(Θ M ) at some selected points t. All the results in the right panels of Figures 7,8,9 & 10 indicate that the convergence occurs in the first interval (t = 0.5) with the rate predicted by the interpolation results for the respective functions (see their plots in the left panels), whereas the order increases in the next intervals (t > 1). Moreover Figure 10 presents the spectral convergence behavior. Finally Figure 6, Figure 11 & Figure 12 show there is no convergence in first interval, in accordance with the interpolation results, but the uniform boundedness holds in all cases.  The numerical experiments indicate that, in the transient phase, e D M t Ψ / ψ is uniformly bounded w.r.t. M when applied to vectors Ψ = ψ(Θ M ) deriving from the discretization of functions ψ, for which the interpolation error is O(M −k ) with k ≥ 0, and it is converging for sufficiently smooth ψ. This supports the effectiveness of the PSD in our context and in particular for the analysis of periodic solutions. Now we consider the second term of both (38) and (39). Let z : R + → R be a continuous function. Then the function w : R + → Y , defined as   Figure 6 for ψ(θ) = e θ 2 − 1. Note that the interpolating polynomial is indistinguishable from the function (left) and that ψ (2) has bounded variation. The convergence rate for t = 0.5 is O(M −2 ) (right).   Figure 10. Same as Figure 6 for ψ(θ) = e −1/θ 2 . Note that the interpolating polynomial is indistinguishable from the function (left) and the spectral convergence for t = 0.5, 1.5, 2.5 (right).  Figure 11. Same as Figure 6 for ψ(θ) = sin 1 θ , if θ < 0, ψ(0) = 0. Note there is no convergence for t = 0.5.
So, to guarantee the convergence of we conjecture that M . The generator of the discrete semigroup is a finite dimensional operator that can be described in terms of a matrix D M , which has been widely studied in the context of PDEs.
In Section 4 we have collected the most relevant results, especially about the spectrum of the matrix D M , of the generator A 0,M , and about their relation with the continuous counterpart A * 0 . We also specified the convergence of the resolvent operators by adapting the proof in [11,14].
About the convergence of the discrete semigroup T 0,M (t) to the continuous semigroup T * 0 (t), very few results are available in the supremum norm, which is the norm of the natural state space of DDEs. In the absence of theoretical results, we have presented a series of numerical tests and compiled some conjectures about the behavior of the semigroup, both in time and as M → ∞.  Figure 14. Same as Figure 13 for the function z built so that it is continuously differentiable with jumps in the second derivative at t = 0.5, 1, 1.5.
To arrive to a rigorous proof, we could start from the representation (42) and consider the properties of the Chebyshev polynomials [55]. Another possibility for proving convergence of the discrete semigroup T 0,M (t) to the continuous semigroup T 0 (t) is to use the convergence results of Theorem 4.1 and to recall that the resolvent operator can be expressed, for Reλ > 0, in terms of the solution operator via Laplace transform as (λI − A * 0 ) −1 (α, ψ) = ∞ 0 e −λs T * 0 (s)(α, ψ) ds, λ / ∈ σ(A * 0 ), We are currently working in this direction. The study of the discrete and continuous shift semigroups is motivated by the PSD of DDEs in * -formulation, which is introduced in this manuscript for the first time. The * -formalism, despite introducing some "complications" due to the embedding into larger spaces, appears to be useful in the context of numerical methods for several reasons. The main advantage is to provide us with the variation-of-constants formula (8), which is particularly useful to relate the nonlinear solution semigroup to the trivial linear semigroup T * 0 (t). This is particularly useful in the prospect of comparing the discrete and continuous nonlinear semigroups. The convergence of the nonlinear semigroups is the next step for the final goal of understanding in which sense the PSD approximates the periodic solutions of (2) and their stability. This is ongoing research of the authors and collaborators announced in the paper [8].
The * framework exploits naturally the pairing between dual spaces and therefore comes with two topologies, the norm topology and the weak* topology. So even if approximations do not necessarily converge with respect to the norm, they still may converge in a weaker sense which, probably, is sufficient for our purposes. We intend to investigate this idea in the near future.
In the case of DDEs, the * -formulation allowed us to give a direct description of the PSD of the operators by means of projections into the discretization space and injections into the bigger space. Moreover, it provides us with a variation-ofconstants formula, which is helpful to relate the nonlinear semigroup to the trivial shift semigroup. We mention also that the * -framework is particularly important also in the PSD of REs and REs/DDEs, where the domain of the infinitesimal generator is described by means of an algebraic, rather than differential, equation. In this context, the * -formulation may allow to overcome the inversion of the nonlinear algebraic equation proposed in [8], treating the condition more efficiently. This is currently under investigation by the authors.
Finally, the PSD approach has been proposed also for DDEs, REs, and REs/DDEs with infinite delay in [37]. Moreover, the analogous theory of * -calculus for infinite delay has been developed in [24]. Similarly to the case with finite delay, the PSD of the shift operator in this context is described by a matrix derived by Laguerre differentiation matrices. Therefore, similar questions to the case of finite delay arise: does the discrete shift operator converge as M → ∞? Is it possible to describe the spectral values, and do they converge to the exact ones? How does the choice of the nodes influence the convergence of the operators and their spectra?
The theoretical aspects are still work-in-progress, and there are still many questions to address. The results and the numerical experiments are encouraging.