Piecewise orthogonal collocation for computing periodic solutions of renewal equations

We extend the use of piecewise orthogonal collocation to computing periodic solutions of renewal equations, which are particularly important in modeling population dynamics. We prove convergence through a rigorous error analysis. Finally, we show some numerical experiments confirming the theoretical results and a couple of applications in view of bifurcation analysis.


Introduction
Including delays within models is often a sound way to describe the relevant phenomena more realistically.Indeed, in various fields of science -such as population dynamics or epidemiology -delays between a cause and the corresponding effects appear rather naturally, which brings the need to resort to delay equations in order to capture adequately the dependence on the past [29,33].
In many applications, the main interest is towards the dynamical analysis of the concerned models, including the computation of invariant sets (such as equilibria and periodic solutions) and the study of their asymptotic stability.Regarding periodic solutions, the piecewise orthogonal collocation method to compute those of Renewal Equations (REs) has been first applied in [13], although the method is not even described therein.Indeed, a formal description appeared first in [4] and [6] where its validity is only shown by means of some numerical experiments.
The aim of the present paper is to give a detailed illustration of the method and provide a rigorous convergence analysis, having in mind the vast presence of REs in the field of population dynamics [12,20,23,26].The convergence analysis follows the main ideas used in [7] for Retarded Functional Differential Equations (RFDEs), which is in turn based on the abstract approach in [28].[7] represents indeed the solution of the long-standing problem of lack of a rigorous proof of convergence for RFDEs.The contribution of the present paper consists in extending the piecewise collocation and its convergence to REs, by tackling the nontrivial challenges due to the inherent differences concerning the class of equations and the relevant spaces of functions.While some of such challenges are mostly technical, the more consistent ones bring the necessity to resort to the theory of resolvents for Volterra Integral Equations (VIEs, see Appendix A).
We conclude this introduction with Subsection 1.1, where we describe the equations of interest and the standard way to formulate the problem of computing periodic solutions as a Boundary Value Problem (BVP).The rest of the paper is divided into three main sections.Section 2 describes the piecewise orthogonal collocation method for REs.Section 3 deals with the theoretical convergence of the method by illustrating how the relevant analysis can be based on the abstract approach in [28].An important part of the proof is collected in Appendix A, in order to not interrupt the reading flow.Finally, Section 4 shows the results of some numerical experiments on REs from population dynamics, also in view of bifurcation analysis.Python demos are freely available at http://cdlab.uniud.it/software.

Renewal equations and boundary value problems
In its most general form, an RE can be written as where, for a positive integer d, F : X → R d is autonomous, in general nonlinear and the state space X is a set of functions from [−τ, 0] to R d for some τ > 0, called the delay.The state of the dynamical system on X at time t associated to (1.1) is denoted by x t , defined as x t (θ) := x(t + θ), θ ∈ [−τ, 0].In particular, we develop our analysis for X := B ∞ ([−τ, 0], R d ), where B ∞ denotes bounded and measurable functions.The elements of B ∞ are considered as functions, not as classes of functions that are equal almost everywhere.Such choice, instead of the classical L1 ([−τ, 0], R d ) [17], is justified by the need of evaluating the functions pointwise in order to deal with collocation.A (non-constant) periodic solution of (1.1) with period ω > 0 1 , if there is any, can be obtained by solving a BVP where the solutions are considered over just one period and the periodicity is imposed to the solution values at the extrema of the period.
Note that this requires to evaluate x at points that fall off the interval [0, ω], due to the delay.In order to deal with this issue, one can exploit the periodicity to bring back the evaluation to the domain [0, ω], which, assuming ω ≥ τ2 , means defining the function x t ∈ X as x t (θ) := x(t + θ), t + θ ∈ [0, ω], x(t + θ + ω), t + θ ∈ [−ω, 0). (1. 2) The periodic BVP can then be formulated as where p is a scalar (usually linear) function defining a so-called phase condition, necessary to remove translational invariance.[18].This is the most natural BVP formulation of the problem following the case of RFDEs, e.g., [9,10,11,19,25,27,30], as well as the most convenient to consider when developing the numerical method.However, we will introduce in Section 3 a slightly different, yet equivalent, formulation in view of the convergence analysis based on [28].
In realistic models, such as those describing structured populations, F has usually the form for some integration kernel K : for some integration kernel k : [0, τ] → R d and some function f : R d → R d .The analysis that follows focuses on (1.4), even though its validity can be also extended to (1.5) (see Remark 3.12 in Section 3).As anticipated previously in the present section, the fundamental differences between the REs we considered and the RFDEs are not limited to their roles in mathematical modeling.Indeed, the theory of REs is typically developed in broader spaces such as L 1 [29,16], resulting, in principle, in a lower degree of differentiability of the relevant solution.

Piecewise orthogonal collocation
This section describes the numerical method used to compute periodic solutions of (1.1), starting from a general right-hand side F which, for the moment, is assumed to be computable without resorting to any further numerical approximations.
Since the period ω is unknown, it is numerically convenient (see, e.g., [19]) to reformulate (1.3) through the map s ω : R → R defined by s ω (t the solution of which is intended to be defined in [0, 1]. (2.1) can be solved numerically through piecewise orthogonal collocation [19].This is particularly useful when adaptive meshes, to better follow the solution profile, might be needed, and is now a standard approach (originally developed for ODEs [8], see MatCont [2]).It means looking for an m-degree piecewise continuous polynomial u in [0, 1] and w ∈ R that satisfy the following system having dimension (1 + Lm) × d + 1: The unknowns are, other than w, those of the form u i,j := u(t i,j ) for (i, j) = (1, 0) and i ∈ {1, . . ., L}, j ∈ {1, . . ., m}3 .
Remark 2.1.Typically, periodic solutions are computed within a continuation framework.This provides a reasonable choice for the phase condition, necessary to ensure the actual wellposedeness of (2.1).
As mentioned at the end of Subsection 1.1, in applications from population dynamics right-hand sides usually feature an integral, therefore cannot be exactly computed in general.This is also the case of (1.4), which reads, once the time has been rescaled, where now Observe that, although the corresponding natural state space is in fact a Banach space of functions defined in [−τ/ω, 0], one could choose spaces of functions defined in [−r, 0] for any 1 ≥ r ≥ τ/ω4 .
Assuming that K can be exactly computed, which is usually the case in applications, the approximation of (2.2) through quadrature as where M is a given approximation level and −τ/ω ≤ η 0 < • • • < η M ≤ 0, can also be exactly computed.Such an approximation corresponds to the secondary discretization introduced in Subsection 3.2 and used in the convergence analysis that follows.The nodes η 0 , . . ., η M and the corresponding weights w 0 , . . ., w M are meant to define a suitable quadrature formula by exploiting possible irregularities in K, meaning that their choice does not need to be made a priori.Moreover, note that the quadrature nodes vary together with ω, since the latter is unknown.In particular, they are completely independent of the collocation nodes mentioned earlier.

Convergence analysis
This section describes the theoretical convergence analysis of the numerical method described in Section 2, following the abstract approach [28].In particular, the convergence analysis that follows applies to the Finite Element Method (FEM), which consists in letting L → ∞ while keeping m fixed and is the classical approach considered in practical applications (e.g., in MatCont [2] or some versions of DDE-Biftool [1,32]).A few comments on the Spectral Element Method (SEM, m → ∞ while keeping L fixed) will follow in Subsection 4.1.
Following the approach for RFDEs in [7], we first reformulate (2.1) as i.e., by imposing the periodicity condition to the states at the extrema of the period rather than to the solution values.In this case the solution x is intended as a map defined in [−1, 1] and there is no need to resort to (1.2).
Although formulations (2.1) and (3.1) are formally different, they are mathematically equivalent and also lead to fundamentally equivalent numerical methods.Indeed, when applying the numerical method described in Section 2 to the problem (3.1) one just introduces redundant variables 5 .
The second step consists in observing that (3.1) fits into the general form of the BVP addressed in [28]: Here the relevant solution v := G(u, α) is assumed to lie in a normed space In the sequel we always use the superscript * to denote quantities relevant to fixed points (which correspond to non-constant periodic solutions).
It follows that (3.1) leads to an instance of (3.4) with G, F and B given respectively by The boundary operator is linear and independent of ω.
The fact that (3.1) can be rewritten as a PAF does not imply that the the convergence framework in [28] can be applied either way.In fact, several assumptions are required.These include theoretical assmuptions, the validity of which depends on the choices of the spaces, as well as on the regularity of the integrand K in the right-hand side (1.4).Subsection 3.1 includes the definitions of such assumptions and their statements as propositions, instanced according to the problems of interest.The other assumptions required concern instead the reduction of the problem to a finite-dimensional one, and will be dealt with similarly in Subsection 3.2.Concerning the proofs of such propositions, we will go through the main points, focusing on the differences with respect to the analogous propositions in the RFDE case [7], with the exception of a more complicated one to which we dedicate the entire Appendix A.

Theoretical assumptions
The hypotheses on the original problem needed to prove the validity of the theoretical assumptions in [28] are: which are measurable with respect to both their arguments; (T4) the map x → D 2 K(r, x) is piecewise continuous for all r ∈ R; (T5) there exist r > 0 and κ ≥ 0 such that for every (v t , ω) ∈ b((v * t , ω * ), r) 7 , uniformly with respect to t ∈ [0, 1].
The first theoretical assumption (AFB, [28, page 534]) concerns the Fréchet-differentiability of the operators F and B appearing in (3.4).Since p is linear, so is B in (3.7), hence the latter is Fréchet-differentiable.The validity of the assumption is thus a direct consequence of the following.
Proposition 3.1.Under (T1), (T2) and (T3), F in (3.6) is Fréchet-differentiable, from the right with respect to ω, at every and Proof.The proof is rather technical and goes as that of [7, Proposition 2.1], therefore we avoid to repeat all the steps for brevity and to better concentrate on the differences.Basically, the expression (3.8), defined through (3.9) and (3.10), is directly proven to satisfy the definition of differentiable function according to [3, Definition 1.1.1].It is worth pointing out that assuming an integral right-hand side such as (1.4), which is anyway typical in applications from population dynamics, is crucial for this proposition in the case of REs.Basically, for the thesis to hold it is required that the right-hand side always lies in a more regular space than U (which is always the case for RFDEs, where U plays the role of the space of the derivatives).This can be observed by looking at the last addend of [7, (2.12)], where the derivative of the state of an element of V appears as a factor.Without any assumption whatsoever on F the same would happen in Proposition 3.1, which would be a problem since V is as regular as U in the present case.
The second theoretical assumption (AG, [28, page 534]) concerns the boundedness of G defined in (3.5).The following proposition concerns its validity, and its proof is an immediate consequence of the definition (3.5).

Proposition 3.2. Under (T2), G is bounded.
The third theoretical assumption (Ax * 1, [28, page 536]) concerns the local Lipschitz continuity of the Fréchet derivative of the fixed point operator Φ in (3.4) at the relevant fixed points.In the sequel (u * , α * , ω * ) ∈ U × A × (0, +∞) is a fixed point of Φ and x * is the corresponding 1-periodic solution of (1.1).With respect to the validity of Assumption Ax * 1 the following holds.Proposition 3.3.Under (T1), (T2), (T3) and (T5), there exist r ∈ (0 Proof.Just as its RFDE counterpart in [7], the proposition can be proved thanks to the fact that u * lies in fact in a more regular subspace of its space U, which is a consequence of the assumption (1.4).
The fourth (and last) theoretical assumption (Ax * 2, [28, page 536]), concerns the well-posedness of a linearized inhomogeneous version of the PAF (3.3).Its validity can be proved under (T1), (T2), (T3) and (T4), together with an additional requirement, which in turn follows from assuming, e.g., the hyperbolicity of the periodic solution 8 of the original problem.It is convenient to introduce the abbreviations Proposition 3.4.Under (T1), (T2), (T3) and (T4), let T * (t, s) : X → X be the evolution operator for the linear homogeneous RE (3.12) Proof.(3.12) can be treated as an initial value problem for v = G(u, α), i.e., for t ∈ [0, 1], imposing then the boundary conditions in (3.12).We can write where, in turn, The first boundary condition in (3.12) gives then

Numerical assumptions
As anticipated, the present subsection deals with the numerical assumptions, which concern the chosen discretization scheme for the numerical method.Such scheme is defined by the primary and the secondary discretizations.
As in [7], the primary discretization consists in reducing the spaces U and A to finite-dimensional spaces U L and A L , given a level of discretization L. This happens by means of restriction operators ρ All of them are linear and bounded.In the following we describe the specific choices we make in this context, based on piecewise polynomial interpolation.
Starting from U, which concerns the interval [0, 1], we choose the uniform outer mesh and inner meshes where 0 < c 1 < • • • < c m < 1 are given abscissae for m a positive integer.Correspondingly, we define whose elements u L are indexed as with components in R d .Finally, we define, for u ∈ U, and, for u L ∈ U L , π + L u L ∈ U as the unique element of the space such that Above Π m is the space of R d -valued polynomials having degree m and, when needed, we represent p ∈ Π + L,m through its pieces as where, for ease of notation, we implicitly set and {ℓ m,i,0 , ℓ m,i,1 , . . ., ℓ m,i,m } is the Lagrange basis relevant to the nodes {t + i,0 } ∪ Ω + L,i .Observe that the latter is invariant with respect to i as long as we fix the abscissae c j , j ∈ {1, . . ., m}, defining the inner meshes (3.16).Indeed, for every i ∈ {1, . . ., L}, where {ℓ m,0 , ℓ m,1 , . . ., ℓ m,m } is the Lagrange basis in [0, 1] relevant to the abscissae c 0 , c 1 , . . ., c m with c 0 := 0. Similarly, for A, which concerns the interval [−1, 0], we choose and Correspondingly, we define with indexing and, for α L ∈ A L , π − L α L ∈ A as the unique element of the space Remark 3.5.It is worth pointing out that more general choices can be made concerning outer and inner meshes.In particular, as already remarked in Section 2, in practical applications adaptive outer meshes represent a standard for RFDEs, see, e.g., [19].As for inner meshes, abscissae including the extrema of [0, 1] can also be considered, paying attention to put the correct constraints at the internal outer nodes, i.e., t ± i for i ∈ {1, . . ., L − 1}.
The secondary discretization consists in replacing F in the first of (3.4) with an operator F M that can be exactly computed, for a given level of discretization M. In particular, we define F M through an approximated version F M of the right-hand side F defined in (1.4) as where Indeed, in realistic applications the integrand function in (1.4) can be exactly computed, as already remarked at the end of Section 2. Correspondingly, Φ M is the operator obtained by replacing F in Φ in (3.4) with its approximated version, i.e., Φ M : (3.32) A secondary discretization for G in (3.4) is instead unnecessary, since it can be evaluated exactly in π + L U L × π − L A L according to (3.17) and (3.26).As for the operator p defining the phase condition in (3.4), we assume that it can be evaluated exactly in π + L U L . 9rom the two discretizations together we can define the discrete version ) of Φ L,M can be found by standard solvers for nonlinear systems of algebraic equations and, as will be shown in Subsection 3.3, its prolongation P L (u * L,M , α * L,M , ω * L,M ) is then considered as an approximation of a fixed point ) is considered as an approximation of the solution v * = G(u * , α * ) of (3.4).
The hypotheses on the discretization method needed to prove the validity of the numerical assumptions in [28]  (N3) the nodes η 0 , . . ., η M , together with the weights w 0 , . . ., w M chosen for the secondary discretization as in (3.31) define an interpolatory quadrature formula which is convergent in The first numerical assumption to be verified in [28] is Assumption AF K B K (page 535).As already observed, B and p are linear functions, thus the proof of its validity is a direct consequence of the following.Proposition 3.6.Under (T1), (T2) and (T3), F M is Fréchet-differentiable, from the right with respect to ω, at every ( v, û, ω) ∈ V × U × (0, +∞) and and Proof.The proposition can be proved as Proposition 3.1, by replacing F in the first of (3.4) with F M in (3.31).
Proof.The proof is substantially the same as that of [7, Proposition 3.7] since the same primary discretization is used; the main difference lies in the spaces involved, and therefore in the norm in which the left-hand side needs to be evaluated.In practice, we can define κ , where Λ m is the Lebesgue constant of the chosen nodes.In particular, unlike the RFDE case, the Lebesgue constant defined by the derivatives of the Lagrange polynomials is not involved.Correspondigly, the last numerical assumption (CS2, page 537), can be seen as the discrete version of Ax * 2 therein, here Proposition 3.4.With respect to its validity, the following holds.Proposition 3.8.Under (T1), (T2), (T3), (T4), (T5), (N1), (N2) and (N3), the operator DΨ L,M (u * , α * , ω * ) is invertible and its inverse is uniformly bounded with respect to both L and M.Moreover, Proof.The proof of this proposition is a bit laborious, just like its counterpart in the RFDE case.The latter has been proved in [7] in several steps, the first of which concerns the invertibility of the operator DΨ L,M (u * , α * , ω * ) defined in (3.33) for L, M large enough and can be proved as [7,Proposition 3.11].The second step concerns the uniform boundedness of DΨ −1 L,M (u * , α * , ω * ) and follows the ideas of [7, Lemma 3.12].The latter is based on [7, Proposition A.8], which states that lim L,M→∞ ω L,M = ω, where ω L,M is the last component of the solution of the discretized version of (3.12).The limit does not necessarily hold in the present case, however it can be proved that |ω L,M − ω| is uniformly bounded.This follows by the fact that ∥ξ * L,M,2 − ξ * 2 ∥ X is in turn uniformly bounded (but not necessarily vanishing) thanks to the choice of U in (T1), where ξ * L,M,2 is the discrete version of ξ * 2 .As a consequence, the error component called ε ω,L,M in [7, Lemma 3.12] cannot be proven to vanish in the present case, but would still be uniformly bounded, and that is enough to complete the second step of the proof.The third and last step consists in proving that Ψ L,M (u * , α * , ω * ) vanishes and goes as the proof of [7, Proposition 3.13].

Final convergence results
From the propositions in the previous subsections we can conclude that our problem of finding a fixed point of Φ in (3.4) satisfies all the assumptions required by [28] under certain hypotheses on the state spaces, the discretization and the regularity of the right-hand side.As a consequence, the relevant FEM converges.Theorem 3.9 ([28, Theorem 2, page 539]).Under (T1), (T2), (T3), (T4), (T5), (N1), (N2) and (N3), there exists a positive integer N such that, for all L, M ≥ N, the operator R L Φ M P L has a fixed point (u ) and v * = G(u * , α * ).Thanks to Proposition 3.8, the error on (v * , ω * ) is determined by the last factor, namely the consistency error.For the latter, thanks to basic results on polynomial interpolation, we can write where Λ m is the Lebesgue constant associated to the collocation nodes and the terms and are called respectively primary and secondary consistency errors.
As for ε L in (3.34), which concerns only the primary discretization, a bound can be obtained from the regularity of u * through the following theorem.Theorem 3.10.Let K ∈ C p (R × R d , R d ) for some integer p ≥ 0.Then, Under (T1), (T2), (N1) and (N2), it holds that u is a periodic solution of (1.1) modulo rescaling of time, and it is bounded by (T2).Thus, if K is continuous, then the periodic extension of v * is continuous in [0, +∞), thus also in [−1, ∞) by periodicity.Using the continuity of K again, we obtain that we immediately have also The whole reasoning can be iterated, proving the first part of the result.
To prove (3.36), we observe first that On the other hand, ε M in (3.35) concerns only the secondary discretization and is therefore absent whenever the latter is not needed.However, concerning our specific problem, according to (3.4) and (3.32), it can be written as and needs to be considered if the integral in (1.4) cannot be exactly computed, in which case (3.39) is basically a quadrature error.Assuming that M varies proportionally to L, one can choose a formula that guarantees at least the same order of the primary consistency error, so that the order of convergence of the final error is in fact the one given by theorem 3.10.
Remark 3.11.In principle, one could discretize the problem by choosing, for each mesh interval, a set of representation nodes used to interpolate which are independent from the collocation nodes.That would mean that the unknowns of the discrete problem are given by the values of the relevant functions at the representation nodes, while the equations need to be satisfied at the collocation nodes.If x r L,M is the vector of the unknowns and Q L : X L → X is the prolongation operator corresponding to the representation nodes (while P L , R L refer to the collocation ones), the problem actually reads R L Q L x r L,M = R L Φ M Q L x r L,M .Thus, the vector x * L,M given by the values of the relevant function at the collocation nodes is the solution of the discrete fixed point problem, in fact, Remark 3.12.The entire convergence analysis can as well be carried out for right-hand sides of the form (1.5).In this case, the different theoretical and numerical assumptions read (T5) there exist r > 0 and κ ≥ 0 such that Moreover, the above can be easily further generalized to the case (3.40)

Results
This section deals with the numerical computation of periodic solutions of some specific REs from the field of population dynamics.In particular, in Subsection 4.1 we will provide experimental proof of the order of convergence (3.36), while in Subsection 4.2 we will show, by means of two examples, how bifurcations may be detected thanks to the computed periodic solutions.
Starting from the exact solution at γ = 4, the branch of periodic orbits is continued up to the first period doubling after the Hopf bifurcation.The continuation is performed using a trivial phase condition by forcing x(0) = σ, and Chebyshev extrema as collocation points.The left panel of Figure 4.1 confirms the O(h m+1 ) behavior (being p = +∞).Given the experimental proof, found in [19], that in the case of RFDEs the order of convergence increases from m to m + 1 when using Gauss-Legendre collocation points, we wonder whether a similar increase can be observed in the case of REs.The above experiment is then replicated using such collocation points, and the right panel of Figure 4.1 confirms the same O(h m+1 ) behavior as the case of Chebyshev extrema.Next, we show that we obtain the same results with an RE that is defined by a right-hand side of the form (3.40) as where α := 1/ 1 0 s 2 e −10s ds = 500e 10 /(e 10 − 61), while γ is the varying parameter.As shown in [31], a Hopf bifurcation occurs when log γ ≈ 1.6553.
Starting from a perturbation of the equilibrium at the Hopf bifurcation point, the branch of periodic orbits is continued up to log γ = 1.75.Given the absence of an exact expression of the true solution, unlike the case (4.1), the error is computed with respect to a reference solution which is in turn computed using L = 1000 and m = 4. Figure 4.2 confirms again the O(h m+1 ) behavior regardless of the chosen collocation nodes.
Concerning the convergence of the SEM, it is not yet clear whether it can also be proved under the general framework used in the current work (see [7,Subsection 4.4] for a brief discussion concerning the RFDE case).However, some numerical experiments run by the authors suggest that the SEM does converge for periodic BVPs defined by REs.In the case of (4.2), figure 4.3 shows a spectral decay of the error as m increases while L = 1 remains fixed, although some numerical instability can be observed when using large (> 35) values of m.
Finally, consider the RE given by where Error on the periodic solution of (4.2) at log γ = 1.75.Left: m = 3 (dashed line) and m = 4 (solid line) using Chebyshev points, compared to straight lines having slope 4 (dashed) and 5 (solid).Right: m = 3 (dashed line) and m = 4 (solid line) using Gauss-Legendre points, compared to straight lines having slope 4 (dashed) and 5 (solid).Observe that the function h is in C 2 (R, R) but not in C 3 (R, R), due to a discontinuity of the third derivative at x = 1/2.Using the method in [31], it can be shown that a Hopf bifurcation occurs at γ ≈ 3.4031.Starting from a perturbation of the equilibrium at the Hopf bifurcation point, the branch of periodic orbits is continued up to γ = 6.The integrals representing the distributed delays are again approximated through Clenshaw-Curtis quadrature, and the experiment is again performed using both Chebyshev and Gauss-Legendre collocation points.The error is computed with respect to a reference solution which is in turn computed using L = 400 and m = 4.
Being the values for m considered (3 and 4) greater than p = 2, the left panel of Figure 4.4 confirms the O(h p+1 ) behavior in the former case, as the right one does in the latter.

Applications
As anticipated, the examples in this section show how the method can be used within a continuation framework in view of a bifurcation analysis.The next RE that we consider to this aim is As shown in [31], a period doubling bifurcation occurs when log γ ≈ 3.8777.In [15] a Floquet theory for REs has been developed and proved valid; in particular, the stability of a periodic solution is determined by the eigenvalues of the monodromy operator of the corresponding linearized equation, according to [15,Corollary 15].Such eigenvalues can, in turn, be computed numerically thanks to the pseudospectral method in [14], which is used in order to obtain Figure 4.5.As shown therein, indeed, two stable periodic solutions can be computed on opposite sides of log γ = 3.8777, one having roughly double minimal period than the other, thus confirming the presence of a period doubling bifurcation.
Period doubling bifurcations also occur in the case (4.1), as shown in [13].In  particular, the second one after the Hopf bifurcation (near which it is not possible to obtain an exact expression of the solution) is detected at γ ≈ 4.497.While, at that point, new stable periodic solutions emerge having double minimal period than the stable old ones, unstable solutions with (roughly) unchanged period also exist, and can be computed using the method that we are proposing.Figure 4.6, indeed, shows that two periodic solutions which are very close to each other can be computed on opposite sides of γ = 4.497, one being stable and the other being unstable, thus confirming the presence of a bifurcation.

Concluding remarks
In the past few decades, piecewise orthogonal collocation has been largely used to compute periodic solutions of various classes of delay equations.The present work aims at giving a rigorous description of such a method in the case of REs, while furnishing a complete theoretical error analysis as well as experimental proofs of its validity.In particular, the theoretical proof is based on the abstract approach in [28] for general BVPs, as it is the case for the work [7] concerning the corresponding proof for RFDEs.The main result concerns the FEM and states that, for smooth problems, the error decays as O(L −(m+1) ) where L is the increasing number of mesh intervals, while m is the constant degree of the piecewise polynomials used.
Given both the convergence analysis for RFDEs in [7] and that in the present paper for REs, we expect to be able to rely on the general approach in [28] also for the case of coupled RFDE/RE systems, motivated by their predominance in population dynamics.The proof is currently work in progress and the expected main challenge is given by the need to resort to the theory of VIEs with measure kernels [21,Chapter 10].
Moreover, it would be interesting to extend the method (and, correspondingly, the convergence analysis) to different and more complex classes of delay equations.The first step that we plan to take in this direction is to try to apply the approach to differential equations with non-constant delays (in particular, state-dependent delays), for which the setting defined in [7] cannot be applied.
T * (1, 0) with respect to (A.1), we get In order to obtain a contradiction with (A.2), we resort to adjoint theory for VIEs (see [21] as a general reference).Indeed, any RE of the form 11 and r ≤ 1 can be written as the VIE where and y s 0 = ψ for some ψ ∈ Y, where we use the notation y s (η) := y(s + η) for η ∈ [0, r].Then, one defines the backward evolution family {V(s, s 0 )} s≤s 0 on Y through V(s, s 0 )ψ = y s .
From the theory of resolvents [21, Chapter 9, Section 3], we can express the solution of (A.3) and that of (A.4) respectively as where R * 0 is the resolvent of (A.3).Given t ∈ R, consider now the pairing [ Observe that such bilinear form is nondegenerate for all t ∈ R whenever K * (and thus K * 0 ) is nontrivial.Indeed, assume by contradiction that there exists ψ ∈ Y such that ψ is nonzero, but [ψ, •] t is constantly 0. By the nondegenerateness of the bilinear form ⟨•, •⟩, this means that the innermost integral is 0 for all α ∈ X and almost all η ∈ [0, 1].If α := x t , where x is the (unique modulo multiplication by constant) 1-periodic solution of the VIE, then such integral is equal to Thus x t+1 is almost everywhere equal to 0. Using periodicity, this means that x is almost everywhere 0, which is only possible if K * is trivial, contradiction.Using similar arguments, one can prove that there is no nonzero α ∈ X such that [•, α] t is constantly zero, after exchanging the order of integration in the definition of [•, •] t .
We claim that the forward monodromy operator and the corresponding backward one are adjoint w.r.t.(A.7), i.e., that (A.8) Indeed, using (A.6), we have As for A, we have where the first equality comes from the definition of g, the second follows from the substitution σ ← t + σ, the third is obtained by exchanging the order of integration between η and σ, the fourth follows from the definition of ψ 0 , the fifth is obtained by exchanging the order of integration between η and β, the sixth is due to the fact that K * 0 (t + σ, t − 1 + η) vanishes for η < σ, the seventh follows from the substitution η ← 1 + η, the eighth is obtained by just renaming the variables, the ninth is due to the fact that K * 0 (t + 1 + β, t + σ) vanishes for β > σ, the tenth is obtained by exchanging the order of integration between β and σ, the eleventh follows from the substitution σ ← σ − t and the last comes from the definition of f .As for B, we have where the first equality comes from the definition of g, the second is obtained by exchanging the order of integration between θ and σ, the third is due to the fact that R * 0 (σ, t − 1 + η) vanishes if σ < t − 1 + η and K * 0 (t + θ, σ) vanishes if σ < t − 1 + θ, the sixth is obtained by changing the order of integration, the seventh follows from the 1-periodicity of R * 0 , the eighth follows from the substitution σ ← σ − t, the ninth is due to the fact that R * 0 (1 + β, σ) vanishes if σ > 1 + β and K * 0 (σ, t + θ) vanishes if θ < σ − t − 1, the tenth follows from the substitutions θ ← θ − t and β ← β − t and the last comes from the definition of f .Eventually, using (A.5), we have Under the assumption that 1 is a simple eigenvalue, both the VIE and its adjoint have a unique 1-periodic solution modulo multiplication by some constant, say x and y respectively.Thus, the associated states y t and x t are respectively the left and right 1-eigenvectors of the operator T(t + 1, t).Again thanks to their uniqueness, we have [y t , x t ] t ̸ = 0 for all t ∈ R.Moreover, the continuity of the map t → [y t , x t ] t and the mean value theorem for definite integrals let us conclude that 1 0 [y t , x t ] t dt ̸ = 0. Finally, observe that are: (N1) the primary discretization of the space U is based on the choices (3.15)-(3.21);(N2) the primary discretization of the space A is based on the choices (3.24)-(3.30);

Figure 4 . 5 :
Figure 4.5: Stable periodic solutions of (4.4) at log γ = 3.83 (top), 3.8777 (middle, period doubling) and 3.9 (bottom), having periods ω = 4, 4 and ≈ 8.003 respectively, computed with L = 20, 20 and 40 and m = 5.Left: representation of two periods (top, middle) and one period (bottom) of the solutions in the scaled interval [0, 2].Right: eigenvalues of the corresponding monodromy operator with respect to the unit circle, all internal in the first and third pictures with the exception of the trivial multiplier 1 (due to linearization, [15, Proposition 10]).
and the operator G : U × A → V represents a linear operator which reconstructs the solution given u and an initial state α, also lying in a Banach space A of functions [−1, 0] → R d .β is a vector of parameters living in a Banach space B. The first line of (3.2) represents the concerned functional equation via the function F : V × U × B → U and the second represents the boundary conditions via B : V × U × B → A × B.
R, where R is the range of the operator I X − T * (1, 0) (see Appendix A for a detailed proof of the latter point).