# Multiscale DeepONet for Nonlinear Operators in Oscillatory Function Spaces
for Building Seismic Wave Responses

Lizuo Liu  Wei Cai  Department of Mathematics, Southern Methodist University,
Dallas, TX 75275, USA

###### Abstract

In this paper, we propose a multiscale DeepONet to represent nonlinear
operator between Banach spaces of highly oscillatory continuous functions. The
multiscale deep neural network (DNN) utilizes a multiple scaling technique to
convert high frequency function to lower frequency functions before using a
DNN to learn a specific range of frequency of the function. The multi-scale
concept is integrated into the DeepONet which is based on a universal
approximation theory of nonlinear operators. The resulting multi-scale
DeepONet is shown to be effective to represent building seismic response
operator which maps oscillatory seismic excitation to the oscillatory building
responses.

###### keywords:

Neural network, universal approximation theory of nonlinear operator.

Suggested Running Head:
Multiscale DeepONet for Operators in Oscillatory Function Spaces

Corresponding Author:
Prof. Wei Cai
Department of Mathematics,
Southern Methodist University,
Dallas, TX 75275
Email<EMAIL_ADDRESS>

##  1 Introduction

Computing operators between physical quantities defined in function spaces
have many applications in forward and inverse problems in scientific and
engineering computations. For example, in wave scattering in inhomogeneous or
random media, the mapping between the media physical properties, which can be
modelled as a random field, and the wave field is a nonlinear operator, which
embodies some of most challenging computational tasks in medical imaging,
geophysical and seismic problems. Another example in earthquake safety of
buildings and structures, the response of structures to the seismic waves
gives rise to another important operator between spaces of highly oscillatory
functions. Due to the infinite dimension of function operators, active
researches have been done recently in efficient representation of both forward
and inverse operators with reduced models and model reductions. It is all
nature that deep neural network has been considered as a machine learning tool
to learn these operators. DNNs have been used for the field of scientific and
engineering computing has shown much promise [1][2][4] [3] [5][6]. Recently,
the DeepOnet[7] was constructed based on the theory of Universal Approximation
of Nonlinear Operators [8] for learning operators between Banach spaces. Other
approach for learning operators include Fourier neural operator [12]. In this
paper,we will further develop the DeepONet for operators between highly
oscillatory function spaces such as those encountered in the building seismic
wave response problems and a multiscale DeepONet based on multiscale DNN
concept will be investigated.

The traditional way to get the responses of buildings is by the finite element
method[9], which could take $O\left(10^{2}\right)$ to $O\left(10^{3}\right)$
seconds for the computation of single case, depending on the grids of the
finite elements. The method we developed uses the multiscale deep neural
network[6] and the DeepONet to learn the mapping between the seismic
excitations and the corresponding displacements of the building at all levels
excited by the seismic waves, which could give the displacements in
$O\left(10^{-2}\right)$ seconds with relative accurate results.

The paper is constructed as follows. In section 2, we introduce the
theoretical basis of the seismic wave response operator and the idea of data
augmentation of the linear differential equation case. In section 3, we give a
short review of the multiscale deep neural network and the integration of
multiscale idea into the DeepONet structure. In section 4, we introduce the
data processing procedure, including the computation of responses and the
elimination of the aliasing of the seismic records. The numerical results are
shown in section 5 and finally, a conclusion is given in section 6.

##  2 Building seismic wave response operator $\cal R$

Building response to seismic waves The problem we want to study is the
dynamics behavior of buildings to the seismic waves[10]. We want to predict
the response of the buildings instantly given the seismic excitations over a
time period $[0,T]$, which defines an operator

$\mathcal{R}:P(t)\mapsto x(t),\quad t\in[0,T]$ (1)

where $P$ is the seismic wave excitation and $x$ is the displacements of the
building excited by the seismic wave excitation $P$. We assume that the
dynamics of a deformed body can be described by the following equation of
motion after an appropriate finite element discretization:

$\mathbf{M}\ddot{x}\left(t\right)=P-F+H-C\dot{x}$ (2)

with initial conditions

$\left\\{\begin{aligned} &x\left(0\right)=0\\\
&x^{\prime}\left(0\right)=0\end{aligned}\right.,$ (3)

where $\mathbf{M}$ is the global mass matrix, $P$ accounts for the global load
vector (nodal load, body force, surface force, etc.), $C$ is the damping
matrix, $H$ is the global hour-glass resisting force vector handling the hour-
glass deformation modes, and $F$ is the assembly of equivalent nodal force
vectors from all the elements

$F=\sum_{i=1}^{n}\int_{V_{i}}\mathbf{B}^{T}\iota dV$ (4)

where $\mathbf{B}$ is the strain-displacement matrix, $\iota$ is the stress
vector, and $V$ is the volume in the current configuration.

The $\mathbf{M},C,F$ and $H$ are given for the operator fitting problem. The
inputs and outputs of the neural network should be the global load vector
$P\left(t\right)$ and the displacements of $x\left(t\right)$ at specific
moments $t_{1},t_{2},\ldots,t_{n}$.

If we only consider $1$ level of the building and only $1$ node in that level,
$x$ in equation $\left(\ref{Solid}\right)$ is a 1D time series, otherwise $x$
is an mD time series if we only set 1 node but considering $m$ levels.

The superposition of the solutions Consider the training dataset
$\\{\left(p_{i},x_{i}\right),i\in\mathcal{I}\\}$, where the seismic wave
excitations $\\{p_{i}\\}$ and the corresponding displacements $\\{x_{i}\\}$
satisfy the equation $\left(\ref{Solid}\right)$, we build every training data
to be fed to the neural network by

$\left\\{\begin{aligned} &\mathcal{P}=\sum_{i\in\iota}w_{i}p_{i},\\\
&\mathcal{X}=\sum_{i\in\iota}w_{i}x_{i}\end{aligned}\right.$ (5)

where $\iota$ is a randomly sampling index subset from the index set
$\mathcal{I}$ whose size is fixed and $w_{i}$’s are the random weights
satisfying

$\sum_{i\in\iota}w_{i}=1$ (6)

Thus, $\mathcal{P}$ and $\mathcal{X}$ satisfy

$\left\\{\begin{aligned}
&\mathbf{M}\ddot{\mathcal{X}}=\mathcal{P}-F+H-C\dot{\mathcal{X}},\\\
&\mathcal{X}\left(0\right)=0,\\\
&\mathcal{X}^{\prime}\left(0\right)=0\end{aligned}\right.$ (7)

by the superposition principle of the nonhomogeneous case, which construct the
general idea of data augmentation for the case that there is a linear mapping
relation between inputs and outputs.

##  3 Multiscale DeepONet for operators in oscillatory function spaces

Solutions for wave propagation and scattering in inhomogeneous media can be
viewed through an operator between the media physical properties such as
permittivity, density, etc and the wave fields through their governing
equations including elasticity wave equations, Maxwell’s equations, and
acoustic wave equations. In high frequency wave problems, the operator will
have either its domain or range involves functions of highly oscillatory
nature.

We will extend the DeepONet in order to represent the operator $R$ between
highly oscillatory functions, namely, the seismic excitations, and the
building responses, both of which are observed to include a wide range of
frequencies. First of all, we define the so-called Tauber-Wiener function to
be used as activation function in DNN.

###### Definition 1

If a function $\sigma:\mathbb{R}\to\mathbb{R}$ $\left(\text{continuous or
discontinuous}\right)$ satisfies that all the linear combinations
$\sum_{i=1}^{I}c_{i}\sigma\left(\lambda_{i}t+\theta_{i}\right),\lambda_{i}\in\mathbb{R},\theta_{i}\in\mathbb{R},c_{i}\in\mathbb{R},i=1,\ldots,I$
are dense in every $C\left[a,b\right]$, then $\sigma$ is called a Tauber-
Wiener $\left(\text{TW}\right)$ function.

The DeepONet is based on the universal approximation for nonlinear operators.
This theory gives a constructive procedure for approximating nonlinear
operator between continuous functions $u(x)$ in $x\in K_{1}\subset\mathcal{X}$
and continuous functions in $y\in K_{2}\subset\mathcal{Y=}$ $\mathbb{R}^{n}$,

$\mathcal{G}:u=u(x)\in V\subset C(K_{1})\rightarrow\mathcal{G(}u)(y)\in
C(K_{2})$ (8)

based on two uniform approximation results:

  * •

Uniform approximation of functions: Given any $\varepsilon_{1},$ functions
$\mathcal{G(}u)(y)$ selected from a compact subset in $C(K_{2})$ can be
uniformed approximated by a two-layered neural network with any Tauber-Wiener
$\left(\text{TW}\right)$ activation function $\sigma_{t}$

$|\mathcal{G(}u)(y)-\sum_{k=1}^{N}c_{k}\cdot\sigma_{t}\left(\omega_{k}\cdot
y+\zeta_{k}\right)|\leq\varepsilon_{1},\text{ \ }\forall\text{ }y\in K_{2}$
(9)

where $c_{k}=c_{k}(u)\circeq c_{k}(\mathcal{G(}u))$, and
$\omega_{k},\zeta_{k}$ are all independent of $y$ and the function
$\mathcal{G(}u)(y)$ being approximated.

  * •

Uniform approximation of functionals by a finite dimensional nonlinear
operator: Given any $\varepsilon_{2},$ functional $c_{k}(u),u\in V$ (compact
subset of $C(K_{1})$) can be approximated by a two-layer neural network with
any Tauber-Wiener $\left(\text{TW}\right)$ activation function $\sigma_{b}$

$|c_{k}(u)-\sum_{i=1}^{M}c_{i}^{k}\sigma_{b}\left(\sum_{j=1}^{m}\xi_{ij}^{k}u\left(x_{j}\right)+\theta_{i}^{k}\right)|\leq\varepsilon_{2},\forall
u\in V,$ (10)

where the cofficients $c_{i}^{k},\xi_{ij}^{k},\theta_{i}^{k}$ and nodes
$\\{x_{j}\\}_{j=1}^{m}$ and $m,M$ are all independent of $u$.

Combining these two uniform approximation, we have an universal approximation
of nonlinear operator when restricted to a compact subset of continuous
function defined on a compact domain $K_{1}$. Namely, given any $\varepsilon,$
we can find $\\{x_{j}\\}_{j=1}^{m}\subset K_{1}\subset\mathcal{X},$coefficents
$c_{i}^{k},\xi_{ij}^{k},\omega_{k},\zeta_{k}$-all independent of continuous
functions $u$ $\in V\subset C(K_{1})$ and $y\in K_{2}\subset$ $\mathbb{R}^{n}$

$|\mathcal{G}(u)(y)-\sum_{k=1}^{N}\sum_{i=1}^{M}c_{i}^{k}\sigma_{b}\left(\sum_{j=1}^{m}\xi_{ij}^{k}u\left(x_{j}\right)+\theta_{i}^{k}\right)\cdot\sigma_{t}\left(\omega_{k}\cdot
y+\zeta_{k}\right)|<\epsilon.$ (11)

Multiscale DeepONet The outer sum of N-DNNs in (11) is identified as trunk in
the DeepONet with activation function $\sigma_{t}$, which is used for
providing uniform approximation of time dependent functions, i.e. building
response in our case, sampled from a compact subsect of $C^{0}[0,T]$; the
inner sum of M DNNs (11) with activation function $\sigma_{b}$ is used for
providing uniform approximation to functionals of infinite dimensions,
selected from a compact subspace of functionals, by m-dimensional nonlinear
operators. As the building response can be highly oscillatory, therefore, we
will apply the multscale DNN concept to generate a multiscale DeepONet as
shown in Fig. 1.

Figure 1: The schematics of the (bFCN-tMS) structure.

In principle, we can also apply the multiscale DNN concept to the branch
network, thus , we have several possible ways to combine the DeepONet and the
multiscale deep neural network. The main criteria that whether we need to set
specific scales for the input is whether there are high frequency components
with respect to the specific variables. From the deep learning viewpoint, The
DeepONet is a network with $m+n$ inputs and only 1 output, where $m$ is the
number of the inputs of the branch net, $n$ is the number of the inputs of the
trunk net. The inputs of the trunk net
$[P\left(t_{1}\right),P\left(t_{2}\right),\ldots P\left(t_{m}\right)]$ are the
evaluations at the scattered sensors $t_{1},t_{2},\ldots,t_{m}$. If we want to
analyse the frequency of response w.r.t every input components to determine
whether the branch net should be set as multiscale scheme, we need to
calculate the fourier transform with respect to each $P\left(t_{i}\right)$. To
do so, we need to interpolate a $m+n$ dimension function
$f:\mathbb{R}^{m+n}\to\mathbb{R}$. If a large size of sensors are used, the
analysis could be troublesome due to the difficulty of interpolating high
dimension data.

The DeepONet approximation for the seismic response operator in (1) for one
location $x(t)$ in the building will be given as for $t\in[0,T]$,

$x(t)=\mathcal{R}(P)(t)\approx\sum_{k=1}^{N}\sum_{i=1}^{M}c_{i}^{k}\sigma_{b}\left(\sum_{j=1}^{m}\xi_{ij}^{k}P\left(t_{j}\right)+\theta_{i}^{k}\right)\cdot\sigma_{t}\left(\omega_{k}\cdot
t+\zeta_{k}\right).$ (12)

Multiscales in DeepONet To simplify this procedure, we consider designing pre-
experiments for different combinations of the DeepONet and the Multiscale
DeepONet. The combinations we are interested in are

  * •

(bMS-tFCN) The branch net is multiscale neural network, the trunk net is the
general fully connected neural network.

  * •

(bFCN-tMS) The branch net is fully connected neural network, the trunk net is
the multiscale neural network.

  * •

(bMS-tMS) The branch net and the trunk net are both multiscale neural network.

  * •

(bFCN-tFCN) The branch net and the trunk net are both fully connect network.

The results can be seen in the section 5.1. We will select the bMS-tFCN
structure based on the results of pre-experiments. The schematics of the bFCN-
tMS structure is shown in Fig. 1.

Loss function Given batch size $\mathcal{N}$ for training process and the
total number of the test records $N$, the loss function is defined as

$\mathcal{L}\left(\theta\right)=\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\frac{1}{\text{max}_{j}|y_{ij}|}\Delta
t\sum_{j=1}^{m}\left(f_{NN}\left(t_{j},\mathbf{u}_{i},\theta\right)-y_{ij}\right)^{2}$
(13)

where $m$ is the length of signal. The penalty is set to be the reciprocal of
the maximum of the absolute value of the response since the discrepancy
regarding the magnitude of responses could be large corresponding to different
seismic records. The larger penalty should applied to assure that the neural
network could predict the response whose magnitude is smaller well.

To evaluate the training process, the mean of the relative L2 error is
considered

$\mathcal{L}\left(\theta\right)=\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\frac{\sqrt{\Delta
t\sum_{j=1}^{m}\left(f_{NN}\left(t_{j},\mathbf{u}_{i},\theta\right)-y_{ij}\right)^{2}}}{\sqrt{\Delta
t\sum_{j=1}^{m}\left(y_{ij}\right)^{2}}}$

The relative L2 error in a complete epoch is defined as

$\mathcal{L}^{\mathcal{R}}_{\text{train}}=\frac{1}{B}\sum_{k=1}^{B}\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\frac{\sqrt{\Delta
t\sum_{j=1}^{m}\left(f_{NN}\left(t_{j},\mathbf{u}_{i},\theta^{\left(k\right)}\right)-y_{ijk}\right)^{2}}}{\sqrt{\Delta
t\sum_{j=1}^{m}\left(y_{ijk}\right)^{2}}}$ (14)

where $B$ is the number of batches, $\theta^{\left(k\right)}$ means the
parameters of neural network at $k$-th batch.

Similarly, we define the relative L2 error for the testing dataset

$\mathcal{L}^{\mathcal{R}}_{\text{test}}=\frac{1}{{N}}\sum_{i=1}^{N}\frac{\sqrt{\Delta
t\sum_{j=1}^{m}\left(f_{NN}\left(t_{j},\mathbf{u}_{i},\theta\right)-y_{ij}\right)^{2}}}{\sqrt{\Delta
t\sum_{j=1}^{m}\left(y_{ij}\right)^{2}}}$ (15)

Note the $N$ is the total number of the test records.

##  4 Data preparation

To compute the responses of the building, we use the package openseespy[9] and
use Newmark’s method to solve equation $\left(\ref{Solid}\right)$ numerically.
We need to solve $N$ differential equations for the nodal displacements $x$,
and those $N$ differential equations could be transformed to several main
equations by expressing the displacements to the combination of the first few
natural vibration modes $\phi_{n}$. The $\phi_{n}$’s should be solved by the
Newmark’s Method. Thus we could obtain $x$ by the inverse transform with
respect to $\phi_{n}$. The details of solving the dynamics are introduced in
Chapter $16$ in [11].

Before solving the response, we need generate simulated seismic records since
the real seismic records are limited to obtain. We use the SeisMonGen[13]
package to generate the simulated seismic records. SeisMoGen is a software
that generate earthquake ground motion simulation based on the wavelet
decomposition under Priestley process assumption. For each specific real
records, the generator can generate 50 random samples.

The real seismic records should be downsampled before simulating the generated
samples by openseespy since the input size of the branch net is always fixed
and the same time step should be considered for different seismic records. The
original idea to re-sample one seismic record $\\{p_{n}\\}$ is retaining every
$K$th sample from the records such that the new records satisfies
$\tilde{f}_{n}=f_{nK}$, where $\tilde{f}_{n}$ is the down sampled sequence.
However, it is not valid to directly down-sample a sequence unless it is known
that the spectrum is $0$ at frequencies at or above the Nyquist frequency
defined by the lower sampling frequency. Since this will cause the phenomenon
called aliasing. Aliasing occurs when a signal is under-sampled. If the signal
sampling rate is small, there will be aliasing introduced.

The aliasing To discuss the relation of aliasing and downsampling, we follow
the notations in Julius O. Simth’s book[14], define

$\displaystyle\operatorname{DOWNSAMPLE}_{L,m}(x)$ $\displaystyle\triangleq
x(mL)$ (16) $\displaystyle m$ $\displaystyle=0,1,2,\ldots,M-1$ $\displaystyle
N$ $\displaystyle=LM$

and

$\displaystyle\operatorname{ALIAS}_{L,m}(x)$
$\displaystyle\triangleq\sum_{l=0}^{L-1}x(m+lM)$ (17) $\displaystyle m$
$\displaystyle=0,1,2,\ldots,M-1$ $\displaystyle N$ $\displaystyle=LM$

There is a downsampling theorem introduced,

######  Theorem 2 (Downsampling theorem)

For all $x\in\mathbf{C}^{N}$,

$\operatorname{DOWNSAMPLE}_{L}\left(x\right)\leftrightarrow\frac{1}{L}\operatorname{ALLAS}_{L}\left(X\right).$

From this theorem, we know that downsampling in the time domain will cause the
amplitudes in the frequency domain overlap, therefore extra frequencies will
pollute the expected frequency bands.

Thus before calculating the response, we need to eradicate aliasing before
downsampling. We will consider to use Butterworth filters to keep the low
frequency part but removing the high frequency part whose frequency is greater
than the half of Nyquist frequency.

Butterworth filters The class of continuous-time Butterworth filters[15] is
that for which the magnitude squared of the frequency response
$B\left(\omega\right)$ is of the form

$\left|B\left(\omega\right)\right|^{2}=\frac{1}{1+\left(\omega/\omega_{c}\right)^{2N}}$
(18)

The parameter $N$ is refered to as the filter order, The higher the order, the
sharper the transition from passband to the stopband. The $\omega_{c}$ is at
which $\left|B\left(\omega\right)\right|$ is at $1/\sqrt{2}$ its values at
$\omega=0$.

The magnitude squared of the frequency response of the class of discrete-time
Butterworth filters is of the form

$\left|B\left(\omega\right)\right|^{2}=\frac{1}{1+\left(\frac{\tan\left(\omega/2\right)}{\tan\left(\omega_{c}/2\right)}\right)^{2N}}$
(19)

The frequency response $B\left(\omega\right)$ satisfies the form

$B\left(s\right)=\frac{\omega_{c}^{N}}{\prod_{p=1}^{N}\left(s+s_{p}\right)}$
(20)

where

$s_{p}=\omega_{c}\exp\left\\{j\left[\frac{\pi\left(2p+1\right)}{2N}\frac{\pi}{2}\right]\right\\}$

is the pole of $B\left(j\omega\right)B^{*}\left(j\omega\right)$.

Amplitude separation in data If we plot the responses, we could conclude that
the high frequency has smaller magnitude, as shown in the Figure 2. The
structure now is

$f\left(P,t\right)=\sum_{i=1}^{n}\epsilon^{i}f_{i}\left(P,t\right)$ (21)

where $f_{i}$’s are DeepONet as shown in the schematics, but the trunk net of
each are multiscale nets with different scales, the schematics as shown in
Figure 3. Based on our observation, the data with larger magnitude contains
lower frequency, the data with smaller magnitude contains higher frequency.
Thus, we should configure the $f_{i}$ with more subnets as $i$ increases.

(a)

(b)

(c)

Figure 2: Several cases that have evident discrepancy in the magnitudes. The
small magnitude response has higher frequency components. Figure 3: The
multiscale DeepONet

We will compare the multiscale MSDeepONet with a large MSDeepONet in section
5, which shows huge improvements.

##  5 Numerical Results

###  5.1 Pre-experiements

Selection of scales - $l*n$ or $l^{n}$ The setting of scales of subnet of
multiscale deep neural network is still an open problem. We will compare two
different settings of scales of the subnets of multiscale deep neural network
in this subsection, one is that we set the scales as equal spaced,
$\\{l,2l,3l,\ldots,\kappa_{up}\\}$, the other is that we set the scale as
exponential spaced $\\{s^{0},s^{1},s^{2},\ldots,\kappa_{up}\\}$, where
$l=\frac{\kappa_{up}}{N},s=\kappa_{up}^{\frac{1}{N}}$ and $N$ is the number of
subnets. We set the branch net as fully connected neural network but the trunk
net as multiscale deep neural network with different scale settings, whose
capability will be proved in the next pre-experiments. In this experiment, a
random case will be selected and an upper bound frequency $\kappa_{up}$ will
be considered, which should be the upper bound of frequency of the
corresponding random case. The number of subnets is $30$,
$\kappa_{up}=60\times 2\pi$. As shown in Figure 4, The loss of the equal
spaced set up decreases faster and will be smaller compared with the
exponentially spaced case. In addition, this case give us another bonus
observation that to learn the oscillating case to a satisfactory accuracy, we
need the MSE loss to be $O\left(10^{-6}\right)$ at least.

Figure 4: The training loss of neural networks with different set up of
scales.

Determining DeepONet structure Like the previous experiment, we use only one
seismic record and the corresponding response for the DeepONet structure
determination.

The multiscale neural network considered as trunk net has $100$ four-layer
subnetworks with $10$ neurons in each layer, whose scales are
$\\{1,1+2\pi,\ldots,1+2k\pi,$ $\ldots,1+200\pi\\}$ respectively, and the
activation is $\sin\left(x\right)$. The output of trunk net is 1000 dimension.
The fully connected network considered as trunk net is also a four-layer
neural network, using $\sin\left(x\right)$ as activation function but with
$1000$ hidden neurons to keep the total number of hidden neurons to be the
same as the multiscale neural network.

The multiscale neural network considered as branch net has $100$ four-layer
subnetworks with $5$ neurons in each layer, whose scales are
$\\{1,1+2\pi,\ldots,1+2k\pi,$ $\ldots,1+200\pi\\}$ respectively, as the same
as the scales of multiscale net considered as trunk net. The activation is
also $\sin\left(x\right)$. The dimension of the output of branch net is as the
same as the dimension of output of trunk net. The fully connected network
considered as trunk net now shoule be a four-layer neural network with $500$
hidden neurons whose activation function is also $\sin\left(x\right)$.

Those four DeepONets with different structures will be trained up to $1500$
epochs. The training results can be shown in the Figure 5. As we can see, both
the bFCN-tMS DeepONet and the bMS-tMS DeepONet capture the frequencies, but
other 2 structures only learn some low-frequency waves after 1500-epoch
training. These results prove that there is no high frequency part that will
affect the learning for the inputs of branch net
$\left[P\left(t_{1}\right),P\left(t_{2}\right),\ldots
P\left(t_{m}\right)\right]$, or at least, we would not gain any improvements
if we use the scales $\\{1,1+2\pi,\ldots,1+2k\pi,$ $\ldots,1+200\pi\\}$ for
these inputs.

This results should be expected since solution $x\left(t\right)$ will contain
high frequency parts with respect to $t$ if the right hand side contains high
frequency part with respect to $t$ from the Equation
$\left(\ref{Solid}\right)$. Thus the trunk net whose input is $t$ should use
multiscale deep neural network to accelerate to capture the high frequency
part. Both multiscale neural net and fully connected neural net could be
considered as the proper choice for branch net, but considering training the
multiscale neural network will be more time-consuming, we select the fully
connected neural net as branch net. We will use the bFCN-tMS structure to
learn. The diagram of the bFCN-tMS structure can be shown in Fig. 1.

(a) (bFCN-tMS)

(b) (bFCN-tFCN)

(c) (bMS-tMS)

(d) (bMS-tFCN)

Figure 5: The output of the training data for different network structures. As
we can see, the bFCN-tMS (a) DeepONet captures the frequency the best, but the
one without multiscale DNN in the trunk of DeepONet bFCN-tFCN (b) and bMS-tFCN
(d) fails to do so in 100 epochs.

###  5.2 Results of multiscale DeepONet with amplitude separation

Predict the responses for a specific floor The multiscale DeepONet in this
section uses 3 sub-DeepONet, as shown in the Figure 3. For each sub-DeepONet,
the branch net is a 4-layer fully connected neural network with 128 hidden
neurons in each layer, whose activation function is Relu. The trunk net is a
multiscale neural network, but with different scales for different order of
$\varepsilon$. The DeepONet multiplying $\varepsilon^{2}$ at the output should
contain more scales since the responses whose amplitude is small contain
higher frequency parts. The DeepONet multiplying $\varepsilon^{0}$ at the
output do not need contain such scales since the responses whose amplitude is
large only contain lower frequency parts. We choose the scale
$\\{1,1+2\pi,\ldots,1+2k\pi,\ldots,1+200\pi\\}$ for the DeepONet multiplying
$\varepsilon^{2}$, the scale $\\{1,1+2\pi,\ldots,1+2k\pi,\ldots,1+100\pi\\}$
for the DeepONet multiplying $\varepsilon$, the scale
$\\{1,1+2\pi,\ldots,1+2k\pi,\ldots,1+20\pi\\}$ for the DeepONet multiplying
$\varepsilon^{0}$. Each subnets of those multiscale deep neural networks are
4-layer neural networks with 8 hidden neurons in each layer. As comparison, we
use a large multiscale DeepONet without amplitude separation trained by the
same datasets and tested by the same testing datasets. The trunk net of
multiscale DeepONet without amplitude separation is a multiscale neural
network which contains 100 four-layer subnets with $24$ hidden neurons in each
layer, whose scale is the same as the multiscale DeepONet with amplitude
seperation. The branch net of the multiscale DeepONet without amplitude
separation is the 4-layer fully connected neural network with $3\times 128$
hidden neurons in each layer, whose activation function is also relu. The
learning rate is $10^{-3}$ and total learning procedure contains 1500 epochs
with 40 batches per epoch.

(a)

(b)

Figure 6: The prediction of the multiscale deepONet. As shown in the figure,
the neural network capture the correct frequencies.

The results of the multi scale DeepONet are shown in the Figure 6, the neural
network predict the correct frequencies. By using the data augmentation, there
is no over-fitting issue and the prediction to the test cases show
satisfactory accuracy. The mean relative L2 error for testing cases at last is
$0.13$. The comparison of evolution of relative L2 error of training dataset
and testing dataset for the multiscale DeepONet with amplitude separation
comparing with what for multiscale DeepONet without separation are shown in
Figure 7. We could conclude that the amplitude separation idea indeed have
some contributions even though the data augmentation is also applied during
training procedure of multiscale DeepONet without amplitude separation.

(a) The relative L2 training error of the multiscale DeepONet with and without
amplitude separation

(b) The relative L2 testing error of the multiscale DeepONet with and without
amplitude separation

Figure 7: The multiscale deepONet with amplitude separation gives better
convergence results.

Prediction of the responses in multiple floors The preceding results show the
multiscale DeepONet with amplitude separation could predict the response of a
specific floor with a satisfactory accuracy. In this section, we show that the
multiscale DeepONet with amplitude separation not only could predict the
response of specific floor, but also could predict the responses for many
floors simultaneously, excited by the same seismic waves. The main difference
of the multiple-floor cases is that there are multiple outputs with different
kinds of frequencies for different outputs. Specifically, the responses of
lower floor contains higher frequencies with smaller amplitudes while the
responses of higher floor contains lower frequencies with larger amplitudes.
Such a amplitudes-frequency disparity can be handled with the proposed
amplitude separation in the multiscale DeepONet.

Now, the difference between the multiscale DeepONet here is that the output of
the neural network will be an $l$-D array. Other settings are the same as what
in the single-floor case. The results are shown in Fig. 8. We could conclude
that the multiscale DeepONet captured the desired frequencies with
satisfactory accuracy from the left figures of Fig. 8 and time responses from
the right figures.

(a) Floor 2

(b) Floor 3

(c) Floor 4

(d) Floor 5

(e) Floor 6

(f) Floor 7

(g) Floor 8

Figure 8: The predictions of the multiscale deepONet for seismic responses of
multiple floors. The neural network capture the correct frequencies (left
column) and time responses (right column).

##  6 Conclusion

In this paper, we proposed a multiscale DeepONet to handle the mappings
between oscillating inputs and outputs. For situations with the amplitude
separation in signals, we also proposed a specific structure to handle large
disparity of scales in amplitudes, resulting in satisfactory results. The
proposed multiscale DeepONet also shows its power when predicting the
responses containing different frequencies for different floors.

Future work The fast convergence of multiscale neural network needs rigorous
mathematical analysis even with lots of evidences of the advantages of
multiscale neural networks. Meanwhile, the theory behind amplitude separation
also needs further studies. The response operator $\mathcal{R}$ function acts
like a lowpass filter and is linear. For nonlinear operators such as the
forward operator in the inverse medium scattering problem, more research will
be done to see the performance of the multiscale DeepONet, a success for these
cases will be very helpful for solving inverse scattering problems.

##  7 Acknowledgement

The authors like to thank Prof. GE. Karniadakis for bringing the attention of
this research project to our attention and Dr. Kamaljyoti Nath for helpful
discussions and assistance during this work.

## References

  * [1] W. N. E, B. Yu, The deep Ritz method: A deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1):112, 2018.
  * [2] J. Han, A. Jentzen, W.N. E, Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences. 2018 Aug 21;115(34):8505-10.
  * [3] X. Jin, S. Cai, H. Li, GE. Karniadakis, NSFnets (Navier-Stokes flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations. Journal of Computational Physics. 2021 Feb 1;426:109951.
  * [4] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019.
  * [5] W. Cai, X.G. Li, and L.Z. Liu. A phase shift deep neural network for high frequency approximation and wave problems. SIAM Journal on Scientific Computing. 2020;42(5):A3285-312.
  * [6] Z.Q. Liu, W. Cai, and Z.Q. John Xu, Multi-scale Deep Neural Network (MscaleDNN) for Solving Poisson-Boltzmann Equation in Complex Domains, arXiv:2007.11207, 2020, Communications in Computational Physics. 2020 Jun;28(5):1970-2001.
  * [7] Lu Lu, Pengzhan Jin, and George Em Karniadakis.  DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193. 2019 Oct 8.
  * [8] Tianping Chen and Hong Chen.  Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems. IEEE Transactions on Neural Networks. 1995 Jul;6(4):911-7.
  * [9] Open System for Earthquake Engineering Simulation - Home Page.
  * [10] Jun-Hong Ding, Xian-Long Jin, Yi-Zhi Guo, and Gen-Guo Li.  Numerical simulation for large-scale seismic response analysis of immersed tunnel.  Engineering structures. 2006 Aug 1;28(10):1367-77.
  * [11] Anil K. Chopra.  Dynamics of Structures.  Pearson, 4th edition edition.
  * [12] Z.Y. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895, 2021.
  * [13] V. L. Nithin, S. Das, and H. B. Kaushik.  Wavelet-based simulation of scenario-specific nonstationary accelerograms and their GMPE compatibility.  Soil Dynamics and Earthquake Engineering. 2017 Aug 1;99:56-67.
  * [14] Julius O. Smith.  Mathematics of the Discrete Fourier Transform (DFT).  http://ccrma.stanford.edu/ jos/mdft/, accessed (date accessed).
  * [15] Alan V. Oppenheim, Alan S. Willsky, and S. Hamid Nawab.  Signals & Systems (2nd Ed.).  Prentice-Hall, Inc.


# (Spectral) Chebyshev collocation methods for solving differential equations

Pierluigi Amodio Dipartimento di Matematica, Università di Bari, Italy,
<EMAIL_ADDRESS>Luigi Brugnano Dipartimento di Matematica e
Informatica “U. Dini”, Università di Firenze, Italy<EMAIL_ADDRESS>Felice Iavernaro Dipartimento di Matematica, Università di Bari, Italy,
<EMAIL_ADDRESS>
###### Abstract

Recently, the efficient numerical solution of Hamiltonian problems has been
tackled by defining the class of energy-conserving Runge-Kutta methods named
Hamiltonian Boundary Value Methods (HBVMs). Their derivation relies on the
expansion of the vector field along the Legendre orthonormal basis.
Interestingly, this approach can be extended to cope with other orthonormal
bases and, in particular, we here consider the case of the Chebyshev
polynomial basis. The corresponding Runge-Kutta methods were previously
obtained by Costabile and Napoli [33]. In this paper, the use of a different
framework allows us to carry out a novel analysis of the methods also when
they are used as spectral formulae in time, along with some generalizations of
the methods.

Keywords: Legendre polynomials, Hamiltonian Boundary Value Methods, HBVMs,
Chebyshev polynomials, Chebyshev collocation methods.

AMS-MSC: 65L06, 65L05.

##  1 Introduction

The numerical solution of Hamiltonian problems has been the subject of many
investigations in the last decades, due to the fact that Hamiltonian problems
are not structurally stable against generic perturbations, like those induced
by a general-purpose numerical method used to approximate their solutions. The
main features of a canonical Hamiltonian system are surely the symplecticness
of the map and the conservation of energy, which cannot be simultaneously
maintained by any given Runge-Kutta or B-series method [32]. Consequently,
numerical methods have been devised in order of either defining a symplectic
discrete map, giving rise to the class of symplectic methods (see, e.g., the
monographs [46, 42, 37, 9], and references therein, the review paper [45], and
related approaches [30]), or being able to conserve the energy, resulting in
the class of energy-conserving methods (see, e.g., [44, 39, 40, 41, 31, 19,
20, 36, 18, 43] and the monographs [15]). In particular, Hamiltonian Boundary
Value Methods (HBVMs) [19, 20, 18, 15, 22] (as well as their generalizations,
see, e.g. [10, 14, 24, 29, 1, 12, 28, 26, 6]) form a special class of Runge-
Kutta methods with a (generally) rank-deficient coefficient matrix. HBVMs, in
turn, can be also viewed as the outcome obtained after a local projection of
the vector field onto a finite-dimensional function space: in particular, the
set of polynomials of a given degree. For this purpose, the Legendre
orthonormal polynomial basis has been considered so far [23] (see also [2]),
and their use as spectral methods in time has been also investigated both
theoretically and numerically [28, 17, 13, 5, 7, 8]. Remarkably, as was
already observed in [23], this idea is even more general and can be adapted to
other finite-dimensional function spaces and/or different bases. Following
this route, in this paper we consider a class of Runge-Kutta methods based on
the use of the Chebyshev polynomial basis.

It turns out that, with this choice, our approach finally leads us back to the
same formulae introduced by Costabile and Napoli by considering the classical
collocation conditions based on the Chebyshev abscissae [33]. Exploiting the
Routh-Hurwitz criterion, in [34] they derived the A-stability property of the
formule up to order 20 while, more recently, they also extended the methods to
cope with $k$-th order problems [35].

In this context, we are mainly interested in a spectral-time implementation of
these formulae, which means that the given approximation accuracy (usually
close to the machine epsilon) is achieved by increasing the order of the
formulae rather than reducing the integration stepsize. An interesting
consequence of such a strategy is that, modulo the effects of round-off
errors, the numerical solution will mimic the exact solution so that, when
applied to a Hamiltonian system the method will inherit both the simplecticity
and the conservation properties.

The advantage of using the Chebyshev basis stems from the fact that all the
entries in the Butcher tableau of the corresponding Runge-Kutta methods can be
given in closed form, thus avoiding the introduction of round-off errors when
numerically computing them (as is the case with the Legendre basis, where the
Gauss-Legendre nodes need to be numerically computed). In this respect, not
only does the analysis provided within the new framework help in deriving the
explicit expression of the coefficients of the methods for any arbitrary high
order, but it is also very useful to discuss the convergence rate when they
are used as spectral methods in time and to derive an efficient implementation
strategy based on the discrete cosine transform: these latter aspects are new,
at the best of our knowledge. Further, also a generalization of the methods is
sketched, similar to that characterizing HBVMs.

With this premise, the structure of the paper is as follows: in Section 2 we
provide a brief introduction to the framework used to derive HBVMs, relying on
the use of the Legendre polynomial basis; in Section 3 the approach is
extended to cope with the Chebyshev polynomial basis; in Section 4 a thorough
analysis of the resulting methods is given, also when they are used as
spectral methods in time; in Section 5 we provide some numerical tests, to
assess the theoretical findings; at last a few conclusions are given in
Section 6.

##  2 Hamiltonian Boundary Value Methods (HBVMs)

In order to introduce HBVMs, let us consider a canonical Hamiltonian problem
in the form

$\dot{y}=J\nabla H(y),\qquad y(0)=y_{0}\in\mathbb{R}^{2m},\qquad
J=\left(\begin{array}[]{cc}&I_{m}\\\ -I_{m}\end{array}\right)=-J^{\top},$ (1)

with $H:\mathbb{R}^{2m}\rightarrow\mathbb{R}$ the Hamiltonian function (or
energy) of the system.111In fact, for isolated mechanical systems, $H$ has the
physical meaning of the total energy. As is easily understood $H$ is a
constant of motion, since

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

being $J$ skew-symmetric. The simple idea, on which HBVMs rely on, is that of
reformulating the previous conservation property in terms of the line integral
of $\nabla H$ along the path defined by the solution $y(t)$:

$H(y(t))=H(y_{0})+\int_{0}^{t}\nabla
H(y(\tau))^{\top}\dot{y}(\tau)\mathrm{d}\tau.$

Clearly, the integral at the right-hand side of the previous equality
vanishes, because the integrand is identically zero by virtue of (1), so that
the conservation holds true for all $t>0$. The solution of (1) is the unique
function satisfying such a conservation property for all $t>0$. However, if we
consider a discrete-time dynamics, ruled by a timestep $h$, there exist
infinitely many paths $\sigma$ such that:

$\displaystyle\sigma(0)$ $\displaystyle=$ $\displaystyle
y_{0},\qquad\sigma(h)=:y_{1}\approx y(h),$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle h\int_{0}^{1}\nabla
H(\sigma(ch))^{\top}\dot{\sigma}(ch)\mathrm{d}c.$ (2)

The path $\sigma$ obviously defines a one-step numerical method that conserves
the energy, since

$H(y_{1})~{}=~{}H(y_{0})+h\int_{0}^{1}\nabla
H(\sigma(ch))^{\top}\dot{\sigma}(ch)\mathrm{d}c~{}=~{}H(y_{0}),$

even though now the integrand is no more identically zero. The methods derived
in this framework have been named line integral methods, due to the fact that
the path $\sigma$ is defined so that the line integral in (2) vanishes. Line
integral methods have been thoroughly analyzed in the monograph [15] (see also
the review paper [16]). Clearly, in the practical implementation of the
methods, the integral is replaced by a quadrature of enough high-order, thus
providing a fully discrete method, even though we shall not consider, for the
moment, this aspect, for which the reader may refer to the above mentioned
references.222This amounts to study HBVMs as continuous-stage Runge-Kutta
methods, as it has been done in [2, 3].

Interestingly enough, after some initial attempts to derive methods in this
class [39, 40, 41, 19, 25], a systematic way for their derivation was found in
[23], which is based on a local Fourier expansion of the vector field in (1).
In fact, by setting

$f(y)=J\nabla H(y),$ (3)

and using hereafter, depending on the needs, either one or the other notation,
one may rewrite problem (1), on the interval $[0,h]$, as:

$\dot{y}(ch)=\sum_{j\geq 0}P_{j}(c)\gamma_{j}(y),\qquad c\in[0,1],\qquad
y(0)=y_{0},$ (4)

where $\\{P_{j}\\}_{j\geq 0}$ is the Legendre orthonormal polynomial basis on
$[0,1]$,

$P_{i}\in\Pi_{i},\qquad\int_{0}^{1}P_{i}(\tau)P_{j}(\tau)\mathrm{d}\tau=\delta_{ij},\qquad
i,j\geq 0,$ (5)

$\Pi_{i}$, hereafter, denotes the vector space of polynomials of degree at
most $i$, $\delta_{ij}$ is the Kronecker symbol, and

$\gamma_{j}(y)=\int_{0}^{1}P_{j}(\tau)f(y(\tau h))\mathrm{d}\tau,\qquad j\geq
0,$ (6)

are the corresponding Fourier coefficients. The solution of the problem is
formally obtained, in terms of the unknown Fourier coefficients, by
integrating both sides of Equation (4):

$y(ch)=y_{0}+h\sum_{j\geq
0}\int_{0}^{c}P_{j}(x)\mathrm{d}x\,\gamma_{j}(y),\qquad c\in[0,1].$ (7)

A polynomial approximation $\sigma\in\Pi_{s}$ can be formally obtained by
truncating the previous series to finite sums:

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

and

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

respectively, with $\gamma_{j}(\sigma)$ defined according to (6), upon
replacing $y$ with $\sigma$. Whichever the degree $s\geq 1$ of the polynomial
approximation, the following result holds true, where the same notation used
above holds.

###### Theorem 1

$H(y_{1})=H(y_{0})$.

Proof In fact, one has:

$\displaystyle
H(y_{1})-H(y_{0})~{}=~{}H(\sigma(h))-H(\sigma(0))~{}=~{}h\int_{0}^{1}\nabla
H(\sigma(ch))^{\top}\dot{\sigma}(ch)\mathrm{d}c$ $\displaystyle=$
$\displaystyle h\int_{0}^{1}\nabla
H(\sigma(ch))^{\top}\sum_{j=0}^{s-1}P_{j}(c)\gamma_{j}(\sigma)\mathrm{d}c~{}=~{}h\sum_{j=0}^{s-1}\left[\int_{0}^{1}P_{j}(c)\nabla
H(\sigma(ch))\mathrm{d}c\right]^{\top}\gamma_{j}(\sigma)$ $\displaystyle=$
$\displaystyle h\sum_{j=0}^{s-1}\left[\int_{0}^{1}P_{j}(c)\nabla
H(\sigma(ch))\mathrm{d}c\right]^{\top}J\left[\int_{0}^{1}P_{j}(c)\nabla
H(\sigma(ch))\mathrm{d}c\right]~{}=~{}0,$

due to the fact that $J$ is skew-symmetric.  $\Box{~{}}$

The next result states that the method has order $2s$.

###### Theorem 2

$y_{1}-y(h)=O(h^{2s+1})$.

Proof See [23, Theorem 1].  $\Box{~{}}$

It must be emphasized that, when the method is used as a spectral method in
time, then the concept of order does not hold anymore. Instead, the following
result can be proved, under regularity assumptions on $f$ in (3).

###### Theorem 3

Let us assume $f(\sigma(t))$ be analytical in a closed ball of radius $r^{*}$
centered at $0$, then for all $h\in(0,h^{*}]$, with $h^{*}<r^{*}$, there exist
$M=M(h^{*})>0$ and $\rho>1$, $\rho\sim h^{-1}$, such that: 333Hereafter,
$|\cdot|$ will devote any convenient vector norm.

$\|\sigma-y\|~{}:=~{}\max_{c\in[0,1]}|\sigma(ch)-y(ch)|~{}\leq~{}hM\rho^{-s}.$
(10)

Proof See [5, Theorem 2].  $\Box{~{}}$

###### Remark 1

When using HBVMs as spectral methods in time, it is more important considering
the measure of the error (10), rather than the error at $h$. In fact, the
timestep used, in such a case, may be large, and even huge, so that the whole
approximation in the interval $(0,h)$ is needed.

###### Remark 2

As previously stated, an important issue, when using numerical methods as
spectral methods, stems from the fact that one needs the relevant coefficients
be computed for very high-order formulae. When some of such coefficients are
evaluated numerically, as is the case for the abscissae of the Gauss-Legendre
formulae, this may introduce errors that, even though small, may affect the
accuracy of the resulting method. It is then useful to have formulae for which
all the involved coefficients are known in closed form, and this motivates the
present paper.

##  3 Chebyshev-Runge-Kutta methods

An interesting way to interpret the approximation (8) is that of looking for
the coefficients $\gamma_{0},\dots,\gamma_{s-1}$, in the polynomial
approximation

$\dot{\sigma}(ch)=\sum_{j=0}^{s-1}P_{j}(c)\gamma_{j},\qquad
c\in[0,1],\qquad\sigma(0)=y_{0}\in\mathbb{R}^{m},$ (11)

to 444For sake of simplicity, hereafter we shall assume $f$ be an analytic
function.

$\dot{y}(ch)=f(y(ch)),\qquad c\in[0,1],\qquad y(0)=y_{0}\in\mathbb{R}^{m},$
(12)

such that the residual function

$r(ch)~{}:=~{}\dot{\sigma}(ch)-f(\sigma(ch)),\qquad c\in[0,1],$ (13)

be orthogonal to $\Pi_{s-1}$. That is,

$0~{}=~{}\int_{0}^{1}P_{i}(c)r(ch)\mathrm{d}c~{}=~{}\int_{0}^{1}P_{i}(c)\left(\dot{\sigma}(ch)-f(\sigma(ch))\right)\mathrm{d}c,\qquad
i=0,\dots,s-1.$ (14)

By considering the orthogonality relations (5), this amounts to require that

$\gamma_{i}=\int_{0}^{1}P_{i}(c)f\left(y_{0}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{i}(x)\mathrm{d}x\gamma_{j}\right)\mathrm{d}c,\qquad
i=0,\dots,s-1,$ (15)

namely (6) with $\sigma$ replacing $y$, and the new approximation given by

$y_{1}~{}:=~{}\sigma(h)~{}\equiv~{}y_{0}+h\sum_{j=0}^{s-1}\int_{0}^{1}P_{j}(x)\mathrm{d}x\gamma_{j}.$
(16)

Approximating the integrals in (15) by a Gauss-Legendre quadrature of order
$2k$, with $k\geq s$, then provides a HBVM$(k,s)$ method [15, 16, 20, 23].

A generalization of the orthogonality requirement (15) consists in considering
a suitable weighting function

$\omega(c)\geq 0,\qquad c\in[0,1],\qquad\int_{0}^{1}\omega(c)\mathrm{d}c=1,$
(17)

and a polynomial basis orthonormal w.r.t. the induced product,

$P_{i}\in\Pi_{i},\qquad\int_{0}^{1}\omega(\tau)P_{i}(\tau)P_{j}(\tau)\mathrm{d}\tau=\delta_{ij},\qquad
i,j\geq 0,$ (18)

then requiring

$0~{}=~{}\int_{0}^{1}\omega(c)P_{i}(c)r(ch)\mathrm{d}c~{}\equiv~{}\int_{0}^{1}\omega(c)P_{i}(c)\left(\dot{\sigma}(ch)-f(\sigma(ch))\right)\mathrm{d}c,\qquad
i=0,\dots,s-1.$ (19)

Consequently, now the coefficients of the polynomial approximation (11) turn
out to solve the following set of equations,

$\gamma_{i}=\int_{0}^{1}\omega(c)P_{i}(c)f\left(y_{0}+h\sum_{j=0}^{s-1}\int_{0}^{c}P_{i}(x)\mathrm{d}x\gamma_{j}\right)\mathrm{d}c,\qquad
i=0,\dots,s-1,$ (20)

in place of (15).555Clearly, (15) is derived by considering the special case
$\omega(c)\equiv 1$. In so doing, we obtain a polynomial approximation
formally still given by (9), with the Fourier coefficients now defined as

$\gamma_{j}~{}\equiv~{}\gamma_{i}(\sigma)~{}:=~{}\int_{0}^{1}\omega(c)P_{i}(c)f(\sigma(ch))\mathrm{d}c,\qquad
i=0,\dots,s-1.$ (21)

Hereafter, we shall consider the weighting function

$\omega(c)=\frac{1}{\pi\sqrt{c(1-c)}},\qquad c\in(0,1),$ (22)

which satisfies (17) and provides the shifted and scaled Chebyshev polynomials
of the first kind, i.e.,666As is usual, hereafter, $T_{j}(x)$, $j\geq 0$,
denote the Chebyshev polynomials of the first kind.

$P_{0}(c)~{}=~{}T_{0}(2c-1)~{}\equiv~{}1,\qquad
P_{j}(c)~{}=~{}\sqrt{2}\,T_{j}(2c-1),\quad j\geq 1,\qquad c\in[0,1],$ (23)

satisfying (18).

For later use, we also report the relations between such polynomials and their
integrals:

$\displaystyle\int_{0}^{c}P_{0}(x)\mathrm{d}x$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\frac{P_{1}(c)}{\sqrt{2}}+P_{0}(c)\right],$
$\displaystyle\int_{0}^{c}P_{1}(x)\mathrm{d}x$ $\displaystyle=$
$\displaystyle\frac{1}{8}\left[P_{2}(c)-\sqrt{2}P_{0}(c)\right],$ (24)
$\displaystyle\int_{0}^{c}P_{j}(x)\mathrm{d}x$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left[\frac{P_{j+1}(c)}{j+1}-\frac{P_{j-1}(c)}{j-1}-\frac{(-1)^{j}2\sqrt{2}P_{0}(c)}{j^{2}-1}\right],\qquad
j\geq 2.$

In particular, for $c=1$ one obtains, by considering that $P_{0}(1)=1$ and
$P_{j}(1)=\sqrt{2}$, for all $j\geq 1$:

$\int_{0}^{1}P_{j}(x)\mathrm{d}x=\left\\{\begin{array}[]{cl}1,&~{}j~{}=~{}0,\\\\[5.69054pt]
0,&~{}j\quad\mbox{odd,}\\\\[5.69054pt]
\frac{\sqrt{2}}{1-j^{2}},&~{}j\quad\mbox{even},\quad j\geq
2.\end{array}\right.$ (25)

###  3.1 Discretization

As is clear, the integrals involved in (20)-(21) cannot be computed exactly,
and need to be approximated by using a numerical method: this latter can be
naturally chosen as the Gauss-Chebysvev interpolatory quadrature formula of
order $2s$ on the interval $[0,1]$, with nodes

$c_{i}=\frac{1+\cos\theta_{i}}{2},\qquad\theta_{i}=\frac{2i-1}{2s}\pi,\qquad
i=1,\dots,s,$ (26)

and weights

$\omega_{i}=\frac{1}{s},\qquad i=1,\dots,s.$

We recall that $P_{s}(c_{i})=0$, $i=1,\dots,s$, due to fact that (see (23)),
$x_{i}=\cos\theta_{i}$, $i=1,\dots,s$ are the roots of $T_{s}(x)$.
Consequently, the Fourier coefficients (20)-(21) will be approximated as:

$\hat{\gamma}_{i}=\frac{1}{s}\sum_{j=1}^{s}P_{i}(c_{j})f\left(y_{0}+h\sum_{\ell=0}^{s-1}\int_{0}^{c_{j}}P_{\ell}(x)\mathrm{d}x\,\hat{\gamma}_{\ell}\right)\mathrm{d}c,\qquad
i=0,\dots,s-1.$ (27)

In so doing, we obtain a new polynomial approximation,

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

in place of $\sigma$. Setting $Y_{j}$ the argument of $f$ in Equation (27),
and taking into account (28), we arrive at the $s$-stage Runge-Kutta method
with stages

$Y_{i}~{}:=~{}u(c_{i}h)~{}=~{}y_{0}+h\sum_{j=1}^{s}\underbrace{\left(\frac{1}{s}\sum_{\ell=0}^{s-1}\int_{0}^{c_{i}}P_{\ell}(x)\mathrm{d}xP_{\ell}(c_{j})\right)}_{=:\,a_{ij}}f(Y_{j}),\qquad
i=1,\dots,s,$ (29)

with the new approximation given by:

$y_{1}~{}:=~{}u(h)=y_{0}+h\sum_{i=1}^{s}\underbrace{\left(\frac{1}{s}\sum_{\ell=0}^{s-1}\int_{0}^{1}P_{\ell}(x)\mathrm{d}xP_{\ell}(c_{i})\right)}_{=:\,b_{i}}f(Y_{i}).$
(30)

Consequently, the abscissae and weights $(c_{i},b_{i})$, $i=1,\dots,s$, and
the entries $a_{ij}$, $i,j=1,\dots,s$, define an $s$-stage Runge-Kutta method.
An explicit expression of the abscissae has been given in (26). Let us now
derive more refined expressions for the weights and the Butcher matrix. For
this purpose, let us define the following matrices:

${\cal P}_{s}=\Big{(}P_{j-1}(c_{i})\Big{)}_{i,j=1,\dots,s},\qquad{\cal
I}_{s}=\left(\int_{0}^{c_{i}}P_{j-1}(x)\mathrm{d}x\right)_{i,j=1,\dots,s},\qquad\Omega=\frac{1}{s}I_{s},\quad\in~{}\mathbb{R}^{s\times
s},$ (31)

with $I_{s}$ the identity matrix, and

$X_{s}=\left(\begin{array}[]{ccccccc}\frac{1}{2}&-\sqrt{2}\beta_{2}&\alpha_{3}&\alpha_{4}&\dots&\alpha_{s-1}&\alpha_{s}\\\\[2.84526pt]
\sqrt{2}\beta_{1}&0&-\beta_{1}&\\\ &\beta_{2}&0&-\beta_{2}&&O\\\
&&\beta_{3}&0&-\beta_{3}\\\ &&&\ddots&\ddots&\ddots\\\
&O&&&\beta_{s-2}&0&-\beta_{s-2}\\\
&&&&&\beta_{s-1}&0\end{array}\right)\in\mathbb{R}^{s\times s},$ (32)

with

$\beta_{j}=\frac{1}{4j},\quad j\geq
1,\qquad\alpha_{j}=(-1)^{j}8\sqrt{2}\beta_{j}\beta_{j-2},\quad j\geq 3.$ (33)

The following results hold true.

###### Lemma 1

With reference to (26) and the matrices (31) and (32), one has:

$({\cal
P}_{s})_{ij}=\left[\sqrt{2}+\delta_{j1}(1-\sqrt{2})\right]\cos(j-1)\theta_{i},\qquad
i,j=1,\dots,s,$ ${\cal I}_{s}={\cal P}_{s}X_{s},\qquad P_{s}^{\top}\Omega{\cal
P}_{s}=I_{s}.$

Proof The first statement follows from (23) and the well-known fact that

$T_{j-1}(\cos\theta_{i})=\cos(j-1)\theta_{i}.$

The second statement follows from (24), by considering that $P_{s}(c_{i})=0$,
$i=1,\dots,s$. Finally, one has, by setting, hereafter,
$e_{i}\in\mathbb{R}^{s}$ the $i$th unit vector:

$\displaystyle e_{i}^{\top}\left({\cal P}_{s}^{\top}\Omega{\cal
P}_{s}\right)e_{j}$ $\displaystyle=$ $\displaystyle\frac{1}{s}({\cal
P}_{s}e_{i})^{\top}({\cal P}_{s}e_{j})$ $\displaystyle=$
$\displaystyle\frac{\left[\sqrt{2}+\delta_{i1}(1-\sqrt{2})\right]\left[\sqrt{2}+\delta_{j1}(1-\sqrt{2})\right]}{s}\sum_{k=1}^{s}\cos(i-1)\theta_{k}\cdot\cos(j-1)\theta_{k}$
$\displaystyle=$ $\displaystyle\delta_{ij},$

due to the fact that, for all $i,j=1,\dots,s$,

$\sum_{k=1}^{s}\cos(i-1)\theta_{k}\cdot\cos(j-1)\theta_{k}~{}=~{}\left\\{\begin{array}[]{cl}0,&i\neq
j,\\\\[2.84526pt] s,&i=j=1,\\\\[2.84526pt] \frac{s}{2},&i=j\geq
1.\,\mbox{~{}$\Box{~{}}$}\end{array}\right.$

Consequently, we can now state the following result.

###### Theorem 4

The Butcher matrix of the Runge-Kutta method (29)-(30) is given by

$A~{}\equiv~{}(a_{ij})~{}:=~{}{\cal I}_{s}{\cal
P}_{s}^{\top}\Omega~{}\equiv~{}{\cal P}_{s}X_{s}{\cal
P}_{s}^{\top}\Omega~{}\equiv~{}{\cal P}_{s}X_{s}{\cal P}_{s}^{-1}.$ (34)

The corresponding weights are given by:

$b_{i}:=\frac{1}{s}\left[1-2\sum_{j=1}^{\lceil
s/2\rceil-1}\frac{\cos\frac{2i-1}{s}j\pi}{4j^{2}-1}\right],\qquad
i=1,\dots,s.$ (35)

Proof From Lemma 1, the last two equalities in (34) easily follow. The first
equality in (34) then follows from (29), by taking into account (31):

$a_{ij}~{}=~{}e_{i}^{\top}Ae_{j}~{}=~{}e_{i}^{\top}{\cal I}_{s}{\cal
P}_{s}^{\top}\Omega e_{j}~{}=~{}\frac{1}{s}(e_{i}^{\top}{\cal
I}_{s})(e_{j}^{\top}{\cal
P}_{s})^{\top}~{}=~{}\frac{1}{s}\sum_{\ell=0}^{s-1}\int_{0}^{c_{i}}P_{\ell}(x)\mathrm{d}xP_{\ell}(c_{j}),$

which coincides with the formula given in (29). The formula (35) of the
weights follows by considering that from (30), taking into account (25) and
(26), one has:

$\displaystyle b_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{s}\sum_{\ell=0}^{s-1}\int_{0}^{1}P_{\ell}(x)\mathrm{d}xP_{\ell}(c_{i})~{}=~{}\left(\int_{0}^{1}P_{0}(c)\mathrm{d}c\,,\dots,\int_{0}^{1}P_{s-1}(c)\mathrm{d}c\right)^{\top}{\cal
P}_{s}^{\top}\Omega e_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{s}\left(1,\,0,\,\frac{-\sqrt{2}}{2^{3}-1},\,0,\,\frac{-\sqrt{2}}{4^{3}-1},\,0,\,\frac{-\sqrt{2}}{6^{3}-1},\,\dots\right)\left(\begin{array}[]{c}1\\\
\sqrt{2}\cos\theta_{i}\\\ \vdots\\\
\sqrt{2}\cos(s-1)\theta_{i}\end{array}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{s}\left[1-2\sum_{j=1}^{\lceil s/2\rceil-1}\frac{\cos
2j\theta_{i}}{4j^{2}-1}\right]~{}=~{}\frac{1}{s}\left[1-2\sum_{j=1}^{\lceil
s/2\rceil-1}\frac{\cos\frac{2i-1}{s}j\pi}{4j^{2}-1}\right].\,\mbox{~{}$\Box{~{}}$}$

Moreover, concerning the weights of the method, the following result holds
true.

###### Theorem 5

The weights $\\{b_{i}\\}_{i=1,\dots,s}$ defined in (35) are all positive.

Proof In fact, from (35), having fixed a given value of $s\geq 1$, one has:

$\displaystyle b_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{s}\left[1-2\sum_{j=1}^{\lceil
s/2\rceil-1}\frac{\cos\frac{2i-1}{s}j\pi}{4j^{2}-1}\right]~{}\geq~{}\frac{1}{s}\left[1-\sum_{j=1}^{\lceil
s/2\rceil-1}\frac{2}{4j^{2}-1}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{s}\left[1-\sum_{j=1}^{\lceil
s/2\rceil-1}\left(\frac{1}{2j-1}-\frac{1}{2j+1}\right)\right]~{}=~{}\frac{1}{s}\left[1-1+\frac{1}{2(\lceil
s/2\rceil-1)+1}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{s}\,\frac{1}{2\lceil
s/2\rceil-1}~{}\geq~{}\frac{1}{s}\,\frac{1}{s}~{}=~{}\frac{1}{s^{2}}.\,\mbox{~{}$\Box{~{}}$}$

In conclusion, we have derived the family of $s$-stage Runge-Kutta method
whose Butcher tableau, with reference to (26), (35), and (31)–(33), is given
by:

$\begin{array}[]{c|c}{\bm{c}}&{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\\\
\hline\cr\\\\[-8.53581pt]
&{\bm{b}}^{\top}\end{array}\quad\equiv\quad\begin{array}[]{c|c}{\bm{c}}&{\cal
P}_{s}X_{s}{\cal P}_{s}^{\top}\Omega\\\ \hline\cr\\\\[-8.53581pt]
&{\bm{b}}^{\top}\end{array}\quad\equiv\quad\begin{array}[]{c|c}{\bm{c}}&{\cal
P}_{s}X_{s}{\cal P}_{s}^{-1}\\\ \hline\cr\\\\[-8.53581pt]
&{\bm{b}}^{\top}\end{array}$ (37)

having set, as is usual,

${\bm{c}}=(c_{1},\,\dots,\,c_{s})^{\top},\qquad{\bm{b}}=(b_{1},\,\dots,\,b_{s})^{\top}.$
(38)

In particular, the last form in (37) can be regarded as a kind of
$W$-transformation [38] of the method.

###### Remark 3

The Runge-Kutta methods in (37)-(38) coincide with the Chebyshev Collocation
Methods derived by Costabile and Napoli in [33]. For sake of brevity, we shall
refer to the $s$-stage method as CCM$(s)$.

###  3.2 Collocation methods

As anticipated above, CCM$(s)$ methods are collocation methods. In fact, from
(27), (29), and (37), by setting

$Y~{}=~{}\left(\begin{array}[]{c}u(c_{1}h)\\\ \vdots\\\
u(c_{s}h)\end{array}\right)~{}\equiv~{}u({\bm{c}}h),\qquad\hat{{\bm{\gamma}}}=\left(\begin{array}[]{c}\hat{\gamma}_{0}\\\
\vdots\\\ \hat{\gamma}_{s-1}\end{array}\right),$

the vector of the stages and of the Fourier coefficients, respectively, one
derives that:

$\dot{u}({\bm{c}}h)={\cal P}_{s}\otimes
I_{m}\hat{{\bm{\gamma}}},\qquad\hat{{\bm{\gamma}}}={\cal
P}_{s}^{\top}\Omega\otimes I_{m}f(u({\bm{c}}h)).$

By combining the last two equations, and taking into account that, by virtue
of Lemma 1, ${\cal P}_{s}{\cal P}_{s}^{\top}\Omega=I_{s}$, one eventually
derives that

$\dot{u}({\bm{c}}h)=f({\bm{c}}h),$

i.e.,

$\dot{u}(c_{i}h)=f(u(c_{i}h)),\qquad i=1,\dots,s.$

Further, the collocation polynomial $u$ satisfies $u(0)=y_{0}$ and
$u(h)=:y_{1}$, as seen above.

###  3.3 Interesting implementation details

An interesting property of the Butcher matrix in the tableau (37) derives from
the fact that the multiplications

$\sqrt{s}\,{\cal P}_{s}^{\top}\Omega\otimes
I_{m}V\equiv\frac{1}{\sqrt{s}}{\cal P}_{s}^{\top}\otimes
I_{m}V,\qquad\mbox{and}\qquad\frac{1}{\sqrt{s}}{\cal P}_{s}\otimes I_{m}Z,$

with $V,Z\in\mathbb{R}^{sm}$ given vectors, amount to the discrete cosine
transform of $V$, ${\tt dct}(V)$, and the inverse discrete cosine transform of
$Z$, ${\tt idct}(Z)$, respectively.777Here, we have used the name of the
Matlab© functions implementing the two transformations. Consequently, a fixed-
point iteration for computing the stages of the CCM$(s)$ method reads, by
setting hereafter the vector
${\bm{e}}=\left(1,\,1,\,\dots,\,1\right)^{\top}\in\mathbb{R}^{s}$:

$Y^{(r+1)}={\bm{e}}\otimes y_{0}+h\,{\tt idct}\left(X_{s}\otimes I_{m}\,{\tt
dct}(f(Y^{(r)}))\right),\qquad r=0,1,\dots,$ (39)

which can be advantageous, for large values of $s$. In fact, by considering
that the matrix $X_{s}$ in (32) is sparse, the complexity of one iteration
amounts to computing the vector field in the stage values, plus $O(ms\log(s))$
flops,888As is usual, 1 flop denotes a basic floating-point operation. whereas
the standard implementation would require $2ms^{2}$ flops.

###  3.4 Generalizations

We observe that a possible generalizations of the collocation methods
described above is that of using, in the discretization procedure of the
integrals involved in (20)-(21), a Gauss-Chebysvev interpolatory quadrature
formula of order $2k$ on the interval $[0,1]$, with nodes and weights

$c_{i}=\frac{1+\cos\theta_{i}}{2},\qquad\theta_{i}=\frac{2i-1}{2k}\pi,\qquad\omega_{i}=\frac{1}{k},\qquad
i=1,\dots,k,$ (40)

for a convenient $k\geq s$. In such a case, for $k>s$, the methods are no more
collocation methods. This is analogous to what happens for HBVM$(k,s)$ methods
[15] that, when $k=s$, reduce to the $s$-stage Gauss-Legendre collocation
methods but, for $k>s$, are no more a collocation methods. By choosing the
absissae (40), the matrices (31)-(32) respectively become:

${\cal
P}_{s}=\Big{(}P_{j-1}(c_{i})\Big{)}_{\scriptsize\begin{array}[]{l}i=1,\dots,k\\\
j=1,\dots,s\end{array}},~{}{\cal
I}_{s}=\left(\int_{0}^{c_{i}}P_{j-1}(x)\mathrm{d}x\right)_{\scriptsize\begin{array}[]{l}i=1,\dots,k\\\
j=1,\dots,s\end{array}}\in\mathbb{R}^{k\times
s},\qquad\Omega=\frac{1}{k}I_{k}\in\mathbb{R}^{k\times k},$ (41)

and (see (32))

$\hat{X}_{s}=\left(\begin{array}[]{ccccccc}\frac{1}{2}&-\sqrt{2}\beta_{2}&\alpha_{3}&\alpha_{4}&\dots&\alpha_{s-1}&\alpha_{s}\\\\[2.84526pt]
\sqrt{2}\beta_{1}&0&-\beta_{1}&\\\ &\beta_{2}&0&-\beta_{2}&&O\\\
&&\beta_{3}&0&-\beta_{3}\\\ &&&\ddots&\ddots&\ddots\\\
&O&&&\beta_{s-2}&0&-\beta_{s-2}\\\ &&&&&\beta_{s-1}&0\\\
\hline\cr&&&&&&\beta_{s}\end{array}\right)\equiv\left(\begin{array}[]{c}X_{s}\\\
\hline\cr 0,\dots,0,\beta_{s}\end{array}\right)\in\mathbb{R}^{s+1\times s}.$

This latter matrix is such that, for $k>s$:

${\cal I}_{s}={\cal P}_{s+1}\hat{X}_{s},$ (42)

where ${\cal P}_{s+1}\in\mathbb{R}^{k\times s+1}$ is defined similarly as in
(41).999Clearly, when $k=s$, (42) reduces to ${\cal I}_{s}={\cal P}_{s}X_{s}$,
according to Lemma 1. Finally, the corresponding Butcher tableau become
(compare with (37)-(38)):

$\begin{array}[]{c|c}{\bm{c}}&{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\\\
\hline\cr\\\\[-8.53581pt]
&{\bm{b}}^{\top}\end{array}\quad\equiv\quad\begin{array}[]{c|c}{\bm{c}}&{\cal
P}_{s+1}\hat{X}_{s}{\cal P}_{s}^{\top}\Omega\\\ \hline\cr\\\\[-8.53581pt]
&{\bm{b}}^{\top}\end{array}\,,\qquad{\bm{c}}=(c_{1},\,\dots,\,c_{k})^{\top},\qquad{\bm{b}}=(b_{1},\,\dots,\,b_{k})^{\top},$

with the weights given by

$b_{i}:=\frac{1}{k}\left[1-2\sum_{j=1}^{\lceil
s/2\rceil-1}\frac{\cos\frac{2i-1}{k}j\pi}{4j^{2}-1}\right],\qquad
i=1,\dots,k,$

in place of (35). Nevertheless, provided that $k\geq s$ holds, the analysis of
such methods is similar to that of CCM$(s)$ methods (all the results in
Section 4 continue formally to hold), so that we shall not consider them
further.

##  4 Analysis of the methods

In this section we carry out the analysis of the CCM$(s)$ method defined by
the Butcher tableau (37)-(38). In particular, we study:

  1. 1.

the symmetry of the method;

  2. 2.

its order of convergence;

  3. 3.

its linear stability;

  4. 4.

its accuracy when used as a spectral method in time.

The analysis of items 2 and 3 uses different arguments from those provided in
[34, 35], whereas the analysis related to items 1 and 4 is novel.

###  4.1 Symmetry

To begin with, let us recall a few symmetry properties of the abscissae (26)
and of the polynomials (23), which we state without proof (they derive from
well known properties of Chebyshev polynomials and abscissae):

$c_{i}=1-c_{s-i+1},\quad i=1,\dots,s,\qquad P_{j}(1-c)=(-1)^{j}P_{j}(c),\quad
c\in[0,1],\quad i=0,\dots,s.$ (43)

For later use, we also set the vector

${\cal
I}_{s}(1)~{}:=~{}\left(\begin{array}[]{c}\int_{0}^{1}P_{0}(x)\mathrm{d}x\\\
\vdots\\\
\int_{0}^{1}P_{s-1}(x)\mathrm{d}x\end{array}\right)~{}\equiv~{}\left(\begin{array}[]{c}1\\\
0\\\ \frac{-\sqrt{2}}{2^{3}-1}\\\ 0\\\ \frac{-\sqrt{2}}{4^{3}-1}\\\ 0\\\
\frac{-\sqrt{2}}{6^{3}-1}\\\ \vdots\end{array}\right),$ (44)

with the last equality following from (25), such that, according to what seen
in the proof of Theorem 4:

${\bm{b}}^{\top}={\cal I}_{s}(1)^{\top}{\cal P}_{s}^{\top}\Omega.$ (45)

We also need to define the following matrices:

$P=\left(\begin{array}[]{ccc}&&1\\\ &\reflectbox{$\ddots$}\\\
1\end{array}\right),\quad D=\left(\begin{array}[]{cccccc}1\\\ &-1\\\ &&1\\\
&&&-1\\\ &&&&\ddots\\\ &&&&&(-1)^{s-1}\end{array}\right)\in\mathbb{R}^{s\times
s}.$ (46)

We can now state the following result.

###### Theorem 6

The method (37)-(38) is symmetric.

Proof Since $P{\bm{c}}={\bm{c}}$, see (43), it is known that the symmetry of
the method (37)-(38) is granted, provided that [37]:

$P{\bm{b}}={\bm{b}},\qquad P\left({\cal I}_{s}{\cal
P}_{s}^{\top}\Omega\right)P={\bm{e}}\,{\bm{b}}^{\top}-{\cal I}_{s}{\cal
P}_{s}^{\top}\Omega.$ (47)

Let us start proving the first equality which, in turn, amounts to show that,
by virtue of (45):

${\bm{b}}^{\top}~{}=~{}{\cal I}_{s}(1)^{\top}{\cal P}_{s}^{\top}\Omega
P~{}=~{}{\bm{b}}^{\top}P.$

In fact, by taking into account (31), (43), and (46), one at first derives
that:

${\cal P}_{s}^{\top}\Omega P=D{\cal P}_{s}^{\top}\Omega.$ (48)

Moreover, by considering that $D{\cal I}_{s}(1)={\cal I}_{s}(1)$, since the
even entries of ${\cal I}_{s}(1)$ are zero (see (44)), one has:

${\bm{b}}^{\top}~{}=~{}{\cal I}_{s}(1)^{\top}{\cal
P}_{s}^{\top}\Omega~{}=~{}{\cal I}_{s}(1)^{\top}D{\cal
P}_{s}^{\top}\Omega~{}=~{}{\cal I}_{s}(1)^{\top}{\cal P}_{s}^{\top}\Omega
P~{}=~{}{\bm{b}}^{\top}P.$

Consequently, the first part of the statement follows. Further, one obtains,
by virtue of (48):

$P({\cal I}_{s}{\cal P}_{s}^{\top}\Omega)P=P({\cal I}_{s}D^{2}{\cal
P}_{s}^{\top}\Omega)P=(P{\cal I}_{s}D)(D{\cal P}_{s}^{\top}\Omega P)=(P{\cal
I}_{s}D){\cal P}_{s}^{\top}\Omega.$

Therefore, from (45) it follows that the second statement in (47) holds true,
provided that

$P{\cal I}_{s}D={\bm{e}}\,{\cal I}_{s}(1)^{\top}-{\cal
I}_{s}\qquad\Longleftrightarrow\qquad(P{\cal
I}_{s}D)_{ij}=\int_{c_{i}}^{1}P_{j-1}(\xi)\mathrm{d}\xi,\quad i,j=1,\dots,s.$

From (31) and (46) one obtains, by virtue of (43):

$\displaystyle(P{\cal I}_{s}D)_{ij}$ $\displaystyle=$
$\displaystyle(-1)^{j-1}\int_{0}^{c_{s-i+1}}P_{j-1}(x)\mathrm{d}x~{}=~{}\int_{0}^{1-c_{i}}(-1)^{j-1}P_{j-1}(x)\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{0}^{1-c_{i}}P_{j-1}(1-x)\mathrm{d}x~{}=~{}-\int_{1}^{c_{i}}P_{j-1}(\xi)\mathrm{d}\xi~{}=~{}\int_{c_{i}}^{1}P_{j-1}(\xi)\mathrm{d}\xi.$

Consequently, the statement is proved.  $\Box{~{}}$

###  4.2 Stability

Because of its symmetry, one concludes that for all $s\geq 1$ the absolute
region of a CCM$(s)$ method coincides with $\mathbb{C}^{-}$, the negative-real
complex plane, provided that all the eigenvalues of $X_{s}$ have positive real
part (in fact, $X_{s}$ is similar to the Butcher matrix, see (37)). We have
numerically verified that $X_{s}$ has all the eigenvalues with positive real
part for all $s\geq 1$ (actually, we stopped at $s=1000$). We can then
conclude that CCM$(s)$ methods are perfectly (or precisely) $A$-stable for all
$s\geq 1$.

###  4.3 Order of convergence

Let us now study the order of convergence of a CCM$(s)$ method. For this
purpose, we need the following preliminary results.

###### Lemma 2

The convergence order of a symmetric method is even.

Proof See [37, Theorem 3.2].  $\Box{~{}}$

###### Lemma 3

Assume that a function $G:[0,h]\rightarrow V$, with $V$ a vector space, admits
a Taylor expansion at 0. Then, with reference to the orthonormal basis
(22)-(23), one has:

$\psi_{j}~{}:=~{}\int_{0}^{1}\omega(c)P_{j}(c)G(ch)\mathrm{d}c~{}=~{}O(h^{j}),\qquad
j=0,1,\dots.$

Proof In fact, one has:

$\displaystyle\psi_{j}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\omega(c)P_{j}(c)G(ch)\mathrm{d}c~{}=~{}\int_{0}^{1}\omega(c)P_{j}(c)\sum_{\ell\geq
0}\frac{G^{(\ell)}(0)}{\ell!}(ch)^{\ell}\mathrm{d}c$ $\displaystyle=$
$\displaystyle\sum_{\ell\geq
0}\frac{G^{(\ell)}(0)}{\ell!}h^{\ell}\int_{0}^{1}\omega(c)P_{j}(c)c^{\ell}\mathrm{d}c~{}=~{}\sum_{\ell\geq
j}\frac{G^{(\ell)}(0)}{\ell!}h^{\ell}\int_{0}^{1}\omega(c)P_{j}(c)c^{\ell}\mathrm{d}c$
$\displaystyle=$ $\displaystyle O(h^{j}).\,\mbox{~{}$\Box{~{}}$}$

###### Lemma 4

With reference to the approximate Fourier coefficients defined in (27), one
has:

$\hat{\gamma}_{j}=O(h^{j}),\qquad j=0,\dots,s-1.$

Proof In fact, from the previous Lemma 3 applied to $G(ch)=f(u(ch))$, we know
that (using the notation (21))

$\gamma_{j}(u)~{}=~{}\int_{0}^{1}\omega(c)P_{j}(c)f(u(ch))\mathrm{d}c~{}=~{}O(h^{j}).$
(49)

On the other hand, the quadrature error,

$\Delta_{j}(h)~{}:=~{}\gamma_{j}(u)-\hat{\gamma}_{j}~{}=~{}O(h^{2s-j}).$ (50)

Consequently,

$\hat{\gamma}_{j}~{}=~{}\gamma_{j}(u)-\Delta_{j}(h)~{}=~{}O(h^{j}),\qquad
j=0,\dots,s-1.\mbox{~{}$\Box{~{}}$}$

Finally, we need to recall some well-known perturbation results for ODE-IVPs.
For this purpose, let us denote $y(t,\xi,\eta)$ the solution of the problem
(compare with (12)),

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

Then:

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

having set

$\Phi(t,\xi,\eta)~{}\equiv~{}\frac{\partial}{\partial\eta}y(t,\xi,\eta),$

the solution of the variational problem

$\dot{\Phi}(t,\xi,\eta)~{}=~{}f^{\prime}(y(t,\xi,\eta))\Phi(t,\xi,\eta),\qquad
t>\xi,\qquad\Phi(\xi,\xi,\eta)=I_{m}.$

We can now state the convergence result.

###### Theorem 7

With reference to the polynomial approximation (28) defined by the CCM$(s)$
method (37)-(38) to the solution of (12), one has:

$\|u-y\|=O(h^{s+1}),\qquad y_{1}-y(h)=O(h^{r+1}),$

with $r=s$ if $s$ is even, or $r=s+1$, otherwise.

Proof By taking into account the above arguments, and using the notation
(49)-(50), one has:

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

Consequently, the first part of the statement holds. The last part of the
statement then follows from Lemma 2, by taking $c=1$.  $\Box{~{}}$

###  4.4 Analysis as spectral method

As was pointed out in the introduction, one main reason for introducing CCMs
is their use as spectral methods in time. This strategy allows for using large
(sometimes huge) timesteps $h$. In this context, the classical notion of order
of convergence, based on the fact that $h\rightarrow 0$, does not apply
anymore, so that a different analysis is needed. We shall here follow similar
steps as those described in [5] for spectral HBVMs.

To begin with, we recall that, for the modified Chebysehv basis (22)-(23) the
following properties hold true:101010The second property has already been used
in the proof of Lemma 3.

$\forall j\geq
0:\quad\|P_{j}\|\leq\sqrt{2},\qquad\int_{0}^{1}\omega(c)P_{j}(c)c^{\ell}\mathrm{d}c=0,\quad\ell=0,\dots,j-1.$
(51)

Moreover, we recall that, with reference to the notation (21):111111Also this
property has already been used, in the proof of Theorem 7.

$f(y(ch))=\sum_{j\geq 0}P_{j}(c)\gamma_{j}(y),\qquad c\in[0,1].$ (52)

We also state the following preliminary result.

###### Lemma 5

With reference to (22) and to the polynomials (23), let us define, for a given
$\rho>1$:

$Q_{j}(\xi)~{}=~{}\int_{0}^{1}\omega(c)\frac{P_{j}(c)}{\xi-c}\mathrm{d}c,\qquad\xi\in{\cal
C}_{\rho}(0),$

having set ${\cal C}_{\rho}(0)=\\{z\in\mathbb{C}:|z|=\rho\\}$. Then

$\|Q_{j}\|_{\rho}~{}:=~{}\max_{\xi\in{\cal
C}_{\rho}(0)}|Q_{j}(\xi)|~{}\leq~{}\sqrt{2}\frac{\rho^{-j}}{(\rho-1)}.$ (53)

Proof By taking into account (51), for $|\xi|=\rho>1$ one has:

$\displaystyle Q_{j}(\xi)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\omega(c)\frac{P_{j}(c)}{\xi-c}\mathrm{d}c~{}=~{}\xi^{-1}\int_{0}^{1}\omega(c)\frac{P_{j}(c)}{1-\xi^{-1}c}\mathrm{d}c$
$\displaystyle=$
$\displaystyle\xi^{-1}\int_{0}^{1}\omega(c)P_{j}(c)\sum_{\ell\geq
0}\xi^{-\ell}c^{\ell}\mathrm{d}c~{}=~{}\xi^{-1}\sum_{\ell\geq
j}\xi^{-\ell}\int_{0}^{1}\omega(c)P_{j}(c)c^{\ell}\mathrm{d}c$
$\displaystyle=$ $\displaystyle\xi^{-j-1}\sum_{\ell\geq
0}\xi^{-\ell}\int_{0}^{1}\omega(c)P_{j}(c)c^{\ell+j}\mathrm{d}c.$

Evaluating the norms, one has, again by virtue of (51):

$\left|\int_{0}^{1}\omega(c)P_{j}(c)c^{\ell+j}\mathrm{d}c\right|~{}\leq~{}\|P_{j}\|\int_{0}^{1}\omega(c)\mathrm{d}c~{}\leq~{}\sqrt{2}.$

Consequently, one obtains:

$|Q_{j}(\xi)|~{}\leq~{}\sqrt{2}\rho^{-j-1}\sum_{\ell\geq
0}\rho^{-\ell}~{}=~{}\sqrt{2}\rho^{-j-1}\frac{1}{1-\rho^{-1}}~{}=~{}\sqrt{2}\frac{\rho^{-j}}{\rho-1}.\,\mbox{~{}$\Box{~{}}$}$

Further, by setting hereafter, for $r>0$,

${\cal B}_{r}(0)=\\{z\in\mathbb{C}:|z|\leq r\\},$

we have the following straightforward result.

###### Lemma 6

Let $g(z)$ be analytical in ${\cal B}_{r^{*}}(0)$, for a given $r^{*}>0$.
Then, for all $h\in(0,h^{*}]$, with $h^{*}<r^{*}$,

$g_{h}(\xi)~{}:=~{}g(\xi h)$ (54)

is analytical in ${\cal B}_{\rho}(0)$, with:

$\rho~{}=~{}\rho(h)~{}:=~{}\frac{r^{*}}{h}~{}\geq~{}\frac{r^{*}}{h^{*}}~{}=:~{}\rho^{*}~{}>~{}1.$
(55)

We can now state the following result.

###### Theorem 8

With reference to (12), let assume that the function

$g(z)~{}:=~{}f(y(z))$ (56)

and the timestep $h$ satisfy the hypotheses of Lemma 6. Then there exists
$\kappa=\kappa(h^{*})$ and $\rho>1$, $\rho\sim h^{-1}$, such that

$|\gamma_{j}(y)|~{}\leq~{}\kappa\rho^{-j},\qquad j\geq 0.$ (57)

Proof By taking into account (21), (52), and (54)-(55), by virtue of Lemma 5
one has:

$\displaystyle\gamma_{j}(y)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\omega(c)P_{j}(c)f(y(ch))\mathrm{d}c~{}\equiv~{}\int_{0}^{1}\omega(c)P_{j}(c)g_{h}(c)\mathrm{d}c$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\omega(c)P_{j}(c)\left[\frac{1}{2\pi i}\int_{{\cal
C}_{\rho}(0)}\frac{g_{h}(\xi)}{\xi-c}\mathrm{d}\xi\right]\mathrm{d}c~{}=~{}\frac{1}{2\pi
i}\int_{{\cal
C}_{\rho}(0)}g_{h}(\xi)\left[\int_{0}^{1}\omega(c)\frac{P_{j}(c)}{\xi-c}\mathrm{d}c\right]\mathrm{d}\xi$
$\displaystyle\equiv$ $\displaystyle\frac{1}{2\pi i}\int_{{\cal
C}_{\rho}(0)}g_{h}(\xi)Q_{j}(\xi)\mathrm{d}\xi.$

Consequently (see (53)),

$|\gamma_{j}(y)|~{}\leq~{}\frac{\rho}{2\pi}\|g_{h}\|_{\rho}\|Q_{j}\|_{\rho}.$

From Lemma 5 we know that

$\|Q_{j}\|_{\rho}~{}\leq~{}\sqrt{2}\frac{\rho^{-j}}{\rho-1}.$

Further,

$\|g_{h}\|_{\rho}~{}=~{}\max_{|z|=\rho}|g_{h}(z)|~{}\leq~{}\max_{|z|\leq\rho^{*}}|g_{h}(z)|~{}\equiv~{}\max_{|z|\leq
r^{*}}|g(z)|~{}=:~{}\|g\|.$

Consequently, one eventually derives:

$|\gamma_{j}(y)|~{}\leq~{}\frac{\rho\sqrt{2}\|g\|}{2\pi(\rho-1)}\rho^{-j}~{}\leq~{}\overbrace{\frac{\rho^{*}\sqrt{2}\|g\|}{2\pi(\rho^{*}-1)}}^{=:\,\kappa}\rho^{-j},$

from which the statement follows by taking into account that (see (55)),
$\rho^{*}=r^{*}/h^{*}$.  $\Box{~{}}$

Finally, we have the following result.

###### Theorem 9

Let us consider the polynomial approximation (9) and (22)-(23), with the
Fourier coefficients given by (21). Then, Theorem 3 continues formally to
hold.

Proof Following similar steps as in the proof of Theorem 7, one has:

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

Consequently, from (51) and Theorem 8, also considering (55), one derives:

$\displaystyle\|\sigma-y\|$ $\displaystyle\leq$ $\displaystyle
h\sqrt{2}\kappa\max_{c,\tau\in[0,1]}\|\Phi(ch,\tau h,\sigma(\tau
h))\|\sum_{j\geq s}\rho^{-j}$ $\displaystyle\leq$ $\displaystyle
h\underbrace{\sqrt{2}\kappa\max_{c,\tau\in[0,1]}\|\Phi(ch^{*},\tau
h^{*},\sigma(\tau
h^{*}))\|\frac{1}{1-(\rho^{*})^{-1}}}_{=:\,M}\rho^{-s}~{}\equiv~{}hM\rho^{-s}.\,\mbox{~{}$\Box{~{}}$}$

###### Remark 4

When a CCM$(s)$ method is used as a spectral method in time, one obviously
assume that the value of $s$ is large enough so that the quadrature errors
falls below the round-off error level. In other words, the polynomials
$\sigma$ and $u$ are undistinguishable, within the round-off error level.

##  5 Numerical tests

We here report a few numerical tests to assess the theoretical findings. They
have been carried out on a 3 GHz Intel Xeon W10 core computer with 64GB of
memory, running Matlab© 2020b. We consider the Kepler problem,

$\dot{q}_{i}=p_{i},\qquad\dot{p}_{i}=\frac{-q_{i}}{(q_{1}^{2}+q_{2}^{2})^{\frac{3}{2}}},\qquad
i=1,2,$ (58)

that, when considering the trajectory starting at

$q_{1}(0)=0.4,\qquad q_{2}(0)=p_{1}(0)=0,\qquad p_{2}(0)=2,$ (59)

has a periodic solution of period $2\pi$. In Table 2 we list the numerical
results obtained by solving this problem on one period by means of the
CCM$(s)$ method, $s=1,\dots,4$, using timestep $h=2\pi/n$, namely, the errorr
err after one period and the corresponding rate of convergence. As one may
see, the listed results agree with what stated in Theorem 7, in particular the
fact that the order of the methods is even, due to their symmetry.

Table 1: Numerical results when solving problem (58)-(59) on the interval
$[0,2\pi]$ with timestep $h=2\pi/n$.

| $s=1$ | $s=2$ | $s=3$ | $s=4$
---|---|---|---|---
$n$ | err | rate | err | rate | err | rate | err | rate
50 | 2.98e+0 | — | 2.24e+0 | — | 7.36e-03 | — | 7.33e-03 | —
100 | 1.66e+0 | 0.8 | 9.45e-01 | 1.2 | 6.15e-04 | 3.6 | 4.46e-04 | 4.0
200 | 5.23e-01 | 1.7 | 2.53e-01 | 1.9 | 4.03e-05 | 3.9 | 2.78e-05 | 4.0
400 | 1.34e-01 | 2.0 | 6.34e-02 | 2.0 | 2.55e-06 | 4.0 | 1.73e-06 | 4.0
800 | 3.35e-02 | 2.0 | 1.58e-02 | 2.0 | 1.60e-07 | 4.0 | 1.08e-07 | 4.0
1600 | 8.38e-03 | 2.0 | 3.96e-03 | 2.0 | 1.00e-08 | 4.0 | 6.77e-09 | 4.0

Table 2: Numerical results when solving problem (58)-(59) on the interval
$[0,20\pi]$ by using the CCM(50) method with timestep $h=2\pi/n$.

$n$ | 3 | 6 | 9 | 12 | 15
---|---|---|---|---|---
period | err | err | err | err | err
1 | 5.04e-12 | 9.14e-14 | 7.37e-14 | 1.27e-13 | 7.44e-14
2 | 9.72e-12 | 7.96e-14 | 1.24e-13 | 3.05e-13 | 5.72e-14
3 | 1.34e-11 | 3.54e-13 | 2.81e-13 | 7.23e-13 | 5.49e-14
4 | 1.90e-11 | 3.00e-13 | 7.80e-13 | 1.27e-12 | 1.49e-13
5 | 2.55e-11 | 4.69e-13 | 1.34e-12 | 2.30e-12 | 2.77e-13
6 | 3.04e-11 | 5.68e-13 | 1.75e-12 | 3.25e-12 | 3.59e-13
7 | 3.44e-11 | 2.13e-13 | 1.73e-12 | 4.23e-12 | 2.37e-13
8 | 3.94e-11 | 2.72e-13 | 1.45e-12 | 5.19e-12 | 2.04e-13
9 | 4.38e-11 | 6.93e-13 | 1.18e-12 | 6.15e-12 | 3.29e-13
10 | 4.77e-11 | 1.54e-12 | 8.38e-13 | 7.01e-12 | 5.00e-13

Next, we consider the use of the methods as a spectral methods in time. In
Figure 1 we plot the values of $|\gamma_{j}|$, $j=0,\dots,29$, for the CCM(30)
method used with timestep $h=2\pi/n$, $n=5,10,15,20$, for the first
integration step. As one may see, the Fourier coefficients decrease
exponentially, with $j$, and the basis of the exponential decreases with $h$
(i.e., they decrease as $\rho^{j}$, with $\rho\sim h^{-1}$). Moreover, when
using $h=\pi/10$ the last Fourier coefficients reach the round-off error
level, thus stagnating.

Figure 1: Plot of $|\gamma_{j}|$, $j=0,\dots,29$, for the CCM$(30)$ method,
for the first integration step of problem (58)-(59) by using the given
timestep $h$.

Now, let us solve problem (58)-(59) on ten periods, by using the CCM(50)
method with timestep $h=2\pi/n$, $n=3,6,9,12,15$. The errors (err), at the end
of each period, are listed in Table 2. As one my see, the errors are uniformly
small, as is expected from a spectral method.

At last, since problem (58)-(59) is Hamiltonian with Hamiltonian

$H(q,p)=\frac{1}{2}\|p\|_{2}^{2}-\frac{1}{\|q\|_{2}},$

let us consider the solution of problem (58)-(59) on the interval $[0,10^{3}]$
with timestep $h=0.1$, by using the CCM(3) and CCM(30) methods. The former
method is only symmetric, whereas the second one is spectral, for the
considered timestep. On the left of Figure 2 is the plot of the Hamiltonian
error for the CCM(3) method: though not energy conserving, the Hamiltonian
error is bounded, since the method is symmetric and the problem reversible. On
the right of the same figure, is the plot of the Hamiltonian error for the
CCM(30) method: in this case, the spectral accuracy reflects on the practical
conservation of the Hamiltonian (as well as of the additional constant of
motions, such as the angular momentum and the Lenz vector, not displayed
here). It must be emphasized that the execution times for running the two
methods (implemented in the same code) is almost the same: 2.9 sec vs. 3.7
sec. Considering that a similar behavior is observed when solving other
classes of problems (see also [5, 17] for a related analysis on HBVMs), it is
then clear that the use of CCMs as spectral methods can be very effective.

Figure 2: Plot of the Hamiltonian error when solving the problem (58)-(59)
with timestep $h=0.1$ one the interval $[0,10^{3}]$ by using the CCM(3) method
(left plot) and the CCM(30) method (right plot).

##  6 Conclusions

In this paper we have re-derived, in a novel framework, the Chebyshev
collocation methods of Costabile and Napoli [33], thus providing a more
comprehensive analysis of the methods, w.r.t. [33, 34], in particular when
used as spectral methods in time. The methods, derived by considering an
expansion of the vector field along the Chebyshev orthonormal polynomial
basis, turn out to be symmetric, perfectly $A$-stable, can have arbitrarily
high order, and the entries of their Butcher tableau can be written in closed
form. Their use as spectral methods in time has been also studied. A few
numerical tests confirm the theoretical achievements. As a further direction
of investigation, the efficient implementation of the methods will be
considered, as well as their application to different kinds of differential
problems.

### Declarations

The authors declare no conflict of interests.

### Acknoledgements

The authors wish to thank the mrSIR project [47] for the financial support.

## References

  * [1] P. Amodio, L. Brugnano, F. Iavernaro. Energy-conserving methods for Hamiltonian Boundary Value Problems and applications in astrodynamics. Adv. Comput. Math. 41 (2015) 881–905. https://doi.org/10.1007/s10444-014-9390-z
  * [2] P. Amodio, L. Brugnano, 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
  * [3] P. Amodio, L. Brugnano, F. Iavernaro. Continuous-Stage Runge-Kutta approximation to Differential Problems. Axioms 11 (2022) 192. https://doi.org/10.3390/axioms11050192
  * [4] P. Amodio, L. Brugnano, F. Iavernaro. Spectrally accurate solutions of nonlinear fractional initial value problems. AIP Conf. Proc. 2116 (2019) 140005. https://doi.org/10.1063/1.5114132
  * [5] P. Amodio, L. Brugnano, 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
  * [6] P. Amodio, L. Brugnano, F. Iavernaro. Arbitrarily high-order energy-conserving methods for Poisson problems. Numer. Algorithms (2022). https://doi.org/10.1007/s11075-022-01285-z
  * [7] P. Amodio, L. Brugnano, F. Iavernaro, C. Magherini. Spectral solution of ODE-IVPs using SHBVMs. AIP Conf. Proc. 2293 (2020) 100002. https://doi.org/10.1063/5.0026404
  * [8] L. Barletti, L. Brugnano, Y. Tang, B. Zhu. Spectrally accurate space-time solution of Manakov systems. J. Comput. Appl. Math. 377 (2020) 112918. https://doi.org/10.1016/j.cam.2020.112918
  * [9] S. Blanes, F. Casas. A Concise Introduction to Geometric Numerical Integration. CRC Press, Boca Raton, FL, USA, 2016.
  * [10] L. Brugnano, M. Calvo, I.J. Montijano, L. Ràndez. Energy preserving methods for Poisson systems. J. Comput. Appl. Math. 236 (2012) 3890–3904. http://doi.org/10.1016/j.cam.2012.02.033
  * [11] L. Brugnano, G. Frasca-Caccia, F. Iavernaro. Efficient implementation of Gauss collocation and Hamiltonian Boundary Value Methods. Numer. Algorithms 65 (2014) 633–650. https://doi.org/10.1007/s11075-014-9825-0
  * [12] L. Brugnano, G. Gurioli, F. Iavernaro, E. Weinmüller. Line integral solution of Hamiltonian systems with holonomic constraints. Appl. Numer. Math. 127 (2018) 56–77. https://doi.org/10.1016/j.apnum.2017.12.014
  * [13] L. Brugnano, G. Gurioli, C. Zhang. Spectrally accurate energy-preserving methods for the numerical solution of the “Good” Boussinesq equation. Numer. Methods Partial Differential Equations 35 (2019) 1343–1362. https://doi.org/10.1002/num.22353
  * [14] L. Brugnano, F. Iavernaro. Line integral methods which preserve all invariants of conservative problems. J. Comput. Appl. Math. 236 (2012) 3905–3919. http://doi.org/10.1016/j.cam.2012.03.026
  * [15] L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems. Chapman et Hall/CRC, Boca Raton, FL, USA, 2016.
  * [16] L. Brugnano, F. Iavernaro. Line Integral Solution of Differential Problems. Axioms 7(2) (2018) 36. https://doi.org/10.3390/axioms7020036
  * [17] L. Brugnano, F. Iavernaro, J.I. Montijano, 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
  * [18] L. Brugnano, F. Iavernaro, T. Susca. Numerical comparisons between Gauss-Legendre methods and Hamiltonian BVMs defined over Gauss points. Monogr. Real Acad. Cienc. Zaragoza 33 (2010) 95–112.
  * [19] L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian BVMs (HBVMs): a family of “drift-free” methods for integrating polynomial Hamiltonian systems. AIP Conf. Proc. 1168 (2009) 715–718. http://doi.org/10.1063/1.3241566
  * [20] L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian Boundary Value Methods (Energy Preserving Discrete Line Integral Methods). JNAIAM J. Numer. Anal. Ind. Appl. Math. 5 (2010) 17–37.
  * [21] L. Brugnano, F. Iavernaro, 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
  * [22] L. Brugnano, F. Iavernaro, D. Trigiante. The lack of continuity and the role of infinite and infinitesimal in numerical methods for ODEs: The case of symplecticity. Appl. Math. Comput. 218 (2012) 8053–8063. https://doi.org/10.1016/j.amc.2011.03.022
  * [23] L. Brugnano, F. Iavernaro, 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
  * [24] L. Brugnano, F. Iavernaro, D. Trigiante. A two-step, fourth-order method with energy preserving properties. Comput. Phys. Commun. 183 (2012) 1860–1868. https://doi.org/10.1016/j.cpc.2012.04.002
  * [25] L. Brugnano, F. Iavernaro, 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
  * [26] L. Brugnano, F. Iavernaro, R. Zhang. Arbitrarily high-order energy-preserving methods for simulating the gyrocenter dynamics of charged particles. J. Comput. Appl. Math. 380 (2020) 112994. https://doi.org/10.1016/j.cam.2020.112994
  * [27] L. Brugnano, J.I. Montijano, L. Rández. High-order energy-conserving Line Integral Methods for charged particle dynamics. J. Comput. Phys. 396 (2019) 209–227. https://doi.org/10.1016/j.jcp.2019.06.068
  * [28] L. Brugnano, J.I. Montijano, 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://doi.org/10.1007/s11075-018-0552-9
  * [29] L. Brugnano, Y. Sun. Multiple invariants conserving Runge-Kutta type methods for Hamiltonian problems. Numer. Algorithms 65 (2014) 611–632. https://doi.org/10.1007/s11075-013-9769-9
  * [30] J.C. Butcher. B-series. Algebraic analysis of numerical methods. Springer Nature, Switzerland, 2021.
  * [31] E. Celledoni, R.I. McLachlan, D.I. McLaren, B. Owren, G.R.W. Quispel, W.M. Wright. Energy-preserving Runge-Kutta methods. M2AN Math. Model. Numer. Anal. 43 (2009) 645–649. https://doi.org/10.1051/m2an/2009020
  * [32] P. Chartier, E. Faou, A. Murua, An algebraic approach to invariant preserving integrators: the case of quadratic and Hamiltonian invariants. Numer. Math. 103 (2006) 575–590. https://doi.org/10.1007/s00211-006-0003-8
  * [33] F. Costabile, A. Napoli. A method for global approximation of the initial value problem. Numer. Algorithms 27  (2001) 119–130.
  * [34] F. Costabile, A. Napoli. Stability of Chebyshev collocation methods. Comput. Math. Appl. 47 (2004) 659–666.
  * [35] F. Costabile, A. Napoli. A class of collocation methods for numerical integration of initial value problems. Comput. Math. Appl. 62 (2011) 3221–3235.
  * [36] E. Hairer. Energy preserving variant of collocation methods. JNAIAM J. Numer. Anal. Ind. Appl. Math. 5 (2010) 73–84.
  * [37] E. Hairer, Ch. Lubich, G. Wanner. Geometric Numerical Integration, 2nd ed. Springer, Berlin, Germany, 2006.
  * [38] E. Hairer, G. Wanner. Solving ordiary differential equations II: Stiff and differential-algebraic problems, second ed. Springer-Verlag, Berlin, 1996.
  * [39] F. Iavernaro, 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
  * [40] F. Iavernaro, 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
  * [41] F. Iavernaro, D. Trigiante. High-order Symmetric Schemes for the Energy Conservation of Polynomial Hamiltonian Problems. JNAIAM J. Numer. Anal. Ind. Appl. Math. 4 (2009) 87–101.
  * [42] B. Leimkuhler, S. Reich. Simulating Hamiltonian Dynamics. Cambridge University Press, Cambridge, UK, 2004.
  * [43] Y.-W. Li, X. Wu. Functionally fitted energy-preserving methods for solving oscillatory nonlinear Hamiltonian systems. SIAM J. Numer. Anal. 54 (2016) 2036–2059. https://doi.org/10.1137/15M1032752
  * [44] R.I. McLachlan, G.R.W. Quispel, N. Robidoux. Geometric integration using discrete gradients. Phil. Trans. R. Soc. Lond. A 357 (1999) 1021–1045.
  * [45] J.M. Sanz-Serna. Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more. SIAM Rev. 58 (2016) 3–33. https://doi.org/10.1137/151002769
  * [46] J.M. Sanz-Serna, M.P. Calvo. Numerical Hamiltonian Problems. Chapman & Hall, London, UK, 1994.
  * [47] https://www.mrsir.it/en/about-us/