A practical approach to computing Lyapunov exponents of renewal and delay equations

We propose a method for computing the Lyapunov exponents of renewal equations (delay equations of Volterra type) and of coupled systems of renewal and delay differential equations. The method consists in the reformulation of the delay equation as an abstract differential equation, the reduction of the latter to a system of ordinary differential equations via pseudospectral collocation, and the application of the standard discrete QR method. The effectiveness of the method is shown experimentally and a MATLAB implementation is provided.


Introduction
A delay equation is a functional equation consisting of "a rule for extending a function of time towards the future on the basis of the (assumed to be) known past" [29].A renewal equation (RE) is a delay equation of Volterra type, i.e., the rule for extension prescribes the value of the unknown function itself, instead of the value of its derivative, as in the case of delay differential equations (DDEs).
The goal of this work is to compute the (dominant) Lyapunov exponents (LEs) of REs and of coupled systems of REs and DDEs (henceforth coupled equations).The usefulness of LEs for measuring the asymptotic exponential behavior of solutions is well known; for example, they can be used to study the average asymptotic stability of solutions, the insurgence of chaotic dynamics and the effects of perturbations on the system, as well as to estimate the entropy or the dimension of attractors.
As for DDEs, recent methods for computing the LEs have been proposed, particularly [20] and [12], which use two different approaches (for other methods, see the references in the cited works).
In [12] the DDE is reformulated as an abstract differential equation and a pseudospectral discretization is applied [13], yielding a system of ordinary differential equations (ODEs); LEs are then computed by using the standard discrete QR method (henceforth DQR) for ODEs proposed in [24,26].In [20], instead, the problem is tackled directly: the DDE is posed in an infinite-dimensional Hilbert space as the state space, the associated family of evolution operators is discretized and the DQR is adapted and applied to the finite-dimensional approximation; for the error analysis, the DQR is raised to infinite dimension and compared to the approximated DQR used for the computations.
As for REs, as far as we know, there are no methods available in the literature for computing LEs.Only a first example of naive computation can be found in [14], where it is done simply to exemplify the versatility of the collocation techniques used therein, without attempting a rigorous formulation and error analysis.
In the present work of computational nature, we develop a practical method following the approach of [12] described above, and based on [36] for the reformulation of REs into abstract differential equations.As in [12], we use the DQR to compute the LEs of the approximating ODE, but, in principle, any method can be used; our choice was motivated by our goal of providing a practical way of computing LEs by using ready-to-use code for a well-known method.
In Section 2 we recall the DQR for linear ODEs.Then, in Section 3 we define the reformulation of REs, DDEs and coupled equations into abstract differential equations and their pseudospectral collocation into ODEs.After describing the implementation choices in Section 4, we present in Section 5 some numerical experiments concerning the convergence of the method for an example RE with many known properties, as well as some examples of computation of LEs of REs and coupled equations.Finally, we present some concluding remarks in Section 6.
The MATLAB codes implementing the method and the scripts to reproduce the experiments of Section 5 are available at http://cdlab.uniud.it/software.

DQR for ODEs
In this section, we first illustrate the DQR to compute the LEs of linear nonautonomous ODEs; in the nonlinear case, one previously linearizes around a reference trajectory in the attractor.Then, we comment on the relevant literature.
Let n be a positive integer and consider the ODE for A : [0, +∞) → R n×n continuous and bounded; also, let Z(t) be the fundamental matrix solution exiting from a given nonsingular matrix Z 0 ∈ R n×n prescribed at time 0 without loss of generality.For any sequence {t k } k∈N of time instants strictly increasing A practical approach to computing Lyapunov exponents of renewal and delay equations Breda, Liessi from t 0 = 0, construct the iterative QR factorization3 starting from Z 0 = Q 0 R 0 and, at each step j = 1, . . ., k, solving the n initial value problems (IVPs) and factorizing the solution at t j as If S(t, s) := Z(t)Z(s) −1 is the state transition matrix associated with (1), then The uniqueness of the QR factorization and (2) give so that, eventually, (upper4 ) LEs are recovered as Above, [R j,j−1 ] i,i denotes the i-th diagonal entry of the j-th triangular factor R j,j−1 .
Obviously, in the implementation ( 5) is truncated to some large T > 0. In the end, each step of the DQR requires the solution of the IVPs (3) and the QR factorization (4).
The above summary was taken mainly from [12], where the DQR is applied to the ODE obtained from the pseudospectral collocation of a given DDE (see Section 3), thus following the original approach of [13] to also address the study of chaotic dynamics.As anticipated in Section 1, the aim of the present work is to extend this procedure to more general classes of delay equations, such as REs and coupled equations.Once the pseudospectral collocation is performed (possibly after linearization, see Remark 4.1 later on), the outcome is an ODE like (1); thus, the DQR applies unchanged, independent of the original delay equation.
The literature on the theory and computation of LEs of ODEs is ample; for a starting reference, see [25], but see also [2] as a reference monograph.QR methods were first proposed in the pioneering works [8,9]; for a complete discussion of the discrete version, see [24].The literature on the computation of LEs of delay equations is mostly of an experimental flavor and, to the best of the authors' knowledge, restricted only to DDEs.As initial references, we can suggest [11,20,31], but see also [22] for a more recent method to reduce DDEs to ODEs by using Galerkin-type projections.Note also that all of these works rely on a Hilbert state space setting to legitimize orthogonal projections, while the technique in [12] is free from this constraint and thus maintains the classical state spaces (typically continuous functions for DDEs and absolutely integrable ones for REs).Beyond the lack of relevant methods, this is part of the motivation for the extension of the approach proposed in [12] to more general delay equations.

Pseudospectral collocation
In this section, we illustrate the use of pseudospectral collocation to reduce delay equations to ODEs, in view of the application of the DQR described in Section 2. For the reader's convenience, we first present, separately, the discretization of an RE in Section 3.1 and that of a DDE in Section 3.2, summarizing from, respectively, [36] and [12] the main aspects for the present objective (for a full treatment, see again [12,36] and also [13]).Eventually, we combine the two approaches in Section 3.3 for a coupled equation.
In what follows, we use the subscripts X and Y to refer, respectively, to REs and DDEs.

Pseudospectral collocation of REs
Let τ > 0 be real and d X > 0 be an integer.Consider the IVP for an RE given by where ϕ ∈ L In [36] an efficient application of pseudospectral collocation to reduce (6) to an IVP for an ODE is proposed based on an equivalent formulation of (6) as an abstract Cauchy problem (ACP) describing the evolution of an integral of the original state x t .In particular, by defining the Volterra integral operator V : L 1 → AC 0 as where AC 0 := AC 0 ([−τ, 0]; R d X ) is the space5 of absolutely continuous functions vanishing at 0, it turns out that ( 6) is equivalent to the ACP A 0,X µ := µ ′ , D(A 0,X ) := {µ ∈ AC 0 : µ = V η for some η ∈ NBV 0 }, where NBV 0 is the space of functions of bounded variation vanishing at 0 and continuous from the right, and q X ∈ NBV 0 is defined as In order to discretize (7), consider a mesh a positive integer.Correspondingly, let P M X ,X : R M X d X → NBV 0 be the interpolation operator on {0} ∪ Ω M X ,X with value 0 at θ 0,X := 0, i.e., where {ℓ 0,X , ℓ 1,X , . . ., ℓ M X ,X } is the Lagrange basis on {0} ∪ Ω M X ,X , and let R M X ,X : NBV 0 → R M X d X be the restriction operator (R M X ,X µ) j := µ(θ j,X ), j = 1, . . ., M X .
Then, the discrete version of ( 7) is given by with 6 and where where Remark 3.1.Instead of NBV 0 , [36] uses the space NBV of functions of bounded variation that vanish at 0 and are continuous from the right on (−τ, 0), but not necessarily at −τ.In that setting, C 0,X := (V ↾ NBV ) −1 is a multi-valued operator, defined as8 C 0,X µ := {η : µ = Vη}, since functions differing only by the jump at −τ are mapped by V to the same element of NBV.
The trivial semigroup {S 0,X (t)} t≥0 on NBV defined as is not strongly continuous.However, its restriction {T 0,X (t)} t≥0 to AC 0 is strongly continuous and A 0,X is its infinitesimal generator.
From this point of view, the semilinear ACP (7) renders a clear separation between the translation along the solutions (through the linear semigroup {T 0,X (t)} t≥0 and its infinitesimal generator A 0,X ) and the rule for extension (basically through the nonlinear right-hand side F of the specific RE), which are the two ingredients of a delay equation.For these and related aspects of the theory of delay equations, see [27,28] for the sun-star (⊙ * ) theory and [29,36] for the more recent twin semigroup theory.

Pseudospectral collocation of DDEs
Let τ > 0 be real and d Y > 0 be an integer.Consider the IVP for a DDE given by where In [12] (but see also [13]) pseudospectral collocation is used to reduce (10) to an IVP for an ODE based on an equivalent formulation of (10) as an ACP describing the evolution of the original state y t , viz.
In order to discretize (11), consider a mesh Then, the discrete version of ( 11) is given by where Remark 3.2.Note that the ACP (11) can be alternatively described as the equivalent semilinear where A 0,Y is the infinitesimal generator of the shift semigroup {T 0,Y (t)} t≥0 defined as and q Y ∈ L ∞ is defined as Equation (13) renders for DDEs the same separation between translation along the solutions and rule for extension as illustrated in Remark 3.1 for REs (see again [28]).The pseudospectral collocation of (13) leads, again, to (12), which can be rewritten equivalently as and all of the others equal to zero.Now, (14) resembles (13).

Pseudospectral collocation of coupled equations
Let τ, d X and d Y be as above.Consider the IVP for a coupled equation given by where For well-posedness, see [27].By combining the approaches of the previous sections, it follows that ( 15) is equivalent to the ACP (u The discrete version of ( 16) is as follows: Note that, now, The total number of approximating ODEs is Remark 3.3.Following Remarks 3.1 and 3.2, it is not difficult to see that ( 16) is equivalent to the semilinear ACP where A 0,X,Y := diag(A 0,X , A 0,Y ) is linear and is nonlinear.The pseudospectral collocation of (18) leads, again, to (17), where, correspondingly, the ODEs for V can be compacted as done in Remark 3.2 for DDEs; see (14).The (numerical) analysis of (18) is current work in progress at CDLab,9 also in view of the corresponding sun-star theory of coupled equations developed in [27] and of the more recent twin semigroup theory of [29].

Implementation
Due to our choice of example equations (see Section 5), in our implementation we considered scalar equations (d X = d Y = 1) of the following kinds: for REs, and for coupled equations, 10 where τ 2 > τ 1 > 0, and f , f 1 , f 2 and g are (possibly) nonlinear functions R → R. Nevertheless, generalizing to other forms of equations is usually fairly straightforward. 11 We implemented the pseudospectral discretization 12 using Chebyshev nodes of type II (extrema) as the meshes {0} ∪ Ω M X ,X and Ω M Y ,Y of points in [−τ, 0] with M X = M Y .To compute the nodes and the corresponding differentiation matrix, we used the cheb routine of [38,Chapter 6].For the interpolation, we used the barycentric Lagrange interpolation formula [10]; the barycentric weights corresponding to Chebyshev extrema are explicitly known and are given therein.For the quadrature of the integrals, we used the Clenshaw-Curtis formula [23,39].We implemented (A 0,X P M X ,X Φ)(θ) as θ)Φ j , computing ℓ ′ j,X as the polynomial interpolating the j-th column of the differentiation matrix, again with the barycentric formula.
In order to apply the DQR described in Section 2 to the approximating ODE, the latter needs to be linearized around a reference solution.The linearization is done explicitly.The solutions are computed by using MATLAB's ode45, which implements the embedded Dormand-Prince (5, 4) method [30,37].For the differential part of (20), the initial value consists of the vector of values of the chosen initial function at Ω M Y ,Y , while, for (19) and the renewal part of (20), the vector D −1 M X ,X u is used, 13 where u is the vector of values of the initial function at Ω M X ,X .
Finally, the DQR for a linear ODE is implemented in the dqr routine of [12], which follows [26].Therein, the IVPs (3) are again solved with the Dormand-Prince (5, 4) pair; however, instead of adapting the step size (initially 0.01) based on the error between the two solutions, the automatic adaptation controls the error between the corresponding LEs.As an initial guess for the fundamental matrix solution, a random matrix is used.The computation is stopped when the specified truncation time T is reached. 10Observe that the unknown of the differential equation is not delayed. 11For instance, one can consider models such as (3.3) in [35], yet with finite age-span, which do not enter class (20).In such cases, the method is implemented following the same strategy described here (discretization of the nonlinear system for computing solutions, linearization and discretization of the linearized system for computing the LEs).However, the authors are currently implementing a general code for a larger class of equations (ideally the most general (15)), which was beyond the scope of this work. 12For more details on pseudospectral methods, see also [38]. 13In (9) the initial value is specified as R M X ,X V ϕ, i.e., the vector representing the polynomial interpolating the exact integral of the initial value ψ.Another approach is to use R M X ,X V P M X ,X R M X ,X ϕ, in which the integral of the polynomial interpolating ψ is used.In our implementation we use neither; our choice is computationally easier and is motivated as follows.
As already noted, in order to represent the integrated state, only the vector U of values at Ω M X ,X is needed, as the value at θ 0,X = 0 is always 0. Computing the derivative of the interpolating polynomial by applying the differentiation matrix to (0, U) T (where the 0 stands for a column vector of d X zeros), we obtain (d M X ,X U, D M X ,X U) T , where d M X ,X ∈ R d X ×M X d X is a row of d X × d X -block entries ℓ ′ j,X (θ 0,X )I d X for j = 1, . . ., M X .Since deriving a polynomial lowers its degree by one, D M X ,X U uniquely determines the derivative of the polynomial represented by U, which motivates our use of D −1 M X ,X u.
Remark 4.1.For REs only, and in particular for the example described in Section 5.1, we experimented also with a different method, based on computing a solution of (19), linearizing the latter around the former and applying the pseudospectral collocation to the resulting linear RE.We computed the solution of the RE with the method described in [34], which is based on the trapezoidal quadrature formula on a uniform grid in [−τ 2 , 0] with the constraint that −τ 1 must be a grid point.Corresponding to a solution x, we considered the linear RE14 See Section 5.1 for a comparison of the approaches.
Remark 4.2.In Remark 4.1 (more precisely in Footnote 14), we observed that, in most cases, REs cannot be linearized.However, in many of those cases, the ODE resulting from the pseudospectral discretization can, in fact, be linearized; for example, the ODEs resulting from ( 19) and (20) can be linearized if f , f 1 , f 2 and g are differentiable.As an example, the linearization of (17) around a solution (U, V) is as follows: where J indicates the Jacobian matrix, JF M X (U(t), In Section 5.1 below, we explicitly show the linearized ODE for an example RE.
We recall that the MATLAB codes implementing the method and the scripts to reproduce the experiments of Section 5 are available at http://cdlab.uniud.it/software.

Results
We present here three example equations: an RE with a quadratic (logistic-like) nonlinearity in Section 5.1, an RE modeling egg cannibalism in Section 5.2 and a simplified version of the Daphnia model with a logistic term for the growth of the resource in Section 5.3.In particular, we use the first example to test the proposed method also from the numerical point of view; we then apply it to the second and third example to compute the exponents.

RE with quadratic nonlinearity
The first equation we study is the RE with quadratic nonlinearity from [14]: i.e., (19) with τ 1 := 1, τ 2 := 3 and f (x) := γ 2 x(1 − x).Its equilibria and their stability properties are known; in particular, its nontrivial equilibrium undergoes a Hopf bifurcation for γ = 2 + π 2 and the branch of periodic solutions arising from there has the analytic expression Observe that the period is 4, independent of γ.Moreover, it is experimentally known that it presents several period-doubling bifurcations, possibly leading to a cascade and, eventually, to chaos [14].
Since several properties of ( 22) are analytically known, we use it to test the effectiveness and efficiency of the method proposed for LE computation, and to compare it to the alternative approach described in Remark 4.1.
For equilibria, the LEs are the real parts of the eigenvalues λ of the infinitesimal generator of the semigroup of solution operators, which are related to the eigenvalues µ of the solution operator that advances the solution by h via µ = e λh . ( For periodic solutions, the LEs are the real parts of the Floquet exponents, which are related to the Floquet multipliers (i.e., the eigenvalues of the monodromy operator) via (24), where µ, λ and h are, respectively, a Floquet multiplier, a Floquet exponent and the period.In both cases, we can thus obtain the LEs by computing the eigenvalues µ of an evolution operator with any time step h for the equilibria (we choose h = τ 2 = 3 for ( 22)) and a time step h equal to the period for periodic solutions (h = 4 for ( 23)), and then computing the real part of log(µ)/h.In order to obtain reference values for our experiments, we compute the spectra of evolution operators with the method of [16], which is based on the pseudospectral collocation of the operator; we use the implementation eigTMNpw of [18,19].
Although computing the solutions of delay equations is not the focus of this work, given that both the main approach and the alternative one of Remark 4.1 involve computing solutions, our first experiment compares the error of the computed solutions with respect to the known periodic solutions of (22).We choose γ = 4 > 2 + π 2 , which corresponds to a stable periodic solution, since the first period-doubling bifurcation is experimentally known to happen at γ ≈ 4.32 [14,16].
Figure 1 shows the errors on the solution of the approximating ODE and on the solution of the original RE (22) with respect to the number of nodes (minus 1) in the grid in [−3, 0], i.e., M X for the pseudospectral discretization and 3r for the trapezoidal method 15 [34]; in both cases, two errors are measured, namely the absolute error at t = 500 and the maximum absolute error on a grid of points in [0, 500] (a uniform grid with step 0.05 for the pseudospectral approach, the time points given by the trapezoidal method for the alternative approach).To solve the approximating ODE given by the pseudospectral discretization, we used ode45 with RelTol = 10 −6 and AbsTol = 10 −7 , which justifies the barrier on the error in Figure 1.
The experiment confirms that the trapezoidal method has order 2, as proved in [34], and that the pseudospectral discretization has infinite order, which is often the case for pseudospectral methods applied to smooth problems [38].Even for rather small values of M X = 3r, the error for the pseudospectral method is several orders of magnitudes smaller than the one for the trapezoidal method. 16n the next experiment, we investigate how the errors on the LEs depend on the choice of M X and of the final time17 T. We choose values of γ corresponding to the stable trivial equilibrium (γ = 0.5), the stable nontrivial equilibrium (γ = 3) and the stable periodic orbit (γ = 4).
Since we are going to use the linearization of the ODE (9) coming from the RE (22), as an example, we show it explicitly here.With reference to Remark 4.2, observe that the right-hand side of ( 22) is not Fréchet-differentiable as a map from L 1 to R, while the right-hand side of the discretized equation is differentiable.The linearization of the approximating ODE around the solution U is given by where JF M X (U(t)) is a row vector with components  22) with γ = 4 with respect to the known periodic solution (23), computed via pseudospectral discretization (solid lines) and directly with the trapezoidal method (dashed lines), measured as the absolute error at t = 500 (•) and as the maximum absolute error on a grid of points in [0, 500] (•), when varying the number of nodes (minus 1) in the grid in [−3, 0], i.e., M X for the pseudospectral discretization and 3r for the trapezoidal method.The exact periodic solution ( 23) is used as the initial value.
In Figures 2 and 3, we can see the absolute errors on the dominant LE increasing either M X or T. The tolerance for dqr is 10 −6 , while those for ode45 are RelTol = 10 −6 and AbsTol = 10 −7 .The reference values are obtained by using eigTMNpw with the default options and 120 as the degree of the collocation polynomials (fixed independently of M X ).
In Figure 2 the final time T = 1000 is fixed and M X increases.We can observe that, apart from very low values of M X , the error reaches a barrier. 18We performed the same experiment with T = 10000 and could make the same observation, although the barrier was smaller by about one order of magnitude.The barrier depends on the error due to the time truncation in (5).Indeed, Figure 3, where M X is constant and T varies, shows that the LEs converge linearly (confirming what is explained in [20]).In Figure 2 the truncation error appears to dominate on the error due to the collocation.
For the dominant nontrivial exponent 19 of the periodic solution, we observe in Figure 3 that, for M X = 8, the error seems to reach a barrier, indicating that more ODEs are necessary to reproduce the properties of the original RE more accurately, as it is reasonable to expect.In other experiments, we observed that, as T increases, error Figure 3. Absolute errors on the dominant LEs of the RE with quadratic nonlinearity (22) for values of γ corresponding to the stable trivial equilibrium (γ = 0.5), the stable nontrivial equilibrium (γ = and the stable periodic orbit (γ = 4).For the last one, both the trivial and the dominant nontrivial exponents are shown.The errors are measured with respect to the exponents computed via eigTMNpw.The exponents are computed for M X ∈ {8, 12, 16, 20} as shown in the legend.
back to chaos (starting at γ ≈ 4.8885).As an example, Figure 4 (second and third row) shows some stable periodic solutions in the branches arising from the first and the second set of bifurcations.Observe that indeed the period approximately doubles at each bifurcation.

Egg cannibalism model
The second equation we consider is the egg cannibalism (toy) model from [15]: with a mat and a max being, respectively, the constant maturation and maximum ages.
Observe that (25) corresponds to (19) with τ 1 := a mat , τ 2 := a max and f (x) := γ 2 xe −x .Also in this case, the equilibria and their stability properties are known, including the occurrence of a Hopf bifurcation for the nontrivial equilibrium at log(γ) = 1 + π 2 , although here the periodic solutions are not explicitly known; again, the presence of period-doubling bifurcations is experimentally known [13][14][15]36].
Figure 5 (top row) presents the diagram of the first two dominant LEs of ( 25) when varying γ, with a mat = 1 and a max = 3.The numerical parameters are M X = 15, T = 1000, a tolerance of 10 −6 for dqr, and RelTol = 10 −6 and AbsTol = 10 −7 for ode45.Similar to the previous example, the dominant exponent is 0 at the Hopf bifurcation, and one exponent remains 0 thereafter.The dominant nontrivial exponent reaches 0 again for the expected period-doubling bifurcations at log(γ) ≈ 3.855 ( [36] finds log(γ) ≈ 3.8777 with M X = 20 and log(γ) ≈ 3.8763 with M X = 40, setting MatCont's tolerances to 10 −10 for Newton's method and 10 −6 for the calculation of the test functions for bifurcation points) and log(γ) ≈ 4.54, and it becomes positive for log(γ) ≥ 4.66, indicating chaos.Figure 5 (bottom row) shows some stable periodic solutions, confirming the (approximate) doubling of the period.

Simplified logistic Daphnia
The third and final equation is a simplified version of the Daphnia model with a logistic term as the growth of the resource, taken from [15]: where b is the birth rate of the consumer population, S is the density of the resource, r and K are, respectively, the growth rate and the carrying capacity of the resource, and   a mat and a max have the same meaning as in Section 5.2.20 Equation ( 26) corresponds to (20) with τ 1 := a mat , τ 2 := a max , f 1 (x) := βx, f 2 (x) := −γx and g(y) := ry(1 − y K ).The equilibria are known, along with the stability properties of the trivial (zero and consumer-free) ones; in particular, the consumer-free equlibrium exchanges stability with the nontrivial equilibrium in a transcritical bifurcation at β = (K(a max − a mat )) −1 ; moreover, when varying β, the nontrivial equilibrium is experimentally known to undergo a Hopf bifurcation [13,15,17].
The diagram of the first two dominant LEs of (26) when varying β is shown in Figure 6.The values of the other model parameters are a mat = 3, a max = 4 and r = K = γ = 1.The numerical parameters are M X = M Y = 15, T = 1000, a tolerance of 10 −6 for dqr, and RelTol = 10 −4 and AbsTol = 10 −5 for ode45.We can observe a spike at β = 1 (albeit not touching the value 0), corresponding to the known transcritical bifurcation, while the Hopf bifurcation seems to happen for γ ≈ 3 ([13] finds γ ≈ 3.0161 with M X = 10 and MatCont's tolerances set to 10 −10 ).We continued the experiment for values of β higher than those shown in Figure 6, but we did not find other bifurcations.
Compared with the previous examples, the diagram seems less accurate (observe the spike corresponding to the transcritical bifurcation and the trivial LE after the Hopf bifurcation).The explanation for this phenomenon is still unknown.In this example, we have increased the tolerances for ode45 in order to improve the computation times; however, for β = 1 the dominant LE differs absolutely from the one computed with RelTol = 10 −6 and AbsTol = 10 −7 (as in the previous examples) by less than 10 −7 and is thus still substantially far from 0 (as a further example, for β = 5 the absolute difference is larger but still less than 10 −5 ).Moreover, as β increases, the periodic solution presents flat regions followed by spikes, which may suggest that the equation is becoming stiff; however, we tried replacing ode45 with MATLAB's ode23s, implementing a modified Rosenbrock formula of order 2 for stiff ODEs [37], with no substantial changes.

Concluding remarks
In this work, we have provided the first method, to our knowledge, for computing LEs of REs and coupled equations.The proposed method appears to be effective when applied to examples with known properties; however, since the nature of our study has been purely computational, further investigation into the method's convergence properties is required.As far as efficiency is concerned, LEs are notoriously expensive to compute [20], and that is true also in this case; the computational cost depends linearly on T, while its relation with M X and M Y is currently unclear from the experiments, even though it is expected to be of order 4, according to (3).
The next step in our research is to tackle the problem following the approach of [20]; as recalled in Section 1, this involves the reformulation of the equation (RE or coupled equation, in our case) in a Hilbert space and the rigorous definition of LEs and the DQR in the new setting.The technique of [12] and of the present work, however, is a rather general approach which can also be applied to other kinds of equations (e.g., certain classes of partial differential equations of interest for mathematical biology [1]), as anticipated in [12] and according to the pragmatic point of view discussed in [14].

Figure 1 .
Figure1.Errors on the solution of the RE with quadratic nonlinearity(22) with γ = 4 with respect to the known periodic solution(23), computed via pseudospectral discretization (solid lines) and directly with the trapezoidal method (dashed lines), measured as the absolute error at t = 500 (•) and as the maximum absolute error on a grid of points in [0, 500] (•), when varying the number of nodes (minus 1) in the grid in [−3, 0], i.e., M X for the pseudospectral discretization and 3r for the trapezoidal method.The exact periodic solution (23) is used as the initial value.

Figure 4 .
Figure 4. Diagram of the first two dominant (in descending order) LEs (top row) and solutions (other rows) of the RE with quadratic nonlinearity (22) when varying γ, computed with M X = 15 and T = 1000.The solutions are computed via pseudospectral discretization with M X = 15, starting from a constant initial value of 0.2.The final time of T = 1000 ensures sufficiently good convergence to the stable periodic solution.For each solution, the last two periods are shown, separated by a vertical dashed line.The values of γ and of the period are given above the diagrams; the values of γ are also marked by vertical dashed lines in the diagram of the LEs.

30 Figure 5 .
Figure 5. Diagram of the first two dominant (in descending order) LEs (top row) and solutions (bottom row) of the egg cannibalism RE (25) with a mat = 1 and a max = 3, when varying γ, computed with M X = 15 and T = 1000.The solutions are computed via pseudospectral discretization with M X = 15, starting from a constant initial value of 0.2.The final time of T = 1000 ensures sufficiently good convergence to the stable periodic solution.For each solution, the last two periods are shown, separated by a vertical dashed line.The values of γ and of the period are given above the diagrams; the values of γ are also marked by vertical dashed lines in the diagram of LEs.

Figure 6 .
Figure 6.Diagram of the first two dominant (in descending order) LEs of the simplified logistic Daphnia equation (26) with a mat = 3, a max = 4 and r = K = γ = 1, when varying β, computed with M X = 15 and T = 1000.