††thanks: These authors contribute equally to this work.††thanks: These
authors contribute equally to this work.

# Topological Origin of Floquet Thermalization
in Periodically Driven Many-body Systems

Hao-Yue Qi  Hefei National Research Center for Physical Sciences at the
Microscale and School of Physical Sciences, University of Science and
Technology of China, Hefei 230026, China  CAS Center for Excellence in Quantum
Information and Quantum Physics, University of Science and Technology of
China, Hefei 230026, China  Yue Wu  Institute for Advanced Study, Tsinghua
University, Beijing 100084, China  Wei Zheng<EMAIL_ADDRESS>Hefei
National Research Center for Physical Sciences at the Microscale and School of
Physical Sciences, University of Science and Technology of China, Hefei
230026, China  CAS Center for Excellence in Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei 230026, China
Hefei National Laboratory, University of Science and Technology of China,
Hefei 230088, China

###### Abstract

Floquet engineering is a powerful manipulation method in modern quantum
technology. However, unwanted heating is the main challenge of Floquet
engineering, therefore the Floquet thermalization has attracting considerable
attentions recently. In this work, we investigate thermalization of
periodically driven many-body systems through the lens of Krylov complexity,
and find a topological origin of different thermalization behaviors. We
demonstrate that If the topology of the Krylov chain is nontrivial, a
periodically driven system will reach a state with finite temperature. When
the Krylov chain is topologically trivial, the system will be heated to
infinite temperature. We further show that the prethermalization can be
understood as the tunnelling process of a quasi-edge mode through the local
gap on Krylov chain. This picture provides a systematically method to obtain
the effective prethermal Hamiltonian.

Quantum systems with periodic driving have attracted much attention in the
decades. First, periodic driving can be used to engineer nontrivial
Hamiltonians. For example, periodic driving has been used to generate
artificial magnetic fields in optical lattice, simulate synthetic lattice
gauge theories with ultracold atomic gases<EMAIL_ADDRESS>, and
engineer topological bands<EMAIL_ADDRESS><EMAIL_ADDRESS>Haldane . Second, periodic driving can induce intrinsically out-of-equilibrium
phenomena without counterpart in equilibrium, such as discrete time-crystals
<EMAIL_ADDRESS>and anomalous Floquet topological states
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Floquet<EMAIL_ADDRESS>.However, since periodic driving breaks
energy conservation, it was expected that under driving, generic many-body
systems would be heated up to featureless infinite temperature ensembles
eventually<EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS>. Fortunately, under high-frequency
driving, there is a class of systems that exhibit the prethermalization
phenomena<EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>REV . Such systems have several approximate conservation quantities. Thus, at
intermediate time scales, systems will first equilibrate to a finite
temperature state governed by a so-called prethermal Hamiltonian. Later,
systems will thermalize to infinite temperature ensembles at much longer time
scales. Prethermalization has been experimentally observed in one-dimensional
Bose gas and spins with dipole interaction<EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS>. However, in what
circumstances, a quantum system will be directly heated to infinite
temperature or exhibit prethermalization behavior is not clear. There is no
simple and clear physical picture of the Floquet thermalization.

Recently, the Krylov complexity of operators has been developed to study the
chaos and thermalization in many-body systems<EMAIL_ADDRESS>. With the
help of Krylov basis, the evolution of an operator under Heisenberg equation
can be mapped to the spreading of a wavepacket in a one-dimensional tight-
binding chain. Krylov complexity has motivated a lot of interesting studies in
closed<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS>and open<EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS>quantum many-body
systems.

In this work, we found topological origin of different thermalization behavior
of periodically driven many-body systems based on Krylov complexity. We first
build time-independent Krylov basis for periodically driven systems, and map
the Floquet dynamics of an operator to the wavepacket evolution on the Krylov
chain. The main results of this work are summarized here:

(1) We construct the unfold Floquet effective Hamiltonian by seeking the zero
mode on Krylov chain. Unlike the effective Hamiltonian obtained by taking the
logarithm of evolution operator in previous works, which has folded spectrum
and does not satisfy the assumption of eigenstate thermalization hypothesis
(ETH), our unfold effective Hamiltonian satisfies ETH and can be used to
estimate local observables at long time scale.

(2) We further find that it is the topology of Krylov chain that determined
the Floquet thermalization behavior. In a finite system, when the Krylov chain
is topologically nontrivial, it sustains a localized zero edge mode. Then the
system will heated to a state with finite temperature. If topology of the
Krylov chain is trivial, the zero mode is extended to the bulk. As a result,
the system will be heated to infinite temperature.

(3) In the thermodynamic limit, the Krylov chain is gapless. Therefore there
is no eigen zero edge mode, and the system will finally be heated to infinite
temperature. However, we find that there could exist a quasi edge mode in high
driving frequency regime, which is protected by the local topology and gap of
the Krylov chain. We further conjecture that this quasi edge mode represents
the Floquet prethermal Hamiltonian. Its life time determines the time scale of
intermediate prethermal state. With this picture in mind, one can
systematically estimate the prethermal Hamiltonian via Krylov basis.

Figure 1: Hopping amplitude $b_{n}$ and zero mode
$|H_{\mathrm{eff}}^{\mathrm{zm}})$ on the Krylov chain. (a)(c) quantum spin
chain with system size $L=14$ and $h_{z}/J=0.809,h_{x}/J=0.9045$. (b)(d) hard-
core Bose-Hubbard model with system size $L=16$ and $J^{\prime}/J=0.2,U/J=1$.
For both two models, the hopping amplitude fluctuates around a non-zero value
site-by-site in high-frequency regime, and the fluctuation becomes smaller
deep in the bulk. The zero mode is an edge mode. In contrast, in the low-
frequency regime, the hopping amplitude is nearly constant. The zero mode
extends over the bulk.

Krylov Chain. In a periodically driven system, $\hat{H}(t+T)=\hat{H}(t)$, the
evolution operator during one driving period is given by
$\hat{U}_{\mathrm{F}}=\mathcal{\hat{T}}\exp\left[-i\int_{0}^{T}dt\hat{H}(t)\right]$,
here we set $\hbar=1$. One can define a time-independent effective Hamiltonian
as

$\hat{U}_{\mathrm{F}}=e^{-i\hat{H}_{\mathrm{eff}}T}.$ (1)

This effective Hamiltonian reproduces the evolution of the driven system in a
stroboscopic fashion. In other words, the state evolution under time-
independent $\hat{H}_{\mathrm{eff}}$ is identical to the one under driven
Hamiltonian $\hat{H}(t)$, only at steps of the driving period $t=nT$. Thus the
effective Hamiltonian well describes the dynamics over time scale much longer
than a single driving period. Numerically, we calculate the effective
Hamiltonian by taking the logarithm of $\hat{U}_{\mathrm{F}}$, as
$\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}=\frac{i}{T}\ln\hat{U}_{\mathrm{F}}$.
Therefore, the long-time evolution of operators can be captured by the
following Heisenberg equation in a stroboscopic way,

$i\partial_{t}\hat{O}=\left[\hat{O},\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}\right].$
(2)

However, the effective Hamiltonian obtained by this method is not unique. It
has a folded energy spectrum. Specifically speaking,
$\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}$ can be diagonalized into
$\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}=\sum\nolimits_{k}E_{k}\left|\psi_{k}\right\rangle\left\langle\psi_{k}\right|$,
and its eigen energies are folded in the first Floquet-Brillouin zone,
$-\pi/T<E_{k}<\pi/T$. As a result, it does not satisfy the assumption of ETH,
and can not predict the expectation values of local observables after a long
time evolution. Certainly, we can construct other equivalent effective
Hamiltonian by shifting the eigenenergies without altering the eigenstates as
$\hat{H}_{\mathrm{eff}}=\sum\nolimits_{k}\left(E_{k}+m_{k}\Omega\right)\left|\psi_{k}\right\rangle\left\langle\psi_{k}\right|$,
where $m_{k}$ is an integer and $\Omega=2\pi/T$ is the driving frequency. Such
Hamiltonians satisfy the definition (1), and thus commute with
$\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}$.

The Krylov basis of operators is based on the Lanczos algorithm. An operator
is denoted by
$\hat{O}=\sum\nolimits_{ij}O_{ij}\left|i\right\rangle\left\langle j\right|$,
where $\\{\left|i\right\rangle\\}$ is a set of orthogonal states and $O_{ij}$
are corresponding matrix elements. We can map the operator $\hat{O}$ to a
vector in the operator space as
$\hat{O}\longrightarrow\left|O\right)=\sum\nolimits_{ij}O_{ij}\left|i\right\rangle\otimes\left|j\right\rangle$.
The corresponding inner product of operators is defined as
$\left(O_{1}|O_{2}\right)=\mathrm{Tr}(\hat{O}_{1}^{{\dagger}}\hat{O}_{2})$. We
further introduce a Liouville super-operator $\mathcal{L}$, acting on an
operator as
$\mathcal{L}\hat{O}\equiv[\hat{H}_{\mathrm{eff}}^{\mathrm{fold}},\hat{O}]$. In
the operator space, it can be mapped into a matrix,
$\mathcal{L}\hat{O}\rightarrow\mathcal{L}\left|O\right)$. Thus the Heisenberg
equation (2) can be expressed into
$i\partial_{t}\left|O\right)=-\mathcal{L}\left|O\right)$. Applying the the
Baker-Campbell-Hausdorff formula, we expand the evolution of operator into
$\left|O(t)\right)=\sum\nolimits_{n}\frac{\left(-it\right)^{n}}{n!}\mathcal{L}^{n}\left|O\right)$.
Here $\left\\{\mathcal{L}^{n}\left|O\right)\right\\}$ forms a complete set of
operator bases, but it is neither normalized nor orthogonal. We then apply the
Gram-Schmidt procedure to obtain a set of orthogonalized and normalized basis,
called Krylov basis as,

$\displaystyle\left|W_{0}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{b_{0}}\left|O\right),$ (3)
$\displaystyle\left|W_{1}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{b_{1}}\mathcal{L}\left|W_{0}\right),$ (4)
$\displaystyle\left|W_{n}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{b_{n}}\left[\mathcal{L}\left|W_{n-1}\right)-b_{n-1}\left|W_{n-2}\right)\right],n\geq
2.$ (5)

Here $b_{n}$ are the so-called Lanczos coefficients. They are matrix elements
of Liouville super-operator in the new basis,
$b_{n}=\left(W_{n}\right|\mathcal{L}\left|W_{n+1}\right)$. We expand the
operator in the Krylov basis as
$\left|O(t)\right)=\sum\nolimits_{n}\varphi_{n}(t)\left|W_{n}\right)$. The
coefficients $\varphi_{n}(t)$ satisfy the following linear equation,

$i\partial_{t}\varphi_{n}=-\left(b_{n-1}\varphi_{n-1}+b_{n}\varphi_{n+1}\right).$
(6)

Therefore, the operator evolution is mapped to a quantum mechanics problem of
one particle hopping on a tight-binding half-infinite chain. Note that the
$b_{n}$ play the role of the hopping amplitude on this half-infinite Krylov
chain. The particle is initially localized at original site, $n=0$. During the
evolution, it hops into the bulk. The Krylov complexity is defined as
$\left(O(t)\right|\mathcal{K}\left|O(t)\right)=\sum\nolimits_{n}n|\varphi_{n}(t)|^{2}$,
where $\mathcal{K=}\sum\nolimits_{n}n\left|W_{n}\right)\left(W_{n}\right|$ is
the Krylov super-operator. It is nothing but the center-of-mass of the wave
function on the Krylov chain.

Models. To illustrate the Floquet thermalization, we study two models with
periodical driving. One is a quantum spin chain, the other is a one-
dimensional hard-core Bose-Hubbard Model. The evolution operators of both
models in one period are

$\hat{U}_{\mathrm{F}}^{\alpha}=e^{-i\hat{H}_{\alpha}T/2}e^{-i\hat{V}_{\alpha}T/2},$
(7)

where $\alpha=\mathrm{spin,boson}$, for spin and Bose-Hubbard models
respectively. The corresponding $\hat{H}_{\alpha}$ are

$\displaystyle\hat{H}_{\mathrm{spin}}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{L}(-J\sigma_{i}^{z}\sigma_{i+1}^{z}+h_{z}\sigma_{i}^{z}),$
(8) $\displaystyle\hat{H}_{\mathrm{boson}}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{L}(J\hat{a}_{i}^{{\dagger}}\hat{a}_{i+1}+J^{\prime}\hat{a}_{i}^{{\dagger}}\hat{a}_{i+2}+\mathrm{H.c.}).$
(9)

And $\hat{V}_{\alpha}$ are given by
$\hat{V}_{\mathrm{spin}}=\sum_{i=1}^{L}h_{x}\sigma_{i}^{x}$,
$\hat{V}_{\mathrm{boson}}=\sum_{i=1}^{L}U\hat{n}_{i}\hat{n}_{i+1}$, where
$\hat{n}_{i}=\hat{a}_{i}^{{\dagger}}\hat{a}_{i}$. Both models are chaotic.

Figure 2: Comparison of zero mode and effective Hamiltonian obtained via
Floquet-Magnus expansion in high-frequency regime on Krylov chain. (a) quantum
spin chain at $TJ=1$, (b) hard-core Bose-Hubbard model at $TJ=1.3$. For both
models, the zero mode is consistent with the Floquet-Magnus expansion. The
simulation parameters are the same as in Fig.1.

Unlike previous studies<EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>, we
choose the seed operator as the average Hamiltonian in one period,

$\hat{O}=\hat{H}^{(0)}=\frac{1}{T}\int\nolimits_{0}^{T}dt\hat{H}(t).$ (10)

For the two models investigated here, we have
$\hat{H}_{\alpha}^{(0)}=(\hat{H}_{\alpha}+\hat{V}_{\alpha})/2$. We numerically
calculate the evolution operator $\hat{U}_{\mathrm{F}}$ and its logarithm,
i.e. $\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}$ with finite system size. Then we
obtain the Krylov basis by the method described above. After that, we expand
super-operator $\mathcal{L}$ on this basis, and obtain the hopping amplitude
on the Krylov chain, $b_{n}$, with varying driving frequencies. Typical
distributions of $b_{n}$ are plotted in Fig.1(a)(b). We find that, for both
two models, the hopping amplitude fluctuate around a non-zero value site-by-
site in high-frequency driving regime. It reminds us the celebrated Su-
Schrieffer-Heeger (SSH) model, which possesses nontrivial topological band
structure and can sustain ingap edge states. The fluctuation is large near the
origin site, but becomes smaller deep in the bulk. In contrast, in the low-
frequency regime, the fluctuation of hopping amplitude is highly suppressed,
and $b_{n}$ is nearly constant.

We then numerically seek the zero mode operator of Liouvillian,
$\mathcal{L}\left|O\right)=0$. That is equivalent to solving the single-
particle eigen wave function with zero energy on the Krylov chain, see
Fig.1(c)(d). We further conclude that the zero mode obtained by this method is
one particular choice of effective Hamiltonian
$\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}$ with unfolded energy spectrum. The
reasons are as follows. First, this effective Hamiltonian commutes with
$\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}$, thus is a zero mode of Liouvillian,
$\mathcal{L}\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}=[\hat{H}_{\mathrm{eff}}^{\mathrm{fold}},\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}]=0$.
Second, we analytically calculate the effective Hamiltonian via Floquet-Magnus
expansion up to the 5th order in the high-frequency regime<EMAIL_ADDRESS>,
and compare it with the zero mode on the Krylov basis. The results are plotted
in Fig.2. Note that for both models, the zero mode is consistent with the
Floquet-Magnus expansion.

We note that the distribution of $\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}$ on the
Krylov chain varies as driving period changes, see Fig.1(c)(d). In short
driving period regime, the effective Hamiltonian is a edge mode, localizing
near the origin of Krylov chain. Thus it possesses small Krylov complexity. In
long driving period regime, $\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}$ extends
over the bulk, and has larger Krylov complexity. The behaviors of the zero
mode can be understood by the band topology of the Krylov chain. Since
distribution of $b_{n}$ does not possess the translation symmetry, it is hard
to investigate the topology in momentum space. We turn to calculate the
topological invariance, the winding number, of Krylov chain via a real-space
approach as<EMAIL_ADDRESS><EMAIL_ADDRESS>;
<EMAIL_ADDRESS>,

$\nu=\frac{1}{2N}\mathrm{Tr}\left(\mathcal{SQ}\left[\mathcal{Q},\mathcal{K}\right]\right).$
(11)

Here
$\mathcal{S=}\sum\nolimits_{n}\left(\left|W_{2n}\right)\left(W_{2n}\right|-\left|W_{2n+1}\right)\left(W_{2n+1}\right|\right)$
is the chiral-symmetry super-operator on the Krylov chain, and $N$ is the
lattice length. One note that $\mathcal{SLS}^{{\dagger}}\mathcal{=-L}$, thus
the eigen values of $\mathcal{L}$ come in pairs. The super-operator
$\mathcal{Q}$ is given by
$\mathcal{Q=}\sum\nolimits_{i}\left(\left|A_{i}\right)\left(A_{i}\right|-\left|A_{i}^{\prime}\right)\left(A_{i}^{\prime}\right|\right)$,
where $\left|A_{i}\right)$ are the eigen modes of $\mathcal{L}$ with positive
eigen values $\lambda_{i}$,
$\mathcal{L}\left|A_{i}\right)=\lambda_{i}\left|A_{i}\right)$, and
$\left|A_{i}^{\prime}\right)$ are chiral-symmetry partners of
$\left|A_{i}\right)$,
$\left|A_{i}^{\prime}\right)=\mathcal{S}\left|A_{i}\right)$. Therefore they
have negative eigen values,
$\mathcal{L}\left|A_{i}^{\prime}\right)=-\lambda_{i}\left|A_{i}^{\prime}\right)$.
We numerically calculate the winding number $\nu$, and compare it with the
overlap $(H^{(0)}|H_{\mathrm{eff}}^{\mathrm{zm}})$ as a function of driving
frequency. The latter represents the weight of zero mode on the origin site,
and can be used to characterize the edge state. If
$(H^{(0)}|H_{\mathrm{eff}}^{\mathrm{zm}})\sim 1$, that indicates the zero mode
is well localized on the edge. The results are plotted in Fig.3(d) and
Fig.4(d). Note that the band topology of the Krylov chain determines the zero
edge mode. In the short driving period regime, winding number $\nu\approx 1$,
that protects a zero edge mode on the boundary. In long driving period regime,
since the gap is closed, the winding number is not well defined, and gives a
non-integer. Thus zero mode is extended to bulk.

Figure 3: Verifying ETH and winding number $\nu$ of Krylov chain in the
quantum spin model. (a-c) Long-time equilibrium values of
$\langle\hat{\sigma}^{x}_{1}(mT)\rangle$,
$\langle\hat{\sigma}^{z}_{1}\hat{\sigma}^{z}_{2}(mT)\rangle$ and
$\langle\hat{H}^{(0)}(mT)\rangle$, as $m\rightarrow\infty$, are compared to
ETH predictions at different driving periods. The dashed line marks the
infinite-temperature value. (d) Local winding number $\nu$ of the Krylov chain
is compared with the overlap $(H^{(0)}|H_{\mathrm{eff}}^{\mathrm{zm}})$, which
exhibits the localization of zero mode on the first Krylov site. We work in
zero momentum sector of positive parity. The initial state is chosen as
$Z_{2}$ state $\mathcal{P}|\uparrow\downarrow\cdots\uparrow\downarrow\rangle$,
where $\mathcal{P}$ is the projector onto this sector. The simulation
parameters are the same as in Fig.1.

Unlike the folded effective Hamiltonian
$\hat{H}_{\mathrm{eff}}^{\mathrm{fold}}$, we found that this zero mode
$\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}$ satisfies the assumption of ETH SM. We
can use it to predict the long-time equilibrium values of local observables as

$\left\langle\hat{O}_{\mathrm{local}}(mT\rightarrow\infty)\right\rangle\approx\mathrm{Tr}\left(\hat{\rho}_{\beta}\hat{O}_{\mathrm{local}}\right),$
(12)

where the canonical ensemble density matrix,
$\hat{\rho}_{\beta}=e^{-\beta\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}}/{\mathrm{Tr}}(e^{-\beta\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}})$.
The inverse temperature $\beta$ is obtained by solving a single equation,
$\left\langle\psi(0)\right|\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}\left|\psi(0)\right\rangle=\mathrm{Tr}(\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}\hat{\rho}_{\beta})$.
Here $\left|\psi(0)\right\rangle$ is a generic initial pure state. In Fig.3
and Fig.4, we plot the comparison of local observables after long time
evolution with the ETH prediction. All these calculations start from the same
initial pure state. Note that, for both models, the long time expectation
values are consistent with the ETH prediction not even in high frequency
limit, but in the whole frequency regime. By comparing the long time
thermalization behavior with $(H^{(0)}|H_{\mathrm{eff}}^{\mathrm{zm}})$ (Fig.3
and Fig.4), we note that when the zero mode is localized on the edge, system
will be driven into an ensemble with finite temperature. If the zero mode is
extended to the bulk, the system will be heated to infinite temperature. Thus
we conclude that it is the band topology of Krylov chain that determines the
long time thermalization of periodically driven many-body systems. When the
Krylov chain is topologically non-trivial, the Floquet heating will be bounded
to finite energy density. When the Krylov chain is topological trivial or
gapless, the heating will be unbounded.

Figure 4: Verifying ETH and winding number $\nu$ of Krylov chain in the hard-
core Bose-Hubbard model. (a-c) Long-time equilibrium values of
$\langle\hat{b}_{1}^{\dagger}\hat{b}_{2}(mT)+h.c.\rangle$,
$\langle\hat{n}_{1}\hat{n}_{2}(mT)\rangle$ and
$\langle\hat{H}^{(0)}(mT)\rangle$, as $m\rightarrow\infty$, are compared to
ETH predictions at different driving periods. The dashed line marks the
infinite-temperature value. (d) Local winding number $\nu$ of the Krylov chain
is compared with the overlap $(H^{(0)}|H_{\mathrm{eff}}^{\mathrm{zm}})$, which
exhibits the localization of zero mode on the first Krylov site. We work in
zero momentum sector of positive parity. The initial state is chosen as ground
state of $\hat{H}_{\mathrm{boson}}$. The simulation parameters are the same as
in Fig.1.

Prethermalization. In the thermodynamic limit, the hopping amplitude $b_{n}$
at large $n$ approaches a constant. Therefore, the Krylov chain is gapless,
and there is no eigen edge mode. As a result, the system will be finally
heated to infinite temperature despite the frequency. However, in the high
frequency limit, the configuration of hopping amplitude near the origin site
is still similar to the SSH model. That indicates, there exists local gap near
the origin. To see this, we plot the local density-of-state
$\rho_{n}\left(\omega\right)=-\frac{1}{\pi}\mathrm{Im}\left\langle
W_{n}\right|\left(\omega-\mathcal{L}+i0^{+}\right)^{-1}\left|W_{n}\right\rangle$
on the Krylov chain in Fig.5. Such local gap will sustain a quasi eigen edge
mode, denoted by $\hat{H}_{\mathrm{eff}}^{\mathrm{quasi}}$. This quasi edge
mode has a finite lifetime, $\tau_{\ast}$. Before this time scale, the initial
wavepacket stays localized near the edge. But after $\tau_{\ast}$, it will
tunnel across the local gap, and spread to the bulk.

The damping of this quasi edge mode in the Krylov chain describes the Floquet
prethermalization. The lifetime of the quasi edge mode $\tau_{\ast}$ is just
the time scale of prethermalization. During this time scale, the system will
thermalize to a quasi steady state with finite temperature given by
$\hat{H}_{\mathrm{eff}}^{\mathrm{quasi}}$. After that, the system will be
continuously heated to infinite temperature. To see that, we plot the
thermalization process of local observables in a long time scale. We find that
their prethermal values are consistent with the prediction of ETH given by
$\hat{H}_{\mathrm{eff}}^{\mathrm{quasi}}$, but away from the prediction of
$\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}$, and $\hat{H}^{(0)}$ see Fig.5(a)(c).
Numerically, we approximate the $\hat{H}_{\mathrm{eff}}^{\mathrm{quasi}}$ by
truncating the zero mode to the length of local gap regime. Therefore, we
conclude that the Floquet prethermalization process is dominated by the local
topological structure on Krylov chain. With this picture in mind, one can
optimize the prethermal Hamiltonian as the proper quasi edge mode
$\hat{H}_{\mathrm{eff}}^{\mathrm{quasi}}$.

Figure 5: Real time thermalization dynamics and local density-of-state on
Krylov chains. (a-b) is the quatnum spin model, with $L=14$ and $TJ=1.2$.
(c-d) is the hard-core Bose-Hubbard model with $L=18$ and $TJ=3.9$. (a)(c) The
real time thermalization dynamics of local observable
$\langle\hat{H}^{(0)}_{\alpha}(mT)\rangle$ in a long time scale. The predicted
values of ETH given by $\hat{H}_{\mathrm{eff}}^{(0)}$,
$\hat{H}_{\mathrm{eff}}^{\mathrm{quasi}}$ and
$\hat{H}_{\mathrm{eff}}^{\mathrm{zm}}$ are shown by green, yellow and red
lines. The dotted gray vertical line marks the value of prethemalization time
scale $\tau_{*}$. (b)(d) Local density-of-state on Krylov chains. Note that
there exists local gap near the origin for both models.

Summary. In summary, we study thermalization of periodically driven many-body
quantum systems via Krylov basis. We found that the zero mode on the Krylov
chain represents the unfolded Floquet effective Hamiltonian. We demonstrate
that the band topology of the Krylov chain dominates the thermalization
behavior. We show that the prethermalization phenomena can be understood as
the damping process of a quasi edge mode protected by the local gap. Our
results bring out the generic connection between Krylov complexity and
thermalization of periodically driven systems that deserves further studies.

## Acknowledgments

We thank H. Zhai and Z. Wang for helpful discussion. This work is supported by
NSFC (Grant No. GG2030007011 and No. GG2030040453) and Innovation Program for
Quantum Science and Technology (No.2021ZD0302004).

## References

  * (1) C. Schweizer, F. Grusdt, M. Berngruber, L. Barbiero, E. Demler, N. Goldman, I. Bloch, and M. Aidelsburger, Floquet approach to Z2 lattice gauge theories with ultracold atoms in optical lattices, Nat. Phys. 15, 1168 (2019)
  * (2) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, T. Esslinger, Experimental realization of the topological Haldane model with ultracold fermions, Nature 515, 237 (2014)
  * (3) W. Zheng and H. Zhai, Floquet topological states in shaking optical lattices, Phys. Rev. A 89, 061603(R) (2014)
  * (4) M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao, Quantum and classical discrete time crystals, Rev. Mod. Phys. 95, 031001 (2023)
  * (5) T. Kitagawa, E. Berg, M. Rudner, and E. Demler, Topological characterization of periodically driven quantum systems, Phys. Rev. B 82, 235114 (2010)
  * (6) M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin, Anomalous Edge States and the Bulk-Edge Correspondence for Periodically Driven Two-Dimensional Systems, Phys. Rev. X, 3, 031005 (2013)
  * (7) K. Wintersperger, C. Braun, F. Nur Ünal, A. Eckardt, M. Di Liberto, N. Goldman, I. Bloch, M. Aidelsburger, Realization of an anomalous Floquet topological system with ultracold atoms, Nat. Phys. 16, 1058 (2020)
  * (8) P. Xu, W. Zheng and H. Zhai, Topological micromotion of Floquet quantum systems, Phys. Rev. B 105, 045139 (2022)
  * (9) A. Lazarides, A. Das, and R. Moessner, Equilibrium states of generic quantum systems subject to periodic driving, Phys. Rev. E 90, 012110 (2014)
  * (10) L. D.lessio and M. Rigol, Long-time Behavior of Isolated Periodically Driven Interacting Lattice Systems, Phys. Rev. X 4, 041048 (2014)
  * (11) H. Kim, T. N. Ikeda, and D. A. Huse, Testing whether all eigenstates obey the eigenstate thermalization hypothesis, Phys. Rev. E 90, 052105 (2014)
  * (12) P. Ponte, A. Chandran, Z. Papić, and D. A. Abanin, Periodically driven ergodic and many-body localized quantum systems, Annals of Physics 353, 196 (2015)
  * (13) D. A. Abanin, W. De Roeck, and F. Huveneers, Exponentially Slow Heating in Periodically Driven Many-Body Systems, Phys. Rev. Lett. 115, 256803 (2015)
  * (14) T. Mori, T. Kuwahara, and K. Saito, Rigorous Bound on Energy Absorption and Generic Relaxation in Periodically Driven Quantum Systems, Phys. Rev. Lett. 116, 120401 (2016)
  * (15) T. Kuwahara, T. Mori, K. Saitoc, Floquet-Magnus theory and generic transient dynamics in periodically driven many-body quantum systems, Annals of Physics 367, 96 (2016)
  * (16) D. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, A Rigorous Theory of Many-Body Prethermalization for Periodically Driven and Closed Quantum Systems. Comm. Math. Phys. 354, 809 (2017).
  * (17) M. Bukov, M. Heyl, D. A. Huse, and A. Polkovnikov, Heating and many-body resonances in a periodically driven two-band system, Phys. Rev. B 93, 155132 (2016)
  * (18) S. A. Weidinger and M. Knap, Floquet prethermalization and regimes of heating in a periodically driven, interacting quantum system, Scientific Reports,
  * (19) E. Canovi, M. Kollar, and M. Eckstein, Stroboscopic prethermalization in weakly interacting periodically driven systems, Phys. Rev. E 93, 012130 (2016)
  * (20) D. A. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, Effective Hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems Phys. Rev. B 95, 014112 (2017)
  * (21) C. Fleckenstein and M. Bukov, Thermalization and prethermalization in periodically kicked quantum spin chains, Phys. Rev. B 103, 144307 (2021)
  * (22) W. W. Ho, T. Mori, D. A. Abanin, E. G. Dalla Torre, Quantum and classical Floquet prethermalization, Annals of Physics 454, 169297 (2023)
  * (23) A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Floquet Prethermalization in a Bose-Hubbard System, Phys. Rev. X 10, 021044 (2020)
  * (24) P. Peng, C. Yin, X. Huang, C. Ramanathan, and P. Cappellaro, Floquet prethermalization in dipolar spin chains, Nat. Phys. 17, 444 (2021)
  * (25) W. Beatrez, O. Janes, A. Akkiraju, A. Pillai, A. Oddo, P. Reshetikhin, E. Druga, M. McAllister, M. Elo, B. Gilbert, D. Suter, and A. Ajoy Floquet Prethermalization with Lifetime Exceeding 90 s in a Bulk Hyperpolarized Solid, Phys. Rev. Lett. 127, 170603 (2021)
  * (26) D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, A universal operator growth hypothesis, Phys. Rev. X 9, 041017 (2019).
  * (27) C. Lv, R. Zhang, and Q, Zhou, Building Krylov complexity from circuit complexity, arXiv:2303.07343
  * (28) R. Zhang and H. Zhai, Universal Hypothesis of Autocorrelation Function from Krylov Complexity, arXiv:2305.02356
  * (29) A. Dymarsky and M. Smolkin, Krylov complexity in conformal field theory, arXiv:2104.09514
  * (30) E. Rabinovici, A. Sánchez-Garrido, R. Shira and J. Sonnerb, Krylov localization and suppression of complexity, JHEP 03, 211 (2022)
  * (31) P. Caputa, J. M. Magan, and D. Patramanis, Geometry of Krylov complexity, Phys. Rev. Res. 4, 013041 (2022)
  * (32) Universal chaotic dynamics from Krylov space, J. Erdmenger, S.-K. Jian, and Z.-Y. Xian, arXiv:2303.12151
  * (33) A. Dymarsky and A. Gorsky, Quantum chaos as delocalization in Krylov space, Phys. Rev. B 102, 085137 (2020)
  * (34) D. J. Yates, A. G. Abanov, and A. Mitra, Lifetime of Almost Strong Edge-Mode Operators in One-Dimensional, Interacting, Symmetry Protected Topological Phases, Phys. Rev. Lett. 124, 206803 (2020)
  * (35) H.-C. Yeh, G. Cardoso, L. Korneev, D. Sels, A. G. Abanov, and A. Mitra, Slowly decaying zero mode in a weakly nonintegrable boundary impurity model, Phys. Rev. B 108, 165143 (2023)
  * (36) P. Zhang and Z. Yu, Dynamical transition of operator size growth in quantum systems embedded in an environment, Phys. Rev. Lett. 130, 250401 (2023)
  * (37) L. Chen, B. Mu, H. Wang, and P. Zhang, Dissecting Quantum Many-body Chaos in the Krylov Space, arXiv:2404.08207
  * (38) D. J. Yates and A. Mitra, Strong and almost strong modes of Floquet spin chains in Krylov subspaces Phys. Rev. B 104, 195121 (2021)
  * (39) D, J. Yates, A. G. Abanov, and A. Mitra Long-lived period-doubled edge modes of interacting and disorder-free Floquet spin chains, Communications Physics 5, 10.1038 (2022).
  * (40) A. A. Nizami and A. W. Shrestha, Krylov construction and complexity for driven quantum systems, arXiv:2305.00256
  * (41) P. Suchsland, R. Moessner,and P. W. Claeys, Krylov complexity and Trotter transitions in unitary circuit dynamics, arXiv:2308.03851
  * (42) H.-C. Yeh and A. Mitra A Universal Model of Floquet Operator Krylov Space, arXiv:2311.15116
  * (43) C. Liu, H. Tang, and H. Zhai, Krylov complexity in open quantum systems, Phys. Rev. Res. 5, 033085 (2023)
  * (44) A. Bhattacharya, P. Nandy, P. P. Nath, and H. Sahua, On Krylov complexity in open systems: an approach via bi-Lanczos algorithm, arXiv:2303.04175
  * (45) N. S. Srivatsa and C. von Keyserlingk, The operator growth hypothesis in open quantum systems, arXiv:2310.15376
  * (46) A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a Floquet-space perspective, New J. Phys. 17, 093039 (2015)
  * (47) A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. (Amsterdam) 321, 2 (2006).
  * (48) R. Bianco and R. Resta, Mapping topological order incoordinate space, Phys. Rev. B 84, 241106(R) (2011).
  * (49) F. Song, S. Yao, and Z. Wang, Non-Hermitian Topological Invariants in Real Space, Phys. Rev. Lett. 123, 246801 (2019)


# Continuous-Stage Runge-Kutta approximation to Differential Problems

Pierluigi Amodio Dipartimento di Matematica, Università di Bari, Via Orabona
4, 70125 Bari, Italy. Luigi Brugnano Dipartimento di Matematica e Informatica
“U. Dini”, Università di Firenze, Viale Morgagni 67/A, 50134 Firenze, Italy.
Felice Iavernaro∗

###### Abstract

In recent years, the efficient numerical solution of Hamiltonian problems has
led to the definition of a class of energy-conserving Runge-Kutta methods
named Hamiltonian Boundary Value Methods (HBVMs). Such methods admit an
interesting interpretation in terms of continuous-stage Runge-Kutta methods,
which is here recalled and revisited for general differential problems.

Keywords: Hamiltonian problems; Hamiltonian Boundary Value Methods; HBVMs;
ODE-IVPs; continuous-stage Runge-Kutta methods; continuous-stage Runge-Kutta-
Nyström methods.

MSC: 65L06; 65P10; 65L05.

##  1 Introduction

Continuous-stage Runge-Kutta methods (csRK, hereafter), introduced by Butcher
[34, 35, 36], have been used in recent years as a useful tool for studying
energy-conserving methods for Hamiltonian problems (see, e.g., [22, 25, 38,
41, 42, 44, 43, 45, 46]). In particular, we shall at first consider methods,
within this latter class, having Butcher tableau in the form:

$\begin{array}[]{c|c}c&a(c,\tau)\\\ \hline\cr&b(c)\end{array},$ (1)

where $c,\tau\in[0,1]$, and

$a(c,\tau):[0,1]\times[0,1]\rightarrow\mathbb{R},\qquad b(c):[0,1]\rightarrow
1,$ (2)

are suitable functions defining the method. Hereafter, we shall use the
notation (2), in place of the more commonly used $a_{c\tau},\,b_{c}$, in order
to make clear that these are functions of the respective arguments. For later
use, we shall also denote

$\dot{a}(c,\tau)=\frac{\mathrm{d}}{\mathrm{d}c}a(c,\tau).$ (3)

When used for solving the ODE-IVP 111In the sequel, we shall assume $f$ to be
analytical.

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

the method (1)-(3) provides an approximation

$\dot{\sigma}(ch)=\int_{0}^{1}\dot{a}(c,\tau)f(\sigma(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1],\qquad\sigma(0)=y_{0},$ (5)

to (4) and, therefore,

$\sigma(ch)=y_{0}+h\int_{0}^{1}a(c,\tau)f(\sigma(\tau h))\mathrm{d}\tau,\qquad
c\in[0,1],$ (6)

with the approximation to $y(h)$ given by

$y_{1}=y_{0}+h\int_{0}^{1}b(c)f(\sigma(ch))\mathrm{d}c.$ (7)

In particular, we shall hereafter consider the natural choiche

$b(c)\equiv 1,\qquad c\in[0,1],$ (8)

though different choices have been also considered [41].222Clearly, because of
consistenty, one requires that $\int_{0}^{1}b(c)\mathrm{d}c=1$. The arguments
studied in this paper strictly follows those in [2] (in turn, inspired by [44,
47, 48]), derived from the energy-conserving methods called Hamiltonian
Boundary Value Methods (HBVMs), which have been the subject of many
investigations [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 19,
20, 17, 21, 22, 23, 26, 24, 25, 27, 28, 29, 30, 32, 31, 33]. Here, we consider
also a relevant generalization, w.r.t. the arguments in [2], as explained
below. With this premise, the structure of the paper is as follows: in Section
2 we recall the basic facts concerning the initial value problem (4); in
Section 3 we generalize the approach to the case of second-order problems,
thus obtaining a continuous-stage Runge-Kutta-Nyström (csRKN, hereafter)
method, whose generalization for $k$th-order problems is also sketched;
Section 4 is devoted to derive relevant families of methods for the previous
classes of problems; at last, in Section 5 a few conclusions are given.

##  2 Approximation of ODE-IVPs

Let us consider, at first, the initial value problem (4). As done in [2], for
our analysis we shall use an expansion of the vector field along the
orthonormal Legendre polynomial basis:

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

where, as usual, $\Pi_{i}$ denotes the vector space of polynomials of degree
$i$, and $\delta_{ij}$ is the Kronecker symbol. Consequently, we can rewrite
(4) as

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

with

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

Integrating side-by-side, and imposing the initial condition, one then obtains
that the solution of (10) is formally given by

$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].$ (12)

For $c=1$, by considering that, by virtue of (9),
$\int_{0}^{1}P_{j}(x)\mathrm{d}x=\delta_{j0}$, and taking into account (4),
one obtains:

$y(h)=y_{0}+h\int_{0}^{1}f(y(ch))\mathrm{d}c\,\equiv\,y_{0}+\int_{0}^{h}\dot{y}(t)\mathrm{d}t,$
(13)

i.e., the Fundamental Theorem of the Calculus. Interestingly, by setting

$a_{\infty}(c,\tau)=\sum_{j\geq 0}\int_{0}^{c}P_{j}(x)\mathrm{d}xP_{j}(\tau),$
(14)

one obtains that (10)-(11) and (12)-(13) can be rewritten, respectively, as:

$\displaystyle\dot{y}(ch)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\dot{a}_{\infty}(c,\tau)f(y(\tau
h))\mathrm{d}\tau,\qquad y(0)=y_{0}\in\mathbb{R}^{m},$ $\displaystyle y(ch)$
$\displaystyle=$ $\displaystyle y_{0}+h\int_{0}^{1}a_{\infty}(c,\tau)f(y(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1],$ (15) $\displaystyle y(h)$
$\displaystyle=$ $\displaystyle y_{0}+h\int_{0}^{1}f(y(ch))\mathrm{d}c.$

Consequently,

$\begin{array}[]{c|c}c&a_{\infty}(c,\tau)\\\ \hline\cr&1\end{array},$ (16)

is the csRK “method” providing the exact solution to the problem.

###  2.1 Vector formulation

For later use, we now cast the formulation of (14) in vector form. For this
purpose, let us introduce the infinite vectors,

${\cal P}_{\infty}(c)=\left(\begin{array}[]{c}P_{0}(c)\\\ P_{1}(c)\\\
\vdots\end{array}\right),\qquad{\cal I}_{\infty}(c)=\int_{0}^{c}{\cal
P}_{\infty}(x)\mathrm{d}x,$ (17)

and the matrix,

$X_{\infty}=\left(\begin{array}[]{cccc}\xi_{0}&-\xi_{1}\\\
\xi_{1}&0&-\xi_{2}\\\ &\xi_{2}&\ddots&\ddots\\\
&&\ddots\end{array}\right),\qquad\xi_{i}=\frac{1}{2\sqrt{|4i^{2}-1|}},\qquad
i=0,1,\dots,$ (18)

also recalling that, by virtue of (9), and due to the well-known relations
between the Legendre polynomials and their integrals,

${\cal
I}_{\infty}(c)^{\top}=P_{\infty}(c)^{\top}X_{\infty},~{}\int_{0}^{1}{\cal
P}_{\infty}(\tau){\cal
P}_{\infty}(\tau)^{\top}\mathrm{d}\tau=I,~{}\int_{0}^{1}{\cal
P}_{\infty}(\tau){\cal I}_{\infty}(\tau)^{\top}\mathrm{d}\tau=X_{\infty},$
(19)

with $I$ the identity operator. Consequently, with reference to (14), one has

$a_{\infty}(c,\tau)={\cal I}_{\infty}(c)^{\top}{\cal P}_{\infty}(\tau)={\cal
P}_{\infty}(c)^{\top}X_{\infty}{\cal P}_{\infty}(\tau),$ (20)

and, in particular, one may regard the Butcher tableau, equivalent to (16)

$\begin{array}[]{c|c}c&{\cal P}_{\infty}(c)^{\top}X_{\infty}{\cal
P}_{\infty}(\tau)\\\ \hline\cr&1\end{array},$ (21)

as the corresponding $W$-transformation [39] of the continuous problem.

It is worth mentioning that, by defining the infinite vector (see (11)),

$\bm{\gamma}:=\left(\begin{array}[]{c}\gamma_{0}(y)\\\ \gamma_{1}(y)\\\
\vdots\end{array}\right)\equiv\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes
I_{m}f(y(\tau h))\mathrm{d}\tau,$ (22)

then (12) can be rewritten as

$y(ch)=y_{0}+h{\cal I}_{\infty}(c)^{\top}\otimes I_{m}\bm{\gamma}$

and, consequently, the vector $\bm{\gamma}$ satisfies the equation

$\bm{\gamma}=\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes I_{m}f(y_{0}+h{\cal
I}_{\infty}(\tau)^{\top}\otimes I_{m}\bm{\gamma})\mathrm{d}\tau,$ (23)

with (compare with (15)),

$y(h)=y_{0}+h\gamma_{0}(y).$ (24)

###  2.2 Polynomial approximation

In order to derive a polynomial approximation $\sigma\in\Pi_{s}$ of (15), it
suffices to truncate the infinite series in (14) after $s$ terms:

$a_{s}(c,\tau)=\sum_{j=0}^{s-1}\int_{0}^{c}P_{j}(x)\mathrm{d}xP_{j}(\tau),$
(25)

so that

$\dot{a}_{s}(c,\tau)=\sum_{j=0}^{s-1}P_{j}(c)P_{j}(\tau),$ (26)

and, therefore, (15) is approximated by

$\displaystyle\dot{\sigma}(ch)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\dot{a}_{s}(c,\tau)f(\sigma(\tau
h))\mathrm{d}\tau,\qquad\sigma(0)=y_{0}\in\mathbb{R}^{m},$
$\displaystyle\sigma(ch)$ $\displaystyle=$ $\displaystyle
y_{0}+h\int_{0}^{1}a_{s}(c,\tau)f(\sigma(\tau h))\mathrm{d}\tau,\qquad
c\in[0,1],$ (27) $\displaystyle y_{1}\,:=\,\sigma(h)$ $\displaystyle=$
$\displaystyle y_{0}+h\int_{0}^{1}f(\sigma(ch))\mathrm{d}c.$

###### Remark 1

As it was shown in [2], the csRK method

$\begin{array}[]{c|c}c&a_{s}(c,\tau)\\\ \hline\cr&1\end{array},$ (28)

with $a_{s}(c,\tau)$ given by (25), is equivalent to the energy-preserving
method, named HBVM$(\infty,s)$, introduced in [22].333In particular, when
$s=1$ one obtains the AVF method [42]; for a related approach, see also [38].

The following properties hold true [25]:

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

i.e., the approximation procedure has order $2s$.

###### Remark 2

If $H:\mathbb{R}^{m}\rightarrow\mathbb{R}$, and

$f(y)=J\nabla H(y),\qquad with\qquad J^{\top}=-J,$

then $H(y_{1})=H(y_{0})$. In the case of Hamiltonian problems, $H$ is the
energy of the system. Consequently, the csRK method (28) is energy-conserving,
as is shown in [25, Theorem 3]).

A corresponding vector formulation of the csRK method (28) can be derived by
replacing the infinite vectors and matrix in (17)-(18) with

${\cal P}_{r}(c)=\left(\begin{array}[]{c}P_{0}(c)\\\ \vdots\\\
P_{r-1}(c)\end{array}\right),\quad r=s,s+1,\qquad{\cal
I}_{s}(c)=\int_{0}^{c}{\cal P}_{s}(x)\mathrm{d}x,$ (30)

and the matrices,

$\hat{X}_{s}=\left(\begin{array}[]{cccc}\xi_{0}&-\xi_{1}\\\
\xi_{1}&0&\ddots\\\ &\ddots&\ddots&-\xi_{s-1}\\\ &&\xi_{s-1}&0\\\
\hline\cr&&&\xi_{s}\end{array}\right)\equiv\left(\begin{array}[]{c}X_{s}\\\
\hline\cr 0\dots 0\,\xi_{s}\end{array}\right)\in\mathbb{R}^{s+1\times s},$
(31)

such that

${\cal I}_{s}(c)^{\top}=P_{s+1}(c)^{\top}\hat{X}_{s},\quad\int_{0}^{1}{\cal
P}_{r}(\tau){\cal
P}_{r}(\tau)^{\top}\mathrm{d}\tau=I_{r},\quad\int_{0}^{1}{\cal
P}_{s}(\tau){\cal I}_{s}(\tau)^{\top}\mathrm{d}\tau=X_{s},$ (32)

with $I_{r}\in\mathbb{R}^{r\times r}$ the identity matrix. Consequently, with
reference to (25), one has

$a_{s}(c,\tau)={\cal I}_{s}(c)^{\top}{\cal P}_{s}(\tau)={\cal
P}_{s+1}(c)^{\top}\hat{X}_{s}{\cal P}_{s}(\tau),$ (33)

and, in particular, one may regard the Butcher tableau, equivalent to (28)

$\begin{array}[]{c|c}c&{\cal P}_{s+1}(c)^{\top}\hat{X}_{s}{\cal
P}_{s}(\tau)\\\ \hline\cr&1\end{array},$ (34)

as the corresponding $W$-transformation of a HBVM$(\infty,s)$ method [2].

Further, according to [23], setting the vector

$\bm{\gamma}_{s}:=\left(\begin{array}[]{c}\gamma_{0}(\sigma)\\\ \vdots\\\
\gamma_{s-1}(\sigma)\end{array}\right),$ (35)

with $\gamma_{j}(\sigma)$ defined as in (11), by formally replacing $y$ by
$\sigma$, one has that such vector satisfies the equation (compare with (23))

$\bm{\gamma}_{s}=\int_{0}^{1}{\cal P}_{s}(\tau)\otimes I_{m}f(y_{0}+h{\cal
I}_{s}(\tau)^{\top}\otimes I_{m}\bm{\gamma}_{s})\mathrm{d}\tau,$ (36)

with (compare with (27) and (24)),

$y_{1}=y_{0}+h\gamma_{0}(\sigma).$ (37)

##  3 Approximation of special second-order ODE-IVPs

An interesting particular case is that of special second order problems,
namely problems in the form

$\ddot{y}(t)=f(y(t)),\qquad t\in[0,h],\qquad
y(0)=y_{0},\,\dot{y}(0)=\dot{y}_{0}\in\mathbb{R}^{m},$ (38)

which, in turn, are a special case of the more general problems

$\ddot{y}(t)=f(y(t),\dot{y}(t)),\qquad t\in[0,h],\qquad
y(0)=y_{0},\,\dot{y}(0)=\dot{y}_{0}\in\mathbb{R}^{m}.$ (39)

Let us study, at first, the special problem (38), which is very important in
many applications (as an example, separable Hamiltonian problems are in such a
form), then discussing problem (39). Setting $\dot{y}(t)=p(t)$, and expanding
the right-hand sides along the Legendre basis gives:

$\displaystyle\dot{y}(ch)$ $\displaystyle=$ $\displaystyle{\cal
P}_{\infty}(c)^{\top}\otimes I_{m}\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes
I_{m}p(\tau h)\mathrm{d}\tau,$ (40) $\displaystyle\dot{p}(ch)$
$\displaystyle=$ $\displaystyle{\cal P}_{\infty}(c)^{\top}\otimes
I_{m}\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes I_{m}f(y(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1].$

###### Remark 3

It is worth mentioning that, with reference to (3) and (20), the two previous
equations can be rewritten as:

$\dot{y}(ch)=\int_{0}^{1}\dot{a}_{\infty}(c,\tau)p(\tau
h)\mathrm{d}\tau,\qquad\dot{p}(ch)=\int_{0}^{1}\dot{a}_{\infty}(c,\tau)f(y(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1].$ (41)

Integrating side by side (40), and imposing the initial condition then gives:

$\displaystyle y(ch)$ $\displaystyle=$ $\displaystyle y_{0}+h{\cal
I}_{\infty}(c)^{\top}\otimes I_{m}\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes
I_{m}p(\tau h)\mathrm{d}\tau,$ $\displaystyle p(ch)$ $\displaystyle=$
$\displaystyle\dot{y}_{0}+h{\cal I}_{\infty}(c)^{\top}\otimes
I_{m}\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes I_{m}f(y(\tau
h))\mathrm{d}\tau$ $\displaystyle=$ $\displaystyle\dot{y}_{0}+h{\cal
I}_{\infty}(c)^{\top}\otimes I_{m}\bm{\gamma},\qquad c\in[0,1],$

with the vector $\bm{\gamma}$ formally still given by (22). Substitution of
the second equation in (3) into the first one, taking into account that 444In
general, $e_{i}$ will denote the infinite vector whose $j$th entry is
$\delta_{ij}$.

$\int_{0}^{1}{\cal
P}_{\infty}(\tau)\mathrm{d}\tau=e_{1}\equiv\left(\begin{array}[]{c}1\\\ 0\\\
\vdots\end{array}\right),\qquad and\qquad{\cal I}_{\infty}(c)^{\top}e_{1}=c,$

and considering (19) and again (22), then gives,

$\displaystyle y(ch)$ $\displaystyle=$ $\displaystyle y_{0}+h{\cal
I}_{\infty}(c)^{\top}\otimes I_{m}\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes
I_{m}\left[\dot{y}_{0}+h{\cal I}_{\infty}(\tau)^{\top}\otimes
I_{m}\bm{\gamma}\right]\mathrm{d}\tau$ $\displaystyle=$ $\displaystyle
y_{0}+ch\dot{y}_{0}+h^{2}{\cal I}_{\infty}(c)^{\top}\otimes
I_{m}\int_{0}^{1}{\cal P}_{\infty}(\tau){\cal
I}_{\infty}(\tau)^{\top}\mathrm{d}\tau\otimes I_{m}\bm{\gamma}$
$\displaystyle=$ $\displaystyle y_{0}+ch\dot{y}_{0}+h^{2}{\cal
I}_{\infty}(c)^{\top}X_{\infty}\otimes I_{m}\bm{\gamma}$ $\displaystyle\equiv$
$\displaystyle y_{0}+ch\dot{y}_{0}+h^{2}\int_{0}^{1}\left[{\cal
I}_{\infty}(c)^{\top}X_{\infty}{\cal P}_{\infty}(\tau)\right]\otimes
I_{m}f(y(\tau h))\mathrm{d}\tau$ (43)

Setting

$\bar{a}_{\infty}(c,\tau)={\cal I}_{\infty}(c)^{\top}X_{\infty}{\cal
P}_{\infty}(\tau)\equiv{\cal P}_{\infty}(c)^{\top}X_{\infty}^{2}{\cal
P}_{\infty}(\tau),$ (44)

one has that (3) can be rewritten as:

$y(ch)=y_{0}+ch\dot{y}_{0}+h^{2}\int_{0}^{1}\bar{a}_{\infty}(c,\tau)f(y(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1].$ (45)

###### Remark 4

We observe that, from the second equation in (40) and (41), and considering
(20), one derives:

$p(ch)=\dot{y}_{0}+h\int_{0}^{1}a_{\infty}(c,\tau)f(y(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1].$ (46)

Moreover, it is worth mentioning that, see (44), (14), and (19),

$\bar{a}_{\infty}(c,\tau)=\int_{0}^{1}{\cal I}_{\infty}(c)^{\top}{\cal
P}_{\infty}(\xi){\cal I}_{\infty}(\xi)^{\top}{\cal
P}_{\infty}(\tau)\mathrm{d}\xi\equiv\int_{0}^{1}a_{\infty}(c,\xi)a_{\infty}(\xi,\tau)\mathrm{d}\xi.$
(47)

From (40) and (45) one obtains that the values at $h$ will be given by (see
(18)):

$\dot{y}(h)~{}\equiv~{}p(h)~{}=~{}\dot{y}_{0}+he_{1}\top\otimes
I_{m}\bm{\gamma}\,=\,\dot{y}_{0}+h\gamma_{0}(y)\,\equiv\,\dot{y}_{0}+h\int_{0}^{1}f(y(ch))\mathrm{d}c,\\\
$ (48)

and

$\displaystyle y(h)$ $\displaystyle=$ $\displaystyle
y_{0}+h\dot{y}_{0}+h^{2}e_{1}^{\top}X_{\infty}\int_{0}^{1}{\cal
P}_{\infty}(c)\otimes I_{m}f(y(ch)\mathrm{d}c$ (49) $\displaystyle=$
$\displaystyle y_{0}+h\dot{y}_{0}+h^{2}e_{1}^{\top}X_{\infty}\int_{0}^{1}{\cal
P}_{\infty}(c)\otimes I_{m}f(y(ch)\mathrm{d}c$ $\displaystyle=$ $\displaystyle
y_{0}+h\dot{y}_{0}+h^{2}\int_{0}^{1}(\xi_{0}e_{1}-\xi_{1}e_{2})^{\top}{\cal
P}_{\infty}(c)\otimes I_{m}f(y(ch)\mathrm{d}c$ $\displaystyle=$ $\displaystyle
y_{0}+h\dot{y}_{0}+h^{2}\int_{0}^{1}[\xi_{0}P_{0}(c)-\xi_{1}P_{1}(c)]f(y(ch)\mathrm{d}c$
$\displaystyle\equiv$ $\displaystyle
y_{0}+h\dot{y}_{0}+h^{2}\int_{0}^{1}(1-c)f(y(ch)\mathrm{d}c.$

In other words, the exact solution of problem (38) is generated by the
following csRKN “method”:

$\begin{array}[]{c|c}c&\bar{a}_{\infty}(c,\tau)\\\ \hline\cr&1-c\\\
\hline\cr&1\end{array}\quad\equiv\quad\begin{array}[]{c|c}c&{\cal
I}_{\infty}(c)^{\top}X_{\infty}{\cal P}_{\infty}(\tau)\\\ \hline\cr&1-c\\\
\hline\cr&1\end{array}\quad\equiv\quad\begin{array}[]{c|c}c&{\cal
P}_{\infty}(c)^{\top}X_{\infty}^{2}{\cal P}_{\infty}(\tau)\\\ \hline\cr&1-c\\\
\hline\cr&1\end{array}.$ (50)

###  3.1 Vector formulation

An interesting alternative formulation of the method (50), akin to (23)-(24)
for first order problems, can be derived by combining (22) and (43):

$\bm{\gamma}=\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes
I_{m}\,f\left(y_{0}+\tau h\dot{y}_{0}+h^{2}{\cal
I}_{\infty}(\tau)^{\top}X_{\infty}\otimes
I_{m}\bm{\gamma}\right)\mathrm{d}\tau,$ (51)

with values at $t=h$ given by (compare with (48)-(49)):

$\dot{y}(h)=\dot{y}_{0}+h\gamma_{0}(y),\qquad
y(h)=y_{0}+h\dot{y}_{0}+h^{2}\left(\xi_{0}\gamma_{0}(y)-\xi_{1}\gamma_{1}(y)\right).$
(52)

###  3.2 The case of the general problem (39)

The arguments used above, can be extended to cope with problem (39) in a
straightforward way. In fact, by following similar steps as before, (45)-(46)
now becomes:

$\displaystyle y(ch)$ $\displaystyle=$ $\displaystyle
y_{0}+ch\dot{y}_{0}+h^{2}\int_{0}^{1}\bar{a}_{\infty}(c,\tau)f(y(\tau
h),\dot{y}(\tau h))\mathrm{d}\tau,$ $\displaystyle\dot{y}(ch)$
$\displaystyle=$
$\displaystyle\dot{y}_{0}+h\int_{0}^{1}a_{\infty}(c,\tau)f(y(\tau
h),\dot{y}(\tau h))\mathrm{d}\tau,\qquad c\in[0,1],$ (53)

with the values at $h$ given by

$\dot{y}_{0}+h\int_{0}^{1}f(y(ch),\dot{y}(ch)))\mathrm{d}c,\qquad
y(h)=y_{0}+h\dot{y}_{0}+h^{2}\int_{0}^{1}(1-c)f(y(ch),\dot{y}(ch))\mathrm{d}c.$
(54)

Consequently, we obtain the general csRKN method

$\displaystyle\begin{array}[]{c|c|c}c&\bar{a}_{\infty}(c,\tau)&a_{\infty}(c,\tau)\\\
\hline\cr&1-c&1\end{array}~{}\equiv~{}\begin{array}[]{c|c|c}c&{\cal
I}_{\infty}(c)^{\top}X_{\infty}{\cal P}_{\infty}(\tau)&{\cal
I}_{\infty}(c)^{\top}{\cal P}_{\infty}(\tau)\\\ \hline\cr&1-c&1\end{array}$
(59) $\displaystyle\equiv$ $\displaystyle\begin{array}[]{c|c|c}c&{\cal
P}_{\infty}(c)^{\top}X_{\infty}^{2}{\cal P}_{\infty}(\tau)&{\cal
P}_{\infty}(c)^{\top}X_{\infty}{\cal P}_{\infty}(\tau)\\\
\hline\cr&1-c&1\end{array},\hskip 85.35826pt$ (62)

in place of (50).

The bad news is that now the system of equations (53) has a doubled size,
w.r.t. (45). On the other hand, the good news is that, upon modifying the
definition of the vector $\bm{\gamma}$ in (22) as follows,

$\bm{\gamma}:=\left(\begin{array}[]{c}\gamma_{0}(y,\dot{y})\\\
\gamma_{1}(y,\dot{y})\\\ \vdots\end{array}\right)\equiv\int_{0}^{1}{\cal
P}_{\infty}(\tau)\otimes I_{m}f(y(\tau h),\dot{y}(\tau h))\mathrm{d}\tau,$
(63)

one has that the equations (53) can be rewritten as

$\displaystyle y(ch)$ $\displaystyle=$ $\displaystyle
y_{0}+ch\dot{y}_{0}+h^{2}{\cal I}_{\infty}(c)^{\top}X_{\infty}\otimes
I_{m}\bm{\gamma},$ $\displaystyle\dot{y}(ch)$ $\displaystyle=$
$\displaystyle\dot{y}_{0}+h{\cal I}_{\infty}(c)^{\top}\otimes
I_{m}\bm{\gamma},\qquad c\in[0,1].$

Consequently, we obtain again a single equation for the vector $\bm{\gamma}$
defined in (63),

$\displaystyle\bm{\gamma}=\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes I_{m}$
$\displaystyle f\left(y_{0}+\tau h\dot{y}_{0}+h^{2}{\cal
I}_{\infty}(\tau)^{\top}X_{\infty}\otimes
I_{m}\bm{\gamma},\,\dot{y}_{0}+h{\cal I}_{\infty}(c)^{\top}\otimes
I_{m}\bm{\gamma}\right)\mathrm{d}\tau.$

Similarly, (48)-(49) respectively become:

$\dot{y}(h)=\dot{y}_{0}+h\gamma_{0}(y,\dot{y}),\quad
y(h)=y_{0}+h\dot{y}_{0}+h^{2}\left(\xi_{0}\gamma_{0}(y,\dot{y})-\xi_{1}\gamma_{1}(y,\dot{y})\right).$
(65)

Remarkably enough, the equations (3.2)-(65) are very similar to (51)-(52).

###  3.3 Polynomial approximation

Polynomial approximations of degree $s$, $\sigma(ch)\approx y(ch)$ and
$\sigma_{1}(ch)\approx p(ch)$, can be obtained by formal substitution of the
matrices ${\cal P}_{\infty}(c),{\cal I}_{\infty}(c)$, and $X_{\infty}$ in
(40)–(49) with the corresponding finite ones, ${\cal P}_{s}(c)$, ${\cal
I}_{s}(c)\equiv{\cal P}_{s+1}(c)\hat{X}_{s}$, and $X_{s}$ defined in
(30)-(31). Consequently, following similar steps as above, one obtains:

$\displaystyle\dot{\sigma}(ch)$ $\displaystyle=$ $\displaystyle{\cal
P}_{s}(c)^{\top}\otimes I_{m}\int_{0}^{1}{\cal P}_{s}(\tau)\otimes
I_{m}\sigma_{1}(\tau h)\mathrm{d}\tau,$ (66)
$\displaystyle\dot{\sigma}_{1}(ch)$ $\displaystyle=$ $\displaystyle{\cal
P}_{s}(c)^{\top}\otimes I_{m}\int_{0}^{1}{\cal P}_{s}(\tau)\otimes
I_{m}f(\sigma(\tau h))\mathrm{d}\tau,\qquad c\in[0,1],$

in place of (40). Similarly, setting (in place of (44))

$\bar{a}_{s}(c,\tau)={\cal I}_{s}(c)^{\top}X_{s}{\cal P}_{s}(\tau)\equiv{\cal
P}_{s+1}(c)^{\top}\hat{X}_{s}X_{s}{\cal P}_{s}(\tau),$ (67)

from (66) one derives:

$\sigma(ch)=y_{0}+ch\dot{y}_{0}+h^{2}\int_{0}^{1}\bar{a}_{s}(c,\tau)f(\sigma(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1],$ (68)

with the new approximations

$y_{1}:=\sigma(h)\qquad and\qquad\dot{y}_{1}:=\sigma_{1}(h)$

given by

$\dot{y}_{1}=\dot{y}_{0}+h\int_{0}^{1}f(\sigma(ch))\mathrm{d}c,\qquad
y_{1}=y_{0}+h\dot{y}_{0}+h^{2}\int_{0}^{1}(1-c)f(\sigma(ch))\mathrm{d}c.$ (69)

###### Remark 5

Similarly as (41), one derives:

$\dot{\sigma}(ch)=\int_{0}^{1}\dot{a}_{s}(c,\tau)\sigma_{1}(\tau
h)\mathrm{d}\tau,\quad\dot{\sigma}_{1}(ch)=\int_{0}^{1}\dot{a}_{s}(c,\tau)f(\sigma(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1].$ (70)

Moreover, (compare with (46)), from the second equation in (66), and
considering (33), one obtains:

$\sigma_{1}(ch)=\dot{y}_{0}+h\int_{0}^{1}a_{s}(c,\tau)f(\sigma(\tau
h))\mathrm{d}\tau,\qquad c\in[0,1].$ (71)

At last, (compare with (47)), from (67), (33), and (32), one has:

$\bar{a}_{s}(c,\tau)=\int_{0}^{1}{\cal I}_{s}(c)^{\top}{\cal P}_{s}(\xi){\cal
I}_{s}(\xi)^{\top}{\cal
P}_{s}(\tau)\mathrm{d}\xi\equiv\int_{0}^{1}a_{s}(c,\xi)a_{s}(\xi,\tau)\mathrm{d}\xi.$
(72)

Summing all up, through (67)–(69) we have defined the following csRKN method:

$\begin{array}[]{c|c}c&\bar{a}_{s}(c,\tau)\\\ \hline\cr&1-c\\\
\hline\cr&1\end{array}\quad\equiv\quad\begin{array}[]{c|c}c&{\cal
I}_{s}(c)^{\top}X_{s}{\cal P}_{s}(\tau)\\\ \hline\cr&1-c\\\
\hline\cr&1\end{array}\quad\equiv\quad\begin{array}[]{c|c}c&{\cal
P}_{s+1}(c)^{\top}\hat{X}_{s}X_{s}{\cal P}_{s}(\tau)\\\ \hline\cr&1-c\\\
\hline\cr&1\end{array}.$ (73)

###### Remark 6

This latter method is equivalent to the csRKN method obtained by applying the
HBVM$(\infty,s)$ method to the special second-order problem (38) [2].

An alternative formulation of the method (73) can be obtained by repeating
arguments similar to those used in Section 3.1. In fact, by setting the vector
(compare with (35)–(37))

$\bm{\gamma}_{s}:=\left(\begin{array}[]{c}\gamma_{j}(\sigma)\\\ \vdots\\\
\gamma_{s-1}(\sigma)\end{array}\right)\equiv\int_{0}^{1}{\cal
P}_{s}(\tau)\otimes I_{m}f(\sigma(\tau h))\mathrm{d}\tau,$ (74)

by virtue of (67)-(68) such a vector satisfies the equation

$\bm{\gamma}_{s}=\int_{0}^{1}{\cal P}_{s}(\tau)\otimes
I_{m}\,f\left(y_{0}+\tau h\dot{y}_{0}+h^{2}{\cal
I}_{s}(\tau)^{\top}X_{s}\otimes I_{m}\bm{\gamma}_{s}\right)\mathrm{d}\tau,$
(75)

with the new approximations given by (compare with (52)):

$\dot{y}_{1}=\dot{y}_{0}+h\gamma_{0}(\sigma),\qquad
y_{1}=y_{0}+h\dot{y}_{0}+h^{2}\left(\xi_{0}\gamma_{0}(\sigma)-\xi_{1}\gamma_{1}(\sigma)\right).$
(76)

For the general problem (39), following similar steps as above, the
generalized csRKN method (62) becomes, with reference to (33):

$\displaystyle\begin{array}[]{c|c|c}c&\bar{a}_{s}(c,\tau)&a_{s}(c,\tau)\\\
\hline\cr&1-c&1\end{array}~{}\equiv~{}\begin{array}[]{c|c|c}c&{\cal
I}_{s}(c)^{\top}X_{s}{\cal P}_{s}(\tau)&{\cal I}_{s}(c)^{\top}{\cal
P}_{s}(\tau)\\\ \hline\cr&1-c&1\end{array}$ (81) $\displaystyle\equiv$
$\displaystyle\begin{array}[]{c|c|c}c&{\cal
P}_{s+1}(c)^{\top}\hat{X}_{s}X_{s}{\cal P}_{s}(\tau)&{\cal
P}_{s+1}(c)^{\top}\hat{X}_{s}X_{s}{\cal P}_{s}(\tau)\\\
\hline\cr&1-c&1\end{array}.\hskip 28.45274pt$ (84)

Also in this case, we can derive a vector formulation, similar to (74)-(75).
In fact, by defining the vector

$\bm{\gamma}_{s}:=\left(\begin{array}[]{c}\gamma_{j}(\sigma,\sigma_{1})\\\
\vdots\\\
\gamma_{s-1}(\sigma,\sigma_{1})\end{array}\right)\equiv\int_{0}^{1}{\cal
P}_{s}(\tau)\otimes I_{m}f(\sigma(\tau h),\sigma_{1}(\tau h))\mathrm{d}\tau,$
(85)

it satisfies the equation

$\displaystyle\bm{\gamma}_{s}=\int_{0}^{1}{\cal P}_{s}(\tau)\otimes I_{m}\,$
$\displaystyle f\left(y_{0}+\tau h\dot{y}_{0}+h^{2}{\cal
I}_{s}(\tau)^{\top}X_{s}\otimes I_{m}\bm{\gamma}_{s},\,\dot{y}_{0}+h{\cal
I}_{s}(\tau)^{\top}\otimes I_{m}\bm{\gamma}_{s}\right)\mathrm{d}\tau,$

with the new approximations given by:

$\dot{y}_{1}=\dot{y}_{0}+h\gamma_{0}(\sigma,\sigma_{1}),\qquad
y_{1}=y_{0}+h\dot{y}_{0}+h^{2}\left(\xi_{0}\gamma_{0}(\sigma,\sigma_{1})-\xi_{1}\gamma_{1}(\sigma,\sigma_{1})\right).$
(87)

###### Remark 7

By comparing the discrete problems (75) and (3.3), one realizes that they have
the same dimension, independently of the fact that the latter one solves the
general problem (39). This fact is even more striking since, as we are going
to sketch in the next section, this will be the case for a general $k$th order
ODE-IVP.

###  3.4 Approximation of general $k$th-order ODE-IVPs

Let us consider the case of a general $k$th order problem,

$\displaystyle y^{(k)}(ch)$ $\displaystyle=$ $\displaystyle
f(y(ch),y^{(1)}(ch),\dots,y^{(k-1)}(ch)),\qquad c\in[0,1],$ $\displaystyle
y^{(i)}(0)$ $\displaystyle=$ $\displaystyle
y^{(i)}_{0}\in\mathbb{R}^{m},\qquad i=0,\dots,k-1.$ (88)

Repeating similar steps as above, by defining the infinite vector
555Hereafter, for sake of brevity, we shall skip the arguments of the Fourier
coefficients $\gamma_{i}$.

$\bm{\gamma}\equiv\left(\begin{array}[]{c}\gamma_{0}\\\ \gamma_{1}\\\
\vdots\end{array}\right):=\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes
I_{m}f(y(\tau h),\dots,y^{(k-1)}(\tau h))\mathrm{d}\tau,$ (89)

one has that it satisfies that equation

$\displaystyle\bm{\gamma}=\int_{0}^{1}{\cal P}_{\infty}(\tau)\otimes I_{m}$
$\displaystyle f\left(\sum_{i=0}^{k-1}\frac{(\tau
h)^{i}}{i!}y_{0}^{(i)}+h^{k}{\cal
I}_{\infty}(\tau)^{\top}X_{\infty}^{k-1}\otimes
I_{m}\bm{\gamma},\dots,y_{0}^{(k-1)}+h{\cal I}_{\infty}(\tau)^{\top}\otimes
I_{m}\bm{\gamma}\right)\mathrm{d}\tau,$

with the values at $t=h$ given by:

$y^{(i)}(h)=\sum_{j=0}^{k-1-i}\frac{h^{j}}{j!}y^{(j+i)}_{0}+h^{k-i}\sum_{j=0}^{k-1-i}b_{j}^{(k-1-i)}\gamma_{j},\qquad
i=0,\dots,k-1,$ (91)

$b_{j}^{(k-1-i)}$, $j=0,\dots,k-1-i$, being the $(j+1)$st entry on the first
row of the matrix

$X_{\infty}^{k-1-i},\qquad i=0,\dots,k-1.$

A polynomial approximation of degree $s$ (resulting, as usual, into an order
$2s$ method) can be derived by formally substituting, in the equation (3.4),
${\cal P}_{\infty}(\tau)$ and ${\cal I}_{\infty}(\tau)$ with ${\cal
P}_{s}(\tau)$ and ${\cal I}_{s}(\tau)$, respectively,666In so doing, the
vector $\bm{\gamma}$ now belongs to $\mathbb{R}^{sm}$. with the new
approximations $y_{1}^{(i)}\approx y^{(i)}(h)$, given by:

$y^{(i)}_{1}\,:=\,\sum_{j=0}^{k-1-i}\frac{h^{j}}{j!}y^{(j+i)}_{0}+h^{k-i}\sum_{j=0}^{k-1-i}b_{j}^{(k-1-i)}\gamma_{j},\qquad
i=0,\dots,k-1,$ (92)

$b_{j}^{(k-1-i)}$, $j=0,\dots,k-1-i$, being now the $(j+1)$st entry on the
first row of the matrix

$X_{s}^{k-1-i},\qquad i=0,\dots,k-1.$

##  4 Discretization

“As is well known, even many relatively simple
integrals cannot be expressed in finite terms of
elementary functions, and thus must be evalu-
ated by numerical methods.”

Dahlquist and Björk [37, page. 521]

As is clear, the integrals involved in the Fourier coefficients (36) need to
be approximated by using a suitable quadrature rule, which we choose as the
interpolatory Gauss-Legendre quadrature of order $2k$, with $k\geq s$, having
weights and abscissae $(b_{i},c_{i})$, $i=1,\dots,k$. In so doing the vector
of the Fourier coefficients (35)-(36) becomes

$\hat{\bm{\gamma}}_{s}:=\left(\begin{array}[]{c}\hat{\gamma}_{0}\\\ \vdots\\\
\hat{\gamma}_{s-1}\end{array}\right),$ (93)

satisfying the equation (compare with (36))

$\hat{\bm{\gamma}}_{s}=\sum_{j=1}^{k}b_{j}{\cal P}_{s}(c_{j})\otimes
I_{m}f(y_{0}+h{\cal I}_{s}(c_{j}))\otimes I_{m}\hat{\bm{\gamma}}_{s})$ (94)

with (compare with (37)),

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

We observe that, by introducing the matrices (see (30))

$\displaystyle{\cal P}_{r}:=\left(\begin{array}[]{c}{\cal
P}_{r}(c_{1})^{\top}\\\ \vdots\\\ {\cal
P}_{r}(c_{k})^{\top}\end{array}\right)\in\mathbb{R}^{k\times r},\quad{\cal
I}_{s}:=\left(\begin{array}[]{c}{\cal I}_{s}(c_{1})^{\top}\\\ \vdots\\\ {\cal
I}_{s}(c_{k})^{\top}\end{array}\right)\in\mathbb{R}^{k\times s},$ (102)
$\displaystyle\Omega=\left(\begin{array}[]{ccc}b_{1}\\\ &\ddots\\\
&&b_{k}\end{array}\right),$ (106)

one has, with reference to (31),

${\cal I}_{s}={\cal P}_{s+1}\hat{X}_{s},\qquad{\cal P}_{s}^{\top}\Omega{\cal
P}_{s}=I_{s},\qquad{\cal P}_{s}^{\top}\Omega{\cal I}_{s}=X_{s}.$ (107)

Further, by also introducing the vector
$\bm{e}=(1,\dots,1)^{\top}\in\mathbb{R}^{k}$, one obtains that (94) can be
rewritten as

$\hat{\bm{\gamma}}_{s}={\cal P}_{s}^{\top}\Omega\otimes I_{m}f(\bm{e}\otimes
y_{0}+h{\cal I}_{s}\otimes I_{m}\hat{\bm{\gamma}}_{s}).$ (108)

Setting

$Y\equiv\left(\begin{array}[]{c}Y_{1}\\\ \vdots\\\
Y_{k}\end{array}\right):=\bm{e}\otimes y_{0}+h{\cal I}_{s}\otimes
I_{m}\hat{\bm{\gamma}}_{s}$ (109)

the vector of the stages of the corresponding Runge-Kutta method, from
(108)-(109) one obtains the stage equation

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

with the new approximation

$y_{1}=y_{0}+h\sum_{i=1}^{k}b_{i}f(Y_{i}).$ (111)

Summing all up, (108)–(111) define the $k$-stage Runge-Kutta method

$\begin{array}[]{c|c}\bm{c}&{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\\\
\hline\cr&\bm{b}^{\top}\end{array}~{},\qquad\bm{b}=(b_{1},\dots,b_{k})^{\top},\qquad\bm{c}=(c_{1},\dots,c_{k})^{\top},$
(112)

named HBVM$(k,s)$ [15, 16, 22]. It is worth mentioning that the Butcher matrix
of the method, with reference to the coefficients of the csRK (33), is given
by:

${\cal I}_{s}{\cal
P}_{s}^{\top}\Omega\equiv\left(\begin{array}[]{c}b_{j}a_{s}(c_{i},c_{j})\end{array}\right)\in\mathbb{R}^{k\times
k}.$

In particular, when $k=s$ one obtains the $s$-stage Gauss collocation method.
However, the use of values $k>s$ (and even $k\gg s$) is useful, in view of
deriving energy-conserving methods for Hamiltonian systems [22, 25, 15, 16].

###### Remark 8

As is clear, the formulation (94)-(95) is computationally more effective than
(110)-(111), having the former (block) dimension $s$, independently of $k$,
which is the dimension of the latter formulation. As observed in [23], this
allows the use of relatively large values of $k$, without increasing too much
the computational cost.

Similar arguments can be repeated in the case of the polynomial approximations
for problems (38), (39), and (88): we here sketch only those concerning the
csRKN (73), providing the correct implementation for the HBVM$(k,s)$ method
for the special second-order problem (38). By formally using the same
approximate Fourier coefficients (93) in place of (74), one has that (75) is
replaced by the following discrete counterpart,

$\hat{\bm{\gamma}}_{s}=\sum_{j=1}^{k}b_{j}{\cal P}_{s}(c_{j})\otimes
I_{m}f(y_{0}+c_{j}h\dot{y}_{0}+h^{2}{\cal I}_{s}(c_{j}))X_{s}\otimes
I_{m}\hat{\bm{\gamma}}_{s})$ (113)

with the new approximations (compare with (76)) given by

$\dot{y}_{1}=\dot{y}_{0}+h\hat{\gamma}_{0},\qquad
y_{1}=y_{0}+h\dot{y}_{0}+h^{2}\left(\xi_{0}\hat{\gamma}_{0}-\xi_{1}\hat{\gamma}_{1}\right).$
(114)

Similarly as in the first order case, (113) can be rewritten as

$\hat{\bm{\gamma}}_{s}={\cal P}_{s}^{\top}\Omega\otimes I_{m}f(\bm{e}\otimes
y_{0}+h\bm{c}\otimes\dot{y}_{0}+h^{2}{\cal I}_{s}X_{s}\otimes
I_{m}\hat{\bm{\gamma}}_{s}),$ (115)

where $\bm{c}$ is the vector of the abscissae. Again, the vector

$Y\equiv\left(\begin{array}[]{c}Y_{1}\\\ \vdots\\\
Y_{k}\end{array}\right):=\bm{e}\otimes
y_{0}+h\bm{c}\otimes\dot{y}_{0}+h^{2}{\cal I}_{s}X_{s}\otimes
I_{m}\hat{\bm{\gamma}}_{s}$ (116)

is the stage vector of a $k$-stage RKN method. In particular, from (115) and
(116) one obtains

$Y=\bm{e}\otimes y_{0}+h\bm{c}\otimes\dot{y}_{0}+h^{2}{\cal I}_{s}X_{s}{\cal
P}_{s}^{\top}\Omega\otimes I_{m}f(Y),$ (117)

and some algebra shows that the new approximations are given by

$\dot{y}_{1}=\dot{y}_{0}+h\sum_{i=1}^{k}b_{i}f(Y_{i}),\qquad
y_{1}=y_{0}+h\dot{y}_{0}+h^{2}\sum_{i=1}^{k}b_{i}(1-c_{i})f(Y_{i}).$ (118)

Consequently, we are speaking about the following RKN method:777Here, $\circ$
denotes the Hadamard, i.e. componentwise, product.

$\begin{array}[]{c|c}\bm{c}&{\cal I}_{s}X_{s}{\cal P}_{s}^{\top}\Omega\\\
\hline\cr&[\bm{b}\circ(\bm{e}-\bm{c})]^{\top}\\\
\hline\cr&\bm{b}^{\top}\end{array}.$ (119)

The Butcher tableau (119) defines the RKN formulation of a HBVM$(k,s)$ method
[2, 23, 15]. Is is worth mentioning that, with reference to (67) and (72)

${\cal I}_{s}X_{s}{\cal
P}_{s}^{\top}\Omega\equiv\left(\begin{array}[]{c}b_{j}\bar{a}_{s}(c_{i},c_{j})\end{array}\right)\equiv\left(\begin{array}[]{c}b_{j}\int_{0}^{1}a_{s}(c_{i},\tau)a_{s}(\tau,c_{j})\mathrm{d}\tau\end{array}\right)\in\mathbb{R}^{k\times
k}.$

###### Remark 9

As observed in Remark 8, also in this case, in consideration that $k\geq s$,
the formulation (114)-(115) of the method is much more efficient than the
usual one, given by (117)-(118), also in view of the use of values of $k\gg
s$, as in the case of separable Hamiltonian problems.

We conclude this section by recalling that:

  * •

very efficient iterative procedures exist for solving the discrete problems
(108) and (115) generated by a HBVM method (see, e.g., [23]);888We mention
that a state-of-art Matlab code is available at the website of the book [15].

  * •

when HBVMs are used as spectral methods in time, i.e., choosing values of $s$
and $k>s$ so that no accuracy improvement can be obtained, for the considered
finite-precision arithmetic and timestep $h$ used [3, 18, 29], then there is
no practical difference between the discrete methods and their continuous-
stage counterparts.

##  5 Conclusions

In this paper we have recalled the basic facts, also reporting some new
insight, on the energy-conserving class of Runge-Kutta methods named HBVMs.
Such method have been here studied within the framework of continuous-stage
Runge-Kutta methods. The extension to second-order problems has been also
studied, providing a natural continuous-stage Runge-Kutta-Nyström formulation
of the methods. Further, also the extension to general $k$th-order problems
has been sketched. The relation with the fully discrete methods has been also
recalled, thus showing the usefulness of using such a framework to study the
fully discrete methods.

#### Conflict of interests.

The authors declare no conflict of interest.

## References

  * [1] Amodio, P.; Brugnano, L.; Iavernaro, F. Energy-conserving methods for Hamiltonian Boundary Value Problems and applications in astrodynamics. Adv. Comput. Math. 2015, 41, 881–905; https://doi.org/10.1007/s10444-014-9390-z
  * [2] Amodio, P.; Brugnano, L.; Iavernaro, F. A note on the continuous-stage Runge-Kutta-(Nyström) formulation of Hamiltonian Boundary Value Methods (HBVMs). Appl. Math. Comput. 2019, 363, 124634; https://doi.org/10.1016/j.amc.2019.124634
  * [3] Amodio, P.; Brugnano, L.; Iavernaro, F. Analysis of Spectral Hamiltonian Boundary Value Methods (SHBVMs) for the numerical solution of ODE problems. Numer. Algorithms 2020, 83, 1489–1508; https://doi.org/10.1007/s11075-019-00733-7
  * [4] Amodio, P.; Brugnano, L.; Iavernaro, F. Arbitrarily high-order energy-conserving methods for Poisson problems. Numer. Algorithms 2022; https://doi.org/10.1007/s11075-022-01285-z
  * [5] Barletti, L.; Brugnano, L.; Frasca-Caccia, G.; Iavernaro, F. Energy-conserving methods for the nonlinear Schroedinger equation. Appl. Math. Comput. 2018, 318, 3–18; https://doi.org/10.1016/j.amc.2017.04.018
  * [6] Barletti, L.; Brugnano, L.; Tang, Y.; Zhu, B. Spectrally accurate space-time solution of Manakov systems. J. Comput. Appl. Math. 2020, 377, 112918; https://doi.org/10.1016/j.cam.2020.112918
  * [7] Brugnano, L.; Calvo, M.; Montijano, J.I.; Ràndez, L. Energy preserving methods for Poisson systems. J. Comput. Appl. Math. 2012, 236, 3890–3904; https://doi.org/10.1016/j.cam.2012.02.033
  * [8] Brugnano, L.; Frasca-Caccia, G.; Iavernaro, F. Efficient implementation of geometric integrators for separable Hamiltonian problems. AIP Conf. Proc. 2013, 1588, 734–737; http://doi.org/10.1063/1.4825598
  * [9] Brugnano, L.; Frasca-Caccia, G.; Iavernaro, F. Hamiltonian Boundary Value Methods (HBVMs) and their efficient implementation. Mathematics in Engineering, Science and Aerospace, 2014 5, 4, 343–411.
  * [10] Brugnano, L.; Frasca-Caccia, G.; Iavernaro, F. Efficient implementation of Gauss collocation and Hamiltonian Boundary Value Methods. Numer. Algorithms 2014, 65, 633–650; https://doi.org/10.1007/s11075-014-9825-0
  * [11] Brugnano, L.; Frasca-Caccia, G.; Iavernaro, F. Energy conservation issues in the numerical solution of the semilinear wave equation. Appl. Math. Comput. 2015, 270, 842–870; https://doi.org/10.1016/j.amc.2015.08.078
  * [12] Brugnano, L.; Frasca-Caccia, G.; Iavernaro, F. Line Integral Solution of Hamiltonian PDEs. Mathematics, 2019, 7(3), 275; https://doi.org/10.3390/math7030275
  * [13] Brugnano, L.; Gurioli, G.; Iavernaro, F.; Weinmüller, E.B. Line integral solution of Hamiltonian systems with holonomic constraints. Appl. Numer. Math. 2018, 127, 56–77; https://doi.org/10.1016/j.apnum.2017.12.014
  * [14] Brugnano, L.; Gurioli, G.; Zhang, C. Spectrally accurate energy-preserving methods for the numerical solution of the ”Good” Boussinesq equation. Numer. Methods Partial Differential Equations 2019, 35, 1343–1362; https://doi.org/10.1002/num.22353
  * [15] Brugnano, L.; Iavernaro, F. Line Integral Methods for Conservative Problems; Chapman et Hall/CRC: Boca Raton, FL, USA, 2016. http://web.math.unifi.it/users/brugnano/LIMbook/
  * [16] Brugnano, L.; Iavernaro, F. Line Integral Solution of Differential Problems. Axioms 2018 , 7(2), 36; https://doi.org/10.3390/axioms7020036
  * [17] Brugnano, L.; Iavernaro, F. Line Integral Methods which preserve all invariants of conservative problems. J. Comput. Appl. Math. 2012, 236, 3905–3919; https://doi.org/10.1016/j.cam.2012.03.026
  * [18] Brugnano, L.; Iavernaro, F.; Montijano, J.I.; Rández, L. Spectrally accurate space-time solution of Hamiltonian PDEs. Numer. Algorithms 2019, 81, 1183–1202; https://doi.org/10.1007/s11075-018-0586-z
  * [19] Brugnano, L.; Iavernaro, F.; Susca, T. Hamiltonian BVMs (HBVMs): implementation details and applications. AIP Conf. Proc., 2009, 1168, 723–726; https://doi.org/10.1063/1.3241568
  * [20] Brugnano, L.; Iavernaro, F.; Susca, T. Numerical comparisons between Gauss-Legendre methods and Hamiltonian BVMs defined over Gauss points. Monogr. Real Acad. Cienc. Zaragoza 2010, 33, 95–112.
  * [21] Brugnano, L.; Iavernaro, F.; Trigiante, D. Hamiltonian BVMs (HBVMs): A family of “drift-free” methods for integrating polynomial Hamiltonian systems. AIP Conf. Proc. 2009, 1168, 715–718; https://doi.org/10.1063/1.3241566
  * [22] Brugnano, L.; Iavernaro, F.; Trigiante, D. Hamiltonian Boundary Value Methods (Energy Preserving Discrete Line Integral Methods). JNAIAM J. Numer. Anal. Ind. Appl. Math. 2010, 5, 17–37;
  * [23] Brugnano, L.; Iavernaro, F.; Trigiante, D. A note on the efficient implementation of Hamiltonian BVMs. J. Comput. Appl. Math. 2011, 236, 375–383; https://doi.org/10.1016/j.cam.2011.07.022
  * [24] Brugnano, L.; Iavernaro, F.; Trigiante, D. The lack of continuity and the role of Infinite and infinitesimal in numerical methods for ODEs: The case of symplecticity. Appl. Math. Comput. 2012, 218, 8053–8063; https://doi.org/10.1016/j.amc.2011.03.022
  * [25] Brugnano, L.; Iavernaro, F.; Trigiante, D. A simple framework for the derivation and analysis of effective one-step methods for ODEs. Appl. Math. Comput. 2012, 218, 8475–8485; https://doi.org/10.1016/j.amc.2012.01.074
  * [26] Brugnano, L.; Iavernaro, F.; Trigiante, D. A two-step, fourth-order method with energy preserving properties. Comput. Phys. Commun. 2012, 183, 1860–1868; http://doi.org/10.1016/j.cpc.2012.04.002
  * [27] Brugnano, L.; Iavernaro, F.; Trigiante, D. Analysis 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. 2015, 20, 650–667; https://doi.org/10.1016/j.cnsns.2014.05.030
  * [28] Brugnano, L.; Iavernaro, F.; Zhang, R. Arbitrarily high-order energy-preserving methods for simulating the gyrocenter dynamics of charged particles. J. Comput. Appl. Math. 2020, 380, 112994; http://doi.org/
  * [29] Brugnano, L.; Montijano, J.I.; Rández, L. On the effectiveness of spectral methods for the numerical solution of multi-frequency highly-oscillatory Hamiltonian problems. Numer. Algorithms 2019, 81, 345–376; https://doi.org/10.1007/s11075-018-0552-9
  * [30] Brugnano, L.; Montijano, J.I.; Rández, L.; High-order energy-conserving Line Integral Methods for charged particle dynamics. J. Comput. Phys. 2019, 396, 209–227; https://doi.org/10.1016/j.jcp.2019.06.068
  * [31] Brugnano, L.; Sun, Y. Multiple invariants conserving Runge-Kutta type methods for Hamiltonian problems. Numer. Algorithms 2014, 65, 611–632; https://doi.org/10.1007/s11075-013-9769-9
  * [32] Brugnano, L.; Gurioli, G.; Sun, Y. Energy-conserving Hamiltonian Boundary Value Methods for the numerical solution of the Korteweg-de Vries equation. J. Comput. Appl. Math. 2019, 351, 117–135; http://doi.org/
  * [33] Brugnano, L.; Zhang, C.; Li, D. A class of energy-conserving Hamiltonian boundary value methods for nonlinear Schrödinger equation with wave operator. Commun. Nonlinear Sci. Numer. Simul. 2018, 60, 33–49; https://doi.org/10.1016/j.cnsns.2017.12.018
  * [34] Butcher, J.C. An algebraic theory of integration methods. Math. Comp. 1972, 26, 79–106; https://doi.org/10.2307/2004720
  * [35] Butcher, J.C. The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods, John Wiley & Sons, Chichester, 1987.
  * [36] Butcher, J.C.; Wanner, G. Runge-Kutta methods: some historical notes. Appl. Numer. Math. 1996, 22, 113–151; https://doi.org/10.1016/S0168-9274(96)00048-7
  * [37] Dahlquist, G.; Björk, Å. Numerical Methods in Scientific Computing; SIAM: Philadelphia, PA, USA, 2008; Volume I.
  * [38] Hairer, E. Energy-preserving variant of collocation methods. JNAIAM J. Numer. Anal. Ind. Appl. Math. 2010, 5, 73–84.
  * [39] Hairer, E.; Wanner, G. Solving Ordinary Differential Equations II, 2nd revised edition. Springer, Heidelberg, 2002.
  * [40] Li, J.; Wu, X. Energy-preserving continuous stage extended Runge-Kutta-Nyström methods for oscillatory Hamiltonian systems. Appl. Numer. Math. 2019, 145, 469–487; https://doi.org/10.1016/j.apnum.2019.05.009
  * [41] Miyatake, Y. ; Butcher, J.C. A characterization of energy preserving methods and the construction of parallel integrators for Hamiltonian systems. SIAM J. Numer. Anal. 2016, 54, No. 3, 1993–2013; https://doi.org/10.1137/15M1020861
  * [42] Quispel, G.R.W. ; McLaren, D.I. A new class of energy-preserving numerical integration methods. J. Phys. A Math. Theor. 2008, 41, 045206; https://doi.org/10.1088/1751-8113/41/4/045206
  * [43] Tang, Q.; Chen, C.M. Continuous finite element methods for Hamiltonian systems. Appl. Math. Mech. 2007, 28, 1071–1080; https://doi.org/10.1007/s10483-007-0809-y
  * [44] Tang, W.; A note on continuous-stage Runge-Kutta methods. Appl. Math. Comput. 2018, 339, 231–241; https://doi.org/10.1016/j.amc.2018.07.044
  * [45] Tang, W.; Sun, Y. Time finite element methods: A unified framework for numerical discretizations of ODEs. Appl. Math. Comput. 2012, 219, 2158–2179; https://doi.org/10.1016/j.amc.2012.08.062
  * [46] Tang, W.; Sun, Y. Construction of Runge-Kutta type methods for solving ordinary differential equations. Appl. Math. Comput. 2014, 234, 179–191; https://doi.org/10.1016/j.amc.2014.02.042
  * [47] Tang, W.; Sun, Y.; Zhang, J. High order symplectic integrators based on continuous-stage Runge-Kutta-Nyström methods. Appl. Math. Comput. 2019, 361, 670–679; https://doi.org/10.1016/j.amc.2019.06.031
  * [48] Tang, W.; Zhang, J. Symplecticity-preserving continuous stage Runge-Kutta-Nyström methods. Appl. Math. Comput. 2018, 323, 204–219; https://doi.org/10.1016/j.amc.2017.11.054