remarkRemark hypothesisHypothesis claimClaim Equations with infinite delay:
pseudospectral approximationF. Scarabel and R. Vermiglio

# Equations with infinite delay:
pseudospectral discretization for numerical stability and bifurcation in an
abstract framework††thanks: Funding: The work of RS is supported by the
Italian Ministry of University and Research (MUR) through the PRIN 2020
project (No. 2020JLWP23) “Integrated Mathematical Approaches to
Socio–Epidemiological Dynamics”, Unit of Udine (CUP: G25F22000430006).

Francesca Scarabel Department of Applied Mathematics, University of Leeds,
Woodhouse, Leeds LS2 9JT, UK<EMAIL_ADDRESS>Rossana Vermiglio
Department of Mathematics and Computer Science, University of Udine, via delle
Scienze 208, 33100 Udine, Italy<EMAIL_ADDRESS>
###### Abstract

We consider nonlinear delay differential and renewal equations with infinite
delay. We extend the work of Gyllenberg _et al_ , _Appl. Math. Comput._ (2018)
by introducing a unifying abstract framework, and derive a finite-dimensional
approximating system via pseudospectral discretization. For renewal equations,
we consider a reformulation in the space of absolutely continuous functions
via integration. We prove the one-to-one correspondence of equilibria between
the original equation and its approximation, and that linearization and
discretization commute. Our most important result is the proof of convergence
of the characteristic roots of the pseudospectral approximation of the
linear(ized) equations when the collocation nodes are chosen as the family of
scaled zeros or extrema of Laguerre polynomials. This ensures that the finite-
dimensional system correctly reproduces the stability properties of the
original linear equation if the dimension of the approximation is large
enough. The result is illustrated with several numerical tests, which also
demonstrate the effectiveness of the approach for the bifurcation analysis of
equilibria of nonlinear equations. The new approach used to prove convergence
also provides the exact location of the spectrum of the differentiation
matrices for the Laguerre zeros and extrema, adding new insights into
properties that are important in the numerical solution of differential
equations by pseudospectral methods.

###### keywords:

Renewal equations, delay differential equations, spectral collocation,
exponentially weighted interpolation, abstract differential equation, linear
stability, differentiation matrix.

34K99, 37M20, 65J99, 65L07, 65P30, 65R20, 92D25

##  1 Introduction

According to a standard definition, “a delay equation is a rule for extending
a function of time towards the future on the basis of the (assumed to be)
known past” [20]. Delay equations (DEs) include renewal equations (REs), which
specify the value of the unknown function in terms of its past values, delay
differential equations (DDEs), which specify the derivative, and systems
coupling both [13, 19, 20].

In this paper we focus on equations with infinite delay (iDEs, including iDDEs
and iREs), which are widely used in mathematical biology, and in particular in
ecology and epidemiology, to describe physiologically structured population
models, in which the individual rates are assumed to depend on a structuring
variable (e.g., age, size, time since infection) that evolves in time [15, 16,
17, 34, 38]. In behavioral epidemic models, the infinite delay can arise to
incorporate the effect of information and waning memory [21, 31]. Other
applications of iDDEs can be found in [30] and the references therein.

Our interest focuses on the long-term behavior of the systems, and in
particular on the existence and stability of equilibria and their bifurcation
analysis (i.e., how these properties change under parameter variation).
Differently from ordinary differential equations (ODEs), DEs generate
infinite-dimensional dynamical systems where the state space is a suitable
function space, hence substantially increasing the complexity of their
theoretical and numerical analysis. One way to overcome theoretically the
infinite dimensionality of the problem is via reduction to finite-dimensional
inertial or center manifolds [23, 28]. However in this paper we take a
different approach, aiming at studying _numerically_ the infinite-dimensional
system by reducing it to a finite-dimensional system of ODEs, whose properties
can be studied with well-established software, see e.g. [5] and references
therein. While in some special cases this reduction can be done exactly, e.g.,
via the linear chain trick [9, 15, 16], for more general kernels one usually
relies on approximations.

Among the various techniques available for this reduction [3, 25, 39],
pseudospectral discretization (PSD) has been successfully applied in the last
decade both to linear(ized) problems for the stability analysis of equilibria,
and to nonlinear ones for the numerical bifurcation in various contexts,
including DEs with bounded delay and PDEs with nonlocal boundary conditions
[1, 4, 5, 6, 7, 18, 10, 11, 36, 37]. For DEs, the main advantage of using PSD
compared to other approximation techniques is the fast (exponential) order of
convergence for the stability of the equilibria, which follows from the fact
that the eigenfunctions of the infinitesimal generator of the linear(ized)
equation are exponential, hence analytic.

To treat the unbounded delay, several approaches are used in the literature,
including domain truncation or suitable change of variables [3, 25, 39], or
even quadrature rules to approximate an iDDE with a DDE with a finite number
of point delays to be studied via established methods (e.g. [5]) or software
like DDE-BIFTOOL [22]. However, these approaches do not exploit some of the
fundamental features of iDEs, in particular the “exponentially fading memory”,
which is the major advantage of the PSD proposed here. In fact, the step from
finite to infinite delay involves several complexities, both from a
theoretical and numerical point of view.

From a theoretical point of view, for iDEs the principle of linearized
stability (PLS) has been proved in a space of exponentially weighted functions
defined on $\mathbb{R}_{-}:=(-\infty,0]$, with weight

(1) $w(\theta):=\mathrm{e}^{\rho\theta},\quad\theta\in\mathbb{R}_{-},$

for some $\rho>0$. Using a modification of the sun-star calculus, Diekmann and
Gyllenberg show that the stability of equilibria of a nonlinear equation is
determined by the part of the spectrum of the infinitesimal generator of the
linearized semigroup in the right-half plane

$\mathbb{C}_{\rho}:=\\{\lambda\in\mathbb{C}\ |\ \Re\lambda>-\rho\\},$

which contains eigenvalues only, that are roots of a characteristic equation
[14, Theorem 3.15 and Theorem 4.7]. Recently, the PLS and the Hopf bifurcation
theorem have been proved for iDDEs by applying the integrated semigroup theory
in [30]. We also refer to [40] for a Hopf result in a space of bounded
uniformly continuous functions. However, in this paper we keep reference to
[14] as it covers iDDEs, iREs and coupled iDDEs and iREs. Finally, we note
that other weight functions can be considered [29], but here we restrict to
the exponential weight (1) for which the PLS has been rigorously proved.

In this theoretical framework and motivated by the accuracy, effectiveness and
flexibility shown by the PSD approach for bounded delay, a detailed
experimental analysis has been carried out using collocation on the scaled
Laguerre nodes in [27] for iDDEs.

Here, we introduce a unifying abstract differential equation (ADE) able to
capture iDDEs and iREs within the same framework. This pragmatic
representation allows us to encompass all the fundamental concepts of iDEs and
the general derivation of the finite-dimensional dynamical system through the
weighted version of the PSD all at once. For iREs, we also propose a new
approach inspired by [20, 37]: via integration, we map the original iRE,
defined in a space of weighted integrable functions, to an equivalent equation
defined in a subspace of absolutely continuous functions on any bounded
subinterval of $\mathbb{R}_{-}$. The main advantages are twofold: point
evaluation, as well as polynomial interpolation and collocation, are well
defined in the new space, and the domain of the infinitesimal generator is
defined by a trivial condition, simplifying the implementation.

Most importantly, here we add the key element to the numerical investigations
in [27]: the rigorous proof that the PSD accurately approximates the stability
properties of the equilibria. To this aim, we first show that there is a one-
to-one correspondence between the equilibria of the ADE and its PSD. Then, we
focus on the linear(ized) problems and prove the convergence of the stability
indicators for the zero equilibrium, i.e., the eigenvalues of the
infinitesimal generator in $\mathbb{C}_{\rho}$. Due to the complexities
introduced by the step from finite to infinite delay and in the absence of
suitable interpolation results on the semi-unbounded interval, the proof
follows a slightly different route compared to the standard ones for finite
delay [7, 37]. The new approach guarantees the convergence of the approximated
eigenvalues by using as collocation nodes the zeros or extrema of standard
Laguerre polynomials suitably scaled. Moreover, the general results for the
Laguerre polynomials presented in the appendix provide an exact description of
the location of the eigenvalues of the differentiation matrix with boundary
condition built on the scaled Laguerre zeros and extrema.

The paper is organized as follows. In Section 2 we recall the basic results of
stability theory for iDEs and introduce their abstract formulation. The
unifying abstract framework is introduced in Section 3, together with the PSD
approach. In Section 4, after deriving the characteristic equations, we show
the one-to-one correspondence of equilibria of the original ADE and its
finite-dimensional approximation, and prove the convergence of the
characteristic roots. To improve readability, the results used in the proofs
are collected in Appendix A for the zeros and extrema of the standard Laguerre
polynomials in $\mathbb{R}_{+}$. Finally, Section 5 illustrates by numerical
experiments the theoretical convergence of the eigenvalues for linear iDEs and
the flexibility of the approach for numerical bifurcation analysis of
equilibria of nonlinear iDEs.

##  2 Abstract formulation of nonlinear iDEs

In this section we recall from [14] some basic results for iDDEs and iREs,
including the PLS for the stability of equilibria. In view of Section 3, we
also show how both iDDEs and iREs (and, by combining the two approaches,
systems coupling both) can be formally written as semilinear ADEs. Although we
typically consider real-valued nonlinear equations and real-valued solutions,
to apply spectral theory we assume that real Banach spaces are embedded into
complex ones and all operators are extended to complex spaces. For further
details on the complexification process we refer to [19, Chapter III.7].

Let $\rho>0$ and consider the weight function (1). For $d\in\mathbb{N}$,
$d>0$, we denote with $|\cdot|$ a norm in $\mathbb{C}^{d}$. Let
$L^{\infty}_{\rho}:=L^{\infty}_{\rho}(\mathbb{R}_{-};\mathbb{C}^{d})$ be the
space of (the equivalence classes of) functions $\varphi$ such that
$|w\varphi|$ is essentially bounded in $\mathbb{R}_{-}$, and
$L^{p}_{\rho}:=L^{p}_{\rho}(\mathbb{R}_{-};\mathbb{C}^{d})$, $1\leq
p<+\infty$, that of $\varphi$ such that $w|\varphi|^{p}$ is integrable,
equipped with the norm

$\|\varphi\|_{p,\rho}:=\|w\varphi\|_{p}=\begin{cases}\left(\int_{-\infty}^{0}w(\theta)|\varphi(\theta)|^{p}\mathop{}\\!\mathrm{d}\theta\right)^{\frac{1}{p}},&1\leq
p<+\infty,\\\\[5.0pt] {\rm
ess}\sup_{\theta\in\mathbb{R}_{-}}w(\theta)|\varphi(\theta)|,&p=\infty.\end{cases}$

Note that $wL^{p}_{\rho}\cong L^{p}$, but we will often write $w$ explicitly
to keep reference to the original state space of the dynamical system.

To ease notation, we write $w\psi(\theta)$ in place of
$w(\theta)\psi(\theta)$. We consider column vectors and use the convention of
writing them as rows, but separating each column block-vector elements with
semicolons. We will also often abuse notation by using a symbol to denote both
a constant function and the value it takes on.

###  2.1 Delay differential equations

A $d$-dimensional nonlinear iDDE is

(2) $\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}(t)=F(y_{t}),\quad
t\geq 0,$

where $y_{t}(\theta):=y(t+\theta)$, $\theta\in\mathbb{R}_{-}$, is the “history
function” at time $t\geq 0.$ We assume that
$F\colon\mathcal{C}_{\rho}\to\mathbb{R}^{d}$ is a continuously Fréchet
differentiable map, where
$\mathcal{C}_{\rho}:=\mathcal{C}_{0,\rho}(\mathbb{R}_{-};\mathbb{R}^{d})$ is
the Banach space of all functions $\psi\colon\mathbb{R}_{-}\to\mathbb{R}^{d}$
such that $w\psi$ is continuous and vanishes at $-\infty$, equipped with the
norm $\|\cdot\|_{\infty,\rho}$.

To represent (2) in the abstract framework, we consider the space

$X:=w\mathbb{C}^{d}\times wL_{\rho}^{\infty}$

equipped with the norm
$\|(w\alpha;w\psi)\|_{X}=\max\\{|\alpha|,\|w\psi\|_{\infty}\\}.$ By defining
the embedding $j\colon\mathcal{C}_{\rho}\to X$ as $j\psi:=(w\psi(0);w\psi)$,
we identify $\mathcal{C}_{\rho}$ with the subspace $Y:=j(\mathcal{C}_{\rho})$
of $X$. Note that, by interpreting $\alpha\in\mathbb{C}^{d}$ as a constant
function in $L_{\rho}^{\infty},$ we have that $X\subset
wL_{\rho}^{\infty}(\mathbb{R}_{-},\mathbb{C}^{\bar{d}})$ with $\bar{d}=2d$.

By defining $u(t):=jy_{t},$ it is easy to verify that $u$ satisfies the
semilinear ADE

(3)
$\frac{\mathop{}\\!\mathrm{d}u}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0}^{\text{DDE}}u(t)+\mathcal{F}^{\text{DDE}}(u(t)),\quad
t\geq 0,$

where the linear unbounded operator $\mathcal{A}_{0}^{\text{DDE}}\colon
D(\mathcal{A}_{0}^{\text{DDE}})\to X$ is

(4)
$\begin{array}[]{l}\mathcal{A}_{0}^{\text{DDE}}(w\alpha;w\psi)=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\theta}(w\alpha;w\psi)-\rho(w\alpha;w\psi)=(0;w\frac{\mathop{}\\!\mathrm{d}\psi}{\mathop{}\\!\mathrm{d}\theta}),\\\\[8.0pt]
D(\mathcal{A}_{0}^{\text{DDE}})=\\{(w\alpha;w\psi)\mid\psi\text{ is Lipschitz
continuous, }\psi(0)=\alpha\\},\end{array}$

and the nonlinear operator $\mathcal{F}^{\text{DDE}}\colon Y\to X$ is

(5)
$\mathcal{F}^{\text{DDE}}(w\alpha;w\psi)=(wF(j^{-1}(w\alpha,w\psi));0)=(wF(\psi);0).$

Here, $\mathcal{A}_{0}^{\text{DDE}}$ captures the “shift” of the solution,
$\mathcal{F}^{\text{DDE}}$ captures the “extension” according to the rule
defined by (2), while both account for “weight multiplication”. Note that the
$\theta$-derivative in (4) should be interpreted as almost everywhere.

Now we consider the nonlinear operator $\mathcal{F}^{\text{DDE}}.$ The
canonical basis $e_{1},\dots,e_{d}$ of $\mathbb{R}^{d}$ determines $d$
linearly independent vectors $\xi_{i}=(e_{i};0)$, $i=1,\dots,d,$ and
$\Xi_{d}:=\text{span}\\{w\xi_{i}\mid i=1,...,d\\}\subset X.$ By defining
$\mathcal{F}_{i}^{\text{DDE}}(w\alpha;w\psi)=F_{i}(\psi)$, we get the finite-
rank representation
$\mathcal{F}^{\text{DDE}}=w\sum_{i=1}^{d}\mathcal{F}_{i}^{\text{DDE}}\xi_{i}$.
Note that $\mathcal{F}^{\text{DDE}}$ inherits the regularity properties of
$F$: if $F$ is continuously Fréchet differentiable, so is
$\mathcal{F}^{\text{DDE}}$.

We are interested in the stability properties of an equilibrium $\bar{\psi}$
of (2), which is a constant solution $\bar{\psi}$ such that $F(\bar{\psi})=0$.
The linearization of (2) at $\bar{\psi}$ is

(6)
$\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}(t)=DF(\bar{\psi})y_{t},\quad
t\geq 0,$

where $DF(\bar{\psi})$ is the Fréchet derivative of $F$ at $\bar{\psi}$.

From (4)–(5) we get that $\bar{u}=(w\bar{\psi};w\bar{\psi})$ is an equilibrium
of (3) and vice versa. Moreover by introducing the linear operator
$\mathcal{L}^{\text{DDE}}(w\alpha;w\psi)=(wDF(\bar{\psi})\psi;0)$, the
abstract formulation of (6) is

(7)
$\frac{\mathop{}\\!\mathrm{d}u}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0}^{\text{DDE}}u(t)+\mathcal{L}^{\text{DDE}}u(t),\quad
t\geq 0.$

By the PLS [14, Theorem 4.7], the stability of $\bar{\psi}$ is determined by
the stability of the zero solution of (6) as follows.

######  Theorem 2.1.

Let $\bar{\psi}$ be an equilibrium of (2) and
$\bar{u}=(w\bar{\psi};w\bar{\psi})$ the corresponding equilibrium of (3). Then

  * •

if all the eigenvalues of
$\mathcal{A}^{\text{DDE}}:=\mathcal{A}_{0}^{\text{DDE}}+\mathcal{L}^{\text{DDE}}$
in $\mathbb{C}_{\rho}$ have negative real part, then $\bar{u}$ is
exponentially stable;

  * •

if at least one eigenvalue of $\mathcal{A}^{\text{DDE}}$ has positive real
part, then $\bar{u}$ is unstable;

  * •

the eigenvalues of $\mathcal{A}^{\text{DDE}}$ in $\mathbb{C}_{\rho}$
correspond to roots of the characteristic equation

$\det(\lambda I_{d}-DF(\bar{\psi})\mathrm{e}^{\lambda\cdot})=0$

in $\mathbb{C}_{\rho}$, which are at most finitely many.

Given an eigenvalue $\lambda\in\mathbb{C}_{\rho}$ of
$\mathcal{A}^{\text{DDE}}$, its eigenfunction is of the form
$(w\alpha;w\psi_{\lambda})$, where
$\psi_{\lambda}(\theta):=\mathrm{e}^{\lambda\theta}\,\alpha$, for
$\alpha\in\mathbb{C}^{d}\setminus\\{0\\}$.

###  2.2 Renewal equations

Consider now the $d$-dimensional nonlinear iRE

(8) $y(t)=F(y_{t}),\quad t>0,$

for $y_{t}\in L^{1}_{\rho}$, where $F\colon L^{1}_{\rho}\to\mathbb{R}^{d}$ is
a continuously Fréchet differentiable map. In particular, this assumption
means that we only consider equations (8) such that their linearization is a
linear Volterra integral equations, ruling out neutral REs with point delays.

To derive the abstract equation, we define the space $X:=wL^{1}_{\rho}$ and
the embedding $j:L^{1}_{\rho}\to X$ given by
$j\varphi(\theta):=w\int_{0}^{\theta}\varphi(s)\mathop{}\\!\mathrm{d}s$, for
$\theta\in\mathbb{R}_{-}$. Note that

$\|j\varphi\|_{X}\leq\int_{-\infty}^{0}\mathrm{e}^{\rho\theta}\int_{\theta}^{0}|\varphi(s)|\mathop{}\\!\mathrm{d}s\mathop{}\\!\mathrm{d}\theta=\int_{-\infty}^{0}\int_{-\infty}^{s}\mathrm{e}^{\rho\theta}\mathop{}\\!\mathrm{d}\theta\
|\varphi(s)|\mathop{}\\!\mathrm{d}s\leq\frac{\|w\varphi\|_{1}}{\rho},$

and $Y:=jL^{1}_{\rho}=\\{w\varphi\in X\mid
w\varphi=j\frac{\mathop{}\\!\mathrm{d}\varphi}{\mathop{}\\!\mathrm{d}\theta}\text{
and }w\frac{\mathop{}\\!\mathrm{d}\varphi}{\mathop{}\\!\mathrm{d}\theta}\in
X\\}$ coincides with the subspace
$AC_{0}:=AC_{0}(\mathbb{R}_{-};\mathbb{C}^{d})$ of all absolutely continuous
functions $w\varphi\in X$ in any closed subset of $\mathbb{R}_{-}$, such that
$w\varphi(0)=0.$ Note that, for $w\varphi\in Y$,
$j^{-1}(w\varphi)=\frac{\mathop{}\\!\mathrm{d}\varphi}{\mathop{}\\!\mathrm{d}\theta}$.

By defining the state $u(t):=jy_{t}$, we can formally derive the ADE

(9)
$\frac{\mathop{}\\!\mathrm{d}u}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0}^{\text{RE}}u(t)+\mathcal{F}^{\text{RE}}(u(t)),\quad
t\geq 0,$

where

(10)
$\displaystyle\mathcal{A}_{0}^{\text{RE}}(w\varphi)=\textstyle{\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\theta}}(w\varphi)-\rho
w\varphi=w\textstyle{\frac{\mathop{}\\!\mathrm{d}\varphi}{\mathop{}\\!\mathrm{d}\theta}},$
$\displaystyle w\varphi\in D(\mathcal{A}_{0}^{\text{RE}})=Y,$ (11)
$\displaystyle\mathcal{F}^{\text{RE}}(w\varphi)=-wF(j^{-1}w\varphi)=-wF(\textstyle{\frac{\mathop{}\\!\mathrm{d}\varphi}{\mathop{}\\!\mathrm{d}\theta}}),$
$\displaystyle w\varphi\in Y.$

Similarly as before, $\mathcal{F}^{\text{RE}}$ has finite rank with range
$\Xi_{d}=w\mathbb{R}^{d}.$

As expected, there is a one-to-one correspondence between the equilibria
$\bar{\varphi}$ of (8), which are constant functions such that
$\bar{\varphi}=F(\bar{\varphi}),$ and the equilibria $\bar{u}=j\bar{\varphi}$
of (9). Moreover, by defining
$\mathcal{L}^{\text{RE}}:=-wDF(\bar{\varphi})j^{-1},$ the abstract formulation
of the linearized RE $y(t)=DF(\bar{\varphi})y_{t}$, $t>0$, gives the
linearized equation

(12)
$\frac{\mathop{}\\!\mathrm{d}u}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0}^{\text{RE}}u(t)+\mathcal{L}^{\text{RE}}u(t),\quad
t\geq 0.$

The following PLS holds [14, Theorem 3.9].

######  Theorem 2.2.

Let $\bar{\varphi}$ be an equilibrium of (8) and $\bar{u}=j\bar{\varphi}$ the
corresponding equilibrium of (11). Then

  * •

if all the eigenvalues of
$\mathcal{A}^{\text{RE}}:=\mathcal{A}_{0}^{\text{RE}}+\mathcal{L}^{\text{RE}}$
in $\mathbb{C}_{\rho}$ have negative real part, then $\bar{u}$ is
exponentially stable;

  * •

if at least one eigenvalue of $\mathcal{A}^{\text{RE}}$ has positive real
part, then $\bar{u}$ is unstable;

  * •

the eigenvalues of $\mathcal{A}^{\text{RE}}$ in $\mathbb{C}_{\rho}$ correspond
to zeros of characteristic equation

$\det(I_{d}-DF(\bar{\varphi})\mathrm{e}^{\lambda\cdot})=0$

in $\mathbb{C}_{\rho},$ which are at most finitely many.

Given an eigenvalue $\lambda\in\mathbb{C}_{\rho}$ of
$\mathcal{A}^{\text{RE}}$, the corresponding eigenfunction is
$w\varphi_{\lambda}$, where
$\varphi_{\lambda}(\theta):=\frac{\mathrm{e}^{\lambda\cdot}-1}{\lambda}\alpha$,
for $\alpha\in\mathbb{C}^{d}\setminus\\{0\\}$. Note also that
$j^{-1}(w\varphi_{\lambda})=\mathrm{e}^{\lambda\cdot}\alpha$.

##  3 A unifying abstract framework

We now propose a unifying abstract framework which is particularly useful to
derive the approximating finite-dimensional system via PSD and to analyse the
convergence in Section 4. We consider a general framework where $X$ is a
Banach space of $\mathbb{C}^{\bar{d}}$-valued functions defined on
$\mathbb{R}_{-}$ equipped with a non-weighted $p$-norm, $1\leq p\leq+\infty,$
denoted by $\|\cdot\|_{X}.$ In other words, $X$ is a closed subspace of
$L^{p}(\mathbb{R}_{-};\mathbb{C}^{\bar{d}}).$

As before, let $\rho>0$ and $w$ be the function (1). Consider the semilinear
ADE

(13)
$\frac{\mathop{}\\!\mathrm{d}u}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0}u(t)+\mathcal{F}(u(t)),\quad
t\geq 0,$

where $\mathcal{A}_{0}\colon D(\mathcal{A}_{0})\to X$ is the linear unbounded
operator

(14)
$\begin{array}[]{l}\mathcal{A}_{0}x=\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}\theta}-\rho
x=w\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\theta}(\frac{x}{w}),\quad
x\in D(\mathcal{A}_{0}),\\\\[8.0pt] D(\mathcal{A}_{0})=\\{x\in
X\mid\mathcal{A}_{0}x\in X\text{ and }x(0)=\mathcal{B}x\\},\end{array}$

with $\mathcal{B}\colon X\to\mathbb{R}^{\bar{d}}$ linear and bounded.

######  Remark 3.1.

The operators $\mathcal{A}_{0}^{\text{DDE}}$ and $\mathcal{A}_{0}^{\text{RE}}$
defined in (4) and (10) are particular instances of (14) with
$\mathcal{B}(w\alpha,w\psi)=(\alpha;\alpha)$ and $\mathcal{B}(w\varphi)=0$,
respectively. In particular, the form of $\mathcal{B}$ characterizes the type
of iDE under investigation.

Note that $D(\mathcal{A}_{0})$ is a subset of the Sobolev space
$W^{p,1}(\mathbb{R}_{-};\mathbb{C}^{\bar{d}}),$ hence it admits one (and only
one) continuous representative [8]. We assume that $\mathcal{F}\colon
D(\mathcal{A}_{0})\to X$ is continuously Fréchet differentiable and admits the
finite-rank representation

(15) $\mathcal{F}(x)=w\sum_{i=1}^{d}\mathcal{F}_{i}(x)\xi_{i}\in\Xi_{d}$

for some $\mathcal{F}_{i}\colon D(\mathcal{A}_{0})\to\mathbb{R}$ sufficiently
smooth, $\xi_{i}$, $i=1,\dots,d$, linearly independent vectors in
$\mathbb{C}^{\bar{d}}$, and $\Xi_{d}:=w\,\text{span}\\{\xi_{i}\mid
i=1,...,d\\}.$ The continuity of the elements of $D(\mathcal{A}_{0})$, where
we assume that $u(t)$ lives for $t\geq 0$, justifies the choice of the PSD
approach, which requires point evaluation to be well defined.

###  3.1 Reduction to finite dimension via PSD

Given a positive integer $N$, called _discretization index_ , we consider a
set of $N$ distinct points in $\mathbb{R}_{-}$,

$\Theta_{N}:=\\{\theta_{1},\dots,\theta_{N}\mid\theta_{N}<\dots<\theta_{1}<0\\}.$

We refer to the set $\Theta_{N}$ as the _collocation nodes_. Then, we
represent every $x\in D(\mathcal{A}_{0})$ by a set of vectors
$X_{j}=x(\theta_{j})$, $j=1,\dots,N$ via the restriction operator

$\mathcal{R}_{N}x=(X_{1};\dots;X_{N})\in\mathbb{C}^{N\bar{d}}.$

The operator $\mathcal{R}_{N}$ is complemented with a prolongation operator
$\mathcal{P}_{N,0}$ that maps $\mathbb{R}^{N\bar{d}}$ to the space
$D(\mathcal{A}_{0})\cap\Pi_{N}^{\bar{d}}$, where $\Pi_{N}^{\bar{d}}$ is the
space of $\mathbb{C}^{\bar{d}}$-valued polynomials of degree at most $N$.
$\mathcal{P}_{N,0}$ is defined such that
$x_{N,0}=\mathcal{P}_{N,0}(X_{1};\dots;X_{N})$ satisfies the implicit
condition $x_{N,0}(0)=\mathcal{B}x_{N,0}$. In practice, we first define the
extended mesh

$\Theta_{N,0}:=\Theta_{N}\cup\\{\theta_{0}=0\\}$

and consider the $N+1$ Lagrange polynomials $\ell_{i,0}$, $i=0,\dots,N$,
associated with $\Theta_{N,0}$. Then, we take

$x_{N,0}=\beta\,\ell_{0,0}+\sum_{i=1}^{N}X_{i}\ell_{i,0},$

with $\beta$ implicitly defined by $x_{N,0}(0)=\mathcal{B}x_{N,0}$. With these
choices, $\mathcal{A}_{0}\mathcal{P}_{N,0}$ maps into $w\Pi_{N-1}^{\bar{d}}$
and the application of $\mathcal{R}_{N}$ is well defined.

Using $\mathcal{R}_{N}$ and $\mathcal{P}_{N,0}$ we can derive the discrete
versions of $\mathcal{A}_{0}$ and $\mathcal{F}$, namely

$\mathcal{A}_{0,N}:=\mathcal{R}_{N}\mathcal{A}_{0}\mathcal{P}_{N,0},\quad\mathcal{F}_{N}:=\mathcal{R}_{N}\mathcal{F}\circ\mathcal{P}_{N,0},$

where $\circ$ denotes composition of nonlinear operators. Finally, the PSD of
(13) is the semilinear ODE

(16)
$\frac{\mathop{}\\!\mathrm{d}U}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0,N}U(t)+\mathcal{F}_{N}(U(t)),\quad
t\geq 0,$

for $U(t)\in\mathbb{C}^{N\bar{d}}$. In the next section, we study how the
equilibria and stability indicators of (13) are approximated by the
counterparts of (16).

##  4 Approximation of equilibria and their stability

An equilibrium $\bar{u}\in D(\mathcal{A}_{0})$ of (13) is a solution of
$\mathcal{A}_{0}\bar{u}+\mathcal{F}(\bar{u})=0$, and, for
$\mathcal{L}:=D\mathcal{F}(\bar{u})$, the linearization of (13) at $\bar{u}$
is

(17)
$\frac{\mathop{}\\!\mathrm{d}u}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0}u(t)+\mathcal{L}u(t),\quad
t\geq 0.$

The linearized operator $\mathcal{L}$ inherits the finite-rank property and
can be written as

$\mathcal{L}x=w\sum_{i=1}^{d}(\mathcal{L}_{i}x)\xi_{i},\quad\text{for all
}x\in D(\mathcal{A}_{0}),$

for $\mathcal{L}_{i}\colon D(\mathcal{A}_{0})\to\mathbb{C}$, $i=1,\dots,d$. We
define $\mathcal{L}^{\xi}\colon D(\mathcal{A}_{0})\to\mathbb{C}^{d}$ by
$\mathcal{L}^{\xi}:=(\mathcal{L}_{1};\dots;\mathcal{L}_{d}).$

Inspired by the PLS for iDEs we assume the following throughout. The local
stability of an equilibrium $\bar{u}$ of the semilinear equation (13) can be
inferred from the stability of the zero solution of the linearized equation
(17), and the position of the eigenvalues of $\mathcal{A}_{0}+\mathcal{L}$ in
$\mathbb{C}_{\rho}$ determines the stability of the zero solution of (17).

######  Proposition 4.1.

Let $\rho=\rho_{1}+\epsilon$ with $\rho_{1},\epsilon>0$. Then,
$\lambda\in\mathbb{C}_{\rho_{1}}$ is an eigenvalue of
$\mathcal{A}_{0}+\mathcal{L}$ if and only if $\lambda$ is a root of the
characteristic equation

(18) $\det(I_{\bar{d}+d}-A(\lambda))=0,$

where $A(\lambda):\mathbb{C}^{\bar{d}+d}\to\mathbb{C}^{\bar{d}+d}$ is the
linear operator

$A(\lambda)(\beta;\gamma)=\left(\mathcal{B}x(\cdot;\lambda,\beta,\gamma);\
\mathcal{L}^{\xi}x(\cdot;\lambda,\beta,\gamma)\right),$

and $x(\cdot;\lambda,\beta,\gamma)$ is the solution of the initial value
problem

(19)
$\begin{cases}\displaystyle{\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}\theta}}=(\lambda+\rho)x+w\sum_{i=1}^{d}\gamma_{i}\xi_{i},\\\\[5.0pt]
x(0)=\beta.\end{cases}$

with $\beta\in\mathbb{C}^{\bar{d}}$ and $\gamma\in\mathbb{C}^{d}$. In this
case, $x(\cdot;\lambda,\beta,\gamma)\in D(\mathcal{A}_{0})$ is the
eigenfunction corresponding to $\lambda$.

######  Proof 4.2.

Let $x=x(\cdot;\lambda,\beta,\gamma)$ be the solution of (19). By variation of
constants,

(20)
$x(\cdot;\lambda,\beta,\gamma)=e(\lambda+\rho)\beta+h(\lambda+\rho)\gamma,$

where $e(\lambda+\rho)\colon\mathbb{C}^{\bar{d}}\to X$ and
$h(\lambda+\rho)\colon\mathbb{C}^{d}\to X$ are defined by

$\displaystyle e(\lambda+\rho)\beta(\theta)$
$\displaystyle:=\mathrm{e}^{(\lambda+\rho)\theta}\beta=w(\theta)\mathrm{e}^{\lambda\theta}\beta,$
$\displaystyle h(\lambda+\rho)\gamma(\theta)$
$\displaystyle:=\textstyle{\left(\sum_{i=1}^{d}\gamma_{i}\xi_{i}\right)\int_{0}^{\theta}\mathrm{e}^{(\lambda+\rho)(\theta-s)}w(s)\,\mathop{}\\!\mathrm{d}s},\quad\theta\in\mathbb{R}_{-}.$

Now, (20) is an eigenfunction of $\mathcal{A}_{0}+\mathcal{L}$ corresponding
to $\lambda$ if and only if there exists $(\beta;\gamma)\neq 0$ such that
$\beta=\mathcal{B}x$ and $\gamma=\mathcal{L}^{\xi}x,$ i.e., if and only if
$\lambda$ is a root of (18).

The computation of (part of) the eigenvalues of $\mathcal{A}_{0}+\mathcal{L}$
in $\mathbb{C}_{\rho}$ is therefore turned into the computation of the roots
of (18). From Section 4, the position of the characteristic roots of (18) in
$\mathbb{C}_{\rho}$ determines the stability of the zero solution of (17).

Regarding the approximating system (16), the PLS holds from the classical
theory of ODEs. To approximate the stability of an equilibrium $\bar{u}$ of
(13), we first need the following results to relate the linearized equations.

######  Theorem 4.3 (One-to-one correspondence of equilibria).

Let $N\geq 1$. If $\overline{u}\in D(\mathcal{A}_{0})$ is an equilibrium of
(13), then $\overline{U}=\mathcal{R}_{N}\overline{u}$ is an equilibrium of
(16). Vice versa, if $\overline{U}$ is an equilibrium of (16), then
$\overline{u}=\mathcal{P}_{N,0}\overline{U}$ is an equilibrium of (13).

Moreover, linearization at equilibrium and PSD commute. In other words, the
linearization of (16) at $\overline{U}$ coincides with the PSD of (17) and
reads

(21)
$\frac{\mathop{}\\!\mathrm{d}U}{\mathop{}\\!\mathrm{d}t}(t)=\mathcal{A}_{0,N}U(t)+\mathcal{L}_{N}U(t),\quad
t\geq 0,$

with $\mathcal{L}_{N}=\mathcal{R}_{N}\mathcal{L}\mathcal{P}_{N,0}$.

######  Proof 4.4.

First, if $\overline{u}\in D(\mathcal{A}_{0})$ is an equilibrium of (16), then
$\overline{u}(0)=\mathcal{B}\overline{u}$. From (15) and the equality
$\mathcal{A}_{0}\overline{u}+\mathcal{F}(\overline{u})=0$, integrating we have

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\theta}\left(\frac{\overline{u}}{w}\right)=-\frac{\mathcal{F}(\overline{u})}{w}=-\sum_{i=1}^{d}\mathcal{F}_{i}(\overline{u})\xi_{i}\in\mathbb{C}^{\bar{d}}\quad\Rightarrow\quad\overline{u}(\theta)=w(\theta)(\mathcal{B}\overline{u}-\theta\sum_{i=1}^{d}\mathcal{F}_{i}(\overline{u})\xi_{i}).$

Since $\overline{u}/w$ is a polynomial of degree one and
$\overline{u}(0)=\mathcal{B}\overline{u}$, defining
$\overline{U}:=\mathcal{R}_{N}\overline{u}$, we have
$\mathcal{P}_{N,0}\overline{U}=\overline{u}$. Therefore

$\mathcal{A}_{0,N}\overline{U}+\mathcal{F}_{N}(\overline{U})=\mathcal{R}_{N}(\mathcal{A}_{0}\overline{u}+\mathcal{F}(\overline{u}))=0,$

hence $\overline{U}$ is an equilibrium of (16). Vice versa, let $\overline{U}$
be an equilibrium of (16) and let
$\overline{u}:=\mathcal{P}_{N,0}\overline{U}\in D(\mathcal{A}_{0}).$ Note that
$(\mathcal{A}_{0}\overline{u})/w$ is a polynomial of degree at most $N-1$, and
the function
$((\mathcal{A}_{0}\overline{u})/w+\sum_{i=1}^{d}\mathcal{F}_{i}(\overline{u})\xi_{i})$
is zero on $\Theta_{N}$. Therefore, it must be zero as a polynomial, and the
polynomial $\overline{u}$ is in fact an equilibrium of (13).

Let us now consider the linearization at $\overline{u}$. Trivially, (21) is
the PSD of (17). On the other hand, for the linearity of $\mathcal{R}_{N}$ and
$\mathcal{P}_{N,0}$ we have
$D(\mathcal{R}_{N}\mathcal{F}\mathcal{P}_{N,0})(\overline{U})=\mathcal{R}_{N}D(\mathcal{F}\mathcal{P}_{N,0})(\overline{U})=(\mathcal{R}_{N}D\mathcal{F}(\mathcal{P}_{N,0}\overline{U})\mathcal{P}_{N,0})\overline{U}=\mathcal{R}_{N}\mathcal{L}\mathcal{P}_{N,0}\overline{U}.$

For (21) we derive a discrete characteristic equation following [7,
Proposition 5.2].

######  Proposition 4.5.

The eigenvalues of $\mathcal{A}_{0,N}+\mathcal{L}_{N}$ in $\mathbb{C}_{\rho}$
are the roots of the discrete characteristic equation

(22) $\det(I_{\bar{d}+d}-A_{N}(\lambda))=0,$

where $A_{N}(\lambda):\mathbb{C}^{\bar{d}+d}\to\mathbb{C}^{\bar{d}+d}$ is the
linear operator

$A_{N}(\lambda)(\beta;\gamma)=(\mathcal{B}x_{N,0}(\cdot;\lambda,\beta,\gamma);\mathcal{L}x_{N,0}(\cdot;\lambda,\beta,\gamma)),$

and $x_{N,0}(\cdot;\lambda,\beta,\gamma)\in w\Pi_{N}^{\bar{d}}$ is the
solution of the collocation equation

(23)
$\begin{cases}\displaystyle{\frac{\mathop{}\\!\mathrm{d}x_{N,0}}{\mathop{}\\!\mathrm{d}\theta}}(\theta_{i})=(\lambda+\rho)x_{N,0}(\theta_{i})+w(\theta_{i})\sum_{i=1}^{d}\gamma_{i}\xi_{i}\,\quad
i=1,\dots,N,\\\\[5.0pt] x_{N,0}(0)=\beta.\end{cases}$

###  4.1 Convergence of the approximation of characteristic roots

We now show that the roots of (22) converge to the roots of (18) as
$N\to\infty$. We first assume that the operator $\mathcal{L}$ can be computed
exactly and so its PSD is given by
$\mathcal{L}_{N}=\mathcal{R}_{N}\mathcal{L}\mathcal{P}_{N,0}$. We later
briefly discuss the impact of quadrature rules.

Here we take $\rho=\rho_{1}+\epsilon$ for some $\rho_{1},\epsilon>0$, and we
choose

(24) $\theta_{j}:=-\frac{t_{j,N}}{2\rho_{1}},\quad j=1,\dots,N,$

where $t_{j,N}$, $j=1,\dots,N$ are either the zeros or the extrema of the
standard Laguerre polynomials in $\mathbb{R}_{+}$ (see Appendix A for
details). We refer to (24) as the _scaled Laguerre zeros_ or _extrema_.

######  Theorem 4.6.

Let $\rho=\rho_{1}+\epsilon$ with $\rho_{1},\,\epsilon>0$, and $B$ be an open
ball in $\mathbb{C}_{\rho_{1}}$. For $\Theta_{N}$ chosen as the scaled
Laguerre zeros or extrema (24), the collocation equation (23) admits a unique
solution $x_{N,0}(\cdot;\lambda,\beta,\gamma)\in w\Pi_{N}^{\bar{d}}$ for all
$N\geq 1$, $\lambda\in B$, $\beta\in\mathbb{C}^{\bar{d}}$ and
$\gamma\in\mathbb{C}^{d}.$ Moreover, let $x(\cdot;\lambda,\beta,\gamma)$ be
the solution of (19) and let

$C(\lambda):=\frac{\lambda}{\lambda+2\rho_{1}},\quad\text{and }\quad
K^{p}_{\epsilon}=\begin{cases}1&\text{ for }p=\infty,\ \epsilon\geq 0,\\\
(\frac{p\epsilon}{2\rho_{1}})^{-1/p}&\text{ for }1\leq p<\infty,\
\epsilon>0.\end{cases}$

The error
$e_{N}=x_{N,0}(\cdot;\lambda,\beta,\gamma)-x(\cdot;\lambda,\beta,\gamma)$
satisfies

$\|e_{N}\|_{X}\leq{C}_{0}(B,\epsilon)D_{N}(B)\,|(\beta;\gamma)|\\\ $

where for the scaled Laguerre zeros

$C_{0}(B,\epsilon):=kK^{p}_{\epsilon}\max_{\lambda\in\overline{B}}\left\\{\frac{\max\\{|\lambda|,1\\}}{\Re\lambda+\rho_{1}}\right\\}\quad\text{
and }\quad D_{N}(B):=\max_{\lambda\in\overline{B}}|C(\lambda)|^{N},$

while for the scaled Laguerre extrema

$C_{0}(B,\epsilon):=kK^{p}_{\epsilon}\,\max_{\lambda\in\overline{B}}\left\\{2+\frac{|\lambda|}{2\rho_{1}(\Re\lambda+\rho_{1})}\right\\}\quad\text{
and }\quad
D_{N}(B):=\max_{\lambda\in\overline{B}}\Big{|}\frac{C(\lambda)^{N}}{1-C(\lambda)^{N+1}}\Big{|},$

where $k$ is a positive constant. Finally, since $|C(\lambda)|<1$,
$\|e_{N}\|_{X}\to 0$ as $N\to\infty$.

######  Proof 4.7.

First we observe that, for the solution of (19), we can write $x=w\varphi\in
L^{p}$ with $\varphi\in L^{p}_{\rho}$ solution of

$\begin{cases}\displaystyle{\frac{\mathop{}\\!\mathrm{d}\varphi}{\mathop{}\\!\mathrm{d}\theta}}=\lambda\varphi+\textstyle{\sum_{i=1}^{d}\gamma_{i}\xi_{i}}&\theta\leq
0,\\\\[5.0pt] \varphi(0)=\beta.\end{cases}$

Second, by defining

$\theta\mapsto
t:=-2\rho_{1}\theta,\quad\lambda\mapsto\mu:=-\frac{\lambda}{2\rho_{1}},\quad\text{and}\quad\sum_{i=1}^{d}\gamma_{i}\xi_{i}\mapsto
c:=-\frac{1}{2\rho_{1}}\sum_{i=1}^{d}\gamma_{i}\xi_{i},$

we transform the problem for $\varphi$ on $\mathbb{R}_{-}$ with
$\Re\lambda>-\rho_{1}$ into the problem

$\begin{cases}\displaystyle{\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}\theta}}=\mu
y+c,&t\geq 0,\\\\[5.0pt] y(0)=\beta,\end{cases}$

for $y$ on $\mathbb{R}_{+}$ with $\Re\mu<\frac{1}{2}.$ Under this mapping,
$\mathrm{e}^{\rho_{1}\theta}=\mathrm{e}^{-t/2},$ while
$w(\theta)=\mathrm{e}^{\rho\theta}=\mathrm{e}^{-(1/2+\delta)t}$ with
$\delta:=\frac{\epsilon}{2\rho_{1}}$, and the collocation nodes $\theta_{j}$
are associated to the zeros or extrema $t_{j}$ of the standard Laguerre
polynomials in $\mathbb{R}_{+}.$ Now the result follows by applying
Proposition A.1 and Proposition A.3 in Appendix A on the positive real line
$\mathbb{R}_{+}$ and mapping the bounds (42) and (47) back to
$\mathbb{R}_{-}.$

Finally, from the error bounds on the collocation equation we are ready to
conclude the error bounds on the characteristic roots in the following
theorem.

######  Theorem 4.8.

Let $\rho=\rho_{1}+\epsilon$, with $\rho_{1},\epsilon>0,$ and
$\lambda^{*}\in\mathbb{C}_{\rho_{1}}$ be a root of (18) with multiplicity $m$.
Let $B$ be an open ball in $\mathbb{C}_{\rho_{1}}$ of center $\lambda^{*}$
such that $\lambda^{*}$ is the unique root of (18) in $B$. For $\Theta_{N}$
chosen as the scaled Laguerre zeros or extrema (24), there exists
$N_{1}:=N_{1}(B)$ such that, for $N\geq N_{1}$, equation (22) has $m$ roots
$\lambda_{i,N}$, $i=1,\dots,m$ (taking into account multiplicities) and

(25)
$\max_{i=1,\dots,m}|\lambda^{*}-\lambda_{i,N}|=O\left(D_{N}(B)^{1/m}\right).$

Hence $\displaystyle{\max_{i=1,\dots,m}}|\lambda^{*}-\lambda_{i,N}|\to 0$ as
$N\to\infty$.

######  Proof 4.9.

The steps follow the lines of [7, Section 5.3.2], to which we refer to further
details. From Theorem 4.6, for all choices of the $p$-norm and for both
choices of the nodes, for all $\lambda\in B$ we can bound

$\displaystyle\|(A(\lambda)-A_{N}(\lambda))(\beta;\gamma)\|_{X}$
$\displaystyle\leq\|(\mathcal{B},\mathcal{L}^{\xi})\|\cdot\|e_{N}\|_{X}\leq\|(\mathcal{B},\mathcal{L}^{\xi})\|\cdot
C_{0}(B,\epsilon)D_{N}(B)|(\beta;\gamma)|,$

with $C_{0}(B,\epsilon)$, $D_{N}(B)$ obtained from Theorem 4.6. Therefore,
considering the operator norm in $\mathbb{C}^{\bar{d}+d}$, following the proof
of [7, Proposition 5.3] we obtain the norm bound

$\|(A(\lambda)-A_{N}(\lambda))\|\leq C_{1}(B)C_{0}(B,\epsilon)D_{N}(B)$

for
$C_{1}(B):=\max\left\\{\left\|(\det)^{\prime}(I_{\bar{d}+d}-A(z)+\Delta)\right\|\colon
z\in\overline{B},\ \Delta\in\mathbb{C}^{\bar{d}+d},\ \|\Delta\|\leq
1\right\\}.$ This allows us to conclude the proof using Rouché’s Theorem from
complex analysis [35, Section 7.7].

The approximations (16) and (21) hold when $\mathcal{F}$ and $\mathcal{L}$ can
be computed exactly. However, such operators often involve integrals, hence
need to be approximated by $\widetilde{\mathcal{F}}$ and
$\widetilde{\mathcal{L}}$ via suitable quadrature rules. In this case, (16),
(21) and (22) still hold with the operators defined by
$\mathcal{F}_{N}=\mathcal{R}_{N}\widetilde{\mathcal{F}}\mathcal{P}_{N,0}$,
$\mathcal{L}_{N}=\mathcal{R}_{N}\widetilde{\mathcal{L}}\mathcal{P}_{N,0}$, and
$A_{N}(\lambda)(\beta;\gamma)=(\mathcal{B}x_{N,0}(\cdot;\lambda,\beta,\gamma);\
\widetilde{\mathcal{L}}x_{N,0}(\cdot;\lambda,\beta,\gamma))$. As a
consequence, the convergence analysis should also account for error

$\widetilde{Q}(B):=\sup\left\\{\frac{\|(\mathcal{L}-\widetilde{\mathcal{L}})x(\cdot;\lambda,\beta,\gamma)\|_{X}}{|(\beta;\gamma)|}\quad\colon\quad\lambda\in\overline{B},\
(\beta;\gamma)\in\mathbb{C}^{\bar{d}+d}\setminus\\{0\\}\right\\}$

and (25) should be replaced by

(26)
$\max_{i=1,\dots,m}|\lambda^{*}-\lambda_{i,N}|=O\left(D_{N}(B)^{1/m}\right)+O\left(\widetilde{Q}(B)^{1/m}\right).$

In practice, it is often convenient to use quadrature formulas associated with
$\Theta_{N}$ or $\Theta_{0,N}$, so that $\mathcal{L}$ is approximated by
$\widetilde{\mathcal{L}}_{N}$ by Gauss–Laguerre (on the scaled Laguerre zeros)
or Gauss–Radau–Laguerre (on the scaled Laguerre extrema plus zero) quadrature
formulas respectively, with corresponding error $\widetilde{Q}_{N}(B)$. In
this case, the convergence $\lambda_{i,N}\to\lambda^{*}$ follows from
$\widetilde{Q}_{N}(B)\to 0$ as $N\to\infty$ and the convergence rate depends
on the regularity of the integrands in suitable Sobolev spaces (see e.g. [32,
Section 5.1] for Gaussian quadrature rules on $\mathbb{R}_{+}$).

##  5 Numerical results

We here present some numerical tests on linear and nonlinear iDEs to
illustrate the error bounds obtained in Theorem 4.8. For linear iDEs we
compute the error on the characteristic roots and show convergence with
exponential order for increasing values of $N$. For nonlinear iDEs, we perform
equilibrium bifurcation analyses using the package MatCont (version 7p1) [12]
for MATLAB to illustrate the flexibility of the approach. The numerical tests
were performed using MATLAB (version 2023a). The standard Laguerre collocation
nodes and quadrature weights are computed using the publicly available MATLAB
suites from [24] and [41]. For ease of implementation, in the case of iDDEs we
use one single entry to discretize the first element of a pair in the space
$X$, hence $(\alpha,\psi)\in X$ is discretized by a vector of dimension $N+1$
(instead of dimension $\bar{d}N=2N$). MATLAB codes to reproduce the tests are
available at https://github.com/scarabel/.

###  5.1 Linear iDEs

Consider the linear iDDE

(27)
$\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}(t)=ay(t)+\int_{0}^{+\infty}k(s)y(t-s)\mathop{}\\!\mathrm{d}s,\quad
t\geq 0,$

for $a\in\mathbb{R}$ and $k\in L^{1}(\mathbb{R}^{+})$ with
$\int_{0}^{\infty}|k(s)|\mathrm{e}^{\rho^{*}s}\mathop{}\\!\mathrm{d}s<\infty$
for some $\rho^{*}>0$. If $a+\int_{0}^{+\infty}k(s)\mathop{}\\!\mathrm{d}s\neq
0,$ then $\bar{y}=0$ is the unique equilibrium of (27). The characteristic
roots in $\mathbb{C}_{\rho}$ are the solutions of
$\lambda=a+\widehat{k}(\lambda)$, where $\widehat{k}$ denotes the Laplace
transform of $k$, and are eigenvalues of the linear operator
$\mathcal{A}^{\text{DDE}}:=\mathcal{A}_{0}^{\text{DDE}}+\mathcal{L}^{\text{DDE}}$,
see (7). We study how they are approximated by the eigenvalues of the PSD
$\mathcal{A}_{N}^{\text{DDE}}$. The parameter values for the numerical tests
are collected in Table 1.

Table 1: Parameter values for the numerical tests on the linear iDDE (27) and iRE (29). |  | kernel | $\mu$ | $a$ | $k_{0}$ | $\sigma$ |  | $\lambda$
---|---|---|---|---|---|---|---|---
(a1) | iDDE | exponential | 2 | 3 | $-6$ | - |  | $\lambda=0$
(a2) | iDDE | exponential | 2 | 3 | $-6$ | - |  | $\lambda=a-\mu$
(b) | iDDE | exponential | 2 | 2 | $-8$ | - |  | $\lambda=\mu\,i$
(c) | iDDE | exponential | 2 | 6 | $-16$ | - |  |  $\lambda=(a-\mu)/2$ (double)
(d) | iDDE | gamma | 4 | 0 | - | 2 |  | $(\lambda-1)(\lambda+\mu)^{2}=\mu^{2}$
(e) | iDDE | gamma | 4 | 0 | - | $\pi$ |  | $\lambda=a+\widehat{k}(\lambda)$
(f) | iRE | non-monotonic | 1 | 1 | 1 | - |  |  (30)
(g) | iRE | non-monotonic | 1.5 | 1 | 3 | - |  |  (30)

For an exponential kernel $k(s)=k_{0}\mathrm{e}^{-\mu s}$,
$\mu\in\mathbb{R}_{+}$, it is easy to show that the characteristic roots exist
in $\\{\Re\lambda>-\mu\\}$ and they can be written explicitly as

$\lambda=(a-\mu\pm\sqrt{(\mu+a)^{2}+4k_{0}})/2.$

Note in particular that: $\lambda=0$ is a root if $k_{0}=-a\mu$; the roots are
real if $k_{0}\geq-(\mu+a)^{2}/4$, and have double multiplicity if
$k_{0}=-(\mu+a)^{2}/4$. In the numerical tests we take the two choices
$\rho_{1}=\mu/2$ and $\rho_{1}=\mu/4$ to compare the impact of the nodes used
for quadrature, and we fix $\rho=\rho_{1}$. To study the convergence of
eigenvalues and eigenvectors, the eigenvector $\psi_{\lambda,N}$ of
$\mathcal{A}_{N}^{\text{DDE}}$ is first normalized such that the first
component equals one to obtain an approximation of
$\psi_{\lambda}=\mathrm{e}^{\lambda\cdot}$. Fig. 1 shows the log-log plot of
the computed errors $|\lambda-\lambda_{N}|$ and
$\|w(\psi_{\lambda}-\psi_{\lambda,N})\|_{\infty,N}$, where
$\|\cdot\|_{\infty,N}$ indicates the discrete supremum vector norm (i.e., the
maximum difference on the collocation nodes), as a function of $N$: for finite
(polynomial) order of convergence we would expect to observe straight
decreasing lines; an exponential (infinite) order of convergence translates
into exponentially decreasing trends [3, Section 2.4]. Fig. 1 illustrates
several features predicted by Theorem 4.8: the error decreases exponentially
as $N$ increases; both larger modulus and larger real part imply slower
convergence (a2, b); the maximal accuracy is halved for double roots (c); and
the order of convergence depends on the quadrature error, so quadrature rules
that are exact on the exponential kernel ($\rho_{1}=\mu/2$) guarantee faster
convergence than non-exact rules ($\rho_{1}=\mu/4$). In the case of
$\lambda=0$ (a1) the theoretical collocation error is zero if the quadrature
is exact, and indeed the errors in both, eigenvalue and eigenvector, fall to
machine precision already for $N=1$ when $\rho_{1}=\mu/2$, while the error
coincides with the quadrature error, and the convergence is exponential, for
$\rho_{1}=\mu/4$. Moreover, for $\lambda=0$, Gauss-Laguerre quadrature rules
and Gauss–Radau–Laguerre rules show a similar trend, while in the other tests
the former seem to guarantee a faster exponential decay compared to the latter
in general.

In tests (d) and (e) we considered a Gamma-distributed kernel
$k(s)=f_{\mu}^{(\sigma)}(s)$, with rate parameter $\mu$ and shape parameter
$\sigma$, i.e.,

(28) $f_{\mu}^{(\sigma)}(s)=\frac{\mu^{\sigma}s^{\sigma-1}\mathrm{e}^{-\mu
s}}{\Gamma(\sigma)},\quad s\geq 0.$

The reference eigenvalues in this case were obtained by solving algebraically
the characteristic equation using the MATLAB solver `fsolve`. In this case,
too, the choice of a suitable $\rho_{1}$ for quadrature rules can be linked to
$\mu$, with $\rho_{1}=\mu/2$ giving faster convergence than $\rho_{1}=\mu/4$.
The case of a non-integer shape parameter in Fig. 1(e) is particularly
interesting, as it illustrates the role of the quadrature error in limiting
the convergence order, as shown by Eq. 26: the errors decay exponentially when
the integrals are computed with the built-in MATLAB integral function (black
lines, tolerance set at $10^{-12}$), while the convergence is of polynomial
order (visible in the numerical errors decaying asymptotically as straight
lines in the log-log diagram) when the integrals are approximated with the
Gaussian quadrature rules associated to the collocation nodes (gray lines).
The polynomial order of convergence in the latter case was estimated as
$N^{\alpha}$ with $\alpha\approx-3.1479$.

Figure 1:  Linear DDE (27): log-log plot of $|\lambda-\lambda_{N}|$ (solid)
and $\|w(\psi_{\lambda}-\psi_{\lambda,N})\|_{\infty,N}$ (dotted), for
parameters specified in Table 1. See main text for further details.

Next, we consider the linear iRE

(29) $y(t)=\int_{0}^{+\infty}k(s)y(t-s)\mathop{}\\!\mathrm{d}s,\quad t>0,$

where $k\in L^{\infty}(\mathbb{R}^{+})$ with ${\rm
ess}\sup_{s\in\mathbb{R}_{+}}|k(s)|\mathrm{e}^{\rho^{*}s}<\infty$ for some
$\rho^{*}>0$. If $\int_{0}^{\infty}k(s)\mathop{}\\!\mathrm{d}s\neq 0$,
$\bar{y}=0$ is the unique equilibrium of (29). The characteristic roots in
$\mathbb{C}_{\rho}$ are the solutions of $1=\widehat{k}(\lambda)$, and are
eigenvalues of the linear operator
$\mathcal{A}^{\text{RE}}:=\mathcal{A}_{0}^{\text{RE}}+\mathcal{L}^{\text{RE}}$,
see (12), with corresponding PSD $\mathcal{A}_{N}^{\text{RE}}$.

Fig. 2 shows the errors in the approximation of $\lambda$ and the
corresponding eigenfunction
$\psi_{\lambda}=j^{-1}(w\varphi_{\lambda})=\mathrm{e}^{\lambda\cdot}$, in the
discrete weighted $\infty$\- and $1$-norm. To obtain a normalization
comparable with the exact eigenfunction, we first applied the differentiation
matrix to the discrete eigenfunction $\varphi_{\lambda,N}$ of
$\mathcal{A}_{N}^{\text{RE}}$ and then normalized it so that the first
component equals one, to obtain an approximation $\psi_{\lambda,N}$ of
$\psi_{\lambda}$. For an exponential kernel, the convergence plots show very
similar trends as for the iDDE, so we did not include them in the paper.
Instead, we considered the non-monotone kernel $k(s)=k_{0}\mathrm{e}^{-\mu
s}(\sin(as)+1)$, for $a,\mu,k_{0}\in\mathbb{R}_{+}$, for which the
characteristic roots (with $\Re\lambda>-\mu$) satisfy

(30) $\lambda^{3}+\lambda^{2}(3\mu-
k_{0})+\lambda(3\mu^{2}-2k_{0}\mu+a^{2}-k_{0}a)+(\mu^{2}+a^{2})(\mu-
k_{0})-k_{0}a\mu=0.$

Parameter values are collected in Table 1. In the numerical tests we
constructed the collocation nodes and quadrature weights for two different
choices, $\rho_{1}=\mu/2$ and $\rho_{1}=\mu/3$, and we fixed
$\rho=\mu>\rho_{1}$. The weighted 1-norm of the error is then computed with
respect to $w(\theta)=\mathrm{e}^{\rho\theta}$. The convergence trends in Fig.
2 highlight similar behavior as for iDDEs, consistent with Theorem 4.8.

Figure 2:  Linear RE (29): log-log plot of $|\lambda-\lambda_{N}|$ (solid) and
$\|\psi_{\lambda}-\psi_{\lambda,N}\|_{1,N}$ (dotted), for parameter values
specified in Table 1.

###  5.2 Numerical bifurcation analysis of nonlinear iDEs

We now consider nonlinear iDEs and perform the numerical bifurcation analysis
of equilibria using the numerical continuation package MatCont applied to the
ODE system (16). For each example below we specify the parameters used to
construct the ODE system and perform the continuation. The bifurcation
diagrams show the output of the numerical continuation. Although Laguerre
zeros appear to ensure better accuracy than Laguerre extrema in the linear
examples (Fig. 1 and 2), numerical instabilities caused by the eigenvalues
with large modulus seem to substantially affect the computations in the
nonlinear case. For this reason we use Laguerre extrema and the corresponding
Gauss–Radau–Laguerre quadrature rules.

Consider the nonlinear iDDE taken from [2]

(31)
$\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}(t)=-\delta_{A}y(t)+b\int_{0}^{+\infty}f_{m/\tau}^{(m)}(s)\mathrm{e}^{-\delta_{J}s-ay(t-s)}y(t-s)\mathop{}\\!\mathrm{d}s,\quad
t\geq 0,$

where $\tau,m,a,b,\delta_{A},\delta_{J}>0$ and $f_{\mu}^{(m)}(s)$ is the Gamma
distribution (28). When $m\in\mathbb{N}$, (31) admits an equivalent ODE
representation via the linear chain trick [27], which we use as reference to
compare the bifurcation diagrams computed via the PSD. If $m$ is large enough,
the positive equilibrium undergoes two Hopf bifurcations when varying $\tau$.
The Hopf bifurcation curve in $(\tau,m)$, which separates the parts of the
plane where the positive equilibrium is either stable or unstable, is shown in
Fig. 3 (top-left panel), showing also that the curves computed for $N=10$ and
$N=20$ are indistinguishable in the figure. Note that, since $y$ appears as an
exponent in (31), to obtain reliable results we used the positive part of the
state vector when computing the right-hand side of the PSD. Fig. 3 shows the
convergence of the observed error in the approximation of the two Hopf
bifurcation points and the corresponding eigenvalues for $m=6.5,\ 7,\ 7.5$.
Note that the discretization approach allows to consider distributions with
non integer shape parameter $m$, offering a very general approach for iDDEs
and iREs.

Figure 3:  Top-left: Hopf bifurcation curve of (31) in $(\tau,m)$ computed
with MatCont, for $\delta_{A}=0.5,\,\delta_{J}=1,\,a=7,\,b=350$, and
$\rho=(\delta_{J}+m/\tau)/4$, for $N=5,10,20$ (stable positive equilibrium
below the curve). Note that the curves for $N=10,\,20$ are indistinguishable.
The other panels show the error in the two Hopf bifurcation points and the
critical characteristic root at Hopf, for $m=6.5,\,7,\,7.5$, with MatCont
tolerance $10^{-10}$. Reference Hopf values are computed from the equivalent
ODE formulation [27].  Figure 4:  Left: Hopf bifurcation curve of (32) in
$(\mu,\gamma/\mu)$, for $\eta:=\beta_{0}\mathrm{e}^{-\mu}/\mu$, computed using
MatCont tolerance and $10^{-6}$, $\rho=\mu/2$. The positive equilibrium exists
above the line $\eta=1$ (‘branching point’, BP). Gray lines are computed with
DDE-BIFTOOL on the equivalent DDE formulation. Right: error in BP and Hopf for
$\mu=2$, with MatCont tolerance $10^{-10}$. Reference values are computed with
$N=50$.

Consider now a nonlinear iRE equivalent to the Nicholson’s blowflies DDE [26]

(32) $y(t)=\beta_{0}\,\mathrm{e}^{-\int_{1}^{\infty}y(t-s)\mathrm{e}^{-\mu
s}\mathop{}\\!\mathrm{d}s}\,\int_{1}^{\infty}y(t-s)\mathrm{e}^{-\mu
s}\mathop{}\\!\mathrm{d}s,\quad t>0,$

for $\beta_{0},\mu\geq 0$, where $y(t)$ represents the total birth rate of a
population where individuals become adult after a maturation delay which is
normalized to 1. Equation (32) admits a unique nontrivial equilibrium
$\overline{y}=(\log\frac{\beta_{0}}{\mu}-\mu)\mu\mathrm{e}^{\mu}$ for
$\beta_{0}>\mu\mathrm{e}^{\mu}$, which undergoes a Hopf bifurcation as
$\beta_{0}$ increases, and, if $\mu$ is large enough, the periodic solution
undergoes a period doubling bifurcation [11]. The stability regions of the
positive equilibrium in a two-parameter plane are shown in Fig. 4 (left),
compared to those obtained by applying DDE-BIFTOOL to the equivalent DDE with
finite delay (for the equivalence see [11]). The numerical error in the
approximation of the transcritical and Hopf bifurcations with respect to
$\beta_{0}$ are illustrated in Fig. 4 (right), for $\mu=4$. For implementation
purposes, the right-hand side of (32) was reformulated integrating by parts,
and then computed with respect to the integrated state
$Y(t):=\int_{0}^{t}y(s)\mathop{}\\!\mathrm{d}s$. The quadrature rules on the
interval $[1,+\infty)$ were obtained by appropriately shifting the standard
Gauss–Radau–Laguerre quadrature nodes. We note that the eigenvalue at the
branching point ($\lambda=0$) is approximated to machine precision already
with $N=1$, consistently with the analysis for linear equations. The
approximation of the bifurcation point loses some accuracy probably due to the
nonlinear solver used by MatCont for evaluating the test functions. The
eigenvalues at the Hopf bifurcation point and the bifurcation point itself
show a slower convergence, consistently with the fact that, at the Hopf
bifurcation point, the eigenvalues are imaginary and therefore have strictly
positive modulus.

##  6 Conclusions

In this paper we have introduced a general abstract framework for the
pseudospectral reduction of nonlinear iDEs to a finite-dimensional system of
ODEs, which can be studied with well-established software for ODEs
(specifically we used MatCont for MATLAB). For iREs we have proposed a novel
approach compared to [27], exploiting the equivalence via integration proposed
for finite delay in [37]. In this way, the PSD returns directly an ODE system
without involving an implicit algebraic condition.

For linear(ized) iDDEs and iREs we have proved the convergence of the PSD
approximation of the characteristic roots to the true ones, answering an
important question raised when the method was first proposed in [27]. The
convergence holds for a general weighted $L^{p}$-norm, therefore also in the
Hilbert space $L^{2}_{\rho}$. The numerical tests illustrate the theoretical
result of convergence of characteristic roots and eigenfunctions in norm
$\|\cdot\|_{\infty,\rho}$ for iDDEs and $\|\cdot\|_{1,\rho}$ for iREs. The
theoretical error bound depends on the modulus, the real part, and the
multiplicity of the target eigenvalue, as well as the discretization index
$N$, the choice of the family of collocation nodes (chosen as the zeros or
extrema of the scaled Laguerre polynomials, orthogonal with respect to a
weight $\mathrm{e}^{\rho_{1}\cdot}$, $\rho_{1}<\rho$) and the quadrature
formula used to approximate the integrals. In particular, approximating
eigenvalues with larger modulus with a given accuracy requires larger
discretization indices. Despite this, since the right half-plane
$\mathbb{C}_{\rho}$ contains a finite number of eigenvalues, our theoretical
result guarantees that all eigenvalues in $\mathbb{C}_{\rho}$ can be
approximated with desired accuracy by choosing $N$ large enough. However,
determining such value of $N$ is a difficult problem, unless one can estimate
_a priori_ the location of the spectrum. To investigate this problem
experimentally, we performed extensive numerical tests, which confirm that the
observed accuracy strongly depends on the specific eigenvalue of interest and
the chosen quadrature formula. For optimal choices, discretization indices of
the order of 20–30 usually guarantee machine precision for eigenvalues with
multiplicity one.

Numerical tests on nonlinear equations demonstrate the flexibility of the
method to study the bifurcation of equilibria. To verify the correctness of
the computed bifurcation diagrams, in this paper we have considered equations
in which the bifurcation diagram can be obtained via an alternative route,
hence allowing to compare the PSD output with reference results: the nonlinear
iDDE (31) admits an equivalent ODE formulation, at least for integer shape
parameters of the kernel, whose bifurcation properties can be studied
independently with MatCont; and the iRE (32) admits an equivalent formulation
as a DDE with finite delay that can be studied with bifurcation software for
DDEs (here we used DDE-BIFTOOL for MATLAB).

Future work includes the proof of convergence of the resolvent operators, that
may be useful in turn to prove the convergence of the solution operators of
the initial value problems and of the periodic solutions. From a computational
point of view, we plan to investigate numerically the impact of different sets
of collocation nodes and use truncated interpolation and quadrature schemes to
improve computation times in the spirit of [33]. While the abstract
formulation encompasses other types of models, the convergence proof in this
paper relies on several assumptions specific to iDEs, including the form of
the boundary condition defined by $\mathcal{B}$ and that the operator
$\mathcal{F}$ maps on a $d$-dimensional space of functions $w\Xi_{d}$, on
which weighted interpolation is exact. One could consider extending the
techniques to more general boundary conditions in the form $\mathcal{B}x=0$.
For more general basis functions of $\Xi_{d}$ or different (e.g., truncated)
interpolation schemes, the analogous of Theorem 4.6 would require the
assumption that the interpolation scheme is convergent on the functions in
$\Xi_{d}$. Similarly, the one-to-one correspondence of equilibria and
linearization (Theorem 4.3) may hold only in an approximate sense, if
interpolation is not exact on the equilibria. Finally, the analysis conducted
in the appendix also provides the exact location of the spectrum of the
reduced differentiation matrix, i.e., the differentiation matrix with zero
boundary condition, associated with the zeros and extrema of the standard
Laguerre polynomials on the positive half-line. Understanding the behavior of
the spectrum of differentiation matrices is a notoriously difficult problem,
and these results add a first significative contribution to the study of such
matrices, which are commonly used in the numerical solution of differential
equations. We aim to perform further explorations in future work.

##  Appendix A General results for a collocation problem on $\mathbb{R}_{+}$

In this section we assume $d=1$ for simplicity and, in order to work with the
family of the standard Laguerre polynomials $\\{L_{N}\\}_{N\geq 0}$,
orthogonal with respect to the weight $\mathrm{e}^{-t}$ on the positive half-
line $\mathbb{R}_{+}:=[0,+\infty)$, we consider the equivalent reformulation
of the collocation problem on $\mathbb{R}_{+},$ and we study the existence and
uniqueness and the convergence of the collocation polynomials for two relevant
choices of the collocation points: the zeros and the extrema of $L_{N}$, for
which we refer for instance to [32, Chapter 2.3.5]. The same approach can be
easily extended to other collocation points.

Given $\beta,c\in\mathbb{C}$ and $\mu\in\mathbb{C}$ with $\Re\mu<1/2$, let
$y(t)=\mathrm{e}^{\mu
t}\beta+\int_{0}^{t}\mathrm{e}^{\mu(t-s)}c\mathop{}\\!\mathrm{d}s$ be the
solution of the initial value problem

(33)
$\begin{cases}\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}=\mu\,y+c,\quad
t\in\mathbb{R}_{+},\\\\[5.0pt] y(0)=\beta.\end{cases}$

Note that, since $\Re\mu<1/2$, for all $\delta\geq 0$ and
$w_{\delta}(t):=\mathrm{e}^{-(1/2+\delta)t},\ t\in\mathbb{R}_{+},$ we have
$w_{\delta}y\in L^{p}(\mathbb{R}_{+},\mathbb{C})$ for $1\leq p\leq+\infty.$

Let $\Pi_{N}$ denote the set of algebraic polynomials of degree at most $N$ on
$\mathbb{R}_{+}$. Given a set of distinct points $\\{t_{i,N}\\}_{i=1}^{N}$ on
$\mathbb{R}_{+},$ the collocation polynomial $p_{N}\in\Pi_{N}$ for (33) is
defined by

(34)
$\begin{cases}\frac{\mathop{}\\!\mathrm{d}p_{N}}{\mathop{}\\!\mathrm{d}t}=\mu\,p_{N}+c,\quad
t\in\\{t_{i,N}\\}_{i=1}^{N},\\\\[5.0pt] p_{N}(0)=\beta.\end{cases}$

Since
$\frac{\mathop{}\\!\mathrm{d}p_{N}}{\mathop{}\\!\mathrm{d}t}-\mu\,p_{N}-c$ is
a polynomial of degree $N$ which vanishes at the collocation points, the
collocation equation (34) admits the following equivalent representation

(35)
$\left\\{\begin{array}[]{l}\frac{\mathop{}\\!\mathrm{d}p_{N}}{\mathop{}\\!\mathrm{d}t}-\mu\,p_{N}=c+k_{N}q_{N},\quad
t\in\mathbb{R}_{+},\\\\[5.0pt] p_{N}(0)=\beta,\end{array}\right.$

where $q_{N}(t)=\prod_{i=1}^{N}(t-t_{i,N})$ and $k_{N}$ is a suitable
constant. By expressing the polynomials in (35) as follows

$\frac{\mathop{}\\!\mathrm{d}p_{N}}{\mathop{}\\!\mathrm{d}t}(t)=\sum\limits_{j=0}^{N-1}a_{j}t^{j},\quad
p_{N}(t)=\beta+\sum\limits_{j=1}^{N}\frac{a_{j-1}}{j}\,t^{j},\quad\text{and
}\quad q_{N}(t)=\sum_{j=0}^{N}\frac{b_{j}}{j!}t^{j},$

from (35) we derive the following recurrence relations

(36) $\left\\{\begin{array}[]{l}a_{0}=\mu\beta+c+k_{N}b_{0}\\\
a_{j}-\frac{\mu}{j}a_{j-1}=k_{N}\frac{b_{j}}{j!},\quad j=1,\ldots,N,\\\
a_{N}=0.\end{array}\right.$

Now, by easy computations, we get from (36)

$a_{j}=\frac{\mu^{j}}{j!}(\mu\beta+c)+k_{N}d_{j},\quad j=1,\ldots,N,$

with $d_{j}$ given by

(37) $d_{j}j!=\sum\limits_{k=0}^{j}\mu^{j-k}\,b_{k},\qquad j=0,1,\dots,N.$

Finally, by requiring that $a_{N}=0,$ we get the following equation for
$k_{N}$

(38) $\mu^{N}(\mu\beta+c)+k_{N}d_{N}N!=0.$

Note that each $d_{j}$ in (37) depends on $\mu$ and, through the polynomial
$q_{N},$ on the collocation points. In the following we focus on two relevant
choices: the zeros and the extrema of the standard Laguerre polynomial
$L_{N}$. For both, we study the existence and uniqueness, and the convergence
w.r.t. the weight function $w_{\delta}$ in $L^{p}(\mathbb{R}_{+},\mathbb{C})$,
namely the behavior of $\|w_{\delta}(p_{N}-y)\|_{p}$ as $N\to\infty.$

Let $\hat{e}_{N}:=w_{\delta}(p_{N}-y).$ From (33) and (35) we easily get

(39)
$\hat{e}_{N}(t)=k_{N}w_{\delta}(t)\int_{0}^{t}\mathrm{e}^{\mu(t-s)}q_{N}(s)\mathrm{d}s,\quad
t\in\mathbb{R}_{+}.$

######  Proposition A.1 (Existence, uniqueness, convergence for Laguerre
zeros).

Let $\mu\in\mathbb{C}$, $\Re\mu<1/2,$ and let $\\{t_{i,N}\\}_{i=1}^{N}$ be the
zeros of the $N$-degree standard Laguerre polynomial $L_{N}.$ Then there
exists a unique polynomial solution $p_{N}$ of (34) for all $N\geq 1.$
Moreover, defining

(40) $C(\mu):=\frac{\mu}{\mu-1},$

we have

(41) $k_{N}q_{N}(t)=-C(\mu)^{N}L_{N}(t)\,(\mu\beta+c),\quad
t\in\mathbb{R}_{+},$

and the following error bound

(42)
$\|\hat{e}_{N}\|_{p}\leq|C(\mu)|^{N}\frac{(|\mu||\beta|+|c|)K_{\delta}^{p}}{(1/2-\Re\mu)},$

where

(43) $K_{\delta}^{p}=\begin{cases}1,&p=\infty,\,\delta\geq 0,\\\
(p\delta)^{-1/p},&1\leq p<\infty,\,\delta>0.\end{cases}$

The convergence $\|\hat{e}_{N}\|_{p}\to 0$ as $N\to\infty$ follows for
$p=\infty$ with $\delta\geq 0,$ and for $\ 1\leq p<+\infty$ with $\delta>0.$

######  Proof A.2.

For the standard Laguerre polynomial $L_{N}$ it holds

$q_{N}(t)=(-1)^{N}N!\,L_{N}(t)=(-1)^{N}N!\,\sum\limits_{j=0}^{N}(-1)^{j}{N\choose
j}\,\frac{t^{j}}{j!},$

see [32, Chapter 2.3.5]. Hence, $b_{j}=(-1)^{N}N!(-1)^{j}{N\choose j}$, and,
from (37),

(44) $d_{N}=(-1)^{N}\sum\limits_{j=0}^{N}\mu^{N-j}(-1)^{j}{N\choose
j}=(1-\mu)^{N}\neq 0$

since $|1-\mu|\geq|1-\Re\mu|>\frac{1}{2}.$ Then the existence and uniqueness
of the collocation solution follows for all $N\geq 1.$ From (38) and (40) we
obtain $k_{N}=-\frac{C(\mu)^{N}}{(-1)^{N}N!}(\mu\beta+c)$ and (41).

Let us study now the error $\hat{e}_{N}$. From (39) and (41) we get

$\hat{e}_{N}(t)=C(\mu)^{N}\,\mathrm{e}^{-\delta
t}\int_{0}^{t}\mathrm{e}^{(\mu-1/2)(t-s)}\mathrm{e}^{-s/2}L_{N}(s)\,\mathrm{d}s\,(\mu\beta+c),\quad
t\in\mathbb{R}_{+}.$

Since $\mathrm{e}^{-s/2}|L_{N}(s)|\leq L_{N}(0)=1$ (see [32, Section 2.3.5]),
we get

$\displaystyle|\hat{e}_{N}(t)|$
$\displaystyle\leq\,|C(\mu)|^{N}\mathrm{e}^{-\delta t}\
\frac{1-\mathrm{e}^{(\Re\mu-1/2)t}}{1/2-\Re\mu}\,(|\mu||\beta|+|c|),\quad
t\in\mathbb{R}_{+}.$

The converge follows since $|C(\mu)|<1$ for $\Re\mu<1/2$ and from the bounds

(45) $\begin{array}[]{ll}\|\mathrm{e}^{-\delta
t}(1-\mathrm{e}^{(\Re\mu-1/2)t})\|_{\infty}\leq 1,&\text{ for all }\delta\geq
0,\\\\[5.0pt] \|\mathrm{e}^{-\delta
t}(1-\mathrm{e}^{(\Re\mu-1/2)t})\|_{p}\leq(p\delta)^{-1/p},&\text{ for all
}\delta>0\text{ and }1\leq p<\infty.\end{array}$

######  Proposition A.3 (Existence, uniqueness, convergence for Laguerre
extrema).

Let $\mu\in\mathbb{C}$, $\Re\mu<1/2,$ and let $\\{t_{i,N}\\}_{i=1}^{N}$ be the
extrema of the $(N+1)$-degree standard Laguerre polynomial, i.e., the zeros of
$\frac{\mathrm{d}}{\mathrm{d}t}L_{N+1}$. Then there exists a unique polynomial
solution $p_{N}$ of (34) for all $N\geq 1.$ Moreover, for $C(\mu)$ defined as
in (40) and

$D_{N}(\mu):=\frac{C(\mu)^{N}}{1-C(\mu)^{N+1}},$

we have that

(46)
$k_{N}q_{N}(t)=D_{N}(\mu)\,\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}L_{N+1}(t)\left(\frac{\mu\beta+c}{\mu-1}\right),\
t\in\mathbb{R}_{+},$

and the following bound

(47)
$\|\hat{e}_{N}\|_{p}\leq|D_{N}(\mu)|\frac{(|\mu||\beta|+|c|)K_{\delta}^{p}}{|\mu-1|}\left(2+\frac{|\mu|}{1/2-\Re\mu}\right)$

with $K_{\delta}^{p}$ defined in (43). The convergence $\|\hat{e}_{N}\|_{p}\to
0$ as $N\to\infty$ follows for $p=\infty$ with $\delta\geq 0,$ and for $\
1\leq p<+\infty$ with $\delta>0.$

######  Proof A.4.

For the polynomial $\frac{\mathrm{d}}{\mathrm{d}t}L_{N+1}$ it holds

$q_{N}(t)=(-1)^{N}N!\,L_{N}^{(1)}(t)=(-1)^{N}N!\sum\limits_{j=0}^{N}\frac{(-1)^{j}}{j!}{N+1\choose
N-j}\,t^{j}$

(see [32, Chapter 2.3.5]). Hence $b_{j}=(-1)^{N}N!(-1)^{j}{N+1\choose j+1}$,
from (37) we get

(48) $\displaystyle
d_{N}=(-1)^{N+1}\mu^{N+1}\sum_{j=0}^{N}(-\mu)^{-(j+1)}{N+1\choose
j+1}=(-1)^{N+1}\left((\mu-1)^{N+1}-\mu^{N+1}\right)\neq 0$

since $\Re\mu<1/2.$ Then the existence and uniqueness of the collocation
solution follows for all $N\geq 1.$ From (38) and from

$\frac{\mu^{N}}{(\mu-1)^{N+1}-\mu^{N+1}}=\frac{C(\mu)^{N}}{(\mu-1)-\mu
C(\mu)^{N}}=\frac{C(\mu)^{N}}{(\mu-1)(1-C(\mu)^{N+1})},$

we obtain

$k_{N}=\frac{\mu^{N}(\mu\beta+c)}{d_{N}N!}=\frac{(-1)^{N+1}}{N!}\left[\frac{C(\mu)^{N+1}}{1-C(\mu)^{N+1}}\beta+\frac{C(\mu)^{N}}{1-C(\mu)^{N+1}}\frac{c}{\mu-1}\right]$

and (46). Let us now consider the error $\hat{e}_{N}$. From (39) and (46) we
get

$\displaystyle\hat{e}_{N}(t)$
$\displaystyle=D_{N}(\mu)\,w_{\delta}(t)\int_{0}^{t}\mathrm{e}^{\mu(t-s)}\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}s}L_{N+1}(s)\mathop{}\\!\mathrm{d}s\,\frac{\mu\beta+c}{\mu-1}$
$\displaystyle=D_{N}(\mu)\,w_{\delta}(t)\Big{(}L_{N+1}(t)-L_{N+1}(0)\mathrm{e}^{\mu
t}+\mu\int_{0}^{t}\mathrm{e}^{\mu(t-s)}L_{N+1}(s)\mathop{}\\!\mathrm{d}s\Big{)}\,\frac{\mu\beta+c}{\mu-1}$

Since $\mathrm{e}^{-s/2}|L_{N+1}(s)|\leq L_{N+1}(0)=1$ (see [32, Chapter
2.3.5]), we get

$|\hat{e}_{N}(t)|\leq|D_{N}(\mu)|\mathrm{e}^{-\delta
t}\left(2+|\mu|\frac{1-\mathrm{e}^{(\Re\mu-1/2)t}}{1/2-\Re\mu}\right)\left(\frac{|\beta||\mu|+|c|}{|\mu-1|}\right)$

and (47) from (45). Since $|D_{N}(\mu)|\to 0$ as $N\to\infty$, the convergence
follows.

As a final important remark, note that from (38) with $\beta=0$ and $c=0$ we
get the equation $d_{N}=0$, whose roots $\mu\in\mathbb{C}$ are the eigenvalues
of the reduced differentiation matrix (with zero boundary condition). In
particular, for Laguerre zeros, the differentiation matrix has only the
eigenvalue $\mu=1$ with multiplicity $N$ (see (44)), while for Laguerre
extrema, the $N$ eigenvalues are aligned along $\Re\mu=\frac{1}{2}$ (see
(48)). This result adds a new insight into the spectral properties of such
matrices, which are very useful in the numerical solution of differential
equations.

## Acknowledgments

The authors are members of CDLab (Computational Dynamics Laboratory,
Department of Mathematics and Computer Science, University of Udine), of INdAM
Research group GNCS, and of UMI Research group “Modellistica socio-
epidemiologica”. F.S. acknowledges the QJMAM Fund for Applied Mathematics for
travel funding that facilitated the completion of this paper.

## References

  * [1] A. Andò, S. De Reggi, D. Liessi, and F. Scarabel, A pseudospectral method for investigating the stability of linear population models with two physiological structures, Math. Bio. Eng., 20(3) (2023), 4493–4515.
  * [2] E. Beretta and D. Breda, Discrete or distributed delay? Effects on stability of population growth, Math. Biosci. Eng., 13 (2016), 19–41.
  * [3] J. P. Boyd, Chebyshev and Fourier Spectral Methods: Second Revised Edition, Dover Books in Mathematics, Dover Publications, Inc., Mineola, New York, second ed., 2001.
  * [4] D. Breda, S. De Reggi, F. Scarabel, R. Vermiglio, and J. Wu, Bivariate collocation for computing R0 in epidemic models with two structures, Computers & Mathematics with Applications, 116 (2022), 15–24.
  * [5] D. Breda, O. Diekmann, M. Gyllenberg, F. Scarabel, and R. Vermiglio, Pseudospectral discretization of nonlinear delay equations: new prospects for numerical bifurcation analysis, SIAM J. Appl. Dyn. Syst., 15 (2016), 1–23.
  * [6] D. Breda, O. Diekmann, S. Maset, and R. Vermiglio, A numerical approach for investigating the stability of equilibria for structured population models, J. Biol. Dyn., 7 (2013), 4–20.
  * [7] D. Breda, S. Maset, and R. Vermiglio, Stability of Linear Delay Differential Equations. A Numerical Approach with MATLAB, SpringerBriefs in Control, Automation and Robotics, Springer, New York, 2015.
  * [8] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations, Springer, New York, 2011.
  * [9] T. Cassidy, P. Gillich, A. Humphries, and C. van Dorp, Numerical methods and hypoexponential approximations for gamma distributed delay differential equations, IMA J. Appl. Math. (2022), hxac027.
  * [10] S. De Reggi, F. Scarabel, and R. Vermiglio, Approximating reproduction numbers: a general numerical method for age-structured models, preprint available at arXiv:2312.13477.
  * [11] B. de Wolff, F. Scarabel, S. Verduyn-Lunel, and O. Diekmann, Pseudospectral approximation of Hopf bifurcation for delay differential equations, SIAM J. Appl. Dyn. Syst., 20(1) (2021), 333–370.
  * [12] A. Dhooge, W. Govaerts, Y. A. Kuznetsov, H. G. E. Meijer, and B. Sautois, New features of the software MatCont for bifurcation analysis of dynamical systems, Math. Comput. Model. Dyn. Syst., 14 (2008), 147–175.
  * [13] O. Diekmann, P. Getto, and M. Gyllenberg, Stability and bifurcation analysis of Volterra functional equations in the light of suns and stars, SIAM J. Math. Anal., 39 (2007), 1023–1069.
  * [14] O. Diekmann and M. Gyllenberg, Equations with infinite delay: Blending the abstract and the concrete, J. Differ. Equ., 252 (2012), 819–851.
  * [15] O. Diekmann, M. Gyllenberg, and J. A. J. Metz, Finite dimensional state representation of linear and nonlinear delay systems, J. Dyn. Differ. Equ., 30 (2018), 1439–1467.
  * [16] O. Diekmann, M. Gyllenberg, and J. A. J. Metz, Finite dimensional state representation of pysiologically structured populations, J. Math. Biol., 80 (2020), 205–273.
  * [17] O. Diekmann, M. Gyllenberg, J. A. J. Metz, S. Nakaoka, and A. M. de Roos, Daphnia revisited: local stability and bifurcation theory for physiologically structured population models explained by way of an example, J. Math. Biol., 61 (2010), 277–318.
  * [18] O. Diekmann, F. Scarabel, and R. Vermiglio, Pseudospectral discretization of delay differential equations in sun-star formulation: results and conjectures, Discr. Cont. Dyn. Syst.-S, 13 (2020), 2575–2602.
  * [19] O. Diekmann, S. A. van Gils, S. M. Verduyn Lunel, and H.-O. Walther, Delay Equations. Functional-, Complex-, and Nonlinear Analysis, vol. 110 of Applied Mathematical Sciences, Springer, New York, 1995.
  * [20] O. Diekmann and S. Verduyn Lunel, Twin semigroups and delay equations, J. Diff. Equ., 286 (2021), 332–410.
  * [21] A. D’Onofrio and P. Manfredi, Behavioral SIR models with incidence-based social-distancing, Chaos, Solitons and Fractals, 159 (2022), 112072.
  * [22] K. Engelborghs, T. Luzyanina, and G. Samaey, DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations, Technical Report TW-330, Department of Computer Science, K.U. Leuven, Leuven, Belgium, 2001.
  * [23] C. Foias, G. R. Sell, and R. Temam, Inertial manifolds for nonlinear evolutionary equations, J. Diff. Equ., 73 (2) (1988), 309–353.
  * [24] W. Gautschi, Gauss–Radau formulae for Jacobi and Laguerre weight functions, Math. Comput. Simulation, 54 (2000), pp. 403–412,  1999 International Symposium on Computational Sciences, to honor John R. Rice (West Lafayette, IN).
  * [25] B.-Y. Guo, L.-L. Wang, and Z.-Q. Wang, Generalized Laguerre interpolation and pseudospectral method for unbounded domains, SIAM J. Num. Anal., 43 (2006), 2567–2589.
  * [26] W. Gurney, S. Blythe, and R. Nisbet, Nicholson’s blowflies revisited, Nature, 287 (1980), 17–21.
  * [27] M. Gyllenberg, F. Scarabel, and R. Vermiglio, Equations with infinite delay: Numerical bifurcation analysis via pseudospectral discretization, Appl. Math. Comp., 333 (2018), 490–505.
  * [28] J. K. Hale and S. M. V. Lunel, Introduction to Functional Differential Equations, vol. 99 of Springer Science & Business Media, 2013.
  * [29] G. Jordan, O. J. Staffans, and R. L. Wheeler, Subspaces of stable and unstable solutions of a functional differential equation in a fading memory space: The critical case, SIAM J. Math. Anal., 18 (1987), 1323–1340.
  * [30] Z. Liu and P. Magal, Functional differential equations with infinite delay in a space of exponentially bounded and uniformly continuous functions, Discr. Cont. Dyn. Syst.-B, 25 (2020), 2271–2292.
  * [31] P. Manfredi and A. D’Onofrio, Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases, Springer, New York, 2013.
  * [32] G. Mastroianni and G. V. Milovanović, Interpolation Processes. Basic Theory and Applications, Springer Monographs in Mathematics, Springer, Berlin, 2008.
  * [33] G. Mastroianni and G. Monegato, Some new applications of truncated Gauss–Laguerre quadrature formulas, Numerical Algorithms, 49 (2008), 283–297.
  * [34] J. A. J. Metz and O. Diekmann, The Dynamics of Physiologically Structured Populations, vol. 68 of Lecture Notes in Biomathematics, Springer, Berlin, Heidelberg, 1986.
  * [35] H. A. Priestley, Introduction to Complex Analysis. Oxford University Press, New York (1990).
  * [36] F. Scarabel, D. Breda, O. Diekmann, M. Gyllenberg, and R. Vermiglio, Numerical bifurcation analysis of physiologically structured population models via pseudospectral approximation, Vietnam J. Math., 49 (2021), 37–67.
  * [37] F. Scarabel, O. Diekmann, and R. Vermiglio, Numerical bifurcation analysis of renewal equations via pseudospectral approximation, J. Comp. Appl. Math., 397 (2021), 113611.
  * [38] F. Scarabel, L. Pellis, N. H. Ogden, and J. Wu, A renewal equation model to assess roles and limitations of contact tracing for disease outbreak control, R. Soc. Open Sci., 8:202091 (2021).
  * [39] J. Shen, T. Tang, and L.-L. Wang, Spectral Methods. Algorithms, Analysis and Applications, vol. 41 of Springer Series in Computational Mathematics, Springer, Berlin, Heidelberg, 2011.
  * [40] O. J. Staffans, Hopf bifurcation for functional and functional differential equations with infinite delay, J. Diff. Equ., 70 (1987), 114–151.
  * [41] J. A. C. Weideman and S. C. Reddy, A MATLAB differentiation matrix suite, ACM T. Math. Software, 26 (2000), 465–519.