Analysis of Energy and QUadratic Invariant Preserving (EQUIP) methods

TL;DR

This paper analyzes EQUIP methods, a class of geometric integrators that conserve energy and quadratic invariants, providing reformulation, refined analysis, and numerical comparisons with Gauss methods.

Contribution

It offers a reformulation and detailed analysis of EQUIP methods, including a practical implementation procedure and comparative numerical tests.

Findings

EQUIP methods conserve Hamiltonian and quadratic invariants.

Numerical tests show EQUIP methods' effectiveness on Hamiltonian problems.

Comparison indicates advantages over traditional Gauss methods.

Abstract

In this paper we are concerned with the analysis of a class of geometric integrators, at first devised in [14, 18], which can be regarded as an energy-conserving variant of Gauss collocation methods. With these latter they share the property of conserving quadratic first integrals but, in addition, they also conserve the Hamiltonian function itself. We here reformulate the methods in a more convenient way, and propose a more refined analysis than that given in [18] also providing, as a by-product, a practical procedure for their implementation. A thorough comparison with the original Gauss methods is carried out by means of a few numerical tests solving Hamiltonian and Poisson problems.

Tables

Table 1. Table 1: Numerical solution of the Kepler problem ( 88 )–( 90 ) by using the EQUIP(6,2) and the 2-stage Gauss (GAUSS2) methods with stepsize h = 2 ​ π / n ℎ 2 𝜋 𝑛 h=2\pi/n over 10 periods.

	EQUIP(6,2)							GAUSS2
$n$	error	rate	$H$-error	$M$-error	$\overline{\alpha }$	rate	iter/	error	rate	$H$-error	rate	$M$-error	iter/
							step						step
20	1.34e-1	–	1.64e-09	3.66e-15	1.51e-3	–	19.6	1.55e0	–	1.95e-3	–	5.75e-15	17.4
30	2.61e-2	4.0	6.10-12	7.88e-15	6.81e-4	2.0	15.6	2.37e-1	4.6	2.27e-4	5.3	1.15e-14	14.2
40	8.36e-3	4.0	1.86e-13	4.22e-16	3.84e-4	2.0	13.6	8.00e-2	3.8	7.65e-5	3.8	1.24e-15	12.8
50	3.45e-3	4.0	1.84e-14	1.58e-15	2.45e-4	2.0	12.5	3.41e-2	3.8	3.28e-5	3.8	2.39e-15	11.7
60	1.67e-3	4.0	2.05e-15	4.70e-16	1.70e-4	2.0	11.8	1.68e-2	3.9	1.61e-5	3.9	3.16e-15	11.3
70	9.01e-4	4.0	1.44e-15	4.28e-16	1.25e-4	2.0	11.4	9.17e-3	3.9	8.83e-6	3.9	4.10e-15	10.6
80	5.29e-4	4.0	1.11e-15	2.46e-15	9.58e-5	2.0	10.8	5.41e-3	3.9	5.22e-6	3.9	1.36e-15	10.3
90	3.31e-4	4.0	2.44e-15	8.11e-16	7.57e-5	2.0	10.5	3.40e-3	4.0	3.27e-6	4.0	3.33e-15	10.0
100	2.18e-4	4.0	9.77e-16	2.98e-15	6.13e-5	2.0	10.2	2.24e-3	4.0	2.16e-6	4.0	6.65e-15	9.7

Equations

$$\dot{y}=f(y),t>0,y(0)=y_0\in \Omega \subseteq {\mathbb{R}}^m.$$

$$C:\Omega \to {\mathbb{R}}^{\nu }$$

$$y(t)\text{mbox}solves(\text{ref}fy)\Rightarrow C(y(t))\equiv C(y_0),\forall t\ge 0.$$

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

$$\nabla C(y{)}^{\top }f(y)=0,\forall y\in \Omega .$$

$$0=C(y(t))-C(y(0))={\int }_0^t\nabla C(y(s){)}^{\top }\dot{y}(s)\mathrm{d}s,\forall t>0,$$

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

$${\cal P}_{s}=\left(\begin{array}[]{ccc}P_{0}(c_{1})&\dots&P_{s-1}(c_{1})\\ \vdots&&\vdots\\ P_{0}(c_{s})&\dots&P_{s-1}(c_{s})\end{array}\right),\qquad\Omega=\left(\begin{array}[]{ccc}b_{1}\\ &\ddots\\ &&b_{s}\end{array}\right),$$

$$P_s(c_i)=0,b_i={\int }_0^1{\ell }_i(x)\mathrm{d}x,{\ell }_i(x)=\prod _{j\ne i}\frac{x-c_j}{c_i-c_j},i=1,\dots ,s,$$

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

$$\begin{array}[]{c|c}{\bm{c}}&A:={\cal P}_{s}X_{s}{\cal P}_{s}^{\top}\Omega\\ \hline\cr\\[-8.53581pt] &{\bm{b}}^{\top}\end{array},\qquad{\bm{b}}=\left(\begin{array}[]{c}b_{1},\,\dots\,,b_{s}\end{array}\right)^{\top},\qquad{\bm{c}}=\left(\begin{array}[]{c}c_{1},\,\dots\,,c_{s}\end{array}\right)^{\top}.$$

$$\Omega A+A^{\top }\Omega -\mathbf{b}{\mathbf{b}}^{\top }=O,$$

$$\dot{y}=f(y)\equiv J\nabla H(y),\qquad J=\left(\begin{array}[]{rr}O_{d}&I_{d}\\ -I_{d}&O_{d}\end{array}\right)=-J^{\top}=-J^{-1},$$

$$H:\Omega \subseteq {\mathbb{R}}^m\to \mathbb{R}$$

$$\begin{array}[]{c|c}{\bm{c}}&A(\alpha):={\cal P}_{s}X_{s}(\alpha){\cal P}_{s}^{\top}\Omega\\ \hline\cr\\[-8.53581pt] &{\bm{b}}^{\top}\end{array},$$

$$X_{s}(\alpha):=X_{s}-\alpha W_{s}\equiv\left(\begin{array}[]{cccc}\xi_{0}&\alpha-\xi_{1}\\ \xi_{1}-\alpha&0&\ddots\\ &\ddots&\ddots&-\xi_{s-1}\\ &&\xi_{s-1}&0\end{array}\right),\quad W_{s}={\bm{e}}_{2}{\bm{e}}_{1}^{\top}-{\bm{e}}_{1}{\bm{e}}_{2}^{\top}=-W_{s}^{\top}.$$

$$||y_1(\alpha )-y(h)||=\alpha O(h^3)+O(h^{2s+1}).$$

$$|\alpha|\leq ch^{2},\qquad\mbox{for a given positive constant $c$}.$$

$$y_1:=y_1(\alpha )\approx y(h),$$

$$\sigma (0)=y_0,\dot{\sigma }(c_{i}h)=f(\sigma (c_{i}h)),i=1,\dots ,s,y_1\equiv y_1(0):=\sigma (h).$$

$$\sigma (0)=y_0,Y_i=\sigma (c_{i}h),i=1,\dots ,s.$$

$$\sigma :[0,h]\cup [h,h+1]\to {\mathbb{R}}^m,$$

$$y_1\equiv y_1(\alpha ):=\sigma (h+1),s.t.C(\sigma (1+h))=C(\sigma (0))\iff C(y_1)=C(y_0).$$

$${\sigma }_1(t)\equiv \sigma (t),t\in [0,h],{\sigma }_2(t)\equiv \sigma (t+h),t\in [0,1],$$

$${\gamma }_j:=\sum _{i=1}^s{b}_i{P}_j(c_i)f(Y_i),j=0,\dots ,s-1,$$

$${\bm{\phi}}_{i}\equiv\left(\begin{array}[]{c}\phi_{i0}\\ \vdots\\ \phi_{i,s-1}\end{array}\right):=X_{s}^{-1}{\bm{e}}_{i},\qquad i=1,2.$$

$${\sigma }_1(ch)=y_0+h\sum _{j=0}^{s-1}{\int }_0^c{P}_j(x)\mathrm{d}x\left[{\gamma }_j-\alpha \left({\varphi }_{2j}{\gamma }_0-{\varphi }_{1j}{\gamma }_1\right)\right],c\in [0,1],$$

$$Y:=\left(\begin{array}[]{c}Y_{1}\\ \vdots\\ Y_{s}\end{array}\right),$$

$$\sigma_{1}({\bm{c}}h):=\left(\begin{array}[]{c}\sigma_{1}(c_{1}h)\\ \vdots\\ \sigma_{1}(c_{s}h)\end{array}\right),$$

$$1=(1,\dots ,1{)}^{\top }\in {\mathbb{R}}^s,$$

Full text

Keywords: Gauss collocation methods, symplectic methods, energy-conserving methods, line integral methods, Hamiltonian problems, Poisson problems.

MSC: 65P10, 65L05.

1 Introduction

The study of numerical methods for differential equations able to retain relevant geometric properties of the solutions, is the core of what nowadays is named geometric integration (see, e.g., the classical monographs [40, 37, 29] and, more recently, [10, 4]). Accordingly, numerical methods able to reproduce specific geometric properties of the continuous dynamical system in the discrete one generated by their application, are often referred to as geometric integrators. In the framework of ODEs, this approach dates back to the work of G. Dahlquist, who devised the linear stability analysis of numerical methods for ODEs in order to assess their behavior when solving dissipative problems. In this respect, $A$-stable methods can be considered as geometric integrators, when dealing with such problems.

For a large class of problems which are of interest in the applications, the geometric properties of the solutions can be taken back to the presence of constants of motion, which are conserved along the evolution. Such problems are, therefore, named conservative problems. In more detail, let us consider an ODE-IVP in the form

[TABLE]

The problem admits $\nu<m$ constants of motion when there exists

[TABLE]

such that

[TABLE]

Hereafter, we shall assume that both $f$ and $C$ are suitably regular.111E.g., analytical. From a local point of view, the conservation property (3) is equivalent to require that

[TABLE]

i.e., since the initial point of the trajectory is arbitrary,

[TABLE]

An alternative way for requiring the conservation of $C(y(t))$ is through the use of a line integral,

[TABLE]

which is clearly equivalent to (4). The reformulation (6) of the invariance property is at the basis of the so-called (discrete) Line Integral Methods (see, e.g., the recent monograph [10]), and here we will exploit this new framework to carry out a theoretical analysis of the methods and to devise a practical implementation strategy.

The methods we are going to consider can be regarded as a suitable modification of the basic Gauss-Legendre collocation methods [24], which we rewrite here by using the well-known $W$-transformation [30]. In more details, introducing the orthonormal Legendre polynomial basis $\{P_{i}(x)\}$ on $[0,1]$ ($\Pi_{i}$ denotes the subspace of polynomials of degree at most $i$)

[TABLE]

the matrices

[TABLE]

with

[TABLE]

the Legendre nodes and weights, and matrix

[TABLE]

the $s$-stage Gauss method is defined by the Butcher tableau

[TABLE]

As is well known (see, e.g., [23]), such methods are algebraically stable, in that

[TABLE]

and, in particular, they preserve any quadratic invariant of (1). In the framework of geometric integration, condition (12) is equivalent to say that such methods are symplectic [36, 38]. The importance of symplectic methods relies to their application to canonical Hamiltonian systems, i.e. problems of the form

[TABLE]

with $m=2d$ and $O_{d}$ and $I_{d}$ the zero and identity matrices of dimension $d$, respectively (hereafter $I_{r}$ will denote the identity matrix of dimension $r$). Moreover,

[TABLE]

is the Hamiltonian function defining the problem, which we assume to be suitably regular. As is well-known, in such a case, the Hamiltonian function is a constant of motion for the dynamical system defined by (13)–(14). By virtue of (12) a symplectic method does preserve quadratic first integrals including the Hamiltonian function itself in case (13) is linear. The capability of conserving quadratic invariants has been recently recognized as one of the main advantages in using symplectic integrators (see [39]).

With this premise, assuming hereafter that only one constant of motion exists for (1) (i.e., $\nu=1$ in (2)), besides possible quadratic invariants,222We observe that the methods could in principle be generalized to cope with the conservation of more invariants (i.e., $\nu>1$), as sketched in [18]. Nevertheless, the arguments become more involved, so that we shall not consider this case, here. the methods we shall study have a Runge-Kutta representation in the form (cf. (11))

[TABLE]

where ${\cal P}_{s}$ and $\Omega$ are the same matrices defined in (8), and (see (10)), by setting ${\bm{e}}_{i}\in\mathbb{R}^{s}$ the $i$th unit vector,

[TABLE]

We observe that, for $\alpha=0$, (15)–(16) coincides with the Butcher tableau (11) of the $s$-stage Gauss collocation method. In the present case, the scalar parameter $\alpha$ is chosen, at each step of application of the method for numerically solving (1), in order to enforce the conservation of the invariant (2) (or (14)) in the numerical solution (see also [18], for the Hamitonian case). To this purpose, exploiting Theorem 5.1 in [31, Chapter IV.5] and the symmetry 333See Theorem 3 below. of the method (15), we see that it has order 2 for any fixed $\alpha\not=0$ and, clearly, order $2s$ when $\alpha=0$. More precisely, denoting by $y_{1}(\alpha)$ the solution yielded by the application of a single step of method (15) to problem (1), we have

[TABLE]

In Section 3 we will show that an $\alpha^{\ast}=O(h^{2s-2})$ actually exists such that $C(y_{1}(\alpha^{\ast}))=C(y_{0})$, which will confer an order $p=2s$ on the resulting method, as is the case for the corresponding Gauss formula. Since $s\geq 2$, we may confine the study of existence and uniqueness of such an $\alpha^{\ast}$ inside a ball of radius $ch^{2}$. Therefore in the sequel we make the following assumption

[TABLE]

The formulae resulting from the choice $\alpha=\alpha^{\ast}$ have been named Energy and QUadratic Invariants Preserving (EQUIP) methods in [18] since, beside preserving the Hamiltonian function, they continue to preserve any quadratic invariant just like the Gauss collocation methods they come from. It must be stressed that this double feature seems to confer great robustness to the methods, as we will see later in the numerical tests. For sake of completeness, we mention that the idea of modifying Runge-Kutta methods, by introducing a suitable parameter to enforce energy conservation, has been also considered, e.g., in [25, 35].

An aspect which has not yet been studied neither in [18] nor in the subsequent research concerns the implementation of the methods. Motivated by this issue, in the present work we exploit the characterization (6) of a constant of motion to show that the methods (15) may be indeed interpreted as discrete line integral methods. This will allow us to carry out a more refined theoretical analysis than that given in [18] and to devise a practical implementation strategy.

The structure of the paper is as follows: in Section 2 we provide full details about the derivation of the parameter $\alpha$; in Section 3 we provide a novel analysis for the corresponding method, besides that given in [18]; in Section 4 we consider the fully discrete method; Section 5 is devoted to the solution of the generated discrete problem; Section 6 provides numerical evidence for the theoretical findings; at last, Section 7 contains a few conclusions and further directions of investigation.

2 Energy and QUadratic Invariants Preserving (EQUIP)

methods

In this section, we derive an expression for the parameter $\alpha$ in (15)–(16), chosen in order to enforce the conservation of $C(y)$. Since we are speaking about a one-step method, it will suffice to consider the very first application of the method, in order to obtain the approximation

[TABLE]

where $h$ is the time step and $y(t)$ is the solution of (1). The derivation of the parameter $\alpha$ will be then derived in the framework of line integral methods, namely methods which enforce the conservation of the invariants through an approximation of the line integral (6). Such methods, at first derived for the numerical solution of Hamiltonian problems [32, 33, 34], eventually resulted in the class of Hamiltonian Boundary Value Methods (HBVMs) [6, 11, 12, 13, 15, 16, 17, 20], and have been then generalized along several directions [1, 2, 5, 7, 9, 19, 21, 41], including the methods here studied.444Despite the common line integral framework to gain the conservation property, the methods that we shall describe here greatly differ from HBVMs. As matter of fact, the latter methods are not symplectic and are conservative only for Hamiltonian problems, whereas the former ones are defined, at each integration step, by a symplectic map and can be conservative for more general problems.

To begin with, we start recalling that, when $\alpha=0$, the method defined by (15)–(16) reduces to the $s$-stage Gauss method, which is defined by the corresponding collocation polynomial $\sigma\in\Pi_{s}$ such that

[TABLE]

Denoting by $Y_{i}$, $i=1,\dots,s$ the stages associated with the method, an equivalent way of defining the collocation polynomial $\sigma(ch)$ is via the interpolation conditions

[TABLE]

When $\alpha\neq 0$, the collocation polynomial is replaced with a piecewise continuous polynomial path, which we continue to denote by $\sigma$,555For sake of notation simplicity, in the sequel we will omit to explicitly highlight the dependence of the continuous polynomial path $\sigma$ and some related quantities on the parameter $\alpha$, if this is clear from the context.

[TABLE]

such that the interpolation conditions (20) are still fulfilled for the stages $Y_{i}$ associated with (15) and, moreover,

[TABLE]

For sake of simplicity, let us denote

[TABLE]

and notice that the interpolation conditions (20) only involve the polynomial $\sigma_{1}$. The following results hold true.

Theorem 1

For fixed $\alpha$ and $h$, and with reference to the polynomial basis (7), let us denote

[TABLE]

where $Y_{i}$ are the stages associated with the Runge-Kutta method (15). Moreover, see (10), we define the vectors

[TABLE]

Then, the polynomial $\sigma_{1}\in\Pi_{s}$ defined as

[TABLE]

satisfies the interpolation conditions (20).

Proof The polynomial $\sigma_{1}$ defined at (25) clearly satisfies the first condition in (20). Setting

[TABLE]

the stage vector of (15)–(16), and

[TABLE]

the stage interpolation conditions in (20) read, in compact notation, $Y=\sigma_{1}({\bm{c}}h)$. Denoting

[TABLE]

one obtains:

[TABLE]

Now, considering the vector form of (23),

[TABLE]

and the property

[TABLE]

from (28) and (25) we deduce that

[TABLE]

This concludes the proof.  $\Box{~{}}$

Concerning the polynomial $\sigma_{2}(t)$ in (22), it will be used to gain the conservation property. In more details, one has:

[TABLE]

Consequently, we can choose a linear polynomial:

[TABLE]

We observe that, from (25) and (34), we obtain: 666Clearly, from (34)–(35) one obtains that $\sigma_{2}(t)\equiv y_{1}$, when $\alpha=0$.

[TABLE]

Concerning the choice of the parameter $\alpha$ such that the conservation condition $C(y_{1})=C(y_{0})$ in (21) is fulfilled, the following result holds true.

Theorem 2

Let us consider the polynomials $\sigma_{1}$ and $\sigma_{2}$ defined at (25) and (34), respectively, and the vectors $\gamma_{j}$ defined at (23). Morever, we define the vectors

[TABLE]

Then $C(y_{1})=C(y_{0})$, provided that:

[TABLE]

Proof By using the framework provided by line integral methods, from (22), (25), (34), and (35), one has:

[TABLE]

from which (37) immediately follows. $\Box{~{}}$

In case of the Hamiltonian problem (13)–(14), the previous result is modified as follows.

Corollary 1

Let us consider the polynomials $\sigma_{1}$ and $\sigma_{2}$ defined at (25) and (34), respectively, and the vectors $\gamma_{j}$ defined at (23). Moreover, with reference to (13), let us define the vectors

[TABLE]

Then $H(y_{1})=H(y_{0})$, provided that:

[TABLE]

Proof The expression (39) follows immediately from (37) by taking into account (13) and setting $C(y)=H(y)$. In fact, from (36) one obtains: $\gamma_{j}(\sigma_{1})=J\rho_{j}(\sigma_{1})$, $j=0,\dots,s-1$, and $\bar{\gamma}(\sigma_{2})=J\bar{\rho}(\sigma_{2})$. $\Box{~{}}$

The method defined by (15)–(16), with $\alpha$ given by (37) (i.e., (39), in the Hamiltonian case (13)), are named [18] $s$-stage EQUIP methods: the name being accounted for by their analysis in the next section. We conclude this section observing that the discrete line integral tool provides us with a new nonlinear system associated with an EQUIP method, where the stage vector $Y$ is replaced by the vector ${\bm{\gamma}}$. In fact, inserting (31) in (29) we obtain that the unknown vectors $\gamma_{j}$ are determined by solving the system

[TABLE]

with $f({\bm{1}}\otimes y_{0}+h({\cal I}_{s}-\alpha{\cal P}_{s}W_{s})\otimes I_{m}\,{\bm{\gamma}})$ denoting the block vector containing the application of the function $f$ to each block entry of the argument. This new representation of an EQUIP method made up of (40), and either (37) or (39), will be exploited in the sequel to derive a number of theoretical results. Compared with [18], where only the existence of an $\alpha$ yielding energy conservation was proved without giving details about how $\alpha$ should be searched, here equation (37) yields an implicit equation the parameter $\alpha$ must satisfy, and will provide us with a practical strategy for the implementation of the methods.

3 Analysis of EQUIP methods

In what follows, we shall obviously assume $s\geq 2$ in (15)–(16).777In general, $s\geq\nu+1$, in case of $\nu$ constants of motion. As already observed, when $\alpha=0$, the $s$-stage EQUIP method (15)–(16) reduces to the $s$-stage Gauss-collocation method. Another relevant property is symmetry.

Theorem 3

For any fixed $\alpha\in\mathbb{R}$ the $s$-stage Runge-Kutta method (15)–(16) is symmetric.

Proof Since there are no redundant stages, the symmetry of the method is equivalent to require (see (27))

[TABLE]

with

[TABLE]

The latter condition is trivial. Concerning the former, one has:

[TABLE]

Since (see (8)) $P\Omega P=\Omega$  and  $W_{s}^{\top}=-W_{s}$, we obtain

[TABLE]

By considering that $P_{j}(c_{i})=(-1)^{j}P_{j}(c_{s-i+1})$ and $W_{s}={\bm{e}}_{2}{\bm{e}}_{1}^{\top}-{\bm{e}}_{1}{\bm{e}}_{2}^{\top}$, one obtains:

[TABLE]

with $G=G^{\top}$. Consequently, matrix ${\cal P}_{s}G{\cal P}_{s}^{\top}$ is symmetric and, therefore:

[TABLE]

A further relevant property of the $s$-stage method (15)–(16) is algebraic stability.

Theorem 4

For any fixed $\alpha\in\mathbb{R}$ the $s$-stage Runge-Kutta method (15)–(16) satisfies (12).

Proof In fact, for any fixed $\alpha\in\mathbb{R}$, from (8)–(10) and (15)–(16) one has:

[TABLE]

As a straightforward consequence, one obtains the following result.

Corollary 2

For any fixed $\alpha\in\mathbb{R}$ the $s$-stage Runge-Kutta method (15)–(16) is symplectic and, therefore, conserves all quadratic invariants of (1).

Proof See [23, 36, 38]. $\Box{~{}}$

Remark 1

The name Energy and QUadratic Invariants Preserving (EQUIP) method derives from the conservation properties associated with Theorem 2 (i.e., Corollary 1, in the Hamiltonian case), and Corollary 2.

We notice that all the properties discussed in Theorems 3 and 4, and in Corollary 2 above, hold true independently of the choice of the parameter $\alpha$. If we let equation (37) (or its variant (39)) come into play, a primary issue is to ascertain the existence and uniqueness of the solution of the nonlinear system associated with the method, namely (37)–(40) (or (39)–(40)). For this purpose, we need some additional preliminary results.

Lemma 1

Let $G:\mathbb{R}\rightarrow V$, vith $V$ a vector space, be a function admitting a Taylor expansion at 0. Then,

[TABLE]

Proof See [17, Lemma 1]. $\Box{~{}}$

Lemma 2

With reference to (23), (36), and (38), one has, under regularity assumptions on $f$ and $C$:

[TABLE]

and

[TABLE]

Proof The thesis easily follows from Lemma 1. $\Box{~{}}$

We need now to discuss the right-hand sides of the expressions at (37) and (39): we shall then study the corresponding numerators and denominators.

Analysis of the numerators in (37) and (39)

Hereafter, let us denote ${\cal N}(h)$ the numerator at the right-hand side of (37) or (39).

Lemma 3

Concerning the numerator in (37), under regularity assumptions on $f$ and $C$, one has

[TABLE]

Proof Following the arguments in [9], let us consider the following expansions along the orthonormal basis (7):

[TABLE]

with $\rho_{j}(\sigma_{1})$ and $\gamma_{j}(\sigma_{1})$ given by (36) and (38), respectively. Because of (5) and (7), one has, then:

[TABLE]

Consequently, by virtue of (41), one obtains

[TABLE]

Then,

[TABLE]

For the particular case of Hamiltonian problems, one has the following simpler result.

Lemma 4

Concerning the numerator in (39), under regularity assumptions on $H$, one has

[TABLE]

Proof In fact, from (41), one has, by taking into account that matrix $J$ is skew-symmetric:

[TABLE]

Analysis of the denominators in (37) and (39)

Before studying the denominators at the right-hand sides of (37) and (39), which we shall hereafter denote by ${\cal D}(h)$, we need to state the following result concerning the leading entries of the vectors defined at (24).

Lemma 5

With reference to (10) and (24), one has (discarding $\phi_{ij}$, when $j\geq s$):

[TABLE]

Moreover, it is quite straightforward, though lengthy, to prove the following expansions, which hold under the assumption (18).

Lemma 6

Let us consider the vectors defined at (23), (36), and (38), and the parameter $\alpha$ defined at either (37) or (39). Moreover, let us set

[TABLE]

Then, under regularity assumptions on the involved functions, and assuming that (18) is satisfied, the following expansions hold true for the following coefficients defined at (23), (36), and (38): 888They are obtained by expanding $f(y)$ and $\nabla C(y)$ at $y=y_{0}$.

[TABLE]

and

[TABLE]

The next results then follow.

Lemma 7

Let us consider the denominator at the right-hand side in (37). Then, under regularity assumptions on $f$ and assuming (18) be satisfied, one has:

[TABLE]

Proof We set  ${\cal D}(h)=D_{1}(h)+D_{2}(h)$,  with

[TABLE]

Then, by taking into account (41) and (48)–(55), we discuss separately the case $s$ odd and $s$ even.

$s$ odd: in such a case, one has, by taking into account that $\phi_{22}\xi_{1}\xi_{2}=-\xi_{1}$,

[TABLE]

Consequently, (56) follows.

$s$ even: in such a case, one has, by taking into account that  $(\xi_{0}^{-1}\xi_{1}-\xi_{0}\xi_{1}^{-1})=-\xi_{0}^{-2}\xi_{1}$  and  $(\xi_{0}^{2}\xi_{1}^{-1}-\xi_{0}^{-2}\xi_{1})=-\xi_{1}$,

[TABLE]

from which (56) follows. $\Box{~{}}$

In the Hamiltonian case, by considering (13)–(14)  (i.e., $Jg(y)\equiv f(y)$), and the expansions (54), one proves the following result in a similar way.

Lemma 8

Let us consider the denominator at the right-hand side in (39). Then, under regularity assumptions on $f$ and assuming (18) be satisfied, one has:

[TABLE]

Existence and uniqueness

As a simple consequence of the previous Lemmas 3–8, one obtains the following result, which provides a different, and more general, proof than what proved in [18, Theorem 3.5] in the Hamiltonian case.

Theorem 5

Assume that $\alpha$ satisfies equation (37) and

[TABLE]

Then,

[TABLE]

The same result holds true in the Hamiltonian case, when $\alpha$ is given by the right-hand side in (39), and

[TABLE]

Proof Depending on the two cases, one obtains $\alpha={\cal N}(h)/{\cal D}(h)$, with ${\cal N}(h)$ and ${\cal D}(h)$ defined according to (43) and (56), when $\alpha$ is given by (37), or (44) and (57), when $\alpha$ is given by (39). In both cases, (59) easily follows, due to the fact that ${\cal N}(h)=O(h^{2s})$ and ${\cal D}(h)=O(h^{2})$. $\Box{~{}}$

Remark 2

We observe that when (58) is not true, then the leading term in (56) vanishes. Practically, this event can be recognized by the slow down or failure of the convergence of the nonlinear iteration later described in Section 5 which, in turn, relies on the result of Theorem 6 below. In such a case, one may set $\alpha=0$, thus performing a single step of the basic $s$-stage Gauss method. The same argument applies, in the Hamiltonian case, when (60) is not satisfied, so that the leading term in (57) vanishes. This strategy will allow to retain the order $2s$ of the underlying $s$-stage Gauss method, as we shall see later in the analysis of the method.

We are now in the right position to prove the existence and uniqueness result. The nonlinear system to be solved at each step of the integration procedure is (see (37)–(39) and (40)) 999We now emphasize, in the first equation in (61), the dependence of the numerator and denominator defining $\alpha$ from the given arguments, rather than from the stepsize $h$.

[TABLE]

Setting $\alpha=0$ in the first equation of (61) leads back to the nonlinear system associated with the Gauss method of order $2s$ that we write as

[TABLE]

Assuming $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ Lipschitz continuous with constant $L>0$, we see that $\Psi_{0}({\bm{\gamma}})$ is Lipschitz continuous in $\mathbb{R}^{sm}$ with constant 101010Hereafter, $\|\cdot\|$ denotes the infinity norm.

[TABLE]

therefore $\Psi_{0}({\bm{\gamma}})$ is a contraction on $\mathbb{R}^{sm}$ for $h$ sufficiently small, and the Banach fixed-point theorem assures the existence of a unique solution $\bar{\bm{\gamma}}$ of (62) which can be determined as the limit of the sequence ${\bm{\gamma}}^{\ell+1}=\Psi_{0}({\bm{\gamma}}^{\ell})$, $\ell=0,1,\dots,$ starting from any ${\bm{\gamma}}_{0}\in\mathbb{R}^{sm}$.

Theorem 6

Under regularity assumptions on the functions $f$, and assuming either (58) or (60) holds true, so that $D(\alpha,{\bm{\gamma}})$ is bounded away from zero in a neighbourhood of $(0,\bar{\bm{\gamma}})$, system (61) admits a unique solution $(\alpha^{\ast},{\bm{\gamma}}^{\ast})$ which can be found as the limit of the sequence

[TABLE]

starting at

[TABLE]

Moreover the solution $(\alpha^{\ast},{\bm{\gamma}}^{\ast})$ satisfies

[TABLE]

Proof In accord with the contraction mapping theorem, we show that constants $h_{0}>0$ and $r>0$ exist such that, for $h\in[0,h_{0})$, the iteration function,

[TABLE]

defined at (64) satisfies the two properties:

• (a)

$F(\alpha,{\bm{\gamma}})$ is a contraction in the closed ball $B((0,\bar{\bm{\gamma}}),r)$ of center $(0,\bar{\bm{\gamma}})$ and radius $r>0$;

• (b)

the sequence $(\alpha^{\ell},{\bm{\gamma}}^{\ell})\in B((0,\bar{\bm{\gamma}}),r)$,  for all  $\ell=0,1,\dots$.

Expanding $\Psi(\alpha,{\bm{\gamma}})$ with respect to the variable $\alpha$ in a neighbourhood of zero yields

[TABLE]

with

[TABLE]

where $f^{\prime}$ stands for the Jacobian matrix of $f$. The function $R(\alpha,{\bm{\gamma}})$ satisfies $R(0,{\bm{\gamma}})=0$ and, since by virtue of Theorem 5 we may assume $|\alpha|\leq c_{1}h^{2s-2}$, we also get

[TABLE]

with $c_{1}$ and $c_{2}$ positive constants independent of $h$, $\alpha$, and ${\bm{\gamma}}$, in a closed ball $B((0,\bar{\bm{\gamma}}),r)$ of center $(0,\bar{\bm{\gamma}})$ and given radius $r>0$. Consequently, in a neighbourhood of $(0,\bar{\bm{\gamma}})$, system (61) may be regarded as a perturbation of system (62), the perturbing term being the function $R(\alpha,{\bm{\gamma}})$.

From (67), it follows that $R(\alpha,{\bm{\gamma}})$ is Lipschitz continuous with constant $K_{1}=O(h^{2s-1})$ and hence $\Psi(\alpha,{\bm{\gamma}})$ is Lipschitz continuous with constant $\mu\leq K_{0}+K_{1}=O(h)$, where $K_{0}$ is defined in (63). This property, together with the Lipschitz continuity of $N(\alpha,{\bm{\gamma}})/D(\alpha,{\bm{\gamma}})$ with constant $O(h^{2s-2})$, assures that $h_{0}>0$ exists such that, for $h\in[0,h_{0})$, the iteration function $F(\alpha,{\bm{\gamma}})$ defined at (66) is indeed a contraction in $B((0,\bar{\bm{\gamma}}),r)$, with constant $\mu$.

To show that system (61) actually admits a solution ${\bm{\gamma}}^{\ast}$ which is $O(h^{2s-1})$-close to $\bar{\bm{\gamma}}$, we set

[TABLE]

and prove that the sequence $(\alpha^{\ell},{\bm{\gamma}}^{\ell})$ defined at (64) is entirely contained in the closed ball $B((0,\bar{\bm{\gamma}}),r)$. We proceed by induction. For $\ell=1$ we get

[TABLE]

and

[TABLE]

Assuming now that $(\alpha^{\ell},{\bm{\gamma}}^{\ell})\in B((0,\bar{\bm{\gamma}}),r)$, we obtain

[TABLE]

and

[TABLE]

Since $(\alpha^{\ast},{\bm{\gamma}}^{\ast})\in B((0,\bar{\bm{\gamma}}),r)$, then clearly $|\alpha^{\ast}|\leq r=O(h^{2s-2})$. Moreover, subtracting $\bar{\bm{\gamma}}=\Psi(0,\bar{\bm{\gamma}})$ from ${\bm{\gamma}}^{\ast}=\Psi(\alpha^{\ast},{\bm{\gamma}}^{\ast})$ we finally obtain

[TABLE]

Order of convergence

Next, we show that method (15)–(16), with $\alpha$ given by (37) (or (39), in case the problem is Hamiltonian), retains the order $2s$ of the underlying $s$-stage Gauss method (i.e., the method corresponding to $\alpha=0$).

Theorem 7

Let $y_{1}\approx y(h)$ be the approximation provided by the $s$-stage EQUIP method (15)–(16), with $\alpha$ given by (37) (respectively, (39) when the problem is Hamiltonian). Then, assuming that the hypotheses of Theorem 6 hold true,

[TABLE]

i.e., the method has order $2s$.

Proof The proof will follow the novel approach defined in [17]. For this purpose, let us denote by $y(t,\theta,\xi)$ the solution of the problem

[TABLE]

Moreover, let $\Phi(t,\theta)$ be the fundamental matrix solution of the corresponding variational problem

[TABLE]

We also recall that, according to Lemma 1,

[TABLE]

By taking into account (20)–(34), let us then define the polynomial

[TABLE]

such that $u(0)=y_{0}$ and $u(h)=y_{1}$. Consqeuently,

[TABLE]

By considering that

[TABLE]

and taking into account (41), (69), and (70), one obtains:

[TABLE]

Due to the fact that $\Gamma_{j}(h)=O(h^{j})$, $\gamma_{j}(u)=O(h^{j})$, $\gamma_{j}-\gamma_{j}(\sigma_{1})=O(h^{2s-j})$, one then concludes that

[TABLE]

Consequently, (68) follows provided that

[TABLE]

By taking into account (25), (70), and (59), one has

[TABLE]

so that:

[TABLE]

having set

[TABLE]

Since $g_{j}$ is a polynomial of degree $s+j$, for $j=0,\dots,s-1$, one has that

[TABLE]

Consequently, by setting

[TABLE]

and taking into account (8), (30), (24), (29), one has:

[TABLE]

By considering that ${\cal P}_{s}^{\top}\Omega{\cal P}_{s}=I_{s}$, one then obtains:

[TABLE]

(73) then follows by taking into account that $\alpha=O(h^{2s-2})$, because of Theorem 6, and, moreover,

[TABLE]

4 Discretization of the involved integrals

We observe that the actual implementation of the $s$-stage EQUIP method (15)–(16) requires the evaluation of the integrals (36) and (38). Consequently, it is not yet an operative method since, quoting e.g., Dahlquist and Björk [27, p. 521], as is well known, even many relatively simple integrals cannot be expressed in finite terms of elementary functions, and thus must be evaluated by numerical methods. For this purpose, one may use a $k$-point Gauss-Legendre quadrature, and we denote EQUIP$(k,s)$ the corresponding method, whose analysis is now completed, w.r.t. that done in the previous section. Hereafter, we shall assume $k\geq s$. Consequently, by setting

[TABLE]

the nodes and weights of the Gauss-Legendre formula of order $2k$, one approximates (36) and (38) by means of

[TABLE]

and

[TABLE]

respectively. It is quite straightforward to prove that, under regularity assumptions on $f$ and $C$,

[TABLE]

As a result, the new approximation is formally still defined by the two polynomials $\sigma_{1}$ and $\sigma_{2}$ defined in (25) and (32)–(34). However, in such formulae, the parameter $\alpha$, defined in (37) or (39), has to be replaced by an approximated one, say $\hat{\alpha}$, computed by using the quadratures (77)–(79). That is,

[TABLE]

or

[TABLE]

respectively. The following result can then be proved.

Lemma 9

Let the parameter $\alpha$ be defined by either (81), when solving problem (1)–(2), or (82), when solving problem (13)–(14). Then, for all $k\geq s$ one has

[TABLE]

with $\alpha$ given by (37) and (39), respectively. In particular, when the problem is Hamiltonian and $k=s$, one has $\hat{\alpha}=0$, so that EQUIP$(s,s)$ reduces to the $s$-stage Gauss-method.

Proof The first part of the proof follows straightforwardly, by taking into account (80), from the line integral formulation used to derive (37) and (39), respectively. In particular, one finds that

[TABLE]

with ${\cal N}(h)=O(h^{2s})$ and ${\cal D}(h)=O(h^{2})$ formally given by (43) and (56), in case of problem (1)–(2), or (44) and (57), in case of problem (13)–(14), respectively.

The fact that $\hat{\alpha}=0$, when the problem is Hamiltonian and $k=s$, follows from (82), by considering that matrix $J$ is skew-symmetric and (see (23) and (79)) $\hat{\gamma}_{j}(\sigma_{1})\equiv\gamma_{j}$. $\Box{~{}}$

Moreover, the following results hold true, whose straightforward proofs, which strictly follow that of Theorem 2, are omitted for sake of brevity.

Theorem 8

With reference to problem (1)–(2) and under regularity assumptions on both $f$ and $C$, and for $k\geq s$, one obtains that the EQUIP$(k,s)$ method defined by the parameter $\hat{\alpha}$ in (81) satisfies:

[TABLE]

Theorem 9

With reference to problem (13)–(14) and under regularity assumptions on $H$, and for $k\geq s$, one obtains that the EQUIP$(k,s)$ method defined by the parameter $\hat{\alpha}$ in (82) satisfies:

[TABLE]

Remark 3

As a consequence of Theorems 8 and 9, one has that an exact conservation can be always gained, in case of a polynomial invariant, by taking $k$ large enough. On the other hand, even in the general case a practical conservation can still be obtained, by taking $k$ large enough so that the quadrature error is within round-off error level. It is worth mentioning that by considering higher values of $k$ does not increase the complexity of the nonlinear problem to be solved, which amounts to compute the $s$ coefficients of the polynomial $\sigma_{1}$, and the parameter $\alpha$.

The previous results allows us to prove that, for all $k\geq s$, the EQUIP$(k,s)$ method still has order $2s$.

Theorem 10

Let $y_{1}\approx y(h)$ be the approximation provided by the EQUIP$(k,s)$ method (15)–(16), with $\alpha=\hat{\alpha}$ given by (81) (respectively, (82) when the problem is Hamiltonian). Then, for all $k\geq s$, and assuming the hypotheses of Theorem 7 hold true,

[TABLE]

i.e., the method retains the order $2s$.

Proof The proof strictly follows that done for Theorem 7, the only difference being that, starting from (74), $\alpha$ has to be formally replaced by $\hat{\alpha}$. Nevertheless, because of (83), the same arguments of that proof apply also in this case, so that (85) follows. $\Box{~{}}$

Remark 4

It is worth noticing that, by taking into account that $\hat{\alpha}=O(h^{2s-2})$, from (31) it follows that the stage order of EQUIP$(k,s)$ methods is $s$, for all $k\geq s$. I.e., the methods also retain the same stage order as that of the underlying $s$-stage Gauss method.

5 The nonlinear iteration

Next, let us duscuss the use of a fixed-point iteration for solving the discrete problem generated by the EQUIP$(k,s)$ method, with $k\geq s$. In principle, this is given by (64), with the only difference that we have to consider $\hat{\alpha}$ in place of $\alpha$. Therefore, by taking into account (83)–(84), convergence can be proved by using the same arguments in the proof of Theorem 6, provided that $k>s$, which we assume hereafter. There are, however, two improvements over the basic iteration (6), as explained below.

First of all, in order to improve the convergence properties of the iteration, we consider the following Gauss-Seidel type variant of (64),

[TABLE]

where we have formally used $\hat{N}$ and $\hat{D}$ in place of $N$ and $D$, respectively, in order to take into account of the quadratures (see (77)–(82)).

Secondly, the quadratures involved in the computation of $\hat{\alpha}$ in place of $\alpha$, potentially introduce a $O(h^{2k+1})$ error in the invariant at each step, according to Theorems 8 and 9. In order to avoid their accumulation, it suffices to compute $\hat{\alpha}$ so that the variation of the invariant, at the current step, matches the opposite of the error until the previous step. This latter error, which we denote $\Delta C$, can be easily computed by evaluating the invariant. As a result, following steps similar to those used in the proof of Theorem 2, one obtains that the iteration (86) becomes:

[TABLE]

In so doing, the errors in the invariant remain uniformly $O(h^{2k+1})$, since their accumulation is prevented, and moreover, the magnitude of $\hat{\alpha}$ remains $O(h^{2(s-1)})$.

6 Numerical results

In this section we consider a few numerical tests concerning the numerical approximation of some conservative problems. In particular, two Hamiltonian problems and two Poisson problems. They are listed below.

The well-known Kepler problem [29, 10], which is Hamiltonian, i.e. in the form (13), with

[TABLE]

It admits also the additional quadratic invariant (i.e., the angular momentum),

[TABLE]

Its solution is periodic of period $T=2\pi$, and resulting into an ellipse of eccentricity $\varepsilon$ in the $q$-plane, when starting at

[TABLE]

In particular, we shall consider the value $\varepsilon=0.5$. 2. 2.

The pendulum problem [10], which is also Hamiltonian with

[TABLE]

Its solution is periodic of period $T\approx 28.57109480185544$, when starting at

[TABLE] 3. 3.

The Poisson problem [5]

[TABLE]

with $y\in\mathbb{R}^{3}$ and

[TABLE]

also admitting, as a further invariant besides $H$, the quadratic Casimir

[TABLE]

The solution turns out to be periodic, of period $T\approx 0.53102669598427$, when choosing

[TABLE] 4. 4.

The Lotka-Volterra problem, which is in the form (93) with $y\in\mathbb{R}^{2}$ and

[TABLE]

The solution turns out to be periodic, of period $T\approx 7.720315563434113$, when choosing

[TABLE]

For such problems, we compare the EQUIP$(k,s)$ methods with the corresponding $s$-stage Gauss collocation methods, $s=2,3$, using a given constant stepsize of the form $h=T/n$, with $T$ the period of the solution and increasing values of $n$. This allows us to easily asses the accuracy of the solution, by checking it at the end of each period. As a measure of the error on the Hamiltonian, by setting $H_{i}$ its numerical value at the $i$th integration step, we consider the mean square value of the error over $N$ steps,

[TABLE]

A similar measure will be used for the invariants (89), $e_{M}$, and (95), $e_{C}$. When needed, by setting $\alpha_{i}$ the computed parameter (81)–(82) for the EQUIP method, we set

[TABLE]

as a measure of its magnitude. The EQUIP methods are implemented via the iteration (87), whereas only the second equation, having fixed $\hat{\alpha}=0$, is used for the Gauss methods. For the EQUIP$(k,s)$ methods, the value of $k$ is chosen in order to obtain a practical energy-conservation, as soon as the stepsize is decreased, for the considered values of $s$: in the present case, $k=6$ is sufficient for our purposes.

We start considering the Kepler problem (88)–(90). In Table 2 we list the obtained results by solving the problem with the EQUIP(6,2) and the 2-stage Gauss (GAUSS2, hereafter) methods. From such results, one deduces that:

• •

both methods are 4-th order, with EQUIP(6,2) one order of magnitude more accurate than GAUSS2, and preserve the quadratic invariant $M$;

• •

for EQUP(6,2), $\bar{\alpha}$ turns out to be $O(h^{2})$ as expected and, moreover, the error on $H$ soon falls within round-off error level, whereas it decreases with order 4 for GAUSS2;

• •

the number of iterations per step (iter/step) is essentially the same for both methods, and it decreases (as is expected) with the stepsize $h$.

Similarly, in Table 2 we list the obtained results by solving the problem with the EQUIP(6,3) and the 3-stage Gauss (GAUSS3, hereafter) methods, showing that:

• •

both methods are 6-th order, with EQUIP(6,3) one order of magnitude more accurate than GAUSS3, and preserve the quadratic invariant $M$;

• •

for EQUP(6,3), $\bar{\alpha}$ turns out to be $O(h^{4})$ as expected and, moreover, the error on $H$ soon falls within round-off error level, whereas it decreases with order 6 for GAUSS3;

• •

the number of iterations per step (iter/step) is essentially the same for both methods, and it decreases (as is expected) with the stepsize $h$.

Next, we briefly study the numerical solution of the pendulum problem (91)–(92) by means of EQUIP$(k,s)$ and the corresponding Gauss method. In particular:

• •

in Table 4 we list the obtained results for the EQUIP(6,2) and GAUSS2 methods, showing that, even though GAUSS2 is able to conserve the Hamiltonian with the prescribed order, nonetheless, the EQUIP(6,2) method is much more accurate in conserving the Hamiltonian, and this reflects into a much more accurate numerical solution;

• •

in Table 4 we list the obtained results for the EQUIP(6,3) and GAUSS3 methods, showing that also in this case, the energy-conserving EQUIP(6,3) is much more reliable than the GAUSS3 method.

We now consider the numerical solution of the Poisson problem (93)–(96) by means of the EQUIP(6,2) and GAUSS2 methods, and the EQUIP(6,3) and GAUSS3 methods, over 50 periods, by using a stepsize $h=T/100$. In Figure 2 are the computed results for the 4th-order (left-plot) and 6th-order (right-plot) methods, all of them exactly conserving the quadratic Casimir (95). As one may see, the EQUIP methods exhibit a linear error growth, whereas the Gauss methods show a quadratic error growth.

A similar result is obtained when numerically solving the Lotka-Volterra problem (97)–(98) by means of the EQUIP(6,2) and GAUSS2 methods, and the EQUIP(6,3) and GAUSS3 methods, over 50 periods, by using a stepsize $h=T/100$. In Figure 2 are the computed results for the 4th-order (left-plot) and 6th-order (right-plot) methods. As one may see, also in this case the EQUIP methods exhibit a linear error growth, whereas the Gauss methods show a quadratic error growth.

7 Conclusions

In this paper we reported a novel analysis for EQUIP($k,s$) methods, a symplectic variant of the $s$-stage Gauss collocation methods with energy-preserving properties. Such methods, which retain the order $2s$ of the underlying Gauss method, were originally devised for Hamiltonian problems and are here extended to cope with general conservative problems. The practical implementation of the methods has been another major focus of the paper, ranging from the discretization of the involved integrals, by means of suitable high-order Gaussian rules, to the solution of the generated discrete problem. In particular, the fixed-point iteration used for the analysis of the methods has been improved, in order to maintain a uniform error level in the invariant. A few numerical tests, on a couple of Hamiltonian and Poisson problems, confirm the theoretical findings and show that such methods can be far more reliable than the original Gauss formulae, due to their better conservation properties.

As a future direction of investigation we mention that, upon devising a Newton-type iteration more efficient than (86), EQUIP($k,s$) methods could be conveniently applied also for numerically solving Hamiltonian PDEs, after a proper space semi-discretization. In such a case, in fact, they will conserve possible quadratic invariants of the solution, besides the (semi)-discrete Hamiltonian (see, e.g., [2]).