Line Integral solution of Hamiltonian PDEs

TL;DR

This paper introduces energy-conserving line integral methods within the Hamiltonian Boundary Value Methods class to numerically solve Hamiltonian PDEs, demonstrating their effectiveness on key equations like wave, Schrödinger, and KdV.

Contribution

It presents a novel approach using line integral methods for Hamiltonian PDEs, emphasizing energy conservation and broad applicability to important equations.

Findings

Effective energy conservation in numerical solutions

Successful application to wave, Schrödinger, and KdV equations

Potential for improved long-term stability in simulations

Abstract

In this paper, we report about recent findings in the numerical solution of Hamiltonian Partial Differential Equations (PDEs), by using energy-conserving line integral methods in the Hamiltonian Boundary Value Methods (HBVMs) class. In particular, we consider the semilinear wave equation, the nonlinear Schr\"odinger equation, and the Korteweg-de Vries equation, to illustrate the main features of this novel approach.

Figures

Tables

Table 1. Table 1: Numerical solution of the sine-Gordon equation ( 92 ) using a time-step h = 100 / n ℎ 100 𝑛 h=100/n .

Gauss 1
$n$	time	$e_u$	rate	$e_H$	rate
2000	2.1	4.61e-02	—	1.98e-03	—
3000	2.8	2.05e-02	2.0	8.79e-04	2.0
4000	3.8	1.15e-02	2.0	4.95e-04	2.0
5000	4.8	7.37e-03	2.0	3.17e-04	2.0
6000	6.1	5.12e-03	2.0	2.20e-04	2.0
Gauss 2
$n$	time	$e_u$	rate	$e_H$	rate
1000	2.4	2.69e-05	—	2.57e-06	—
1500	3.3	5.28e-06	4.0	5.05e-07	4.0
2000	3.9	1.67e-06	4.0	1.59e-07	4.0
2500	5.3	6.83e-07	4.0	6.53e-08	4.0
3000	6.6	3.29e-07	4.0	3.15e-08	4.0
HBVM(4,1)
$n$	time	$e_u$	rate	$e_H$
1000	2.4	1.37e-02	—	7.11e-15
1500	3.3	6.15e-03	2.0	1.07e-14
2000	4.2	3.47e-03	2.0	8.88e-15
2500	5.7	2.22e-03	2.0	1.07e-14
3000	7.0	1.55e-03	2.0	8.88e-15
HBVM(4,2)
$n$	time	$e_u$	rate	$e_H$
1000	2.8	2.11e-05	—	1.07e-14
1500	4.0	4.18e-06	4.0	8.88e-15
2000	4.6	1.32e-06	4.0	1.07e-14
2500	6.0	5.42e-07	4.0	7.11e-15
3000	7.0	2.61e-07	4.0	8.88e-15
SHBVM
$n$	time	$k$	$s$	$e_u$	$e_H$
50	2.7	22	20	2.87e-12	3.55e-15
75	1.6	20	18	3.61e-13	7.11e-15
100	1.3	15	12	3.53e-13	3.55e-15

Table 2. Table 2: Numerical solution of the nonlinear Schrödinger equation ( 94 ) using a time-step h = 20 / n ℎ 20 𝑛 h=20/n .

Gauss 1
$n$	time	$e_{uv}$	rate	$e_1$	$e_2$	$e_H$	rate
400	19.1	4.93e-01	—	1.60e-14	2.01e-16	9.58e-04	—
600	26.3	2.32e-01	1.9	4.75e-14	6.70e-16	1.79e-04	4.1
800	33.2	1.31e-01	2.0	2.38e-14	1.39e-16	5.55e-05	4.1
1000	40.7	8.41e-02	2.0	2.22e-14	5.20e-17	2.25e-05	4.0
Gauss 2
$n$	time	$e_{uv}$	rate	$e_1$	$e_2$	$e_H$	rate
400	41.3	1.71e-03	—	3.13e-14	3.12e-17	5.20e-08	—
600	60.4	3.41e-04	4.0	1.91e-14	2.43e-17	2.14e-09	7.9
800	72.1	1.08e-04	4.0	1.91e-14	2.78e-17	2.18e-10	7.9
1000	84.7	4.44e-05	4.0	2.00e-14	3.82e-17	3.69e-11	8.0
HBVM(2,1)
$n$	time	$e_{uv}$	rate	$e_1$	rate	$e_2$	rate	$e_H$
400	36.3	5.23e-01	—	1.40e-04	—	4.25e-06	—	3.55e-15
600	50.7	2.45e-01	1.9	2.64e-05	4.1	7.96e-07	4.1	3.55e-15
800	60.1	1.38e-01	2.0	8.23e-06	4.1	2.47e-07	4.1	4.00e-15
1000	74.5	8.89e-02	2.0	3.35e-06	4.0	1.01e-07	4.0	4.00e-15
HBVM(4,2)
$n$	time	$e_{uv}$	rate	$e_1$	rate	$e_2$	rate	$e_H$
400	43.2	1.74e-03	—	6.66e-09	—	1.83e-10	—	4.44e-15
600	61.4	3.47e-04	4.0	2.70e-10	7.9	7.48e-12	7.9	3.55e-15
800	77.0	1.10e-04	4.0	2.74e-11	8.0	7.61e-13	7.9	4.00e-15
1000	92.3	4.52e-05	4.0	4.63e-12	8.0	1.29e-13	8.0	3.55e-15
SHBVM
$n$	time	$k$	$s$	$e_{uv}$	$e_1$	$e_2$	$e_H$
50	55.0	20	18	3.13e-11	1.35e-14	6.59e-17	4.88e-15
75	53.6	16	14	2.27e-11	1.33e-14	7.29e-17	3.11e-15
100	61.6	14	12	2.47e-11	1.40e-14	6.25e-17	3.11e-15

Table 3. Table 3: Numerical solution of the Korteweg–de Vries equation ( 96 ) using a time-step h = 10 / n ℎ 10 𝑛 h=10/n .

Gauss 1
$n$	time	$e_u$	rate	$e_H$	rate
10000	5.1	1.10e+00	—	9.38e-07	—
20000	8.6	2.75e-01	2.0	5.85e-08	4.0
30000	13.8	1.22e-01	2.0	1.16e-08	4.0
40000	17.6	6.88e-02	2.0	3.67e-09	4.0
50000	20.7	4.40e-02	2.0	1.50e-09	4.0
Gauss 2
$n$	time	$e_u$	rate	$e_H$	rate
10000	9.4	9.33e-05	—	7.30e-12	—
20000	18.1	5.84e-06	4.0	1.99e-13	5.2
30000	25.3	1.15e-06	4.0	2.84e-13	***
40000	31.3	3.65e-07	4.0	6.54e-13	***
50000	39.5	1.50e-07	4.0	7.25e-13	***
HBVM(2,1)
$n$	time	$e_u$	rate	$e_H$
10000	9.4	9.65e-01	—	6.39e-14
20000	16.6	2.42e-01	2.0	5.68e-14
30000	21.8	1.08e-01	2.0	7.11e-14
40000	29.0	6.05e-02	2.0	6.39e-14
50000	33.6	3.87e-02	2.0	6.39e-14
HBVM(3,2)
$n$	time	$e_u$	rate	$e_H$
10000	12.8	8.49e-05	—	5.68e-14
20000	24.4	5.32e-06	4.0	5.68e-14
30000	34.6	1.05e-06	4.0	6.39e-14
40000	43.0	3.32e-07	4.0	5.68e-14
50000	54.5	1.36e-07	4.0	7.11e-14
SHBVM
$n$	time	$k$	$s$	$e_u$	$e_H$
400	7.2	20	18	1.31e-11	4.26e-14
600	7.4	16	14	3.70e-12	4.26e-14
800	8.6	14	12	4.75e-12	4.26e-14

Equations

$$\dot{y}=J\nabla H(y)=:f(y),y(0)=y_0\in {\mathbb{R}}^{2m},$$

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

$$H(y(t))-H(y_0)={\int }_0^t\nabla H(y(\tau ){)}^{\top }\dot{y}(\tau )\mathrm{d}\tau ={\int }_0^t\nabla H(y(\tau ){)}^{\top }J\nabla H(y(\tau ))\mathrm{d}\tau =0,$$

$$\sigma (0)=y_0,\sigma (h)=:y_1,y_1\approx y(h),$$

$$H(y_1)-H(y_0)$$

$$={\int }_0^h\nabla H(\sigma (t){)}^{\top }\dot{\sigma }(t)\mathrm{d}t=h{\int }_0^1\nabla H(\sigma (ch){)}^{\top }\dot{\sigma }(ch)\mathrm{d}c=0,$$

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

$$\dot{\sigma }(ch)=\sum _{j=0}^{s-1}{P}_j(c){\gamma }_j(\sigma ),c\in [0,1],$$

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

$${\int }_0^1\nabla H(\sigma (ch){)}^{\top }\dot{\sigma }(ch)\mathrm{d}c$$

$${\gamma }_j(\sigma )=J{\int }_0^1{P}_j(c)\nabla H(\sigma (ch))\mathrm{d}c\equiv {\int }_0^1{P}_j(c)f(\sigma (ch))\mathrm{d}c,$$

$${\gamma }_j(\sigma )={\int }_0^1{P}_j(c)f\left(y_0+h\sum _{i=0}^{s-1}{\int }_0^c{P}_i(\tau )\mathrm{d}\tau {\gamma }_i(\sigma )\right)\mathrm{d}c,j=0,\dots ,s-1.$$

$${\gamma }_j(\sigma )\approx \sum _{\ell =1}^k{b}_{\ell }{P}_j(c_{\ell })f(\sigma (c_{\ell }h))=:{\hat{\gamma }}_j,j=0,\dots ,s-1,$$

$${\hat{\gamma }}_j=\sum _{\ell =1}^k{b}_{\ell }{P}_j(c_{\ell })f\left(y_0+h\sum _{i=0}^{s-1}{\int }_0^{c_{\ell }}{P}_i(\tau )\mathrm{d}\tau {\hat{\gamma }}_i\right),j=0,\dots ,s-1,$$

$$H(y_{1})-H(y_{0})=\left\{\begin{array}[]{cc}0,&\mbox{if}\quad H\in\Pi_{\nu}\quad\mbox{with}\quad\nu\leq 2k/s,\\[5.69054pt] O(h^{2k+1}),&\mbox{otherwise}.\end{array}\right.$$

$$Y_i$$

$$y_1=y_0+h{\hat{\gamma }}_0\equiv y_0+h\sum _{i=1}^k{b}_{i}f(Y_i).$$

$$\begin{array}[]{c|c}\bm{c}&{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\\ \hline\cr&\bm{b}^{\top}\end{array}$$

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

$${\cal P}_{s}=\left(\begin{array}[]{ccc}P_{0}(c_{1})&\dots&P_{s-1}(c_{1})\\ \vdots&&\vdots\\ P_{0}(c_{k})&\dots&P_{s-1}(c_{k})\end{array}\right),\leavevmode\nobreak\ {\cal I}_{s}=\left(\begin{array}[]{ccc}\int_{0}^{c_{1}}P_{0}(x)\mathrm{d}x&\dots&\int_{0}^{c_{1}}P_{s-1}(x)\mathrm{d}x\\ \vdots&&\vdots\\ \int_{0}^{c_{k}}P_{0}(x)\mathrm{d}x&\dots&\int_{0}^{c_{k}}P_{s-1}(x)\mathrm{d}x\end{array}\right)\leavevmode\nobreak\ \in\mathbb{R}^{k\times s}.$$

$$Y=\bm{e}\otimes y_{0}+h{\cal I}_{s}{\cal P}_{s}^{\top}\Omega\otimes I_{2m}f(Y),\quad Y=\left(\begin{array}[]{c}Y_{1}\\ \vdots\\ Y_{k}\end{array}\right),\,f(Y)=\left(\begin{array}[]{c}f(Y_{1})\\ \vdots\\ f(Y_{k})\end{array}\right),\,\bm{e}=\left(\begin{array}[]{c}1\\ \vdots\\ 1\end{array}\right)\in\mathbb{R}^{k},$$

$$F(\hat{\mathbf{\gamma }}):=\hat{\mathbf{\gamma }}-{\mathcal{P}}_s^{\top }\Omega \otimes I_{2m}f\left(\mathbf{e}\otimes y_0+h{\mathcal{I}}_s\otimes I_{2m}\hat{\mathbf{\gamma }}\right)=0,$$

$$\hat{\bm{\gamma}}=\left(\begin{array}[]{c}\hat{\gamma}_{0}\\ \vdots\\ \hat{\gamma}_{s-1}\end{array}\right),\qquad\hat{\gamma}_{i}\in\mathbb{R}^{2m},\quad i=0,\dots,s-1.$$

$$\ddot{q}=\nabla U(q),q(0)=q_0,\text{leavevmode}\dot{q}(0)=p_0\in {\mathbb{R}}^m,$$

$$Q=\mathbf{e}\otimes q_0+h{\mathcal{I}}_s{\mathcal{P}}_s^{\top }\Omega \otimes I_{m}P,P=\mathbf{e}\otimes p_0+h{\mathcal{I}}_s{\mathcal{P}}_s^{\top }\Omega \otimes I_m\nabla U(Q),$$

$$Q=\left(\begin{array}[]{c}Q_{1}\\ \vdots\\ Q_{k}\end{array}\right),\qquad P=\left(\begin{array}[]{c}P_{1}\\ \vdots\\ P_{k}\end{array}\right),\qquad\nabla U(Q)=\left(\begin{array}[]{c}\nabla U(Q_{1})\\ \vdots\\ \nabla U(Q_{k})\end{array}\right).$$

$${\cal P}_{s}^{\top}\Omega{\cal I}_{s}=X_{s}\equiv\left(\begin{array}[]{cccc}\xi_{0}&-\xi_{1}\\ \xi_{1}&0&\ddots\\ &\ddots&\ddots&-\xi_{s-1}\\ &&\xi_{s-1}&0\end{array}\right),\quad\xi_{i}=\left(2\sqrt{|4i^{2}-1|}\right)^{-1},\quad i=0,\dots,s-1,$$

$$Q=\mathbf{e}\otimes q_0+h\mathbf{c}\otimes p_0+h^2{\mathcal{I}}_s{X}_s{\mathcal{P}}_s^{\top }\Omega \otimes I_m\nabla U(Q).$$

$$\bar{\bm{\gamma}}\equiv\left(\begin{array}[]{c}\bar{\gamma}_{0}\\ \vdots\\ \bar{\gamma}_{s-1}\end{array}\right)={\cal P}_{s}^{\top}\Omega\otimes I_{m}\nabla U(Q),\qquad\bar{\gamma}_{i}\in\mathbb{R}^{m},\quad i=0,\dots,s-1,$$

$$G(\bar{\mathbf{\gamma }}):=\bar{\mathbf{\gamma }}-{\mathcal{P}}_s^{\top }\Omega \otimes I_m\nabla U\left(\mathbf{e}\otimes q_0+h\mathbf{c}\otimes p_0+h^2{\mathcal{I}}_s{X}_s\otimes I_m\bar{\mathbf{\gamma }}\right)=0.$$

$$q_1=q_0+h{p}_0+h^2\left({\xi }_0{\bar{\gamma }}_0-{\xi }_1{\bar{\gamma }}_1\right),p_1=p_0+h{\bar{\gamma }}_0,$$

Full text

Keywords: Hamiltonian problems; energy-conserving methods; Hamiltonian Boundary Value Methods; HBVMs; line integral methods; spectral methods; Hamiltonian PDEs; semilinear wave equation; nonlinear Schrödinger equation; Korteweg–de Vries equation.

MSC: 65P10, 65M70, 65M20, 65L05, 65L06

1 Introduction

The numerical solution of ordinary differential equations (ODE) problems, though researched for over sixty years, is still a very active field of investigation, following a number of trends, such as:

1. a)

the search for methods suited for specific relevant classes of problems; 2. b)

their efficient implementation on a computer; 3. c)

the extension of existing methods to cope with wider classes of problems.

Point a) is particularly interesting, since it is nowadays well understood that relevant classes of problems do possess specific geometric properties in their solutions and, often, one is interested in reproducing such properties in the discrete solution obtained by a numerical method. As matter of fact, the term Geometric Integration has been coined to denote the study of numerical methods able to preserve such properties. These latter methods, in turn, are named geometric integrators. As an example, when dealing with dissipative problems, $A$-stable methods are geometric integrators, since they retain the asymptotic stability of equilibria. Nevertheless, when stability results by first approximation do not apply, things become much more involved. This is the case, for example, of Hamiltonian problems, i.e., problems in the form

[TABLE]

with $J=-J^{\top}$ and $H$ a scalar function (which we shall hereafter assume to be suitably regular), called the Hamiltonian or energy. Due to the skew-symmetry of $J$, this latter function turns out to be conserved along the solution of (1). In fact, one has:

[TABLE]

Hamiltonian problems are very important in the applications and, for this reason, their numerical simulation has been the subject of many researches: we refer the reader, e.g., to the monographs [67, 62, 53, 19, 6] and references therein. In particular, numerical methods able to conserve $H$ are geometric integrators, referred to as energy-conserving methods.

Point b) is also paramount: in fact, no numerical method can be really useful, if it cannot be efficiently implemented on a computer. Therefore, a particular care has to be devoted to devise robust implementation techniques, in order to make the studied methods suitable for solving a wide class of problems. In particular, the availability of efficient Newton-type procedures for solving the discrete problems generated by the methods turns out to be central, when numerically solving the Hamiltonian problems described at the next point.

At last, point c) is one of the main focuses of the present paper. In fact, according to [67, page 157], one effective way of solving Hamiltonian PDEs is to discretize, at first, the space variable(s). In so doing, under appropriate space discretizations, one obtains a large-size Hamiltonian problem, which can be then solved by using a suitable geometric integrator. In particular, for sake of simplicity and brevity, in this paper we shall deal with initial-boundary value problems in one space dimension, equipped with periodic boundary conditions, even though the arguments can be extended to cope with higher space dimensions, as is sketched in Section 3.3. As was anticipated above, the numerical solution of the Hamiltonian problems arising from the space discretization of Hamiltonian PDEs will require the use of effective Newton-type procedures, in order to avoid severe step-size limitations.

With these premises, the present paper is devoted to report about recent findings in the numerical solution of Hamiltonian PDEs by using Hamiltonian Boundary Value Methods (HBVMs), a class of energy-conserving Runge-Kutta methods for Hamiltonian problems. The novelty in their use stems from the fact that they provide effective and arbitrarily high-order energy-conserving methods for the time integration of the Hamiltonian semi-discrete problems obtained from Hamiltonian PDEs. In fact, low order methods have been mainly considered for this purpose, so far (see, e.g., [41, 61, 50, 64, 51, 59]). Further approaches can be found in [47, 45, 48, 63, 44, 43, 49, 52, 46]. In more details, the structure of this paper is as follows:

• •

in Section 2 we recall the main facts about HBVMs, also sketching their efficient blended implementation;

• •

in Section 3 we describe the space discretization of the semilinear wave equation, and the efficient solution of the resulting Hamiltonian ODE problem via HBVMs. For this equation we shall provide full details, whereas the whole procedure will be only sketched for the subsequent equations;

• •

in Section 4 we see that the same approach can be used for the nonlinear Schrödinger equation;

• •

in Section 5 we consider, instead, the Korteweg–de Vries equation;

• •

Section 6 contains some numerical tests, aimed at showing the effectiveness of the proposed approach;

• •

at last, a few conclusions are given in Section 7.

2 Hamiltonian Boundary Value Methods (HBVMs)

HBVMs are energy-conserving methods derived within the framework of (discrete) line integral methods, initially proposed in [56, 57, 54, 55, 58], and later refined in [23, 24, 22, 29, 30, 31]. The approach has also been extended along several directions [24, 25, 10, 18, 27, 28, 39, 14], including Hamiltonian BVPs [1], constrained Hamiltonian problems [15], highly-oscillatory problems [38, 21, 2], and Hamiltonian PDEs [13, 3, 17, 21, 40, 16]. We also refer to the review paper [20] and to the monograph [19].

The basic idea line integral methods rely on is that the conservation of an invariant can be recast as the vanishing of a corresponding line-integral. In the case of the Hamiltonian $H$ for (1), one has:

[TABLE]

due to the fact that the integrand is identically zero. Consequently, $H(y(t))=H(y_{0})$, for all $t\geq 0$. Nevertheless, when dealing with a discrete time dynamics, ruled by a time-step $h>0$, one can consider a path $\sigma:[0,h]\rightarrow\mathbb{R}^{2m}$ such that

[TABLE]

and

[TABLE]

but without requiring the integrand to be identically zero. In such a case, there are infinitely many paths satisfying (2)–(3), each providing a corresponding line integral method. In particular, we here consider a polynomial path, which we expand along the orthonormal Legendre basis:

[TABLE]

where, as is usual, $\Pi_{i}$ is the set of polynomials of degree $i$ and $\delta_{ij}$ is the Kronecker symbol. In order to obtain a path $\sigma\in\Pi_{s}$ satisfying (2)-(3), let us then consider the expansion

[TABLE]

in terms of the $s$ unknown vector coefficients $\{\gamma_{j}(\sigma)\}$. In order to fulfill (2), integrating both sides of (5) and taking into account that (see (4)) $\int_{0}^{1}P_{j}(c)\mathrm{d}c=\delta_{j0}$, one obtains:

[TABLE]

Taking into account (5), condition (3) becomes

[TABLE]

which is satisfied by choosing (see (1)):

[TABLE]

because of the skew-symmetry of matrix $J$. Therefore, this specific energy-conserving line integral method is defined by the polynomial path $\sigma$, whose coefficients satisfy the following set of $s$ nonlinear vector equations, derived from (6) and (7):

[TABLE]

Moreover, it can be proved that $\sigma(h)-y(h)=O(h^{2s+1})$, i.e., the approximation procedure has order $2s$ [30, Theorem 1] (see also [20]). However, this procedure does not yet provide a numerical method since, quoting e.g. Dahlquist and Björk [42, page 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.” In particular, since we are dealing with a polynomial approximation, we consider the Gaussian interpolatory quadrature rule, based at the zeros $0<c_{1}<\dots<c_{k}<1$ of $P_{k}$, whose weights we denote, respectively, by $b_{1},\dots,b_{k}$, which is well-known to have order $2k$.444 I.e., it is exact for polynomial integrands up to order $2k-1$. Consequently, with reference to (7), we obtain the approximation

[TABLE]

where, for sake of brevity, we continue to denote $\sigma$ the polynomial approximation. The new discrete problem is then given by

[TABLE]

which remarkably has, alike (8), dimension $s$, independently of $k$.

Definition 1**.**

The discrete problem (10) defines a HBVM$(k,s)$ method. The limit as $k\rightarrow\infty$, given by (8), defines a HBVM$(\infty,s)$ formula.

It is possible to prove the following result [19, 20].

Theorem 1**.**

For all $k\geq s$, by using the $k$ Gauss-Legendre abscissae, a HBVM$(k,s)$ method is symmetric and of order $2s$. Moreover, it reduces to the $s$-stage Gauss collocation method, when $k=s$. Concerning energy-conservation when applied for solving (1), one has:

[TABLE]

It is worth mentioning that, because of (11), by choosing $k$ large enough one can either obtain:

• •

an exact conservation of energy, when $H$ is a polynomial;

• •

a practical conservation of energy, otherwise. In fact, in such a case, it is enough that the energy error falls within the round-off error level.

2.1 Runge-Kutta form of HBVM$(k,s)$

It is possible to see that, actually, a HBVM$(k,s)$ method is a $k$-stage Runge-Kutta method. In fact, by setting in (9) $Y_{\ell}:=\sigma(c_{\ell}h)$, $\ell=1,\dots,k$, one obtains:

[TABLE]

with the new approximation given by

[TABLE]

It can be readily seen that (12)-(13) define the $k$-stage Runge-Kutta method with Butcher tableau

[TABLE]

with

[TABLE]

and

[TABLE]

For this Runge-Kutta method, the stage equation (12) has (block) dimension $k$ and is given by

[TABLE]

having set, in general, $I_{r}\in\mathbb{R}^{r\times r}$ the identity matrix. Nonetheless, the equivalent discrete problem (10), whose dimension is $s$ independently of $k$, turns out to be given by:

[TABLE]

where

[TABLE]

Once (17) is solved, according to (13) the new approximation is given by  $y_{1}=y_{0}+h\hat{\gamma}_{0}$.

2.2 Special second-order problems

Sometimes, the problem (1) assumes the form of a special second-order problem,

[TABLE]

for which, setting $p=\dot{q}$,  y=\left(\begin{array}[]{c}q\\ p\end{array}\right)  and  $H(y)\equiv H(q,p)=\frac{1}{2}p^{\top}p-U(q)$. In such a case, the dimension of the blocks of the discrete problem can be halved. In fact, by using (14) for solving (19), one sees that the stage equations for $q$ and $p$ are respectively given by:

[TABLE]

having set

[TABLE]

Plugging the second equation in (20) into the first one, considering that ${\cal I}_{s}{\cal P}_{s}^{\top}\Omega\bm{e}=\bm{c}$ and, moreover,

[TABLE]

one then obtains:

[TABLE]

Setting (compare with (18))

[TABLE]

and taking into account (22), one then obtains the new discrete problem (compare with (17)):

[TABLE]

Once it has been solved, it can be seen that the new approximations are given by (see, e.g., [19, Chapter 4]):

[TABLE]

where $\xi_{0}$ and $-\xi_{1}$ are the nonzero entries on the first row of matrix $X_{s}$ defined in (21).

2.3 Blended iteration

The efficient solution of the discrete problem (17) has been studied in a series of papers [26, 11, 12, 19]. We here recall the main facts about the so called blended implementation of HBVMs, which represents a Newton-type iteration for solving (17). This approach, at first sketched in [9], has then been analyzed in [32] and developed in [34, 35, 36]. It has been then implemented in the Fortran codes BiM [33] and BiMD [37], for the numerical solution of stiff ODE-IVPs and linearly implicit DAEs: both codes can be retrieved at [72]; the latter code is also available at the Test Set for IVP Solvers [71]. The blended implementation of HBVMs has then been considered in [26] and implemented in the Matlab function hbvm available at the url [73]. We also mention that, more recently, this approach has been also considered for RKN methods [70].

Let us then consider the simplified Newton iteration for solving (17) which, by taking into account (21) amounts to solving the following set of linear systems:

[TABLE]

with $f^{\prime}(y_{0})$ the Jacobian of $f$ evaluated at $y_{0}$. This iteration, though straightforward and very effective, requires, however, the factorization of a $2ms\times 2ms$ matrix, which can be cumbersome, when $s$ and/or $m$ are large. To get rid of this problem, by considering that matrix $X_{s}$ is nonsingular one at first considers the following equivalent formulation of (25), having set $\rho_{s}$ a positive, and for the moment unspecified, parameter:

[TABLE]

The next step is to consider the blending of the two equivalent formulations (25) and (26) with weights $\theta_{s}$ and $I_{s}\otimes I_{2m}-\theta_{s}$, respectively, where:

[TABLE]

In so doing, one obtains a new linear system, whose coefficient matrix has the inverse which can be approximated by $\theta_{s}$. Skipping the details (for which we refer to [26], see also [19, 20]), one then obtains the following blended iteration for solving (17):

[TABLE]

which only requires to factor the matrix $\Sigma$ in (27), having the same size as that of the continuous problem. Concerning the choice of the parameter $\rho_{s}$, as is shown in [32], the optimal choice, based on a linear convergence analysis, turns out to be:

[TABLE]

where, as is usual, $\sigma(X_{s})$ is the spectrum of $X_{s}$.

In the case of the special second-order problem (19), the simplified Newton iteration for solving (24) becomes:

[TABLE]

with $\nabla^{2}U(q_{0})$ the Hessian of $U$ evaluated at $q_{0}$. Consequently, similar steps as above can be repeated, via the following formal substitutions:

[TABLE]

As a result, the blended iteration for solving (24) is given by:

[TABLE]

with the parameter $\rho_{s}$ still given by (29) and

[TABLE]

Consequently, also in such a case, one has only to factor a matrix having the same size as that of the continuous problem.

2.4 Blended iteration for semilinear problems

Once more, we stress that the availability of a Newton-type iteration for solving (17) is paramount, in order to avoid severe step-size limitations, when such a problem is derived from the space discretization of Hamiltonian PDEs. In fact, in such a case, the resulting ODE problem turns out to be in the form

[TABLE]

with the dimension and the norm of matrix $A$ tending to infinity, as the space discretization is made more and more accurate, whereas $\|g\|$ remains bounded, if the solution is bounded. Consequently, one can consider a constant approximation of the Jacobian of the right-hand side of (33), given by the matrix $A$ of the linear term. As a result, the matrix $\Sigma$ defined in (27) becomes

[TABLE]

which is constant for all time-steps and, consequently, it needs to be factored only once.

Similarly, when problem (19) is in the form

[TABLE]

with $A^{2}$ symmetric and semi-positive definite, and $\|A^{2}\|\gg\|g\|$, one can approximate the matrix $\Sigma$ in (32) as

[TABLE]

which, also in this case, is constant for all time-steps and needs to be factored only once.

We end this subsection by stressing that, for the problems that we shall consider in the sequel, matrix $A$ in (34), or matrix $A^{2}$ in (36), has a block structure with diagonal blocks. As a result, the corresponding blended iterations (28) and (31) are computationally inexpensive. Moreover, the linear algebra can be made still more efficient, as is done in the Matlab function hbvm available at [73], by considering a matrix formulation of the iteration [19, 68].

2.5 HBVMs as spectral methods in time

To conclude this quick introduction to HBVMs, we mention their use as spectral methods in time, which has been the subject of recent investigations [38, 21, 2]. We mention that the use of Runge-Kutta methods as spectral methods in time has been considered previously in [7, 4, 5, 69] (see also [30]). In more details, if we consider the expansion of the right-hand side of (1), on the interval $[0,h]$, along the Legendre basis (4), one has:

[TABLE]

where $\gamma_{j}(y)$ is defined according to (7), by formally replacing $\sigma$ by $y$. On the other hand, the polynomial approximation $\sigma$ defined in (5) is obtained by truncating the previous series after $s$ terms. However, by considering that

[TABLE]

one has that

[TABLE]

the more regular $f(y)$, the faster the convergence to 0 of $\|\gamma_{j}(y)\|_{2}$, as $j\rightarrow\infty$. Consequently, when using a finite precision arithmetic with machine epsilon $\varepsilon$, if one truncates the expansion (37) when the Fourier coefficient $\gamma_{s}(y)$ is negligible, w.r.t. the previous ones, then one obtains that (37) and (5) become indistinguishable, in the used finite precision arithmetic. A straightforward criterion for this to happen, considered in [38, 21], is to require that

[TABLE]

with $tol\approx\varepsilon$. Moreover, the analysis in [2] shows that one could even use $tol\approx\sqrt{\varepsilon}$ in (38), still obtaining full machine accuracy at $t=h$. At last (see (9)), by choosing $k$ large enough, one may obtain full machine accuracy in the approximation of $\gamma_{j}(\sigma)$ by means of $\hat{\gamma}_{j}$, $j=0,\dots,s-1$. As a result, the use of HBVMs as spectral methods in time (which we shall denote by SHBVMs, as an abbreviation for spectral HBVMs) usually requires the use of relatively large values of $s$ and $k$. This, in turn, is not a big issue; in fact:

• •

on one hand, we have the availability of the blended iteration (28) (or (31)), whose computational cost is mildly affected by such parameters, also considering the approximation (34) (or (36));

• •

on the other hand, SHBVMs will allow the use of relatively large time-steps.

Summing all up, overall SHBVMs will result to be extremely effective and competitive, as is testified by the numerical tests reported in Section 6 (see also [2, 38, 21]).

3 The semilinear wave equation

The first Hamiltonian PDE that we consider is the semilinear wave equation:

[TABLE]

with $f^{\prime}$ the derivative of $f$. The problem (39) is completed by prescribing periodic boundary conditions. Hereafter, we shall assume the solution to be suitably regular, as a periodic function in space. Moreover, for sake of brevity, we shall omit the arguments of the involved functions, when not necessary. By setting $v=u_{t}$, one obtains that (39) is a Hamiltonian PDE, with Hamiltonian functional

[TABLE]

so that, by setting

[TABLE]

the vector of the functional derivatives of ${\cal H}$, with

[TABLE]

one has:

[TABLE]

which is formally in the form (1). As in the ODE case, also now one has the conservation of the Hamiltonian.

Theorem 2**.**

Assuming that the solution of (39) is suitably smooth in space, the Hamitonian (40) is conserved, when periodic boundary conditions are prescribed.

Proof.

In fact, from (39) and (40), and taking into account that $v=u_{t}$, one has:

[TABLE]

because of the periodic boundary conditions. ∎

In order to numerically solve (39), according to what sketched in the introduction, we at first discretize the space variable, with the aim of obtaining a corresponding Hamiltonian ODE problem. For this purpose, we consider the following orthonormal basis on the interval $[a,b]$, which takes into account of the periodic boundary conditions [13, 19, 3, 16, 17, 20]:

[TABLE]

In fact, for all allowed $i,j$, one has:

[TABLE]

Consequently, for suitable time dependent coefficients $\alpha_{0}(t),\alpha_{1}(t),\beta_{1}(t),\dots$, one has:

[TABLE]

The infinite expansion (46) can be cast in vector form, by defining the infinite-dimensional vectors

[TABLE]

as

[TABLE]

By also introducing the infinite matrix

[TABLE]

and considering that (45) can be written in matrix form as

[TABLE]

the identity operator, we then prove the following result.

Theorem 3**.**

With reference to (42)–(50), problem (39) can be rewritten as the special second-order problem

[TABLE]

By setting $\bm{p}=\dot{\bm{q}}$, this latter problem is Hamiltonian, with Hamiltonian function

[TABLE]

which turns out to be equivalent to the Hamiltonian functional (40).

Proof.

Problem (51) is clearly Hamiltonian, w.r.t. the Hamiltonian (52), since

[TABLE]

Let us then show that:

• •

(51) is equivalent to (39);

• •

(52) is equivalent to (40).

Concerning the first point, we observe that

[TABLE]

with an obvious meaning of $\bm{\omega}"(x)$, so that (39) can be rewritten as

[TABLE]

Multiplying both sides by $\bm{\omega}(x)$, then integrating in space from $a$ to $b$, and taking into account (50), give us (51).

Concerning the second point, the statement easily follows by considering that $v(x,t)=u_{t}(x,t)=\bm{\omega}(x)^{\top}\bm{p}$, and

[TABLE]

Moreover, by defining the matrix

[TABLE]

where matrix $J_{2}$ is that defined in (41), one has:

[TABLE]

so that, by taking again into account (50), one obtains:

[TABLE]

The proof is completed by considering that, from (48), $f(\bm{\omega}(x)^{\top}\bm{q}(t))=f(u(x,t))$. ∎

3.1 Discretization

In order to solve problem (51) on a computer, the infinite expansion (46) must be truncated at a convenient index $N$. In so doing, (47) and (49) respectively become

[TABLE]

so that (48) continues formally to hold, even though now $u$ is no more the solution of (39).555Observe that, for sake of brevity, we continue to use the same notation $\bm{\omega},\bm{q}$, and $D$ used for the infinite expansion, though, hereafter, they will refer to the finite counterparts (54). Nevertheless, in the spirit of Fourier-Galerkin methods [8], by requiring the residual obtained by plugging $u$ into (39) be orthogonal to the functional subspace666The same Fourier-Galerkin procedure will be used for the Hamiltonian PDEs studied in the following sections.

[TABLE]

which, for fixed $t$, contains $u$, one obtains a finite set of $2N+1$ ODEs, formally still given by (51), for which Theorem 3 continues formally to hold, with the only exception that now $H(\bm{q},\bm{p})$ is no more equivalent to the Hamiltonian functional (40), but only yields an approximation to it. Nevertheless, it is well known that, under suitable regularity assumptions on $f$ and the initial data $u_{0}$ and $v_{0}$, this truncated version converges exponentially to the original functional (40), as $N\rightarrow\infty$ (this phenomenon is usually referred to as spectral accuracy), as well as the truncated version of $u$ converges to the infinite expansion (46).

The resulting finite-dimensional semi-discrete problem (51), which is still Hamiltonian, is the one we will solve by using line-integral methods. Actually, it is not yet ready to be solved, since the integral $\int_{a}^{b}\bm{\omega}(x)f^{\prime}(\bm{\omega}(x)^{\top}\bm{q})\mathrm{d}x$, appearing in it, needs to be evaluated. For this purpose, since we are dealing with an integrand which is periodic in space, a composite trapezoidal rule based at the abscissae

[TABLE]

can be considered. We refer, e.g., to [13, Theorems 5, 6], for a proper choice of the number, $m+1$, of points in (55), able to preserve the property of spectral accuracy.

Problem (51) can then be solved by using a HBVM$(k,s)$ method, for which the accuracy results of Theorem 1 hold true. In particular, concerning the conservation of the semi-discrete Hamiltonian (52), the next result holds true, which follows from (11).

Theorem 4**.**

If a HBVM$(k,s)$ method is used with time-step $h$ for solving (51), one has

[TABLE]

having set $\bm{q}_{1}\approx\bm{q}(h)$ and $\bm{p}_{1}\approx\bm{p}(h)$ the new approximations.

3.2 The nonlinear iteration

In light of what previously stated, in order to obtain a spectral accuracy in space a suitably large value of $N$ in (54) has to be considered (a practical criterion for its choice will be sketched in Section 6). Consequently, the special second-order problem (51) is semilinear, with a bounded nonlinear term,777When the solution is bounded. and the linear term given by $-D^{2}\bm{q}$. On the other hand, both the size (i.e., $2N+1$) and the norm (i.e., $\left(\frac{2\pi N}{b-a}\right)^{2}$) of the matrix $D^{2}$ tend to infinity, as $N\rightarrow\infty$. Consequently, when using a HBVM$(k,s)$ method for solving (51), the blended iteration (29), (31)-(32) can be conveniently used, to get rid of the large norm of the linear term, with matrix $\Sigma$ approximated as in (36). As a result, it turns out to be given by

[TABLE]

which is a diagonal matrix and, therefore, $\Sigma^{-1}$ can be cheaply computed and stored. Consequently, the complexity of the blended iteration turns out to be comparable with that of an explicit method, though not suffering from the step-size restrictions of this latter. As a matter of fact, the use of an explicit method usually would require $h\|D\|<1$, i.e., $h=O(N^{-1})$, which may be restrictive, when $N\gg 1$.

3.3 Extension to higher space dimensions

For completeness, in this section we sketch the generalization of most of the previous arguments to the case where the space domain of the wave equation is, for sake of simplicity, the square $[a,b]^{2}:=[a,b]\times[a,b]$:

[TABLE]

As before, the problem (58) is completed by prescribing periodic boundary conditions. In such a case, the Hamiltonian functional becomes, by setting as usual  $v=u_{t}$,

[TABLE]

Because of the periodic boundary conditions, we can again consider the orthonormal basis (42)–(44) in each space dimension, thus obtaining the expansion (for sake of brevity, let us set $s_{0}\equiv 0$)

[TABLE]

involving the additional time-dependent coefficients $\eta_{0}(t),\mu_{1}(t),\eta_{1}(t),\dots$. With reference to the infinite-dimensional vectors in (47), and defining the vectors

[TABLE]

the infinite expansion (60) can be cast in vector form as

[TABLE]

Consequently, by taking into account (49)-(50), one obtains that (compare with Theorem 3) problem (58) can be recast as the infinite set of second order ODEs:

[TABLE]

By setting $\bm{p}_{1}\otimes\bm{p}_{2}=\dot{\bm{q}}_{1}\otimes\dot{\bm{q}}_{2}$, one then obtains that problem (3.3) is Hamiltonian, with Hamiltonian function

[TABLE]

This latter function, in turn, is equivalent to the Hamiltonian functional (59), via the expansion (61). Then, as done in the one dimensional case, the vectors $\bm{\omega}(x),\bm{\omega}(y),\bm{q}_{1}(t),\bm{q}_{2}(t)$, are truncated after $2N+1$ terms, for a convenient large value of $N$, so that (3.3) becomes a Hamiltonian set of $(2N+1)^{2}$ ODEs, with Hamiltonian (63). This problem can be solved by adapting the arguments previously explained in the one dimensional case, even though now the complexity is clearly increased. Remarkably enough, however, the diagonal structure of the Jacobian of the linear term in (3.3), i.e., $-(D\otimes D)^{2}$, is still preserved.888Evidently, this property holds true whichever is the dimension of the considered space domain.

4 The nonlinear Schrödinger equation

We now consider the nonlinear Schrödinger equation, which is very important in many applications (see, e.g., the introduction in [3]). In real variables, it takes the form,

[TABLE]

$f^{\prime}$ being the derivative of a suitably regular function $f$. The problem is completed with periodic boundary conditions and, hereafter, we shall assume the initial functions to be suitably regular (as periodic functions), in order to guarantee a suitably smooth solution. Such an equation can be written in the form (41), with $\nabla{\cal H}$ the vector of the functional derivatives of the Hamiltonian functional

[TABLE]

This latter functional is conserved, because of the periodic boundary conditions [3, Theorem 1]. Additional conserved (quadratic) functionals are the mass and the momentum [3, Theorem 2], respectively given by:

[TABLE]

In order to obtain a space discretization which takes into account of the periodic boundary conditions, we consider again the expansion along the Fourier basis (42)–(45), for $u$ and $v$. The expansion for $u$ is formally still given by (46). Similarly, that for $v$ will be given by:

[TABLE]

for suitable time dependent coefficients, $\eta_{0}(t),\eta_{1}(t),\mu_{1}(t),\dots$. By using the infinite vectors (47) and

[TABLE]

we can cast the expansions of $u$ and $v$ in vector form, respectively, as (48) and

[TABLE]

As a consequence, the following result holds true, whose proof is similar to that of Theorem 2 (see also [3, Section 2]).

Theorem 5**.**

With reference to (42)–(50), problem (64) can be rewritten as the infinite-dimensional Hamiltonian ODE problem

[TABLE]

This latter problem is Hamiltonian w.r.t. the Hamiltonian

[TABLE]

which turns out to be equivalent to the Hamiltonian functional (65). Moreover, the two quadratic invariants (66) can be respectively rewritten as

[TABLE]

where $\bar{D}$ is the matrix defined in (53).

As in the case of the nonlinear wave equation, in order to solve problem (70) on a computer, one needs to truncate the infinite expansions (46) and (67) at a convenient index $N$. In so doing, the infinite vectors and matrices (47), (49), and (68) become those in (54) and

[TABLE]

respectively. As a result, one eventually arrives again at the finite-dimensional Hamiltonian ODE problem (70), having dimension $4N+2$, with the Hamiltonian and the invariants still given by (71) and (72), respectively. Again, spectral accuracy is expected, if the solution is regular enough in space (as a periodic function). Finally, we mention that also in this case the integrals in space can be computed by means of a composite trapezoidal rule, based at the abscissae (55), for a suitably large value of $m$.

Again, we can use an HBVM$(k,s)$ method for solving (70). Concerning the conservation of the Hamiltonian, the following straightforward result follows from (11).

Theorem 6**.**

If a HBVM$(k,s)$ method is used with time-step $h$ for solving (70), one has

[TABLE]

having set $\bm{q}_{1}\approx\bm{q}(h)$ and $\bm{p}_{1}\approx\bm{p}(h)$ the new approximations.

4.1 The nonlinear iteration

Following the same arguments discussed in the previous section, in order to obtain a spectral accuracy in space a suitably large value of $N$ in (54) and (73) has to be considered (we remind that a practical criterion for its choice will be given in Section 6). Consequently, the Hamiltonian problem (70) is semilinear, with a bounded nonlinear term, if the solution is bounded, and the linear term given by (see (41))

[TABLE]

On the other hand, the norm of the matrix $D^{2}$ is $\left(\frac{2\pi N}{b-a}\right)^{2}$, and tends to infinity, as $N\rightarrow\infty$, as well as its size. Consequently, when using a HBVM$(k,s)$ method for solving (70), the blended iteration (27)–(29) can be conveniently used, to get rid of the large norm of the linear term, with matrix $\Sigma$ approximated as in (34) and, in the present context, given by

[TABLE]

which is a block matrix with diagonal blocks and, therefore, $\Sigma^{-1}$ can be cheaply computed and stored. As matter of fact, one has [3, Theorem 5]:

[TABLE]

which is again a block matrix with diagonal blocks (actually, two vectors are enough to store it). As a consequence, also in the present case the complexity of the blended iteration turns out to be comparable with that of an explicit method, though not suffering from step-size restrictions. As a matter of fact, the use of an explicit method would require $h\|D^{2}\|<1$, i.e., $h=O(N^{-2})$, which may be very restrictive, when $N\gg 1$.

5 The Korteweg–de Vries (KdV) equation

The last Hamiltonian PDE that we consider is the Korteweg–de Vries equation, recently investigated in [16] by using line integral methods,

[TABLE]

with $\alpha\beta\neq 0$, and coupled with periodic boundary conditions. As usual, we shall assume that $u_{0}$ is smooth enough, as a periodic function, so that $u(x,t)$ turns out to be suitably regular, as a periodic function in space [60]. Equation (76) can be written in Hamiltonian form as

[TABLE]

with ${\cal H}[u]$ the Hamiltonian functional

[TABLE]

and

[TABLE]

its functional derivative.999Actually, it can be seen that there is a further Hamiltonian formulation of (76) [65], so that the PDE has a so called bi-Hamiltonian structure. Because of the periodic boundary conditions, the Hamiltonian functional (77) turns out to be conserved. Another conserved functional is given by

[TABLE]

as it can be readily shown. In order to obtain a space discretization, we consider an expansion along the usual orthonormal basis (42)–(45), which provides us with an expression formally still given by (46). However, because of the conservation of the functional (78), one obtains that

[TABLE]

Consequently, the expansion (46) now becomes

[TABLE]

In order to put this expansion in vector form, let us introduce the infinite vectors

[TABLE]

In so doing, we can rewrite (79) as:

[TABLE]

Consequently, the conservation of (78) is automatically granted. Moreover, by defining the infinite matrix

[TABLE]

such that

[TABLE]

and similarly for the higher derivatives, and considering that

[TABLE]

one verifies that (76) can be rewritten as the infinite dimensional ODE problem (see [16, Lemma 3] for full details)

[TABLE]

For this problem, the following result holds true [16, Theorem 1].

Theorem 7**.**

Problem (85) is in the form (1) with (see (41))

[TABLE]

and the Hamiltonian given by

[TABLE]

This latter is equivalent to the Hamiltonian functional (77), via (81) and (83)-(84).

As done before, in order for the problem (85) to be solvable on a computer, the infinite expansion in (79) must be truncated to a convenient index $N$. In so doing, one still formally retrieves the vector formulation (81), where now the vectors

[TABLE]

are hereafter used in place of (80). Similarly, by replacing matrix (82) with

[TABLE]

one obtains a set of $2N$ Hamiltonian equations, formally still given by (85), with the Hamiltonian $H$ also formally given by (86). As for the Hamiltonian PDEs previously studied, spectral accuracy is expected, as $N\rightarrow\infty$, upon regularity assumptions on $u_{0}$.101010Concerning the integrals appearing in (85) and (86), they can be exactly computed via a composite trapezoidal rule based at the abscissae (55), by choosing $m>3N$ [16].

Having got the finite dimensional Hamiltonian ODE problem (85), we can use a HBVM$(k,s)$ method for its time integration. Concerning energy conservation, the following result easily follows from (11).

Theorem 8**.**

A HBVM$(k,s)$ method used for solving (85) is energy-conserving, for all $k\geq 3s/2$.

5.1 The nonlinear iteration

Also in this case, problem (85) is semilinear. However, it is worth observing that the Hessian of the Hamiltonian $H$ in (86) is given by

[TABLE]

with $u(x,t)$ given by the expansion (81). Consequently, by considering the constant approximation $u(x,t)\equiv\hat{u}_{0}$ (due to the conservation of (78)), and taking into account (84), one obtains the constant approximate (diagonal) Hessian

[TABLE]

Therefore, the blended iteration (27)–(29) can be conveniently used, by considering the resulting approximated matrix (see (41) and (88))

[TABLE]

which is a block matrix with diagonal blocks. Moreover, one has [16, Theorem 3]:

[TABLE]

which can be easily computed (once for all) and stored (in fact, only two vectors of length $N$ are needed). Consequently, the complexity of the blended iteration turns out to be comparable with that of an explicit method, though not suffering from its step-size restrictions which, for the present problem, would require $h=O(N^{-3})$.

6 Numerical tests

In this section, we report a few numerical tests, aimed at assessing the effectiveness of HBVMs for solving the previously studied Hamiltonian PDEs. In particular, the spectral version of HBVMs (SHBVMs) will be recognized to be very promising. In more details, we shall compare the following methods:

• •

the symplectic $s$-stage Gauss methods, $s=1,2$;

• •

the energy-conserving HBVM$(k,s)$ methods, $s=1,2$, and $k$ suitably chosen;

• •

the SHBVM method.

The comparisons will be quite fair, since the same Matlab function 111111It is a modification of the function hbvm available at [73]. implements all methods. All numerical tests have been done on a 2.8 GHz Intel Core i7 computer with 16GB of memory, running Matlab 2017b.

To begin with, let us define the criterion used for getting spectral accuracy in space, i.e., for a correct choice of $N$ in (54), (73), (87), and (88). In more details, $N$ has been chosen in order to fulfil both the two following requirements:

• •

a good approximation of the initial condition. This is achieved by requiring

[TABLE]

with $\varepsilon$ the machine epsilon, for problems (51) and (70), or

[TABLE]

for problem (85);

• •

a good approximation of the Hamiltonian. This is achieved by computing the initial value $H(\bm{q}_{0},\bm{p}_{0})=:H_{0}$ of the semi-discrete Hamiltonian (i.e., (52), or (71), or (86)) for consecutive values of $N$, and checking that the absolute value of the difference, $\Delta H_{0}$, satisfies:

[TABLE]

6.1 The semilinear wave equation

We consider the so called sine-Gordon equation [13, Section 7] with a breather soliton solution,

[TABLE]

where we choose $\gamma=1.5$. Its solution, depicted in the upper plot in Figure 1, is:

[TABLE]

In the lower plot of Figure 1 there are the graphs of $E_{0}$ and $\Delta H_{0}$, as defined in (89) and (91), respectively. From such plots, one infers that the choice $N=250$ is adequate to obtain spectral accuracy in space. In Table 1 we list the obtained numerical results by solving the resulting semi-discrete problem (51) with time-step $h=100/n$. In more details: the execution time (in sec), the maximum solution and Hamiltonian errors, $e_{u}$ and $e_{H}$, respectively, and the rate of convergence, where appropriate; for the SHBVM method, we also list the used values of $k$ and $s$, the latter obtained by using $tol\approx\sqrt{\varepsilon}$ in (38) and $k$ suitably larger than $s$. From the obtained results, one sees that:

• •

the higher-order methods perform better than the lower-order ones;

• •

the energy-conserving methods are slightly more efficient than the symplectic ones, when the largest time-steps are used;

• •

the spectral method turns out to be the most effective one, and uses much larger time-steps.

6.2 The nonlinear Schrödinger equation

We consider the so called focusing equation,121212The de-focusing case is obtained when the sign of the coupling term is reversed.

[TABLE]

where the initial conditions at $t=0$ are taken from the known solution,

[TABLE]

depicted in the upper plot of Figure 2, plus (approximate) boundary conditions. In the lower plot of the same figure, there are the plots of $E_{0}$ and $\Delta H_{0}$, as defined in (89) and (91), respectively. From such plots, one infers that the choice $N=600$ is adequate to obtain spectral accuracy in space. For this problem, the symplectic $s$-stage Gauss methods conserve the quadratic invariants (72), whereas the HBVM$(2s,s)$ methods are energy conserving (according to Theorem 6, since $f(x)=x^{2}$). For the SHBVM method, we use $tol\approx 10^{-1}\sqrt{\varepsilon}$ in (38). In Table 2 we list the numerical results obtained by solving the resulting semi-discrete problem (70) with time-step $h=20/n$: besides the execution time (in sec), we list the maximum solution, mass, momentum, and Hamiltonian errors, $e_{uv}$, $e_{1}$, $e_{2}$, and $e_{H}$, respectively, along with the rate of convergence, where appropriate; for the SHBVM method, we also list the used values of $k$ and $s$. We observe that a kind of super-convergence occurs in the invariants (twice the convergence order of the solution) for the Gauss and HBVM methods. In this case, the symplectic and energy-conserving methods turn out to be almost equivalent, with the higher-order methods more efficient than the lower-order ones. However, the SHBVM method outperform all of them, being able to use much larger time-steps, and having a uniformly small error in both the solution and the invariants (which are all conserved within the round-off error level).

6.3 The Korteweg–de Vries equation

This example is adapted from [16, Example 2]:

[TABLE]

equipped with periodic boundary conditions and the initial condition obtained from the known cnoidal wave solution,

[TABLE]

Here $\mathrm{cn}:=\mathrm{cn}(z|m)$ is the Jacobi elliptic function with modulus $m$, $K(m)$ is the complete elliptic integral of the first kind, and the following parameters have been used:

[TABLE]

The initial part of the solution (97) is depicted in the upper plot of Figure 3, whereas in the lower plot one may find $E_{0}$ and $\Delta H_{0}$, as defined in (90) and (91), respectively, versus $N$ . From the latter plots, one infers that the choice $N=50$ is adequate to obtain spectral accuracy in space. By recalling the result of Theorem 8 for HBVMs, in Table 3 we list the numerical results obtained by solving the resulting semi-discrete problem (85) with time-step $h=10/n$, in terms of: execution time (in sec); maximum solution and Hamiltonian errors, $e_{u}$ and $e_{H}$, respectively; rate of convergence, where appropriate.131313Also in this case, for the Gauss method a super-convergence occurs in the Hamiltonian error. For the SHBVM method, we also list the used values of $k$ and $s$, the latter obtained by using $tol\approx 10^{-1}\sqrt{\varepsilon}$ in (38) and $k$ suitably larger than $s$. From the obtained results, one sees that the energy-conserving and symplectic methods are almost equivalent, with the higher-order methods performing better than the lower-order ones. Also in this case, however, the spectral method turns out to be the most effective, being able to use much larger time-steps, with uniformly small solution and Hamiltonian errors.

6.4 A few remarks

From the obtained results, we can draw a few conclusions, which we report in the sequel.

Energy-conservation.

When the conservation of energy is not an issue, the performance of energy-conserving HBVMs seems to be comparable with that of the symplectic Gauss formulae of the same order. Clearly, things may change when energy-conservation is an important feature (see, e.g., the example in [13, Section 7]).

Order of the methods.

From the numerical results, one clearly sees that the second-order methods are outperformed by higher-order HBVMs and/or Gauss methods. In particular, for problems (94) and (96), the second-order HBVM(2,1) method is exactly energy-conserving, and can be regarded as a high-performance implementation of the AVF method in [66]. Despite this, its performance is not comparable with that of the higher-order methods.

Spectral methods in time.

The obtained numerical results further confirm what recently observed in [38, 21, 2], i.e., that the use of HBVMs as spectral methods in time is a very promising way of getting very high-performance ODE solvers, due to the effectiveness of the underlying blended iteration described in Section 2.

7 Conclusions

In this paper we have reviewed the basic facts concerning the use of energy-conserving line integral methods for efficiently solving Hamiltonian PDEs. This has been done by performing, at first, a suitable space discretization, along a Fourier orthonormal basis, thus obtaining a corresponding high-dimensional Hamiltonian problem. In particular, we have studied the semilinear wave equation, the nonlinear Schrödinger equation, and the Korteweg–de Vries equation in one dimension. It is worth mentioning, however, that: as sketched in Section 3.3, the used space discretization can be straightforwardly extended to the case of more space dimensions; additional Hamiltonian PDEs have been considered in [40, 17]. In the future, we plan to further investigate Hamiltonian PDEs within the same framework.