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

TL;DR

The paper tackles the challenge of computing Lyapunov exponents for renewal (Volterra-type delay) equations and their couplings with delay differential equations. It introduces a practical pipeline that reformulates REs as abstract differential equations, discretizes via pseudospectral collocation to yielding finite-dimensional ODEs, and then applies the standard discrete QR method to compute the dominant exponents, with a MATLAB implementation provided. The authors validate the approach on three renewal/delay models, recovering known bifurcation structure (Hopf, period-doubling) and chaotic regimes, and compare against alternative discretizations to demonstrate accuracy and efficiency. The work fills a gap in LE computation for REs, broadening numerical analysis tools for delay systems and enabling quantitative stability and complexity assessments in renewal-type dynamics. The methodology is poised to impact analyses in mathematical biology and other fields where renewal-type delays arise, offering a practical, implementable route to LE estimation beyond classical DDEs.

Abstract

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.

Figures (6)

• Figure 1: Errors on the solution of the RE with quadratic nonlinearity \ref{['quadRE']} with $\gamma=4$ with respect to the known periodic solution \ref{['quadREpersol']}, computed via pseudospectral discretization (solid lines) and directly with the trapezoidal method (dashed lines), measured as the absolute error at $t=500$ ($\bullet$) and as the maximum absolute error on a grid of points in $[0, 500]$ ($\circ$), 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 \ref{['quadREpersol']} is used as the initial value.
• Figure 2: Absolute errors on the dominant LEs of the RE with quadratic nonlinearity \ref{['quadRE']} for 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$). 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 final time for dqr is $T=1000$.
• Figure 3: Absolute errors on the dominant LEs of the RE with quadratic nonlinearity \ref{['quadRE']} for 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$). 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 \in \{8, 12, 16, 20\}$ as shown in the legend.
• Figure 4: Diagram of the first two dominant (in descending order) LEs (top row) and solutions (other rows) of the RE with quadratic nonlinearity \ref{['quadRE']} when varying $\gamma$, 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 $\gamma$ and of the period are given above the diagrams; the values of $\gamma$ are also marked by vertical dashed lines in the diagram of the LEs.
• Figure 5: Diagram of the first two dominant (in descending order) LEs (top row) and solutions (bottom row) of the egg cannibalism RE \ref{['cannibRE']} with $a_{\text{mat}}=1$ and $a_{\text{max}}=3$, when varying $\gamma$, 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 $\gamma$ and of the period are given above the diagrams; the values of $\gamma$ are also marked by vertical dashed lines in the diagram of LEs.

Theorems & Definitions (5)

• Remark 3.1
• Remark 3.2
• Remark 3.3
• Remark 4.1
• Remark 4.2