Isospectral Property of Hamiltonian Boundary Value Methods (HBVMs) and their blended implementationWork developed within the project “Numerical methods and software for differential equations”.

# Isospectral Property of Hamiltonian Boundary Value Methods (HBVMs) and their blended implementation††thanks: Work developed within the project “Numerical methods and software for differential equations”.

L. Brugnano Luigi Brugnano Dipartimento di Matematica, Università di Firenze, Viale Morgagni 67/A, 50134 Firenze (Italy).
22email: luigi.brugnano@unifi.itFelice Iavernaro Dipartimento di Matematica, Università di Bari, Via Orabona 4, 70125 Bari (Italy).
44email: felix@dm.uniba.itDonato Trigiante Dipartimento di Energetica, Università di Firenze, Via Lombroso 6/17, 50134 Firenze (Italy).
66email: trigiant@unifi.it
F. Iavernaro Luigi Brugnano Dipartimento di Matematica, Università di Firenze, Viale Morgagni 67/A, 50134 Firenze (Italy).
22email: luigi.brugnano@unifi.itFelice Iavernaro Dipartimento di Matematica, Università di Bari, Via Orabona 4, 70125 Bari (Italy).
44email: felix@dm.uniba.itDonato Trigiante Dipartimento di Energetica, Università di Firenze, Via Lombroso 6/17, 50134 Firenze (Italy).
66email: trigiant@unifi.it
D. Trigiante Luigi Brugnano Dipartimento di Matematica, Università di Firenze, Viale Morgagni 67/A, 50134 Firenze (Italy).
22email: luigi.brugnano@unifi.itFelice Iavernaro Dipartimento di Matematica, Università di Bari, Via Orabona 4, 70125 Bari (Italy).
44email: felix@dm.uniba.itDonato Trigiante Dipartimento di Energetica, Università di Firenze, Via Lombroso 6/17, 50134 Firenze (Italy).
66email: trigiant@unifi.it
February 6, 2010.
###### Abstract

One main issue, when numerically integrating autonomous Hamiltonian systems, is the long-term conservation of some of its invariants, among which the Hamiltonian function itself. Recently, a new class of methods, named Hamiltonian Boundary Value Methods (HBVMs) has been introduced and analysed BIT (), which are able to exactly preserve polynomial Hamiltonians of arbitrarily high degree. We here study a further property of such methods, namely that of having, when cast as Runge-Kutta methods, a matrix of the Butcher tableau with the same spectrum (apart the zero eigenvalues) as that of the corresponding Gauss-Legendre method, independently of the considered abscissae. Consequently, HBVMs are always perfectly -stable methods. Moreover, this allows their efficient blended implementation, for solving the generated discrete problems.

###### Keywords:
polynomial Hamiltonian energy preserving methods extended collocation methods Hamiltonian Boundary Value Methods HBVMs block Boundary Value Methods blended implicit methods Runge-Kutta methods blended iteration
###### Msc:
65P10 65L05 65L06 65L80 65H10
journal:

## 1 Introduction

Hamiltonian problems are of great interest in many fields of application, ranging from the macro-scale of celestial mechanics, to the micro-scale of molecular dynamics. They have been deeply studied, from the point of view of the mathematical analysis, since two centuries. Their numerical solution is a more recent field of investigation, which has led to define symplectic methods, i.e., the simplecticity of the discrete map, considering that, for the continuous flow, simplecticity implies the conservation of . However, the conservation of the Hamiltonian and simplecticity of the flow cannot be satisfied at the same time unless the integrator produces the exact solution (see (HLW, , page 379)). More recently, the conservation of energy has been approached by means of the definition of the discrete line integral, in a series of papers IP1 (); IP2 (); IT1 (); IT2 (); IT3 (), leading to the definition of Hamiltonian Boundary Value Methods (HBVMs) BIS (); BIT2 (); BIT (); BIT1 (), which are a class of methods able to preserve, for the discrete solution, polynomial Hamiltonians of arbitrarily high degree (and then, a practical conservation of any sufficiently differentiable Hamiltonian). In more details, in BIT (), HBVMs based on Lobatto nodes have been analysed, whereas in BIT1 () HBVMs based on Gauss-Legendre abscissae have been considered. In the last reference, it has been actually shown that both formulae are essentially equivalent to each other, in the sense that they share the same order (twice the number of fundamental stages) and stability properties, and both methods provide the very same numerical solution, when the number of the so called silent stages tends to infitiny.111Actually, they both provide the same numerical solution also when the Hamiltonian is a polynomial and the number of silent stages is high enough to ensure the conservation property of the Hamiltonian function itself (see (2)). In this paper this conclusion if further supported, since we prove that all such methods, when cast as Runge-Kutta methods, have the corresponding matrix of the tableau, whose nonzero eigenvalues coincide with those of the corresponding Gauss-Legendre formula (isospectral property of HBVMs).

This property can be used to define an efficient iteration for solving the discrete problems generated by the methods, via their blended implementation. Indeed, after posing HBVMs in block BVM form, they can be recast in the framework of blended implicit methods, which have been studied in a series of papers B00 (); BM02 (); BM04 (); BM07 (); BM08 (); BM09 (); BM09a (); BMM06 (); BT2 () (see also C. Magherini’s PhD Thesis M04 ()). The latter methods have been successfully implemented in the two computational codes BiM and BiMD BIMcode (); the latter code is also included in the current release of the “Test Set for IVP Solvers” TestSet ().

With this premise, the structure of the paper is the following: in Section 2 the basic facts about HBVMs are recalled; in Section 3 we state the main result of this paper, concerning the isospectral property; in Section 4 the discrete problem to be actually solved is defined; in Section 5 it is shown that a corresponding blended iteration can be devised for its efficient solution, which can be tuned by choosing a free parameter; in Section 6 the optimal choice of the free parameter is done, on the basis of the isospectral property of HBVM, by using a linear analysis of convergence; finally, in Section 7 a few concluding remarks are given.

## 2 Hamiltonian Boundary Value Methods

The arguments in this section are worked out starting from the arguments used in BIT (); BIT1 () to introduce and analyse HBVMs. We consider canonical Hamiltonian problems in the form

 ˙y=J∇H(y),y(t0)=y0∈R2m, (1)

where is a skew-symmetric constant matrix, and the Hamiltonian is assumed to be sufficiently differentiable. The key formula which HBVMs rely on, is the line integral and the related property of conservative vector fields:

 H(y1)−H(y0)=h∫10˙σ(t0+τh)T∇H(σ(t0+τh))dτ, (2)

for any , where is any smooth function such that

 σ(t0)=y0,σ(t0+h)=y1. (3)

Here we consider the case where is a polynomial of degree , yielding an approximation to the true solution in the time interval . The numerical approximation for the subsequent time-step, , is then defined by (3).

After introducing a set of distinct abscissae

 0

we set

 Yi=σ(t0+cih),i=1,…,s, (5)

so that may be thought of as an interpolation polynomial, interpolating the fundamental stages , , at the abscissae (4). We observe that, due to (3), also interpolates the initial condition .

###### Remark 1

Sometimes, the interpolation at is explicitly required. In such a case, the extra abscissa is formally added to (4). This is the case, for example, of a Lobatto distribution of the abscissae BIT ().

Let us consider the following expansions of and for :

 ˙σ(t0+τh)=s∑j=1γjPj(τ),σ(t0+τh)=y0+hs∑j=1γj∫τ0Pj(x)dx, (6)

where is a suitable basis of the vector space of polynomials of degree at most and the (vector) coefficients are to be determined. We shall consider an orthonormal basis of polynomials on the interval , i.e.:

 ∫10Pi(t)Pj(t)dt=δij,i,j=1,…,s, (7)

where is the Kronecker symbol, and has degree . Such a basis can be readily obtained as

 Pi(t)=√2i−1^Pi−1(t),i=1,…,s,

with the shifted Legendre polynomial, of degree , on the interval .

###### Remark 2

From the properties of shifted Legendre polynomials (see, e.g., AS () or the Appendix in BIT ()), one readily obtains that the polynomials satisfy the three-terms recurrence relation:

 P1(t) ≡ 1,P2(t)=√3(2t−1), Pj+2(t) = (2t−1)2j+1j+1√2j+32j+1Pj+1(t)−jj+1√2j+32j−1Pj(t),j≥1.

We shall also assume that is a polynomial, which implies that the integrand in (2) is also a polynomial so that the line integral can be exactly computed by means of a suitable quadrature formula. It is easy to observe that in general, due to the high degree of the integrand function, such quadrature formula cannot be solely based upon the available abscissae : one needs to introduce an additional set of abscissae , distinct from the nodes , in order to make the quadrature formula exact:

 ∫10˙σ(t0+τh)T∇H(σ(t0+τh))dτ= s∑i=1βi˙σ(t0+cih)T∇H(σ(t0+cih))+r∑i=1^βi˙σ(t0+^cih)T∇H(σ(t0+^cih)),

where , , and , , denote the weights of the quadrature formula corresponding to the abscissae , i.e.,

 βi = ∫10⎛⎝s∏j=1,j≠it−cjci−cj⎞⎠(r∏j=1t−^cjci−^cj)dt,i=1,…,s, (9) ^βi = ∫10(s∏j=1t−cj^ci−cj)⎛⎝r∏j=1,j≠it−^cj^ci−^cj⎞⎠dt,i=1,…,r.
###### Remark 3

In the case considered in the previous Remark 1, i.e. when is formally added to the abscissae (4), the first product in each formula in (9) ranges from to . Moreover, also the range of becomes . However, for sake of simplicity, we shall not consider this case further.

According to IT2 (), the right-hand side of (2) is called discrete line integral, while the vectors

 ^Yi≡σ(t0+^cih),i=1,…,r, (10)

are called silent stages: they just serve to increase, as much as one likes, the degree of precision of the quadrature formula, but they are not to be regarded as unknowns since, from (6) and (10), they can be expressed in terms of linear combinations of the fundamental stages (5).

###### Definition 1

The method defined by substituting the quantities in (6) into the right-hand side of (2), and by choosing the unknown coefficients in order that the resulting expression vanishes, is called Hamiltonian Boundary Value Method with steps and degree , in short HBVM(,), where BIT ().

In the sequel, we shall see that HBVMs may be expressed through different, though equivalent, formulations: some of them can be directly implemented in a computer program, the others being of more theoretical interest.

Because of the equality (2), we can apply the procedure directly to the original line integral appearing in the left-hand side. With this premise, by considering the first expansion in (6), the conservation property reads

 s∑j=1γTj∫10Pj(τ)∇H(σ(t0+τh))dτ=0,

which, as is easily checked, is certainly satisfied if we impose the following set of orthogonality conditions:

 γj=∫10Pj(τ)J∇H(σ(t0+τh))dτ,j=1,…,s. (11)

Then, from the second relation of (6) we obtain, by introducing the operator

 L(f;h)σ(t0+ch)= σ(t0)+hs∑j=1∫c0Pj(x)dx∫10Pj(τ)f(σ(t0+τh))dτ,c∈[0,1],

that is the eigenfunction of relative to the eigenvalue :

 σ=L(J∇H;h)σ. (13)
###### Definition 2

Equation (13) is the Master Functional Equation (MFE) defining  BIT1 ().

###### Remark 4

Some further details are in order to better elucidate the role of the MFE in devising our methods. First of all we observe that, by definition, the MFE intrinsically brings, with its polynomial solutions , the conservation property of the Hamiltonian function: indeed (13) is equivalent to (2) under the choice (11).

This means that, when searching for its solutions, one should always take care of the precise dimension of the polynomial vector space, say , is intended to belong to: the higher is , the higher must be the number of silent stages (and hence the number of steps ) to guarantee that (2) be satisfied. This explains the way the solutions of the MFE depends on .

It is also clear that, assuming the same kind of distribution for all the the nodes (see later), (2) will be satisfied starting from a suitable number of steps on. This implies that, for all , HBVM will define the same polynomial of degree , such that .

To practically compute , we set (see (5) and (6))

 Yi=σ(t0+cih)=y0+hs∑j=1aijγj,i=1,…,s, (14)

where

 aij=∫ci0Pj(x)dx,i,j=1,…,s.

Inserting (11) into (14) yields the final formulae which define the HBVMs class based upon the orthonormal basis :

 Yi=y0+hs∑j=1aij∫10Pj(τ)J∇H(σ(t0+τh))dτ,i=1,…,s. (15)

We recall once again that we are working under the assumption (2), namely that the Hamiltonian is a polynomial and that we are considering a sufficient number of additional abscissae such that the line integral and its discrete counterpart do coincide. This implies that we can replace the integrals appearing in (15) by sums representing the associated quadrature formulae introduced in (2), without introducing any discretization error.

This leads back to express the HBVM method in terms of the fundamental stages and the silent stages (see (10)). By using the notation

 f(y)=J∇H(y), (16)

we obtain

 Yi=y0+hs∑j=1aij(s∑l=1βlPj(cl)f(Yl)+r∑l=1^βlPj(^cl)f(ˆYl)),i=1,…,s. (17)

We again stress that the silent stages may be removed from (17) by observing that, for example,

 ^Yl=s∑i=1ℓ(t0+^cih)Yi,l=1,…,r,

where are the cardinal Lagrange polynomials defined on the nodes , .

From the above discussion it is clear that formulae (15) also make sense in the non-polynomial case. In fact, supposing to choose the abscissae so that the sums in (17) converge to an integral as , the resulting formula is again (15).

###### Definition 3

Formula (15) is named -HBVM of degree or HBVMBIT1 ().

This implies that HBVMs may be as well applied in the non-polynomial case since, in finite precision arithmetic, HBVMs are undistinguishable from their limit formulae (15), when a sufficient number of silent stages is introduced. The aspect of having a practical exact integral, for large enough, was already stressed in BIS (); BIT (); BIT1 (); IP1 (); IT2 ().

### 2.1 Runge-Kutta formulation of HBVMs

On the other hand, we emphasize that, in the non-polynomial case, (15) becomes an operative method only after that a suitable strategy to approximate the integrals appearing in it is taken into account. In the present case, if one discretizes the Master Functional Equation (2)–(13), HBVM are then obtained, essentially by extending the discrete problem (17) also to the silent stages (10). In order to simplify the exposition, we shall use (16) and introduce the following notation:

 {ti}={ci}∪{^ci}, {ωi}={βi}∪{^βi}, yi=σ(t0+tih), fi=f(σ(t0+tih)),i=1,…,k.

The discrete problem defining the HBVM then becomes,

 yi=y0+hs∑j=1∫ti0Pj(x)dxk∑ℓ=1ωℓPj(tℓ)fℓ,i=1,…,k. (18)

By introducing the vectors

 y=(yT1,…,yTk)T,e=(1,…,1)T∈Rk,

and the matrices

 Ω=diag(ω1,…,ωk),Is, Ps∈Rk×s, (19)

whose th entry are given by

 (Is)ij=∫ti0Pj(x)dx,(Ps)ij=Pj(ti), (20)

we can cast the set of equations (18) in vector form as

 y=e⊗y0+h(IsPTsΩ)⊗I2mf(y),

with an obvious meaning of . Consequently, the method can be seen as a Runge-Kutta method with the following Butcher tableau:

 t1⋮tkIsPTsΩω1… ωk (21)
###### Remark 5

We observe that, because of the use of an orthonormal basis, the role of the abscissae and of the silent abscissae is interchangeable, within the set . This is due to the fact that all the matrices , , and depend on all the abscissae , and not on a subset of them, and they are invariant with respect to the choice of the fundamental abscissae .

Hereafter, we shall consider a Gauss distribution of the abscissae , so that the resulting HBVM method BIT1 ():

• has order for all ;

• is symmetric and perfectly -stable (i.e., its stability region coincides with the left-half complex plane, BT ());

• reduces to the Gauss-Legendre method of order , when ;

• exactly preserves polynomial Hamiltonian functions of degree , provided that

 k≥νs2.

## 3 The Isospectral Property

We are now going to prove a further additional result, related to the matrix appearing in the Butcher tableau (21), corresponding to HBVM, i.e., the matrix

 A=IsPTsΩ∈Rk×k,k≥s, (22)

whose rank is . Consequently it has a -fold zero eigenvalue. In this section, we are going to discuss the location of the remaining eigenvalues of that matrix.

Before that, we state the following preliminary result, whose proof can be found in (HW, , page 79).

###### Lemma 1

The eigenvalues of the matrix

 Xs=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝12−ξ1ξ10⋱⋱⋱−ξs−1ξs−10⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (23)

with

 ξj=12√(2j+1)(2j−1),j≥1, (24)

coincide with those of the matrix in the Butcher tableau of the Gauss-Legendre method of order .

We also need the following preliminary result, whose proof derives from the properties of shifted-Legendre polynomials (see, e.g., AS () or the Appendix in BIT ()).

###### Lemma 2

With reference to the matrices in (19)–(20), one has

 Is=Ps+1^Xs,

where

 ^Xs=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝12−ξ1ξ10⋱⋱⋱−ξs−1ξs−10ξs⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

with the defined by (24).

The following result then holds true.

###### Theorem 1 (Isospectral Property of HBVMs)

For all and for any choice of the abscissae such that the quadrature defined by the weights is exact for polynomials of degree , the nonzero eigenvalues of the matrix in (22) coincide with those of the matrix of the Gauss-Legendre method of order .

###### Proof

For , the abscissae have to be the Gauss-Legendre nodes, so that HBVM reduces to the Gauss Legendre method of order , as outlined at the end of Section 2.

When , from the orthonormality of the basis, see (7), and considering that the quadrature with weights is exact for polynomials of degree (at least) , one easily obtains that (see (19)–(20))

 PTsΩPs+1=(Is 0),

since, for all  ,  and  :

 (PTsΩPs+1)ij=k∑ℓ=1ωℓPi(tℓ)Pj(tℓ)=∫10Pi(t)Pj(t)dt=δij.

By taking into account the result of Lemma 2, one then obtains:

 APs+1 = IsPTsΩPs+1=Is(Is 0)=Ps+1^Xs(Is 0)=Ps+1(^Xs 0) ≡ Ps+1⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝12−ξ10ξ10⋱⋮⋱⋱−ξs−1⋮ξs−100ξs0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

with the defined according to (24). Consequently, one obtains that the columns of constitute a basis of an invariant (right) subspace of matrix , so that the eigenvalues of are eigenvalues of . In more detail, the eigenvalues of are those of (see (23)) and the zero eigenvalue. Then, also in this case, the nonzero eigenvalues of coincide with those of , i.e., with the eigenvalues of the matrix defining the Gauss-Legendre method of order

## 4 Solving the discrete problem

We shall now consider some computational aspects concerning HBVM. In more details, we now show how its cost depends essentially on , rather than on , in the sense that the nonlinear system to be solved, for obtaining the discrete solution, has (block) dimension . This has been already shown in BIT (), but here we derive the same result in a slightly more compact way, which will allow us to easily introduce blended HBVMs in the next section.

In order to simplify the notation, we shall fix the fundamental stages at , since we have already seen that, due to the use of an orthonormal basis , they could be in principle chosen arbitrarily, among the abscissae . With this premise, we have, from (15),

 yi=y0+hs∑j=1aijk∑ℓ=1ωℓPj(tℓ)fℓ,i=1,…,s. (25)

This equation is now coupled with that defining the silent stages, i.e., from (6) and (10),

 yi=y0+s∑j=1γj∫ti0Pj(t)dt,i=s+1,…,k. (26)

Let us now partition the matrices in (19)–(20) into and , containing the entries defined by the fundamental abscissae and the silent abscissae, respectively. Similarly, we partition the vector into , containing the fundamental stages, and containing the silent stages and, accordingly, let and be the diagonal matrices containing the corresponding entries in matrix . Finally, let us define the vectors , , and .

Consequently, we can rewrite (25) and (26), as

 y1 = e⊗y0+hIs1(PTs1 PTs2)(Ω1Ω2)⊗I2m(f(y1)f(y2)), (27) y2 = u⊗y0+hIs2⊗I2mγ, (28)

 y1=e⊗y0+hIs1⊗I2mγ,

thus giving

 y2 = (u−Is2I−1s1e)⊗y0+Is2I−1s1⊗I2my1 (29) ≡ ^u⊗y0+A1⊗I2my1,

in place of (28), where, evidently, . By setting

 B1=Is1PTs1Ω1∈Rs×s,B2=Is1PTs2Ω2∈Rs×k−s, (30)

substitution of (29) into (27) then provides, at last, the system of (block) size to be actually solved:

 F(y1) ≡ y1−e⊗y0−h[B1⊗I2mf(y1)+ B2⊗I2mf(^u⊗y0+A1⊗I2my1)]=0.

By using the simplified Newton method for solving (4), and setting

 C=B1+B2A1∈Rs×s, (32)

one obtains the iteration:

 (Is⊗I2m−hC⊗J0)δ(n) = −F(y(n)1)≡ψ(n)1, (33) y(n+1)1 = y(n)1+δ(n),n=0,1,…,

where is the Jacobian of evaluated at . Because of the result of Theorem 1, the following property of matrix holds true.

###### Theorem 2

The eigenvalues of matrix in (32) coincide with those of matrix (23), i.e., with the eigenvalues of the matrix of the Butcher array of the Gauss-Legendre method of order .

Proof Assuming, as usual for simplicity, that the fundamental stages are the first ones, one has that the discrete problem

 y=(eu)⊗y0+hA⊗I2mf(y),

which defines the Runge-Kutta formulation of the method, is equivalent, by virtue of (27), (29), (30), to

 (IsOs×r−A1Ir)⊗I2m(y1y2)= (e^u)⊗y0+h(B1B2Or×sOr×r)⊗I2m([]cf(y1)f(y2)),

where, as usual, .222As observed in IP2 (); IT3 (), such formulation fits the framework of block BVMs. Consequently, the eigenvalues of the matrix defined in (22) coincide with those of the matrix pencil

 ( (IsOs×r−A1Ir), (B1B2Or×sOr×r) ),

that is

 μ∈σ(A)  ⇔  μ(IsOs×r−A1Ir)(uv)=(B1B2Or×sOr×r)(uv),

for some nonzero vector . By setting , one obtains the zero eigenvalues of the pencil. For the remaining (nonzero) ones, it must be , so that:

 μu=(B1u+B2v)=(B1u+B2A1u)=Cu  ⇔  μ∈σ(C). □
###### Remark 6

From the result of Theorem 2, it follows that the spectrum of doesn’t depend on the choice of the fundamental abscissae, within the nodes . On the contrary, its condition number does: the latter appears to be minimized when the fundamental abscissae are symmetrically distributed, and approximately evenly spaced, in the interval . As a practical “rule of thumb”, the following algorithm appears to be almost optimal:333We plan to investigate this aspect further, in a forthcoming paper.

1. let the abscissae be chosen according to a Gauss-Legendre distribution of nodes;

2. then, let us consider equidistributed nodes in , say ;

3. select, as the fundamental abscissae, those nodes, among the , which are the closest ones to the ;

4. define matrix in (32) accordingly.

Clearly, for the above algorithm to provide a unique solution (resulting in a symmetric choice of the fundamental abscissae), the difference has to be even which, however, can be easily accomplished.

In order to give evidence of the effectiveness of the above algorithm, in Figure 2 we plot the condition number of matrix , for , and . As one can see, the condition number of turns out to be nicely bounded, for increasing values of , which makes the implementation (that we are going to analyze in the next section) effective also when finite precision arithmetic is used. For comparison, in Figure 2 there is the same plot, obtained by fixing the fundamental abscissae as the first ones.

## 5 Blended HBVMs

The solution of problem (33) is now cast into the framework of blended implicit methods B00 (); BM02 (); BM04 (); BM07 (); BM08 (); BMM06 (); BT2 (); M04 () as below described. First of all, we observe that, since is nonsingular, we can recast problem (33) in the equivalent form

 γ(C−1⊗I2m−hIs⊗J0)δ(n)=−γC−1⊗I2mF(y(n)1)≡ψ(n)2, (34)

where is a free parameter to be chosen later. Let us now introduce the weight (matrix) function

 θ=Is⊗Φ−1,Φ=I2m−hγJ0∈R2m×2m, (35)

and the blended formulation of the system to be solved,

 Mδ(n) ≡ [θ(Is⊗I2m−hC⊗J0)+ (36) (I−θ)γ(C−1⊗I2m−hIs⊗J0)]δ(n) = θψ(n)1+(I−θ)ψ(n)2≡ψ(n).

The latter system has still the same solution as the previous ones, since it is obtained as the blending, with weights and , of the two equivalent forms (33) and (34). For iteratively solving (36), we use the corresponding blended iteration, formally given by B00 (); BM02 (); BM04 (); BM07 (); BM08 (); BMM06 (); BT2 (); M04 ():

 δ(n,ℓ+1)=δ(n,ℓ)−θ(Mδ(n,ℓ)−ψ(n)),ℓ=0,1,…. (37)
###### Remark 7

A nonlinear variant of the iteration (37) can be obtained, by setting

 y(n,ℓ+1)=y(n,ℓ)+δ(n,ℓ),ψ(n,ℓ)1=−F(y(n,ℓ)1),

and similarly defined, as:

 δ(n,ℓ+1)=δ(n,ℓ)−θ(Mδ(n,ℓ)−ψ(n,ℓ)),ℓ=0,1,…. (38)
###### Remark 8

We emphasize that, for actually performing the iteration (35)–(37), as well as (38), one has to factor only the matrix in (35), which has the same size as that of the continuous problem, due to the (block) diagonal structure of  .

We end this section by observing that the above iteration (37) depends on a free parameter . It will be chosen in order to optimize the convergence properties of the iteration, according to a linear analysis of convergence, which is sketched in the next section.

## 6 Linear analysis of convergence

The linear analysis of convergence for the iterations (37) is carried out by considering the usual scalar test equation (see, e.g., BM09 () and the references therein),

 y′=λy,R(λ)<0.

By setting, as usual , the two equivalent formulations (33) and (34) become, respectively (omitting, for sake of brevity, the upper index ),

 (Is−qC)δ=ψ1,γ(C−1−qIs)δ=ψ2.

Moreover,

 θ≡θ(q)=(1−γq)−1Is, (39)

and the blended iteration (37) becomes

 δ(ℓ+1)=(Is−θ(q)M(q))δ(ℓ)+θ(q)ψ(q), (40)

with

 M(q) = θ(q)(Is−qC)+(Is−θ(q))γ(C−1−qIs), (41) ψ(q) = θ(q)ψ1+(Is−θ(q))ψ2.

Consequently, the iteration will be convergent if and only if the spectral radius, say , of the iteration matrix,

 Z(q)=Is−θ(q)M(q), (42)

is less than 1. The set

 Γ={q∈C:ρ(q)<1}

is the region of convergence of the iteration. The iteration is said to be (see BM09 () for details):

• -convergent,   if ;

• -convergent,   if it is -convergent and, moreover,  ,  as  .

For the iteration (40) one verifies that (see (39), (41), and (42))

 Z(q)=q(1−γq)2C−1(C−γIs)2, (43)

which is the null matrix at and at . Consequently, the iteration will be -convergent (and, therefore, -convergent), provided that maximum amplification factor,

 ρ∗≡maxR(q)=0ρ(q) ≤1.

From (43) one has that 444Hereafter, will denote the spectrum of matrix .

 μ∈σ(C) ⇔ q(μ−γ)2μ(1−γq)2∈σ(Z(q)).

By taking into account that

 maxR(q)=0|q||(1−γq)2|=12γ,

one then obtains that

 ρ∗=maxμ∈σ(C)|μ−γ|22γ|μ|.

For the Gauss-Legendre methods (and, then, for any matrix having the same spectrum), it can be shown that BM02 (); BMM06 () the choice

 γ=|μmin|≡minμ∈σ(C)|μ|, (44)

minimizes , which turns out to be given by

 ρ∗=1−cosφmin <1,φmin=Arg(μmin). (45)

In Table 1, we list the optimal value of the parameter , along with the corresponding maximum amplification factor , for various values of , which confirm that the iteration (40) is -convergent.

###### Remark 9

We then conclude that the blended iteration (37) turns out to be -convergent for HBVM methods, for all and .

## 7 Conclusions

In this paper, computational aspects related to the efficient implementation of HBVM methods with steps and degree () (in short, HBVM), have been recast in the framework of blended implicit methods. In more details, we have seen that the discrete problem generated by HBVM amounts to a nonlinear system of (block) dimension . Its efficient solution can be obtained by considering the blended formulation of the discrete problem, for which an efficient diagonal splitting can be easily defined. Consequently, to implement the nonlinear iteration, only one matrix having the same size as that of the continuous problem has to be factored. The free parameter, on which the blended iteration depends on, can be easily chosen, because of the isospectral property of HBVM, resulting in an -convergent iteration. Last, but not least, also the conditioning of the discrete problem depends only on , and it appears to tend to a (nicely) bounded limit, as grows. We plan, in the future, to implement HBVMs in blended formulation (in short, blended HBVMs), in a computational code for numerically solving Hamiltonian problems.

## References

• (1) M. Abramovitz, I.A. Stegun. Handbook of Mathematical Functions. Dover, 1965.
• (2) L. Brugnano. Blended block BVMs (BVMs): A family of economical implicit methods for ODEs. J. Comput. Appl.Math. 116 (2000) 41–62.
• (3) L. Brugnano, F. Iavernaro, T. Susca. Hamiltonian BVMs (HBVMs): implementation details and applications. “Proceedings of ICNAAM 2009”, AIP Conf. Proc. 1168 (2009) 723–726.
• (4) L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian BVMs (HBVMs): a family of “drift-free” methods for integrating polynomial Hamiltonian systems. “Proceedings of ICNAAM 2009”, AIP Conf. Proc. 1168 (2009) 715–718.
• (5) L. Brugnano, F. Iavernaro, D. Trigiante. Analisys of Hamiltonian Boundary Value Methods (HBVMs): a class of energy-preserving Runge-Kutta methods for the numerical solution of polynomial Hamiltonian dynamical systems. BIT (2009), submitted. (arXiv:0909.5659)
• (6) L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian Boundary Value Methods (Energy Preserving Discrete Line Integral Methods). Jour. of Numer. Anal., Industr. and Appl. Math. (2009) submitted. (arXiv:0910.3621)
• (7) L. Brugnano, C. Magherini. Blended implementation of block implicit methods for ODEs. Appl. Numer. Math. 42 (2002) 29–45.
• (8) L. Brugnano, C. Magherini. The BiM code for the numerical solution of ODEs. J. Comput. Appl. Math. 164–165 (2004) 145–158.
• (9) L. Brugnano, C. Magherini. Blended implicit methods for solving ODE and DAE problems, and their extension for second order problems. J. Comput. Appl. Math. 205 (2007) 777–790.
• (10) L. Brugnano, C. Magherini. Blended General Linear Methods based on Generalized BDF. AIP Conf. Proc. 1048 (2008) 871–874.
• (11) L. Brugnano, C. Magherini. Recent Advances in Linear Analysis of Convergence for Splittings for Solving ODE problems. Appl. Numer. Math. 59 (2009) 542–557.
• (12) L. Brugnano, C. Magherini. Blended General Linear Methods based on Boundary Value Methods in the GBDF family. Journal of Numerical Analysis, Industrial and Applied Mathematics 4, 1-2 (2009) 23–40.
• (13) L  Brugnano, C. Magherini, F. Mugnai. Blended implicit methods for the numerical solution of DAE problems. J. Comput. Appl. Math. 189 (2006) 34–50.
• (14) L. Brugnano, D. Trigiante. Solving Differential Problems by Multistep Initial and Boundary Value Methods. Gordon and Breach, Amsterdam, 1998.
• (15) L. Brugnano, D. Trigiante. Block implicit methods for ODEs, in: D. Trigiante (Ed.), Recent Trends in Numerical Analysis. Nova Science Publ. Inc., New York, 2001, pp. 81–105.
• (16) E. Hairer, C. Lubich, G. Wanner. Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, 2 ed., Springer, Berlin, 2006.
• (17) E. Hairer, G. Wanner. Solving Ordinary Differential Equations II, 2 ed., Springer, Berlin, 1996.
• (18) F. Iavernaro, B. Pace. -Stage Trapezoidal Methods for the Conservation of Hamiltonian Functions of Polynomial Type. AIP Conf. Proc. 936 (2007) 603–606.
• (19) F. Iavernaro, B. Pace. Conservative Block-Boundary Value Methods for the Solution of Polynomial Hamiltonian Systems. AIP Conf. Proc. 1048 (2008) 888–891.
• (20) F. Iavernaro, D. Trigiante. Discrete conservative vector fields induced by the trapezoidal method. J. Numer. Anal. Ind. Appl. Math. 1 (2006) 113–130.
• (21) F. Iavernaro, D. Trigiante. State-dependent symplecticity and area preserving numerical methods. J. Comput. Appl. Math. 205 no. 2 (2007) 814–825.
• (22) F. Iavernaro, D. Trigiante. High-order symmetric schemes for the energy conservation of polynomial Hamiltonian problems. J. Numer. Anal. Ind. Appl. Math. 4,1-2 (2009) 87–101.
• (23) C. Magherini. Numerical Solution of Stiff ODE-IVPs via Blended Implicit Methods: Theory and Numerics. PhD thesis, Dipartimento di Matematica “U. Dini”, Università degli Studi di Firenze, September 2004 (Available at the url BIMcode ()).
• (24) Codes BiM/BiMD Homepage:  http://www.math.unifi.it/~brugnano/BiM/index.html
• (25) Test Set for IVP Solvers (release 2.4):  http://pitagora.dm.uniba.it/~testset/
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters