DOI HERE 00 Advance Access Publication Date: Day Month Year Paper Published by
Oxford University Press on behalf of the Institute of Mathematics and its
Applications. All rights reserved.

Nicholas Hale.

[*]Corresponding author<EMAIL_ADDRESS>
0Year 0Year 0Year

# A spectral collocation method for functional and delay differential
equations

Nicholas Hale  Department of Mathematical Sciences, Stellenbosch University,
Stellenbosch, 7600, South Africa

(2024; 2023; Date; Date; Date)

###### Abstract

A framework for Chebyshev spectral collocation methods for the numerical
solution of functional and delay differential equations (FDEs and DDEs) is
described. The framework combines interpolation via the barycentric resampling
matrix with a multidomain approach used to resolve isolated discontinuities
propagated by non-smooth initial data. Geometric convergence is demonstrated
for several examples of linear and nonlinear FDEs and DDEs with various delay
types, including discrete, proportional, continuous, and state-dependent
delay. The framework is a natural extension of standard spectral collocation
methods based on polynomial interpolants and can be readily incorporated into
existing spectral discretisations, such as in Chebfun/Chebop, allowing the
automated and efficient solution of a wide class of nonlinear functional and
delay differential equations.

###### keywords:

Delay differential equation; spectral collocation; pseudospectral; barycentric
interpolation

##  1 Introduction

Delay differential equations (DDEs) arise in a variety of application areas,
notably population dynamics (Gopalsamy, 1992; Kuang, 1993), epidemiology
(Brauer and Castillo-Chávez, 2001; Rihan, 2021), control theory (Kyrychko and
Hogan, 2010; Michiels and Niculescu, 2007), and more (Erneux, 2009). Below we
review some basic concepts and refer the reader to any of the many excellent
introductory texts on DDEs and FDEs for details (Driver, 2012; Hale, 1971;
Kolmanovskii and Myshkis, 1999; Michiels and Niculescu, 2007)).

The most elementary DDEs are first-order initial value problems of the form

$\begin{array}[]{lllll}y^{\prime}(t)&=&f(t,y(t),y(\tau(t))),&0\leq t\leq T,\
\quad\tau(t)\leq t,\\\ y(t)&=&\phi(t),&t\leq 0,\end{array}$ (1.1)

with $\phi$, $f$, and $\tau$ given. More complex models may involve multiple
delays, systems of equations, higher-order derivatives, DDE boundary value
problems, DDE eigenvalue problems, or partial delay differential equations.
The most common types of delay are: discrete delay, where $\tau(t)=t-p$,
$p>0$; proportional delay, where $\tau(t)=qt$, $0<q<1$; continuous delay,
where $y^{\prime}(t)=f(t,y(t),\int^{t}_{0}y(\tau)f(t-\tau)\,d\tau)$; or some
collection/combination of the above. Furthermore, one may consider neutral-
type DDEs, whereby $y^{\prime}(t)=f(t,y(t),y(\tau(t)),y^{\prime}(\tau(t)))$,
or state-dependent delays, where $\tau=\tau(t,y(\tau))$. DDEs are a special
case of functional differential equations (FDEs), which relax the condition
that $\tau(t)\leq t$. For simplicity, we present most of the methodology in
this work in the context of DDEs. However the assumption $\tau(t)\leq t$ is
not required in the framework we propose, making it equally applicable for
DDEs and FDEs.

Like ordinary differential equations (ODEs), most DDEs and FDEs arising in
practical settings do not have closed-form solutions, so some kind of
approximation is required. Standard numerical methods for the approximating
DDE solutions typically involve Runge–Kutta-like solvers equipped with a
suitable interpolation to evaluate the delay terms (the method of steps), or
rewriting the DDE as a system of ODEs with additional variables to represent
the delayed terms. For an overview of such methods, we refer the reader to
(Bellen and Zennaro, 2013), (Engelborghs et al., 2002), or (Shampine and
Thompson, 2009), and references therein. A major difficulty arising in DDEs
and FDEs not present in ODEs is the propagation of discontinuities arising
from non-smooth initial data. Locating and resolving these discontinuities is
necessary in order to compute an accurate and efficient approximation, even
for low order methods (Baker et al., 1995).

For smooth ODEs, spectral collocation methods (also known as pseudospectral
methods) typically exhibit geometric convergence as the number of degrees of
freedom is increased, making them particularly efficient when a high degree of
accuracy is required (Boyd, 2001; Fornberg, 1996; Trefethen, 2000). Moreover,
spectral collocation methods extend naturally to boundary value problems and
differential-eigenvalue problems. Given the power of spectral methods in
solving ODEs, it is natural to apply them to DDEs and FDEs, and there are
numerous examples of this in the literature. We mention a few key papers here,
but more detailed overviews of the development of spectral methods in these
areas are given by Meng et al. (2019) and Jafari et al. (2021).

One of the first spectral methods for DDEs was introduced by Fox et al.
(1971), where a Chebyshev-tau approach was used to solve the linear pantograph
equation. Later, a Legendre-tau approach for DDEs with discrete delay was
presented by Ito et al. (1991). More recently there have been developments in
using spectral methods to solve continuous delay DDEs (Brunner and Hu, 2005,
2007), nonlinear DDEs (Ali et al., 2009), neutral and state-dependent DDEs
(Sedaghat et al., 2012), DDE differential eigenvalue problems (Liu and Ortiz,
1989), and periodic DDEs (Engelborghs et al., 2001; Butcher and Bobrenkov,
2011). In addition, spectral collocation methods have proven useful for
computing characteristic roots of DDEs (Breda et al., 2005; Wu and Michiels,
2012). The work most closely related to that which we present here is that of
Trif (2011), whose MATLAB software Chebpack (Trif, 2012) implements a
Chebyshev-tau approach allowing the ready solution of a variety linear
neutral-type DDEs with discrete and pantograph delay, including higher-order
DDEs and DDE eigenvalue problems. The aim of the present work is to describe a
unifying framework based on Chebyshev spectral collocation and the barycentric
resampling matrix (see Section 2.1) that can solve at once a variety and
combination of the types of DDE and FDE problems discussed above, and
furthermore in such a way that it can be easily incorporated in Chebfun
(Driscoll et al., 2014) to utilise its convenient syntax, nonlinear, and
adaptive solves. Rather than attempt to prove convergence results for some
limited class of problems, we shall demonstrate by numerous examples the
efficacy and generality of the proposed framework.

To that end, the outline of this paper is as follows. In Section 2 we provide
an overview of some necessary preliminaries, specifically the barycentric
interpolation formula and how it leads to the pseudospectral differentiation
and barycentric resampling matrices; our main tools in solving DDEs via
spectral collocation. In Section 3 we show how these matrices are combined to
solve some simple linear and nonlinear DDEs of various delay types. In Section
4 we discuss what small modifications are required so as to apply the method
to FDEs. In Section 5 we provide details of how the approach can be
incorporated in the Chebfun/Chebop framework for the automated solution of
DDEs/FDEs and demonstrate some more challenging examples solved via this
software. In Section 6 we discuss several possible extensions before
concluding Section 7.

##  2 Preliminaries

###  2.1 Barycentric interpolation formula

Given data $\\{(t_{1},y_{1}),\ldots,(t_{n},y_{n})\\}$ with $\\{t_{k}\\}$
finite and distinct, the unique polynomial interpolant of degree at most $n-1$
may be written in barycentric form as

$p(t)={\displaystyle\sum_{k=1}^{n}\frac{w_{k}y_{k}}{t-t_{k}}}\bigg{/}{\displaystyle\sum_{k=1}^{n}\frac{w_{k}}{t-t_{k}}},$
(2.1)

where the barycentric weights are given by $w_{k}=C_{n}\prod_{\ell\neq
k}(t_{k}-t_{\ell})^{-1}$ with $C_{n}$ an arbitrary nonzero constant (Berrut
and Trefethen, 2004, (4.2)).111$C_{n}$ can be chosen to normalise the weights,
for instance, such that $\max_{k}|w_{k}|=1$. Constant aside, any other choice
of weights will give rise to a rational interpolant which may or may not be
free of poles in the interpolation interval. Stability and conditioning of
(2.1) depends on the choice of nodes, $\\{t_{k}\\}$, which we discuss in
Section 2.2. For now, let us denote by ${\bm{t}}$ and $\bm{y}$ the
$\bm{t}=[t_{1},t_{2},\ldots,t_{n}]^{\top}$ and
$\bm{y}=[y_{1},y_{2},\ldots,y_{n}]^{\top}$. Since (2.1) is linear in the data
$\\{y_{k}\\}$, it can be expressed as $p(t)=P(t;\bm{t})\bm{y}$, where

$P(t;\bm{t})=\left[\begin{array}[]{ccccccc}\ell_{1}(t)&\ell_{2}(t)&\cdots&\ell_{n}(t)\end{array}\right]\
\ \text{and}\ \ \
\ell_{k}(t)={\displaystyle\frac{w_{k}}{t-t_{k}}}\bigg{/}{\displaystyle\sum_{k=1}^{n}\displaystyle\frac{w_{k}}{t-t_{k}}}.$
(2.2)

Here $P(t;\bm{t})=[\ell_{1}(t),\dots,\ell_{n}(t)]$ is a quasimatrix, an
$\infty$-by-$n$ matrix, whose columns contain the degree $n-1$ Lagrange
polynomials for the points $\\{t_{k}\\}$.

By differentiating the columns of $P(t;\bm{t})$ and evaluating the result at
$t=\bm{t}$, one obtains the pseudospectral differentiation matrix,
$D=P^{\prime}(\bm{t};\bm{t})$ (Berrut and Trefethen, 2004, Section 9.3). Code
for constructing such matrices is given in Table 1. Alternatively, given a
different set of points, say $\\{\tau_{j}\\}_{j=1}^{m}$,
$P(\bm{\tau};\bm{t})\bm{y}$ evaluates the degree $n$ polynomial interpolant of
the data $\\{(t_{k},y_{k})\\}_{k=1}^{n}$ at the $m$ points
$\bm{\tau}=[\tau_{1},\ldots,\tau_{m}]^{\top}$. Driscoll and Hale (2016) use
this barycentric resampling matrix with $m<n$ to downsample discretised linear
operators in spectral collocation, leaving space to apply boundary conditions
or other side constraints. Here we will use such matrices to evaluate delay
terms in DDEs. For instance, a delay term $y(qt-p)$ in a DDE can be
discretized as something like $P(q\bm{t}-p;\bm{t})\bm{y}$.222Here we assume
$p$ and $q$ are such that $q\bm{t}-p$ is contained within the interpolation
interval of (2.1). In general, one may need to subdivide $[0,T]$ and consider
multiple interpolants in order to achieve this. See Section 3.2 for details.
Concise Matlab code for constructing $P(\bm{\tau};\bm{t})$ is given in Table
1.

⬇

1function D = diffmat(t, w)

2 D = w.’./(w.*(t-t.’));

3 n = -

4 (t); ii = 1:n+1:n^2;

5 D(ii) = 0; D(ii) = -sum(D,2);

6end

⬇

1function P = barymat(tau, t, w)

2 P = w.’./(tau-t.’);

3 P = P./sum(P, 2);

4 P(isnan(P)) = 1;

5end

Table 1: Left: Matlab code to generate the pseudospectral differentiation
matrix, $D$. Similar code is available in a variety of software libraries,
such as diffmat in Chebfun (Driscoll et al., 2014) and poldiff in DMSUITE
(Weideman and Reddy, 2000). Right: Matlab code (modified from (Driscoll and
Hale, 2016)) to generate a barycentric resampling matrix, $P(\tau;\bm{t})$.

###  2.2 Interpolation at Chebyshev points

It is well known that interpolation at arbitrary nodes may be ill-conditioned.
For example, interpolating at uniformly-spaced nodes leads to the Runge
phenomenon and exponential ill-conditioning (Trefethen, 2019, Chapter 13).
Suitable nodes for polynomial interpolation cluster towards the endpoints of
the interpolation interval. Of such points, the Chebyshev–Gauss–Lobatto nodes,

$t_{k}=-\cos\left({(k-1)\pi}/(n-1)\right),\quad k=1,\ldots,n,$ (2.3)

(which hereafter we simply call Chebyshev points) on the interval $[-1,1]$ are
the most widely used in spectral methods for nonperiodic problems (Boyd, 2001;
Fornberg, 1996; Trefethen, 2000), and we focus on these in this manuscript.
However, the approach we describe is equally valid in the case of collocation
with, say, Gauss–Legendre or Gauss–Jacobi nodes, or even a rational spectral
method defined by some alternative choice of the barycentric weights (Berrut
and Klein, 2014). The barycentric weights for polynomial interpolation at the
Chebyshev points are known explicitly (Berrut and Trefethen, 2004):

$w_{k}=\left\\{\begin{array}[]{llll}1/2,&k=1,\\\
(-1)^{k+1},&k=2,\ldots,n-1,\\\ (-1)^{n+1}/2,&k=n.\\\ \end{array}\right.$ (2.4)

With these points and weights the barycentric formula (2.1) is numerically
stable and the interpolation problem is well-conditioned (Higham, 2004). If
the $\\{y_{k}\\}$ are sampled from a function $y(t)$ analytic in a
neighbourhood of $[-1,1]$, then the Chebyshev polynomial interpolant and its
derivatives converge geometrically to $y(t)$ as $n$ is increased. In
particular,
$\|p^{(\nu)}-y^{(\nu)}\|_{L^{\infty}_{[-1,1]}}=\mathcal{O}(\rho^{-n})$ as
$n\rightarrow\infty$, where $\rho$ is the radius of the largest Bernstein
ellipse with foci at $\pm 1$ that lies within the domain of analyticity of $y$
(Trefethen, 2019, Theorem 21.1).

Problems defined on finite intervals other than $[-1,1]$ can be handled by an
affine transformation. For instance, if solving on a domain $[0,T]$ then the
points $t_{k}$ are mapped by

$t_{k}\mapsto\tfrac{T}{2}(t_{k}+1),$ (2.5)

and the barycentric weights remain unchanged.

###  2.3 Spectral collocation for ODEs

The differentiation matrix, $D$, as computed by Table 1a, is the primary tool
for spectral collocation (Trefethen, 2000, Chapter 11). Higher-order
differentiation matrices can be computed by taking powers of $D$, although
there are more accurate and more efficient approaches of calculating higher-
order matrices directly (Baltensperger and Trummer, 2003), (Tee, 2006,
(3.28)). Non-constant coefficients can be incorporated by multiplying by a
diagonal matrix containing the values of the coefficient function evaluated at
the grid points. Nonlinear problems can be solved by a Newton–Raphson
iteration, and piecewise smooth problems by domain decomposition (both of
which are discussed below). Boundary conditions can be enforced in a number of
ways. The most common approach is boundary bordering (sometimes called row
replacement), whereby rows of the discretized operator are removed and
replaced with the boundary constraints (Boyd, 2001, Section 6.4), (Trefethen,
2000, p.135). An alternative approach is to project the discretized operator
onto a lower degree subspace via the barycentric resampling matrix and use the
boundary constraints to ‘square-up’ the problem (Driscoll and Hale, 2016). For
simplicity we employ the more traditional boundary bordering approach in this
work, but any reasonable method of enforcing the boundary conditions
(including rectangular projection) is equally applicable.

######  Example 1.

For comparison with the DDEs that follow, we briefly recall how a simple ODE
can be solved via spectral collocation. In particular, consider the elementary
differential equation

$y^{\prime}(t)=-y(t),\quad y(0)=1,$ (2.6)

with solution $y(t)=e^{-t}$. To solve, the $n\times n$ Chebyshev
differentiation matrix is constructed using the code in Table 1a and the ODE
discretized as $A\bm{y}=\bm{b}$ with $A=D+I$. Boundary conditions are enforced
by boundary bordering, i.e., removing the first rows of $A$ and $\bm{b}$ and
replacing them with the boundary constraint (here the first row of the
identity matrix and the scalar 1, respectively). The resulting square system
is then solved for $\bm{y}$, the approximate solution on the Chebyshev grid
$\bm{t}$. Concise code for solving (2.6) on the interval $[0,1]$ and a plot
demonstrating geometric convergence of the discretisation as $n$ is increased
are given in Figure 1.

⬇

1n = 12; y0 = 1;

2[t,~,w] = chebpts(n, [0,1]);

3z = zeros(n-1, 1);

4D = diffmat(t, w);

5I = eye(n); B = I(1,:);

6A = D + I; b = z;

7y = [B;A(2:end,:)]\\[y0;b(2:end)]

Figure 1: Left: MATLAB code for solving the ODE (2.6) using Chebyshev spectral
collocation. The line [t,~,w] = chebpts(n, [0,1]) (from Chebfun) returns the
n-point Chebyshev grid on the interval $[0,1]$ and the associated barycentric
weights (2.4). Centre: Exact and approximate solution using a 12-point
discretisation. Right: Geometric convergence of the solution as $n$ is
increased down to around machine precision. Here, and throughout (unless
otherwise stated), the reported error is the $L^{\infty}$ error at the
interpolation points.

##  3 Spectral collocation for DDEs

We now demonstrate how the Chebyshev spectral collocation method for ODEs can
be readily extended to solve DDE problems. The examples in this section are
intentionally simple, so as to clearly demonstrate the proposed method. In
particular, in most cases we choose a specific right hand-side so that the
solution is known (typically some form of exponential), allowing us to easily
compute the error in the approximation. More interesting and practical
examples are given in Section 5.

###  3.1 Proportional delay

We start by considering proportional (sometimes called pantograph-type) delay
of the form $y(qt)$, $q<1$ on the interval $[0,1]$.333Any other finite domain
can be dealt with via an affine change of variables, and discrete delays
introduced by such a change of variables can be dealt with in the manner
described in the next section. In this setting $qt\in[0,1]$ and no history
function or domain decomposition is required. Such a delay term can therefore
be implemented directly by the resampling matrix $P(q\bm{t};\bm{t})$.

######  Example 2.

Consider the DDE

$y^{\prime}(t)=-y(t)-y(t/2)+e^{-t/2},\quad y(0)=1,$ (3.1)

with solution $y=e^{-t}$.444A similar DDE appears in (Ali et al., 2009), where
it is solved using a Legendre collocation method. A stiff version of the
problem is considered by El-Safty and Hussien (1998), where a small parameter
is introduced in front of the $u^{\prime}$ term. Discretisation of the
differential and identity operators proceeds as in Example 1 and the delay
term is discretized using the barycentric resampling matrix as described
above. This results now in a system $A\bm{y}=\bm{b}$ where
$A=D+I+P(\bm{t}/2;\bm{t})$. As before, boundary conditions can be enforced by
row replacement. Code is provided in Figure 2 (left) and, as in the ODE case,
we observe geometric convergence as the discretization size is increased
(right). Multiple delays of proportional type, e.g., $y(q_{1}t)$, $y(q_{2}t)$,
can by included by incorporating additional resampling matrices,
$P(q_{1}\bm{t};\bm{t})$, $P(q_{2}\bm{t};\bm{t})$. Extensions to systems of
DDEs and/or nonlinear problems proceed as in the ODE case (Driscoll and Hale,
2016, Section 4).

⬇

1n = 12; y0 = 1;

2[t,~,w] = chebpts(n, [0,1]);

3D = diffmat(t, w);

4I = eye(n); B = I(1,:);

5P = barymat(t/2, t, w);

6A = D + I + P; b = -exp(t/2);

7u = [B;A(2:end,:)]\\[y0;b(2:end)];

Figure 2: Left: MATLAB code for solving the DDE (3.1) using Chebyshev spectral
collocation. The code is a simple modification of that from Figure 1 for
solving the ODE (2.6). Centre: Solution using 12-point Chebyshev grid. Right:
Convergence as $n$ is increased. As in the ODE case, geometric convergence to
around the level of machine precision is observed.

###  3.2 Discrete delay and domain decomposition

We next consider discrete delay. For simplicity we start only with problems
with trivial history functions, i.e., $\phi(t)=0$ for all $t<0$. Non-trivial
history functions can be dealt with in the manner described in the next
section.

A challenge arising with discrete delays (or with proportional delays on
intervals not including the origin) is that discontinuities caused by non-
smooth initial data will propagate into the solution whenever the delay
becomes active (Baker et al., 1995), and the rapid convergence of the spectral
collocation method will be lost (Trefethen, 2000, Chapter 4). If the delay (or
delays) are short relative to the length of the solution interval, then many
discontinuities will be introduced and the use of a spectral method is
unlikely to be beneficial. However, if only a few discontinuities are
introduced with locations known a priori, then in these can be dealt with by a
domain decomposition/multidomain approach. In particular, if the solution is
piecewise analytic on each subdomain, then geometric convergence may be
maintained, as in the case of ODEs (Driscoll and Fornberg, 1998). We describe
this approach first for a simple problem leading to only one singularity in
the domain of interest before outlining the more general approach. A similar
approach is employed by (Borgioli et al., 2020) in the context of stability
analysis of linear time-periodic delayed dynamical systems (see also Example
11).

######  Example 3.

Consider the DDE

$y^{\prime}(t)=-y(t)-y(t-\tfrac{1}{2}),\quad y(t<0)=0,\quad y(0)=1,$ (3.2)

with solution

$y=\left\\{\begin{array}[]{llll}e^{-t},&0\leq t\leq\tfrac{1}{2},\\\
e^{-t+1/2}(\tfrac{1}{2}-t+{e^{-1/2}}),&\tfrac{1}{2}\leq t\leq
1.\end{array}\right.$ (3.3)

From the delay term $y(t-\tfrac{1}{2})$, we know a priori that the solution
will contain a jump in its first derivative at $t=\tfrac{1}{2}$ when the delay
becomes active (and its second, third, fourth, $\ldots$ derivatives at
$t=1,\tfrac{3}{2},2,\ldots$ if we were to integrate further – see below). As
such, we discretise on a piecewise domain with two Chebyshev grids,
$\bm{t}_{L}\in[0,\tfrac{1}{2}]$ and $\bm{t}_{R}\in[\tfrac{1}{2},1]$, and seek
an approximate solution $\bm{y}=[\bm{y}_{L};\bm{y}_{R}]$ on each grid. A
discretized differential operator is constructed for each subinterval and a
continuity condition $[-\bm{e}_{n_{L}}^{\top},\bm{e}_{1}^{\top}]\bm{y}=0$
enforced so that the solution is continuous at $t=\tfrac{1}{2}$.555For higher-
order problems, continuity of additional derivatives is enforced at the
interface by including appropriate rows of suitable differentiation matrices.
Combined with the initial condition and boundary bordering, this results in a
square system that can be solved for $\bm{y}$. Excluding the delay term,
details of such domain decomposition are described in detail in (Driscoll and
Hale, 2016, Section 4.3).

To incorporate the delay term, we first note that $y(t-\tfrac{1}{2})=0$ for
$t<\tfrac{1}{2}$, and so this term may be neglected in the first subinterval.
For any $t$ in the second subinterval, $[\tfrac{1}{2},1]$, we note that
$t-\tfrac{1}{2}\in[0,\tfrac{1}{2}]$, thus the solution on the first
subinterval is interpolated to the appropriate points by the resampling matrix
$P(\bm{t}_{R}-\tfrac{1}{2},\bm{t}_{L})$.666If the width and the discretization
size of the two intervals is the same, then this will simply be the identity
matrix, but in general the resampling matrix will be nontrivial. See Remark 1
for further discussion. This gives rise to the block linear system in the left
panel of Figure 3, where again we use boundary bordering to
incorpomichiels2020rate the boundary and continuity conditions. For this
example, and indeed any initial-value problem DDE, the resulting system is
block lower-triangular and one can solve first for $\bm{y}_{L}$ and then for
$\bm{y}_{R}$. In the context of solving DDEs, this is referred to as the
method of steps. Whilst this structure could certainly be exploited here, FDEs
and BVP-DDEs can be implemented using the same framework as described above,
and these will not result in block lower-triangular systems. For simplicity
and generality, we therefore simply solve the resulting block discretisation
all at once using standard direct methods and defer exploration of more
efficient domain decomposition strategies for future work. MATLAB code for
reproducing the solution to (3.2) (and all other examples in this paper) is
available at the online repository: Hale (2024).

$\left[\\!\\!\\!\begin{array}[]{cccc}\bm{e}_{1}^{\top}\\\\[1.00006pt]
\overline{D_{L}}+\overline{I}&\\\
-\bm{e}_{n_{L}}^{\top}&\\!\\!\\!\bm{e}_{1}^{\top}\\\\[1.00006pt]
\overline{P}(\bm{t}_{R}-\tfrac{1}{2};\bm{t}_{L})&\\!\\!\\!\overline{D_{R}}+\overline{I}\end{array}\\!\\!\\!\right]\\!\left[\\!\\!\\!\begin{array}[]{cccc}\bm{y}_{L}\\\
\bm{y}_{R}\end{array}\\!\\!\\!\right]=\left[\\!\\!\begin{array}[]{cccc}1\\\
\bm{0}\\\ 0\\\ \bm{0}\end{array}\\!\\!\\!\right]$

Figure 3: Left: Block linear system resulting from the spectral collocation
discretisation of (3.2). The overline notation is used to indicate that the
first row of the corresponding matrix is removed in order to implement
boundary/continuity constraints via boundary bordering. The system is block-
lower triangular, and hence one could first solve for $\bm{y}_{L}$ and then
$\bm{y}_{R}$. However, this will not be the case in context of FDEs as we
discuss shortly. Centre: Solution using 10- and 11-point Chebyshev grids on
$[0,1/2]$ and $[1/2,1]$, respectively. Right: Convergence as $n$ is increased,
where $n$ and $n+1$ are the discretisation sizes of the left and right
subintervals, respectively. Using a multidomain discretisation on
$[0,1/2]\cup[1/2,1]$ so that the solution is analytic on each subdomain
results in geometric convergence (solid line). If only a single Chebyshev grid
on $[0,1]$ is used, then algebraic convergence is observed (dashed line) due
to the discontinuity in the derivative of the solution at $t=1/2$.

The other two panels of Figure 3 show the computed piecewise smooth solution
$\bm{y}=[\bm{y}_{L};\bm{y}_{R}]$ (centre) and the convergence of the method as
the discretization size is increased (right – solid line). Since the solution
is piecewise analytic on the two subintervals, geometric convergence is
obtained, here using an $n$-point discretisation on the left subinterval and
an ($n+1$)-point discretisation on the right. Finally, the dashed line in the
third figure shows the result of applying a similar approach to discretise
this DDE but without introducing the breakpoint at $t=\tfrac{1}{2}$. As
expected, the jump in $y^{\prime}$ at this point results in only algebraic
convergence of the numerical solution. As with most numerical methods for DDEs
of this type, knowledge of the discontinuity locations is vital.

######  Remark 1.

The choice of using a different discretisation size on each subinterval in the
example above is certainly not necessary. We do so to demonstrate that this
offers no difficulty to the method. For this example, choosing an $n$-point
grid on both subintervals would result in
$\overline{P}(\bm{t}_{R}-\tfrac{1}{2};\bm{t}_{L})$ simplifying to an identity
matrix, thus providing a sparser linear system. However, the ability to use a
more general interpolation operator not only allows for differing refinement
on each of the subintervals, but is necessary when subintervals are not of a
uniform size, which is typically the case when there are two or more discrete
delays present or the delay is nonlinear – see Example 5.

Consider now a more general situation and suppose that the propagation of
singularities requires that the solution interval $[0,T]$ is decomposed into
$m$ subintervals with breakpoints $0=T_{0}<T_{1}<\ldots<T_{m}=T$, and each
subinterval $[T_{k-1},T_{k}]$ is discretised using an $n_{k}$-point Chebyshev
grid, $\bm{t}_{k}$, for $k=1,\ldots,m$. An evaluation/delay operator
$y(\tau(t))$ can then be realised as a block matrix

$P=\left[\begin{array}[]{cccc}P_{1,1}&P_{1,2}&\cdots&P_{1,m}\\\
P_{2,1}&P_{2,2}&\cdots&P_{2,m}\\\ \vdots&\vdots&\ddots&\vdots\\\
P_{m,1}&P_{m,2}&\cdots&P_{m,m}\\\ \end{array}\right],$

where each $P_{j,k}$ is an $n_{j}\times n_{k}$ matrix defined by

$[P_{j,k}]_{\ell,:}=\left\\{\begin{array}[]{llll}P(\bm{[\tau}_{j}]_{\ell},\bm{t}_{k}),&T_{k}<[\bm{\tau}_{j}]_{\ell}\leq
T_{k+1},\quad\text{where }\bm{\tau}_{j}=\tau(\bm{t}_{j}).\\\
\bm{0}^{\top},&\rm{otherwise}.\end{array}\right.$

Note that, as in Example 3, if $\tau(t)\leq t$ then $P_{j,k}=0$ for $k>j$,
i.e., $E$ is block lower-triangular. However, for FDEs this will not be the
case in general.

######  Example 4.

We repeat Example 3, but now solve on the domain $[0,2]$. The propagated
singularities appear at $t=0.5,1,1.5$, and so we decompose $[0,2]$ into 4
uniform subintervals. The results are shown in Figure 4.

$\left[\\!\\!\\!\begin{array}[]{ccccccccccc}\bm{e}_{1}^{\top}\\\
\overline{D_{1}}+\overline{I}&\\\
-\bm{e}_{n_{1}}^{\top}&\\!\\!\\!\bm{e}_{1}^{\top}\\\
\overline{P_{2,1}}&\\!\\!\\!\overline{D_{2}}+\overline{I}\\\
&-\bm{e}_{n_{2}}^{\top}&\bm{e}_{1}^{\top}\\\
&\overline{P_{3,2}}&\overline{D_{3}}+\overline{I}&\\\
&&-\bm{e}_{n_{3}}^{\top}&\\!\\!\\!\bm{e}_{1}^{\top}\\\
&&\overline{P_{4,3}}&\\!\\!\\!\overline{D_{4}}+\overline{I}\end{array}\\!\\!\\!\right]\\!\left[\\!\\!\\!\begin{array}[]{cccc}\bm{y}_{1}\\\
\bm{y}_{2}\\\ \bm{y}_{3}\\\
\bm{y}_{4}\end{array}\\!\\!\\!\right]=\left[\\!\\!\begin{array}[]{cccc}1\\\
\bm{0}\\\ {0}\\\ \bm{0}\\\ 0\\\ \bm{0}\\\ {0}\\\
\bm{0}\end{array}\\!\\!\\!\right]$

Figure 4: Left: Block discretisation of the DDE on the interval $[0,2]$ with
‘breakpoints’ at $0.5,1,1.5$. Here $P_{j,k}$ is given by
$P(\bm{t}_{j}-\tfrac{1}{2};\bm{t}_{k})$. Centre: Solution using 10-, 11-, 12-,
13-point grids on the respective subintervals. Right: Geometric convergence as
discretisation size is increased ($n$, $n+1$, $n+2$, $n+3$-point
dicretisations on each subinterval, respectively.)

Although the solution may look smooth at $t=1$ and $t=1.5$, there are
discontinuities in the 2nd and 3rd derivatives of the solution at these
points, respectively. However, since the solution is piecewise analytic
between these points, we again observe geometric convergence in the numerical
solution when using 4 subintervals with the discontinuity locations as
breakpoints. Note again that the system is block lower triangular, so one can
emulate the method of steps in solving first for $\bm{y}_{1}$, then
$\bm{y}_{2}$, and so on.

###  3.3 Non-trivial history functions

In many applications the DDE history function, $\phi(t)$, is constant and
consistent with the initial condition so that $y(t\leq 0)=y_{0}$. In this
case, one can simply consider an equivalent delay
$\tau(t)\mapsto\max(\tau(t),0)$ and enforce a constant history function in a
similar manner as described in the previous section. More generally, at least
for state-independent delays, history functions can be dealt with by rewriting
the DE in a piecewise-defined form with a constant history function. In
particular, a DDE in the form of (1.1) can be rewritten as

$y^{\prime}(t)=\left\\{\begin{array}[]{lll}f(t,y(t),\phi(t)),&\tau(t)<0,\\\
f(t,y(t),y(\tau(t)),&0\leq\tau(t),\\\ \end{array}\right.$ (3.4)

with $y(0)=y_{0}$. Multiple and neutral-type delays can be dealt with in a
similar manner. To ensure exponential convergence, a multidomain approach with
breakpoints included whenever these delays become active can be employed, as
described in the previous section.

###  3.4 Nonlinear delay

The approaches described above can be applied to arbitrary smooth delays
functions $\tau(t)\leq t$. In such instances, locating the propagated
discontinuities can be challenging, but in principle it can be done (Neves and
Feldstein, 1976). With these locations known, the domain can be subdivided and
multidomain approach used as described above.

######  Example 5.

We demonstrate on a simple example with the nonlinear delay
$\tau(t)=t^{2}-1/4$. In particular, consider the DDE

$y^{\prime}(t)=-y(t)-y(t^{2}-1/4),\quad y(t<0)=0,\quad y(0)=1.$ (3.5)

Here, as in the previous example, a discontinuity in the derivative of the
solution is introduced at $t=1/2$, where $\tau(t)=0$. There will then be an
additional discontinuity introduced in the second derivative at
$t=\sqrt{3}/2$, where $\tau(t)=1/2$. In this manner, all discontinuities can
be traced out, even for non-constant delays. For this example, if we solve for
$0\leq t\leq 1$, these are the only two discontinuities required. (The next
occurs at $t=\tfrac{1}{2}\sqrt{2\sqrt{3}+1}\approx 1.05$.) Since the
subintervals $[0,1/2]$, $[1/2,\sqrt{3}/2]$, and $[\sqrt{3}/2,1]$ are of
differing widths, even if the same discretisation size is used on each the
resampling matrices are not trivial. We therefore discretise using the
approach outlined in Section 3.2. The left panel of Figure 5 shows the result
of solving with an $n=12$-point discretisation on each subinterval. The centre
panel demonstrates geometric convergence as $n$ is increased.777Using the
method of steps, an explicit (but lengthy) solution to (3.5) can be expressed
in terms of the integral of an error function. We omit the details. The final
panel shows the sparsity plot of the $n=12$ discretisation.

Figure 5: Left: Analytic and computed solution to (3.5) on a 12-point grid.
Centre: Geometric convergence as the discretisation size is increased. Right:
Sparsity plot of the 12-point Chebyshev spectral discretisation. Here, even
though the same value of $n$ is used on each subinterval, the interpolation
operator is non-trivial (because the intervals are of differing widths),
resulting in dense off-diagonal blocks of the discretisation.

###  3.5 Nonlinear problems & state-dependent delay

In spectral methods for ODEs, nonlinear problems are typically solved by a
Newton–Raphson iteration, or some variant thereof–see (Driscoll and Hale,
2016, Section 4.2) for details. DDEs allow for an additional form of
nonlinearity, whereby the nonlinearity occurs in the delay. This is known as
state-dependent delay (Hartung et al., 2006).

Suppose a term of the form $y(g(y))$ appears in the DDE, where $g$ is some
differentiable function defined on the range of $y$. To linearise, we make the
following observation:

$\frac{d}{dy}\bigg{(}P(g(y))y\bigg{)}=P(g(y))+g^{\prime}(y)P^{\prime}(g(y))y=P(g(y))+g^{\prime}(y)P(g(y))\frac{dy}{dt},$
(3.6)

where we have omitted the $\bm{t}$ dependence of $P(\cdot;\bm{t})$ and the $t$
dependence of $y(t)$ for brevity. Thus the discrete Jacobian of the
interpolation operator, $\frac{d}{d\bm{y}}P(g(\bm{y});\bm{t})$ , can be
computed by

$\frac{d}{d\bm{y}}P(g(\bm{y}))\bm{y}=P(g(\bm{y}))+\text{diag}\Big{(}g^{\prime}(\bm{y})\circ\big{(}P(g(\bm{y}))D\bm{y}\big{)}\Big{)}.$
(3.7)

This can then be incorporated in a Newton–Raphson iteration to solve state-
dependent problems.

A detailed description of the existence and uniqueness of solutions to state-
dependent DDEs is beyond the scope of this paper. Furthermore, the propagation
of discontinuities is significantly more complex than in the constant delay
case (Neves and Feldstein, 1976). We put aside these challenges for now, and
instead consider below a reasonably simple case for which a known (smooth)
solution exists. A less contrived example of state-dependent delay is
demonstrated in Section 5.

######  Example 6.

Consider the DDE

$y^{\prime}(t)=-y(y(t))+f(t),\quad t\in[0,1],\qquad y(0)=0.$ (3.8)

We employ the method of manufactured solutions, taking
$f(t)=\cos(t)+\sin(\sin(t))$ so that $y(t)=\sin(t)$. Figure 6 demonstrates the
result of applying the approach outlined above to the DDE (3.8), using an
initial guess $y(t)=t$ in the Newton–Raphson iteration. The table on the left
shows quadratic convergence of the iterations in both the norm of the residual
and the norm of the Newton update for an $n=12$-point discretisation. As with
previous examples, the centre panel shows the resulting computed solution and
the rightmost panel demonstrates geometric convergence as $n$ is increased.

    Newton iteration:
    iter  norm(res, inf)  norm(dy, inf)
      0   0.71407355247   0.26232516612
      1   0.05480002458   0.01314905164
      2   0.00016794991   0.00002292528
      3   0.00000000051   0.00000000004

Figure 6: Computed solution to the state-dependent delay problem, (3.8). Left:
Quadratic convergence of the Newton iteration applied to the $n=12$
discretisation with initial guess $y=t$. Centre: Computed solution taking
$n=12$. Right: Geometric convergence as $n$ is increased.

######  Remark 2.

A careful reader may have noticed that in the above Example we carefully chose
$y(t)=\sin(t)$ so that the state-dependent delay satisfied $\tau(t,y(t))\leq
t$, and was indeed delay differential equation. However, as with previous
examples, this is not a requirement of the method, and a state-dependent
problem with $\tau(t,y(t))\not<t$ is considered in Section 4.

###  3.6 Continuous delay

Continuous (sometimes called distributed) delays result in Volterra-type
integro-differential equations. Driscoll (2010) presents a detailed
description of a Chebyshev spectral collocation method for such equations, and
we refer the reader to that work for further information. We include an
example here both for completeness and to demonstrate that, since both methods
are natural extensions to standard Chebyshev spectral methods for ODEs, they
can be effortlessly combined to solve problems with a combination of both
discrete, proportional, and continuous delay types.

######  Example 7.

Consider the DDE below, which combines both proportional and continuous delay:

$y^{\prime}(t)+\tfrac{1}{2}y(t/2)=\int_{0}^{t}y(s)e^{-(t-s)^{2}}\,ds,\qquad
y(0)=1.$ (3.9)

The code in the left panel of Figure 7 demonstrates compact MATLAB code for
solving this problem, combining the resampling approach of Section 3.1 with
code from (Driscoll, 2010) for discertising the Volterra operator. The
computed solution is depicted in the centre panel and geometric convergence as
the discretisation size is increased is demonstrated on the right.

⬇

1K = @(s,t) exp(-(t-s).^2);

2n = 18; dom = [0 1];

3t = chebpts(n, dom);

4D = diffmat(n, dom);

5P = barymat(t/2, t);

6V = K(t, t’).*cumsummat(n, dom);

7A = D + P/2 - V ; rhs = 0*t;

8A(1,:) = eye(1,n); rhs(1) = 1;

9y = A\rhs;

Figure 7: Left: MATLAB code for solving the combined proportional and
continuous DDE (3.9). The Chebfun function cumsummat(n,dom) computes an n-by-n
spectral collocation matrix representing indefinite integration. Centre:
Computed solution on an 18-point grid. Right: Geometric convergence as the
discretisation size, $n$, is increased. Here the $n=30$ discretisation is used
a proxy for the unknown true solution.

###  3.7 Other types of DDE problems

We have described above how the barycentric resampling approach can be
incorporated into standard spectral collocation methods to solve the most
common types of DDE problems. The approach can be further generalised to
neutral-type DDEs, with a discrete operator $P(\bm{\tau},\bm{t})D$
representing $y^{\prime}(\tau(t))$, as well as higher-order problems, DDE
boundary-value problems, systems of DDEs, DDE eigenvalue problems, and partial
delay different equations using natural extensions of the techniques described
in (Driscoll and Hale, 2016, Sections 4–6.) for ODEs. Examples of some of
these are given in Section 5 below.

##  4 Spectral collocation for functional differential equations

Unlike time-marching methods, the approach introduced in Section 3 did not
require $\tau(t)\leq t$, making it equally applicable to more general function
differential equations (FDEs). As with DDEs, one must account for propagation
of discontinuities caused by non-smooth initial data, as well as the potential
need for an additional ‘history’ function for $t>T$ (Hale, 1971). For the sake
of brevity we avoid these complications here, which can be treated in a manner
similar to those described in the previous section, and instead consider two
examples for which $\tau(t)\in[0,T]$.

######  Example 8.

Consider the following FDE:

$y^{\prime}(t)=-y(t)-y(1-t^{2})+f(t),\quad t\in[0,1],\qquad y(0)=1,$ (4.1)

where here $f(t)$ is chosen as $e^{t^{2}-1}$ so that (4.1) has the exact
solution $y(t)=e^{-t}$. Clearly the $1-t^{2}$ in the argument of the second
term on the right hand-side does not satisfy $1-t^{2}<t$ for all $0\leq t\leq
1$, making this an FDE but not a DDE. However, as shown in Figure 8,
discretisation proceeds almost identically to the DDE (3.1) and geometric
convergence is again observed as the discretisation size is increased.

⬇

1n = 12; y0 = 1;

2[t,~,w] = chebpts(n, [0,1]);

3D = diffmat(t, w);

4I = eye(n); B = I(1,:);

5P = barymat(1-t.^2, t, w);

6A = D + I + P; b = exp(t.^2-1);

7u = [B;A(2:end,:)]\\[y0;b(2:end)];

Figure 8: Left: MATLAB code for solving the FDE (4.1) using Chebyshev spectral
collocation. The only change from solving the DDE (3.1) is the first argument
of barymat on line 6 and the definition of the righ hand-side vector, b.
Centre: Solution using 12-point Chebyshev grid. Right: Convergence as $n$ is
increased. As in the DDE case, geometric convergence to around the level of
machine precision is observed.

######  Example 9.

Here we consider a FDE variant of the DDE in Example 6 in which $y(t)\not<t$,
namely888This FDE was the original motivation for this work after a student in
the author’s undergraduate class became confused when solving the much more
simple ODE $y^{\prime}=-y^{2}$, $y(0)=1$ and asked if one could solve such a
“crazy-wicked” equation as (4.2)!

$y^{\prime}(t)=-y(y(t)),\qquad y(0)=1.$ (4.2)

From the initial condition it follows that $y(t)<0$ and it turns out that
$y(t)\in[0,1]$ for $t\in[0,1]$, meaning that no history functions or domain
subdivisions are required here. Again, the fact that this is an FDE rather
than a DDE has almost no impact on the discretisation, and the implementation
proceeds almost identically to that of Example 6. Figure 9 shows the
convergence of the Newton iteration, the computed solution, and its
convergence.

    Newton iteration:
    iter    norm(res)        norm(dy)
      0   1.00000000000   1.075290658380
      1   0.25000000000   0.159726357356
      2   0.00686128071   0.002791677486
      3   0.00000843021   0.000005995919
      4   0.00000000002   0.000000000006

Figure 9: Computed solution to the state-dependent FDE (4.2). Left: Quadratic
convergence of the Newton iteration applied to the $n=12$ discretisation with
initial guess $y=1$. Centre: Computed solution taking $n=12$. Right: Geometric
convergence as $n$ is increased, taking the $n=30$ solution as a proxy for the
unknown exact solution.

##  5 Chebfun & more examples

Chebfun is an open-source software system for numerical computation with
functions (Driscoll et al., 2014). Its Chebop class and associated methods
allow Chebfun to solve a variety of linear and nonlinear ODE boundary-value
problems through the use of adaptive spectral collocation, automatic
differentiation, and robust Newton–Raphson iteration (Driscoll et al., 2008;
Birkisson and Driscoll, 2012). Chebfun can also solve integral equations
(Driscoll, 2010), differential eigenvalue problems (Driscoll et al., 2008), as
well as initial-value PDEs through the use of method-of-lines type
discretisations (Trefethen et al., 2017, Chapter 22).

The framework outlined in this manuscript uses the same underlying Chebyshev
spectral collocation discretisation as Chebop uses for ODEs, allowing the
natural extension of Chebop to solve linear FDEs. Furthermore, Chebop’s
automatic differentiation capabilities are easily modified to allow
linearisation of state-dependant delays using (3.7), thus also allowing the
solution of nonlinear FDEs. A detailed description of the software engineering
involved in this extension is beyond the scope of this work, but the syntax,
usability, and versatility of using Chebfun/Chebop to solve DDEs is
demonstrated in various examples below. In addition, Chebfun implementations
of the examples from Section 3 can be found at the repository (Hale, 2024).

###  5.1 Chebfun example 1 (proportional and continuous delay):

Our first Chebfun example is a linear DDE with both proportional and
continuous delay, namely

$y^{\prime}(t)=\frac{1}{100}(qt-t-10)y(qt)+\frac{1}{100}(t+20)e^{-1}+\frac{1}{100}\int_{0}^{t}y(s)\,ds+\frac{1}{1000}\int_{0}^{qt}(t-\tau)y(\tau)\,d\tau,\quad
t\in[0,20],$ (5.1)

with $y(0)=e^{-1}$, with exact solution $y(t)=e^{\frac{1}{10}t-1}$. This
example appears in (Wei and Chen, 2012), where it is solved using a Legendre-
based spectral collocation method. The short code in Figure 10 demonstrates
the ease with which this DDE can be solved in Chebfun. Here the computed
solution, taking $q=\frac{1}{2}$, is returned as a length 14 chebfun and has a
relative error of around $1.29\times 10^{-15}$.

⬇

1 dom = [0 20]; q = 1/2;

2 N = chebop(@(t,y) diff(y) - 1/100*(q*t-t-10)*y(q*t) - ...

3 1/100*(t+20)*exp(-1) - 1/100*cumsum(y) - ...

4 1/1000*feval(volt(@(x,s)(x/q-s), y), q*t), dom);

5 N.lbc = exp(-1);

6 y = N\0;

7 sol = chebfun(@(t) exp(t/10-1), dom);

8 err = norm(y-sol, inf)./norm(sol, inf)

Figure 10: Left: Chebfun code to solve the DDE (5.1). Right: Computed solution
and relative error.

###  5.2 Chebfun example 2 (nonlinear, neutral type):

Our next example is a nonlinear DDE of neutral type:

$y^{\prime}(t)=2\cos(2t)y(t/2)^{2\cos(t)}+\log(y^{\prime}(t/2))-\log(2\cos(t))-\sin(t),\quad
y(0)=1.$ (5.2)

This DDE was introduced by Jackiewicz (1984) and appears in the MATLAB example
ddex5 where it solved using the built-in MATLAB routine for neutral DDEs,
ddensd (MathWorks, 2012). The DDE has two consistent initial slopes, $s=2$ and
$s=-W(-2e^{-2})=0.406\ldots$, where here $W$ is the Lambert W function, each
corresponding to a different solution (Shampine, 2008). Chebfun code for
computing both solutions is given in the listings a the top of Figure 11. When
$s=2$, the DDE has an explicit solution $y=e^{\sin(2t)}$, and both the Chebfun
and ddensd show good agreement with this (with relative errors of around
$10^{-14}$ and $10^{-5}$, respectively), as shown in the left panel of Figure
11. For $s=-W(-2e^{-2})$ there is no known closed form solution. In the right
panel of Figure 11 we see that the numerical solututions computed by Chebfun
and ddesnd differ considerably, even over a relatively short time interval.
This author believes the cause of this is that the DDE is highly ill-
conditioned in this case, but that this is mitigated by the fully implicit
(i.e., global) nature and high accuracy of the interpolation approach
implemented in Chebfun. We defer further investigation to future work.

⬇

1N = chebop(@(t,y) diff(y) - (2*cos(2*t)*y(t/2)^(2*cos(t)) + ...

2 log(deriv(y,t/2)) - log(2*cos(t)) - sin(t)), [0,.1]);

3N.lbc = 1;

4for s = [2, -lambertw(-2/exp(2))]

5 N.init = chebfun(@(t) 1+s*t, [0,.1]);

6 u = N\0;

7 figure, plot(u), title([’s = ’ num2str(s, 14)])

8end

Figure 11: Top: Chebfun code to solve (5.2). The Chebfun function deriv(y,t/2)
evaluates $y^{\prime}$ at $t/2$. Bottom left: Agreement of Chebfun and ddesnd
solutions when $s=2$. Bottom right: Disagreement of Chebfun and ddesnd
solutions when $s=W(-2e^{-2})$.

###  5.3 Chebfun example 3 (functional, state-dependent)

Here we simply show how Example 4.2 can be implemented as a three-line code in
Chebfun. The solution is depicted in Figure 6.

⬇

1N = chebop(@(t,y) diff(y) + y(y), [0 1]);

2N.lbc = 1;

3y = N\0;

Figure 12: Chebfun code to solve (4.2).

###  5.4 Chebfun example 4 (unknown delay)

Here we revisit Example 2 from Section 3.1. However, rather than specify the
proportional delay, $p$, we consider this an unknown which must be solved in
order to satisfy the a two-point DDE boundary-value problem

$y^{\prime}(t)=-y(t)-y(pt)+e^{-t/2},\qquad y(0)=1,\ y(1)=0.25.$ (5.3)

Here we employ Chebop’s functionality for solving systems of equations
containing unknown parameters without the need to introduce extra equations
into the system. With $p$ unknown, the problem is nonlinear, and Chebfun’s
Newton–Raphson iteration with default initial guess and parameters converges
in around 5 iterations. The computed solution (represented by a degree 10
chebfun) and the determined value of $p$ are shown in Figure 13.

⬇

1N = chebop(@(t,y,p) diff(y) + y + y(p*t) - exp(-t/2), [0 1]);

2N.lbc = @(y,p) y - 1;

3N.rbc = @(y,p) y - 0.25;

4yp = N\0;

5plot(yp{1}), title([’p = ’ num2str(yp{2}, 10)])

Figure 13: Left: Chebfun code to solve (5.3). Right: Computed solution.

###  5.5 Chebfun example 5 (eigenvalue problem)

As our final example in this section we consider a delay differential
eigenvalue problem.999The most common eigenvalue problem associated with
functional or delay differential equations is that of obtaining the
characteristic roots of the Jacobian at critical points to determine their
stability properties. This is not what are doing here. Chebfun solves
differential eigenvalue problems by computing an increasingly large Chebyshev
spectral discretisation of the differential operator, as described in
(Driscoll and Hale, 2016, Section 5), and looking for the ‘smoothest’
eigenvalues. With the barycentric resampling approach, this is easily adapted
to the delay or functional case. Here we consider the DDE eigenvalue problem

$y^{\prime\prime}(t)=\lambda y(t/2),\qquad y(0)=y(1)=0.$ (5.4)

Chebfun code is given in Figure 14, in addition to the computed eigenvalues
and eigenfunctions.

⬇

1A = chebop(@(y) diff(y,2), [0 1]);

2A.bc = ’dirichlet’;

3B = chebop(@(t,y) -y(t/2), [0 1]);

4[V, D] = eigs(A, B);

5fprintf(’%10.10g\n’, diag(D))

6plot(V)

⬇

113.05485001

2169.7286494

31398.543635

49480.135715

557516.69512

6326708.2900

Figure 14: Left: Chebfun code. Centre: Computed eigenvalues: Right:
Eigenfunctions.

##  6 Extensions

###  6.1 Laguerre and Hermite spectral methods

Laguerre and Hermite spectral collocation methods employ weighted-polynomial
bases, $p(t)e^{-t/2}$ and $p(t)e^{-t^{2}/2}$, respectively. The approach
described in Sections 3 and 4 can therefore be naturally extended to
functional and delay differential defined on semi-infinite and infinite
domains. For example, if $\bm{t}=[t_{1},\ldots,t_{n}]^{\top}$ is an $n$-point
Gauss–Laguerre (or Gauss–Laguerre–Lobatto) grid and
$\bm{w}=[w_{1},\ldots,w_{n}]^{\top}$ are the associated barycentric weights ,
then a functional operator $y(qt+p)$ can be discretised as

${\rm diag}(e^{-(q\bm{t}+p)/2})P(q\bm{t}+p;\bm{t}){\rm
diag}(e^{\bm{t}/2})\bm{y}.$

This approach has recently been employed in (Hale et al., 2024) in the
numerical investigation of the nonlinear functional differential equation

$u^{\prime}(t)+u(t)=u^{2}(\alpha t),\quad\alpha>0,\ t>0,\qquad u(0)=1.$

###  6.2 Fourier spectral methods

Thus far we have focussed on spectral methods based on polynomial
interpolation. However, similar ideas can be applied when solving periodic
DDEs using trigonometric interpolants. In particular, the degree $n$
trigonometric polynomial interpolant of an arbitrary function at $2n$
equidistant nodes, $\theta_{k}=k\pi/n$, $k=0,\ldots,2n-1$ on the interval
$[0,2\pi]$ may be expressed in barycentric form as (Henrici, 1979)

$t_{n}(\theta)=\sum_{k=0}^{2n-1}(-1)^{k}\cot\frac{\theta-\theta_{k}}{2}y_{k}\bigg{/}\sum_{k=0}^{2n-1}(-1)^{k}\cot\frac{\theta-\theta_{k}}{2},$
(6.1)

and a trigonometric barycentric resampling matrix constructed as in the left
hand-side of Figure 2. This matrix can be combined with the trigonometric
differentiation matrix (see (Trefethen, 2000, Chapter 1)) to solve DDEs on
periodic domains in a similar manner to that described in Section 3. An
advantage here is that the assumption of periodicity prevents the introduction
of discontinuities from non-smooth initial data.

######  Example 10.

Let us first consider a simple linear periodic DDE, namely

$u^{\prime\prime}(t)+\sin(t)u^{\prime}(t-\pi/\sqrt{2})+\cos(t)u(t-\pi/2)=1,\qquad
t\in[0,2\pi],$ (6.2)

with $u(t)$ $2\pi$-periodic. Code to solve this periodic problem is show in
the right hand-side of Table 2. Figure 15 shows the computed solution with
$n=32$ (left) and the convergence (centre) of the approach as the
discretisation size is increased. In this case, we do not have a closed form
expression for this solution and so take the $n=32$ solution as a proxy. The
Fourier coefficients of the $n=32$ solution are shown in the rightmost panel.
As in the case of the non-periodic problems, we observed geometric convergence
in the number of degrees of freedom, in this case down to an accuracy of
around $10^{-12}$.

⬇

1function P = barymattrig(z, t)

2 T = t(end)-2*t(1)+t(2);

3 P = cot((z-t.’)*pi/T);

4 P(:,2:2:end) = -P(:,2:2:end);

5 P = P./sum(P,2);

6 P(isnan(P)) = 1;

7end

⬇

1t = linspace(0, 2*pi, 2*n+1)’; t(end) = [];

2D1 = diffmat(2*n, 1, ’periodic’, [0 2*pi]);

3D2 = diffmat(2*n, 2, ’periodic’, [0 2*pi]);

4P1 = barymattrig(t-pi/sqrt(2),t);

5P2 = barymattrig(t-pi/2,t);

6A = D2+diag(sin(t))*P1*D1+diag(cos(t))*P2;

7u = A\ones(2*n,1);

Table 2: Left: Matlab code for constructing the trigonometric barycentric
resampling matrix for an even number of interpolation nodes t. For an odd
number of points, simply replace cot with csc. Right: Matlab code for solving
(6.2). Here we use diffmat from Chebfun to compute the first- and second-order
trigonometric differentiation matrices. TODO: This code needs updating.

Figure 15: Left: Solution using two 32-point uniform grid. Centre: Convergence
Right: Fourier coefficients.

######  Example 11.

A more interesting application is the use of these methods in the accurate
computation of limit cycles for autonomous DDE initial value problems, which
can subsequently be used in the computation of the relevant Floquet
multipliers (Luzyanina and Engelborghs, 2002). For example, consider the
following Lotka–Volterra predator-prey model with a delayed Holling type-II
response:

$\begin{array}[]{llll}x^{\prime}(t)&=&x(t)-\displaystyle\frac{x(t)^{2}}{K}-\frac{x(t)y(t-s)}{1+x(t)},\\\
y^{\prime}(t)&=&-\gamma
y(t)+\delta\displaystyle\frac{x(t)y(t-s)}{1+x(t)},\end{array},\qquad t>0.$
(6.3)

For certain parameter choices, one can show that the solution exhibits a
stable limit cycle. Here we aim to compute this limit cycle accurately using a
Fourier spectral method. The procedure is as follows:

  * (a)

For some suitable initial condition $(x_{0},y_{0})$ in the neighbourhood of
the equilibrium point
$(x_{e},y_{e})=(\tfrac{\alpha}{1-\alpha},(1-\tfrac{x_{e}}{K})(1+x_{e}))$,
solve (6.3) for some time to obtain an approximation to the limit cycle and
its unknown period, $T$.

  * (b)

Rescale (6.3) to a periodic problem on $[0,1]$ by an affine transformation of
$t$ in the unknown frequency, $1/T$.

  * (c)

Use the Fourier spectral collocation method with baryentric resampling to
implement the delay and estimates of the limit cycle and its period from (b)
as an initial guess, solve the periodic nonlinear DDE for $x$, $y$, and $T$.

  * (d)

Finally, rescale time to obtain the limit cycle solution on $[0,T]$.

A similar approach is employed by (Borgioli et al., 2020) in the context of
stability analysis of linear time-periodic delayed dynamical systems with
piecewise-smooth coefficients. In Figure we show the result of this procedure
when solving (6.3) with $K=1.4$, $\gamma=2/15$, $\delta=1$, and $s=1$, for
which the period of the limit cycle is found to be 30.83847284. Checking the
residual of the solution with Chebfun gives a forward error of around $6\times
10^{-14}$.

Figure 16: Left: Limit cycle solution computed using 129-point grid. Centre:
Limit cycle. Right: Fourier coefficients of solution.

###  6.3 Rational spectral methods

In addition to polynomial and trigonometric interpolants, the approach
described in Sections 3 and 4 can be applied to spectral methods based on
rational interpolants defined by their barycentric weights, such as those
described in (Berrut and Klein, 2014) and (Kosloff and Tal-Ezer, 1993) for
uniform or near-uniformly-spaced nodes, respectively, or those in (Hale and
Tee, 2009) for nodes clustered to resolve local features in the solution.
Furthermore, the advent of the AAA algorithm for rational approximation has
sparked a renewed interest in numerical methods based on rational functions
(Nakatsukasa et al., 2018, 2023). The AAA algorithm involves a barycentric
representation, and as such any spectral collocation method based on AAA ideas
for BVP ODEs could be naturally extended to DDEs and FDEs.

###  6.4 Functional equations

The focus of this work has been on functional differential equations, but the
area of functional equations (FEs) is in itself of interest, and there is
scope for applying the methods described in this manuscript in such a setting.

######  Example 12.

Consider the Schröder-type functional equation (Schröder, 1870)

$u(f(t))=\lambda u(t),\qquad f(t)=\lambda\sin(t).$ (6.4)

For $0<\lambda<1$, the theorem of Koenigs (1884) states that there is a unique
solution holomorphic in a neighbourhood of the origin satisfying
$u(0)=0,u^{\prime}(0)=1$. Here will will approximate this solution on a finite
interval by interpolating at scaled Chebyshev points and employing the
barycentric resampling matrix as described in previous sections.

Figure 17 shows code for implementing the above on the domain $[0,\pi]$ with
$\lambda=\tfrac{1}{2}$ and compares the result with the iterated composition
solution $u(t)=\lim_{n\rightarrow\infty}\lambda^{-n}f^{\circ n}(t)$, where
$f^{\circ n}(t):=\underbrace{f(f(f(\ldots)))}_{n\text{ times}}.$ As with the
FDE examples of previous sections, we observe geometric convergence as the
number of interpolation points is increased. The interpolation approach has
the advantage that it may be applied to more general FEs where iterated
solutions are not readily available.

⬇

1lam = 0.5; f = @(t) lam*sin(t);

2n = 30; [t, ~, v] = chebpts(n, [0,pi]);

3P = barymat(f(t), t, w); I = eye(n)

4L = P - lam*I;

5[~, ~, V] = svd(L); u = V(:,end);

6D = diffmat(t, w); u = u./(D(1,:)*u);

Figure 17: Left: MATLAB code to solve (6.4) using the barycentric resampling
appraoch. The final line normalises to ensure that $u^{\prime}(0)=0$. Centre:
Agreement of the computed solutions using interpolation and iterated
composition. Right: Geometric convergence of the interpolation solution as the
discretisation size is increased.

##  7 Conclusion

The barycentric resampling matrix was used to create a framework enabling the
solution of various delay and functional differential equations—including
pantograph, discrete, continuous, and state-dependent delay—via Chebyshev
spectral collocation. Furthermore, this framework was shown to be application
neutral-type equations, higher-order DDEs/FDEs, and DDE eigenvalue problems.
The close connection to standard Chebyshev spectral methods for ODEs enables
the ready implementation of the method to existing code bases, as was
demonstrated here in the case of Chebfun/Chebop. Numerous examples were
presented, demonstrating both rapid (typically exponential) convergence as
well as convenient implementation. Furthermore, the method was shown to be
readily extendable to methods based on other classical orthogonal polynomials,
trigonometric polynomials, or rational functions. We close with a few
additional observations and comments.

In the coded examples given, boundary bordering was used to implement boundary
conditions. If rectangular projection is used, a la (Driscoll and Hale, 2016),
then one could explicitly construct the discrete evaluation operation
operators to return values on the projected grid, thus avoiding an addition
matrix multiplication. Formulas for the explicit construction of the
rectangular differentiation matrices are also known (Xu and Hale, 2016).

Here we employed spectral collocation, which typically gives rise to dense
linear systems, requiring $O(n^{3})$ floating point operations to solve an
$n$-point discretisation using dense linear algebra. By utilising the
connection to discrete cosine transform (DCT) and the FFT, Chebyshev
differentiation matrices can be applied in $\mathcal{O}(n\log(n))$ operations,
making iterative methods like MINRES/GMRES (with suitable preconditioners)
attractive (Canuto and Quarteroni, 1985). In general, the interpolation
operator/barycentric resampling matrix requires $O(n^{2})$ operations to apply
directly, making an iterative solver less attractive in the DDE/FDE setting.
However, this can be mitigated by the use of a non-uniform DCT (Dutt and
Rokhlin, 1993; Potts and Steidl, 2003).

An alternative to spectral collocation methods are spectral methods based on
coefficient expansions (e.g., Galerkin or tau methods). The ultraspherical
spectral method is such an example, which gives rise to sparse (almost-banded)
matrices when applied to the ODEs with polynomial coefficients (Olver and
Townsend, 2013). One could introduce an interpolation operator to these
discretisations, similarly to that in Section 3. This would destroy the
sparseness of the discretisation, but other beneficial properties, such as
improved conditioning, would likely be maintained, making such an approach
attractive. We leave further investigations for future work.

## Acknowledgements

The author is grateful to André Weideman and Toby Driscoll for useful
discussions.

## References

  * Ali et al. [2009] I. Ali, H. Brunner, and T. Tang.  A spectral method for pantograph-type delay differential equations and its convergence analysis.  _J. Comput. Math._ , 27(2-3):254–265, 2009.  ISSN 0254-9409,1991-7139.
  * Baker et al. [1995] C. T. H. Baker, C. A. H. Paul, and D. R. Willé.  Issues in the numerical solution of evolutionary delay differential equations.  _Adv. Comput. Math._ , 3(3):171–196, 1995.  ISSN 1019-7168,1572-9044.  10.1007/BF02988625.  URL https://doi.org/10.1007/BF02988625.
  * Baltensperger and Trummer [2003] R. Baltensperger and M. R. Trummer.  Spectral differencing with a twist.  _SIAM Journal on Scientific Computing_ , 24(5):1465–1487, 2003.  10.1137/S1064827501388182.  URL https://doi.org/10.1137/S1064827501388182.
  * Bellen and Zennaro [2013] A. Bellen and M. Zennaro.  _Numerical methods for delay differential equations_.  Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, 2013.  ISBN 978-0-19-967137-3.  First paperback reprint of the 2003 original [MR1997488].
  * Berrut and Klein [2014] J.-P. Berrut and G. Klein.  Recent advances in linear barycentric rational interpolation.  _Journal of Computational and Applied Mathematics_ , 259:95–107, 2014.
  * Berrut and Trefethen [2004] J.-P. Berrut and L. N. Trefethen.  Barycentric Lagrange interpolation.  _SIAM Rev._ , 46(3):501–517, 2004.  ISSN 0036-1445.  10.1137/S0036144502417715.  URL https://doi.org/10.1137/S0036144502417715.
  * Birkisson and Driscoll [2012] A. Birkisson and T. A. Driscoll.  Automatic Fréchet differentiation for the numerical solution of boundary-value problems.  _ACM Trans. Math. Softw._ , 38(4), aug 2012.  ISSN 0098-3500.  10.1145/2331130.2331134.  URL https://doi.org/10.1145/2331130.2331134.
  * Borgioli et al. [2020] F. Borgioli, D. Hajdu, T. Insperger, G. Stepan, and W. Michiels.  Pseudospectral method for assessing stability robustness for linear time-periodic delayed dynamical systems.  _Internat. J. Numer. Methods Engrg._ , 121(16):3505–3528, 2020.  ISSN 0029-5981,1097-0207.  10.1002/nme.6368.  URL https://doi.org/10.1002/nme.6368.
  * Boyd [2001] J. P. Boyd.  _Chebyshev and Fourier spectral methods_.  Dover Publications, Inc., Mineola, NY, second edition, 2001.  ISBN 0-486-41183-4.
  * Brauer and Castillo-Chávez [2001] F. Brauer and C. Castillo-Chávez.  _Mathematical models in population biology and epidemiology_ , volume 40 of _Texts in Applied Mathematics_.  Springer-Verlag, New York, 2001.  ISBN 0-387-98902-1.  10.1007/978-1-4757-3516-1.  URL https://doi.org/10.1007/978-1-4757-3516-1.
  * Breda et al. [2005] D. Breda, S. Maset, and R. Vermiglio.  Pseudospectral differencing methods for characteristic roots of delay differential equations.  _SIAM Journal on Scientific Computing_ , 27(2):482–495, 2005.
  * Brunner and Hu [2005] H. Brunner and Q. Hu.  Optimal superconvergence orders of iterated collocation solutions for volterra integral equations with vanishing delays.  _SIAM Journal on Numerical Analysis_ , 43(5):1934–1949, 2005.
  * Brunner and Hu [2007] H. Brunner and Q. Hu.  Optimal superconvergence results for delay integro-differential equations of pantograph type.  _SIAM journal on numerical analysis_ , 45(3):986–1004, 2007.
  * Butcher and Bobrenkov [2011] E. A. Butcher and O. A. Bobrenkov.  On the chebyshev spectral continuous time approximation for constant and periodic delay differential equations.  _Communications in Nonlinear Science and Numerical Simulation_ , 16(3):1541–1554, 2011.
  * Canuto and Quarteroni [1985] C. Canuto and A. Quarteroni.  Preconditioned minimal residual methods for chebyshev spectral calculations.  _Journal of Computational Physics_ , 60(2):315–337, 1985.
  * Driscoll [2010] T. A. Driscoll.  Automatic spectral collocation for integral, integro-differential, and integrally reformulated differential equations.  _Journal of Computational Physics_ , 229(17):5980–5998, 2010.
  * Driscoll and Fornberg [1998] T. A. Driscoll and B. Fornberg.  A block pseudospectral method for Maxwell’s equations. I. One-dimensional case.  _J. Comput. Phys._ , 140(1):47–65, 1998.  ISSN 0021-9991.  10.1006/jcph.1998.5883.  URL https://doi.org/10.1006/jcph.1998.5883.
  * Driscoll and Hale [2016] T. A. Driscoll and N. Hale.  Rectangular spectral collocation.  _IMA J. Numer. Anal._ , 36(1):108–132, 2016.  ISSN 0272-4979.  10.1093/imanum/dru062.  URL https://doi.org/10.1093/imanum/dru062.
  * Driscoll et al. [2008] T. A. Driscoll, F. Bornemann, and L. N. Trefethen.  The chebop system for automatic solution of differential equations.  _BIT Numerical Mathematics_ , 48:701–723, 2008.  10.1007/s10543-008-0198-4.  URL http://www.springerlink.com/content/46471515240111m2/?p=864e78b520294229958456c91ab618fc&pi=5.
  * Driscoll et al. [2014] T. A. Driscoll, N. Hale, and L. N. Trefethen.  _Chebfun Guide_.  Pafnuty Publications, 2014.  URL www.chebfun.org.
  * Driver [2012] R. D. Driver.  _Ordinary and delay differential equations_ , volume 20.  Springer Science & Business Media, 2012.
  * Dutt and Rokhlin [1993] A. Dutt and V. Rokhlin.  Fast Fourier transforms for nonequispaced data.  _SIAM J. Sci. Comput._ , 14(6):1368–1393, 1993\.  ISSN 1064-8275.  10.1137/0914081.  URL https://doi.org/10.1137/0914081.
  * El-Safty and Hussien [1998] A. El-Safty and M. A. Hussien.  Chebyshev solution for stiff delay differential equations.  _Int. J. Comput. Math._ , 68(3-4):323–335, 1998\.  ISSN 0020-7160.  10.1080/00207169808804699.  URL https://doi.org/10.1080/00207169808804699.
  * Engelborghs et al. [2001] K. Engelborghs, T. Luzyanina, K. I. Hout, and D. Roose.  Collocation methods for the computation of periodic solutions of delay differential equations.  _SIAM Journal on Scientific Computing_ , 22(5):1593–1609, 2001.
  * Engelborghs et al. [2002] K. Engelborghs, T. Luzyanina, and D. Roose.  Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL.  _ACM Trans. Math. Software_ , 28(1):1–21, 2002\.  ISSN 0098-3500,1557-7295.  10.1145/513001.513002.  URL https://doi.org/10.1145/513001.513002.
  * Erneux [2009] T. Erneux.  _Applied delay differential equations_ , volume 3 of _Surveys and Tutorials in the Applied Mathematical Sciences_.  Springer, New York, 2009.  ISBN 978-0-387-74371-4.
  * Fornberg [1996] B. Fornberg.  _A practical guide to pseudospectral methods_ , volume 1 of _Cambridge Monographs on Applied and Computational Mathematics_.  Cambridge University Press, Cambridge, 1996.  ISBN 0-521-49582-2.  10.1017/CBO9780511626357.  URL https://doi.org/10.1017/CBO9780511626357.
  * Fox et al. [1971] L. Fox, D. F. Mayers, J. R. Ockendon, and A. B. Tayler.  On a functional differential equation.  _IMA Journal of Applied Mathematics_ , 8(3):271–307, 1971.
  * Gopalsamy [1992] K. Gopalsamy.  _Stability and oscillations in delay differential equations of population dynamics_ , volume 74 of _Mathematics and its Applications_.  Kluwer Academic Publishers Group, Dordrecht, 1992.  ISBN 0-7923-1594-4.  10.1007/978-94-015-7920-9.  URL https://doi.org/10.1007/978-94-015-7920-9.
  * Hale [1971] J. K. Hale.  Functional differential equations.  In _Analytic theory of differential equations_ , pages 9–22. Springer, 1971.
  * Hale [2024] N. Hale.  Software repository containing for “A spectral collocation method for functional and delay differential equations”, 2024.  URL https://github.com/nickhale/ddes/.
  * Hale and Tee [2009] N. Hale and T. W. Tee.  Conformal maps to multiply slit domains and applications.  _SIAM J. Sci. Comput._ , 31(4):3195–3215, 2009\.  ISSN 1064-8275,1095-7197.  10.1137/080738325.  URL https://doi.org/10.1137/080738325.
  * Hale et al. [2024] N. Hale, E. Thomann, and J. Weideman.  Numerical and asymptotic solutions to a nonlinear functional differential equation.  _Submitted to the IMA Journal of Applied Mathematics_ , 2024.
  * Hartung et al. [2006] F. Hartung, T. Krisztin, H.-O. Walther, and J. Wu.  Functional differential equations with state-dependent delays: theory and applications.  In _Handbook of differential equations: ordinary differential equations. Vol. III_ , Handb. Differ. Equ., pages 435–545. Elsevier/North-Holland, Amsterdam, 2006.  ISBN 978-0-444-52849-0; 0-444-52849-0.  10.1016/S1874-5725(06)80009-X.  URL https://doi.org/10.1016/S1874-5725(06)80009-X.
  * Henrici [1979] P. Henrici.  Barycentric formulas for interpolating trigonometric polynomials and their conjugates.  _Numer. Math._ , 33(2):225–234, jun 1979.  ISSN 0029-599X.  10.1007/BF01399556.  URL https://doi.org/10.1007/BF01399556.
  * Higham [2004] N. J. Higham.  The numerical stability of barycentric Lagrange interpolation.  _IMA J. Numer. Anal._ , 24(4):547–556, 2004.  ISSN 0272-4979.  10.1093/imanum/24.4.547.  URL https://doi.org/10.1093/imanum/24.4.547.
  * Ito et al. [1991] K. Ito, H. T. Tran, and A. Manitius.  A fully-discrete spectral method for delay-differential equations.  _SIAM Journal on Numerical Analysis_ , 28(4):1121–1140, 1991.  10.1137/0728060.  URL https://doi.org/10.1137/0728060.
  * Jackiewicz [1984] Z. Jackiewicz.  One-step methods of any order for neutral functional differential equations.  _SIAM J. Numer. Anal._ , 21(3):486–511, 1984\.  ISSN 0036-1429.  10.1137/0721036.  URL https://doi.org/10.1137/0721036.
  * Jafari et al. [2021] H. Jafari, M. Mahmoudi, and M. Noori Skandari.  A new numerical method to solve pantograph delay differential equations with convergence analysis.  _Advances in Difference Equations_ , 2021(1):1–12, 2021.
  * Koenigs [1884] G. Koenigs.  Recherches sur les intégrales de certaines équations fonctionnelles.  _Ann. Sci. École Norm. Sup. (3)_ , 1:3–41, 1884.  ISSN 0012-9593.  URL http://www.numdam.org/item?id=ASENS_1884_3_1__S3_0.
  * Kolmanovskii and Myshkis [1999] V. Kolmanovskii and A. Myshkis.  _Introduction to the theory and applications of functional-differential equations_ , volume 463 of _Mathematics and its Applications_.  Kluwer Academic Publishers, Dordrecht, 1999.  ISBN 0-7923-5504-0.  10.1007/978-94-017-1965-0.  URL https://doi.org/10.1007/978-94-017-1965-0.
  * Kosloff and Tal-Ezer [1993] D. Kosloff and H. Tal-Ezer.  A modified Chebyshev pseudospectral method with an $O(N^{-1})$ time step restriction.  _J. Comput. Phys._ , 104(2):457–469, 1993.  ISSN 0021-9991,1090-2716.  10.1006/jcph.1993.1044.  URL https://doi.org/10.1006/jcph.1993.1044.
  * Kuang [1993] Y. Kuang.  _Delay differential equations with applications in population dynamics_ , volume 191 of _Mathematics in Science and Engineering_.  Academic Press, Inc., Boston, MA, 1993.  ISBN 0-12-427610-5.
  * Kyrychko and Hogan [2010] Y. N. Kyrychko and S. J. Hogan.  On the use of delay equations in engineering applications.  _J. Vib. Control_ , 16(7-8):943–960, 2010.  ISSN 1077-5463,1741-2986.  10.1177/1077546309341100.  URL https://doi.org/10.1177/1077546309341100.
  * Liu and Ortiz [1989] K.-M. Liu and E. L. Ortiz.  Numerical solution of ordinary and partial functional-differential eigenvalue problems with the tau method.  _Computing_ , 41(3):205–217, 1989.
  * Luzyanina and Engelborghs [2002] T. Luzyanina and K. Engelborghs.  Computing Floquet multipliers for functional differential equations.  _Internat. J. Bifur. Chaos Appl. Sci. Engrg._ , 12(12):2977–2989, 2002.  ISSN 0218-1274,1793-6551.  10.1142/S0218127402006291.  URL https://doi.org/10.1142/S0218127402006291.
  * MathWorks [2012] MathWorks.  ddensd: Solve delay differential equations (ddes) of neutral type, 2012\.  URL https://www.mathworks.com/help/matlab/ref/ddensd.html.
  * Meng et al. [2019] T. Meng, Z. Wang, and L. Yi.  An h–p version of the chebyshev spectral collocation method for nonlinear delay differential equations.  _Numerical Methods for Partial Differential Equations_ , 35(2):664–680, 2019.
  * Michiels and Niculescu [2007] W. Michiels and S.-I. Niculescu.  _Stability and stabilization of time-delay systems: an eigenvalue-based approach_.  SIAM, 2007.
  * Nakatsukasa et al. [2018] Y. Nakatsukasa, O. Sète, and L. N. Trefethen.  The aaa algorithm for rational approximation.  _SIAM Journal on Scientific Computing_ , 40(3):A1494–A1522, 2018.
  * Nakatsukasa et al. [2023] Y. Nakatsukasa, O. Sete, and L. N. Trefethen.  The first five years of the aaa algorithm.  _arXiv preprint arXiv:2312.03565_ , 2023.
  * Neves and Feldstein [1976] K. W. Neves and A. Feldstein.  Characterization of jump discontinuities for state dependent delay differential equations.  _J. Math. Anal. Appl._ , 56(3):689–707, 1976\.  ISSN 0022-247X.  10.1016/0022-247X(76)90033-0.  URL https://doi.org/10.1016/0022-247X(76)90033-0.
  * Olver and Townsend [2013] S. Olver and A. Townsend.  A fast and well-conditioned spectral method.  _siam REVIEW_ , 55(3):462–489, 2013.
  * Potts and Steidl [2003] D. Potts and G. Steidl.  Fast summation at nonequispaced knots by NFFTs.  _SIAM J. Sci. Comput._ , 24(6):2013–2037, 2003\.  ISSN 1064-8275,1095-7197.  10.1137/S1064827502400984.  URL https://doi.org/10.1137/S1064827502400984.
  * Rihan [2021] F. A. Rihan.  _Delay differential equations and applications to biology_.  Forum for Interdisciplinary Mathematics. Springer, Singapore, 2021.  ISBN 978-981-16-0625-0; 978-981-16-0626-7.  10.1007/978-981-16-0626-7.  URL https://doi.org/10.1007/978-981-16-0626-7.
  * Schröder [1870] E. Schröder.  Ueber iterirte Functionen.  _Math. Ann._ , 3(2):296–322, 1870.  ISSN 0025-5831,1432-1807.  10.1007/BF01443992.  URL https://doi.org/10.1007/BF01443992.
  * Sedaghat et al. [2012] S. Sedaghat, Y. Ordokhani, and M. Dehghan.  Numerical solution of the delay differential equations of pantograph type via chebyshev polynomials.  _Communications in Nonlinear Science and Numerical Simulation_ , 17(12):4815–4830, 2012.
  * Shampine [2008] L. F. Shampine.  Dissipative approximations to neutral DDEs.  _Appl. Math. Comput._ , 203(2):641–648, 2008\.  ISSN 0096-3003,1873-5649.  10.1016/j.amc.2008.05.010.  URL https://doi.org/10.1016/j.amc.2008.05.010.
  * Shampine and Thompson [2009] L. F. Shampine and S. Thompson.  Numerical solution of delay differential equations.  In _Delay differential equations_ , pages 245–271. Springer, New York, 2009.  ISBN 978-0-387-85594-3.
  * Tee [2006] T. W. Tee.  _An adaptive rational spectral method for differential equations with rapidly varying solutions_.  PhD thesis, University of Oxford, 2006.
  * Trefethen [2000] L. N. Trefethen.  _Spectral methods in MATLAB_ , volume 10 of _Software, Environments, and Tools_.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000.  ISBN 0-89871-465-6.  10.1137/1.9780898719598.  URL https://doi.org/10.1137/1.9780898719598.
  * Trefethen [2019] L. N. Trefethen.  _Approximation Theory and Approximation Practice, Extended Edition_.  SIAM, 2019.
  * Trefethen et al. [2017] L. N. Trefethen, A. Birkisson, and T. A. Driscoll.  _Exploring ODEs_.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017\.  10.1137/1.9781611975161.  URL https://epubs.siam.org/doi/abs/10.1137/1.9781611975161.
  * Trif [2011] D. Trif.  Matrix based operatorial approach to differential and integral problems.  _MATLAB-A Ubiquitous Tool for the Practical Engineer_ , 2011.
  * Trif [2012] D. Trif.  Chebpack, 2012.  URL https://www.mathworks.com/matlabcentral/fileexchange/32227-chebpack.
  * Wei and Chen [2012] Y. Wei and Y. Chen.  Legendre spectral collocation methods for pantograph volterra delay-integro-differential equations.  _Journal of Scientific Computing_ , 53:672–688, 2012.
  * Weideman and Reddy [2000] J. A. C. Weideman and S. C. Reddy.  A MATLAB differentiation matrix suite.  _ACM Trans. Math. Software_ , 26(4):465–519, 2000\.  ISSN 0098-3500.  10.1145/365723.365727.  URL https://doi.org/10.1145/365723.365727.
  * Wu and Michiels [2012] Z. Wu and W. Michiels.  Reliably computing all characteristic roots of delay differential equations in a given right half plane using a spectral method.  _Journal of Computational and Applied Mathematics_ , 236(9):2499–2514, 2012.
  * Xu and Hale [2016] K. Xu and N. Hale.  Explicit construction of rectangular differentiation matrices.  _IMA Journal of Numerical Analysis_ , 36(2):618–632, 2016.  10.1093/imanum/drv013.