# A new framework for polynomial approximation to differential equations 111
Cite as: L. Brugnano, G. Frasca-Caccia, F. Iavernaro, V. Vespri. A new
framework for polynomial approximation to differential equations. Adv. Comput.
Math. 48, 76 (2022) https://doi.org/10.1007/s10444-022-09992-w

Luigi Brugnano222Università di Firenze, Italy<EMAIL_ADDRESS>orcidID: 0000-0002-6290-4107) Gianluca Frasca-Caccia333Università di Salerno,
Italy<EMAIL_ADDRESS>Felice Iavernaro444Università di Bari, Italy
<EMAIL_ADDRESS>orcidID: 0000-0002-9716-7370) Vincenzo
Vespri555Università di Firenze, Italy<EMAIL_ADDRESS>
###### Abstract

In this paper we discuss a framework for the polynomial approximation to the
solution of initial value problems for differential equations. The framework
is based on an expansion of the vector field along an orthonormal basis, and
relies on perturbation results for the considered problem. Initially devised
for the approximation of ordinary differential equations, it is here further
extended and, moreover, generalized to cope with constant delay differential
equations. Relevant classes of Runge-Kutta methods can be derived within this
framework.

Keywords: ordinary differential equations, delay differential equations,
orthogonal polynomials, local Fourier expansion, polynomial approximations,
Runge-Kutta methods.

MCS: 65L05, 65L03, 65L06, 65P10.

##  1 Introduction

In this paper, we shall deal with the definition of a framework to discuss
polynomial approximations to the solution of initial value problems for
ordinary differential equations (ODEs),

$\dot{y}(t)=f(t,y(t)),\qquad t\in[t_{0},T],\qquad
y(t_{0})=y_{0}\in\mathbb{R}^{m},$ (1)

and delay differential equations (DDEs) in the form,

$\displaystyle\dot{y}(t)$ $\displaystyle=$ $\displaystyle
f(t,y(t),y(t-\tau)),\qquad t\in[t_{0},T],\qquad y(t_{0})=y_{0},$ (2)
$\displaystyle y(t)$ $\displaystyle=$
$\displaystyle\phi(t),\qquad\qquad\qquad\qquad t\in[t_{0}-\tau,t_{0}),$

where $\tau>0$ is a constant delay and, usually, $y_{0}=\phi(t_{0})$. In the
sequel, we shall always assume that $f$ and $\phi$ are suitably regular in
their respective arguments. As is well known, the two problems are related in
many ways but, at the same time, have quite different features, which reflect
on their numerical solution. We refer, e.g., to the comprehensive monograph
[27], concerning (1), and [5] (see also [20]) for (2).

In more detail, in this paper we shall fully develop a novel framework for
deriving numerical methods for solving (1), which is then extended to cope
with (2).

The framework we are interested in relies on a local expansion of the vector
field in (1) along an orthonormal basis. Such basis will be, in the present
case, the Legendre polynomial basis $\\{P_{j}\\}_{j\geq 0}$:

$P_{j}\in\Pi_{j},\qquad\int_{0}^{1}P_{i}(x)P_{j}(x)\mathrm{d}x=\delta_{ij},\qquad
i,j=0,1,\dots,$ (3)

where, as is usual, $\Pi_{j}$ is the vector space of polynomials of degree
$j$, and $\delta_{ij}$ is the Kronecker symbol. The idea is actually not new:
early use of this approach are, for example, Hulme [29, 30], Bottasso [7], and
Betsch and Steinmann [6]; it is also at the basis of the energy-conserving
class of Runge-Kutta methods named HBVMs [12] (see also the monograph [9] and
the review paper [10]).

The approach that we shall pursue has been initially devised in [14], where
the target was problem (1), and its potentialities have been disclosed by
using HBVMs as spectral methods in time for efficiently solving highly
oscillatory problems [19] and, subsequently, Hamiltonian PDEs [11]. A
corresponding error analysis is given in [3]. Moreover, this allows deriving a
formulation of HBVMs as continuous-stage Runge-Kutta methods [1, 2].

Starting from this background, in this work we carry out a complete
perturbation analysis of problems (1) and (2), and set up a unique and
comprehensive framework to deal with the numerical solution of both problems
by exploiting the same discretization procedure. In more detail, the truncated
Fourier series may be interpreted as a projection of the differential problem
onto a finite dimensional vector space, leading to a new, numerically easy-to-
handle, differential problem. This latter may be regarded as a perturbation of
the original one, so that the perturbation analysis turns out to be crucial to
understand how the solutions of the two problems are related. At the best of
our knowledge, the perturbation results for problem (2) are new, and provide a
powerful general tool of analysis. That the same framework may cover problems
of different nature constitute, in our opinion, a specific advancement in this
field, and reveals its potentialities to deal with other classes of problems
(which will be the subject of future investigations).

With this premise, the paper is organized as follows: Section 2 concerns the
result pertaining to the ODE case; Section 3 contains the corresponding
results for the DDE case; Section 4 contains some numerical tests for the DDE
case, involving methods which are new in this setting; at last, some
concluding remarks and possible developments are reported in Section 5.

##  2 The ODE case

Without loss of generality, we shall consider problem (1) in the simpler form:

$\dot{y}(t)=f(y(t)),\qquad t\in[t_{0},T],\qquad
y(t_{0})=y_{0}\in\mathbb{R}^{m}.$ (4)

Having fixed the mesh

$t_{n}=t_{0}+nh,\qquad n=0,\dots,N,\qquad h=\frac{T-t_{0}}{N},$ (5)

we formally set, for $n=1,\dots,N$:

$\hat{\sigma}_{n}(ch):=y(t_{n-1}+ch)\equiv\hat{\sigma}(t_{n-1}+ch),\qquad
c\in[0,1],$ (6)

the restriction of the solution of problem (4) to the time interval
$[t_{n-1},t_{n}]$ (the function $\hat{\sigma}(t)\equiv y(t)$ is introduced for
notational purposes). Consequently, $\hat{\sigma}_{n}$ satisfies the
differential equation

$\dot{\hat{\sigma}}_{n}(ch)=\sum_{j\geq
0}P_{j}(c)\gamma_{j}(\hat{\sigma}_{n}),\qquad
c\in[0,1],\qquad\hat{\sigma}_{n}(0)=y(t_{n-1}),$ (7)

so that,

$\hat{\sigma}_{n}(ch)=y(t_{n-1})+h\sum_{j\geq
0}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\gamma_{j}(\hat{\sigma}_{n}),\qquad
c\in[0,1],$ (8)

and, by virtue of (3),

$y(t_{n})=y(t_{n-1})+h\gamma_{0}(\hat{\sigma}_{n})\equiv\hat{\sigma}(t_{n}),$
(9)

where, in general, for any suitably regular function
$z:[0,h]\rightarrow\mathbb{R}^{m}$,

$\gamma_{j}(z):=\int_{0}^{1}P_{j}(\zeta)f(z(\zeta h))\mathrm{d}\zeta.$ (10)

We now look for a piecewise polynomial approximation $\sigma(t)$, to the
solution of (4), such that, setting for $n=1,\dots,N$,

$\sigma_{n}(ch)\equiv\sigma(t_{n-1}+ch),\qquad c\in[0,1],$ (11)

its restriction to the time interval $[t_{n-1},t_{n}]$, $\sigma_{n}\in\Pi_{s}$
and satisfies the differential equation:

$\dot{\sigma}_{n}(ch)=\sum_{j=0}^{s-1}P_{j}(c)\gamma_{j}(\sigma_{n}),\qquad
c\in[0,1],\qquad\sigma_{n}(0)=y_{n-1},$ (12)

obtained by truncating the infinite series in (7) to a finite sum, with

$y_{n}:=\sigma_{n}(h)\equiv\sigma(t_{n}).$ (13)

Consequently, $\sigma_{n}$ can be formally written as:

$\sigma_{n}(ch)=y_{n-1}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\gamma_{j}(\sigma_{n}),\qquad
c\in[0,1],$ (14)

and (compare with (9)),

$y_{n}=y_{n-1}+h\gamma_{0}(\sigma_{n}).$ (15)

###  2.1 Preliminary results

We here provide a few preliminary results, which will be needed to derive the
main ones in the following subsections. Some of them are taken from [14] but
we also report them here, for sake of completeness.

###### Theorem 1

Let $G:[0,h]\rightarrow V$, with $V$ a vector space, admit a Taylor expansion
at 0. Then, for all $j\geq 0$:

$\int_{0}^{1}P_{j}(\zeta)G(\zeta h)\mathrm{d}\zeta=O(h^{j}).$

Proof By virtue of (3), one has:

$\displaystyle\int_{0}^{1}P_{j}(\zeta)G(\zeta h)\mathrm{d}\zeta$
$\displaystyle=$ $\displaystyle\int_{0}^{1}P_{j}(\zeta)\sum_{k\geq
0}\frac{G^{(k)}(0)}{k!}(\zeta h)^{k}\mathrm{d}\zeta=\sum_{k\geq
0}\frac{G^{(k)}(0)}{k!}h^{k}\int_{0}^{1}P_{j}(\zeta)\zeta^{k}\mathrm{d}\zeta$
$\displaystyle=$ $\displaystyle\sum_{k\geq
j}\frac{G^{(k)}(0)}{k!}h^{k}\int_{0}^{1}P_{j}(\zeta)\zeta^{k}\mathrm{d}\zeta=O(h^{j}).\,\framebox{}$

As a straightforward consequence, setting $G(\zeta h):=f(z(\zeta h))$, the
following result is proved.

###### Corollary 1

With reference to (10), one has: $\gamma_{j}(z)=O(h^{j})$.

Let us denote by

$y(t)\equiv y(t,\xi,\eta)$ (16)

the solution of the problem (compare with (4)):

$\dot{y}(t)=f(y(t)),\qquad t\in[\xi,T],\qquad y(\xi)=\eta\in\mathbb{R}^{m}.$
(17)

Hereafter, for sake of brevity, we may use either one of the two notations in
(16), depending on the needs. The following theorem contains standard
perturbation results w.r.t. all the arguments (see, e.g., [27, Section I.14]).

###### Theorem 2

With reference to the solution (16) of problem (17), one has:

$a)~{}\frac{\partial}{\partial t}y(t)=f(y(t)),\qquad
b)~{}\frac{\partial}{\partial\eta}y(t)=\Phi(t,\xi),\qquad
c)~{}\frac{\partial}{\partial\xi}y(t)=-\Phi(t,\xi)f(\eta),$

where $\Phi(t,\xi)$ is the solution of the variational problem

$\dot{\Phi}(t,\xi)=F(y(t))\Phi(t,\xi),\qquad
t\in[\xi,T],\qquad\Phi(\xi,\xi)=I\in\mathbb{R}^{m\times m},$

having set

$F(y)=\frac{\partial}{\partial y}f(y).$ (18)

From this theorem, the following result readily follows where, hereafter,
$|\cdot|$ will denote any convenient vector norm.

###### Corollary 2

With reference to (17), and assuming that $\xi\in[t_{n-1},t_{n}]$, one has:

$y(t,\xi,\eta+\delta\eta)=y(t,\xi,\eta)+\Phi(t,\xi)\delta\eta+(t-\xi)O(|\delta\eta|^{2}),\qquad
t\in[t_{n-1},t_{n}].$

###  2.2 Main results (ODE case)

With reference to (5)–(15), we are now in the position of stating the results
concerning the approximation error at the grid points,

$y(t_{n})-y_{n}\equiv\hat{\sigma}_{n}(h)-\sigma_{n}(h),\qquad n=1,\dots,N,$
(19)

and, more in general, on each subinterval $[t_{n-1},t_{n}]$:

$\delta\sigma_{n}(ch):=\hat{\sigma}_{n}(ch)-\sigma_{n}(ch)\equiv\hat{\sigma}(t_{n-1}+ch)-\sigma(t_{n-1}+ch),\qquad
c\in[0,1].$ (20)

For the first step of the approximation procedure, the following theorem holds
true, the proof being similar to that of [14, Theorem 1].

###### Theorem 3

With reference to (19) and (20), one has:

$y(t_{1})-y_{1}=O(h^{2s+1}),\qquad\|\delta\sigma_{1}\|:=\max_{c\in[0,1]}|\hat{\sigma}_{1}(ch)-\sigma_{1}(ch)|=O(h^{s+1}).$

Proof By virtue of Corollary 1 and Theorem 2, one has:

$\displaystyle\hat{\sigma}_{1}(ch)-\sigma_{1}(ch)~{}=~{}y(t_{0}+ch,t_{0},y_{0})-y(t_{0}+ch,t_{0}+ch,\sigma_{1}(ch))$
$\displaystyle=$ $\displaystyle
y(t_{0}+ch,t_{0},\sigma_{1}(0))-y(t_{0}+ch,t_{0}+ch,\sigma_{1}(ch))=\int_{ch}^{0}\frac{\mathrm{d}}{\mathrm{d}t}y(t_{0}+ch,t_{0}+t,\sigma_{1}(t))\mathrm{d}t$
$\displaystyle=$
$\displaystyle\int_{ch}^{0}\left[\left.\frac{\partial}{\partial\xi}y(t_{0}+ch,\xi,\sigma_{1}(t))\right|_{\xi=t_{0}+t}+\left.\frac{\partial}{\partial\eta}y(t_{0}+ch,t_{0}+t,\eta)\right|_{\eta=\sigma_{1}(t)}\dot{\sigma}_{1}(t)\right]\mathrm{d}t$
$\displaystyle=$
$\displaystyle\int_{0}^{ch}\Phi(t_{0}+ch,t_{0}+t)\left[f(\sigma_{1}(t))-\dot{\sigma}_{1}(t)\right]\mathrm{d}t$
$\displaystyle=$ $\displaystyle h\int_{0}^{c}\Phi(t_{0}+ch,t_{0}+\zeta
h)\left[f(\sigma_{1}(\zeta h))-\dot{\sigma}_{1}(\zeta
h)\right]\mathrm{d}\zeta$ $\displaystyle=$ $\displaystyle
h\int_{0}^{c}\Phi(t_{0}+ch,t_{0}+\zeta h)\left[\sum_{j\geq
0}P_{j}(\zeta)\gamma_{j}(\sigma_{1})-\sum_{j=0}^{s-1}P_{j}(\zeta)\gamma_{j}(\sigma_{1})\right]\mathrm{d}\zeta$
$\displaystyle=$ $\displaystyle h\sum_{j\geq
s}\left[\int_{0}^{c}P_{j}(\zeta)\Phi(t_{0}+ch,t_{0}+\zeta
h)\mathrm{d}\zeta\right]\underbrace{\gamma_{j}(\sigma_{1})}_{=\,O(h^{j})}.$

Consequently, the second part of the statement follows for $c\in(0,1)$,
whereas, when $c=1$ one deduces, by virtue of Theorem Theorem 1:

$\displaystyle y(t_{1})-y_{1}$ $\displaystyle\equiv$
$\displaystyle\hat{\sigma}_{1}(h)-\sigma_{1}(h)\,=\,h\sum_{j\geq
s}\overbrace{\underbrace{\left[\int_{0}^{1}P_{j}(\zeta)\underbrace{\Phi(t_{1},t_{0}+\zeta
h)}_{=:\,G(\zeta
h)}\mathrm{d}\zeta\right]}_{=O(h^{j})}\gamma_{j}(\sigma_{1})}^{=\,O(h^{2j})}~{}=~{}O(h^{2s+1}).\,\framebox{}$

For the remaining steps, the following result holds true.

###### Theorem 4

With reference to (19) and (20), for $n=1,\dots,N$ one has:

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+O(h^{2s+1}),\quad\|\delta\sigma_{n}\|:=\max_{c\in[0,1]}|\delta\sigma_{n}(ch)|=O(h^{s+1}).$

Proof By induction on $n$. For $n=1$ the statement follows from the previous
Theorem 3. Assuming it true for $n-1$, for $n$ one has:

$\displaystyle\hat{\sigma}_{n}(ch)-\sigma_{n}(ch)$ $\displaystyle=$
$\displaystyle
y(t_{n-1}+ch,t_{n-1},\hat{\sigma}_{n}(0))-y(t_{n-1}+ch,t_{n-1}+ch,\sigma_{n}(ch))$
$\displaystyle=$
$\displaystyle\underbrace{y(t_{n-1}+ch,t_{n-1},\hat{\sigma}_{n}(0))-y(t_{n-1}+ch,t_{n-1},\sigma_{n}(0))}_{=:\,E_{n,1}(ch)}~{}+$
$\displaystyle\underbrace{y(t_{n-1}+ch,t_{n-1},\sigma_{n}(0))-y(t_{n-1}+ch,t_{n-1}+ch,\sigma_{n}(ch))}_{=:\,E_{n,2}(ch)}.$

By using similar arguments as those used in the proof of Theorem 3, one
deduces that

$E_{n,2}(ch)=\left\\{\begin{array}[]{cc}O(h^{s+1}),&c\in(0,1),\\\\[5.69054pt]
O(h^{2s+1}),&c=1.\end{array}\right.$

Moreover, considering that, by the induction hypothesis,

$\delta\sigma_{n}(0)=\hat{\sigma}_{n}(0)-\sigma_{n}(0)=y(t_{n-1})-y_{n-1}=(n-1)\,O(h^{2s+1}),$

from Corollary 2, one has:

$\displaystyle E_{n,1}(ch)$ $\displaystyle=$
$\displaystyle\underbrace{\Phi(t_{n-1}+ch,t_{n-1})}_{=\,I+O(ch)}\delta\sigma_{n}(0)+ch\,O(|\delta\sigma_{n}(0)|^{2})$
$\displaystyle=$ $\displaystyle y(t_{n-1})-y_{n-1}+(n-1)O(ch^{2s+2}).$

Consequently, for $c=1$ one obtains the first part of the statement, whereas
the second part follows by taking $c\in(0,1)$.

###### Remark 1

We observe that the two equivalent equations (see (10), (12), and (14)):

$\dot{\sigma}_{n}(ch)=\sum_{j=0}^{s-1}P_{j}(c)\int_{0}^{1}P_{j}(\zeta)f(\sigma_{n}(\zeta
h))\mathrm{d}\zeta,\qquad c\in[0,1],\qquad\sigma_{n}(0)=y_{n-1},$

and

$\sigma_{n}(ch)=y_{n-1}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{j}(x)\mathrm{d}x\int_{0}^{1}P_{j}(\zeta)f(\sigma_{n}(\zeta
h))\mathrm{d}\zeta,\qquad c\in[0,1],$ (21)

define a so called HBVM$(\infty,s)$ method on the interval $[t_{n-1},t_{n}]$
(equation (21) is named Master Functional Equation in [12]. See also [9, 10]).
Consequently, such method defines an order $2s$ approximation procedure for
all $s\geq 1$, which can be also recast as a continuous-stage Runge-Kutta
method [1]. In particular, the case $s=1$ corresponds to the so called AVF
method [40]; the case $s\geq 1$ has been also considered in [26].

An interesting question concerns the difference between the Fourier
coefficients of the solution (8)–(10) and those of the polynomial
approximation (14) on the interval $[t_{n-1},t_{n}]$. The next result
clarifies the issue.

###### Theorem 5

With reference to (8), (10), and (14), for all $n=1,\dots,N$ one has:

$\delta\gamma_{j}^{n}:=\gamma_{j}(\hat{\sigma}_{n})-\gamma_{j}(\sigma_{n})=O(h^{2s-j}),\qquad
j=0,\dots,s-1.$

Proof First of all, from (3), (8), (14), and Theorem 4 we know that:

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+h[\gamma_{0}(\hat{\sigma}_{n})-\gamma_{0}(\sigma_{n})]=y(t_{n-1})-y_{n-1}+O(h^{2s+1}).$

Consequently, from the last equality one derives:

$\delta\gamma_{0}^{n}=\gamma_{0}(\hat{\sigma}_{n})-\gamma_{0}(\sigma_{n})=O(h^{2s}).$

Further, by taking into account (18) and (20), one obtains:

$\displaystyle O(h^{2s})$ $\displaystyle=$
$\displaystyle\gamma_{0}(\hat{\sigma}_{n})-\gamma_{0}(\sigma_{n})~{}=~{}\int_{0}^{1}\left[f(\hat{\sigma}_{n}(\zeta
h))-f(\sigma_{n}(\zeta h))\right]\mathrm{d}\zeta$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\underbrace{\int_{0}^{1}F(\sigma_{n}(\zeta
h)+c\,\delta\sigma_{n}(\zeta h))\mathrm{d}c}_{=:\,G(\zeta
h)}\delta\sigma_{n}(\zeta h)\mathrm{d}\zeta~{}=~{}\int_{0}^{1}G(\zeta
h)\delta\sigma_{n}(\zeta h)\mathrm{d}\zeta.$

Now, considering that $P_{0}(x)\equiv 1$ and, for all $\zeta\in[0,1]$,

$\displaystyle\int_{0}^{\zeta}P_{j}(x)\mathrm{d}x$ $\displaystyle=$
$\displaystyle\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta),\qquad j\geq 1,$
$\displaystyle with\qquad\xi_{j}=\left(2\sqrt{4j^{2}-1}\right)^{-1},$

one has:

$\displaystyle\delta\sigma_{n}(\zeta h)$ $\displaystyle=$
$\displaystyle\hat{\sigma}_{n}(\zeta h)-\sigma_{n}(\zeta h)$ $\displaystyle=$
$\displaystyle
y(t_{n-1})-y_{n-1}\,+\,h\sum_{j=0}^{s-1}\int_{0}^{\zeta}P_{j}(x)\mathrm{d}x\,\delta\gamma_{j}^{n}\,+\,h\sum_{j\geq
s}\int_{0}^{\zeta}P_{j}(x)\mathrm{d}x\,\gamma_{j}(\hat{\sigma}_{n})$
$\displaystyle=$ $\displaystyle y(t_{n-1})-y_{n-1}\,+\,\zeta
h\delta\gamma_{0}^{n}+h\sum_{j=1}^{s-1}\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\delta\gamma_{j}^{n}\,+$
$\displaystyle h\sum_{j\geq
s}\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\gamma_{j}(\hat{\sigma}_{n}).$

Consequently, from Theorem 1 and Corollary 1, one obtains:

$\displaystyle O(h^{2s})~{}=~{}\int_{0}^{1}G(\zeta h)\delta\sigma_{n}(\zeta
h)\,\mathrm{d}\zeta~{}=~{}\underbrace{\int_{0}^{1}G(\zeta
h)\mathrm{d}\zeta}_{=\,O(1)}\underbrace{[y(t_{n-1})-y_{n-1}]}_{=\,(n-1)\,O(h^{2s+1})}$
$\displaystyle+~{}h\underbrace{\int_{0}^{1}G(\zeta
h)\zeta\mathrm{d}\zeta}_{=\,O(1)}\,\underbrace{\delta\gamma_{0}^{n}}_{=\,O(h^{2s})}~{}+~{}h\sum_{j=1}^{s-1}\int_{0}^{1}G(\zeta
h)\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\mathrm{d}\zeta\,\delta\gamma_{j}^{n}$
$\displaystyle+\underbrace{h\sum_{j\geq s}\overbrace{\int_{0}^{1}G(\zeta
h)\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\mathrm{d}\zeta}^{=\,O(h^{j-1})}\,\overbrace{\gamma_{j}(\hat{\sigma}_{n})}^{=\,O(h^{j})}}_{=\,O(h^{2s})},$

from which,

$O(h^{2s})\,=\,h\sum_{j=1}^{s-1}\overbrace{\int_{0}^{1}G(\zeta
h)\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\mathrm{d}\zeta}^{=\,O(h^{j-1})}\delta\gamma_{j}^{n}$

follows and, therefore, one concludes that $\delta\gamma_{j}^{n}=O(h^{2s-j})$,
$j=0\dots,s-1$.

###  2.3 Conservative/dissipative problems

An interesting case [25, 9, 36] is that when problem (4) is in the form

$\dot{y}(t)=S\nabla H(y(t)),\qquad t\in[t_{0},T],\qquad
y(t_{0})=y_{0}\in\mathbb{R}^{m},$ (23)

with $S\in\mathbb{R}^{m\times m}$ either a skew-symmetric matrix,
$S^{\top}=-S$, or a negative semidefinite matrix, $S\leq 0$, whereas $\nabla
H$ is the gradient of a scalar function usually called the Hamiltonian. As is
clear:

  * •

when $S^{\top}=-S$:

$\frac{\mathrm{d}}{\mathrm{d}t}H(y(t))=\nabla H(y(t))^{\top}\dot{y}(t)=\nabla
H(y(t))^{\top}S\nabla H(y(t))=0,$

so that $H$ is a conserved quantity, and the problem is said to be
conservative;

  * •

when $S\leq 0$:

$\frac{\mathrm{d}}{\mathrm{d}t}H(y(t))=\nabla H(y(t))^{\top}\dot{y}(t)=\nabla
H(y(t))^{\top}S\nabla H(y(t))\leq 0,$

and the problem is said to be dissipative.

The next result shows that this behavior is preserved by the approximations
(12)–(15), upon observing that in this case (10) can be conveniently rewritten
as

$\gamma_{j}(z)=S\int_{0}^{1}P_{j}(\zeta)\nabla H(z(\zeta
h))\mathrm{d}\zeta=:S\beta_{j}(z).$ (24)

###### Theorem 6

With reference to (12)–(15) applied for approximating problem (23), for all
$n=1,\dots,N$ one has:

  * •

$H(y_{n})=H(y_{n-1})$, when $S^{\top}=-S$;

  * •

$H(y_{n})\leq H(y_{n-1})$, when $S\leq 0$.

Proof In fact, by considering that $y_{n}=\sigma_{n}(h)$,
$y_{n-1}=\sigma_{n}(0)$, and taking into account (24), one has:

$\displaystyle
H(y_{n})-H(y_{n-1})~{}=~{}H(\sigma_{n}(h))-H(\sigma_{n}(0))=\int_{0}^{h}\frac{\mathrm{d}}{\mathrm{d}t}H(\sigma_{n}(t))\mathrm{d}t$
$\displaystyle=$ $\displaystyle h\int_{0}^{1}\nabla
H(\sigma_{n}(ch))^{\top}\dot{\sigma}_{n}(ch)\mathrm{d}c=h\int_{0}^{1}\nabla
H(\sigma_{n}(ch))^{\top}\sum_{j=0}^{s-1}P_{j}(c)\gamma_{j}(\sigma_{n})\mathrm{d}c$
$\displaystyle=$ $\displaystyle
h\sum_{j=0}^{s-1}\underbrace{\left[\int_{0}^{1}\nabla
H(\sigma_{n}(ch))P_{j}(c)\mathrm{d}c\right]^{\top}}_{=\,\beta_{j}(\sigma_{n})^{\top}}S\beta_{j}(\sigma_{n})=h\sum_{j=0}^{s-1}\beta_{j}(\sigma_{n})^{\top}S\beta_{j}(\sigma_{n})~{}=:~{}\Delta
H_{n}.$

Consequently, if $S^{\top}=-S$, then $\Delta H_{n}=0$, whereas $\Delta
H_{n}\leq 0$, when $S\leq 0$.

###### Remark 2

According to Remark 1, one then obtains that HBVM$(\infty,s)$ methods can
preserve the conservative/dissipative feature of problem (23).

###  2.4 Discretization and Runge-Kutta formulation

Quoting Dahlquist and Björk [22, p. 521] “as is well known, even many
relatively simple integrals cannot be expressed in finite terms of elementary
functions, and thus must be evaluated by numerical methods.” In this context,
this quite obvious statement means that the approximation procedure defined by
(12) and (10) does not yet provide a “true” numerical method. In fact, the
integrals defining the Fourier coefficients,

$\gamma_{j}(\sigma_{n})=\int_{0}^{1}P_{j}(\zeta)f(\sigma_{n}(\zeta
h))\mathrm{d}\zeta,\qquad j=0,\dots,s-1,$ (25)

need to be numerically approximated by using a quadrature formula. Since we
are dealing with a polynomial approximation, it is quite natural to do this by
using an interpolatory quadrature with abscissae and weights $(c_{i},b_{i})$,
$i=1,\dots,k$ (we shall always assume $k$ distinct abscissae):

$\gamma_{j}(\sigma_{n})=\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(\sigma_{n}(c_{i}h))+\Delta_{j}(h),$
(26)

where $\Delta_{j}(h)$ is the quadrature error. The following straightforward
result holds true.

###### Theorem 7

If the quadrature $(c_{i},b_{i})$, $i=1,\dots,k$ has order $q$, i.e., it is
exact for polynomial integrands of degree $q-1$, then

$\Delta_{j}(h)=O(h^{q-j}),\qquad j=0,\dots,s-1.$

###### Remark 3

As is well known, since the quadrature (26) is based at $k$ (distinct)
abscissae, $q\in\\{k,\dots,2k\\}$: the lower limit is obtained by a generic
choice of the abscissae, whereas the upper one is achieved by placing them at
the zeros of $P_{k}(c)$.

When using a quadrature, clearly the Fourier coefficients (25) may be not
exactly evaluated anymore. This implies that we are actually computing a
possibly different piecewise polynomial approximation $u(t)$ such that
(compare with (10)–(15)), for all $n=1,\dots,N$:

$\displaystyle u_{n}(ch)$ $\displaystyle\equiv$ $\displaystyle
u(t_{n-1}+ch),\qquad c\in[0,1],$ (27) $\displaystyle\dot{u}_{n}(ch)$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{s-1}P_{j}(c)\hat{\gamma}_{j}(u_{n}),\qquad
c\in[0,1],\qquad u_{n}(0)=y_{n-1},$ (28)

with (see (26))

$\displaystyle\hat{\gamma}_{j}(u_{n})$ $\displaystyle:=$
$\displaystyle\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(u_{n}(c_{i}h))\equiv\gamma_{j}(u_{n})-\Delta_{j}(h),\quad
j=0,\dots,s-1,$ (29) $\displaystyle u_{n}(ch)$ $\displaystyle=$ $\displaystyle
y_{n-1}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\hat{\gamma}_{j}(u_{n}),\qquad
c\in[0,1],$ (30) $\displaystyle u_{n}(h)$ $\displaystyle=:$ $\displaystyle
y_{n}\equiv y_{n-1}+h\hat{\gamma}_{0}(u_{n}).$ (31)

Actually, (29)–(31) define the $n$th integration step, by using a timestep
$h$, performed with the $k$ stage Runge-Kutta method having stages:

$Y_{i}^{n}:=u_{n}(c_{i}h),\qquad i=1,\dots,k.$ (32)

In fact, evaluating (30) at the abscissae $c_{1},\dots,c_{k}$, and
substituting in it the $s$ approximate Fourier coefficients (29), one obtains,
after rearranging terms,

$\displaystyle Y_{i}^{n}$ $\displaystyle=$ $\displaystyle
y_{n-1}+h\sum_{\ell=1}^{k}\underbrace{b_{\ell}\sum_{j=0}^{s-1}\int_{0}^{c_{i}}P_{j}(x)\mathrm{d}x\,P_{j}(c_{\ell})}_{=:\,a_{i\ell}}f(Y_{\ell}^{n}),\qquad
i=1,\dots,n,$ (33) $\displaystyle y_{n}$ $\displaystyle=$ $\displaystyle
y_{n-1}+h\sum_{i=1}^{k}b_{i}f(Y_{i}^{n}).$ (34)

In other words, we have derived the $k$-stage Runge-Kutta method with
abscissae and weights $(c_{i},b_{i})$, $i=1,\dots,k$, and Butcher matrix
$A=\begin{pmatrix}a_{i\ell}\end{pmatrix}\in\mathbb{R}^{k\times k}$. Next
theorem puts the Butcher tableau in a more compact form [9].

###### Theorem 8

The Butcher tableau of the Runge-Kutta method (33)–(34) is given by

$\begin{array}[]{c|c}\bm{c}&{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\\\
\hline\cr&\\\\[-9.95845pt] &\bm{b}^{\top}\end{array}$ (35)

where

$\bm{b}=\begin{pmatrix}b_{1}&\dots&b_{k}\end{pmatrix}^{\top},\qquad\bm{c}=\begin{pmatrix}c_{1}&\dots&c_{k}\end{pmatrix}^{\top},\qquad\Omega=\mathrm{diag}(\bm{b}),$

and

${\cal P}_{s}=\begin{pmatrix}P_{0}(c_{1})&\dots&P_{s-1}(c_{1})\\\
\vdots&&\vdots\\\ P_{0}(c_{k})&\dots&P_{s-1}(c_{k})\end{pmatrix},\qquad{\cal
I}_{s}=\begin{pmatrix}\int_{0}^{c_{1}}P_{0}(x)\mathrm{d}x&\dots&\int_{0}^{c_{1}}P_{s-1}(x)\mathrm{d}x\\\
\vdots&&\vdots\\\
\int_{0}^{c_{k}}P_{0}(x)\mathrm{d}x&\dots&\int_{0}^{c_{k}}P_{s-1}(x)\mathrm{d}x\end{pmatrix}.$

It is possible to derive an alternative formulation of the Runge-Kutta method
(35). In fact, using the relation (2.2) between the integrals of the Legendre
polynomials and the polynomials themselves, and considering that

$\int_{0}^{c}P_{0}(x)\mathrm{d}x=\xi_{1}P_{1}(c)+\xi_{0}P_{0}(c),\qquad\xi_{0}=\frac{1}{2},$

one has that ${\cal I}_{s}={\cal P}_{s+1}\hat{X}_{s}$, where

${\cal P}_{s+1}=\begin{pmatrix}P_{0}(c_{1})&\dots&P_{s}(c_{1})\\\
\vdots&&\vdots\\\
P_{0}(c_{k})&\dots&P_{s}(c_{k})\end{pmatrix},\qquad\hat{X}_{s}=\begin{pmatrix}\xi_{0}&-\xi_{1}\\\
\xi_{1}&0&\ddots\\\ &\ddots&\ddots&-\xi_{s-1}\\\ &&\xi_{s-1}&0\\\
\hline\cr&&&\xi_{s}\end{pmatrix}=:\begin{pmatrix}X_{s}\\\ \hline\cr
0\,\dots\,0\,\xi_{s}\end{pmatrix}.$

Consequently, the Butcher tableau (35) can be rewritten as

$\begin{array}[]{c|c}\bm{c}&{\cal P}_{s+1}\hat{X}_{s}{\cal
P}_{s}^{\top}\Omega\\\ \hline\cr&\\\\[-9.95845pt]
&\bm{b}^{\top}\end{array}~{}.$

When the quadrature (26) has order $q\geq 2s$, it is quite straightforward to
prove that

${\cal P}_{s}^{\top}\Omega{\cal P}_{s}=I_{s},\qquad{\cal
P}_{s}^{\top}\Omega{\cal
P}_{s+1}=[I_{s}~{}\bm{0}]\in\mathbb{R}^{s\times(s+1)},$

where in general, hereafter, $I_{r}\in\mathbb{R}^{r\times r}$ is the identity
matrix (when the dimension of the identity matrix is not explicitly indicated,
it will be easily deducible from the context). Consequently,

${\cal P}_{s}^{\top}\Omega\left[{\cal P}_{s+1}\hat{X}_{s}{\cal
P}_{s}^{\top}\Omega\right]{\cal P}_{s}=X_{s},$

which can be regarded as a generalization of the $W$-transformation in [28,
Theorem 5.6, p. 79]. In addition to this, when $q\geq 2s$ also the following
results hold true (for sake of brevity, we do not discuss the case $q<2s$,
since it has no practical interest).

###### Theorem 9

With reference to (4)–(10) and (27)–(31), and assuming that the quadrature
formula (26) has order $q\geq 2s$, one has:

$y(t_{1})-y_{1}=O(h^{2s+1}),\qquad\max_{c\in[0,1]}|\hat{\sigma}_{1}(ch)-u_{1}(ch)|=O(h^{s+1}).$

###### Theorem 10

With reference to (4)–(10) and (27)–(31), and assuming that the quadrature
formula (26) has order $q\geq 2s$, for $n=1,\dots,N$ one has:

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+O(h^{2s+1}),\qquad\max_{c\in[0,1]}|\hat{\sigma}_{n}(ch)-u_{n}(ch)|=O(h^{s+1}).$

###### Theorem 11

With reference to (8)–(10) and (29)–(31), and assuming that the quadrature
formula (26) has order $q\geq 2s$, for all $n=1,\dots,N$ one has:

$\delta\hat{\gamma}_{j}^{n}:=\gamma_{j}(\hat{\sigma}_{n})-\hat{\gamma}_{j}(u_{n})=O(h^{2s-j}),\qquad
j=0,\dots,s-1.$

Concerning the case of conservative/dissipative problems in the form (23), the
result of Theorem 6 modifies as follows.

###### Theorem 12

With reference to (27)–(31) applied for approximating problem (23), and
assuming that the quadrature formula (26) has order $q\geq 2s$, for all
$n=1,\dots,N$ one has:

  * •

if $H$ is a polynomial of degree not larger than $q/s$, then the result of
Theorem 6 continues to hold;

  * •

differently,

    * –

$H(y_{n})=H(y_{n-1})+O(h^{q+1})$, when $S^{\top}=-S$,

    * –

$H(y_{n})\leq H(y_{n-1})+O(h^{q+1})$, when $S\leq 0$.

We here provide only the proof of Theorem 9 (see also [14, Theorem 4]), since
those of Theorem 10, Theorem 11, and Theorem 12 can be similarly obtained by
slightly adapting the corresponding proofs of Theorem 4, Theorem 5, and
Theorem 6, respectively.

Proof (of Theorem 9) By taking into account the result of Theorem 7, one has:

$\displaystyle\hat{\sigma}_{1}(ch)-u_{1}(ch)=y(t_{0}+ch,t_{0},y_{0})-y(t_{0}+ch,t_{0}+ch,u_{1}(ch))$
$\displaystyle=$ $\displaystyle
y(t_{0}+ch,t_{0},u_{1}(0))-y(t_{0}+ch,t_{0}+ch,u_{1}(ch))=\int_{ch}^{0}\frac{\mathrm{d}}{\mathrm{d}t}y(t_{0}+ch,t_{0}+t,u_{1}(t))\mathrm{d}t$
$\displaystyle=$
$\displaystyle\int_{ch}^{0}\left[\left.\frac{\partial}{\partial\xi}y(t_{0}+ch,\xi,u_{1}(t))\right|_{\xi=t_{0}+t}+\left.\frac{\partial}{\partial\eta}y(t_{0}+ch,t_{0}+t,\eta)\right|_{\eta=u_{1}(t)}\dot{u}_{1}(t)\right]\mathrm{d}t$
$\displaystyle=$
$\displaystyle\int_{0}^{ch}\Phi(t_{0}+ch,t_{0}+t)\left[f(u_{1}(t))-\dot{u}_{1}(t)\right]\mathrm{d}t$
$\displaystyle=$ $\displaystyle h\int_{0}^{c}\Phi(t_{0}+ch,t_{0}+\zeta
h)\left[f(u_{1}(\zeta h))-\dot{u}_{1}(\zeta h)\right]\mathrm{d}\zeta$
$\displaystyle=$ $\displaystyle h\int_{0}^{c}\Phi(t_{0}+ch,t_{0}+\zeta
h)\left[\sum_{j\geq
0}P_{j}(\zeta)\gamma_{j}(u_{1})-\sum_{j=0}^{s-1}P_{j}(\zeta)\hat{\gamma}_{j}(u_{1})\right]\mathrm{d}\zeta$
$\displaystyle=$ $\displaystyle h\int_{0}^{c}\Phi(t_{0}+ch,t_{0}+\zeta
h)\left[\sum_{j\geq
0}P_{j}(\zeta)\gamma_{j}(u_{1})-\sum_{j=0}^{s-1}P_{j}(\zeta)\left(\gamma_{j}(u_{1})-\Delta_{j}(h)\right)\right]\mathrm{d}\zeta$
$\displaystyle=$ $\displaystyle h\sum_{j\geq
s}\left[\int_{0}^{c}P_{j}(\zeta)\Phi(t_{0}+ch,t_{0}+\zeta
h)\mathrm{d}\zeta\right]\underbrace{\gamma_{j}(u_{1})}_{=\,O(h^{j})}$
$\displaystyle\,+\,h\sum_{j=0}^{s-1}\left[\int_{0}^{c}P_{j}(\zeta)\Phi(t_{0}+ch,t_{0}+\zeta
h)\mathrm{d}\zeta\right]\underbrace{\Delta_{j}(h)}_{=\,O(h^{q-j})}.$

Consequently, the second part of the statement follows by considering that,
for $c\in(0,1)$, this quantity is

$O(h^{s+1})+O(h^{q-s+2})=O(h^{s+1}),$

since $q\geq 2s$, whereas, when $c=1$ one deduces, by virtue of Theorem 1, and
considering that $t_{1}=t_{0}+h$:

$\displaystyle y(t_{1})-y_{1}$ $\displaystyle\equiv$
$\displaystyle\hat{\sigma}_{1}(h)-u_{1}(h)\,=\,h\sum_{j\geq
s}\overbrace{\underbrace{\left[\int_{0}^{1}P_{j}(\zeta)\underbrace{\Phi(t_{1},t_{0}+\zeta
h)}_{=:\,G(\zeta
h)}\mathrm{d}\zeta\right]}_{=\,O(h^{j})}\gamma_{j}(u_{1})}^{=\,O(h^{2j})}$
$\displaystyle\,+\,h\sum_{j=0}^{s-1}\overbrace{\underbrace{\left[\int_{0}^{1}P_{j}(\zeta)\underbrace{\Phi(t_{1},t_{0}+\zeta
h)}_{=:\,G(\zeta
h)}\mathrm{d}\zeta\right]}_{=\,O(h^{j})}\Delta_{j}(h)}^{=\,O(h^{q})}\,=\,O(h^{2s+1})+O(h^{q+1})~{}=~{}O(h^{2s+1}).\,\framebox{}$

###### Remark 4

When the $k$ abscissae are placed at the zeros of $P_{k}(c)$, and $k\geq s$,
one obtains a HBVM$(k,s)$ method, whose order is $2s$ [12, 9, 10]. It is worth
mentioning that the HBVM$(s,s)$ method is nothing but the $s$-stage Gauss-
Legendre collocation method. Moreover, the HBVM$(k,1)$ methods correspond to
the second-order Runge-Kutta methods described in [21]. Different choices of
the quadrature have been also considered in [31, 32, 33, 15].

###  2.5 Solving the discrete problems

Sometimes, the number of stages $k$ of the Runge-Kutta method (35) can be much
larger than the degree $s$ of the underlying polynomial approximation
(29)–(31). This is the case, for example, of HBVM$(k,s)$ methods when used as
energy-conserving methods [12, 9, 10] (see also Theorem 12 in Section 2.3). In
such a case, it is clear that the usual implementation of the Runge-Kutta
method leads to the solution of a discrete problem having (block) dimension
$k$. Nevertheless, for sake of completeness we now recall how the discrete
problem to be solved can be actually recast so as to have (block) dimension
$s$, independently of $k$ [13]. This clearly allows for relatively large
values of $k$, thus making possible the use of an arbitrarily high-order
quadrature (26). Let us then consider the first integration step of the method
for solving (4) with timestep $h$, (thus, we can skip the index $n$ of the
step). Setting
$\bm{1}=\begin{pmatrix}1\,,\,\dots\,,\,1\end{pmatrix}^{\top}\in\mathbb{R}^{k}$,
and $Y$ the stage vector of (block) dimension $k$, one obtains that the stage
equation for (35) is given by:

$Y=\bm{1}\otimes y_{0}+h{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\otimes
I_{m}\,f(Y),$ (36)

with an obvious meaning of $f(Y)$. However, we observe that [13]

${\cal P}_{s}^{\top}\Omega\otimes
I_{m}\,f(Y)=:\hat{\bm{\gamma}}\equiv\begin{pmatrix}\hat{\gamma}_{0}(u_{1})\\\
\vdots\\\ \hat{\gamma}_{s-1}(u_{1})\end{pmatrix},$

i.e., the (block) vector with the $s$ coefficients of the polynomial
approximation $u_{1}(ch)$ (see (29)–(30)). Consequently, (36) can be rewritten
as

$Y=\bm{1}\otimes y_{0}+h{\cal I}_{s}\otimes I_{m}\,\hat{\bm{\gamma}}.$

By combining the last two equations one eventually obtains:

$\hat{\bm{\gamma}}={\cal P}_{s}^{\top}\Omega\otimes
I_{m}\,f\left(\bm{1}\otimes y_{0}+h{\cal I}_{s}\otimes
I_{m}\,\hat{\bm{\gamma}}\right),$ (37)

which is a discrete problem, equivalent to (36), having (block) dimension $s$,
independently of $k$. Once this equation has been solved, the new
approximation is derived, according to (31), as

$y_{1}=y_{0}+h\hat{\gamma}_{0}(u_{1}).$

It is also worth mentioning that very effective nonlinear iterations have been
devised for solving (37) [13, 8, 9] (the most effective being that derived
from the so called blended iteration introduced in [17], see also [18]).

##  3 The DDE case

As for the ODE case, also for DDEs we shall consider, without loss of
generality, the simpler problem

$\displaystyle\dot{y}(t)$ $\displaystyle=$ $\displaystyle
f(y(t),y(t-\tau)),\qquad t\in[t_{0},T],\qquad y(t_{0})=y_{0},$ (38)
$\displaystyle y(t)$ $\displaystyle=$
$\displaystyle\phi(t),\qquad\qquad\qquad\quad\,t\in[t_{0}-\tau,t_{0}),$

in place of (2) where, usually, $y_{0}=\phi(t_{0})$. Moreover, we shall
suppose that both the timestep $h$ defining the discrete mesh (5) and the
width of the integration interval, $T-t_{0}$, are commensurable with the
delay:

$\tau=\nu h,\qquad T-t_{0}=K\tau,\qquad K,\nu\in\mathbb{N},$ (39)

so that the discrete mesh is now given by:

$t_{n}=t_{0}+nh,\qquad n=-\nu,\dots,N\equiv K\nu.$ (40)

On one hand, similarly as done in the ODE case, let us denote, for notational
purposes, by $\hat{\sigma}(t)\equiv y(t)$ the solution of (38), and

$\hat{\sigma}_{n}(ch):=\hat{\sigma}(t_{n-1}+ch),\qquad c\in[0,1],\qquad
n=1-\nu,\dots,N,$ (41)

its restriction to the time interval $[t_{n-1},t_{n}]$. Consequently,

$\hat{\sigma}_{n}(ch)\equiv\phi(t_{n-1}+ch),\qquad c\in[0,1],\qquad
n=1-\nu,\dots,0,$ (42)

whereas, for $n=1,\dots,N$, one has (compare with (7)–(10)):

$\dot{\hat{\sigma}}_{n}(ch)=\sum_{j\geq
0}P_{j}(c)\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu}),\qquad
c\in[0,1],\qquad\hat{\sigma}_{n}(0)=y(t_{n-1}),$ (43)

so that,

$\hat{\sigma}_{n}(ch)=y(t_{n-1})+h\sum_{j\geq
0}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu}),\qquad
c\in[0,1],$ (44)

and

$y(t_{n})=y(t_{n-1})+h\gamma_{0}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu})\equiv\hat{\sigma}(t_{n}),$
(45)

where, in general, for any suitably regular functions
$z,w:[0,h]\rightarrow\mathbb{R}^{m}$,

$\gamma_{j}(z,w):=\int_{0}^{1}P_{j}(\zeta)f(z(\zeta h),w(\zeta
h))\mathrm{d}\zeta.$ (46)

On the other hand, we shall look for a piecewise approximation to
$\hat{\sigma}(t)$, i.e. $\sigma(t)$, such that (compare with (11)–(15))

$\sigma_{n}(ch):=\sigma(t_{n-1}+ch),\qquad c\in[0,1],\qquad n=1-\nu,\dots,N,$
(47)

denotes its restriction to the time interval $[t_{n-1},t_{n}]$. Consequently,
one has:

$\sigma_{n}(ch)\equiv\phi(t_{n-1}+ch),\qquad c\in[0,1],\qquad
n=1-\nu,\dots,0,$ (48)

whereas, for $n=1,\dots,N$, $\sigma_{n}\in\Pi_{s}$ satisfies the differential
equation

$\dot{\sigma}_{n}(ch)=\sum_{j=0}^{s-1}P_{j}(c)\gamma_{j}(\sigma_{n},\sigma_{n-\nu}),\qquad
c\in[0,1],\qquad\sigma_{n}(0)=y_{n-1},$ (49)

so that,

$\sigma_{n}(ch)=y_{n-1}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\gamma_{j}(\sigma_{n},\sigma_{n-\nu}),\qquad
c\in[0,1],$ (50)

and

$y_{n}=y_{n-1}+h\gamma_{0}(\sigma_{n},\sigma_{n-\nu})=:\sigma_{n}(h),$ (51)

with $\gamma_{j}(\sigma_{n},\sigma_{n-\nu})$ defined according to (46). In the
sequel, we shall discuss the accuracy of the approximations:

$\displaystyle y(t_{n})-y_{n}$ $\displaystyle\equiv$
$\displaystyle\hat{\sigma}_{n}(h)-\sigma_{n}(h),$ (52)
$\displaystyle\delta\sigma_{n}(ch)$ $\displaystyle:=$
$\displaystyle\hat{\sigma}_{n}(ch)-\sigma_{n}(ch),\quad c\in(0,1),\qquad
n=1,\dots,N.$

For this purpose, some preliminary results are given in the next section.

###  3.1 Preliminary results

We start with the generalization of Corollary 1 to the present setting.

###### Corollary 3

With reference to (46), one has: $\gamma_{j}(z,w)=O(h^{j})$.

Proof Immediate from Theorem 1, by setting $G(\zeta h):=f(z(\zeta h),w(\zeta
h))$.

We also need perturbation results corresponding to those of Theorem 2 for
ODEs. For this purpose, it is sufficient to discuss them for a local problem
defined on two contiguous time subintervals of width $\tau$: the former
containing the memory, the latter containing the solution to be computed.
Without loss of generality, we shall then fix the reference interval
$[t_{0}-\tau,t_{0}+\tau]$, where we consider the following problem, defined
for a generic $\xi\in[t_{0},t_{0}+\tau]$:

$\displaystyle\dot{y}(t)$ $\displaystyle=$ $\displaystyle
f(y(t),y(t-\tau)),\qquad t\in[t_{0},t_{0}+\tau],\qquad
y(\xi)=\eta\in\mathbb{R}^{m},$ (53) $\displaystyle y(t)$ $\displaystyle=$
$\displaystyle\phi(t),\qquad\qquad\qquad\quad\,t\in[t_{0}-\tau,t_{0}).$

Problem (53) defines a generalization of the localized one associated to (38)
(obtained for $\xi=t_{0}$ and $\eta=y_{0}$), and we shall denote its solution
by

$y(t)\equiv y(t,\xi,\eta,\phi;t_{0}),$ (54)

in order to emphasize its dependence on the first four parameters, whereas the
last one refers to the time subinterval. We shall also use the following
notation:

$F_{1}(z,w)=\frac{\partial}{\partial z}f(z,w),\qquad
F_{2}(z,w)=\frac{\partial}{\partial w}f(z,w).$ (55)

###### Remark 5

It is clear that the function $\phi$ in (53) represents the memory term of the
equation, and it is a known function. The same will happen in the subsequent
reference interval, $[t_{0},t_{0}+2\tau]$, obtained by shifting to the right
the previous one by $\tau$, once the solution of (53) has been computed, and
so forth.

To begin with, let us state the following straightforward result, whose proof
is omitted for brevity.

###### Theorem 13

The solution (54) of problem (53) is defined on the whole time interval
$[t_{0},t_{0}+\tau]$, independently of the point $\xi\in[t_{0},t_{0}+\tau]$
where the condition $\eta$ is given.

The following result then holds true (compare with Theorem 2).

###### Theorem 14

With reference to the solution (54) of problem (53), one has:

$a)~{}\frac{\partial}{\partial t}y(t)=f(y(t),y(t-\tau)),\quad
b)~{}\frac{\partial}{\partial\eta}y(t)=\Phi(t,\xi;t_{0}),\quad
c)~{}\frac{\partial}{\partial\xi}y(t)=-\Phi(t,\xi;t_{0})f(y(\xi),y(\xi-\tau)),$

where $\Phi(t,\xi;t_{0})$ satisfies (see (55)):

$\displaystyle\dot{\Phi}(t,\xi;t_{0})$ $\displaystyle=$ $\displaystyle
F_{1}(y(t),y(t-\tau))\Phi(t,\xi;t_{0}),\qquad t\in[t_{0},t_{0}+\tau],$
$\displaystyle\Phi(\xi,\xi;t_{0})$ $\displaystyle=$ $\displaystyle
I\in\mathbb{R}^{m\times m},$ (56) $\displaystyle\Phi(t,\xi;t_{0})$
$\displaystyle=$ $\displaystyle O\in\mathbb{R}^{m\times m},\qquad
t\in[t_{0}-\tau,t_{0}).$

Proof The statement $a)$ clearly follows from (53). From the same equation one
also derives that, for $t\in[t_{0},t_{0}+\tau]$,

$\displaystyle\qquad\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial}{\partial\eta}y(t)\right)$
$\displaystyle=$
$\displaystyle\frac{\partial}{\partial\eta}\dot{y}(t)~{}=~{}\frac{\partial}{\partial\eta}f(y(t),y(t-\tau))$
$\displaystyle=$ $\displaystyle
F_{1}(y(t),y(t-\tau))\frac{\partial}{\partial\eta}y(t)\,+\,F_{2}(y(t),y(t-\tau))\frac{\partial}{\partial\eta}y(t-\tau).$

Moreover, at $t=\xi$,

$\frac{\partial}{\partial\eta}y(\xi)=\frac{\partial}{\partial\eta}\eta=I,$

and, for $t\in[t_{0}-\tau,t_{0})$,

$\frac{\partial}{\partial\eta}y(t)=\frac{\partial}{\partial\eta}\phi(t)=O.$

This latter equality implies that, for $t\in[t_{0},t_{0}+\tau)$, the term
$F_{2}(y(t),y(t-\tau))\frac{\partial}{\partial\eta}y(t-\tau)$ in (3.1)
vanishes, thus reducing to the first equation in (56), so that $b)$ eventually
follows. Finally, by virtue of Theorem 13, let $t^{*}$ be a generic point in
the interval $[t_{0},t_{0}+\tau]$, and denote

$y^{*}=y(t^{*},\xi,\eta,\phi;t_{0}).$

Consequently, since $\xi\in[t_{0},t_{0}+\tau]$ as well, one has:

$\eta=y(\xi,t^{*},y(t^{*},\xi,\eta,\phi;t_{0}),\phi;t_{0}),$

so that we eventually arrive at the identity

$y^{*}=y(t^{*},\xi,y(\xi,t^{*},y^{*},\phi;t_{0}),\phi;t_{0}).$

By taking into account the results of the previous points $a)$ and $b)$, one
derives:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\xi}y(t^{*},\xi,y(\xi,t^{*},y^{*},\phi;t_{0}),\phi;t_{0})$
$\displaystyle=$
$\displaystyle\frac{\partial}{\partial\xi}y(t^{*},\xi,\eta,\phi;t_{0})+\frac{\partial}{\partial\eta}y(t^{*},\xi,\eta,\phi;t_{0})\left.\frac{\partial}{\partial
t}y(t,t^{*},y^{*},\phi;t_{0})\right|_{t=\xi}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\xi}y(t^{*},\xi,\eta,\phi;t_{0})+\Phi(t^{*},\xi;t_{0})f(y(\xi),y(\xi-\tau)).$

The statement $c)$ then follows, by taking into account that $t^{*}$ is
generic.

One main difference with the ODE case, stems from the fact that now (54) also
depends on the memory term $\phi$, which is a functional parameter.
Consequently, we now look for a Frechét derivative such that, for any
perturbation $\delta\phi\in C([t_{0}-\tau,t_{0}])$ and
$t\in[t_{0}-\tau,t_{0}+\tau]$:

$\lim_{\varepsilon\rightarrow
0}\frac{y(t,\xi,\eta,\phi+\varepsilon\delta\phi;t_{0})-y(t,\xi,\eta,\phi;t_{0})}{\varepsilon}=\int_{t_{0}-\tau}^{t_{0}}\frac{\delta}{\delta\phi(\zeta)}y(t)\delta\phi(\zeta)\mathrm{d}\zeta,$
(58)

where

$\frac{\delta}{\delta\phi(\zeta)}y(t)\,:\,(t,\zeta)\in[t_{0}-\tau,t_{0}+\tau]\times[t_{0}-\tau,t_{0})\rightarrow\begin{pmatrix}\frac{\delta}{\delta\phi_{j}(\zeta)}y_{i}(t)\end{pmatrix}\in\mathbb{R}^{m\times
m},$ (59)

is the functional derivative of (54) (see, e.g., [24, Appendix A]), with
$y_{i}$ and $\phi_{j}$ the respective entries of $y$ and $\phi$. For later
use, we recall that, for a given $\hat{t}\in[t_{0}-\tau,t_{0})$ and
$i,j=1,\dots,m$,

$\displaystyle~{}~{}~{}\frac{\delta}{\delta\phi_{j}(\hat{t}\,)}y_{i}(t)$
$\displaystyle\equiv$ $\displaystyle\lim_{\varepsilon\rightarrow
0}\frac{y_{i}(t,\xi,\eta,\phi+\varepsilon
e_{j}\delta_{\hat{t}};t_{0})-y_{i}(t,\xi,\eta,\phi;t_{0})}{\varepsilon}$
$\displaystyle=$
$\displaystyle\int_{t_{0}-\tau}^{t_{0}}\frac{\delta}{\delta\phi_{j}(\zeta)}y_{i}(t)\delta_{\hat{t}}(\zeta)\mathrm{d}\zeta,$

with $e_{j}\in\mathbb{R}^{m}$ the $j$th unit vector and, hereafter,
$\delta_{\hat{t}}(t)$ is the Dirac delta function centered at $\hat{t}$. The
following result holds true.

###### Lemma 1

With reference to (58) and (59), for any fixed
$t\in[\xi,t_{0}+\tau]\subseteq[t_{0},t_{0}+\tau]$ one has:

$\frac{\delta}{\delta\phi(\zeta)}y(t)=O\in\mathbb{R}^{m\times
m},\qquad\forall\zeta\in[t_{0}-\tau,\xi-\tau)\cup(t-\tau,t_{0}).$

Proof Having fixed $t\in[\xi,t_{0}+\tau]$, it follows that
$\forall\zeta\in[t_{0}-\tau,\xi-\tau)\cup(t-\tau,t_{0})$, setting as usual
$\delta_{\zeta}$ the Dirac delta centered at $\zeta$, one has:

$y(t,\xi,\eta,\phi+\varepsilon\delta_{\zeta};t_{0})=y(t,\xi,\eta,\phi;t_{0}),\qquad\forall\varepsilon\in\mathbb{R}.$

In fact, by virtue of (53), the solution (54) is independent of the values of
$\phi$ outside the interval $[\xi-\tau,t-\tau]$. Consequently, by taking into
account (3.1), it follows that:

$\frac{\delta}{\delta\phi(\zeta)}y(t)=\lim_{\varepsilon\rightarrow
0}\frac{y(t,\xi,\eta,\phi+\varepsilon\delta_{\zeta};t_{0})-y(t,\xi,\eta,\phi;t_{0})}{\varepsilon}=O.\,\framebox{}$

Taking into account Lemma 1, the following result provides a more practical
characterization of the functional derivative (58)–(59). Figure 1 displays the
location of the most relevant points and subintervals involved in Theorem 15.

Figure 1: Relevant time subintervals for Theorem 15.

###### Theorem 15

With reference to the solution (54) of problem (53), for any
$\hat{t}\in(\xi-\tau,t_{0})$ one has:

$\frac{\delta}{\delta\phi(\hat{t}\,)}y(t)=\Psi(t,\hat{t}\,;t_{0}),$ (61)

where $\Psi(t,\hat{t}\,;t_{0})$ satisfies (see (55)):

$\displaystyle\dot{\Psi}(t,\hat{t}\,;t_{0})$ $\displaystyle=$ $\displaystyle
F_{1}(y(t),y(t-\tau))\Psi(t,\hat{t}\,;t_{0}),\qquad
t\in(\hat{t}+\tau,t_{0}+\tau],$
$\displaystyle\quad\Psi(\hat{t}+\tau,\hat{t}\,;t_{0})$ $\displaystyle=$
$\displaystyle F_{2}(y(\hat{t}+\tau),y(\hat{t}\,)),$ (62)
$\displaystyle\Psi(t,\hat{t};t_{0})$ $\displaystyle=$
$\displaystyle\delta_{\hat{t}}(t)I,\qquad t\in[t_{0}-\tau,\hat{t}+\tau),$

with $O,I\in\mathbb{R}^{m\times m}$ and $\delta_{\hat{t}}(t)$ the Dirac delta
function.

Proof In fact, for $t\in[t_{0},t_{0}+\tau]$ one has, by virtue of (53):

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\frac{\delta}{\delta\phi(\hat{t}\,)}y(t,\xi,\eta,\phi;t_{0})~{}=~{}\frac{\delta}{\delta\phi(\hat{t}\,)}\dot{y}(t)~{}=~{}\frac{\delta}{\delta\phi(\hat{t}\,)}f(y(t),y(t-\tau))$
$\displaystyle=$ $\displaystyle
F_{1}(y(t),y(t-\tau))\frac{\delta}{\delta\phi(\hat{t}\,)}y(t)+F_{2}(y(t),y(t-\tau))\frac{\delta}{\delta\phi(\hat{t}\,)}y(t-\tau),$

i.e., using the notation (61),

$\displaystyle\qquad\dot{\Psi}(t,\hat{t}\,;t_{0})$ $\displaystyle=$
$\displaystyle
F_{1}(y(t),y(t-\tau))\Psi(t,\hat{t}\,;t_{0})+F_{2}(y(t),y(t-\tau))\Psi(t-\tau,\hat{t}\,;t_{0}).$
(63)

Moreover, at $t=\xi$,

$\displaystyle\Psi(\xi,\hat{t}\,;t_{0})$ $\displaystyle=$
$\displaystyle\frac{\delta}{\delta\phi(\hat{t}\,)}y(\xi)=\frac{\delta}{\delta\phi(\hat{t}\,)}\eta=O,$
(64)

since the condition $y(\xi)=\eta$ is independent of the history $\phi$.
Further, taking into account (3.1), for all $i,j=1,\dots,m$ and for
$t\in[t_{0}-\tau,t_{0})$, one has:

$\frac{\delta}{\delta\phi_{j}(\hat{t}\,)}y_{i}(t)=\frac{\delta}{\delta\phi_{j}(\hat{t}\,)}\phi_{i}(t)=\lim_{\varepsilon\rightarrow
0}\frac{[\phi_{i}(t)+\varepsilon\delta_{ij}\delta_{\hat{t}}(t)]-\phi_{i}(t)}{\varepsilon}=\delta_{ij}\delta_{\hat{t}}(t),$

with $\delta_{ij}$ the Kronecker delta. Consequently,

$\displaystyle\Psi(t,\hat{t}\,;t_{0})$ $\displaystyle=$
$\displaystyle\frac{\delta}{\delta\phi(\hat{t}\,)}y(t)=\delta_{\hat{t}}(t)I,\qquad
t\in[t_{0}-\tau,t_{0}).$ (65)

From (63)–(65), one then derives, considering that (see Figure 1)
$t_{0}-\tau\leq\xi-\tau<\hat{t}:$

$\displaystyle\quad\Psi(t,\hat{t}\,;t_{0})$ $\displaystyle=$
$\displaystyle\int_{\xi}^{t}\dot{\Psi}(\zeta,\hat{t}\,;t_{0})\mathrm{d}\zeta=\left\\{\begin{array}[]{cc}O,&t\in[t_{0},\hat{t}+\tau),\\\\[5.69054pt]
F_{2}(y(\hat{t}+\tau),y(\hat{t}\,)),&t=\hat{t}+\tau.\end{array}\right.$ (68)

From (65) and (68) the last two equations in (62) follow. Consequently, from
(63), one obtains

$\dot{\Psi}(t,\hat{t}\,;t_{0})=F_{1}(y(t),y(t-\tau))\Psi(t,\hat{t}\,;t_{0}),\qquad
t\in(\hat{t}+\tau,t_{0}+\tau],$

which completes the proof of (62).

As a straightforward consequence, the following result holds true, which
guarantees the regularity of $\Psi$ w.r.t. its first two arguments (again, for
sake of clarity, refer to Figure 1).

###### Corollary 4

With reference to the solution (54) of problem (53), and considering (55),
(56), and (62), for any $\hat{t}\in(\xi-\tau,t_{0})$ one has:

$\Psi(t,\hat{t}\,;t_{0})=\Phi(t,\hat{t}+\tau;t_{0})F_{2}(y(\hat{t}+\tau),y(\hat{t}\,)),\qquad
t\in[\hat{t}+\tau,t_{0}+\tau].$ (69)

Finally, the following result holds true (compare with Corollary 2 of the ODE
case).

###### Corollary 5

With reference to the solution (54) of problem (53), and considering (56) and
(69), for any $\delta\phi\in C([t_{0}-\tau,t_{0}])$ one has:

$\displaystyle y(t,\xi,\eta+\delta\eta,\phi+\delta\phi;t_{0})$
$\displaystyle=$ $\displaystyle
y(t,\xi,\eta,\phi;t_{0})+\Phi(t,\xi;t_{0})\delta\eta+\int_{\xi-\tau}^{t-\tau}\Psi(t,\zeta;t_{0})\delta\phi(\zeta)\mathrm{d}\zeta$
$\displaystyle+\,(t-\xi)\,O(|\delta\eta|+\|\delta\phi\|)^{2},\qquad\qquad
t\in[\xi,t_{0}+\tau],$

with $\|\delta\phi\|=\max_{\zeta\in[\xi-\tau,t-\tau]}|\delta\phi(\zeta)|$.

Proof The statement follows from Theorem 14, part $b)$, and Theorem 15, by
taking into account (58) and the result of Lemma 1.

###  3.2 Main results (DDE case)

We are now in the position of discussing the accuracy of the approximations
(52). To begin with, the following result holds true.

###### Theorem 16

With reference to (38)–(52), for $n=1,\dots,\nu$ one has:

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+O(h^{2s+1}),\quad\|\delta\sigma_{n}\|:=\max_{c\in[0,1]}|\delta\sigma_{n}(ch)|=O(h^{s+1}).$

Proof The statement follows from Theorem 4 by considering that, for
$n=1,\dots,\nu$, $t\in[t_{0},t_{0}+\tau]$ in (38), so that
$y(t-\tau)\equiv\phi(t-\tau)$, which is a known function, thus obtaining an
ODE.

This result allows us to state the following one, which generalizes that of
Theorem 5 to the present case.

###### Theorem 17

With reference to (44), (46), (50), and (52) if for $n\geq 1$ one has:

$y(t_{r})-y_{r}=y(t_{r-1})-y_{r-1}+O(h^{2s+1}),\qquad r=1,\dots,n,$

then

$\delta\gamma_{j}^{n}:=\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu})-\gamma_{j}(\sigma_{n},\sigma_{n-\nu})=O(h^{2s-j}),\qquad
j=0,\dots,s-1.$

Proof The proof is by generalized induction. For $n=1,\dots,\nu$ the statement
follows from Theorem 5 and Theorem 16 since, in this case,

$\hat{\sigma}_{n-\nu}(ch)\equiv\sigma_{n-\nu}(ch)\equiv\phi(t_{n-1}+ch-\tau),\qquad
c\in[0,1],$

so that we are dealing with an ODE. Assume now it true up to $n-1$, and prove
for $n$. By hypothesis, and from (45) and (51), we know that

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+O(h^{2s+1})\equiv
y(t_{n-1})-y_{n-1}+h\delta\gamma_{0}^{n},$

so that $\delta\gamma_{0}^{n}=O(h^{2s})$ follows. Then, by taking into account
(55), it follows that:

$\displaystyle O(h^{2s})$ $\displaystyle=$
$\displaystyle\delta\gamma_{0}^{n}\,=\,\gamma_{0}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu})-\gamma_{0}(\sigma_{n},\sigma_{n-\nu})$
$\displaystyle=$ $\displaystyle\int_{0}^{1}\left[f(\hat{\sigma}_{n}(\zeta
h),\hat{\sigma}_{n-\nu}(\zeta h))-f(\sigma_{n}(\zeta h),\sigma_{n-\nu}(\zeta
h))\right]\mathrm{d}\zeta$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\left[\,F_{1}(\sigma_{n}(\zeta
h)+c\,\delta\sigma_{n}(\zeta h),\sigma_{n-\nu}(\zeta
h)+c\,\delta\sigma_{n-\nu}(\zeta h))\delta\sigma_{n}(\zeta h)~{}+\right.$
$\displaystyle\qquad\quad~{}\left.F_{2}(\sigma_{n}(\zeta
h)+c\,\delta\sigma_{n}(\zeta h),\sigma_{n-\nu}(\zeta
h)+c\,\delta\sigma_{n-\nu}(\zeta h))\delta\sigma_{n-\nu}(\zeta
h)\,\right]\mathrm{d}c\,\mathrm{d}\zeta$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\underbrace{\int_{0}^{1}F_{1}(\sigma_{n}(\zeta
h)+c\,\delta\sigma_{n}(\zeta h),\sigma_{n-\nu}(\zeta
h)+c\,\delta\sigma_{n-\nu}(\zeta h))\mathrm{d}c}_{=:\,G_{1}(\zeta
h)}\,\delta\sigma_{n}(\zeta h)\,\mathrm{d}\zeta~{}+$
$\displaystyle\int_{0}^{1}\underbrace{\int_{0}^{1}F_{2}(\sigma_{n}(\zeta
h)+c\,\delta\sigma_{n}(\zeta h),\sigma_{n-\nu}(\zeta
h)+c\,\delta\sigma_{n-\nu}(\zeta h))\mathrm{d}c}_{=:\,G_{2}(\zeta
h)}\,\delta\sigma_{n-\nu}(\zeta h)\,\mathrm{d}\zeta$ $\displaystyle=$
$\displaystyle\int_{0}^{1}G_{1}(\zeta h)\delta\sigma_{n}(\zeta
h)\mathrm{d}\zeta\,+\,\int_{0}^{1}G_{2}(\zeta h)\delta\sigma_{n-\nu}(\zeta
h)\mathrm{d}\zeta.$

Let us discuss in detail the term

$\int_{0}^{1}G_{1}(\zeta h)\delta\sigma_{n}(\zeta h)\,\mathrm{d}\zeta,$

since the remaining one,

$\int_{0}^{1}G_{2}(\zeta h)\delta\sigma_{n-\nu}(\zeta
h)\,\mathrm{d}\zeta=O(h^{2s}),$

is similarly discussed, by taking into account the induction hypothesis. By
virtue of (2.2), one has:

$\displaystyle\delta\sigma_{n}(\zeta h)$ $\displaystyle=$
$\displaystyle\hat{\sigma}_{n}(\zeta h)-\sigma_{n}(\zeta h)$ $\displaystyle=$
$\displaystyle
y(t_{n-1})-y_{n-1}\,+\,h\sum_{j=0}^{s-1}\int_{0}^{\zeta}P_{j}(x)\mathrm{d}x\,\delta\gamma_{j}^{n}\,+\,\sum_{j\geq
s}\int_{0}^{\zeta}P_{j}(x)\mathrm{d}x\,\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu})$
$\displaystyle=$ $\displaystyle y(t_{n-1})-y_{n-1}\,+\,\zeta
h\delta\gamma_{0}^{n}+h\sum_{j=1}^{s-1}\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\delta\gamma_{j}^{n}\,+$
$\displaystyle h\sum_{j\geq
s}\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu}).$

Consequently, from Theorem 1 and Corollary 3, one obtains:

$\displaystyle O(h^{2s})~{}=~{}\int_{0}^{1}G_{1}(\zeta
h)\delta\sigma_{n}(\zeta
h)\,\mathrm{d}\zeta~{}=~{}\underbrace{\int_{0}^{1}G_{1}(\zeta
h)\mathrm{d}\zeta}_{=\,O(1)}\underbrace{[y(t_{n-1})-y_{n-1}]}_{=\,(n-1)\,O(h^{2s+1})}$
$\displaystyle+\,h\underbrace{\int_{0}^{1}G_{1}(\zeta
h)\zeta\mathrm{d}\zeta}_{=\,O(1)}\,\underbrace{\delta\gamma_{0}^{n}}_{=\,O(h^{2s})}\,+\,h\sum_{j=1}^{s-1}\int_{0}^{1}G_{1}(\zeta
h)\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\mathrm{d}\zeta\,\delta\gamma_{j}^{n}$
$\displaystyle+\underbrace{h\sum_{j\geq s}\overbrace{\int_{0}^{1}G_{1}(\zeta
h)\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\mathrm{d}\zeta}^{=\,O(h^{j-1})}\,\overbrace{\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu})}^{=\,O(h^{j})}}_{=\,O(h^{2s})},$

from which,

$O(h^{2s})=h\sum_{j=1}^{s-1}\underbrace{\int_{0}^{1}G_{1}(\zeta
h)\left[\xi_{j+1}P_{j+1}(\zeta)-\xi_{j}P_{j-1}(\zeta)\right]\mathrm{d}\zeta}_{=\,O(h^{j-1})}\delta\gamma_{j}^{n}$

follows and, therefore, one concludes that $\delta\gamma_{j}^{n}=O(h^{2s-j})$,
$j=0\dots,s-1$.

As a consequence, the following result can be stated.

###### Theorem 18

With reference to (38)–(52), for $n=1,\dots,N\equiv K\nu$ one has:

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+O(h^{2s+1}),\quad\|\delta\sigma_{n}\|:=\max_{c\in[0,1]}|\delta\sigma_{n}(ch)|=O(h^{s+1}).$

Proof The proof is done by induction on groups of $\nu$ consecutive steps. For
the first $\nu$ steps, the statement follows from Theorem 16. Assume now, by
induction, that it holds true up to $t_{k\nu}=t_{0}+k\nu h$, and let us prove
for $n=k\nu+1,\dots,(k+1)\nu$. For this purpose, for $k=1,\dots,K-1$ let us
set:

$\phi_{k}(t)\equiv\sigma(t),\qquad\hat{\phi}_{k}(t)\equiv\hat{\sigma}(t)\equiv
y(t),\qquad t\in[t_{(k-1)\nu},t_{k\nu}).$

Assuming, again, true the statement for $n-1$, and using the notation (54),
one has:

$\displaystyle\delta\sigma_{n}(ch)\,=\,\hat{\sigma}_{n}(ch)-\sigma_{n}(ch)$
$\displaystyle=$ $\displaystyle
y(t_{n-1}+ch,t_{n-1},y(t_{n-1}),\hat{\phi}_{k};t_{k\nu})-y(t_{n-1}+ch,t_{n-1}+ch,\sigma_{n}(ch),\phi_{k};t_{k\nu})$
$\displaystyle=$
$\displaystyle\underbrace{y(t_{n-1}+ch,t_{n-1},\overbrace{\sigma_{n}(0)}^{=\,y_{n-1}},\phi_{k};t_{k\nu})-y(t_{n-1}+ch,t_{n-1}+ch,\sigma_{n}(ch),\phi_{k};t_{k\nu})}_{=:\,E_{n,1}^{(k)}(ch)}$
$\displaystyle+~{}\underbrace{y(t_{n-1}+ch,t_{n-1},y(t_{n-1}),\hat{\phi}_{k};t_{k\nu})-y(t_{n-1}+ch,t_{n-1},\overbrace{y_{n-1}}^{=\,\sigma_{n}(0)},\phi_{k};t_{k\nu})}_{=:\,E_{n,2}^{(k)}(ch)}.$

From Theorem 16, it follows that

$\displaystyle E_{n,1}^{(k)}(h)$ $\displaystyle=$
$\displaystyle\sigma_{n}(0)-\sigma_{n}(0)+O(h^{2s+1})\,=\,O(h^{2s+1}),$ (70)
$\displaystyle\|E_{n,1}^{(k)}\|$ $\displaystyle:=$
$\displaystyle\max_{c\in[0,1]}|E_{n,1}^{(k)}(ch)|=O(h^{s+1}).$

Moreover, from Corollary 5, and considering that $h\nu=\tau$, one has:

$\displaystyle E_{n,2}^{(k)}(ch)$ $\displaystyle=$
$\displaystyle\overbrace{\Phi(t_{n-1}+ch,t_{n-1};t_{k\nu})}^{=\,I+O(ch)}\overbrace{\left[y(t_{n-1})-y_{n-1}\right]}^{=\,\delta\sigma_{n}(0)}$
$\displaystyle+\,h\int_{0}^{c}\Psi(t_{n-1}+ch,t_{n-1}+\zeta
h-\tau;t_{k\nu})\delta\sigma_{n-\nu}(\zeta h)\mathrm{d}\zeta$
$\displaystyle+\,ch\,O(|\delta\sigma_{n}(0)|+\|\delta\sigma_{n-\nu}\|)^{2}.$

By considering that

$|\delta\sigma_{n}(0)|=(n-1)O(h^{2s+1}),\qquad\|\delta\sigma_{n-\nu}\|=O(h^{s+1}),$

one eventually derives

$\|E_{n,2}^{(k)}\|:=\max_{c\in[0,1]}|E_{n,2}^{(k)}(ch)|=O(h^{s+2}),$

from which the second part of the statement follows, by taking into account
(70). Moreover, when $c=1$ then $t_{n-1}+h=t_{n}$ and, by virtue of Theorem 1
and Theorem 17, one obtains:

$\displaystyle\int_{0}^{1}\overbrace{\Psi(t_{n},t_{n-1}+\zeta
h-\tau;t_{k\nu})}^{=:\,G(\zeta h)}\delta\sigma_{n-\nu}(\zeta
h)\mathrm{d}\zeta$ $\displaystyle=$ $\displaystyle\int_{0}^{1}G(\zeta
h)\left[\sum_{j=0}^{s-1}P_{j}(\zeta)\delta\gamma_{j}^{n-\nu}+\sum_{j\geq
s}P_{j}(\zeta)\gamma_{j}(\hat{\sigma}_{n-\nu},\hat{\sigma}_{n-2\nu})\right]$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{s-1}\underbrace{\int_{0}^{1}P_{j}(\zeta)G(\zeta
h)\mathrm{d}\zeta}_{=\,O(h^{j})}\underbrace{\delta\gamma_{j}^{n-\nu}}_{=\,O(h^{2s-j})}+\sum_{j\geq
s}\underbrace{\int_{0}^{1}P_{j}(\zeta)G(\zeta
h)\mathrm{d}\zeta}_{=\,O(h^{j})}\underbrace{\gamma_{j}(\hat{\sigma}_{n-\nu},\hat{\sigma}_{n-2\nu})}_{=\,O(h^{j})}~{}=~{}O(h^{2s}).$

Consequently,

$E_{n,2}^{(k)}(h)=y(t_{n-1})-y_{n-1}+O(h^{2s+1}),$

and also the first part of the statement follows.

###  3.3 Discretization

The discretization issue proceeds as in the ODE case. In fact, also in the DDE
case, the Fourier coefficients (see (46) and (49)),

$\gamma_{j}(\sigma_{n},\sigma_{n-\nu})=\int_{0}^{1}P_{j}(\zeta)f(\sigma_{n}(\zeta
h),\sigma_{n-\nu}(\zeta h))\mathrm{d}\zeta,\qquad j=0,\dots,s-1,$

need to be approximated by using a (interpolatory) quadrature rule of order
$q$, thus providing a possibly different piecewise approximation $u(t)$,

$u(t)\equiv\phi(t),\quad t<t_{0},\qquad u_{n}(ch):=u(t_{n-1}+ch),\quad
c\in[0,1],\quad n=1-\nu,\dots,N,$

such that, for $n\geq 1$:

$\dot{u}_{n}(ch)=\sum_{j=0}^{s-1}P_{j}(c)\hat{\gamma}_{j}(u_{n},u_{n-\nu}),\qquad
c\in[0,1],\qquad u_{n}(0)=y_{n-1}.$ (71)

Consequently,

$u_{n}(ch)=y_{n-1}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\hat{\gamma}_{j}(u_{n},u_{n-\nu}),\qquad
c\in[0,1],$ (72)

and

$y_{n}=y_{n-1}+h\hat{\gamma}_{0}(u_{n},u_{n-\nu})=:u_{n}(h),$ (73)

where (see (46)),

$\hat{\gamma}_{j}(u_{n},u_{n-\nu}):=\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(u_{n}(c_{i}h),u_{n-\nu}(c_{i}h))=\gamma_{j}(u_{n},u_{n-\nu})-\Delta_{j}(h),$
(74)

with $(c_{i},b_{i})$ the abscissae and weights of the quadrature, and
$\Delta_{j}(h)=O(h^{q-j})$ the quadrature error, where $q$ is the order of the
quadrature.

Formulae (71) and (46) form a subclass of the so called natural continuous RK
methods for DDEs (see [5, Sec. 6.2]). As a consequence, their convergence
properties could be as well derived by more classical approaches such as
Bellman’s method of steps, which is an analytic procedure specific for DDEs.
In the present context, the main goal is to show how the framework based on
the perturbation theory applied to the truncated Fourier expansion is easily
adapted to cope with DDEs, therefore we will pursue this route of
investigation. A further strength of this approach is the possibility of
analyzing the convergence properties of the truncated Fourier approximations
when these are used as spectral methods in time. In this regard, the analysis
for the ODE case has been addressed in [3], while a spectral implementation of
the methods for DDEs has been considered in [16].

By using standard arguments (which we omit, as done in the ODE case), we can
derive the following results, representing the corresponding counterparts of
Theorem 17 and Theorem 18, respectively.

###### Theorem 19

With reference to (44), (46), (71)–(73), and assuming that the quadrature
formula (74) has order $q\geq 2s$, if for $n\geq 1$ one has:

$y(t_{r})-y_{r}=y(t_{r-1})-y_{r-1}+O(h^{2s+1}),\qquad r=1,\dots,n,$

then

$\delta\hat{\gamma}_{j}^{n}:=\gamma_{j}(\hat{\sigma}_{n},\hat{\sigma}_{n-\nu})-\hat{\gamma}_{j}(u_{n},u_{n-\nu})=O(h^{2s-j}),\qquad
j=0,\dots,s-1.$

###### Theorem 20

With reference to (44), (46), (71)–(73), and assuming that the quadrature
formula (74) has order $q\geq 2s$, for $n=1,\dots,N\equiv K\nu$ one has:

$y(t_{n})-y_{n}=y(t_{n-1})-y_{n-1}+O(h^{2s+1}),\qquad\max_{c\in[0,1]}|\hat{\sigma}_{n}(ch)-u_{n}(ch)|=O(h^{s+1}).$

###### Remark 6

It is worth mentioning that the result of Theorem 20 states that the super-
convergence order $2s$ at the mesh-points $t_{n}$ is obtained, even though
possibly different Runge-Kutta methods are used at each integration step,
provided that they define a polynomial approximation of degree $s$. This, in
turn, represents a generalization of the results in [4] for collocation
methods.

We conclude this section, by recalling that the considerations in Remark 4
continue to hold in the DDE case and by observing that, concerning the
implementation of the resulting Runge-Kutta method used for solving problem
(38), the arguments in Section 2.5, mutatis mutandis, apply as well.

##  4 Numerical tests

In this section we report a few numerical tests for the DDE case. In fact, in
the ODE case, HBVMs have been extensively used as energy-conserving methods
for Hamiltonian systems (see, e.g., [2, 9, 10, 11, 19]). We show that, under
some circumstances, their use can be advantageous also in the DDE case.
Hereafter, we consider a class of DDEs defined by a Hamiltonian function

$H:(q,p)\in\mathbb{R}^{m}\times\mathbb{R}^{m}\rightarrow\mathbb{R},$ (75)

through the equations

$\displaystyle\dot{q}(t)$ $\displaystyle=$ $\displaystyle\quad
H_{p}(q(t),p(t))+\alpha H_{p}(q(t-\tau),p(t-\tau)),$ $\displaystyle\dot{p}(t)$
$\displaystyle=$ $\displaystyle-\left[H_{q}(q(t),p(t))+\alpha
H_{q}(q(t-\tau),p(t-\tau))\right],\quad$

with $\alpha$ a real parameter, $\tau>0$ the delay, and $H_{q}$ and $H_{p}$
the partial derivatives of $H$ w.r.t. $q$ and $p$, respectively. The problem
is completed by the initial conditions

$q(t)=\phi(t),\qquad p(t)=\psi(t),\qquad t\in[-\tau,0].$ (77)

The introduction of such a kind of delay Hamiltonian system is partly inspired
by the problem of looking for periodic orbits of DDEs, which has been attacked
by many authors in the past (see, e.g., [23, 34, 37, 38, 39, 35, 41]). In this
respect, the first two examples below show an attractive periodic orbit with
integer period lying on a level set of the Hamiltonian function (75) which is,
therefore, a constant of motion once the periodic orbit has been approached.
In the third example we are instead interested in simulating the correct
qualitative behavior of a dissipative Hamiltonian delay problem in the phase
space when the dynamics takes place in a neighborhood of a separatrix. Taking
aside a theoretical discussion of problem (75)-(LABEL:eq:qpdde) which would go
beyond the scopes of the present work, we infer its properties for the three
considered instances by preliminarily applying a high order integrator with
very small stepsize, in order to get a very accurate numerical solution that
will be taken as a reference trajectory in the phase space.

For all the three problems, we show that a very accurate approximation of the
Hamiltonian function allows us to reproduce the correct geometric features of
the solution in the discrete setting. To the best of our knowledge, this is
the first instance of the use of HBVMs in the context of DDEs displaying
geometric properties. For comparison purposes, we also solve the problems with
the classical Gauss collocation integrator of the same order. The numerical
tests have been implemented in Matlab (R2020b) on a 3 GHz Intel Xeon W10 core
computer with 64GB of memory.

#### Problem 1

With reference to (75)–(77), the first problem is defined as follows:

$m=1,\qquad
H(q,p)=\frac{1}{4}\left(q^{4}+p^{4}\right),\qquad\alpha=10^{-1},\qquad\tau=1,\qquad\phi(t)\equiv\sqrt{2},\qquad\psi(t)\equiv
0.$ (78)

We solve this problem by using the following methods:

  * •

HBVM(2,2) (i.e., the 2-stage Gauss method),

  * •

HBVM(4,2).

Both methods are fourth-order, according to Theorem 10, with HBVM(4,2) energy-
conserving once a periodic orbit of integer period is eventually reached (see
Theorem 12). Problem (78) possesses an attracting periodic orbit with period
$T=2\tau=2$ which suggests using a stepsize $h$ equal to a submultiple of
$\tau$, in order to mimic a corresponding discrete periodic solution. As we
are going to see, unlike the 2-stage Gauss collocation method, the
conservation property of HBVM(4,2) results in a precise resolution of this
task. We solve the problem on the interval $[0,2\cdot 10^{3}]$ by using a
timestep $h=\tau/5=0.2$. Figure 2 summarizes the obtained results.

  * •

In the upper row of the figure are the plots of the numerical Hamiltonian,
$H(q_{n},p_{n})$, from which one deduces that both methods quite soon reach a
stationary behavior.

  * •

To better discern the asymptotic behavior of the two numerical solutions, the
central pictures show the plots of
$\left|H(q_{n},p_{n})-H(q_{n-1},p_{n-1})\right|$ for the two methods. From
these plots one infers that, while the stationary value of the Hamiltonian is
constant for the HBVM(4,2) method, it is only approximately constant for the
HBVM(2,2) method, with oscillations having amplitude of order $10^{-2}$.

  * •

The bottom row contains the plots of the numerical trajectory in the phase
plane for both methods, relative to the interval $[2\cdot 10^{2},2\cdot
10^{3}]$ (i.e., after the transient phase). For both methods, the solution
seems to repeat every 10 points (i.e., with period $T=2\tau=2$). However, the
points obtained by the HBVM(2,2) method are not actually periodic, whereas
they are (within to machine precision and independently of the used stepsize
$h$) for the HBVM(4,2) method. To confirm this, in Table 1 we list the last 20
points of the trajectories computed after each period $T=10h$ and lying inside
the two small circles highlighted in the plots. As one may see, only the first
4 digit of the points of the trajectory computed by the HBVM(2,2) method are
retained, whereas the points computed by the HBVM(4,2) method differ at most
on the last digit.

Table 1: The last 20 points of the trajectories inside the circles in the plots on the bottom line of Figure 2. HBVM(2,2) | HBVM(4,2)
---|---
$q$ | $p$ | $q$ | $p$
1.344913051657652 | 1.924341608176171 | 1.364023296679203 | 1.918490612087558
1.344895222079097 | 1.924347768593204 | 1.364023296679201 | 1.918490612087558
1.344877390115245 | 1.924353929533504 | 1.364023296679201 | 1.918490612087558
1.344859555766308 | 1.924360090996891 | 1.364023296679201 | 1.918490612087558
1.344841719032502 | 1.924366252983184 | 1.364023296679200 | 1.918490612087558
1.344823879914032 | 1.924372415492206 | 1.364023296679201 | 1.918490612087558
1.344806038411116 | 1.924378578523775 | 1.364023296679200 | 1.918490612087558
1.344788194523964 | 1.924384742077714 | 1.364023296679201 | 1.918490612087559
1.344770348252789 | 1.924390906153841 | 1.364023296679202 | 1.918490612087558
1.344752499597800 | 1.924397070751980 | 1.364023296679202 | 1.918490612087558
1.344734648559217 | 1.924403235871946 | 1.364023296679200 | 1.918490612087559
1.344716795137255 | 1.924409401513561 | 1.364023296679199 | 1.918490612087559
1.344698939332121 | 1.924415567676645 | 1.364023296679200 | 1.918490612087559
1.344681081144034 | 1.924421734361018 | 1.364023296679199 | 1.918490612087559
1.344663220573212 | 1.924427901566499 | 1.364023296679200 | 1.918490612087558
1.344645357619866 | 1.924434069292906 | 1.364023296679200 | 1.918490612087558
1.344627492284208 | 1.924440237540062 | 1.364023296679200 | 1.918490612087558
1.344609624566458 | 1.924446406307785 | 1.364023296679200 | 1.918490612087559
1.344591754466830 | 1.924452575595894 | 1.364023296679201 | 1.918490612087558
1.344573881985545 | 1.924458745404208 | 1.364023296679203 | 1.918490612087558

Figure 2: Numerical results for problem (78) solved by using HBVM(2,2), left
plots, and HBVM(4,2), right plots, using a timestep $h=0.2$ (see the text for
details).

#### Problem 2

The second example is similar in nature to the previous one but considers a
non-polynomial Hamiltonian function with two degrees of freedom. With
reference to (75)–(77), it is defined by:

$\displaystyle m=2,\qquad
H(q,p)=\frac{1}{4}\left(q_{1}^{4}+q_{2}^{4}+p_{1}^{4}+p_{2}^{4}\right)+\frac{\pi}{2}\left(\frac{1}{\|q\|_{2}^{2}}+\frac{2}{\|p\|_{2}^{2}}\right),$
(79) $\displaystyle\alpha=5\cdot
10^{-2},\qquad\tau=1,\qquad\phi(t)\equiv(0.1,~{}1)^{\top},\qquad\psi(t)\equiv(1,~{}0.2)^{\top}.$

Again, we have experienced the existence of a periodic orbit with period
$T=2\tau=2$. We solve this problem on the interval $[0,10^{3}]$ with timestep
$h=\tau/10=0.1$, by using the following methods:

  * •

HBVM(2,2) (i.e., the 2-stage Gauss method),

  * •

HBVM(10,2).

Both methods are fourth-order, the latter being practically energy-conserving,
for the given timestep, in the event that a periodic orbit is reached.

Also in this case, the conservation property of HBVM(10,2) turns out to be
crucial in reproducing a discrete orbit with period precisely equal to 2,
while a small phase drift affects the solution yielded by the 2-stage Gauss
collocation method. Figure 3, which is similar to Figure 2, summarizes the
obtained results.

  * •

In the upper row of the figure are the plots of the numerical Hamiltonian,
namely $H(q_{n},p_{n})$: for both methods it seems to reach a stationary
behavior.

  * •

The second row shows the plots of
$\left|H(q_{n},p_{n})-H(q_{n-1},p_{n-1})\right|$ for the two methods. From
these plots one infers that, while the stationary value of the Hamiltonian is
constant (up to round-off) for the HBVM(10,2) method, it is only approximately
constant for the HBVM(2,2) method, with oscillations having amplitude of order
$10^{-3}$.

  * •

The bottom row contains the plots of the numerical trajectory in the
$q_{1}-q_{2}$ plane for both methods, relative to the interval
$[10^{2},10^{3}]$ (i.e., after the transient phase). For both methods, the
solution seems to repeat every 20 points (i.e., with period $T=2\tau=2$).
However, only the points obtained by the HBVM(10,2) method are actually
periodic. To confirm this, in Table 2 we list the last 20 points of the
trajectories computed after each period $T=20h$ and lying inside the two small
circles displayed in the plots. As one may see, the Gauss collocation method
only retain the first 5 digits after each period, whereas the points computed
by the HBVM(10,2) method differ at most on the last digit.

Table 2: The last 20 points of the trajectories inside the circles in the plots on the bottom line of Figure 3. HBVM(2,2) | HBVM(10,2)
---|---
$q_{1}$ | $q_{2}$ | $q_{1}$ | $q_{2}$
1.500006047618583 | 1.868403720200248 | 1.595245320422993 | 1.813631211153069
1.500014079966090 | 1.868399733165645 | 1.595245320422992 | 1.813631211153067
1.500022113284264 | 1.868395745553114 | 1.595245320422991 | 1.813631211153069
1.500030147573164 | 1.868391757362593 | 1.595245320422994 | 1.813631211153068
1.500038182832831 | 1.868387768594017 | 1.595245320422991 | 1.813631211153067
1.500046219063319 | 1.868383779247324 | 1.595245320422993 | 1.813631211153069
1.500054256264677 | 1.868379789322453 | 1.595245320422994 | 1.813631211153069
1.500062294436967 | 1.868375798819335 | 1.595245320422991 | 1.813631211153067
1.500070333580236 | 1.868371807737912 | 1.595245320422991 | 1.813631211153068
1.500078373694533 | 1.868367816078119 | 1.595245320422993 | 1.813631211153069
1.500086414779919 | 1.868363823839889 | 1.595245320422992 | 1.813631211153067
1.500094456836429 | 1.868359831023165 | 1.595245320422991 | 1.813631211153069
1.500102499864130 | 1.868355837627883 | 1.595245320422994 | 1.813631211153068
1.500110543863066 | 1.868351843653977 | 1.595245320422991 | 1.813631211153067
1.500118588833287 | 1.868347849101385 | 1.595245320422993 | 1.813631211153069
1.500126634774853 | 1.868343853970044 | 1.595245320422994 | 1.813631211153069
1.500134681687813 | 1.868339858259893 | 1.595245320422991 | 1.813631211153067
1.500142729572213 | 1.868335861970863 | 1.595245320422991 | 1.813631211153068
1.500150778428111 | 1.868331865102895 | 1.595245320422993 | 1.813631211153069
1.500158828255558 | 1.868327867655927 | 1.595245320422992 | 1.813631211153067
1.500166879054607 | 1.868323869629889 | 1.595245320422991 | 1.813631211153069
1.500174930825299 | 1.868319871024725 | 1.595245320422994 | 1.813631211153068

Figure 3: Numerical results for problem (79) solved by using HBVM(2,2), left
plots, and HBVM(10,2), right plots, using a timestep $h=0.1$ (see the text for
details).

#### Problem 3

For the last problem, we are no more interested in periodic trajectories.
Instead, we consider a delay Hamiltonian problem with dissipation. This can be
achieved by choosing a negative value of the parameter $\alpha$ in
(LABEL:eq:qpdde). With reference to (75)–(77), the selected parameters are:

$m=1,\quad H(q,p)=\frac{1}{2}p^{2}-\cos
q,\quad\alpha=-10^{-5},\quad\tau=1,\quad\phi(t)\equiv 0,\quad\psi(t)\equiv
1.99999.$ (80)

This problem is a dissipative delay-variant of the nonlinear pendulum, with
the initial condition chosen close to the separatrix (the level set
$H(q,p)=1$) between the two different regimes of the pendulum: librations
around the straight-down stationary position, and rotations. For the given
initial conditions, the pendulum should undergo damped oscillations with a
decreasing trend of the Hamiltonian function $H(q_{n},p_{n})$. Consequently,
when using relatively large stepsizes, it is fundamental to reproduce the
correct dissipation of the Hamiltonian along the numerical trajectory.

We solve this problem on the interval $[0,500]$, with a timestep
$h=\tau/2=0.5$, by using the following methods:

  * •

HBVM(2,2) (i.e., the 2-stage Gauss method),

  * •

HBVM(10,2).

Figure 4 summarizes the obtained results.

  * •

In the upper row of the figure are the plots of the numerical Hamiltonian,
$H(q_{n},p_{n})$, from which one deduces that both methods have a dissipation
trend of the energy $H$. Nevertheless, for HBVM(2,2) the values of the
Hamiltonian becomes quite larger than 1 in the initial part of the trajectory
and undergoes fictitious oscillations which cause the numerical solution to
escape the correct region of the phase space where the dynamics should take
place, as we are going to see. This is not the case for the HBVM(10,2) method,
whose numerical Hamiltonian decreases in the correct way, thus remaining
always smaller than 1.

  * •

The central pictures show the numerical solution in the phase space. As one
may see, the numerical solution provided by HBVM(2,2) “jumps” twice, before
being trapped into an invariant region. This means that the pendulum undergoes
two complete rotations until it looses enough energy and begins oscillating
around the rest position. On the contrary, the numerical solution obtained by
using HBVM(10,2) always remains in the correct region.

  * •

The bottom row contains the plots of the numerical solution w.r.t. time,
confirming that the numerical solution provided by the HBVM(2,2) method
“jumps” twice, whereas that obtained by the HBVM(10,2) method does not.

Figure 4: Numerical results for problem (80) solved by using HBVM(2,2), left
plots, and HBVM(10,2), right plots, using a timestep $h=0.5$ (see the text for
details).

##  5 Conclusions

In this paper we have fully developed a thorough approach for obtaining
polynomial approximations to the solution of initial value ODE and DDE
problems. It allows us to derive a wide class of Runge-Kutta methods, whose
properties are easily discussed within the framework, as well as their actual
implementation. Some numerical tests, concerning the numerical simulation of
solutions of certain DDE problems of Hamiltonian type, confirm this. The
present approach leaves room for generalizations along several directions: in
particular to different kind of problems, besides the ones considered here.
Another relevant direction of investigation consists in looking for
approximations belonging to functional subspaces different than polynomials:
that is, by considering orthonormal functional bases different from (3). Both
directions will be the subject of future investigations.

## Funding and conflicts of interests

The authors have no affiliation with any organization with a direct or
indirect financial interest in the subject matter discussed in the manuscript.
The authors acknowledge the financial support from the mrSIR crowdfunding
[42].

## References

  * [1] P. Amodio, L. Brugnano, and F. Iavernaro. A note on the continuous-stage Runge-Kutta-(Nyström) formulation of Hamiltonian Boundary Value Methods (HBVMs). Appl. Math. Comput., 363 (2019) 124634. https://doi.org/10.1016/j.amc.2019.124634
  * [2] P. Amodio, L. Brugnano, and F. Iavernaro. Continuous-Stage Runge-Kutta Approximation to Differential Problems. Axioms, 11 (2022) 192. https://doi.org/10.3390/axioms11050192
  * [3] P. Amodio, L. Brugnano, and F. Iavernaro. Analysis of Spectral Hamiltonian Boundary Value Methods (SHBVMs) for the numerical solution of ODE problems. Numer. Algorithms, 83 (2020) 1489–1508. https://doi.org/10.1007/s11075-019-00733-7
  * [4] A. Bellen. One step collocation for delay differential equations. J. Comput. Appl. Math., 10 (1984) 275–283. https://doi.org/10.1016/0377-0427(84)90039-6
  * [5] A. Bellen and M. Zennaro. Numerical Methods for Delay Differential Equations. Clarendon Press, Oxford, 2003.
  * [6] P. Betsch and P. Steinmann. Conservation properties of a time FE method. I. Time-stepping schemes for $N$-body problems. Internat. J. Numer. Methods Engrg., 49 (2000) 599–638. https://doi.org/10.1002/1097-0207(20001020)49:5<599::AID-NME960>3.0.CO;2-9
  * [7] C.L. Bottasso. A new look at finite elements in time: a variational interpretation of Runge-Kutta methods. Appl. Numer. Math., 25 (1997) 355–368. https://doi.org/10.1016/S0168-9274(97)00072-X
  * [8] L. Brugnano, G. Frasca-Caccia, and F. Iavernaro. Efficient implementation of Gauss collocation and Hamiltonian Boundary Value Methods. Numer. Algorithms, 65 (2014) 633–650. http://doi.org/10.1007/s11075-014-9825-0
  * [9] L. Brugnano and F. Iavernaro. Line Integral Methods for Conservative Problems. Chapman et Hall/CRC, Boca Raton, FL, 2016. https://doi.org/10.1201/b19319
  * [10] L. Brugnano and F. Iavernaro. Line Integral Solution of Differential Problems. Axioms, 7(2) (2018) 36. https://doi.org/10.3390/axioms7020036
  * [11] L. Brugnano, F. Iavernaro, J.I. Montijano, and L. Rández. Spectrally accurate space-time solution of Hamiltonian PDEs. Numer. Algorithms, 81 (2019) 1183–1202. https://doi.org/10.1007/s11075-018-0586-z
  * [12] L. Brugnano, F. Iavernaro, and D. Trigiante. Hamiltonian Boundary Value Methods (Energy Preserving Discrete Line Integral Methods). JNAIAM J. Numer. Anal. Ind. Appl. Math., 5, no. 1-2 (2010) 17–37.
  * [13] L. Brugnano, F. Iavernaro, and D. Trigiante. A note on the efficient implementation of Hamiltonian BVMs. J. Comput. Appl. Math., 236 (2011) 375–383. https://doi.org/10.1016/j.cam.2011.07.022
  * [14] L. Brugnano, F. Iavernaro, and D. Trigiante. A simple framework for the derivation and analysis of effective one-step methods for ODEs. Appl. Math. Comput., 218 (2012) 8475–8485. https://doi.org/10.1016/j.amc.2012.01.074
  * [15] L. Brugnano, F. Iavernaro, and D. Trigiante. Analisys of Hamiltonian Boundary Value Methods (HBVMs): a class of energy-preserving Runge-Kutta methods for the numerical solution of polynomial Hamiltonian systems. Commun. Nonlinear Sci. Numer. Simul., 20 (2015) 650-667. https://doi.org/10.1016/j.cnsns.2014.05.030
  * [16] L. Brugnano, F. Iavernaro, and P. Zanzottera. A multiregional extension of the SIR model, with application to the COVID-19 spread in Italy. Math. Meth. Appl. Sci., 44 (2021) 4414–4427. https://doi.org/10.1002/mma.7039
  * [17] L. Brugnano and C. Magherini. Blended implementation of block implicit methods for ODEs. Appl. Numer. Math., 42 (2002) 29–45. https://doi.org/10.1016/S0168-9274(01)00140-4
  * [18] L. Brugnano and C. Magherini. Recent advances in linear analysis of convergence for splittings for solving ODE problems. Appl. Numer. Math., 59 (2009) 542–557. https://doi.org/10.1016/j.apnum.2008.03.008
  * [19] L. Brugnano, J.I. Montijano, and L. Rández. On the effectiveness of spectral methods for the numerical solution of multi-frequency highly-oscillatory Hamiltonian problems. Numer. Algorithms, 81 (2019) 345–376. http://dx.doi.org/10.1007/s11075-018-0552-9
  * [20] H. Brunner. Collocation Methods for Volterra Integral and Related Functional Equations. Cambridge University Press, Cambridge, 2004.
  * [21] E. Celledoni, R.I. McLachlan, D. McLaren, B. Owren, G.R.W. Quispel, and W.M. Wright. Energy preserving Runge-Kutta methods. M2AN Math. Model. Numer. Anal., 43 (2009) 645–649. https://doi.org/10.1051/m2an/2009020
  * [22] G. Dahlquist and Å. Björk. Numerical Methods in Scientific Computing, Volume I. SIAM, Philadelphia, 2008.
  * [23] J.G. Dos Reis and R.L. Baroni. On the existence of periodic solutions for autonomous retarded functional-differential equations on $R^{2}$. Proc. Roy. Soc. Edinburgh Sect., A 102 (1986) 259–262. https://doi.org/10.1017/S0308210500026342
  * [24] E. Engel and R.M. Dreizler. Density Functional Theory, an advanced course. Springer, Berlin, 2011.
  * [25] D. Furihata and T. Matsuo. Discrete Variational Derivative Method: A Structure-Preserving Numerical Method for Partial Differential Equations. Chapman and Hall/CRC, Boca Raton, FL, 2010.
  * [26] E. Hairer. Energy-preserving variants of collocation methods. JNAIAM J. Numer. Anal. Ind. Appl. Math., 5, no. 1-2 (2010) 73–84.
  * [27] E. Hairer, S.P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I, nonstiff problems. Second revised edition (3rd printing). Springer, Heidelberg, 2008.
  * [28] E. Hairer and G. Wanner. Solving Ordinary Differential Equations I, nonstiff problems. Second revised edition. Springer, Heidelberg, 2002.
  * [29] B.L. Hulme. One-step piecewise polynomial Galerkin methods for initial value problems. Math. Comp., 26 (1972) 415–426. https://doi.org/10.1090/S0025-5718-1972-0321301-2
  * [30] B.L. Hulme. Discrete Galerkin and related one-step methods for ordinary differential equations. Math. Comp., 26 (1972) 881–891. https://doi.org/10.1090/S0025-5718-1972-0315899-8
  * [31] F. Iavernaro and B. Pace. $s$-Stage trapezoidal methods for the conservation of Hamiltonian functions of polynomial type. AIP Conf. Proc. 936 (2007) 603–606. https://doi.org/10.1063/1.2790219
  * [32] F. Iavernaro and B. Pace. Conservative block-boundary value methods for the solution of polynomial Hamiltonian systems. AIP Conf. Proc. 1048 (2008) 888–891. https://doi.org/10.1063/1.2991075
  * [33] F. Iavernaro, and D. Trigiante. High-order symmetric schemes for the energy conservation of polynomial Hamiltonian problems. JNAIAM J. Numer. Anal. Ind. Appl. Math., 4, no. 1-2 (2009) 87–111.
  * [34] J.L. Kaplan and J.A. Yorke. Ordinary differential equations which yield periodic solutions of differential delay equations. J. Math. Anal. Appl., 48 (1974) 317–324. https://doi.org/10.1016/0022-247X(74)90162-0
  * [35] J. Mallet-Paret and R.D. Nussbaum. Stability of periodic solutions of state-dependent delay-differential equations. J. Differential Equations, 250 (2011) 4085–4103. https://doi.org/10.1016/j.jde.2010.10.023
  * [36] Y. Miyatake and J.C. Butcher. A characterization of energy-preserving methods and the construction of parallel integrators for Hamiltonian systems. SIAM J. Numer. Anal., 54 (2016) 1993–2013. https://doi.org/10.1137/15M1020861
  * [37] R.D. Nussbaum. Periodic solutions of some nonlinear, autonomous functional differential equations. Bull. Amer. Math. Soc., 79 (1973) 811–814. https://doi.org/10.1016/0022-0396(73)90053-3
  * [38] R.D. Nussbaum. Periodic solutions of some nonlinear, autonomous functional differential equations. II. J. Differential Equations, 14 (1973) 360–394. https://doi.org/10.1090/S0002-9904-1973-13330-0
  * [39] R.D. Nussbaum. Uniqueness and nonuniqueness for periodic solutions of $x^{\prime}(t)=-g(x(t-1))$. J. Differential Equations, 34 (1979) 25–54. https://doi.org/10.1016/0022-0396(79)90016-0
  * [40] G.R.W. Quispel and D.I. McLaren. A new class of energy-preserving numerical integration methods. J. Phys. A, 41 (2008) 045206. https://doi.org/10.1088/1751-8113/41/4/045206
  * [41] H.-O. Walther. Existence of a non-constant periodic solution of a nonlinear autonomous functional differential equation representing the growth of a single species population. J. Math. Biol., 1 (1975) 227–240. https://doi.org/10.1007/BF01273745
  * [42] https://www.mrsir.it/en/about-us/