{'fraction_non_alphanumeric': 0.060285172998304916, 'fraction_numerical': 0.01952338219164423, 'mean_word_length': 3.927869496855346, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical computation of periodic solutions of renewal equations from population dynamics

We describe a piecewise collocation method for computing periodic solutions of renewal equations, obtained as an extension of the corresponding method in [K. Engelborghs et al., SIAM J. Sci. Comput., 22 (2001), pp. 1593–1609] for retarded functional differential equations. Then, we rigorously prove its convergence under the abstract framework proposed in [S. Maset, Numer. Math., 133 (2016), pp. 525– 555], as previously done in [A.A. and D.B., SIAM J. Numer. Anal., 58 (2020), pp. 3010–3039] for general retarded functional differential equations. Finally, we show some numerical experiments on models from populations dynamics which confirm the order of convergence obtained theoretically, as well as a few 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 the dependence on the past adequately [28,30].
In many applications, the main interest is towards the dynamical analysis of the models, which includes the computation of invariant sets (such as equilibria and periodic solutions) and the study of their asymptotic stability. Such analysis is harder for delay equations than for Ordinary Differential Equations (ODEs), the main general reason being the infinite dimension of the dynamical systems generated by the former. As a result, concerning the aforementioned aspects of the dynamical analysis, while the theory behind the subject of periodic solutions of Retarded Functional Differential Equations (RFDEs) is completely established [15,21], the current literature on systems involving also Renewal Equations (REs) is nowhere near as developed. For instance, a Floquet theory for REs which allows to study the local asymptotic stability of their periodic solutions through the principle of linearized stability was established only recently in [12]. Concerning the existing relevant computational tools, the pseudospectral method proposed in [9] allows to perform a numerical bifurcation analysis of nonlinear delay equations, and is based on approximating the equations via collocation, meaning to study the resulting ODE system through numerical packages for stability and bifurcation analyses for ODEs (e.g., the widespread continuation-based software package MatCont [2]). This is what makes it an instance of what is called pragmatic-pseudospectral approach in [10], where the validity of the method has been substantiated through its application on a class of REs. However, while the convergence of the method in [9] concerning the approximation of equilibria and their stability has been also proved theoretically, that concerning periodic solutions was only conjectured, both in [9] and [10]. The pragmatic-pseudospectral approach is opposed to the so-called expert-pseudospectral one, first used in [10] for REs, which relies more directly on the Principle of Linearized Stability: it requires to compute (an approximation of) the sought solution and then linearize the system around it in order to investigate its local stability. Some progress has been made nevertheless as far as RFDEs are concerned. With regard to software tools, some already allow to study various aspects of RFDEs defined by discrete delays without any prior manual discretization. The main example is perhaps DDE-Biftool [1], in which the computation of periodic solutions is based on the piecewise orthogonal collocation method already proposed in [17], although a theoretical proof of the convergence of the method for the general form of an RFDE was presented much later and only recently in [6], based on the abstract approach in [27] to compute numerical solutions of general Boundary Value Problems (BVPs, see also [25,26]). An extension of this collocation method to compute periodic solutions of REs has been first applied in [10] in the framework of the expert-pseudospectral approach, although the method is not even described therein. Indeed, a formal description appeared first in [4] and [5] where its validity was only shown by means of some numerical experiments.
The aim of the present paper is to discuss this method for computing periodic solutions of REs in detail for the first time, having in mind their vast presence in the field of population dynamics [8,19,22,24]. Moreover, we prove its validity rigorously by means of a theoretical convergence analysis, inspired by the successful approach of [6] for RFDEs, by tackling the nontrivial challenges due to the inherent differences concerning the class of equations and the relevant spaces of functions.
We conclude this introduction with Subsection 1.1, where we describe the equa-tions of interest and the standard way to formulate the problem of computing periodic solutions as a BVP. The rest of the paper is divided into three main sections. Section 2 describes the piecewise orthogonal collocation method in practice. Section 3 deals with the theoretical convergence of the method by illustrating how the relevant analysis can be based on the abstract approach in [27]. The key points of the proofs of the propositions therein, which feature the most interesting differences with respect to the proof in [6] for RFDEs, are presented in Appendix A. Finally, Section 4 shows the results of some numerical experiments on REs from population dynamics. In particular, those in Subsection 4.1 confirm the order of convergence of the method obtained theoretically in Section 3, while those in Subsection 4.2 are meant to show related bifurcations by using the method within a continuation framework. 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 classically (see, e.g., [14]) defined as X := L 1 ([−τ, 0], 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 A periodic solution of (1.1) with period ω > 0, 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, ω]. Formally, this means defining a periodic state x t ∈ X as 1 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, which is necessary in order to remove translational invariance. An example of phase condition is the trivial one, of the form x k (0) =x for some k ∈ {1, . . . , d}, wherex is fixed. Otherwise, an integral phase condition is of the form for some k ∈ {1, . . . , d} wherex is a given reference solution [16].
In realistic models, such as those describing structured populations, F is usually an integral function, e.g., of the form (1.5) 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.13 in Section 3).

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., [17]) to reformulate (1.3) through the map s ω : R → R defined by the solution of which is intended to be defined in [0, 1] and represents a 1-periodic function.
(2.1) can be solved numerically through polynomial collocation. This would mean looking for an m-degree polynomial u in [0, 1] and w ∈ R such that for given collocation points 0 < θ 1 < · · · < θ m < 1, where u is defined as in (1.2). However, following [17], the method can be improved using piecewise orthogonal collocation. 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, see MatCont [2]). In this case the numerical solution u in [0, 1] is a piecewise continuous polynomial obtained by solving the following system having dimension for a given mesh 0 = t 0 < · · · < t L = 1 and collocation points 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} 2 .
Remark 2.1. Typically, periodic solutions are computed within a continuation framework. In that case, concerning the choice of a suitable phase condition p, either somex or somex, defined as in Section 1, is available: indeed, if the continuation starts close to a Hopf bifurcation (see, e.g., [23, Sections 3.4 and 3.5]), a coordinate of the equilibrium giving rise to it is a natural choice forx. Alternatively, one could use an integral phase condition wherex is a cycle with period 2π/β and β is in turn (the absolute value of) the imaginary part of the conjugate pair determining the Hopf bifurcation. This cycle is intended to represent an approximation of a periodic solution corresponding to a value of the parameter obtained by slightly perturbing the Hopf one, and a reasonable guess for the amplitude is given by √ α, where α is the real part of the aforementioned conjugate pair at the perturbed value of the parameter. On the other hand, at the subsequent continuation steps,x can be defined as a component of the periodic solution computed at the previous continuation step.
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 x t ∈ X := L 1 ([−1, 0], R d ). 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 r ≥ τ ω , thus also for r = 1 (see footnote before (1.2)). This is, in fact, necessary in cases when ω might vary (recall that it is an unknown) while the spaces need to be fixed, as is required for the forthcoming analysis.
Assuming that the integrand K can be exactly computed, which is usually the case in applications, the approximation of (2.2) 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 main ideas of the theoretical convergence analysis of the numerical method described in Section 2, following the abstract approach [27], intended for neutral functional differential equations. To be more precise, two methods can be associated to a piecewise collocation strategy: the Finite Element Method (FEM), which consists in letting L → ∞ while keeping m fixed, and the Spectral Element Method (SEM), consisting in letting m → ∞ while keeping L fixed. The convergence analysis that follows only applies to the former, which is anyway the classical approach considered in practical applications (e.g., in MatCont [2] or some versions of DDE-Biftool [1]).

Remark 3.1.
It is not yet clear whether the convergence of the SEM can also be proved under the general framework used in the current work (see [6,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, although some numerical instability is observed (see Figure 4.3 in Section 4).
The first step consists in reformulating (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 therefore there is no need to resort to the periodic state (1.2) in the formulation.
Although formulations (2.1) and (3.1) are formally different, they 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. However, this reformulation was used in [6] in the RFDE case since it was necessary in order to apply the approach [27]. Therefore, the same formulation will be considered here for uniformity reasons 3 .
The second step consists in observing that (3.1) fits into the general form of the BVP addressed in [27], that is (3.2) and the relevant solution v := G(u, α) lies in a normed space of functions V, u lies in a Banach space of functions U, and the operator G : U × A → V represents a linear operator which reconstructs the solution given a function u and an initial value/state α 4 , also lying in a Banach space A. β is a vector of possible parameters, usually varying together with the solution and living in a Banach space B. The first line of (3.2) represents the functional equation of neutral type via the function F : V × U × B → U, the second line represents the boundary condition via B : V × U × B → A × B. The latter usually includes a proper boundary condition on the solution (the component in A) and a further condition imposing the necessary constraints on the parameters (the component in B). Note that in the case of a neutral differential equation the space U would play the role of the space of the derivative of the sought solution. In the case of an RE, however, no derivatives are involved and U can play the role of the space of the solution.
In [27], (3.2) is then translated into a fixed-point problem, the so-called Problem in Abstract Form (PAF) which consists in finding In the sequel we always use the superscript * to denote quantities relevant to fixed points.
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 [27] 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, in Appendix A we will go through the main points, focusing on the differences with respect to the analogous propositions in the RFDE case [6].

Theoretical assumptions
The hypotheses on the original problem needed to prove the validity of the theoretical assumptions are collected below, where B ∞ denotes bounded and measurable functions.
(T3) K : R × R d → R d is continuous and has partial derivatives D 1 K, D 2 K.
(T5) There exist r > 0 and κ ≥ 0 such that Note that the choice of working with the spaces B ∞ instead of the classical L 1 spaces is justified by the need of evaluating the functions pointwise in order to deal with collocation.
The first theoretical assumption (AFB, [27, page 534]) concerns the differentiability of the operators F and B appearing in (3.4) in the sense of Fréchet. 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. is Fréchet-differentiable, from the right with respect to ω, at every (v,û,ω) ∈ V × U × (0, +∞) and and (3.10) The second theoretical assumption (AG, [27, page 534]) concerns the boundedness of the Green operator G defined in (3.5). The following proposition concerns its validity. The third theoretical assumption (Ax * 1, [27, 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.

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. The former is defined as explained in [6] for the RFDE case, which is reported below for convenience.
The primary discretization consists in reducing the spaces U and A to finitedimensional spaces U L and A L , given a level of discretization L. This happens by means of restriction operators ρ + L : U → U L , ρ − L : A → A L and prolongation operators π + L : U L → U, π − L : A L → A, which extend respectively to and 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 with components in R d . Finally, we define, for u ∈ U, ρ + L u := (u(0), u(t + 1,1 ), . . . , u(t + 1,m ), . . . , u(t + L,1 ) . . . , u(t + L,m )) T ∈ U L (3. 17) and, for u L ∈ U L , π + L u L ∈ U as the unique element of the space 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.14). Indeed, for every i ∈ {1, . . . , L}, and Correspondingly, we define and, for ψ L ∈ A L , π − L ψ L ∈ A as the unique element of the space Elements in Π − L,m are represented in the same way as those of Π + L,m by suitably adapting both (3.20) and (3.21).

Remark 3.6.
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., [17]. 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 (3.29) where − τ ω = α 0 < · · · < α M = 0. 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.15) and (3.24). As for the operator p defining the phase condition in (3.4), we assume that it can be evaluated exactly in π + L U L . This is indeed true in the case of integral phase conditions if the piecewise quadrature is based on the mesh of the primary discretization, which is the standard approach in practical applications.
From the two discretizations together we can define the discrete version can be found by standard solvers for nonlinear systems of algebraic equations and, as will be shown in Subsection 3.3, its pro- The hypotheses on the discretization method needed to prove the validity of the numerical assumptions in [27] are collected below.
(N1) The primary discretization of the space U is based on the choices (3.13)-(3.19).
(N2) The primary discretization of the space A is based on the choices (3.22)-(3.28).
(N3) The nodes α 0 , . . . , α M , together with the weights w 0 , . . . , w M chosen for the secondary discretization as in (3.29) define an interpolatory quadrature formula which is convergent in The first numerical assumption to be verified in [27] is Assumption AF K B K (page 535). As already observed, B and p are linear functions, thus its validity is a direct consequence of the following.
For the remaining numerical assumptions, it is useful to define Ψ L,M : The second numerical assumption in [27] is CS1 (page 536), which is somehow the discrete version of Ax * 1 therein, here Proposition 3.4. With respect to its validity, the following holds.  Correspondigly, the last numerical assumption (CS2, page 537), can be seen as the discrete version of Ax * 2 therein, here Proposition 3.5. With respect to its validity, the following holds. Proposition 3.9. Under (T1), (T2), (T3), (T5), (N1), (N2), (N3) and (N4), the operator DΨ L,M (u * , ψ * , ω * ) is invertible and its inverse is uniformly bounded with respect to both L and M. Moreover, with r 1 and κ as in Proposition 3.8.

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 [27] under certain hypotheses on the state spaces, the discretization and the regularity of the right-hand side. As a consequence, the relevant FEM converges.
Thanks to Proposition 3.9, 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.32), which concerns only the primary discretization, a bound can be obtained from the regularity of u * through the following theorem.
, on the other hand, 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.30), it can be written as and needs to be considered if the integral in (1.4) cannot be exactly computed, in which case (3.35) 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 magnitude of the final error is in fact the one given by theorem 3.11.

Remark 3.12.
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 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.13. 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

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.34), while in Subsection 4.2 we will show, by means of two examples, how bifurcations may be shown thanks to the computed periodic solutions.

Numerical tests
The first RE that we consider is for which, as shown in [10], the exact expression of the periodic solution between a Hopf bifurcation (at γ = 2 + π/2) and the first period doubling (at γ ≈ 4.327) is The integral representing the distributed delay was approximated through a Clenshaw-Curtis quadrature [31] rescaled to the interval [−3, −1]. Starting from the exact solution at γ = 4, the branch of periodic orbits was continued up to the first period doubling after the Hopf bifurcation. The continuation was performed using trivial phase condition withx = σ, and Chebyshev extrema as collocation points. The left plot of Figure 4.1 confirms the O(h m ) behavior. Given the experimental proof, found in [17], that the order of convergence increases when using Gauss-Legendre collocation points, the same behavior can be expected in the case of REs. The above experiment was then replicated using such collocation points, and the right plot of Figure 4.1 confirms the O(h m+1 ) behavior.
Note that a theoretical proof of such convergence property of Gauss-Legendre nodes is still missing (indeed, all theoretical results mentioned in [17] only concern 10 −1. ODEs). Moreover, there is also no proof that the order of convergence cannot be higher than the one obtained in Theorem 3.11 in some cases, as the following example makes clear. The second RE that we consider is defined by a right-hand side of the form (3.36) as where α := 1 1 0 s 2 e −10s ds = 500e 10 e 10 − 61 , while γ is the varying parameter. As shown in [29], 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 was continued up to log γ = 1.75. The integrals representing the distributed delays were again approximated through Clenshaw-Curtis, and the continuation was performed using both Chebyshev and Gauss-Legendre collocation points. Given the absence of an exact expression of the true solution, unlike the case (4.1), the error was computed with respect to a reference solution which was in turn computed using L = 1000 and m = 4. The left plot of

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 [29], a period doubling bifurcation occurs when log γ ≈ 3.8777. In [12] 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 [12,Corollary 15]. Such eigenvalues can, in turn, be computed numerically thanks to the pseudospectral method in [11], which was used in order to obtain Period doubling bifurcations also occur in the case (4.1) with quadratic right-hand side, as shown in [10]. 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.5, 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 methods have 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 [27] general BVPs, as it was the case for the work [6] concerning the corresponding proof for RFDEs. The main result concerns the FEM and states that the error decays as O(L −m ) where L is the increasing number of mesh intervals, while m is the constant degree of the piecewise polynomials used.
The proof we provided considers the periodic BVP formulation (3.1), which features an infinite-dimensional periodicity condition. It is worth remarking that the only reason why we decided to focus on (3.1) instead of the equivalent formulation (2.1) is uniformity with the corresponding proof in [6] for RFDEs, in which case we had encountered some difficulties with the finite-dimensional periodicity formulation. However, we do not expect the same issues to arise in the RE case and, in fact, we plan to work in the near future in order to prove the opposite.  Given both the convergence analysis for RFDEs in [6] and that in the present paper for REs, we expect to be able to rely on the general approach in [27] also for the case of coupled RFDE/RE systems, motivated by their predominance in population dynamics. The proof is currently work in progress and planned to be completed soon.
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 [6] to differential equations with non-constant delays (in particular, state-dependent delays), for which the setting defined in [6] cannot be applied.

A Key points in the proofs of the results in Section 3
In the following we provide some insight into the proofs of Propositions 3.2-3.5 in Subsection 3.1 and 3.7-3.9. in Subsection 3.2, i.e., those needed in order to prove the convergence of the method in [27] for the problem (3.1). Such proofs are entirely based on those for the RFDE case presented in [6], therefore we restrict ourselves to discussing the key differences with respect to the latter, while avoiding the unnecessary repetitions. Such differences are mostly attributable to the lower regularity of the space A with respect to the RFDE case, as well as to the integral form (1.4) of the right-hand side.
The proof of Proposition 3.2 goes as that of [6, Proposition 2.1], meaning that 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. 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 [6, (2.12)], where the derivative of a state of an element of V appears as a factor. Without any assumption whatsoever on F the same would happen in Proposition 3.2, which would be a problem since V is as regular as U in the present case.
The proof of Proposition 3.3 is as immediate as that of [6,Proposition 2.2]. In this case, the norm of the operator turns out to be exactly 1.
Proposition 3.4, just as its RFDE counterpart, can be proved thanks to the fact that u * lies in fact in a more regular subspace of its space U, which is again a consequence of the assumption (1.4).
A little more care is required for Proposition 3.5.
Moving on to the numerical assumptions, Proposition 3.7 can be proved as Proposition 3.2, by replacing F in the first of (3.4) with F M in (3.29). As for Proposition 3.8, the proof is substantially the same as that of [6, 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, an estimate for κ is obtained from the Lebesgue constant of the chosen nodes without recurring to the corresponding constant defined by the derivatives of the Lagrange polynomials.
The proof of Proposition 3.9 is more laborious, just like its counterpart in the RFDE case. The latter has been proved in [6] in several steps, the first of which concerns the invertibility of the operator DΨ L,M (u * , ψ * , ω * ) defined in (3.31) for L, M large enough and can be proved as [6,Proposition 3.11]. The second step concerns the uniform boundedness of DΨ −1 L,M (u * , ψ * , ω * ) and follows the ideas of [6, Lemma 3.12]. The latter is based on [6, 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 [6, 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 [6, Proposition 3.13].

A.1 The nongeneric case
This subsection of the Appendix completes the proof of Propositon 3.5 by showing that k 1 introduced in the proof cannot be 0. Assume for a contradiction that ξ * 1 ∈ R, i.e., that M * (0) 1 = M * 1 ∈ R, and define Y := L ∞ ([0, 1]; R d ) 6 . Consider the standard bilinear form ·, · : Y × X → R defined as In order to obtain a contradiction with (A.5), we resort to adjoint theory for Volterra Integral Equations (VIEs, see [20] as a general reference). Indeed, any RE of the form 6 Here elements in R d are intended as row vectors. with x t 0 = ϕ for some ϕ ∈ X := L 1 ([−1, 0]; R d ) 7 and r ≤ 1 can be written as the VIE Existence and uniqueness [20,Chapter 9] allow to define the forward evolution family 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 [20, Chapter 9] we can express the solution of (A.6) and that of (A.7) respectively as and where R * 0 is the resolvent of (A.6). 7 In our case K * (t, σ − t) := ω * D 2 K(ω * (σ − t), v * (σ)) for t ≥ t 0 and σ ∈ [t − r, t] with r := τ/ω * ≤ 1. 8 It is enough to consider the integrals at the right-hand side of the VIE and of its adjoint over all R, by taking into account the definition of K * 0 .
Given t ∈ R, consider now the pairing [·, ·] t : Y × X → R defined as 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) 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.10), i.e., that Indeed, using (A.9) 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 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.8) we have which proves (A.11).
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 which contradicts (A.5) thanks to the fact that y 0 is the left 1-eigenvector of T * (1, 0).

A.2 Other proofs
This subsection of the Appendix collects the rest of the proofs, which are rather technical, of the propositions stated in Section 3.
Proof of Proposition 3.2. The thesis holds once that (3.8) is proved according to [3, Definition 1.1.1], i.e., u ∈ U, v ∈ V and ω > 0, As for the left-hand side, by using (3.6), the choice of V in (T2) leads to evaluate for t ∈ [0, 1]. Let us define and By (A.13), in order to prove the thesis we need to evaluate A(t) + B(t) + C(t). First of all, observe that With reference to (A.14), this gives where the first equality follows by the fundamental theorem of calculus and the third equality is given by (T3). By (A.18) one obtains and, substituting into (A.14), which in turn gives both and by substituting into (A.15).
The following technical lemma will be frequently used, starting from the proof of Proposition 3.4.
Proof. By (T2), u * and ψ * are bounded. By (1.4) and (T3), F(ψ * ) has a bounded derivative. By (3.5) and the first component of (3.4), u * has a bounded derivative, therefore it is Lipschitz continuous with Lipschitz constant u * ′ ∞ . By the second component of (3.4), ψ * = u * 1 , thus it has a bounded derivative as well and therefore it is Lipschitz continuous with Lipschitz constant ψ * ′ ∞ . Eventually, v * is Lipschitz continuous being it the continuous junction of two Lipschitz continuous functions.
Proof of Proposition 3.4. In this proof, the notations v := G(u, ψ), v * := G(u * , ψ * ), v := G(u, ψ) will be used for brevity. The thesis is equivalent to the existence of r > 0 and κ ≥ 0 such that for all (u, ψ, ω) ∈ B((u * , ψ * , ω * ), r) and all (u, ψ, ω) ∈ U × A × (0, +∞). Since both G and p are linear, it is sufficient to prove the statement with respect to the first component of DΦ, i.e., the one in U. Then, by defining and through (3.9) and (3.10), the goal is to bound for all t ∈ [0, 1] given the choice of U in (T2). As for (A.23), it can be rewritten as for The plan is to bound every single term A i and B i , i = 1, 2, to eventually get the desired bound for P. Then the same will be done for R in (A.25) and S in (A.26), while a bound for Q in (A.24) can be obtained more straightforwardly. Note that all quantities with the superscript * are constant since related to the fixed point. Clearly B 1 can be rewritten as B 1,1 + B 1,2 , where As for B 1,1 , from τ ≤ ω (T4) gives As for the first addend in the right-hand side above, follows directly from (3.5) again. Thus, As for B 1,2 , it is possible to define since the map t → v * t is uniformly continuous (Lemma A.1) and so is DK at the pair (ω * , v * t ) corresponding to the fixed point, under (T4). This leads to the bound which in turn leads to Moreover, B 2 R d ←X×R ≤ τκ 1 ω * (u, ψ, ω) U×A ×R . Eventually, note that every product A i B j for i, j ∈ {1, 2} in the last member of (A.27) contains at least a factor of index 1, which is bounded by some constant times Therefore, for each product there exist κ i,j such that holds for m ∈ {1, 2}. Thus, for r P := min{r/2, 1}, where r is as in (T4), by virtue of (A.28) there exists κ P := max i,j κ i,j ≥ 0 such that for all (u, ψ, ω) ∈ B((u * , ψ * , ω * ), r P ). Since every product contains exactly one A-term and one B-term, the constant κ P can be defined as A bound for the term Q in (A.24) can be retrieved as done for B 1 . Indeed, under (T3) one can define κ D := max since K is uniformly continuous at the pair (v * t , ω * ) corresponding to the fixed point. This leads, as done for B 1 , to the bound Then, by (A.28), for r Q := r/2 and holds for all (u, ψ, ω) ∈ B((u * , ψ * , ω * ), r Q ).
The term R in (A.25), can be written in the form where R(t) is the most subtle term in proving the proposition. The sought bounds can be obtained thanks to the fact that the fixed point lies indeed in a more regular subspace than A × A × [0, +∞) (Lemma A.1, indeed already used for B 1,2 ). In particular, as observed while computing B 1,2 , D 1 can be bounded as while |D 2 | is bounded by a constant, namely |D 2 | ≤ τ ω * . As for E 1 , note that by (T3) the kernel K is Lipschitz-continuous with respect to the second argument with some constant κ N . Moreover, v * is Lipschitz-continuous with some constant κ * by Lemma A.1. Thus Eventually, note that every product D j E k for j, k ∈ {1, 2} in the last member of (A.31) contains a factor of index 1, which are always bounded by some constant times (u, ψ, ω) − (u * , ψ * , ω * ) U×A ×R , while the factors of index 2 are simply bounded by constants. Arguing as done above for P, for r R := min{1, r/2} there exists κ R ≥ 0 such that for all (u, ψ, ω) ∈ B((u * , ψ * , ω * ), r R ). In particular, κ R can be defined as The term S in (A.26), can be written in the form The terms F 1 and F 2 can be bounded as the terms B 1 and B 2 , once observed that the extra factor θ in the integrand function satisfies Eventually, note that every product A j F k for j, k ∈ {1, 2} in the last member of (A.32) contains a factor of index 1, which are always bounded by some constant times (u, ψ, ω) − (u * , ψ * , ω * ) U×A ×R , while the factors of index 2 are simply bounded by constants. Arguing as done above for P and R, for r S := min{1, r/2} there exists κ S ≥ 0 such that for all (u, ψ, ω) ∈ B((u * , ψ * , ω * ), r R ). In particular, κ S can be defined as The thesis eventually follows by choosing r = min{r P , r Q , r R , r S } = min{1, r/2}, κ = κ P + κ Q + κ R + κ S .
Proof of Proposition 3.7. The proof goes as the one of Proposition 3.2, after replacing F with F M .
The proof of Proposition 3.9 is not immediate and needs to be separated into several steps, which require intermediate lemmas. The main step concerns the invertibility of DΨ L,M (u * , ψ * , ω * ). In principle, one could attempt to prove it by resorting to the Banach's perturbation lemma, thus showing first that This in turn would require (3.30) and (3.31). The latter cannot hold for all u ∈ U given the choice of U in (T2). Therefore, in the following we prove the invertibility of DΨ L,M (u * , ψ * , ω * ) directly by following the lines of the proof of Proposition 3.5. To this aim, we first need to show that the initial value problem for is well-posed, and thus defining an associated evolution operator T * L,M (t, s) : X → X is meaningful (for t, s ∈ [0, 1] and t ≥ s). In the sequel it is also convenient to use the abbreviations for ψ ∈ X has a unique solution x L,M .
Proof. Use x = G(u, ψ) according to (3.5). By virtue of (A.34), (A.35) becomes Well-posedness is thus equivalent to the invertibility of I U − π + L ρ + L K * ,+ M : U → U, for which we resort to the Banach's perturbation lemma, since the invertibility of I U − K * ,+ : U → U is guaranteed by the well-posedness of the corresponding initial value problem under (T4). The thesis follows by [7, (A.7) in Lemma A.6]. where ψ L,M is the eigenfunction associated to µ L,M normalized as ψ L,M X = 1.
Proof. We give the proof for s = 0, the extension to s ∈ (0, 1) being straightforward.
Correspondingly, thanks to Lemma A.2, let G(u L,M , ψ) be the solution of (A.33) exiting from the same ψ, where u L,M satisfies The relevant evolution operators are defined, for t ∈ [0, 1], respectively by By recalling that G in (3.5) is linear we get T where r + L,M := (π + L ρ + L K * ,+ M − K * ,+ )u and r − L,M := (π + L ρ + L K * ,− M − K * ,− )ψ. We already showed in the proof of Lemma A.2 that I U − π + L ρ + L K * ,+ M is invertible through the Banach's perturbation lemma. By the latter it is also possible to show that (I U − π + L ρ + L K * ,+ M ) −1 U←U ≤ 2 (I U − K * ,+ ) −1 U←U holds for L and M sufficiently large. Now (A.38) follows since both r + L,M U ≤ π + L ρ + L K * ,+ M − K * ,+ U←U u U and r − L,M U ≤ π + L ρ + L K * ,− M − K * ,− U←A ψ A vanish by the formulation of [7, Lemma A.6] for REs. The second part follows from standard results on spectral approximation of linear operators. In particular, (A.36) implies strongly-stable convergence of µI Y − T * L,M (t, 0) to µI Y − T * (t, 0) for every finite eigenvalue µ of T * (t, 0) [ The first boundary condition in (A.39) gives then The first addend in the right-hand side above vanishes thanks to (N3), [ As for the difference ξ * L,M,2 − ξ * 2 X , it cannot vanish as L and M grow, since u 0 is not necessarily continous. However, it can be observed that  0). Indeed, the right-hand side above vanishes under (N2) thanks to [ The thesis follows from the uniform boundedness of |k 2 − k L,M,2 | with respect to both L and M.
Next, we need to show that the inverse of DΨ L,M (u * , ψ * , ω * ) is bounded uniformly with respect to L and M, which in turn follows from (A.52). We thus need to show that (u L,M , ψ L,M , ω L,M ) U×A×B is bounded uniformly with respect to both L and M. To this aim we prove that (u L,M , ψ L,M , ω L,M ) is related to the solution of the collocation of (the secondary discretization of) an equivalent version of (3.12) according to the primary discretization under (N1) and (N2). Indeed, we first need to rearrange the terms of (3.12) to give a proper sense to the collocation problem since, in general, u is not continuous therein (because of u 0 ), while the range of π + L ρ + L contains only continuous functions. Consider then Note that existence and uniqueness of (e + L,M , e − L,M ) follows already from Propositions A.4 and 3.5, so that the invertibility of the operator is already proved. Anyway, if we manage to prove that lim L,M→∞ then we can apply the Banach's perturbation lemma to recover the bound On the other hand, ε − L,M does not necessarily vanish but is anyway bounded uniformly with respect to L and M since u 0 is not necessarily continuous, even though it is bounded. Consequently, it is not difficult to argue that ε − L,M A ≤ ψ 0 A + 2 u 0 U by taking into account for possible jumps in u 0 . It is left to prove that ε ω,L,M is uniformly bounded. From (A.60) we have in which the third addend at the right-hand side is uniformly bounded thanks to Proposition A.5, which also proves that ω L,M is uniformly bounded. Therefore, the first and the second addends vanish, thanks to the same arguments adopted for (A.65) under (N3) and (N4) and also thanks to [7, (A.2) in Lemma A.1].
In the proof of Proposition A.5, it is also shown that ψ L,M is uniformly bounded. Eventually, we obtain that (u L,M , ψ L,M , ω L,M ) U×A×B is bounded uniformly with respect to both L and M thanks to (A.58) and Proposition 3.5.
Proof of Proposition 3.9. Thanks to Lemma A.6 and to the fact that r 1 and κ in Proposition 3.8 are independent of L and M, it remains to prove that Ψ L,M (u * , ψ * , ω * ) U×A×B vanishes. We have Ψ L,M (u * , ψ * , ω * ) U×A×B ≤ (I U×A×B − P L R L )(u * , ψ * , ω * ) U×A×B since Φ(u * , ψ * , ω * ) = (u * , ψ * , ω * ). The second addend in the right-hand side above vanishes under (N3) and thanks to the result corresponding to [7,Corollary A.3] for REs. The first addend vanishes as well by Lemma A.1, which shows in particular that u * and ψ * are continuous.