Continuous-Stage Runge–Kutta Approximation to Differential Problems

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. In this review paper, we recall this aspect and extend it to higher-order differential problems.

1. Introduction

Continuous-stage Runge–Kutta methods (csRK methods, hereafter), introduced by Butcher [1,2,3], have been used in recent years as a useful tool for studying structure-preserving methods for Hamiltonian problems (see, e.g., [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]). In particular, we shall at first consider methods, within this latter class, having Butcher tableau in the form:

$$\begin{array}{cc}c& a(c,\tau )\\ & b\left(c\right)\end{array},$$

(1)

where

$$a:[0,1]\times [0,1]\to \mathbb{R},b:[0,1]\to [0,\infty ),$$

(2)

are suitable functions defining the method. Hereafter, we shall use the notation (1), 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

$$\dot{y}\left(t\right)=f\left(y\left(t\right)\right),t\in [0,h],y\left(0\right)=y_0\in {\mathbb{R}}^m,$$

(4)

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

$$\dot{\sigma }\left(ch\right)={\int }_0^1\dot{a}(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1],\sigma \left(0\right)=y_0,$$

(5)

to (4), and therefore,

$$\sigma \left(ch\right)=y_0+h{\int }_0^{1}a(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1],$$

(6)

with the approximation to $y\left(h\right)$ given by

$$y_1=y_0+h{\int }_0^{1}b\left(c\right)f\left(\sigma \left(ch\right)\right)\mathrm{d}c.$$

(7)

In particular, we shall hereafter consider the natural choice

$$b\left(c\right)\equiv 1,c\in [0,1],$$

(8)

though different choices have been also considered [9]. In such a case, because of consistency, one requires that

$${\int }_0^{1}b\left(c\right)\mathrm{d}c=1.$$

In this review paper, we consider the csRKs derived from the energy-conserving methods called Hamiltonian Boundary Value Methods (HBVMs), which have been the subject of many investigations in recent years (see, e.g., [4,5,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37]). In particular, the used arguments strictly follow those in [21] (in turn, inspired by [12,38,39]). However, we also consider an interesting generalization, with respect to the arguments in [21], 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, thereby obtaining a continuous-stage Runge–Kutta-Nyström method (csRKN method, hereafter), whose generalization for kth-order problems is also sketched out; Section 4 is devoted to derive relevant families of methods for the previous classes of problems, and we report some numerical testing; 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 [21], for our analysis we shall use an expansion of the vector field along the orthonormal Legendre polynomial basis:

$$P_i\in {\Pi }_i,{\int }_0^1{P}_i\left(x\right)P_j\left(x\right)\mathrm{d}x={\delta }_{ij},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, assuming for sake of simplicity that f is analytical, we can rewrite (4) as

$$\dot{y}\left(ch\right)=\sum _{j\ge 0}{P}_j\left(c\right){\gamma }_j\left(y\right),c\in [0,1],y\left(0\right)=y_0\in {\mathbb{R}}^m,$$

(10)

with

$${\gamma }_j\left(y\right)={\int }_0^1{P}_j\left(\tau \right)f\left(y\left(\tau h\right)\right)\mathrm{d}\tau ,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\left(ch\right)=y_0+h\sum _{j\ge 0}{\int }_0^c{P}_j\left(x\right)\mathrm{d}x{\gamma }_j\left(y\right),c\in [0,1].$$

(12)

For $c=1$, by considering that, by virtue of (9), ${\int }_0^1{P}_j\left(x\right)\mathrm{d}x={\delta }_{j0}$, and taking into account (4), one obtains:

$$y\left(h\right)=y_0+h{\int }_0^{1}f\left(y\left(ch\right)\right)\mathrm{d}c\equiv y_0+{\int }_0^h\dot{y}\left(t\right)\mathrm{d}t,$$

(13)

i.e., the fundamental theorem of the calculus. Interestingly, by setting

$$a_{\infty }(c,\tau )=\sum _{j\ge 0}{\int }_0^c{P}_j\left(x\right)\mathrm{d}x{P}_j\left(\tau \right),$$

(14)

one obtains that (10)–(13) can be rewritten (see (3)), respectively, as:

$$\begin{array}{ccc}\hfill \dot{y}\left(ch\right)& =& {\int }_0^1{\dot{a}}_{\infty }(c,\tau )f\left(y\left(\tau h\right)\right)\mathrm{d}\tau ,y\left(0\right)=y_0\in {\mathbb{R}}^m,\hfill \\ \hfill y\left(ch\right)& =& y_0+h{\int }_0^1{a}_{\infty }(c,\tau )f\left(y\left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1],\hfill \\ \hfill y\left(h\right)& =& y_0+h{\int }_0^{1}f\left(y\left(ch\right)\right)\mathrm{d}c.\hfill \end{array}$$

(15)

Consequently,

$$\begin{array}{cc}c& a_{\infty }(c,\tau )\\ & 1\end{array}$$

(16)

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

2.1. Vector Formulation

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

$${\mathcal{P}}_{\infty }\left(c\right)=\left(\begin{array}{c}{P}_0\left(c\right)\\ P_1\left(c\right)\\ ⋮\end{array}\right),{\mathcal{I}}_{\infty }\left(c\right)={\int }_0^c{\mathcal{P}}_{\infty }\left(x\right)\mathrm{d}x,$$

(17)

and the infinite matrix,

$$X_{\infty }=\left(\begin{array}{cc}{\xi }_0& -{\xi }_1\\ {\xi }_1& 0& -{\xi }_2\\ & {\xi }_2& \ddots & \ddots \\ & & \ddots \end{array}\right),{\xi }_i=\frac{1}{2\sqrt{|4i^2-1|}},i=0,1,\dots ,$$

(18)

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

$${\mathcal{I}}_{\infty }{\left(c\right)}^{\top }=P_{\infty }{\left(c\right)}^{\top }{X}_{\infty },{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right){\mathcal{P}}_{\infty }{\left(\tau \right)}^{\top }\mathrm{d}\tau =I,{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right){\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }\mathrm{d}\tau =X_{\infty },$$

(19)

with I being the identity operator. Consequently, (14) can be rewritten as

$$a_{\infty }(c,\tau )={\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{\mathcal{P}}_{\infty }\left(\tau \right)={\mathcal{P}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)$$

(20)

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

$$\begin{array}{cc}c& {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)\\ & 1\end{array},$$

(21)

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

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

$$\mathbf{\gamma }:=\left(\begin{array}{c}{\gamma }_0\left(y\right)\\ {\gamma }_1\left(y\right)\\ ⋮\end{array}\right)\equiv {\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f\left(y\left(\tau h\right)\right)\mathrm{d}\tau$$

(22)

where, as is usual, ⊗ denotes the Kronecker product, then (12) can be rewritten as

$$y\left(ch\right)=y_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m\mathbf{\gamma }$$

and, consequently, the vector $\mathbf{\gamma }$ satisfies:

$$\mathbf{\gamma }={\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f(y_0+h{\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }\otimes I_m\mathbf{\gamma })\mathrm{d}\tau ,$$

(23)

with (compare with (15))

$$y\left(h\right)=y_0+h{\gamma }_0\left(y\right).$$

(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\left(x\right)\mathrm{d}x{P}_j\left(\tau \right),$$

(25)

so that

$${\dot{a}}_s(c,\tau )=\sum _{j=0}^{s-1}{P}_j\left(c\right)P_j\left(\tau \right),$$

(26)

and therefore, (15) is approximated by

$$\begin{array}{ccc}\hfill \dot{\sigma }\left(ch\right)& =& {\int }_0^1{\dot{a}}_s(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,\sigma \left(0\right)=y_0\in {\mathbb{R}}^m,\hfill \\ \hfill \sigma \left(ch\right)& =& y_0+h{\int }_0^1{a}_s(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1],\hfill \\ \hfill y_1:=\sigma \left(h\right)& =& y_0+h{\int }_0^{1}f\left(\sigma \left(ch\right)\right)\mathrm{d}c.\hfill \end{array}$$

(27)

Remark 1.

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

$$\begin{array}{cc}c& a_s(c,\tau )\\ & 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 [4] (in particular, when $s=1$ one obtains the AVF method [10]; see also [6] for a related approach).

The following property holds true ([5] Theorem 1):

$$y_1-y\left(h\right)=O\left(h^{2s+1}\right);$$

(29)

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

Remark 2.

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

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

then  $H\left(y_1\right)=H\left(y_0\right)$. 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 [5] 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

$${\mathcal{P}}_r\left(c\right)=\left(\begin{array}{c}{P}_0\left(c\right)\\ ⋮\\ P_{r-1}\left(c\right)\end{array}\right),r=s,s+1,{\mathcal{I}}_s\left(c\right)={\int }_0^c{\mathcal{P}}_s\left(x\right)\mathrm{d}x,$$

(30)

and the matrices

$${\widehat{X}}_s=\left(\begin{array}{cc}{\xi }_0& -{\xi }_1\\ {\xi }_1& 0& \ddots \\ & \ddots & \ddots & -{\xi }_{s-1}\\ & & {\xi }_{s-1}& 0\\ & & & {\xi }_s\end{array}\right)\equiv \left(\begin{array}{c}{X}_s\\ 0\dots 0{\xi }_s\end{array}\right)\in {\mathbb{R}}^{s+1\times s},$$

(31)

such that

$${\mathcal{I}}_s{\left(c\right)}^{\top }=P_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s,{\int }_0^1{\mathcal{P}}_r\left(\tau \right){\mathcal{P}}_r{\left(\tau \right)}^{\top }\mathrm{d}\tau =I_r,{\int }_0^1{\mathcal{P}}_s\left(\tau \right){\mathcal{I}}_s{\left(\tau \right)}^{\top }\mathrm{d}\tau =X_s,$$

(32)

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

$$a_s(c,\tau )={\mathcal{I}}_s{\left(c\right)}^{\top }{\mathcal{P}}_s\left(\tau \right)={\mathcal{P}}_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s{\mathcal{P}}_s\left(\tau \right),$$

(33)

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

$$\begin{array}{cc}c& {\mathcal{P}}_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s{\mathcal{P}}_s\left(\tau \right)\\ & 1\end{array},$$

(34)

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

Further, according to [32], setting the vector

$${\mathbf{\gamma }}_s:=\left(\begin{array}{c}{\gamma }_0\left(\sigma \right)\\ ⋮\\ {\gamma }_{s-1}\left(\sigma \right)\end{array}\right),$$

(35)

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

$${\mathbf{\gamma }}_s={\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_{m}f(y_0+h{\mathcal{I}}_s{\left(\tau \right)}^{\top }\otimes I_m{\mathbf{\gamma }}_s)\mathrm{d}\tau ,$$

(36)

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

$$y_1=y_0+h{\gamma }_0\left(\sigma \right).$$

(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}\left(t\right)=f\left(y\left(t\right)\right),t\in [0,h],y\left(0\right)=y_0,\dot{y}\left(0\right)={\dot{y}}_0\in {\mathbb{R}}^m,$$

(38)

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

$$\ddot{y}\left(t\right)=f(y\left(t\right),\dot{y}\left(t\right)),t\in [0,h],y\left(0\right)=y_0,\dot{y}\left(0\right)={\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 we discuss (39). Setting $\dot{y}\left(t\right)=p\left(t\right)$, and expanding the right-hand sides along the Legendre basis, gives:

$$\begin{array}{ccc}\hfill \dot{y}\left(ch\right)& =& {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}p\left(\tau h\right)\mathrm{d}\tau ,\hfill \\ \hfill \dot{p}\left(ch\right)& =& {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f\left(y\left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1].\hfill \end{array}$$

(40)

Remark 3.

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

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

(41)

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

$$\begin{array}{ccc}\hfill y\left(ch\right)& =& y_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}p\left(\tau h\right)\mathrm{d}\tau ,\hfill \\ \hfill p\left(ch\right)& =& {\dot{y}}_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f\left(y\left(\tau h\right)\right)\mathrm{d}\tau \hfill \\ & =& {\dot{y}}_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m\mathbf{\gamma },c\in [0,1],\hfill \end{array}$$

(42)

with the vector $\mathbf{\gamma }$ formally still given by (22). Setting, in general, $e_i$ the infinite vector whose jth entry is ${\delta }_{ij}$, substitution of the second equation in (42) into the first one, taking into account that

$${\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\mathrm{d}\tau =e_1\equiv \left(\begin{array}{c}1\\ 0\\ ⋮\end{array}\right),{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{e}_1=c,$$

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

$$\begin{array}{ccc}\hfill y\left(ch\right)& =& y_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_m\left[{\dot{y}}_0+h{\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }\otimes I_m\mathbf{\gamma }\right]\mathrm{d}\tau \hfill \\ \hfill & =& y_0+ch{\dot{y}}_0+h^2{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right){\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }\mathrm{d}\tau \otimes I_m\mathbf{\gamma }\hfill \\ \hfill & =& y_0+ch{\dot{y}}_0+h^2{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }\otimes I_m\mathbf{\gamma }\hfill \\ \hfill & \equiv & y_0+ch{\dot{y}}_0+h^2{\int }_0^1\left[{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)\right]\otimes I_{m}f\left(y\left(\tau h\right)\right)\mathrm{d}\tau .\hfill \end{array}$$

(43)

By setting

$${\overline{a}}_{\infty }(c,\tau )={\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)\equiv {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }^2{\mathcal{P}}_{\infty }\left(\tau \right),$$

(44)

one has that (42) can be rewritten as:

$$y\left(ch\right)=y_0+ch{\dot{y}}_0+h^2{\int }_0^1{\overline{a}}_{\infty }(c,\tau )f\left(y\left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1].$$

(45)

Remark 4.

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

$$p\left(ch\right)={\dot{y}}_0+h{\int }_0^1{a}_{\infty }(c,\tau )f\left(y\left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1].$$

(46)

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

$${\overline{a}}_{\infty }(c,\tau )={\int }_0^1{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{\mathcal{P}}_{\infty }\left(\xi \right){\mathcal{I}}_{\infty }{\left(\xi \right)}^{\top }{\mathcal{P}}_{\infty }\left(\tau \right)\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:

$$\dot{y}\left(h\right)\equiv p\left(h\right)={\dot{y}}_0+h{e}_1^{\top }\otimes I_m\mathbf{\gamma }={\dot{y}}_0+h{\gamma }_0\left(y\right)\equiv {\dot{y}}_0+h{\int }_0^{1}f\left(y\left(ch\right)\right)\mathrm{d}c,$$

(48)

and (see (18))

$$\begin{array}{ccc}\hfill y\left(h\right)& =& y_0+h{\dot{y}}_0+h^2{e}_1^{\top }{X}_{\infty }{\int }_0^1{\mathcal{P}}_{\infty }\left(c\right)\otimes I_{m}f\left(y\left(ch\right)\right)\mathrm{d}c\hfill \\ \hfill & =& y_0+h{\dot{y}}_0+h^2{e}_1^{\top }{X}_{\infty }{\int }_0^1{\mathcal{P}}_{\infty }\left(c\right)\otimes I_{m}f\left(y\left(ch\right)\right)\mathrm{d}c\hfill \\ \hfill & =& y_0+h{\dot{y}}_0+h^2{\int }_0^1{({\xi }_0{e}_1-{\xi }_1{e}_2)}^{\top }{\mathcal{P}}_{\infty }\left(c\right)\otimes I_{m}f\left(y\left(ch\right)\right)\mathrm{d}c\hfill \\ \hfill & =& y_0+h{\dot{y}}_0+h^2{\int }_0^1[{\xi }_0{P}_0\left(c\right)-{\xi }_1{P}_1\left(c\right)]f\left(y\left(ch\right)\right)\mathrm{d}c\hfill \\ \hfill & \equiv & y_0+h{\dot{y}}_0+h^2{\int }_0^1(1-c)f\left(y\left(ch\right)\right)\mathrm{d}c.\hfill \end{array}$$

(49)

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

$$\begin{array}{cc}c& {\overline{a}}_{\infty }(c,\tau )\\ & 1-c\\ & 1\end{array}\equiv \begin{array}{cc}c& {\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)\\ & 1-c\\ & 1\end{array}\equiv \begin{array}{cc}c& {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }^2{\mathcal{P}}_{\infty }\left(\tau \right)\\ & 1-c\\ & 1\end{array}.$$

(50)

3.1. Vector Formulation

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

$$\mathbf{\gamma }={\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f\left(y_0+\tau h{\dot{y}}_0+h^2{\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }{X}_{\infty }\otimes I_m\mathbf{\gamma }\right)\mathrm{d}\tau ,$$

(51)

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

$$\dot{y}\left(h\right)={\dot{y}}_0+h{\gamma }_0\left(y\right),y\left(h\right)=y_0+h{\dot{y}}_0+h^2\left({\xi }_0{\gamma }_0\left(y\right)-{\xi }_1{\gamma }_1\left(y\right)\right).$$

(52)

3.2. The Case of the General Second-Order Problem

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) and (46) now become:

$$\begin{array}{ccc}\hfill y\left(ch\right)& =& y_0+ch{\dot{y}}_0+h^2{\int }_0^1{\overline{a}}_{\infty }(c,\tau )f(y\left(\tau h\right),\dot{y}\left(\tau h\right))\mathrm{d}\tau ,\hfill \\ \hfill \dot{y}\left(ch\right)& =& {\dot{y}}_0+h{\int }_0^1{a}_{\infty }(c,\tau )f(y\left(\tau h\right),\dot{y}\left(\tau h\right))\mathrm{d}\tau ,c\in [0,1],\hfill \end{array}$$

(53)

with the values at h given by

$${\dot{y}}_0+h{\int }_0^{1}f(y\left(ch\right),\dot{y}\left(ch\right)))\mathrm{d}c,y\left(h\right)=y_0+h{\dot{y}}_0+h^2{\int }_0^1(1-c)f(y\left(ch\right),\dot{y}\left(ch\right))\mathrm{d}c.$$

(54)

Consequently, we obtain the general csRKN method

$$\begin{array}{ccc}\hfill \begin{array}{ccc}c& {\overline{a}}_{\infty }(c,\tau )& a_{\infty }(c,\tau )\\ & 1-c& 1\end{array}& \equiv & \begin{array}{ccc}c& {\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)& {\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{\mathcal{P}}_{\infty }\left(\tau \right)\\ & 1-c& 1\end{array}\hfill \\ \hfill & \equiv & \begin{array}{ccc}c& {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }^2{\mathcal{P}}_{\infty }\left(\tau \right)& {\mathcal{P}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }{\mathcal{P}}_{\infty }\left(\tau \right)\\ & 1-c& 1\end{array},\hfill \end{array}$$

(55)

in place of (50).

The bad news is that now the system of Equations (53) has a doubled size, with respect to (45). On the other hand, the good news is that, upon modifying the definition of the vector $\mathbf{\gamma }$ in (22) as follows,

$$\mathbf{\gamma }:=\left(\begin{array}{c}{\gamma }_0(y,\dot{y})\\ {\gamma }_1(y,\dot{y})\\ ⋮\end{array}\right)\equiv {\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f(y\left(\tau h\right),\dot{y}\left(\tau h\right))\mathrm{d}\tau ,$$

(56)

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

$$\begin{array}{ccc}\hfill y\left(ch\right)& =& y_0+ch{\dot{y}}_0+h^2{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }{X}_{\infty }\otimes I_m\mathbf{\gamma },\hfill \\ \hfill \dot{y}\left(ch\right)& =& {\dot{y}}_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m\mathbf{\gamma },c\in [0,1].\hfill \end{array}$$

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

$$\mathbf{\gamma }={\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f\left(y_0+\tau h{\dot{y}}_0+h^2{\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }{X}_{\infty }\otimes I_m\mathbf{\gamma },{\dot{y}}_0+h{\mathcal{I}}_{\infty }{\left(c\right)}^{\top }\otimes I_m\mathbf{\gamma }\right)\mathrm{d}\tau .$$

(57)

Similarly, the new approximations in (54) become:

$$\dot{y}\left(h\right)={\dot{y}}_0+h{\gamma }_0(y,\dot{y}),y\left(h\right)=y_0+h{\dot{y}}_0+h^2\left({\xi }_0{\gamma }_0(y,\dot{y})-{\xi }_1{\gamma }_1(y,\dot{y})\right).$$

(58)

Remarkably enough, the Equations (57) and (58) are very similar to (51) and (52), respectively.

3.3. Polynomial Approximation

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

$$\begin{array}{ccc}\hfill \dot{\sigma }\left(ch\right)& =& {\mathcal{P}}_s{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_m{\sigma }_1\left(\tau h\right)\mathrm{d}\tau ,\hfill \\ \hfill {\dot{\sigma }}_1\left(ch\right)& =& {\mathcal{P}}_s{\left(c\right)}^{\top }\otimes I_m{\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_{m}f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1],\hfill \end{array}$$

(59)

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

$${\overline{a}}_s(c,\tau )={\mathcal{I}}_s{\left(c\right)}^{\top }{X}_s{\mathcal{P}}_s\left(\tau \right)\equiv {\mathcal{P}}_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s{X}_s{\mathcal{P}}_s\left(\tau \right),$$

(60)

from (59) one derives:

$$\sigma \left(ch\right)=y_0+ch{\dot{y}}_0+h^2{\int }_0^1{\overline{a}}_s(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1],$$

(61)

with the new approximations,

$$y_1:=\sigma \left(h\right)and{\dot{y}}_1:={\sigma }_1\left(h\right),$$

given by

$${\dot{y}}_1={\dot{y}}_0+h{\int }_0^{1}f\left(\sigma \left(ch\right)\right)\mathrm{d}c,y_1=y_0+h{\dot{y}}_0+h^2{\int }_0^1(1-c)f\left(\sigma \left(ch\right)\right)\mathrm{d}c.$$

(62)

Remark 5.

Similarly to (41), one derives:

$$\dot{\sigma }\left(ch\right)={\int }_0^1{\dot{a}}_s(c,\tau ){\sigma }_1\left(\tau h\right)\mathrm{d}\tau ,{\dot{\sigma }}_1\left(ch\right)={\int }_0^1{\dot{a}}_s(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1].$$

(63)

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

$${\sigma }_1\left(ch\right)={\dot{y}}_0+h{\int }_0^1{a}_s(c,\tau )f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,c\in [0,1].$$

(64)

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

$${\overline{a}}_s(c,\tau )={\int }_0^1{\mathcal{I}}_s{\left(c\right)}^{\top }{\mathcal{P}}_s\left(\xi \right){\mathcal{I}}_s{\left(\xi \right)}^{\top }{\mathcal{P}}_s\left(\tau \right)\mathrm{d}\xi \equiv {\int }_0^1{a}_s(c,\xi )a_s(\xi ,\tau )\mathrm{d}\xi .$$

(65)

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

$$\begin{array}{cc}c& {\overline{a}}_s(c,\tau )\\ & 1-c\\ & 1\end{array}\equiv \begin{array}{cc}c& {\mathcal{I}}_s{\left(c\right)}^{\top }{X}_s{\mathcal{P}}_s\left(\tau \right)\\ & 1-c\\ & 1\end{array}\equiv \begin{array}{cc}c& {\mathcal{P}}_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s{X}_s{\mathcal{P}}_s\left(\tau \right)\\ & 1-c\\ & 1\end{array}.$$

(66)

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) [21].

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

$${\mathbf{\gamma }}_s:=\left(\begin{array}{c}{\gamma }_j\left(\sigma \right)\\ ⋮\\ {\gamma }_{s-1}\left(\sigma \right)\end{array}\right)\equiv {\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_{m}f\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,$$

(67)

by virtue of (60) and (61), such a vector satisfies the equation

$${\mathbf{\gamma }}_s={\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_{m}f\left(y_0+\tau h{\dot{y}}_0+h^2{\mathcal{I}}_s{\left(\tau \right)}^{\top }{X}_s\otimes I_m{\mathbf{\gamma }}_s\right)\mathrm{d}\tau ,$$

(68)

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

$${\dot{y}}_1={\dot{y}}_0+h{\gamma }_0\left(\sigma \right),y_1=y_0+h{\dot{y}}_0+h^2\left({\xi }_0{\gamma }_0\left(\sigma \right)-{\xi }_1{\gamma }_1\left(\sigma \right)\right).$$

(69)

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

$$\begin{array}{ccc}\hfill \begin{array}{ccc}c& {\overline{a}}_s(c,\tau )& a_s(c,\tau )\\ & 1-c& 1\end{array}& \equiv & \begin{array}{ccc}c& {\mathcal{I}}_s{\left(c\right)}^{\top }{X}_s{\mathcal{P}}_s\left(\tau \right)& {\mathcal{I}}_s{\left(c\right)}^{\top }{\mathcal{P}}_s\left(\tau \right)\\ & 1-c& 1\end{array}\hfill \\ \hfill & \equiv & \begin{array}{ccc}c& {\mathcal{P}}_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s{X}_s{\mathcal{P}}_s\left(\tau \right)& {\mathcal{P}}_{s+1}{\left(c\right)}^{\top }{\widehat{X}}_s{\mathcal{P}}_s\left(\tau \right)\\ & 1-c& 1\end{array}.\hfill \end{array}$$

(70)

Additionally, we can derive a vector formulation similar to (67) and (68). In fact, by defining the vector

$${\mathbf{\gamma }}_s:=\left(\begin{array}{c}{\gamma }_j(\sigma ,{\sigma }_1)\\ ⋮\\ {\gamma }_{s-1}(\sigma ,{\sigma }_1)\end{array}\right)\equiv {\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_{m}f(\sigma \left(\tau h\right),{\sigma }_1\left(\tau h\right))\mathrm{d}\tau ,$$

(71)

it turns out that it satisfies the equation

$${\mathbf{\gamma }}_s={\int }_0^1{\mathcal{P}}_s\left(\tau \right)\otimes I_{m}f\left(y_0+\tau h{\dot{y}}_0+h^2{\mathcal{I}}_s{\left(\tau \right)}^{\top }{X}_s\otimes I_m{\mathbf{\gamma }}_s,{\dot{y}}_0+h{\mathcal{I}}_s{\left(\tau \right)}^{\top }\otimes I_m{\mathbf{\gamma }}_s\right)\mathrm{d}\tau ,$$

(72)

with the new approximations given by:

$${\dot{y}}_1={\dot{y}}_0+h{\gamma }_0(\sigma ,{\sigma }_1),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).$$

(73)

Remark 7.

By comparing the discrete problems (68) and (72), 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 out in the next section, this will be the case for a general kth order ODE-IVP.

3.4. Approximation of General kth-Order ODE-IVPs

Let us consider the case of a general kth order problem:

$$y^{\left(k\right)}\left(ch\right)=f(y\left(ch\right),y^{\left(1\right)}\left(ch\right),\dots ,y^{(k-1)}\left(ch\right)),c\in [0,1],y^{\left(i\right)}\left(0\right)=y_0^{\left(i\right)}\in {\mathbb{R}}^m,i=0,\dots ,k-1.$$

(74)

Hereafter, for sake of brevity, we shall skip the arguments of the Fourier coefficients ${\gamma }_i$. Then, by repeating similar steps as above, by defining the infinite vector

$$\mathbf{\gamma }\equiv \left(\begin{array}{c}{\gamma }_0\\ {\gamma }_1\\ ⋮\end{array}\right):={\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f(y\left(\tau h\right),\dots ,y^{(k-1)}\left(\tau h\right))\mathrm{d}\tau ,$$

(75)

one has that it satisfies the equation:

$$\begin{array}{ccc}\hfill {\displaystyle \mathbf{\gamma }=}& & {\int }_0^1{\mathcal{P}}_{\infty }\left(\tau \right)\otimes I_{m}f\left(\sum _{i=0}^{k-1}\frac{{\left(\tau h\right)}^i}{i!}{y}_0^{\left(i\right)}+h^k{\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }{X}_{\infty }^{k-1}\otimes I_m\mathbf{\gamma },\dots ,y_0^{(k-1)}+h{\mathcal{I}}_{\infty }{\left(\tau \right)}^{\top }\otimes I_m\mathbf{\gamma }\right)\mathrm{d}\tau ,\hfill \end{array}$$

(76)

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

$$y^{\left(i\right)}\left(h\right)=\sum _{j=0}^{k-1-i}\frac{h^j}{j!}{y}_0^{(j+i)}+h^{k-i}\sum _{j=0}^{k-1-i}{b}_j^{(k-1-i)}{\gamma }_j,i=0,\dots ,k-1,$$

(77)

$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},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 (76), ${\mathcal{P}}_{\infty }\left(\tau \right)$, ${\mathcal{I}}_{\infty }\left(\tau \right)$, and $X_{\infty }$, with ${\mathcal{P}}_s\left(\tau \right)$, ${\mathcal{I}}_s\left(\tau \right)$, and $X_s$, respectively, (consequently, the vector $\mathbf{\gamma }$ now belongs to ${\mathbb{R}}^{sm}$), with the new approximations $y_1^{\left(i\right)}\approx y^{\left(i\right)}\left(h\right)$, given by:

$$y_1^{\left(i\right)}:=\sum _{j=0}^{k-1-i}\frac{h^j}{j!}{y}_0^{(j+i)}+h^{k-i}\sum _{j=0}^{k-1-i}{b}_j^{(k-1-i)}{\gamma }_j,i=0,\dots ,k-1,$$

(78)

$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},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 evaluated by numerical methods.”

Dahlquist and Björk [41] p. 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\ge s$, having weights and abscissae $(b_i,c_i)$, $i=1,\dots ,k$. In so doing, the vector of the Fourier coefficients (35) becomes

$${\widehat{\mathbf{\gamma }}}_s:=\left(\begin{array}{c}{\widehat{\gamma }}_0\\ ⋮\\ {\widehat{\gamma }}_{s-1}\end{array}\right),$$

(79)

satisfying the Equation (compare with (36))

$${\widehat{\mathbf{\gamma }}}_s=\sum _{j=1}^k{b}_j{\mathcal{P}}_s\left(c_j\right)\otimes I_{m}f(y_0+h{\mathcal{I}}_s\left(c_j\right)\otimes I_m{\widehat{\mathbf{\gamma }}}_s)$$

(80)

with (compare with (37)),

$$y_1=y_0+h{\widehat{\gamma }}_0.$$

(81)

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

$${\mathcal{P}}_r:=\left(\begin{array}{c}{\mathcal{P}}_r{\left(c_1\right)}^{\top }\\ ⋮\\ {\mathcal{P}}_r{\left(c_k\right)}^{\top }\end{array}\right)\in {\mathbb{R}}^{k\times r},{\mathcal{I}}_s:=\left(\begin{array}{c}{\mathcal{I}}_s{\left(c_1\right)}^{\top }\\ ⋮\\ {\mathcal{I}}_s{\left(c_k\right)}^{\top }\end{array}\right)\in {\mathbb{R}}^{k\times s},\Omega =\left(\begin{array}{c}{b}_1\\ & \ddots \\ & & b_k\end{array}\right),$$

(82)

one has, with reference to (31),

$${\mathcal{I}}_s={\mathcal{P}}_{s+1}{\widehat{X}}_s,{\mathcal{P}}_s^{\top }\Omega {\mathcal{P}}_s=I_s,{\mathcal{P}}_s^{\top }\Omega {\mathcal{I}}_s=X_s.$$

(83)

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

$${\widehat{\mathbf{\gamma }}}_s={\mathcal{P}}_s^{\top }\Omega \otimes I_{m}f(\mathit{e}\otimes y_0+h{\mathcal{I}}_s\otimes I_m{\widehat{\mathbf{\gamma }}}_s).$$

(84)

By setting

$$Y\equiv \left(\begin{array}{c}{Y}_1\\ ⋮\\ Y_k\end{array}\right):=\mathit{e}\otimes y_0+h{\mathcal{I}}_s\otimes I_m{\widehat{\mathbf{\gamma }}}_s,$$

(85)

the vector of the stages of the corresponding Runge–Kutta method, from (84) and (85) one obtains the stage equation

$$Y=\mathit{e}\otimes y_0+h{\mathcal{I}}_s{\mathcal{P}}_s^{\top }\Omega \otimes I_{m}f\left(Y\right),$$

(86)

with the new approximation

$$y_1=y_0+h\sum _{i=1}^k{b}_{i}f\left(Y_i\right).$$

(87)

Summing all up, (84)–(87) define the k-stage Runge–Kutta method:

$$\begin{array}{cc}\mathit{c}& {\mathcal{I}}_s{\mathcal{P}}_s^{\top }\Omega \\ & {\mathit{b}}^{\top }\end{array},\mathit{b}={(b_1,\dots ,b_k)}^{\top },\mathit{c}={(c_1,\dots ,c_k)}^{\top },$$

(88)

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

$${\mathcal{I}}_s{\mathcal{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 [4,5,29,30]. In fact, it can be proved (see, e.g., [5] Corollary 3) that a HBVM$(k,s)$ method, with $k\ge s$:

• Has order $2s$;
• Is energy-conserving, for all polynomial Hamiltonians of degree not larger that $2k/s$;
• For general (and suitably regular Hamiltonians), the Hamiltonian error is $O\left(h^{2k+1}\right)$. In this case, however, by using a value of k large enough so that the Hamiltonian error falls within the round-off error level, the method turns out to be practically energy-conserving.

Remark 8.

As is clear, the formulation (80) and (81) is computationally more effective than (86) and (87), having the former (block) dimension s, independently of k, which is the dimension of the latter formulation. As observed in [32], 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 (74): here we sketch only those concerning the csRKN (66), providing the correct implementation for a HBVM$(k,s)$ method for the special second-order problem (38). By formally using the same approximate Fourier coefficients (79) in place of (67), one has that (68) is replaced by the following discrete counterpart,

$${\widehat{\mathbf{\gamma }}}_s=\sum _{j=1}^k{b}_j{\mathcal{P}}_s\left(c_j\right)\otimes I_{m}f(y_0+c_{j}h{\dot{y}}_0+h^2{\mathcal{I}}_s\left(c_j\right)X_s\otimes I_m{\widehat{\mathbf{\gamma }}}_s)$$

(89)

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

$${\dot{y}}_1={\dot{y}}_0+h{\widehat{\gamma }}_0,y_1=y_0+h{\dot{y}}_0+h^2\left({\xi }_0{\widehat{\gamma }}_0-{\xi }_1{\widehat{\gamma }}_1\right).$$

(90)

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

$${\widehat{\mathbf{\gamma }}}_s={\mathcal{P}}_s^{\top }\Omega \otimes I_{m}f(\mathit{e}\otimes y_0+h\mathit{c}\otimes {\dot{y}}_0+h^2{\mathcal{I}}_s{X}_s\otimes I_m{\widehat{\mathbf{\gamma }}}_s),$$

(91)

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

$$Y\equiv \left(\begin{array}{c}{Y}_1\\ ⋮\\ Y_k\end{array}\right):=\mathit{e}\otimes y_0+h\mathit{c}\otimes {\dot{y}}_0+h^2{\mathcal{I}}_s{X}_s\otimes I_m{\widehat{\mathbf{\gamma }}}_s$$

(92)

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

$$Y=\mathit{e}\otimes y_0+h\mathit{c}\otimes {\dot{y}}_0+h^2{\mathcal{I}}_s{X}_s{\mathcal{P}}_s^{\top }\Omega \otimes I_{m}f\left(Y\right),$$

(93)

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\left(Y_i\right),y_1=y_0+h{\dot{y}}_0+h^2\sum _{i=1}^k{b}_i(1-c_i)f\left(Y_i\right).$$

(94)

Consequently, we are speaking about the following k-stages RKN method:

$$\begin{array}{cc}\mathit{c}& {\mathcal{I}}_s{X}_s{\mathcal{P}}_s^{\top }\Omega \\ & {[\mathit{b}\circ (\mathit{e}-\mathit{c})]}^{\top }\\ & {\mathit{b}}^{\top }\end{array},$$

(95)

where  ∘  denotes the Hadamard, i.e., componentwise, product. The Butcher tableau (95) defines the RKN formulation of a HBVM$(k,s)$ method [21,29,32]. It is worth mentioning that, with reference to (60) and (65), the Butcher matrix can be written in terms of the corresponding csRKN method (66) as follows:

$${\mathcal{I}}_s{X}_s{\mathcal{P}}_s^{\top }\Omega \equiv \left(\begin{array}{c}{b}_j{\overline{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\ge s$, the formulation (90) and (91) of the method is much more efficient than the usual one given by (93) and (94), in view of the use of values of $k\gg s$, as in the case of separable Hamiltonian problems, as the numerical example below will show.

It must be emphasized that:

• Very efficient iterative procedures exist for solving the discrete problems (84) and (91) generated by a HBVM method (see, e.g., [32]). As matter of fact, we mention that the state-of-art Matlab code hbvm is available at the website of the book [29];
• When HBVMs are used as spectral methods in time, i.e., choosing values of s and $k>s$ so that no further accuracy improvement can be obtained, for the considered finite-precision arithmetic and timestep h used [22,31,35], then there is no practical difference between the discrete methods and their continuous-stage counterparts.

Numerical Tests

We now report a few numerical examples comparing energy-conserving methods with respect to symplectic methods. For this purpose, we used the same Matlab code hbvm mentioned above: in fact, for $k=s$, it also provides a very efficient implementation of symplectic s-stage Gauss collocation methods. All tests have been executed on a 3 GHz Intel Xeon W10 core computer with 64GB of memory, using Matlab^© R2020b.

Symplectic methods usually perform very well when solving Hamiltonian problems (see, e.g., the monographs [42,43,44,45] or [46,47,48,49]), but in some circumstances (e.g., when relatively large timesteps are used) energy-conserving methods may be more reliable. As an example of this fact, we consider the separable Hamiltonian problem described by the following polynomial Hamiltonian:

$$H(q,p)=\frac{1}{2}\left(p_1^2+p_2^2\right)+\frac{1}{2}\left(\beta q_1^2+q_2^2\right)+\alpha {(q_1-\gamma q_2)}^{2n},q,p\in {\mathbb{R}}^2.$$

(96)

In particular, we consider the following values of the parameters:

$$\alpha =\beta =n=5,\gamma =2.48.$$

(97)

At first, we numerically integrated the trajectory starting at $q\left(0\right)={(1,1)}^{\top }$, $p\left(0\right)={(0,0)}^{\top }$, with a timestep $h=5·{10}^{-3}$, on the interval $[0,250]$, by using the symplectic 2-stage Gauss method (i.e., HBVM(2,2)), and the HBVM(10,2) method, which is energy-conserving: in both cases, the methods were used as RKN methods. A reference solution has been computed by using a spectral HBVM [22,31,35]. The two methods required comparable execution times (8.5 s vs. 12.8 s, respectively), and the obtained numerical results are comparable as well, as is shown in Figure 1 where we have plotted: the Hamiltonian error; the phase–space plot in the $q_1-p_1$ plane; the phase–space plot in the $q_2-p_2$ plane. The Gauss method, though not energy-conserving, exhibits a quasi-conservation of the numerical Hamiltonian. Subsequently, we integrated the same orbit with a timestep $h={10}^{-2}$. Again, the two methods required comparable execution times (6.5 s vs. 8.7 s, respectively) but the obtained numerical results are quite different, as is shown in Figure 2, where we plot, as before: the Hamiltonian error; the phase–space plot in the $q_1-p_1$ plane; the phase–space plot in the $q_2-p_2$ plane. As one may see, now the energy-conserving method is qualitatively more accurate than the symplectic 2-stage Gauss method, which does not exhibit anymore the quasi-conservation of the energy.

5. Conclusions

In this paper, we have reviewed the basic facts and reported some new insights on the energy-conserving class of Runge–Kutta methods named HBVMs. Methods have been here studied within the framework of continuous-stage Runge–Kutta methods. The extension to second-order problems has been also recalled, providing a natural continuous-stage Runge–Kutta-Nyström formulation of the methods. Further, also the extension to general kth-order problems has been sketched out. The relation with the fully discrete methods has been also recalled, thereby showing the usefulness of using such a framework to study the fully discrete methods.