Line Integral solution of Hamiltonian PDEs

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.


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: a) the search for methods suited for specific relevant classes of problems; b) their efficient implementation on a computer; 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 forṁ y = J∇H(y) =: f (y), with J = −J 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: d dt H(y) = ∇H(y) ẏ = ∇H(y) J∇H(y) = 0.
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 [6,19,53,62,67] 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,50,51,59,61,64]). Further approaches can be found in [43-49, 52, 63]. 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.

Hamiltonian Boundary Value Methods (HBVMs)
HBVMs are energy-conserving methods derived within the framework of (discrete) line integral methods, initially proposed in [54][55][56][57][58], and later refined in [22][23][24][29][30][31]. The approach has also been extended along several directions [10,14,18,24,25,27,28,39], including Hamiltonian BVPs [1], constrained Hamiltonian problems [15], highly-oscillatory problems [2,21,38], and Hamiltonian PDEs [3,13,16,17,21,40]. 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: ∇H(y(τ )) J∇H(y(τ ))dτ = 0, due to the fact that the integrand is identically zero. Consequently, H(y(t)) = H(y 0 ), for all t ≥ 0. Nevertheless, when dealing with a discrete time dynamics, ruled by a time-step h > 0, one can consider a path σ : [0, h] → R 2m such that and 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: where, as is usual, Π i is the set of polynomials of degree i and δ ij is the Kronecker symbol.
In order to obtain a path σ ∈ Π s satisfying (2)-(3), let us then consider the expansioṅ in terms of the s unknown vector coefficients {γ j (σ)}. In order to fulfill (2), integrating both sides of (5) and taking into account that (see (4)) 1 0 P j (c)dc = δ j0 , one obtains: Taking into account (5), condition (3) becomes which is satisfied by choosing (see (1)): because of the skew-symmetry of matrix J. Therefore, this specific energy-conserving line integral method is defined by the polynomial path σ, whose coefficients satisfy the following set of s nonlinear vector equations, derived from (6) and (7): Moreover, it can be proved that σ(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 < · · · < c k < 1 of P k , whose weights we denote, respectively, by b 1 , . . . , b k , which is well-known to have order 2k. 1 Consequently, with reference to (7), we obtain the approximation where, for sake of brevity, we continue to denote σ the polynomial approximation. The new discrete problem is then given bŷ which remarkably has, alike (8), dimension s, independently of k.
1 I.e., it is exact for polynomial integrands up to order 2k − 1.
It is possible to prove the following result [19,20].
Theorem 1. For all k ≥ 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: 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.

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 := σ(c h), = 1, . . . , k, one obtains: with the new approximation given by It can be readily seen that (12)-(13) define the k-stage Runge-Kutta method with Butcher tableau with and For this Runge-Kutta method, the stage equation (12) has (block) dimension k and is given by having set, in general, I r ∈ R r×r the identity matrix. Nonetheless, the equivalent discrete problem (10), whose dimension is s independently of k, turns out to be given by: Once (17) is solved, according to (13) the new approximation is given by y 1 = y 0 + hγ 0 .

Special second-order problems
Sometimes, the problem (1) assumes the form of a special second-order problem, for which, setting p =q, y = q p and H(y) ≡ H(q, p) = 1 2 p 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: having set Plugging the second equation in (20) into the first one, considering that I s P s Ωe = c and, moreover, one then obtains: Setting (compare with (18)) and taking into account (22), one then obtains the new discrete problem (compare with (17)): Once it has been solved, it can be seen that the new approximations are given by (see, e.g., [19,Chapter 4]): where ξ 0 and −ξ 1 are the nonzero entries on the first row of matrix X s defined in (21).

Blended iteration
The efficient solution of the discrete problem (17) has been studied in a series of papers [11,12,19,26]. 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: with f (y 0 ) the Jacobian of f evaluated at y 0 . This iteration, though straightforward and very effective, requires, however, the factorization of a 2ms × 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 ρ s a positive, and for the moment unspecified, parameter: The next step is to consider the blending of the two equivalent formulations (25) and (26) with weights θ s and I s ⊗ I 2m − θ s , respectively, where: In so doing, one obtains a new linear system, whose coefficient matrix has the inverse which can be approximated by θ s . Skipping the details (for which we refer to [26], see also [19,20]), one then obtains the following blended iteration for solving (17): which only requires to factor the matrix Σ in (27), having the same size as that of the continuous problem. Concerning the choice of the parameter ρ s , as is shown in [32], the optimal choice, based on a linear convergence analysis, turns out to be: where, as is usual, σ(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: with ∇ 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: As a result, the blended iteration for solving (24) is given by: with the parameter ρ s still given by (29) and Consequently, also in such a case, one has only to factor a matrix having the same size as that of the continuous problem.

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 forṁ 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 Σ defined in (27) becomes which is constant for all time-steps and, consequently, it needs to be factored only once. Similarly, when problem (19) is in the form with A 2 symmetric and semi-positive definite, and A 2 g , one can approximate the matrix Σ in (32) as 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].

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 [2,21,38]. We mention that the use of Runge-Kutta methods as spectral methods in time has been considered previously in [4,5,7,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: where γ j (y) is defined according to (7), by formally replacing σ by y. On the other hand, the polynomial approximation σ defined in (5) is obtained by truncating the previous series after s terms. However, by considering that the more regular f (y), the faster the convergence to 0 of γ j (y) 2 , as j → ∞. Consequently, when using a finite precision arithmetic with machine epsilon ε, if one truncates the expansion (37) when the Fourier coefficient γ 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 [21,38], is to require that with tol ≈ ε. Moreover, the analysis in [2] shows that one could even use tol ≈ √ ε 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 γ j (σ) by means ofγ j , j = 0, . . . , 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,21,38]).

The semilinear wave equation
The first Hamiltonian PDE that we consider is the semilinear wave equation: with f 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 so that, by setting the vector of the functional derivatives of H, with one has: 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: 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 [3,13,16,17,19,20]: In fact, for all allowed i, j, one has: Consequently, for suitable time dependent coefficients α 0 (t), α 1 (t), β 1 (t), . . . , one has: The infinite expansion (46) can be cast in vector form, by defining the infinite-dimensional vectors By also introducing the infinite matrix and considering that (45) can be written in matrix form as 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 By setting p =q, this latter problem is Hamiltonian, with Hamiltonian function which turns out to be equivalent to the Hamiltonian functional (40).
Concerning the first point, we observe that with an obvious meaning of ω"(x), so that (39) can be rewritten as Multiplying both sides by ω(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) = ω(x) p, and Moreover, by defining the matrix where matrix J 2 is that defined in (41), one has: so that, by taking again into account (50), one obtains: The proof is completed by considering that, from (48), f (ω(x) q(t)) = f (u(x, t)).

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 so that (48) continues formally to hold, even though now u is no more the solution of (39). 2 Nevertheless, in the spirit of Fourier-Galerkin methods [8], by requiring the residual obtained by plugging u into (39) be orthogonal to the functional subspace 3 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(q, 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 → ∞ (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 b a ω(x)f (ω(x) q)dx, 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 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).
having set q 1 ≈ q(h) and p 1 ≈ p(h) the new approximations.

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, 4 and the linear term given by −D 2 q. On the other hand, both the size (i.e., 2N + 1) and the norm i.e., 2πN b−a 2 of the matrix D 2 tend to infinity, as N → ∞. Consequently, when using a HBVM(k, s) method for 2 Observe that, for sake of brevity, we continue to use the same notation ω, q, and D used for the infinite expansion, though, hereafter, they will refer to the finite counterparts (54). 3 The same Fourier-Galerkin procedure will be used for the Hamiltonian PDEs studied in the following sections. 4 When the solution is bounded.
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 Σ approximated as in (36). As a result, it turns out to be given by which is a diagonal matrix and, therefore, Σ −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 1.

Extension to higher space dimensions
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 , 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 ≡ 0) involving the additional time-dependent coefficients η 0 (t), µ 1 (t), η 1 (t), . . . . With reference to the infinite-dimensional vectors in (47), and defining the vectors the infinite expansion (60) can be cast in vector form as 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: By setting p 1 ⊗ p 2 =q 1 ⊗q 2 , one then obtains that problem (62) is Hamiltonian, with Hamiltonian function 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 ω(x), ω(y), q 1 (t), q 2 (t), are truncated after 2N + 1 terms, for a convenient large value of N , so that (62) 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 (62), i.e., −(D ⊗ D) 2 , is still preserved. 5

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, f 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 ∇H the vector of the functional derivatives of the Hamiltonian functional 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: 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: for suitable time dependent coefficients, η 0 (t), η 1 (t), µ 1 (t), . . . . By using the infinite vectors (47) and we can cast the expansions of u and v in vector form, respectively, as (48) and v(x, t) = ω(x) p(t).
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 infinitedimensional Hamiltonian ODE probleṁ This latter problem is Hamiltonian w.r.t. the Hamiltonian which turns out to be equivalent to the Hamiltonian functional (65). Moreover, the two quadratic invariants (66) can be respectively rewritten as whereD 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 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 having set q 1 ≈ q(h) and p 1 ≈ p(h) the new approximations.

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)) On the other hand, the norm of the matrix , and tends to infinity, as N → ∞, 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 Σ approximated as in (34) and, in the present context, given by which is a block matrix with diagonal blocks and, therefore, Σ −1 can be cheaply computed and stored. As matter of fact, one has [3, Theorem 5]: 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 1.

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, with αβ = 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 with H[u] the Hamiltonian functional and its functional derivative. 6 Because of the periodic boundary conditions, the Hamiltonian functional (77) turns out to be conserved. Another conserved functional is given by 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), Consequently, the expansion (46) now becomes In order to put this expansion in vector form, let us introduce the infinite vectors (80) In so doing, we can rewrite (79) as: Consequently, the conservation of (78) is automatically granted. Moreover, by defining the infinite matrix such that and similarly for the higher derivatives, and considering that one verifies that (76) can be rewritten as the infinite dimensional ODE problem (see [16,Lemma 3] for full details) For this problem, the following result holds true [16, Theorem 1].
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 are hereafter used in place of (80). Similarly, by replacing matrix (82) with 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 → ∞, upon regularity assumptions on u 0 . 7 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).

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 with u(x, t) given by the expansion (81). Consequently, by considering the constant approximation u(x, t) ≡û 0 (due to the conservation of (78)), and taking into account (84), one obtains the constant approximate (diagonal) Hessian Therefore, the blended iteration (27)- (29) can be conveniently used, by considering the resulting approximated matrix (see (41) and (88)) which is a block matrix with diagonal blocks. Moreover, one has [16,Theorem 3]: 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 ).

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 8 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 with ε the machine epsilon, for problems (51) and (70), or for problem (85); • a good approximation of the Hamiltonian. This is achieved by computing the initial value H(q 0 , 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, ∆H 0 , satisfies:

The semilinear wave equation
We consider the so called sine-Gordon equation [13, Section 7] with a breather soliton solution, where we choose γ = 1.5. Its solution, depicted in the upper plot in Figure 1, is: In the lower plot of Figure 1 there are the graphs of E 0 and ∆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 ≈ √ ε 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.

The nonlinear Schrödinger equation
We consider the so called focusing equation, 9 where the initial conditions at t = 0 are taken from the known solution, 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 ∆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 ≈ 10 −1 √ ε in (38). In Table 2 we list the numerical results obtained by solving 9 The de-focusing case is obtained when the sign of the coupling term is reversed. 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 energyconserving 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).

The Korteweg-de Vries equation
This example is adapted from [16, Example 2]: equipped with periodic boundary conditions and the initial condition obtained from the known cnoidal wave solution, Here cn := 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: = 10 −2 , m = 0.9, a = 192m K 2 (m), ν = 64 (2m−1)K 2 (m), x 0 = 1/2.
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 ∆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. 10 For the SHBVM method, we also list the used values of k and s, the latter obtained by using tol ≈ 10 −1 √ ε 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.

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 energyconservation 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 secondorder 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 [2,21,38], 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.

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 [17,40]. In the future, we plan to further investigate Hamiltonian PDEs within the same framework.