Line Integral Solution of Differential Problems

In recent years, the numerical solution of differential problems, possessing constants of motion, has been attacked by imposing the vanishing of a corresponding line integral. The resulting methods have been, therefore, collectively named (discrete) line integral methods, where it is taken into account that a suitable numerical quadrature is used. The methods, at first devised for the numerical solution of Hamiltonian problems, have been later generalized along several directions and, actually, the research is still very active. In this paper we collect the main facts about line integral methods, also sketching various research trends, and provide a comprehensive set of references.


Introduction
The numerical solution of differential problems in the form ẏ(t) = f (y(t)), t ≥ 0, is needed in a variety of applications.In many relevant instances, the solution has important geometric properties and the name geometric integrator has been coined to denote a numerical method able to preserve them (see, e.g., the monographs [1][2][3][4]).Often, the geometric properties of the vector field are summarized by the presence of constants of motion, namely functions of the state vector which are conserved along the solution trajectory of (1).For this reason, in such a case one speaks about a conservative problem.For sake of simplicity, let us assume, for a while, that there exists only one constant of motion, say C(y(t)) ≡ C(y 0 ), ∀t ≥ 0, ∀y 0 ∈ D, along the solution y(t) of (1).Hereafter, we shall assume that both f : D → R m and C : D → R are suitably smooth (e.g., analytical).In order for the conservation property (2) to hold, one requires that d dt C(y(t)) = ∇C(y(t)) ẏ(t) = ∇C(y(t)) f (y(t)) = 0.
In fact, since y satisfies (1), one obtains that the integrand is given by ∇C(y(τ)) f (y(τ)) ≡ 0, because of (3).On the other hand, if one is interested in obtaining an approximation to y, ruled by a discrete-time dynamics with time-step h, one can look for any path σ joining y 0 to y 1 ≈ y(h), i.e., and such that Definition 1.The path σ satisfying ( 5) and ( 6) defines a line integral method providing an approximation y 1 to y(h) such that C(y 1 ) = C(y 0 ).
Obviously, the process is then repeated on the interval [h, 2h], starting from y 1 , and so on.We observe that the path σ now provides the vanishing of the line integral in (6) without requiring the integrand be identically zero.This, in turn, allows much more freedom during the derivation of such methods.In addition to this, it is important to observe that one cannot, in general, directly impose the vanishing of the integral in (6) since, in most cases, the integrand function does not admit a closed form antiderivative.Consequently, in order to obtain a ready to use numerical method, the use of a suitable quadrature rule is mandatory.
Since we shall deal with polynomial paths σ, it is natural to look for an interpolatory quadrature rule defined by the abscissae and weights (c i , b i ), i = 1, . . ., k.In order to maximize the order of the quadrature, i.e., 2k, we place the abscissae at the zeros of the kth shifted and scaled Legendre polynomial P k (i.e., P k (c i ) = 0, i = 1, . . ., k).Such polynomials provide an orthonormal basis for functions in L 2 [0, 1]: ∀i, j = 0, 1, . . ., (7) with δ ij denoting the Kronecker delta.Consequently, (6) becomes Definition 2. The path σ satisfying ( 5) and ( 8) defines a discrete line integral method providing an approximation y 1 to y(h) such that C(y 1 ) ≈ C(y 0 ), within the accuracy of the quadrature rule.
As is clear, if C is such that the quadrature in (8) is exact, then the method reduces to the line integral method satisfying (5) and (6), exactly conserving the invariant.In the next sections we shall make the above statements more precise and operative.
With these premises, the paper is organized as follows: in Section 2 we shall deal with the numerical solution of Hamiltonian problems; Poisson problems are then considered in Section 3; constrained Hamiltonian problems are studied in Section 4; Hamiltonian partial differential equations (PDEs) are considered in Section 5; highly oscillatory problems are briefly discussed in Section 6; at last, Section 7 contains some concluding remarks.

Hamiltonian Problems
A canonical Hamiltonian problem is in the form ẏ = J∇H(y), y(0 with H the Hamiltonian function and, in general, I r hereafter denoting the identity matrix of dimension r.Because of the skew-symmetry of J, one readily verifies that H is a constant of motion for (9): For isolated mechanical systems, H has the physical meaning of the total energy, so that it is often referred to as the energy.When solving (9) numerically, it is quite clear that this conservation property becomes paramount to get a correct simulation of the underlying phenomenon.The first successful approach in the numerical solution of Hamiltonian problems has been the use of symplectic integrators.The characterization of a symplectic Runge-Kutta method c A b is based on the following algebraic property of its Butcher tableau [29,30] (see also [31]) which is tantamount to the conservation of any quadratic invariant of the continuous problem.Under appropriate assumptions, symplectic integrators provide a bounded Hamiltonian error over long time intervals [2], whereas generic numerical methods usually exhibit a drift in the numerical Hamiltonian.Alternatively, one can look for energy conserving methods (see, e.g., [32][33][34][35][36][37][38][39]).We here sketch the line integral solution to the problem.According to (5) and (6) where the coefficients γ j ∈ R 2m are at the moment unknown.Integrating term by term, and imposing the initial condition, yields the following polynomial of degree s: By defining the approximation to y(h) as where we have taken into account the orthonormality conditions (7), so that 1 0 P j (x)dx = δ j0 , one then obtains that the conditions (5) are fulfilled.In order to fulfil also (6) with C = H, one then requires, by taking into account ( 11) This latter equation is evidently satisfied, due to the skew-symmetry of J, by setting Consequently, (12) becomes which, according to ([12], Definition 1), is the master functional equation defining σ.Consequently, the conservation of the Hamiltonian is assured.Next, we discuss the order of accuracy of the obtained approximation, namely the difference σ(h) − y(h): this will be done in the next section, by using the approach defined in [18].

Local Fourier Expansion
By introducing the notation one has that (9) can be written, on the interval [0, h], as with γ j (y) defined according to (16), by formally replacing σ with y.Similarly, with the polynomial σ in (15) satisfying, by virtue of ( 11) and ( 14), the differential equation: The following result can be proved (see [18], Lemma 1).
Lemma 1.Let g : [0, h] → V, with V a vector space, admit a Taylor expansion at 0. Then, for all j ≥ 0 : Moreover, let us denote by y(t, t, ỹ) the solution of the ODE-IVPs and by Φ(t, t) the fundamental matrix solution of the variational problem associated to it.We recall that We are now in the position to prove the result concerning the accuracy of the approximation (13).
) (in other words, the polynomial σ defines an approximation procedure of order 2s).

Hamiltonian Boundary Value Methods
Quoting Dahlquist and Björk [40], p. 521, as is well known, even many relatively simple integrals cannot be expressed in finite terms of elementary functions, and thus must be evaluated by numerical methods.In our framework, this obvious statement means that, in order to obtain a numerical method from (15), we need to approximate the integrals appearing in that formula by means of a suitable quadrature procedure.In particular, as anticipated above, we shall use the Gauss-Legendre quadrature of order 2k, whose abscissae and weights will be denoted by (c i , b i ) (i.e., P k (c i ) = 0, i = 1, . . ., k).Hereafter, we shall obviously assume k ≥ s.In so doing, in place of (15), one obtains a (possibly different) polynomial, where (see ( 16)): with ∆ j (h) the quadrature error.Considering that the quadrature is exact for polynomial integrands of degree 2k − 1, one has: As a consequence, u ≡ σ if H is a polynomial of degree ν ≤ 2k/s.In such a case, H(u(h)) − H(u(0)) ≡ H(σ(h)) − H(σ(0)) = 0, i.e., the energy is exactly conserved.Differently, one has, by virtue of (22) and Lemma 1: The result of Theorem 1 continues to hold for u.In fact, by using arguments similar to those used in the proof of that theorem, one has, by taking into account (22) and that k ≥ s: Actually, by observing that in (21) only the values of u at the abscissae are needed, one obtains, by setting Y i := u(c i h), and rearranging the terms: with the new approximation given by Consequently, we are speaking about the k-stage Runge-Kutta method with Butcher tableau given by: c and The next result summarizes the properties of HBVMs sketched above, where we also take into account that the abscissae are symmetrically distributed in the interval [0, 1] (we refer to [1,18]  We conclude this section by showing that, for HBVM(k, s), whichever is the value k ≥ s considered, the discrete problem to be solved has dimension s, independently of k.This fact is of paramount importance, in view of the use of relatively large values of k, which are needed, in order to gain a (at least practical) energy conservation.In fact, even for non polynomial Hamiltonians, one obtains a practical energy conservation, once the O(h 2k+1 ) Hamiltonian error falls, by choosing k large enough, within the round-off error level.
By taking into account the stage Equation (24), and considering that Y i = u(c i h), one has that the stage vector can be written as where is the block vector (of dimension s) with the coefficients (22) of the polynomial u in (21).By combining ( 29) and (31), one then obtains the discrete problem, equivalent to (24), having (block) dimension s, independently of k.Once this has been solved, one verifies that the new approximation ( 25) turns out to be given by (compare also with ( 13)): Next section will concern the efficient numerical solution of the discrete problem generated by the application of a HBVM(k, s) method.In fact, a straightforward fixed-point iteration, may impose severe stepsize limitations.On the other hand, the application of the simplified Newton iteration for solving (34) reads, by considering that (see ( 27) and ( 28)) and setting f (y 0 ) the Jacobian of f evaluated at y 0 : This latter iteration, in turn, needs to factor a matrix whose size is s times larger than that of the continuous problem.This can represent an issue, when large-size problems are to be solved and/or large values of s are considered.

Blended Implementation of HBVMs
We here sketch the main facts concerning the so called blended implementation of HBVMs, a Newton-like iteration alternative to (36), which only requires to factor a matrix having the same size as that of the continuous problem, thus resulting into a much more efficient implementation of the methods [1,16].This technique derives from the definition of blended implicit methods, which have been at first considered in [41,42], and then developed in [43][44][45].Suitable blended implicit methods have been implemented in the Fortran codes BIM [46], solving stiff ODE-IVPs, and BIMD [47], also solving linearly implicit DAEs up to index 3.The latter code is also available on the Test Set for IVP Solvers [48] (see also [49]), and turns out to be among the most reliable and efficient codes currently available for solving stiff ODE-IVPs and linearly implicit DAEs.It is worth mentioning that, more recently, the blended implementation of RKN-type methods has been also considered [50].
Let us then consider the iteration (36), which requires the solution of linear systems in the form By observing that matrix X s defined at (35) is nonsingular, we can consider the equivalent linear system where ρ s is a positive parameter to be determined.For this purpose, let be the Jordan canonical form of f (y 0 ).For simplicity, we shall assume that Λ is diagonal, and let λ be any of its diagonal entries.Consequently, the two linear systems ( 37) and (38), projected in the invariant subspace corresponding to that entry, respectively become again being equivalent to each other (i.e., having the same solution x ∈ R s ).We observe that the coefficient matrix of the former system is I s + O(q), when q ≈ 0, whereas that of the latter one is −ρ s q(I s + O(q −1 )), when |q| 1.Consequently, one would like to solve the former system, when q ≈ 0, and the latter one, when |q| 1.This can be done automatically by considering a weighting function θ(q) such that then considering the blending of the two equivalent systems (39) with weights θ(q) and I s − θ(q), respectively: In particular, ( 40) can be accomplished by choosing Consequently, one obtains that M(q) = I s + O(q), q ≈ 0, and M(q) = −ρ s q(I s + O(q −1 )), |q| 1.
As a result, one can consider the following splitting for solving the problem: The choice of the scalar parameter ρ s is then made in order to optimize the convergence properties of the corresponding iteration.According to the analysis in [42], we consider where, as is usual, σ(X s ) denotes the spectrum of matrix X s .A few values of ρ s are listed in Table 1.Coming back to the original iteration (36), one has that the weighting function (42) now becomes which requires to factor only the matrix having the same size as that of the continuous problem, thus obtaining the following blended iteration for HBVMs: It is worth mentioning that: • in the special case of separable Hamiltonian problems, the blended implementation of the methods can be made even more efficient, since the discrete problem can be cast in terms of the generalized coordinates only (see [16] or ( [1], Chapter 4)); • the coding of the blended iteration becomes very high-performance by considering a matrix formulation of (45) (see, e.g., ([1], Chapter 4.2.2) or [51]).As matter of fact, it has been actually implemented in the Matlab code hbvm, which is freely available on the internet at the url [52].
In order to give evidence of the usefulness of energy conservation, let us consider the solution of the well-known pendulum problem, with Hamiltonian When considering the trajectory starting at [1,53] one obtains a periodic solution of period T ≈ 28.57109480185544.In Table 2 we list the obtained results when solving the problem over 10 periods, by using a stepsize h = T/n, with HBVM(6,3) and HBVM(3,3) (i.e., the symplectic 3-stage Gauss collocation method), for increasing values of n.
As one may see, even though both methods are 6th order accurate, nevertheless, HBVM(6,3) becomes (practically) energy-conserving as soon as n ≥ 40, whereas HBVM(3,3) does not.One clearly sees that, for the problem at hand, the energy-conserving method is pretty more accurate than the non conserving one.

Energy and QUadratic Invariants Preserving (EQUIP) Methods
According to Theorem 2, when k = s, HBVM(s, s) reduces to the s-stage Gauss method.For such a method, one has, with reference to ( 26)-( 28) and ( 35) with k = s, so that the Butcher matrix in (26) becomes the W-transformation [54] of the s-stage Gauss method, i.e., the Butcher matrix is A = P s X s P −1 s .Moreover, since A = P s X s P s Ω, the method is easily verified to be symplectic, because of (10).In fact, by setting in general e i ∈ R s the ith unit vector, and with reference to the vector e defined in ( 30) with k = s, one has: We would arrive at the very same conclusion if we replace X s by In particular, if the parameter α is small enough, matrix αV will act as a perturbation of the underlying Gauss formula and the question is whether it is possible to choose α, at each integration step, such that the resulting integrator may be energy conserving.In order for V = O, we need to assume, hereafter, s ≥ 2. In particular, by choosing it is possible to show [26] that the scalar parameter α can be chosen, at each integration step, such that, when solving the Hamiltonian problem ( 9) with a sufficiently small stepsize h: the method retains the order 2s of the original s-stage Gauss method, This fact is theoretically intriguing, since this means that we have a kind of state-dependent Runge-Kutta method, defined by the Butcher tableau which is energy conserving and is defined, at each integration step, by a symplectic map, given by a small perturbation of that of the underlying s-stage Gauss method.Consequently, EQUIP methods do not infringe the well-known result about the nonexistence of energy conserving symplectic numerical methods [55,56].Since the symplecticity condition ( 10) is equivalent to the conservation of all quadratic invariants of the problem, these methods have been named Energy and QUadratic Invariants Preserving (EQUIP) methods [26,57].It would be interesting to study the extent to which the solutions generated by an EQUIP method inherit the good long time behavior of the associated Gauss integrator with reference to the nearly conservation property of further non-quadratic first integrals.This investigation would likely involve a backward error analysis approach, similar to that done in [2], and up to now remains an open question.
For a thorough analysis of such methods we refer to [26,53].In the next section, we sketch their line-integral implementation when solving Poisson problems, a wider class than that of Hamiltonian problems.

Poisson Problems
Poisson problems are in the form When B(y) ≡ J as defined in (9), then one retrieves canonical Hamiltonian problems.As in that case, since B(y) is skew-symmetric, then H, still referred to as the Hamiltonian, is conserved: Moreover, any scalar function C(y) such that ∇C(y) B(y) = 0 is also conserved, since: C is called a Casimir function for (50).Consequently, all possible Casimirs and the Hamiltonian H are conserved quantities for (50).In the sequel, we show that EQUIP methods can be conveniently used for solving such problems.As before, the scalar parameter α in (49) is selected in such a way that the Hamiltonian H is conserved.Moreover, since the Butcher matrix in ( 49) satisfies (10), then all quadratic Casimirs turn out to be conserved as well.The conservation of all quadratic invariants, in turn, is an important property as it has been observed in [58].
Let us then sketch the choice of the parameter α to gain energy conservation (we refer to [53] for full details).By setting . . .
one has that the Butcher matrix in ( 49) can be written as Consequently, by denoting f (y) = B(y)∇H(y), and setting Y i := u(c i h), i = 1, . . ., s, the stages of the method, one obtains that the polynomial u(ch) is given by with the (block) vectors γj formally defined as in ( 22) with k = s.In vector form, one has then (compare with (31)): We observe that, from (51), one also obtains: Nevertheless, the new approximation, still given by now differs from Consequently, in order to define a path joining y 0 to y 1 , to be used for imposing energy-conservation by zeroing a corresponding line-integral, we can consider the polynomial path made up by u plus As a result, by considering that w(1) = y 1 , w(0) = u(h), u(0) = y 0 , we shall choose α such that (see ( 51)-( 56)): In more details, by defining the vectors and resorting to the usual line integral argument, it is possible to prove the following result ([53], Theorem 2).
Theorem 3. (57) holds true, provided that As in the case of HBVMs, however, we shall obtain a practical numerical method only provided that the integrals in (58) are suitably approximated by means of a quadrature which, as usual, we shall choose as the Gauss-Legendre formula of order 2k.In so doing, one obtains an EQUIP(k, s) method.
The following result can be proved ( [53], Theorem 8).Theorem 4.Under suitable regularity assumptions on both B(y) and H(y), one has that for all k ≥ s, the EQUIP(k, s) method has order 2s, conserves all quadratic invariants and, moreover, We observe that, for EQUIP(k, s), even a not exact conservation of the energy may result in a much better error growth, as the next example shows.We consider the Lotka-Volterra problem [53], which is in the form (50) with By choosing the following parameters and initial values, one obtains a periodic solution of period T ≈ 7.720315563434113.If we solve the problems ( 60) and (61) with the EQUIP(6,3) and the 3-stage Gauss methods with stepsize h = T/50 over 100 periods, we obtain the error growths, in the numerical solution, depicted in Figure 1.As one may see, the EQUIP(6,3) method (which only approximately conserves the Hamiltonian), exhibits a linear error growth; on the contrary, the 3-stage Gauss method (which exhibits a drift in the numerical Hamiltonian) has a quadratic error growth.Consequently, there is numerical evidence that EQUIP methods can be conveniently used for numerically solving Poisson problems (a further example can be found in [53]).Error growth over 100 periods when solving problems ( 60) and (61) by using the EQUIP(6,3) method (blue solid line) and the 3-stage Gauss method (red dashed line) with stepsize h = T/50 ≈ 0.15.The dotted lines show the linear and quadratic error growths.

Constrained Hamiltonian Problems
In this section, we report about some recent achievements concerning the line integral solution of constrained Hamiltonian problems with holonomic contraints [23].This research is at a very early stage and, therefore, it is foreseeable that new results will follow in the future.
To begin with, let us consider the separable problem defined by the Hamiltonian with M a symmetric and positive definite matrix, subject to the ν < m holonomic constraints We shall assume that the points are regular for the constraints, so that ∇g(q) has full column rank and, therefore, the ν × ν matrix ∇g(q) M −1 ∇g(q) is nonsingular.By introducing the vector of the Lagrange multipliers λ ∈ R ν , problems ( 62) and ( 63) can be equivalently cast in Hamiltonian form by defining the augmented Hamiltonian Ĥ(q, p, λ) = H(q, p) + λ g(q), ( thus obtaining the equivalent constrained problem: subject to the initial conditions q(0) = q 0 , p(0) = p 0 , such that g(q 0 ) = 0, ∇g(q 0 ) M −1 p 0 = 0. (66) We observe that the first requirement (g(q 0 ) = 0) obviously derives from the given constraints.The second, in turn, derives from 0 = ġ(q) = ∇g(q) q = ∇g(q) M −1 p, which has to be satisfied by the solution of (65).These latter constraints are sometimes referred to as the hidden constraints.A formal expression of the vector of the Lagrange multipliers can be obtained by further differentiating the previous expression, thus giving ∇g(q) M −1 ∇g(q) λ = ∇ 2 g(q) M −1 p, M −1 p − ∇g(q) M −1 ∇U(q), which is well defined, because of the assumption that ∇g(q) M −1 ∇g(q) is nonsingular.Consequently, where the notation λ(q, p) means that the vector λ is a function of q and p.It is easily seen that both the two Hamiltonians ( 62) and ( 64) are conserved along the solution of the problems (65) and (66), and assume the same value.For numerically solving the problem, we shall consider a discrete mesh with time-step h, i.e., t n = nh, n = 0, 1, . . ., looking for approximations such that, starting from (q n , p n ), one arrives at (q n+1 , p n+1 ) by choosing λ n in order for: Consequently, the new approximation conserves the Hamiltonan and exactly satisfies the constraints, but only approximately the hidden contraints.In particular, we shall consider a piecewise constant approximation of the vector of the Lagrange multipliers λ, i.e., λ n is assumed to be constant on the interval [t n , t n+1 ].In other words, we consider the sequence of problems (compare with (65)), for n = 0, 1, . . .: where the constant vector λ n is chosen in order to satisfy the constraints at t n+1 .That is, such that the new approximations, defined as satisfy (68).The reason for choosing λ n as a constant vector stems from the following result, which concerns the augmented Hamiltonian (64).
Proof.In fact, denoting by Ĥq (q, p, λ) the gradient of Ĥ with respect to the q variables, and similarly for Ĥp , the usual line integral argument provides: Consequently, one obtains that g(q n+1 ) = g(q n ) ⇔ H(q n+1 , p n+1 ) = H(q n , p n ), i.e., energy conservation is equivalent to satisfy the constraints.
In order to fulfil the constraints, we shall again resort to a line integral argument.For this purpose, we need to define the Fourier coefficients (along the Legendre basis ( 7)) of the functions appearing at the right-hand sides in (69): so that, in particular, We also need the following result.
We are now in the position of deriving a formal expression for the constant approximation to the vector of the Lagrange multipliers through the usual line integral approach: Because of Lemma 2, one then obtains that (see (35)), The following result can be proved [23].
Theorem 6.The Equation ( 73) is consistent with (67) and is well defined for all sufficiently small h > 0.
At this point, in order to obtain a numerical method, the following two steps need to be done: 1. truncate the infinite series in (72) to finite sums, say up to j = s − 1.In so doing, the expression of λ n changes accordingly, since in (73) the infinite sums will consequently arrive up to j = s − 1; 2. approximate the integrals in (71), for j = 0, . . ., s − 1.As usual, we shall consider a Gauss-Legendre formula of order 2k, with k ≥ s, thus obtaining corresponding approximations which we denote, for n = 0, 1, . . ., In so doing, it can be seen that one retrieves the usual HBVM(k, s) method defined in Section 2, applied for solving (69), coupled with the following equation for λ n , representing the approximation of (73): The following result can be proved [23].
Theorem 7.For all sufficiently small stepsizes h, the HBVM(k, s) method coupled with (74), used for solving (65) and (66) over a finite interval, is well defined and symmetric.It provides a sequence of approximations (q n , p n , λ n ) such that (see (67)): Moreover, in the case where λ(q(t), p(t)) ≡ λ, for all t ≥ 0, then It is worth mentioning that, even in the case where the vector of the Lagrange multipliers is not constant, so that all HBVM(k, s) are second-order accurate for all s ≥ 1, the choice s > 1 generally provides a much smaller solution error.
We refer to [23] for a number of examples of application of HBVMs to constrained Hamiltonian problems.We here only provide the application of HBVM(4,4) (together with (74) with s = 4) for solving the so called conical pendulum problem.In more details, let us have a pendulum of unit mass connected to a fixed point (the origin) by a massless rod of unit length.The initial conditions are such that the motion is periodic of period T and takes place in the horizontal plane q 3 = z 0 .Normalizing the acceleration of gravity, the augmented Hamiltonian is given by Ĥ(q, p, λ) = 1 2 p p + e 3 q + λ(q q − 1) ≡ H(q, p) + λg(q), where p = q and e i ∈ R 3 is the ith unit vector.Choosing q(0) = 2 − 1 2 (e 1 − e 3 ), p(0) = 2 − 1 4 e 2 , results in Since the augmented Hamiltonian Ĥ is quadratic and λ is constant, any HBVM(s, s) method, coupled with (74), is energy and constraint conserving, and of order 2s.If we apply the HBVM(4,4) method for solving the problem over 10 periods by using the stepsize h = T/N, N = 10, 20, 40, it turns out the λ n ≡ λ, g(q n ) = 0, ∇g(q n ) p n = 0, and H(q n , p n ) = H(q 0 , p 0 ) within the round-off error level, for all n = 0, 1, . . ., 10N.On the other hand, the corresponding solution errors, after 10 periods, turn out to be given by 4.9944 × 10 −8 , 1.9676 × 10 −10 , 7.3944 × 10 −13 , thus confirming the order 8 of convergence of the resulting method.

Hamiltonian PDEs
Quoting [4], p. 187, the numerical solution of time dependent PDEs may often be conceived as consisting of two parts.First the spatial derivatives are discretized by finite differences, finite elements, spectral methods, etc. to obtain a system of ODEs, with t as the independent variable.Then this system of ODEs is integrated numerically.If the PDEs are of Hamiltonian type, one may insist that both stages preserve the Hamiltonian structure.Thus the space discretization should be carried out in such a way that the resulting system of ODEs is Hamiltonian (for a suitable Poisson bracket).This approach has been systematically used for solving a number of Hamiltonian PDEs by using HBVMs [1,20,22,59,60] and this research is still under development.We here sketch te main facts for the simplest possible example, provided by the semilinear wave equation, with f the derivative of f , coupled with the initial conditions and periodic boundary conditions.We shall assume that f , φ 0 , φ 1 are regular enough (the last two functions, as periodic functions).By setting (hereafter, when not necessary, we shall avoid the arguments of the functions) , the problem can be cast into Hamiltonian form as i.e., where δH δz is the vector of the functional derivatives [2,3,22] of the Hamiltonian functional Because of the periodic boundary conditions, this latter functional turns out to be conserved.In fact, by considering that because of the periodicity in space.Consequently, Moreover, since u(x, t) is periodic for x ∈ [a, b], we can expand it in space along the following slight variant of the Fourier basis, so that, for all allowed i, j: In so doing, for suitable time dependent coefficients γ j (t), η j (t), one obtains the expansion: having introduced the infinite vectors the identity operator, and introducing the infinite matrix (see (77)) so that ω (x) = Dω(x) and ω (x) = D 2 ω(x) = −D Dω(x), one then obtains that (79) can be rewritten as the infinite system of ODEs: subject to the initial conditions (see (76)) The following result is readily established.
Theorem 8. Problem (87) is Hamiltonian, with Hamiltonian This latter is equivalent to the Hamiltonian functional (80), via the expansion (83).
Proof.The first part of the proof is straightforward.Concerning the second part, one has, by virtue of (85): In order to obtain a computational procedure, one needs to truncate the infinite series in (83) to a finite sum, i.e., u(x, t) ≈ c 0 (x)γ 0 (t) + N ∑ j=1 c j (x)γ j (t) + s j (x)η j (t) =: u N (x, t). (90) Clearly, such an approximation no longer satisfies the wave Equation (79).Nevertheless, in the spirit of Fourier-Galerkin methods [61], by requiring that the residual be orthogonal to the functional subspace V N = span {c 0 (x), c 1 (x), s 1 (x), . . ., c N (x), s N (x)} containing the approximation u(x, t) for all times t ≥ 0, one obtains a finite dimensional ODE system, formally still given by ( 87) and (88), upon replacing the involved infinite vectors and matrices, previously defined in (84) and (86), with the following finite dimensional ones (of dimension 2N + 1): . . .
Moreover, the result of Theorem 8 continues formally to hold for the finite dimensional problem, with the sole exception that now the Hamiltonian (89) only yields an approximation to the continuous functional (80).Nevertheless, it is well known that, under suitable regularity assumptions on f and the initial data, this truncated version converges exponentially to the original functional (80), as N → ∞ (this phenomenon is usually referred to as spectral accuracy).Since problem (87) is Hamiltonian, one can use HBVM(k, s) methods for solving it.It is worth mentioning that, in so doing, the blended implementation of the methods can be made extremely efficient, by considering that: • an accurate approximation in space usually requires the use of large values of N; • as a consequence, in most cases one has in a suitable neighbourhood of the solution.
Consequently, the Jacobian matrix of (87) can be approximated by the linear part alone, i.e., as This implies that matrix Σ involved in the definition of Θ in ( 44) becomes (all the involved matrices are diagonal and have dimension 2N + 1): where ρ s is the parameter defined in (43), and h is the used time-step.As a result, Σ: • is constant for all time steps, so that it has to be computed only once; • it has a block diagonal structure (and, in particular, D is positive definite).
The above features make the resulting blended iteration (45) inexpensive, thus allowing the use of large values of N and large time-steps h.As an example, let us solve the sine-Gordon equation [22], with γ > 0, whose solution is, by considering the value γ = 1.5, In such a case, the value N = 300 in (90) and ( 91) is sufficient to obtain an error in the spatial semi-discretization comparable with the round-off error level.Then, we solve in time the semi-discrete problem (87), of dimension 2N + 2 = 602, by using the HBVM(k, s) methods.In Table 3 we list the obtained maximum errors in the computed solution, by using a time-step h = 1, along with the corresponding Hamiltonian errors and execution times, for various choices of (k, s) (in particular, for k = s, s = 1, 2, 3, we have the s-stage Gauss collocation methods, which are symplectic but not energy conserving) (all numerical tests have been performed on a laptop with a 2.2 GHz dual core i7 processor, 8 GB of memory, and running Matlab R2017b).From the listed results, one sees that HBVM(2s, s) methods become energy-conserving for s ≥ 4 and the error decreases until the round-off error level, for s = 10, with a computational time 10 times larger than that of the implicit mid-point rule (obtained for k = s = 1) which, however, has a solution error 10 13 times larger.

Highly Oscillatory Problems
The Hamiltonian system of ODEs (87) is a particular instance of problems in the form where, without loss of generality, A is a symmetric and positive definite matrix and f is a scalar function such that in a neighbourhood of the solution.Moreover, hereafter we consider the 2-norm, so that ω equals the largest eigenvalue of A (in general, any convenient upper bound would suffice).The problem is Hamiltonian, with Hamiltonian Problems in the form (94) satisfying (95) are named (multi-frequency) highly oscillatory problems, since they are ruled by the linear part, possessing large (possibly different) complex conjugate eigenvalues.Such problems have been investigated since many years, starting from the seminal papers [62,63], and we refer to the monograph [64] for more recent findings.A common feature of the methods proposed so far, however, is that of requiring the use of time-steps h such that hω < 1, either for stability and/or accuracy requirements.We here sketch the approach recently defined in [27], relying on the use of HBVMs, which will allow the use of stepsizes h without such a restriction.
To begin with, and for analysis purposes, let us recast problem (94) in first order form, by setting with J 2 the 2 × 2 skew-symmetric matrix defined in (77).As is usual in this context, one considers at first the linear part of (97), The solution of (98), y(t) = e J 2 ⊗A t y 0 , on the interval [0, h] is readily seen to be given, according to the local Fourier expansion described in Section 2.1, by: However, when using a finite precision arithmetic with machine epsilon u, the best we can do is to approximate where hereafter .
= means "equal within round-off error level", provided that γj (y) In fact, further terms in the infinite series in (99) would provide a negligible contribution, in the used finite precision arithmetic.By defining the function with J j+ 1 2 (•) the Bessel functions of the first kind, it can be shown ( [27], Criterion 1) that the requirement (101) is accomplished by requiring In so doing, one implicitly defines a function ϕ u , depending on the machine epsilon u, such that The plot of such a function is depicted in Figure 2 for the double precision IEEE: as one may see, the function is well approximated by the line [27] 24 + 0.7 • ωh. (105) Next, one consider the whole problem (97), whose solution, on the interval [0, h], can be written as As observed before, when using a finite precision arithmetic with machine epsilon u, the best we can do is to approximate (106) with a polynomial: provided that Assuming the ansatz f (y(ch)) ∼ e νJ 2 ⊗A ch , for a suitable ν ≥ 1 (essentially, one requires that ∇ f is well approximated by a polynomial of degree ν), the requirement (109) is accomplished by choosing ( [27], Criterion 2) s = ϕ u (νωh), where ϕ u is the same function considered in (104).By taking into account that ν ≥ 1 and that ϕ u is an increasing funcion (see Figure 2), one then obtains that s ≥ s 0 .
For solving it, we shall use the SHBVM method with parameters: and a time-step h = 20/N, for various values of N, as specified in Table 4.In that table we also list the corresponding: • maximum absolute error in the computed solution, e q ; • the maximum relative error in the numerical Hamiltonian, e H ; • the value of ωh (which is always much greater than 1); • the parameters (s 0 , s, k), computed according to (104), (110), and (111), respectively; • the execution time (in sec).
As before, the numerical tests have been performed on a laptop with a 2.2 GHz dual core i7 processor, 8 GB of memory, and running Matlab R2017b.(27,45,47) 3.59 2000 10.0 2.89 × 10 −10 3.33 × 10 −16 (26,44,46) 3. 46 As one may see, the method is always energy conserving and the error, as expected, is uniformly small.Also the execution times are all very small (of the order of 3.5 s).It must be emphasized that classical methods, such as the Gautschi or the Deuflhard method, would require the use of much smaller time-steps, and much larger execution times (we refer to [27] for some comparisons).

Conclusions
In this paper, we have reviewed the main facts concerning the numerical solution of conservative problems within the framework of (discrete) line integral methods.Relevant instances of line integral methods are provided by the energy-conserving Runge-Kutta methods named Hamiltonian Boundary Value Methods (HBVMs) and the Energy and QUadratic Invariants Preserving (EQUIP) methods, providing efficient geometric integrators for Hamiltonian problems.It is worth mentioning that HBVMs can be easily adapted to efficiently handle Hamiltonian BVPs [19], whereas EQUIP methods can be also used for numerically solving Poisson problems.In this paper, we have also reviewed some active research trends, concerning the application of HBVMs for numerically solving constrained Hamiltonian problems, Hamiltonian PDEs, and highly-oscillatory problems.The effectiveness of HBVMs is emphasized by the availability of an efficient nonlinear iteration for solving the generated discrete problems, relying on the so called blended implementation of the methods.Matlab software implementing HBVMs for Hamiltonian problems is available at the url [52].

Definition 3 .
The polynomial u defined at(21) defines a Hamiltonian Boundary Value Method (HBVM) with parameters k and s, in short HBVM(k, s).

Figure 1 .
Figure 1.Error growth over 100 periods when solving problems(60) and (61) by using the EQUIP(6,3) method (blue solid line) and the 3-stage Gauss method (red dashed line) with stepsize h = T/50 ≈ 0.15.The dotted lines show the linear and quadratic error growths.

Figure 2 .
Figure 2. Plot of the function ϕ u defined in (104) versus ωh (blue plus line), for the double precision IEEE, together with its linear approximation (105) (red solid line).
for full details).

Table 3 .
Solution error e u and Hamiltonian error e H when solving problem (93) with stepsize h = 1 and N = 300 in (90), by using HBVM(k, s) methods.

Table 4 .
Duffing problem (112)-(115) solved by the SHBVM(k, s, s 0 ) method with parameters (116) and time-step h = 20/N: e q is the maximum absolute error in the computed solution; e H is the maximum relative error in the numerical Hamiltonian.