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

1.   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" ^[1]. 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 ^[2] and ^[3], which use two different approaches (for other methods, see the references in the cited works).

In ^[3] the DDE is reformulated as an abstract differential equation and a pseudospectral discretization is applied ^[4], 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 ^[5,6]. In ^[2], 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 ^[7], 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 ^[3] described above, and based on ^[8] for the reformulation of REs into abstract differential equations. As in ^[3], 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.

2.   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\colon[0, +\infty)\to\mathbb{R}^{n\times n}$ continuous and bounded; also, let $Z(t)$ be the fundamental matrix solution exiting from a given nonsingular matrix $Z_{0}\in\mathbb{R}^{n\times n}$ prescribed at time $0$ without loss of generality. For any sequence $\{t_{k}\}_{k\in\mathbb{N}}$ of time instants strictly increasing from $t_{0} = 0$, construct the iterative QR factorization¹

¹In what follows, a QR factorization of a nonsingular matrix is intended as the unique one with positive diagonal elements.

starting from $Z_{0} = Q_{0}R_{0}$ and, at each step $j = 1, \dots, 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 (2.1), then

The uniqueness of the QR factorization and (2.2) give

so that, eventually, (upper²) LEs are recovered as

² Lower exponents come either as $\liminf$ or as upper exponents of the adjoint system. Note, however, that for regular ODEs in the sense of Lyapunov (see, e.g., [9, Definition 3.5.1]) the LEs exist as exact limits, and such quantities are meaningful for stability statements in the original nonlinear system, thus avoiding the so-called Perron effect (see, e.g., ^[10]). Nevertheless, an in-depth study of the theoretical requirements of the original delay equation guaranteeing the regularity of the ODE obtained by pseudospectral collocation is beyond the scope of the present work, as well as the extension of the regularity concept itself to the infinite-dimensional case of delay equations (extension that the authors, to the best of their knowledge, are not aware of; see, anyway, ^[11] and the references therein).

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 (2.5) is truncated to some large $T > 0$. In the end, each step of the DQR requires the solution of the IVPs (2.3) and the QR factorization (2.4).

The above summary was taken mainly from ^[3], 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 ^[4] 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 (2.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 ^[12], but see also ^[9] as a reference monograph. QR methods were first proposed in the pioneering works ^[13,14]; for a complete discussion of the discrete version, see ^[6]. 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 ^[2,15,16], but see also ^[17] 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 ^[3] 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 ^[3] to more general delay equations.

3.   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, ^[8] and ^[3] the main aspects for the present objective (for a full treatment, see again ^[3,8] and also ^[4]). 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.

3.1.   Pseudospectral collocation of REs

Let $\tau > 0$ be real and $d_{X} > 0$ be an integer. Consider the IVP for an RE given by

where $\phi\in L^{1}:= L^{1}([-\tau, 0];\mathbb{R}^{d_{X}})$, $F\colon L^{1}\to\mathbb{R}^{d_{X}}$ and $x_{t}$, defined as $x_{t}(\theta):= x(t+\theta)$ for $\theta\in[-\tau, 0]$, denotes the history or state function (so that $x_{0} = \phi$ represents the initial state). If $F$ is globally Lipschitz, the IVP (3.1) has a unique solution on $[-\tau, +\infty)$ [18, Theorem 3.8].

In ^[8] an efficient application of pseudospectral collocation to reduce (3.1) to an IVP for an ODE is proposed based on an equivalent formulation of (3.1) 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 $\mathcal{V}\colon L^{1}\to AC_{0}$ as

where $AC_{0}:= AC_{0}([-\tau, 0];\mathbb{R}^{d_{X}})$ is the space³ of absolutely continuous functions vanishing at $0$, it turns out that (3.1) is equivalent to the ACP

³Hereinafter, we do not indicate the domain and codomain of a function space when clear from the context.

through $u(t) = \mathcal{V}x_{t}$. Above, $\mathcal{A}_{0, X}\colon\mathcal{D}(\mathcal{A}_{0, X})\subset NBV_{0}\to NBV_{0}$ is defined as $(\mathcal{V} _{{\upharpoonright}{NBV_{0}}})^{-1}$, i.e.,

where $NBV_{0}$ is the space of functions of bounded variation vanishing at $0$ and continuous from the right, and $q_{X}\in NBV_{0}$ is defined as

In order to discretize (3.2), consider a mesh $\Omega_{M_{X}, X}$ of $M_{X}$ points $-\tau\leq\theta_{M_{X}, X} < \dots < \theta_{1, X} < 0$ with $M_{X}$ a positive integer. Correspondingly, let $P_{M_{X}, X}\colon\mathbb{R}^{M_{X}d_{X}}\to NBV_{0}$ be the interpolation operator on $\{0\}\cup\Omega_{M_{X}, X}$ with value $0$ at $\theta_{0, X}:= 0$, i.e.,

where $\{\ell_{0, X}, \ell_{1, X}, \dots, \ell_{M_{X}, X}\}$ is the Lagrange basis on $\{0\}\cup\Omega_{M_{X}, X}$, and let $R_{M_{X}, X}\colon NBV_{0}\to\mathbb{R}^{M_{X}d_{X}}$ be the restriction operator

Then, the discrete version of (3.2) is given by

with $U(t)\in\mathbb{R}^{M_{X}d_{X}}$ that approximates the integrated state $\mathcal{V}x_{t}$ according to $U_{j}(t)\approx(\mathcal{V}x_{t})(\theta_{j, X})$, $j = 1, \dots, M_{X}$⁴, and where $D_{M_{X}, X}:= R_{M_{X}, X}\mathcal{A}_{0, X}P_{M_{X}, X}\in\mathbb{R}^{M_{X}d_{X}\times M_{X}d_{X}}$ has $d_{X}\times d_{X}$-block entries

⁴Note then that $\mathcal{A}_{0, X}P_{M_{X}, X}U(t)\approx x_{t}$.

where $I_{d_{X}}$ is the identity on $\mathbb{R}^{d_{X}}$, $F_{M_{X}}:= F\circ\mathcal{A}_{0, X}P_{M_{X}, X}$ and $\mathbf{1}_{M_{X}, X}\in\mathbb{R}^{M_{X}d_{X} \times d_{X}}$ has all $d_{X}\times d_{X}$-block entries $I_{d_{X}}$ ⁵.

⁵Note that $\mathbf{1}_{M_{X}, X}$ discretizes $-q_{X}$.

Remark 3.1. Instead of $NBV_{0}$, ^[8] uses the space $NBV$ of functions of bounded variation that vanish at $0$ and are continuous from the right on $(-\tau, 0)$, but not necessarily at $-\tau$. In that setting, $\mathcal{C}_{0, X}:= (\mathcal{V} _{{\upharpoonright}{NBV}})^{-1}$ is a multi-valued operator, defined as⁶ $\mathcal{C}_{0, X}\mu := \{\eta:\mu = V\eta\}$, since functions differing only by the jump at $-\tau$ are mapped by $\mathcal{V}$ to the same element of $NBV$. The trivial semigroup $\{S_{0, X}(t)\}_{t\geq0}$ on $NBV$ defined as

⁶It turns out that $D(\mathcal{C}_{0, X}) = D(\mathcal{A}_{0, X})$ (see (3.3)).

is not strongly continuous. However, its restriction $\{T_{0, X}(t)\}_{t\geq0}$ to $AC_{0}$ is strongly continuous and $\mathcal{A}_{0, X}$ is its infinitesimal generator.

From this point of view, the semilinear ACP (3.2) renders a clear separation between the translation along the solutions (through the linear semigroup $\{T_{0, X}(t)\}_{t\geq0}$ and its infinitesimal generator $\mathcal{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 ^[18,19] for the sun–star ($\odot\ast$) theory and ^[1,8] for the more recent twin semigroup theory.

3.2.   Pseudospectral collocation of DDEs

Let $\tau > 0$ be real and $d_{Y} > 0$ be an integer. Consider the IVP for a DDE given by

where $\psi\in C:= C([-\tau, 0];\mathbb{R}^{d_{Y}})$, $G\colon C\to\mathbb{R}^{d_{Y}}$ and $y_{t}$ is defined as $x_{t}$ in Section 3.1 (so, again, $y_{0} = \psi$ represents the initial state). If $G$ is globally Lipschitz, the IVP (3.5) has a unique solution on $[-\tau, +\infty)$ [20, Section 2.2].

In ^[3] (but see also ^[4]) pseudospectral collocation is used to reduce (3.5) to an IVP for an ODE based on an equivalent formulation of (3.5) as an ACP describing the evolution of the original state $y_{t}$, viz.

through $v(t) = y_{t}$. Above, $\mathcal{A}_{Y}\colon\mathcal{D}(\mathcal{A}_{Y})\subset C\to C$ is defined as

In order to discretize (3.6), consider a mesh $\Omega_{M_{Y}, Y}$ of $1+M_{Y}$ points $-\tau = \theta_{M_{Y}, Y} < \theta_{M_{Y}-1, Y} < \dots < \theta_{1, Y} < \theta_{0, Y}:= 0$ with $M_{Y}$ a positive integer. Correspondingly, let $P_{M_{Y}, Y}\colon\mathbb{R}^{(1+M_{Y})d_{Y}}\to C$ be the interpolation operator on $\Omega_{M_{Y}, Y}$, i.e.,

where $\{\ell_{0, Y}, \ell_{1, Y}, \dots, \ell_{M_{Y}, Y}\}$ is the Lagrange basis on $\Omega_{M_{Y}, Y}$, and let $R_{M_{Y}, Y}:C\to\mathbb{R}^{(1+M_{Y})d_{Y}}$ be the restriction operator

Then, the discrete version of (3.6) is given by

for $V(t)\in\mathbb{R}^{(1+M_{Y})d_{Y}}$ that approximates the state $y_{t}$ according to $V_{j}(t)\approx(y_{t})(\theta_{j, Y})$, $j = 0, 1, \dots, M_{Y}$, and where $D_{M_{Y}, Y}$ has $d_{Y}\times d_{Y}$-block entries

where $I_{d_{Y}}$ is the identity on $\mathbb{R}^{d_{Y}}$ and $G_{M_{Y}}:= G\circ P_{M_{Y}, Y}$.

Remark 3.2. Note that the ACP (3.6) can be alternatively described as the equivalent semilinear ACP

where $\mathcal{A}_{0, Y}$ is the infinitesimal generator of the shift semigroup $\{T_{0, Y}(t)\}_{t\geq0}$ defined as

and $q_{Y}\in L^{\infty}$ is defined as

Equation (3.8) 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 ^[19]). The pseudospectral collocation of (3.8) leads, again, to (3.7), which can be rewritten equivalently as

where $D_{0, M_{Y}, Y}$ is as $D_{M_{Y}, Y}$ but with an additional $d_{Y}$-block row of zeros; also, $\mathbf{1}_{M_{Y}, Y}\in\mathbb{R}^{(1+M_{Y})d_{Y}\times d_{Y}}$ has the first $d_{Y}\times d_{Y}$-block equal to $I_{d_{Y}}$ and all of the others equal to zero. Now, (3.9) resembles (3.8).

3.3.   Pseudospectral collocation of coupled equations

Let $\tau$, $d_{X}$ and $d_{Y}$ be as above. Consider the IVP for a coupled equation given by

where $\phi\in L_{1}$, $\psi\in C$, $F\colon L^{1}\times C\to\mathbb{R}^{d_{X}}$ and $G\colon L^{1}\times C\to\mathbb{R}^{d_{Y}}$. For well-posedness, see ^[18].

By combining the approaches of the previous sections, it follows that (3.10) is equivalent to the ACP

with $\mathcal{B}_{X, Y}\colon\mathcal{D}(\mathcal{B}_{X, Y})\subset NBV_{0}\times Y\to NBV_{0}\times Y$ defined as

through $(u(t), v(t)) = (\mathcal{V}x_{t}, y_{t})$.

The discrete version of (3.11) is as follows:

Note that, now, $F_{M_{X}}:= F\circ(\mathcal{A}_{0, X}P_{M_{X}, X}, P_{M_{Y}, Y})$ and $G_{M_{Y}}:= G\circ(\mathcal{A}_{0, X}P_{M_{X}, X}, P_{M_{Y}, Y})$. The total number of approximating ODEs is $M_{X}d_{X}+(1+M_{Y})d_{Y}$, which becomes $(2M+1)d$ if $d_{X} = d_{Y} = d$ and $M_{X} = M_{Y} = M$.

Remark 3.3. Following Remarks 3.1 and 3.2, it is not difficult to see that (3.11) is equivalent to the semilinear ACP

where $\mathcal{A}_{0, X, Y}:= {\rm{diag}}(\mathcal{A}_{0, X}, \mathcal{A}_{0, Y})$ is linear and

is nonlinear. The pseudospectral collocation of (3.13) leads, again, to (3.12), where, correspondingly, the ODEs for $V$ can be compacted as done in Remark 3.2 for DDEs; see (3.9). The (numerical) analysis of (3.13) is current work in progress at CDLab⁷ , also in view of the corresponding sun–star theory of coupled equations developed in ^[18] and of the more recent twin semigroup theory of ^[1].

⁷http://cdlab.uniud.it/

4.   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⁸, where $\tau_2 > \tau_1 > 0$, and $f$, $f_1$, $f_2$ and $g$ are (possibly) nonlinear functions $\mathbb{R}\to\mathbb{R}$. Nevertheless, generalizing to other forms of equations is usually fairly straightforward⁹.

⁸Observe that the unknown of the differential equation is not delayed.

⁹For instance, one can consider models such as (3.3) in ^[21], yet with finite age-span, which do not enter class (4.2). 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 (3.10)), which was beyond the scope of this work.

We implemented the pseudospectral discretization¹⁰ using Chebyshev nodes of type Ⅱ (extrema) as the meshes $\{0\}\cup\Omega_{M_X, X}$ and $\Omega_{M_Y, Y}$ of points in $[-\tau, 0]$ with $M_X = M_Y$. To compute the nodes and the corresponding differentiation matrix, we used the $\texttt{cheb}$ routine of [22, Chapter 6]. For the interpolation, we used the barycentric Lagrange interpolation formula ^[23]; 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 ^[24,25]. We implemented $(\mathcal{A}_{0, X}P_{M_{X}, X}\Phi)(\theta)$ as $\sum_{j = 1}^{M_{X}}\ell_{j, X}'(\theta)\Phi_{j}$, computing $\ell_{j, X}'$ as the polynomial interpolating the $j$-th column of the differentiation matrix, again with the barycentric formula.

¹⁰For more details on pseudospectral methods, see also ^[22].

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 $\texttt{ode45}$, which implements the embedded Dormand–Prince $(5, 4)$ method ^[26,27]. For the differential part of (4.2), the initial value consists of the vector of values of the chosen initial function at $\Omega_{M_Y, Y}$, while, for (4.1) and the renewal part of (4.2), the vector $D_{M_X, X}^{-1}u$ is used¹¹, where $u$ is the vector of values of the initial function at $\Omega_{M_X, X}$.

¹¹In (3.4) the initial value is specified as $R_{M_{X}, X}\mathcal{V}\phi$, i.e., the vector representing the polynomial interpolating the exact integral of the initial value $\psi$. Another approach is to use $R_{M_{X}, X}\mathcal{V}P_{M_{X}, X}R_{M_{X}, X}\phi$, in which the integral of the polynomial interpolating $\psi$ 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 $\Omega_{M_X, X}$ is needed, as the value at $\theta_{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}\in\mathbb{R}^{d_{X}\times M_{X}d_{X}}$ is a row of $d_{X}\times d_{X}$-block entries $\ell'_{j, X}(\theta_{0, X})I_{d_{X}}$ for $j = 1, \dots, 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_{M_X, X}^{-1}u$.

Finally, the DQR for a linear ODE is implemented in the $\texttt{dqr}$ routine of ^[3], which follows ^[5]. Therein, the IVPs (2.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.

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 (4.1), 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 ^[28], which is based on the trapezoidal quadrature formula on a uniform grid in $[-\tau_2, 0]$ with the constraint that $-\tau_1$ must be a grid point. Corresponding to a solution $\bar{x}$, we considered the linear RE¹²

¹²In general, (4.3) may actually not be the linearization of (4.1) around $\bar{x}$ in $L^1$. Indeed, the right-hand side of the equation is not Fréchet-differentiable unless $f$ is affine. See [18, Section 3.5] for details, in particular with respect to studying the stability of equilibria; the extension of the results therein is an open problem.

See Section 5.1 for a comparison of the approaches.

Remark 4.2. In Remark 4.1 (more precisely in Footnote 12), 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 (4.1) and (4.2) can be linearized if $f$, $f_1$, $f_2$ and $g$ are differentiable. As an example, the linearization of (3.12) around a solution $(\overline{U}, \overline{V})$ is as follows:

where $J$ indicates the Jacobian matrix, $JF_{M_{X}}(\overline{U}(t), \overline{V}(t))$ is a $d_X \times (M_Xd_X+(1+M_Y)d_Y)$ matrix and $JG_{M_{Y}}(\overline{U}(t), \overline{V}(t))$ is a $d_Y \times (M_Xd_X+(1+M_Y)d_Y)$ matrix. 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.

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

5.1.   RE with quadratic nonlinearity

The first equation we study is the RE with quadratic nonlinearity from ^[7]:

i.e., (4.1) with $\tau_1:= 1$, $\tau_2:= 3$ and $f(x) := \frac{\gamma}{2}x(1-x)$. Its equilibria and their stability properties are known; in particular, its nontrivial equilibrium undergoes a Hopf bifurcation for $\gamma = 2 + \frac{\pi}{2}$ and the branch of periodic solutions arising from there has the analytic expression

Observe that the period is $4$, independent of $\gamma$. Moreover, it is experimentally known that it presents several period-doubling bifurcations, possibly leading to a cascade and, eventually, to chaos ^[7].

Since several properties of (5.1) 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 $\lambda$ of the infinitesimal generator of the semigroup of solution operators, which are related to the eigenvalues $\mu$ of the solution operator that advances the solution by $h$ via

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 (5.3), where $\mu$, $\lambda$ 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 $\mu$ of an evolution operator with any time step $h$ for the equilibria (we choose $h = \tau_2 = 3$ for (5.1)) and a time step $h$ equal to the period for periodic solutions ($h = 4$ for (5.2)), and then computing the real part of $\log(\mu)/h$. In order to obtain reference values for our experiments, we compute the spectra of evolution operators with the method of ^[29], which is based on the pseudospectral collocation of the operator; we use the implementation $\texttt{eigTMNpw}$ of ^[30,31].

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 (5.1). We choose $\gamma = 4 > 2+\frac{\pi}{2}$, which corresponds to a stable periodic solution, since the first period-doubling bifurcation is experimentally known to happen at $\gamma \approx 4.32$ ^[7,29].

Figure 1 shows the errors on the solution of the approximating ODE and on the solution of the original RE (5.1) 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^{13 [28]}; 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 $\texttt{ode45}$ with $\texttt{RelTol} = 10^{-6}$ and $\texttt{AbsTol} = 10^{-7}$, which justifies the barrier on the error in Figure 1.

¹³ As noted in Remark 4.1, $-\tau_1 = -1$ must be a grid point in $[-\tau_2, 0] = [-3, 0]$; we thus choose the number $r$ of nodes (minus $1$) in $[-\tau_1, 0]$ as the discretization parameter for the trapezoidal method ^[28], resulting in $3r+1$ nodes in total.

The experiment confirms that the trapezoidal method has order $2$, as proved in ^[28], and that the pseudospectral discretization has infinite order, which is often the case for pseudospectral methods applied to smooth problems ^[22]. 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¹⁴.

¹⁴As already noted, computing the solutions is not the focus of the present work. Admittedly, there are other more sophisticated methods in the literature: see, e.g., ^[32,33,34]. However, in most cases, they are not readily applicable to (4.1) and (4.2), due to the discontinuity in the integration kernel at $-\tau_1$, when the integral is considered on the interval $[-\tau_2, 0]$; in other cases, the implementation of the method is not available and is not as straightforward as ^[28]. It is worth mentioning that pseudospectral methods for computing periodic solutions of REs and coupled equations, exhibiting the usual infinite order of convergence, are available in ^[35,36,37]; they are based on solving the corresponding boundary value problem, so they cannot be straightforwardly adapted to computing generic solutions.

In the next experiment, we investigate how the errors on the LEs depend on the choice of $M_X$ and of the final time ¹⁵ $T$. We choose values of $\gamma$ corresponding to the stable trivial equilibrium ($\gamma = 0.5$), the stable nontrivial equilibrium ($\gamma = 3$) and the stable periodic orbit ($\gamma = 4$).

¹⁵Observe that the $\texttt{dqr}$ routine does not stop exactly at $T$, but at the first time step exceeding $T$, i.e., it does not refine the final step.

Since we are going to use the linearization of the ODE (3.4) coming from the RE (5.1), as an example, we show it explicitly here. With reference to Remark 4.2, observe that the right-hand side of (5.1) is not Fréchet-differentiable as a map from $L^1$ to $\mathbb{R}$, while the right-hand side of the discretized equation is differentiable. The linearization of the approximating ODE around the solution $\overline{U}$ is given by

where $JF_{M_{X}}(\overline{U}(t))$ is a row vector with components

In Figures 2 and 3, we can see the absolute errors on the dominant LE increasing either $M_X$ or $T$. The tolerance for $\texttt{dqr}$ is $10^{-6}$, while those for $\texttt{ode45}$ are $\texttt{RelTol} = 10^{-6}$ and $\texttt{AbsTol} = 10^{-7}$. The reference values are obtained by using $\texttt{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¹⁶. We 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 (2.5). Indeed, Figure 3, where $M_X$ is constant and $T$ varies, shows that the LEs converge linearly (confirming what is explained in ^[2]). In Figure 2 the truncation error appears to dominate on the error due to the collocation.

¹⁶Note in Figure 2 that for $\gamma = 3$ the error is initially below the final barrier; we do not have an explanation of this phenomenon, but in Figure 3 for $\gamma = 3$ and $M_X = 8$ we can observe rather erratic behavior of the error when varying $T$, despite the fact that it has a linearly decreasing bound. Moreover, we remark that the value of $M_X$ required to reach the error barrier for a given $T$ depends, in general, on the specific equation and the specific values of its parameters.

For the dominant nontrivial exponent¹⁷ 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 barriers are reached also for increasing values of $M_X$.

¹⁷The trivial exponent $0$ is always an LE for a periodic solution due to the translation invariance.

We showed here the results of the main approach only. We performed the same experiments also with the alternative approach of Remark 4.1, using, for the linearization, both the exact solutions (which are known for the chosen values of $\gamma$) and the numerical solutions computed by using the trapezoidal method. With the exact solutions, we obtained almost exactly the same values: this means that the solution of the ODE is a good enough approximation of the solution of the RE, and that exchanging the linearization and the collocation does not influence the results. However, when using the numerical solutions obtained via the trapezoidal method, the errors on the LEs were higher: in the experiment shown in Figure 2, the errors in the periodic case were one order of magnitude larger, while in the experiment of Figure 3 the errors reached barriers ranging between $10^{-2}$ and $10^{-1}$. For these reasons, we henceforth use only the main approach as the more practical one.

Figure 4 (compare with [7, Figure 2.3]) presents, on the top row, the diagram of the first two dominant LEs of (5.1) when varying $\gamma$, computed with $M_X = 15$ and $T = 1000$, following previous experimental considerations. We can observe that the dominant LE is $0$ at the expected Hopf bifurcation ($\gamma = 2+\frac{\pi}{2}$), after which one LE is always $0$ since periodic solutions appear. The dominant nontrivial exponent reaches $0$ again for the expected period-doubling bifurcations at $\gamma\approx 4.32, 4.49, 4.53$, and it becomes positive for $\gamma\geq4.55$, indicating the insurgence of a chaotic regime, which is compatible with what was obtained in ^[7]. Finally, for $\gamma \in [4.86, 4.9]$ other stability islands appear, corresponding to a branch of stable periodic solutions (appearing at $\gamma\approx4.8665$) and the corresponding cascade of period-doubling bifurcations (at $\gamma\approx4.8795, 4.8860$) leading back to chaos (starting at $\gamma\approx4.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.

5.2.   Egg cannibalism model

The second equation we consider is the egg cannibalism (toy) model from ^[38]:

with $a_{\text{mat}}$ and $a_{\text{max}}$ being, respectively, the constant maturation and maximum ages. Observe that (5.4) corresponds to (4.1) with $\tau_1:= a_{\text{mat}}$, $\tau_2:= a_{\text{max}}$ and $f(x) := \frac{\gamma}{2}x \mathrm{e}^{-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(\gamma) = 1+\frac{\pi}{2}$, although here the periodic solutions are not explicitly known; again, the presence of period-doubling bifurcations is experimentally known ^[4,7,8,38].

Figure 5 (top row) presents the diagram of the first two dominant LEs of (5.4) when varying $\gamma$, with $a_{\text{mat}} = 1$ and $a_{\text{max}} = 3$. The numerical parameters are $M_X = 15$, $T = 1000$, a tolerance of $10^{-6}$ for $\texttt{dqr}$, and $\texttt{RelTol} = 10^{-6}$ and $\texttt{AbsTol} = 10^{-7}$ for $\texttt{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(\gamma)\approx 3.855$ (^[8] finds $\log(\gamma)\approx3.8777$ with $M_X = 20$ and $\log(\gamma)\approx3.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(\gamma)\approx 4.54$, and it becomes positive for $\log(\gamma)\geq4.66$, indicating chaos. Figure 5 (bottom row) shows some stable periodic solutions, confirming the (approximate) doubling of the period.

5.3.   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 ^[38]:

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_{\text{mat}}$ and $a_{\text{max}}$ have the same meaning as in Section 5.2¹⁸. Equation (5.5) corresponds to (4.2) with $\tau_1:= a_{\text{mat}}$, $\tau_2:= a_{\text{max}}$, $f_1(x) := \beta x$, $f_2(x) := -\gamma x$ and $g(y) := ry(1-\frac{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 $\beta = (K(a_{\text{max}}-a_{\text{mat}}))^{-1}$; moreover, when varying $\beta$, the nontrivial equilibrium is experimentally known to undergo a Hopf bifurcation ^[4,38,39].

¹⁸Note that the second term on the right-hand side of the equation for $S$ in (5.5) can be rewritten as $\gamma b(t)/\beta$, rendering that equation an ODE. This does not bring any simplification since the integral term is computed for the first equation anyway. We prefer to keep the form of (5.5) because it comes from a more general class where the fertility and consumption rates are not constant, but, rather, are functions of the age and the size of the individuals (see (1) in ^[38]).

The diagram of the first two dominant LEs of (5.5) when varying $\beta$ is shown in Figure 6. The values of the other model parameters are $a_{\text{mat}} = 3$, $a_{\text{max}} = 4$ and $r = K = \gamma = 1$. The numerical parameters are $M_X = M_Y = 15$, $T = 1000$, a tolerance of $10^{-6}$ for $\texttt{dqr}$, and $\texttt{RelTol} = 10^{-4}$ and $\texttt{AbsTol} = 10^{-5}$ for $\texttt{ode45}$. We can observe a spike at $\beta = 1$ (albeit not touching the value $0$), corresponding to the known transcritical bifurcation, while the Hopf bifurcation seems to happen for $\gamma\approx3$ (^[4] finds $\gamma\approx3.0161$ with $M_X = 10$ and MatCont's tolerances set to $10^{-10}$). We continued the experiment for values of $\beta$ 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 $\texttt{ode45}$ in order to improve the computation times; however, for $\beta = 1$ the dominant LE differs absolutely from the one computed with $\texttt{RelTol} = 10^{-6}$ and $\texttt{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 $\beta = 5$ the absolute difference is larger but still less than $10^{-5}$). Moreover, as $\beta$ increases, the periodic solution presents flat regions followed by spikes, which may suggest that the equation is becoming stiff; however, we tried replacing $\texttt{ode45}$ with MATLAB's $\texttt{ode23s}$, implementing a modified Rosenbrock formula of order $2$ for stiff ODEs ^[27], with no substantial changes.

6.   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 ^[2], 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 (2.3).

The next step in our research is to tackle the problem following the approach of ^[2]; 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 ^[3] 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 ^[40]), as anticipated in ^[3] and according to the pragmatic point of view discussed in ^[7].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors are members of INdAM research group GNCS and of UMI research group "Modellistica socio-epidemiologica". This work was partially supported by the Italian Ministry of University and Research (MUR) through the PRIN 2020 project (No. 2020JLWP23) "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (CUP: E15F21005420006), and by the GNCS 2023 project "Sistemi dinamici e modelli di evoluzione: tecniche funzionali, analisi qualitativa e metodi numerici" (CUP: E53C22001930001).

Conflict of interest

The authors declare there is no conflict of interest.