On a Class of Conjugate Symplectic Hermite-Obreshkov One-Step Methods with Continuous Spline Extension

The class of A-stable symmetric one-step Hermite–Obreshkov (HO) methods introduced by F. Loscalzo in 1968 for dealing with initial value problems is analyzed. Such schemes have the peculiarity of admitting a multiple knot spline extension collocating the differential equation at the mesh points. As a new result, it is shown that these maximal order schemes are conjugate symplectic, which is a benefit when the methods have to be applied to Hamiltonian problems. Furthermore, a new efficient approach for the computation of the spline extension is introduced, adopting the same strategy developed for the BS linear multistep methods. The performances of the schemes are tested in particular on some Hamiltonian benchmarks and compared with those of the Gauss–Runge–Kutta schemes and Euler–Maclaurin formulas of the same order.


Introduction
We are interested in the numerical solution of the Cauchy problem, that is the first order Ordinary Differential Equation (ODE), y (t) = f(y(t)), t ∈ [t 0 , t 0 + T], associated with the initial condition: where f : IR m → IR m , m ≥ 1, is a C R−1 , R ≥ 1, function on its domain and y 0 ∈ IR m is assigned.Note that there is no loss of generality in assuming that the equation is autonomous.In this context, here, we focus on one-step Hermite-Obreshkov (HO) methods ( [1], p. 277).Unlike Runge-Kutta schemes, a high order of convergence is obtained with HO methods without adding stages.Clearly, there is a price for this because total derivatives of the f function are involved in the difference equation defining the method, and thus, a suitable smoothness requirement for f is necessary.Multiderivative methods have been considered often in the past for the numerical treatment of ODEs, for example also in the context of boundary value methods [2], and in the last years, there has been a renewed interest in this topic, also considering its application to the numerical solution of differential algebraic equations; see, e.g., [3][4][5][6][7][8].Here, we consider the numerical solution of Hamiltonian problems which in canonical form can be written as follows: with: where q and p are the generalized coordinates and momenta, H : IR 2 → IR is the Hamiltonian function and I stands for the identity matrix of dimension .Note that the flow ϕ t : y 0 → y(t) associated with the dynamical system (3) is symplectic; this means that its Jacobian satisfies: ∂ϕ t (y) ∂y J ∂ϕ t (y) ∂y = J, ∀ y ∈ IR 2 . ( A one-step numerical method Φ h : IR 2 → IR 2 with stepsize h is symplectic if the discrete flow y n+1 = Φ h (y n ), n ≥ 0, satisfies: Recently, the class of Euler-Maclaurin HO methods for the solution of Hamiltonian problems has been analyzed in [9,10] where, despite the non-existence results of symplectic multiderivative methods shown in [11], the conjugate symplecticity of the methods was proven.Two numerical methods Φ h , Ψ h are conjugate to each other if there exists a global change of coordinate χ h , such that: with χ h = y + O(h) uniformly for y varying in a compact set and • denoting a composition operator [12].If one method is conjugate to a symplectic method is said to be conjugate symplectic, this is a less strong requirement than symplecticity, which allows the numerical solution to have the same long-time behavior of a symplectic method.Observe that the conjugate symplecticity here refers to a property of the discrete flow of the two numerical methods; this should be not confused with the group of conjugate symplectic matrices, the set of matrices M ∈ C 2 that satisfy M H J M = J, where H means Hermitian conjugate [13].
In this paper, we consider the symmetric one-step HO methods, which were analyzed in [14,15] in the context of spline applications.We call them BSHO methods, since they are connected to B-Splines, as we will show.BSHO methods have a formulation similar to that of the Euler-Maclaurin formulas, and the order two and four schemes of the two families are the same.As a new result, we prove that BSHO methods are conjugate symplectic schemes, as is the case for the Euler-Maclaurin methods [9,10], and so, both families are suited to the context of geometric integration.BSHO methods are also strictly related to BS methods [16,17], which are a class of linear multistep methods also based on B-splines suited for addressing boundary value problems formulated as first order differential problems.Note that also BS methods were firstly studied in [14,15], but at that time, they were discarded in favor of BSHO methods since; when used as initial value methods, they are not convergent.In [16,17], the same schemes have been studied as boundary value methods, and they have been recovered in particular in connection with boundary value problems.As for the BSHO methods, the discrete solution generated by a BS method can be easily extended to a continuous spline collocating the differential problem at the mesh points [18].The idea now is to rely on B-splines with multiple inner knots in order to derive one-step HO schemes.The inner knot multiplicity is strictly connected to the number of derivatives of f involved in the difference equations defining the method and consequently with the order of the method.The efficient approach introduced in [18] dealing with BS methods for the computation of the collocating spline extension is here extended to BSHO methods, working with multiple knots.Note that we adopt a reversed point of view with respect to [14,15] because we assume to have already available the numerical solution generated by the BSHO methods and to be interested in an efficient procedure for obtaining the B-spline coefficients of the associated spline.
The paper is organized as follows.In Section 2, one-step symmetric HO methods are introduced, focusing in particular on BSHO methods.Section 3 is devoted to proving that BSHO methods are conjugate symplectic methods.Then, Section 4 first shows how these methods can be revisited in the spline collocation context.Successively, an efficient procedure is introduced to compute the B-spline form of the collocating spline extension associated with the numerical solution produced by the R-th BSHO, and it is shown that its convergence order is equal to that of the numerical solution.Section 6 presents some numerical results related to Hamiltonian problems, comparing them with those generated by Euler-Maclaurin and Gauss-Runge-Kutta schemes of the same order.

One-Step Symmetric Hermite-Obreshkov Methods
Let t i , i = 0, . . ., N, be an assigned partition of the integration interval [t 0 , t 0 + T], and let us denote by u i an approximation of y(t i ).Any one-step symmetric Hermite-Obreshkov (HO) method can be written as follows, clearly setting u 0 := y 0 , where h n := t n+1 − t n and where u r , for j ≥ 1, denotes the total (j − 1)-th derivative of f with respect to t computed at u r , Note that u r ≈ y (j) (t r ), and on the basis of (1), the analytical computation of the j-th derivative y (j) involves a tensor of order j.For example, y (2) (t) = df dt (y(t)) = ∂f ∂y (y(t)) f(y(t)) (where ∂f ∂y becomes the Jacobian m × m matrix of f with respect to y when m > 1).As a consequence, it is u (2) r = ∂f ∂y (u r ) f(u r ).We observe that the definition in (8) implies that only u n+1 is unknown in (7), which in general is a nonlinear vector equation in IR m with respect to it.
For example, the one-step Euler-Maclaurin [1] formulas of order 2s with s ∈ IN, s ≥ 1, (where the b 2i denote the Bernoulli numbers, which are reported in Table 2) belong to this class of methods.These methods will be referred to in the following with the label EMHO (Euler-Maclaurin Hermite-Obreshkov).
Here, we consider another class of symmetric HO methods that can be obtained by defining as follows the polynomial P 2R , appearing in ([1], Lemma 13.3), the statement of which is reported in Lemma 1.
Lemma 1.Let R be any positive integer and P 2R be a polynomial of exact degree 2R.Then, the following one-step linear difference equation, defines a multiderivative method of order 2R.
Referring to the methods obtainable by Lemma 1, if in particular the polynomial P 2R is defined as in (10), then we obtain the class of methods in which we are interested here.They can be written as in (7) with, which are reported in Table 1, for R = 1, . . ., 5. In particular, for R = 1 and R = 2, we obtain the trapezoidal rule and the Euler-Maclaurin method of order four, respectively.These methods were originally introduced in the spline collocation context, dealing in particular with splines with multiple knots [14,15], as we will show in Section 4. We call them BSHO methods since we will show that they can be obtained dealing in particular with the standard B-spline basis.The stability function of the R-th one-step symmetric BSHO method is the rational function corresponding to the (R, R)-Padé approximation of the exponential function, as is that of the same order Runge-Kutta-Gauss method ( [19], p. 72).It has been proven that methods with this stability function are A-stable ( [19], Theorem 4.12).For the proof of the statement of the following corollary, which will be useful in the sequel, we refer to [15] and {u i } N i=0 denotes the related numerical solution produced by the R-th one-step symmetric BSHO method in (7)- (11), it is: u

Conjugate Symplecticity of the Symmetric One-Step BSHO Methods
Following the lines of the proof given in [10], in this section, we prove that one-step symmetric BSHO methods are conjugate symplectic schemes.The following lemma, proved in [20], is the starting point of the proof, and it makes use of the B-series integrator concept.On this concern, referring to [12] for the details, here, we just recall that a B-series integrator is a numerical method that can be expressed as a formal B-series, that is it has a power series in the time step in which each term is a sum of elementary differentials of the vector field and where the number of terms is allowed to be infinite.Lemma 2. Assume that Problem (1) admits a quadratic first integral Q(y) = y T Sy (with S denoting a constant symmetric matrix) and that it is solved by a B-series integrator Φ h (y).Then, the following properties, where all formulas have to be interpreted in the sense of formal series, are equivalent: We observe that Lemma 2 is used in [21] to prove the conjugate symplecticity of symmetric linear multistep methods.With similar arguments, we prove the following theorem.
Theorem 1.The map u 1 = Φ h (u 0 ) associated with the one-step method ( 7)-( 11) admits a B-series expansion and satisfies Property (a) of Lemma 2.
Proof.By defining the two characteristic polynomials of the trapezoidal rule: and the shift operator E(u n ) := u n+1 , the R-th method described in (7) reads, We now consider a function v(t), a stepsize h and the shift operator E h (v(t)) := v(t + h), and we look for a continuous function v(t) that satisfies (12) in the sense of formal series (a series where the number of terms is allowed to be infinite), using the relation By multiplying both side of the previous equation by Dρ(e hD ) −1 , we obtain: that is, Now, since Bernoulli numbers define the Taylor expansion of the function z/(e z − 1) and b 0 = 1, b 1 = −1/2 and b j = 0 for the other odd j, we have: Thus, we can write that: With some algebra, we arrive at the following relation, with: Observe that δ j = 0 for j = 1, . . ., R − 1, since the method is of order 2R (see [12], Theorem 3.1, page 340).Therefore, we derive the modified initial value differential equation associated with the numerical scheme by coupling (15) with the initial condition v(t 0 ) = y 0 .Thus, the one-step symmetric BSHO methods are B-series integrators.The proof of Lemma 2 Property (a) follows exactly the same steps of the analogous proof in Theorem 1 of [10] and in [12] (Theorem 4.10, page 591).
In Table 2, we report the coefficients δ R for R ≤ 5 and the corresponding Bernoulli numbers.We can observe that the truncation error in the modified initial value problem is smaller than the one of the EMHO methods of the same order, which is equal to b i /i! (see [10]).The conjugate symplecticity property of a numerical scheme makes it suitable for the solution of Hamiltonian problems, since a conjugate symplectic method has the same long-time behavior of a symplectic one.A well-known pair of conjugate symplectic methods is composed by the trapezoidal and midpoint rules.Observe that the trapezoidal rule belongs to both the classes BSHO and EMHO of multiderivative methods, and its characteristic polynomial plays an important role in the proof of their conjugate symplecticity.

The Spline Extension
A (vector) Hermite polynomial of degree 2R + 1 interpolating both u n and u n+1 respectively at t n and t n+1 together with assigned derivatives u n+1 , k = 1, . . ., R, can be computed using the Newton interpolation formulas with multiple nodes.On the other hand, in his Ph.D. thesis [15], Loscalzo proved that a polynomial of degree 2R verifying the same conditions exists if and only if (7) is fulfilled with the β coefficients defined as in (11).Note that, since the polynomial of degree 2R + 1 fulfilling these conditions is always unique and its principal coefficient is given by the generalized divided difference u[t n , . . ., t n , t n+1 , . . ., t n+1 ] of order 2R + 1 associated with the given R-order Hermite data, the n-th condition in (7) holds iff this coefficient vanishes.If all the conditions in ( 7) are fulfilled, it is possible to define a piecewise polynomial, the restriction to [t n , t n+1 ] of which coincides with this polynomial, and it is clearly a C R spline of degree 2R with breakpoints at the mesh points.Now, when the definition given in ( 8) is used together with the assumption u 0 = y 0 , the conditions in ( 7) become a multiderivative one-step scheme for the numerical solution of (1).Thus, the numerical solution u n , n = 0, . . ., N it produces and the associated derivative values defined as in ( 8) can be associated with the above-mentioned 2R degree spline extension.Such a spline collocates the differential equation at the mesh points with multiplicity R, that is it verifies the given differential equation and also the equations y (j) (t) = d (j−1) (f•y) dt j−1 (t), j = 2, . . ., R at the mesh points.This piecewise representation of the spline is that adopted in [15].Here, we are interested in deriving its more compact B-spline representation.Besides being more compact, this also allows us to clarify the connection between BSHO and BS methods previously introduced in [16][17][18].For this aim, let us introduce some necessary notation.Let S 2R , be the space of C R 2R-degree splines with breakpoints at t i , i = 0, . . ., N, where t 0 < • • • < t N = t 0 + T. Since we relate to the B-spline basis, we need to introduce the associated extended knot vector: where: which means that all the inner breakpoints have multiplicity R in T and both t 0 and t N have multiplicity 2R + 1.The associated B-spline basis is denoted as B i , i = −2R, . . ., (N − 1)R and the dimension of S 2R as D, with D := (N + 1)R + 1.
The mentioned result proven by Loscalzo is equivalent to saying that, if the β coefficients are defined as in (11), any C R spline of degree 2R with breakpoints at the mesh points fulfills the relation in (7), where u (j) n denotes the j-th spline derivative at t n .In turn, this is equivalent to saying that such a relation holds for any element of the B-spline basis of S 2R .Thus, setting α := (−1 ; 1) T ∈ IR 2  and β (i) considering the local support of the B-spline basis, we have that (α; β (1) ; ...; β (R) ), where the punctuation mark ";" means vertical catenation (to make a column-vector), can be also characterized as the unique solution of the following linear system, G (n) (α; β (1) ; . . .; where e 2R+2 = (0; . . .; 0; 1) ∈ IR 2R+2 and: R+1 defined as, (19) where B (j) i denotes the j-th derivative of B i .Note that the last equation in (17), 2β In order to prove the non-singularity of the matrix G (n) , we need to introduce the following definition, Definition 1.Given a non-decreasing set of abscissas Θ := {θ i } M i=0 , we say that a function g 1 agrees with another function g 2 at Θ if g (j) 2 (θ i ), j = 0, . . ., m i − 1, i = 0, . . ., M, where m i denotes the multiplicity of θ i in Θ.
Proof.Observe that the restriction to I n = [t n , t n+1 ] of the splines in S 2R generates Π 2R since there are no inner knots in I n .Then, restricting to I n , Π 2R can be also generated by the B-splines of S 2R not vanishing in I n , that is from B (n−2)R , . . ., B nR .Since the polynomial in Π 2R agreeing with a given function in: is unique, it follows that also the corresponding (2R + 1) × (2R + 1) matrix collocating the spline basis active in I n is nonsingular.Such a matrix is the principal submatrix of G (n)T of order 2R + 1.Thus now, considering that the restriction to I n of any function in S 2R is a polynomial of degree 2R, we prove by reductio ad absurdum that the last row of G (n) cannot be a linear combination of the other rows.In fact, in the opposite case, there would exist a polynomial P of degree 2R such that P(t n ) = P(t n+1 ) = 0, P (t n ) = P (t n+1 ) = −1, and P (j) (t n ) = P (j) (t n+1 ) = 0, j = 2, . . ., R. Considering the specific interpolation conditions, this P does not fulfill the n-th condition in (7).This is absurd, since Loscalzo [15] has proven that such a condition is equivalent to requiring degree reduction for the unique polynomial of degree less than or equal to 2R + 1, fulfilling R + 1 Hermite conditions at both t n and t n+1 .
Note that this different form for defining the coefficient of the R-th BSHO scheme is analogous to that adopted in [17] for defining a BS method on a general partition.However, in this case, the coefficients of the scheme do not depend on the mesh distribution, so there is no need to determine them solving the above linear system.On the other hand, having proven that the matrix G (n) is nonsingular will be useful in the following for determining the B-spline form of the associated spline extension.
Thus, let us now see how the B-spline coefficients of the spline in S 2R associated with the numerical solution generated by the R-th BSHO can be efficiently obtained, considering that the following conditions have to be imposed, Now, we are interested in deriving the B-spline coefficients c i , i = −2R, . . ., (N − 1)R, of s 2R , Relying on the representation in (21), all the conditions in (20) can be re-written in the following compact matrix form, where c = (c −2R ; . . .; c (N−1)R ) ∈ IR mD , with c j ∈ IR m , I m is the identity matrix of size m × m, D is the dimension of the spline space previously introduced and where: A := (A 1 ; A 2 ; . . .; A R+1 ) , with each A being a (R + 1)-banded matrix of size (N + 1) × D (see Figure 1) with entries defined as follows: The following theorem related to the rectangular linear system in (22) ensures that the collocating spline s 2R is well defined.
Theorem 2. The rectangular linear system in (22) has always a unique solution, if the entries of the vector on its right-hand side satisfy the conditions in (7) with the β coefficients given in (11).
Proof.The proof is analogous to that in [18] (Theorem 1), and it is omitted.
We now move to introduce the strategy adopted for an efficient computation of the B-spline coefficients of s 2R .

Efficient Spline Computation
Concerning the computation of the spline coefficient vectors: the unique solution of ( 22) can be computed with several different strategies, which can have very different computational costs and can produce results with different accuracy when implemented in finite arithmetic.Here, we follow the local strategy used in [18].Taking into account the banded structure of A i , i = 1, . . ., R + 1, we can verify that ( 22) implies the following relations, where u = (u 0 ; . . .; u N ), c (i) := (c (i−3)R ; . . .; c (i−1)R ) ∈ IR m (2R+1) , i = 1, . . ., N and: As a consequence, we can also write that, where ĉ(i) := (c (i) ; 0) ∈ IR m (2R+2) .Now, for all integers r < 2R + 2, we can define other R + 1 auxiliary vectors α(R) i,r , β(R) l,i,r , l = 1, . . ., R ∈ IR 2 , defined as the solution of the following linear system, where e r is the r-th unit vector in IR 2R+2 (that is the auxiliary vectors define the r-th column of the inverse of G (i) ).Then, we can write, From this formula, considering (25), we can conclude that: Thus, solving all the systems (26) for i = 1, . . ., N, r = r 1 (i), . . ., r 2 (i), with: all the spline coefficients are obtained.Note that, with this approach, we solve D auxiliary systems, the size of which does not depend on N, using only N different coefficient matrices.Furthermore, only the information at t i−1 and t i is necessary to compute c (i−3)R+r−1 .Thus, the spline can be dynamically computed at the same time the numerical solution is advanced at a new time value.This is clearly of interest for a dynamical adaptation of the stepsize.
In the following subsection, relying on its B-spline representation, we prove that the convergence order of s 2R to y is equal to that of the numerical solution.This result was already available in [15] (see Theorem 4.2 in the reference), but proven with different longer arguments.

Spline Convergence
Let us assume the following quasi-uniformity requirement for the mesh, where M l and M u are positive constants not depending on h, with M l ≤ 1 and M u ≥ 1.Note that this requirement is a standard assumption in the refinement strategies of numerical methods for ODEs.We first prove the following result, that will be useful in the sequel.
Proposition 2. If y ∈ S 2R and so in particular if y is a polynomial of degree at most 2R, then: where y n := y(t n ), y n := d j y d j t (t n ), j = 1, . . ., R, n = 0, . . ., N, and the spline extension s 2R coincides with y.
Proof.The result follows by considering that the divided difference vanishes and, as a consequence, the local truncation error of the methods is null.
Then, we can prove the following theorem (where for notational simplicity, we restrict to m = 1), the statement of which is analogous to that on the convergence of the spline extension associated with BS methods [18].In the proof of the theorem, we relate to the quasi-interpolation approach for function approximation, the peculiarity of which consists of being a local approach.For example, in the spline context considered here, this means that only a local subset of a given discrete dataset is required to compute a B-spline coefficient of the approximant; refer to [22] for the details.Theorem 3. Let us assume that the assumptions on f done in Corollary 1 hold and that (28) holds.Then, the spline extension s 2R approximates the solution y of (1) with an error of order O(h 2R ) where h := max i=0,...,N−1 h i .Proof.Let s 2R denote the spline belonging to S 2R obtained by quasi-interpolating y with one of the rules introduced in Formula (5.1) in [22] by point evaluation functionals.From [22] (Theorem 5.2), under the quasi-uniformity assumption on the mesh distribution, we can derive that such a spline approximates y with maximal approximation order also with respect to all the derivatives, that is, where K is a constant depending only on R, M l and M u .
On the other hand, by using the triangular inequality, we can state that: Thus, we need to consider the first term on the right-hand side of this inequality.On this concern, because of the partition of unity property of the B-splines, we can write: where c := (c −2R ; . . .; c (N+1)R+1 ) and c := (c −2R ; . . .; c (N+1)R+1 ).
Now, for any function g ∈ C 2R [t 0 , t 0 + T], we can define the following linear functionals, where: and the vector ( α(R) i,r ; β(R) 1,i,r ; . . .; β(R) R,i,r ) has been defined in the previous section.Considering from Proposition 2 that s 2R , as well as any other spline belonging to S 2R can be written as follows, 29), we can deduce that: ) is defined in (26) as the r-th column of the inverse of the matrix G (i) .On the other hand, the entries of such nonsingular matrix do not depend on h, but because of the locality of the B-spline basis and of the R-th multiplicity of the inner knots, only on the ratios h j /h j+1 , j = i − 1, i, which are uniformly bounded from below and from above because of (28).
Thus, there exists a constant C depending on M l , M u and R such that G (i) −1 ≤ C, which implies that the same is true for any one of the mentioned coefficient vectors.From the latter, we deduce that for all indices, we find: On the other hand, taking into account the result reported in Corollary 1 besides (29), we can easily derive that w (i)

Approximation of the Derivatives
The computation of the derivative u (j) n , j ≥ 2, from the corresponding u n is quite expensive, and thus, usually, methods not requiring derivative values are preferred.Therefore, as well as for any other multiderivative method, it is of interest to associate with BSHO methods an efficient way to compute the derivative values at the mesh points.We are exploiting a number of possibilities, such as: • using generic symbolic tools, if the function f is known in closed form; • using a tool of automatic differentiation, like ADiGator, a MATLAB Automatic Differentiation Tool [23]; • using the Infinity Computer Arithmetic, if the function f is known as a black box [6,7,10]; • approximating it with, for example, finite differences.
As shown in the remainder of this section, when approximate derivatives are used, we obtain a different numerical solution, since the numerical scheme for its identification changes.In this case, the final formulation of the scheme is that of a standard linear multistep method, being still derived from (7) with coefficients in (11), but by replacing derivatives of order higher than one with their approximations.In this section, we just show the relation of these methods with a class of Boundary Value Methods (BVMs), the Extended Trapezoidal Rules (ETRs), linear multistep methods used with boundary conditions [24].Similar relations have been found in [25] with HO and the equivalent class of the super-implicit methods, which require the knowledge of functions not only at past, but also at future time steps.The ETRs can be derived from BSHO when the derivatives are approximated by finite differences.Let us consider the order four method with R = 2.In this case, the first derivative of f could be approximated using central differences: =: f i and u =: f i , is: after the approximation becomes: rearranging, we recover the ETR of order four: With similar arguments for the method of order six, R = 3, by approximating the derivatives with the order four finite differences: and: and rearranging, we obtain the sixth order ETR method: This relation allows us to derive a continuous extension of the ETR schemes using the continuous extension of the BSHO method, just substituting the derivatives by the corresponding approximations.Naturally, a change of the stepsize will now change the coefficients of the linear multistep schemes.
Observe that BVMs have been efficiently used for the solution of boundary value problems in [26], and the BS methods are also in this class [16].
It has been proven in [21] that symmetric linear multistep methods are conjugate symplectic schemes.Naturally, in the context of linear multistep methods used with only initial conditions, this property refers only to the trapezoidal method, but when we solve boundary value problems, the correct use of a linear multistep formula is with boundary conditions; this makes the corresponding formulas stable, with a region of stability equal to the left half plane of C (see [24]).The conjugate symplecticity of the methods is the reason for their good behavior shown in [27,28] when used in block form and with a sufficiently large block for the solution of conservative problems.
Remark 1.We recall that, even when approximated derivatives are used, the numerical solution admits a C R 2R-degree spline extension verifying all the conditions in (22), where all the u (j) n , j ≥ 2 appearing on the right-hand side have to be replaced with the adopted approximations.The exact solution of the rectangular system in ( 22) is still possible, since (7) with coefficients in (11) is still verified by the numerical solution u n , n = 0, . . ., N, by its derivatives u (1) n = f(u n ), n = 0, . . ., N and by the approximations of the higher order derivatives.The only difference in this case is that the continuous spline extension collocates at the breakpoints of just the given first order differential equation.

Numerical Examples
The numerical examples reported here have two main purposes: the first is to show the good behavior of BSHO methods for Hamiltonian problems, showing both the linear growth of the error for long time computation and the conservation of the Hamiltonian.To this end, we compare the methods with the symplectic Gauss-Runge-Kutta methods and with the conjugate symplectic EMHO methods.On the other hand, we are interested in showing the convergence properties of the spline continuous extensions.Observe that the availability of a continuous extension of the same order of the method is an important property.In fact for high order methods, especially for superconvergent methods like the Gauss ones, it is very difficult to find a good continuous extension.The natural continuous extension of these methods does not keep the same order of accuracy, without adding extra stages [29].Observe also that a good continuous extension is an important tool, for example for the event location.
We report results of our experiments for BSHO methods of order six and eight.We recall that the order two BSHO method corresponds to the well-known trapezoidal rule, the property of conjugate symplecticity of which is well known (see for example [12]) and the continuous extension by the B-spline of which has been already developed in [18].The order four BSHO belongs also to the EMHO class, and it has been analyzed in detail in [10].

Kepler Problem
The first example is the classical Kepler problem, which describes the motion of two bodies subject to Newton's law of gravitation.This problem is a completely integrable Hamiltonian nonlinear dynamical system with two degrees of freedom (see, for details, [30]).The Hamiltonian function: describes the motion of the body that is not located in the origin of the coordinate systems.This motion is an ellipse in the q 1 -q 2 plane, the eccentricity e of which is set using as starting values: and with period µ := 2π.The first integrals of this problem are: the total energy H, the angular momentum: M(q 1 , q 2 , p 1 , p 2 ) := q 1 p 2 − q 2 p 1 .
Only three of the four first integrals are independent, so, for example, A 1 can be neglected.
As in [10], we set e = 0.6 and h = µ/200, and we integrate the problem over 10 3 periods.Setting y := (q 1 , q 2 , p 1 , p 2 ), the error y j − y 0 1 in the solution is computed at specific times fixed equal to multiples of the period, that is at t j = 2πj, with j = 1, 2, . . .; the errors in the invariants have been computed at the mesh points t n = πn, n = 1, 3, 5 . ... Figure 2 reports the obtained results for the sixth and eighth order BSHO (dotted line, BSHO6, BSHO8), the sixth order EMHO (solid lines, EMHO6) and the sixth and eighth order Gauss-Runge-Kutta (GRK) (dashed lines, GRK6, GRK8) methods.In the top-left picture, the absolute error of the numerical solution is shown; the top-right picture shows the error in the Hamiltonian function; the error in the angular momentum is drawn in the bottom-left picture, while the bottom-right picture concerns the error in the second component of the Lenz vector.As expected from a symplectic or a conjugate symplectic integrator, we can see a linear drift in the error y j − y 0 1 as the time increases (top left plot).As well as for the other considered methods, we can see that BSHO methods guarantee a near conservation of the Hamiltonian function, of the second component of the Lenz vector and of the angular momentum (other pictures).This latter quadratic invariant is precisely conserved (up to machine precision) by GRK methods due to their symplecticity property.We observe also that, as expected, the error for the BSHO6 method is 3  10 of the error of the EMHO6 method.
To check the convergence behavior of the continuous extensions, we integrated the problem over 10 periods starting with stepsize h = µ/N, N = 100.We computed a reference solution using the order eight method with a halved stepsize, and we computed the maximum absolute error on the doubled grid.The results are reported in Table 3 for the solution and the first derivative and clearly show that the continuous extension respects the theoretical order of convergence.

Non-Linear Pendulum Problem
As a second example, we consider the dynamics of a pendulum under the influence of gravity.This dynamics is usually described in terms of the angle q that the pendulum forms with its stable rest position: q + sin q = 0, where p = q is the angular velocity.The Hamiltonian function associated with ( 31) is: An initial condition (q 0 , p 0 ) such that |H(q 0 , p 0 )| < 1 gives rise to a periodic solution y(t) = (q(t), p(t)) corresponding to oscillations of the pendulum around the straight-down stationary position.In particular, starting at y 0 = (q 0 , 0) , the period of oscillation may be expressed in terms of the complete elliptical integral of the first kind as: For the experiments, we choose q 0 = π/2; thus, the period µ is equal to 7.416298709205487.We use the sixth and eighth order BSHO and GRK methods and the sixth order EMHO method with stepsize h = µ/20 to integrate the problem over 2 • 10 4 periods.Setting y = (q, p), again, the errors y j − y 0 in the solution are evaluated at times that are multiples of the period µ, that is for t j = µj, with j = 1, 2, . . .; the energy error H(y n ) − H(y 0 ) has been computed at the mesh points t n = 11hn, n = 1, 2, . ... Figure 3 reports the obtained results.In the left plot, we can see that, for all the considered methods, the error in the solution grows linearly as time increases.A near conservation of the energy function is observable in both pictures on the right.The amplitudes of the bounded oscillations are similar for both methods, confirming the good long-time behavior properties of BSHO methods for the problem at hand.To check the convergence behavior of the continuous extensions, we integrated the problem over 10 periods starting with stepsize h = µ/N, N = 10.We computed a reference solution using the order eight method with a halved stepsize, and we compute the maximum absolute error on the doubled grid.The results are reported in Table 4 for the solution and the first derivative and clearly show, also for this example, that the continuous extension respects the theoretical order of convergence.

Conclusions
In this paper, we have analyzed the BSHO schemes, a class of symmetric one-step multi-derivative methods firstly introduced in [14,15] for the numerical solution of the Cauchy problem.As a new result, we have proven that these are conjugate symplectic schemes, thus suited to the context of geometric integration.Moreover, an efficient approach for the computation of the B-spline form of the spline extending the numerical solution produced by any BSHO method has been presented.The spline associated with the R-th BSHO method collocates the differential equation at the mesh points with multiplicity R and approximates the solution of the considered differential problem with the same accuracy O(h 2R ) characterizing the numerical solution.The relation between BSHO schemes and symmetric linear multistep methods when the derivatives are approximated by finite differences has also been pointed out.
Future related work will consist in studying the possibility of associating with the BSHO schemes a dual quasi-interpolation approach, as already done dealing with the BS linear multistep methods in [16,18,31].

Figure 2 .
Figure 2. Kepler problem: results for the sixth (BSHO6, red dotted line) and eighth (BSHO8, purple dotted line) order BSHO methods, sixth order Euler-Maclaurin method (EMHO6, blue solid line) and sixth (Gauss-Runge-Kutta (GRK6), yellow dashed line) and eighth (GRK8-green dashed line) order Gauss methods.(Top-left) Absolute error of the numerical solution; (top-right) error in the Hamiltonian function; (bottom-left) error in the angular momentum; (bottom-right) error in the second component of the Lenz vector.

Figure 3 .
Figure 3. Nonlinear pendulum problem: results for the Hermite-Obreshkov method of order six and eight (BSHO6, red, and BSHO8, purple dotted lines), for the sixth order Euler-Maclaurin (EMHO6, blue solid line) and Gauss methods (GRK6, yellow, and GRK8, green dashed lines) applied to the pendulum problem.(Left) plot: absolute error of the numerical solution; (upper-right) and (bottom-right) plots: error in the Hamiltonian function for the sixth order and eighth order integrators, respectively.

, Corollary 1. Let
us assume that f ∈ C 2R+1 (D), where D := {y ∈ IR m | ∃t ∈ [t 0 , t 0 + T] such that y − y(t) 2 ≤ L b }, with L b > 0.Then, there exists a positive constant h b such that if max

Table 2 .
Coefficients of the modified differential equations and Bernoulli numbers.

Table 3 .
Kepler problem: maximum absolute error of the numerical solution and its derivative computed for 10 periods.