A note on the continuous-stage Runge-Kutta-(Nystr\"om) formulation of Hamiltonian Boundary Value Methods (HBVMs)

TL;DR

This paper explores the continuous-stage Runge-Kutta(-Nyström) formulation of Hamiltonian Boundary Value Methods (HBVMs), providing deeper insight into their structure and properties for energy-conserving numerical solutions of Hamiltonian problems.

Contribution

It introduces a continuous-stage formulation of HBVMs, enhancing understanding of their structure and potential advantages in solving Hamiltonian systems.

Findings

Provides a new perspective on HBVMs as continuous-stage methods

Deepens theoretical understanding of energy-conserving numerical methods

Facilitates potential improvements in numerical algorithms for Hamiltonian problems

Abstract

In recent years, the class of energy-conserving methods named Hamiltonian Boundary Value Methods (HBVMs) has been devised for numerically solving Hamiltonian problems. In this short note, we study their natural formulation as continuous-stage Runge-Kutta(-Nystr\"om) methods, which allows a deeper insight in the methods.

Equations

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

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

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

$$\begin{array}[]{c|c}c&a_{c\tau}\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}~{}.$$

$$\ddot{q}(t)=f(q(t)),t\in [0,h],q(0)=q_0,\dot{q}(0)=p_0\in {\mathbb{R}}^m,$$

$$P_j\in {\Pi }_j,{\int }_0^1{P}_i(c)P_j(c)\mathrm{d}c={\delta }_{ij},\forall i,j\ge 0,$$

$$\dot{y}(ch)=\sum _{j\ge 0}{P}_j(c){\gamma }_j(y),c\in [0,1],{\gamma }_j(y)={\int }_0^1{P}_j(\tau )f(y(\tau h))\mathrm{d}\tau ,j\ge 0,$$

$$y(ch)=y_0+h\sum _{j\ge 0}{\int }_0^c{P}_j(x)\mathrm{d}x{\gamma }_j(y),c\in [0,1].$$

$${\cal P}_{\infty}(c):=\left(\begin{array}[]{c}P_{0}(c)\\[2.84526pt] P_{1}(c)\\ \vdots\end{array}\right),\quad{\cal I}_{\infty}(c):=\int_{0}^{c}{\cal P}_{\infty}(x)\mathrm{d}x\equiv\left(\begin{array}[]{c}\int_{0}^{c}P_{0}(x)\mathrm{d}x\\[5.69054pt] \int_{0}^{c}P_{1}(x)\mathrm{d}x\\ \vdots\end{array}\right),\quad{\bm{\gamma}}(y)=\left(\begin{array}[]{c}\gamma_{0}(y)\\[2.84526pt] \gamma_{1}(y)\\ \vdots\end{array}\right),$$

$$\dot{y}(ch)={\mathcal{P}}_{\infty }(c{)}^{\top }\otimes I_m\mathbf{\gamma }(y),c\in [0,1],\mathbf{\gamma }(y)={\int }_0^1{\mathcal{P}}_{\infty }(\tau )\otimes I_{m}f(y(\tau h))\mathrm{d}\tau ,$$

$$y(ch)=y_0+h{\mathcal{I}}_{\infty }(c{)}^{\top }\otimes I_m\mathbf{\gamma }(y),c\in [0,1].$$

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

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

$${\int }_0^1{\mathcal{P}}_{\infty }(c){\mathcal{I}}_{\infty }(c{)}^{\top }\mathrm{d}c=X_{\infty }.$$

$$y(ch)$$

$$y_1$$

$$y(ch)$$

$$y_1$$

$$\begin{array}[]{c|c}c&{\cal I}_{\infty}(c)^{\top}{\cal P}_{\infty}(\tau)\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}\qquad\equiv\qquad\begin{array}[]{c|c}c&{\cal P}_{\infty}(c)^{\top}X_{\infty}{\cal P}_{\infty}(\tau)\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}\qquad=:\qquad\begin{array}[]{c|c}c&a_{c\tau}^{(\infty)}\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}~{}.$$

$$a_{c\tau }^{(\infty )}$$

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

$$\begin{array}[]{c|c}c&{\cal I}_{s}(c)^{\top}{\cal P}_{s}(\tau)\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}\qquad=:\qquad\begin{array}[]{c|c}c&a_{c\tau}^{(s)}\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}~{},$$

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

$$\sigma (ch)=y_0+h\sum _{j=0}^{s-1}{\int }_0^c{P}_j(x)\mathrm{d}x{\int }_0^1{P}_j(\tau )f(\sigma (\tau h))\mathrm{d}\tau ,c\in [0,1],$$

$${\cal I}_{s}(c)^{\top}={\cal P}_{s+1}(c)^{\top}\left(\begin{array}[]{ccccc}\xi_{0}&-\xi_{1}\\ \xi_{1}&0&-\xi_{2}\\ &\xi_{2}&\ddots&\ddots\\ &&\ddots&\ddots&-\xi_{s-1}\\ &&&\xi_{s-1}&0\\ \hline\cr&&&&\xi_{s}\end{array}\right)=:{\cal P}_{s+1}(c)^{\top}\hat{X}_{s}\equiv{\cal P}_{s+1}(c)^{\top}\left(\begin{array}[]{c}X_{s}\\ \hline\cr 0,\dots,0,\xi_{s}\end{array}\right),$$

$$a_{c\tau }^{(s)}$$

$$\begin{array}[]{c|c}c&{\cal I}_{s}(c)^{\top}{\cal P}_{s}(\tau)\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}\qquad\equiv\qquad\begin{array}[]{c|c}c&{\cal P}_{s+1}(c)^{\top}\hat{X}_{s}{\cal P}_{s}(\tau)\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}\qquad\equiv\qquad\begin{array}[]{c|c}c&a_{c\tau}^{(s)}\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}~{},$$

$$\begin{array}[]{c|c}c&{\cal P}_{s}(c)^{\top}X_{s}{\cal P}_{s}(\tau)\\[5.69054pt] \hline\cr\\[-5.69054pt] &1\end{array}~{},$$

$${\int }_0^1{P}_j(\tau )f(\sigma (\tau h))\mathrm{d}\tau$$

$$\begin{array}[]{c|c}{\bm{c}}&{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\\[5.69054pt] \hline\cr\\[-5.69054pt] &{\bm{b}}^{\top}\end{array}\qquad\equiv\qquad\begin{array}[]{c|c}{\bm{c}}&{\cal P}_{s+1}\hat{X}_{s}{\cal P}_{s}^{\top}\Omega\\[5.69054pt] \hline\cr\\[-5.69054pt] &{\bm{b}}^{\top}\end{array}\qquad=:\qquad\begin{array}[]{c|c}{\bm{c}}&A=(a_{ij})\\[5.69054pt] \hline\cr\\[-5.69054pt] &{\bm{b}}^{\top}\end{array}\,,$$

$${\bm{b}}=\left(\begin{array}[]{ccc}b_{1},&\dots,&b_{k}\end{array}\right)^{\top},\qquad{\bm{c}}=\left(\begin{array}[]{ccc}c_{1},&\dots,&c_{k}\end{array}\right)^{\top},$$

Full text

Keywords: continuous-stage Runge-Kutta methods, Runge-Kutta-Nyström methods, Hamiltonian Boundary Value Methods, HBVMs.

MSC: 65L05, 65P10.

1 Introduction

The numerical solution of Hamiltonian problems has been recently tackled by defining energy-conserving methods, which can be regarded as continuous-stage Runge-Kutta (RK, hereafter) methods (e.g., [34, 18, 30, 36]). In their simplest (and most effective) form,111For a more general form which, however, we shall not consider here, we refer, e.g., to [33]. continuous-stage RK methods are “methods” that, when applied for solving an initial value problem for ODEs (ODE-IVP, hereafter), which we assume without loss of generality in the form

[TABLE]

with $f$ analytical, define an approximating function $u:[0,h]\rightarrow\mathbb{R}^{m}$ such that

[TABLE]

with $a_{c\tau}:[0,1]\times[0,1]\rightarrow\mathbb{R}$, and a corresponding approximation to $y(h)$,

[TABLE]

As is usual, this procedure can be summarized by the following (generalized) Butcher tableau,

[TABLE]

We observe that (2)-(3) is not yet an actual numerical method, due to the fact that the involved integrals need to be conveniently approximated by means of quadrature rules. In so doing, one obtains “usual” RK methods.222I.e., having discrete stages. Nevertheless, (2)-(3) can be useful for purposes of analysis [29, 20, 35, 36, 38, 37, 32] since, essentially, it allows to discuss all Runge-Kutta methods derived by using different quadratures for approximating the involved integrals. In particular, the papers [38, 37] have inspired the present note, where we provide the continuous-stage RK formulation of Hamiltonian Boundary Value Methods (HBVMs) [17, 16, 18, 21, 24, 12, 14], a class of energy-conserving methods for Hamiltonian problems, which have been developed along several directions [4, 13, 22, 23, 26, 8, 9], including Hamitonian BVPs [1], highly-oscillatory problems [25, 2], Hamiltonian PDEs [7, 3, 15, 6, 10, 27, 11], and also considering their efficient implementation [19, 5]. Here, we shall also consider the continuous formulation of such methods when applied for solving special second-order problems [12], i.e., problems in the form

[TABLE]

where, for the sake of brevity, we shall again assume $f$ to be analytical.

With these premises, the structure of the paper is as follows: in Section 2 we study the case of first order ODE problems; Section 3 is devoted to study the case where one solves special second-order problems; at last, a few concluding remarks are drawn in Section 4.

2 The framework

Generalizing the arguments in [20], let us consider the orthonormal Legendre polynomial basis $\{P_{j}\}_{j\geq 0}$ on the interval $[0,1]$:

[TABLE]

where $\Pi_{j}$ is the set of polynomials of degree $j$. Then, the ODE-IVP (1) can be written, by expanding the right-hand side along the Legendre basis, as

[TABLE]

from which, integrating side by side, one obtains the following formal expression for the solution of (1):

[TABLE]

The above equations can be cast in vector form by introducing the infinite vectors

[TABLE]

respectively as:

[TABLE]

and

[TABLE]

Moreover, by considering that

[TABLE]

with $I$ the identity operator and

[TABLE]

one also obtains that

[TABLE]

Setting  $y_{1}\equiv y(h),$  we can cast (10) as:

[TABLE]

which, by virtue of (11)–(13), can be also written as

[TABLE]

In other words, we are speaking about the application of the following continuous-stage RK method to problem (1) :

[TABLE]

As is clear, by virtue of (11)-(12), the coefficients of this “continuous-stage RK method”, providing the exact solution of (1), are given by

[TABLE]

2.1 Polynomial approximation

In order to obtain a polynomial approximation $\sigma\in\Pi_{s}$ to $y$, let us now introduce the truncated vectors

[TABLE]

in place of the corresponding infinite ones in (8). In so doing, we replace (16) with the continuous-stage RK method

[TABLE]

whose coefficients are now polynomials of degree $s$. Consequently, by setting now $y_{1}\equiv\sigma(h)$ the approximation to $y(h)$, one obtains:

[TABLE]

The following straightforward result holds true.

Theorem 1

The continuous-stage RK method (19)-(20) coincides with the HBVM$(\infty,s)$ method in [18].333In particular when $s=1$ one retrieves the AVF method in [34].

Proof In fact, from (18), one has that (20) is equivalent to

[TABLE]

which, according to [18, Definition 1] (see also [14]), is the Master Functional Equation defining a HBVM$(\infty,s)$ method.  $\Box{~{}}$

Furthermore, by considering that (see (18) and (12))

[TABLE]

one easily obtains that (compare with (17))

[TABLE]

As a result, from (19) and (22), one obtains that

[TABLE]

which is clearly equivalent to (20).

Remark 1

We observe that, in a sense, (24) can be regarded as a continuous extension of the $W$-transformation in [31, Section IV.5]. Moreover, by considering, in place of (24), the following Butcher tableau,

[TABLE]

one obtains the continuous extension of the low-rank symplectic methods in [28].

2.1.1 Discretization

We conclude this section by recalling that [18, 20, 12] for the polynomial $\sigma$ defined in (20)-(21), one has $\sigma(h)-y(h)=O(h^{2s+1})$.444One could obtain the result also by using the symplifying assumptions for continuous-stage RK methods [35, 36, 33]. Moreover, by approximating the integrals

[TABLE]

appearing in (21) by means of a Gauss-Legendre formula of order $2k$, one obtains a HBVM$(k,s)$ method, which retains the order $2s$ of the approximation defined by (21), for all $k\geq s$. In particular, when $k=s$, one obtains the $s$-stage Gauss-Legendre collocation method. As a result, the Butcher tableau of a HBVM$(k,s)$ method turns out to be given by

[TABLE]

with $\hat{X}_{s}$ the matrix defined in (22),

[TABLE]

the vectors containing the weights and abscissae of the quadrature, respectively,555Any quadrature is in principle allowed, provided that it is enough accurate.

[TABLE]

and

[TABLE]

In particular, from (2.1) one obtains that the entries of matrix $A$ in (25) are given by

[TABLE]

3 Second order problems

Inspired by [38, 37] (see also [12]), we now consider the case of special second order problems, i.e., ODE-IVPs in the form (4). By setting $p(t)\equiv\dot{q}(t)$, one then obtains the following equivalent system of first order ODEs,

[TABLE]

HBVMs have been considered for numerically solving this problem [19]. We can then consider the use of HBVM$(\infty,s)$, too. To begin with, by applying same steps as above, one then obtains that (29) can be formally written as

[TABLE]

Simplifying the expressions, integrating side by side, and imposing the initial conditions, then gives

[TABLE]

Substituting the second equation in the first one, and taking into account (11)-(12), then gives, setting e_{1}=\left(\begin{array}[]{ccc}1,&0,&\dots{}\end{array}\right)^{\top} and considering that ${\cal I}_{\infty}(c)e_{1}=c$,

[TABLE]

where, by considering that (see (12))

[TABLE]

and taking into account (8), we have set :

[TABLE]

Moreover, by setting $q_{1}\equiv q(h)$ and (see (29)) $\dot{q}_{1}\equiv p(h)$, one obtains

[TABLE]

and, by also considering that  $f(q(\tau h))=\sum_{j\geq 0}P_{j}(\tau)\int_{0}^{1}P_{j}(\xi)f(q(\xi h))\mathrm{d}\xi$, $\tau\in[0,1]$,

[TABLE]

Next, by taking into account (11), one obtains:

[TABLE]

In conclusion, we can summarize the above procedure as follows (see (3)):

[TABLE]

In other words, we are speaking about the application of the following “continuous-stage Runge-Kutta-Nyström (RKN, hereafter) method” for solving problem (29), i.e., (4) :

[TABLE]

which provides the exact solution of the problem.

3.1 Polynomial approximation

As done for first order problems, also in this case we can consider a polynomial approximation $\sigma\in\Pi_{s}$ to $q$. This is done by resorting to the same finite vectors and matrices defined in (18) and (22), resulting into the following continuous-stage RKN method:

[TABLE]

which defines the application of the HBVM$(\infty,s)$ method for solving (4). One has, then,

[TABLE]

It is well-known [12, 14] that $q_{1}-q(h)=\dot{q}_{1}-\dot{q}(h)=O(h^{2s+1})$.666Also in this case, one could derive the result through the simplifying assumptions for continuous-stage RKN methods [37].

Remark 2

We observe, however, that in order for (32) to hold, one must have $s\geq 2$. Conversely, one would obtain $\bar{b}_{\xi}\equiv 1$, in place of  $\bar{b}_{\xi}=1-\xi$.

Moreover, considering that (compare with (30))

[TABLE]

one obtains:

[TABLE]

in place of (3).

3.1.1 Discretization

We conclude this section by recalling that, by approximating the integrals appearing in (36) by means of a Gauss-Legendre formula of order $2k$, one obtains a HBVM$(k,s)$ method, which retains the order $2s$ of the approximation defined by (36), for all $k\geq s$.777In particular, when $k=s$, one obtains the RKN method induced by the $s$-stage Gauss collocation method, $s\geq 2$. The Butcher tableau of this $k$-stage RKN method turns out to be given by:

[TABLE]

with $\circ$ the Hadamard (i.e., componentwise) product, and the same matrices and vectors defined in (22) and (26)–(28). As in the case of first order problems, one has that the entries of the Butcher matrix $\bar{A}$ in (39) are given by (see (38))

[TABLE]

for all $k\geq s$ and $s\geq 2$.

4 Conclusions

In this paper, we have studied the formulation of the class of energy-conserving methods named Hamiltonian Boundary Value Methods (HBVMs) as continuous-stage RK methods. When applied for solving special second-order problems, such methods also provide a class of continuous-stage RKN methods, whose derivation has been provided in full details. The formulation of HBVMs as continuous-stage RK/RKN methods, in turn, is interesting by itself, even though the efficient implementation and analysis of the methods is better addressed, in our opinion, in their original formulation (see, e.g., the monograph [12] or the review paper [14].)