The Lack of Continuity and the Role of Infinite and Infinitesimal in Numerical Methods for ODEs: the Case of Symplecticity
aa r X i v : . [ m a t h . NA ] O c t The Lack of Continuity and the Role of Infinite andInfinitesimal in Numerical Methods for ODEs: the Caseof Symplecticity ✩ Luigi Brugnano a , Felice Iavernaro b , Donato Trigiante c a Dipartimento di Matematica, Universit`a di Firenze, Viale Morgagni 67/A, 50134 Firenze(Italy). b Dipartimento di Matematica, Universit`a di Bari, Via Orabona 4, 70125 Bari (Italy). c Dipartimento di Energetica, Universit`a di Firenze, Via Lombroso 6/17, 50134 Firenze(Italy).
Abstract
When numerically integrating canonical Hamiltonian systems, the long-termconservation of some of its invariants, among which the Hamiltonian functionitself, assumes a central role. The classical approach to this problem has led tothe definition of symplectic methods, among which we mention Gauss-Legendrecollocation formulae. Indeed, in the continuous setting, energy conservationis derived from symplecticity via an infinite number of infinitesimal contacttransformations . However, this infinite process cannot be directly transferred tothe discrete setting. By following a different approach, in this paper we describea sequence of methods, sharing the same essential spectrum (and, then, thesame essential properties), which are energy preserving starting from a certainelement of the sequence on, i.e., after a finite number of steps.
Keywords: polynomial Hamiltonian, energy preserving methods, symplecticmethods, Hamiltonian Boundary Value Methods, HBVMs, Runge-Kuttacollocation methods
1. Introduction
In order to make easier the following considerations, it is necessary to stressthe differences between a continuous problem, for example (1), and a discreteproblem obtained by applying to it a whatever numerical method. The maindifference between them, very often underrated by many authors, is the lack ✩ Work developed within the project “Numerical methods and software for differential equa-tions”.
Email addresses: [email protected] (Luigi Brugnano), [email protected] (Felice Iavernaro), [email protected] (Donato Trigiante)
Preprint submitted to Elsevier April 9, 2018 f continuity in time. The latter does not affect many aspects such as, e.g.,the study of the qualitative behavior of solutions around asymptotically stablecritical points (stability analysis). In fact, the respective theories, with minorchanges, are very similar (see, e.g., [27]). As a consequence, many tools, alreadydevised in the continuous analysis, can be transferred to the discrete analysisalmost unchanged. This is the case, for example, of the linearization aroundasymptotically stable critical points, which has been extensively used in thenumerical analysis of methods for differential problems (linear stability analy-sis). There are, however, other aspects for which continuity plays an essentialrole. For example, two solutions of the continuous problem need to stay awayfrom each other while, in the discrete case, they may interlace (without havingcommon points, of course). This fact has many mathematical and even physicalimplications (see, e.g., [26]). In this paper we shall deal with another case inwhich continuity plays an essential role. It regards the role of symplecticitywhich is central in discussing energy conservation in continuous Hamiltonianproblems, while it is less crucial in the energy conservation of discrete problems.This depends on the interplay between infinitesimal contact transformationsand the need of infinite processes (number of iterations) which cannot be opera-tively used in Numerical Analysis. This question has been already discussed inprevious papers (see, e.g., [7]). Here, after a rapid introduction to the subject,we shall focus on a particular aspect, although very important, which concernsa property of the numerical methods, in the general class of collocation meth-ods, which comes out by no more requiring symplecticity but still providingconservation of the Hamiltonian functions on a subset of points of the mesh.Clearly, this permits to avoid the drift of energy experienced when using manynumerical methods proposed in the recent literature.With this premise, the structure of the paper is the following: in Section 2 werecall the basic facts about canonical Hamiltonian problems and the approachesused for their numerical solution; one of them, resulting in the recently intro-duced class of
Hamiltonian Boundary Value Methods (HBVMs) , is then sketchedin Section 3; in Section 4 we state the main result of this paper, concerning the isospectral property of such methods; in Section 5 such property is further gen-eralized to study the existing connections between HBVMs and Runge-Kuttacollocation methods; a few concluding remarks are finally given in Section 6.
2. Canonical Hamiltonian problems
Canonical Hamiltonian problems are in the form˙ y = J ∇ H ( y ) , y ( t ) = y ∈ R m , (1)where J is a skew-symmetric constant matrix, the Hamiltonian H ( y ) is assumedto be sufficiently differentiable, and the state vector splits into two blocks, y =( q T , p T ) T , q, p ∈ R m where, for mechanical systems, q denote the positions and p the (generalized) momenta. Such problems are of great interest in many fields2f application, ranging from the macro-scale of celestial mechanics, to the micro-scale of molecular dynamics. They have been deeply studied, from the point ofview of the mathematical analysis, since two centuries. Their numerical solutionis a more recent field of investigation, where the main difficulty in dealing withthem numerically stems from the fact that the meaningful isolated critical pointsof such systems are only marginally stable: neighboring solution curves do noteventually approach the equilibrium point either in future or in past times. Thisimplies that the geometry around them critically depends on perturbations ofthe linear part. Consequently, the use of a linear test equation, which essentiallycaptures the geometry of the linear part, whose utility has been enormous insettling the dissipative case, cannot be of any utility in the present case.It is then natural to look for other properties of Hamiltonian systems that canbe imposed on the discrete methods in order to make them effective. The firstproperty which comes to mind is the symplecticity of the flow ϕ t := y y ( t )associated with (1). This property can be described either in geometric form(invariance of areas, volumes, etc.) or in analytical form: (cid:18) ∂ϕ t ∂y (cid:19) T J (cid:18) ∂ϕ t ∂y (cid:19) = J. In one way or the other, it essentially consists in moving infinitesimally on thetrajectories representing the solutions. Infinitesimally means retaining only thelinear part of the infinitesimal time displacement δt . It can be shown that thisproduces new values of the variables q + δq, p + δp which leave unchanged thevalue of the Hamiltonian H ( q + δq, p + δp ) = H ( q, p ) ( Infinitesimal ContactTransformation (ICT), see [13, p. 386]). Consequently, since the composition ofsuch infinitesimal transformations maintains the invariance, so does an infinitenumber of them.It is not surprising that the first numerical attempts to design conservativemethods have tried to transfer similar arguments to discrete methods, i.e. todesign symplectic integrators [24, 12] (see also the monographs [25, 21, 14] formore details on the subject). A backward error analysis has shown that sym-plecticity seems somehow to improve the long-time behavior properties of thenumerical solutions. Indeed, for a symplectic method of order r , implementedwith a constant stepsize h , the following estimation reveals how the numericalsolution y n may depart from the manifold H ( y ) = H ( y ) of the phase space,which contains the continuous solution itself: H ( y n ) − H ( y ) = O ( nh e − h h h r ) , (2)where h > is sufficiently small and h ≤ h . Relation (2) implies that a lineardrift of the energy with respect to time t = nh may arise. However, due to thepresence of the exponential, such a drift will not appear as far as nh ≤ e h h : thiscircumstance is often referred to by stating that symplectic methods conservethe energy function on exponentially long time intervals (see, for example, [14,Theorem 8.1, p. 367]). This is clearly a surrogate of the definition of stability in3hat the “good behaviour” of the numerical solution is not extended on infinite time intervals.As matter of fact, since symplecticity requires an infinite sequence of in-finitesimal contact transformations it cannot be transferred “sic et simpliciter” to the discrete methods, simply because infinite processes are not permitted inNumerical Analysis. A more efficient approach would require to design meth-ods which avoid the necessity of using infinite processes while preserving theconstant of motion, i.e., yielding a numerical solution belonging to the manifold H ( y ) = H ( y ). In this paper we consider a recently introduced class of methodsof any high order that provide energy conservation. More specifically, with anygiven order, one can associate an infinite sequences of methods, differing fromeach other for the number of internal mesh points that cover the same timewindow, say [ t i , t i + h ]: the more points we include, the better the conservationproperties of the method. However, we show that, at least for polynomial Hamil-tonian functions, such a process of increasing the number of internal points isnot infinite. In fact we show that there exists a finite value of new added pointsstarting from which the method become conservative, whatever is the stepsize h used.The evolution of the approaches to the problem, i.e. to get efficient energy-conserving methods, has been slow. As a matter of fact, the first unsuccessfulattempts to derive energy-preserving Runge-Kutta methods culminated in thewrong general feeling that such methods could not even exist [20]. One of thefirst successful attempts to solve the problem, outside the class of Runge-Kuttamethods, is represented by discrete gradient methods (see [22] and referencestherein) which are second order accurate. Purely algebraic approaches havebeen also introduced (see, e.g., [11]), without presenting any energy preserv-ing method. A further approach was considered in [23], where the averagedvector field method was proposed and shown to conserve the energy function,although it is only second order accurate. As was recently outlined in [10], ap-proximating the integral appearing in such method by means of a quadratureformula (based upon polynomial interpolation) yields a family of second orderRunge-Kutta methods. These latter methods represent an instance of energy-preserving Runge-Kutta methods for polynomial Hamiltonian problems: theirfirst appearance may be found in [16], under the name of s -stage trapezoidalmethods . Additional examples of fourth and sixth-order conservative Runge-Kutta methods (for polynomial Hamiltonians of suitable degree) were presentedin [17, 19]. All such energy-conserving methods have been derived by means ofthe new concept of discrete line integral .The evolution of this idea eventually led to the definition of HamiltonianBoundary Value Methods (HBVMs) [2, 3, 4], which is a wide class of methodsable to preserve, for the discrete solution, polynomial Hamiltonians of arbitrar-ily high degree (and then, a practical conservation of any sufficiently differen-tiable Hamiltonian). In more details, in [3] HBVMs defined at Lobatto nodeshave been analysed, whereas in [4] HBVMs defined at Gauss-Legendre abscis-sae have been considered. In the last reference, it has been actually shownthat both formulae are essentially equivalent to each other, since the order and4tability properties of the methods turn out to be independent of the abscissaedistribution, and both methods are equivalent, when the number of the so called silent stages tends to infinity . In this paper this conclusion is further supported,since we prove that HBVMs, when cast as Runge-Kutta methods, are such thatthe corresponding matrix of the tableau has the nonzero eigenvalues coincidentwith those of the corresponding Gauss-Legendre formula ( isospectral property of HBVMs). This property will be also used to further analyse the existingconnections between HBVMs and Runge-Kutta collocation methods.
3. Hamiltonian Boundary Value Methods
The arguments in this section are worked out starting from those used in [3,4] to introduce and analyse HBVMs. Starting from the canonical Hamiltonianproblems (1), the key formula which HBVMs rely on, is the line integral andthe related property of conservative vector fields: H ( y ) − H ( y ) = h Z ˙ σ ( t + τ h ) T ∇ H ( σ ( t + τ h ))d τ, (3)for any y ∈ R m , where σ is any smooth function such that σ ( t ) = y , σ ( t + h ) = y . (4)Here we consider the case where σ ( t ) is a polynomial of degree s , yielding anapproximation to the true solution y ( t ) in the time interval [ t , t + h ]. Thenumerical approximation for the subsequent time-step, y , is then defined by(4). After introducing a set of s distinct abscissae,0 < c , . . . , c s ≤ , (5)we set Y i = σ ( t + c i h ) , i = 1 , . . . , s, (6)so that σ ( t ) may be thought of as an interpolation polynomial, interpolatingthe fundamental stages Y i , i = 1 , . . . , s , at the abscissae (5). We observe that,due to (4), σ ( t ) also interpolates the initial condition y . Remark 1.
Sometimes, the interpolation at t is explicitly required. In such acase, the extra abscissa c = 0 is formally added to (5). This is the case, forexample, of a Lobatto distribution of the abscissae [3]. Let us consider the following expansions of ˙ σ ( t ) and σ ( t ) for t ∈ [ t , t + h ]:˙ σ ( t + τ h ) = s X j =1 γ j P j ( τ ) , σ ( t + τ h ) = y + h s X j =1 γ j Z τ P j ( x ) d x, (7)where { P j ( t ) } is a suitable basis of the vector space of polynomials of degreeat most s − { γ j } are to be determined. We5hall consider an orthonormal polynomial basis on the interval [0 ,
1] (though, inprinciple, different bases could be also considered [4, 16, 19, 17]): Z P i ( t ) P j ( t )d t = δ ij , i, j = 1 , . . . , s, (8)where δ ij is the Kronecker symbol, and P i ( t ) has degree i −
1. Such a basis canbe readily obtained as P i ( t ) = √ i − P i − ( t ) , i = 1 , . . . , s, (9)with ˆ P i − ( t ) the shifted Legendre polynomial, of degree i −
1, on the interval[0 ,
1] (see, e.g., [1]). We shall also assume that H ( y ) is a polynomial, whichimplies that the integrand in (3) is also a polynomial so that the line integralcan be exactly computed by means of a suitable quadrature formula. It is easyto observe that in general, due to the high degree of the integrand function,such quadrature formula cannot be solely based upon the available abscissae { c i } : one needs to introduce an additional set of abscissae { ˆ c , . . . , ˆ c r } , distinctfrom the nodes { c i } , in order to make the quadrature formula exact: Z ˙ σ ( t + τ h ) T ∇ H ( σ ( t + τ h ))d τ = (10) s X i =1 β i ˙ σ ( t + c i h ) T ∇ H ( σ ( t + c i h )) + r X i =1 ˆ β i ˙ σ ( t + ˆ c i h ) T ∇ H ( σ ( t + ˆ c i h )) , where β i , i = 1 , . . . , s , and ˆ β i , i = 1 , . . . , r , denote the weights of the quadratureformula defined at the abscissae { c i } ∪ { ˆ c i } . Then, according to [3, 4], we givethe following definition. Definition 1.
The method defined by the polynomial σ ( t ) , determined by sub-stituting the quantities in (7) into the right-hand side of (10) , and by choosingthe unknown coefficient { γ j } in order that the resulting expression vanishes, iscalled Hamiltonian Boundary Value Method with k steps and degree s , in short HBVM( k , s ) , where k = s + r . According to [18], the right-hand side of (10) is called discrete line integral associated with the map defined by the HBVM( k , s ) method, while the vectorsˆ Y i ≡ σ ( t + ˆ c i h ) , i = 1 , . . . , r, (11)are called silent stages : they just serve to increase, as much as one likes, thedegree of precision of the quadrature formula, but they are not to be regardedas unknowns since, from (7) and (11), they can be expressed in terms of linearcombinations of the fundamental stages (6).Because of the equality (10), we can apply the procedure described in Def-inition 1 directly to the original line integral appearing in the left-hand side.6ith this premise, by considering the first expansion in (7), the conservationproperty reads s X j =1 γ Tj Z P j ( τ ) ∇ H ( σ ( t + τ h ))d τ = 0 , which, as is easily checked, is satisfied if we impose the following set of orthog-onality conditions: γ j = Z P j ( τ ) J ∇ H ( σ ( t + τ h ))d τ, j = 1 , . . . , s. (12)Then, from the second relation of (7) we obtain, by introducing the operator L ( f ; h ) σ ( t + ch ) = (13) σ ( t ) + h s X j =1 Z c P j ( x )d x Z P j ( τ ) f ( σ ( t + τ h ))d τ, c ∈ [0 , , that σ is the eigenfunction of L ( J ∇ H ; h ) relative to the eigenvalue λ = 1: σ = L ( J ∇ H ; h ) σ. (14)According to [4], (14) is called the Master Functional Equation defining σ : itcharacterizes HBVM( k, s ) methods, for all k ≥ s . Indeed, such methods areuniquely defined by the polynomial σ , of degree s , the number of steps k beingonly required to obtain the exact quadrature formula (10).To practically compute σ , we set (see (6) and (7)) Y i = σ ( t + c i h ) = y + h s X j =1 a ij γ j , i = 1 , . . . , s, (15)where a ij = Z c i P j ( x )d x, i, j = 1 , . . . , s. Inserting (12) into (15) yields the final formulae which define the HBVMs classbased upon the orthonormal basis { P j } : Y i = y + h s X j =1 a ij Z P j ( τ ) J ∇ H ( σ ( t + τ h )) d τ, i = 1 , . . . , s. (16)For sake of completeness, we report the nonlinear system associated withthe HBVM( k, s ) method, in terms of the fundamental stages { Y i } and the silentstages { ˆ Y i } (see (11)), by using the notation f ( y ) = J ∇ H ( y ) . (17)7t represents the discrete counterpart of (16), which may be directly retrieved byevaluating, for example, the integrals in (16) by means of the (exact) quadratureformula introduced in (10): Y i = (18)= y + h s X j =1 a ij s X l =1 β l P j ( c l ) f ( Y l ) + r X l =1 ˆ β l P j (ˆ c l ) f ( b Y l ) ! , i = 1 , . . . , s. From the above discussion it is clear that, in the non-polynomial case, supposingto choose the abscissae { ˆ c i } so that the sums in (18) converge to an integral as r ≡ k − s tends to infinity , the resulting formula is (16), which has been named ∞ -HBVM of degree s or HBVM ( ∞ , s ) in [4]. This implies that HBVMs may beas well applied in the non-polynomial case since, in finite precision arithmetic,HBVMs are undistinguishable from their limit formulae (16), when a sufficientnumber of silent stages is introduced, so that a practical energy conservation isobtained, for k large enough [3, 4, 16, 19]. On the other hand, we emphasizethat, in the non-polynomial case, (16) becomes an operative method only afterthat a suitable strategy to approximate the integrals appearing in it is taken intoaccount. In the present case, if one discretizes the Master Functional Equation (13)–(14), HBVM( k, s ) are then obtained, essentially by extending the discreteproblem (18) also to the silent stages (11). In more details, by using (17) andintroducing the following notation: { τ i } = { c i } ∪ { ˆ c i } , { ω i } = { β i } ∪ { ˆ β i } ,y i = σ ( t + τ i h ) , f i = f ( σ ( t + τ i h )) , i = 1 , . . . , k, the discrete problem defining the HBVM( k, s ) method becomes, y i = y + h s X j =1 Z τ i P j ( x )d x k X ℓ =1 ω ℓ P j ( τ ℓ ) f ℓ , i = 1 , . . . , k. (19)By defining the vectors y = ( y T , . . . , y Tk ) T and e = (1 , . . . , T ∈ R k , and thematrices Ω = diag( ω , . . . , ω k ) , I s , P s ∈ R k × s , (20)whose ( i, j )th entry are given by( I s ) ij = Z τ i P j ( x )d x, ( P s ) ij = P j ( τ i ) , (21)we can cast the set of equations (19) in vector form as y = e ⊗ y + h ( I s P Ts Ω) ⊗ I m f ( y ) , f ( y ). Consequently, the method can be regardedas a Runge-Kutta method with the following Butcher tableau: τ ... τ k I s P Ts Ω ω . . . ω k (22)In particular, when a Gauss distribution of the abscissae { τ , . . . , τ k } is con-sidered, it can be proved that the resulting HBVM( k, s ) method [4] (see also[3, 5]): • has order 2 s for all k ≥ s ; • is symmetric and perfectly A -stable (i.e., its stability region coincides withthe left-half complex plane, C − [9]); • reduces to the Gauss-Legendre method of order 2 s , when k = s ; • exactly preserves polynomial Hamiltonian functions of degree ν , providedthat k ≥ νs . (23)
4. The Isospectral Property
We are now going to prove a further additional result, related to the matrixappearing in the Butcher tableau (22), i.e., the matrix A = I s P Ts Ω ∈ R k × k , k ≥ s, (24)whose rank is s . Consequently it has a ( k − s )-fold zero eigenvalue. In thissection, we are going to discuss its essential spectrum , i.e., the location of theremaining s nonzero eigenvalues of that matrix. Before that, we state a coupleof preliminary results: their proofs follow, respectively, from [15, Theorem 5.6,p. 83] and from the properties of shifted Legendre polynomials (see, e.g., [1] orthe Appendix in [3]). Lemma 1.
The eigenvalues of the matrix X s = − ξ ξ . . .. . . . . . − ξ s − ξ s − , (25) with ξ j = 12 p (2 j + 1)(2 j − , j ≥ , (26) coincide with those of the matrix in the Butcher tableau of the Gauss-Legendremethod of order s . emma 2. With reference to the matrices in (20)–(21), one has I s = P s +1 ˆ X s , where ˆ X s = − ξ ξ . . .. . . . . . − ξ s − ξ s − ξ s , with the ξ j defined by (26). The following result then holds true.
Theorem 1 (Isospectral Property of HBVMs).
For all k ≥ s and for anychoice of the abscissae { τ i } such that the quadrature defined by the weights { ω i } is exact for polynomials of degree s − , the nonzero eigenvalues of the matrix A in (24) coincide with those of matrix (25), characterizing the Gauss-Legendremethod of order s . Proof For k = s , the abscissae { τ i } have to be the s Gauss-Legendre nodes, sothat HBVM( s, s ) reduces to the Gauss Legendre method of order 2 s , as alreadyoutlined at the end of Section 3. When k > s , from the orthonormality of thebasis, see (8), and considering that the quadrature with weights { ω i } is exactfor polynomials of degree (at least) 2 s −
1, one obtains that (see (20)–(21)) forall i = 1 , . . . , s and j = 1 , . . . , s + 1, (cid:0) P Ts Ω P s +1 (cid:1) ij = k X ℓ =1 ω ℓ P i ( t ℓ ) P j ( t ℓ ) = Z P i ( t ) P j ( t )d t = δ ij , and, therefore, P Ts Ω P s +1 = ( I s ) . By taking into account the result of Lemma 2, one then obtains: A P s +1 = I s P Ts Ω P s +1 = I s ( I s ) = P s +1 ˆ X s ( I s ) = P s +1 (cid:16) ˆ X s (cid:17) = P s +1 − ξ ξ − ξ s − ... ξ s − ξ s ≡ P s +1 e X s , (27)with the { ξ j } defined according to (26). Consequently, one obtains that thecolumns of P s +1 constitute a basis of an invariant (right) subspace of matrix10 , so that the eigenvalues of e X s are eigenvalues of A . In more detail, theeigenvalues of e X s are those of X s (see (25)) and the zero eigenvalue. Then,also in this case, the nonzero eigenvalues of A coincide with those of X s , i.e.,with the eigenvalues of the matrix defining the Gauss-Legendre method of order2 s . ✷ It turns out that such methods, in the form here presented (i.e., havingchosen the polynomial basis (9)), can be regarded as a generalization of Gaussmethods, in the sense that, they share the same nonzero spectrum, for all k ≥ s . In the limit k → ∞ , the same essential spectrum is retained by the limitoperator (see (14)). In the case of a polynomial Hamiltonian, such sequence ofmethods starts to be energy-preserving for a finite value of k . Moreover, eventhough for a general Hamiltonian the method becomes energy-preserving in thelimit k → ∞ , nevertheless, when using finite precision arithmetic, the limitis practically obtained for a finite value of k , namely as soon as full machineprecision accuracy is achieved.
5. HBVMs and Runge-Kutta collocation methods
By using the previous results and notations, we now further elucidate theexisting connections between HBVMs and Runge-Kutta collocation methods.Our starting point is a generic collocation method with k stages, defined by thetableau τ ... τ k A ω . . . ω k (28)where, for i, j = 1 , . . . , k : A = ( α ij ) ≡ (cid:18)Z τ i ℓ j ( x )d x (cid:19) , ω j = Z ℓ j ( x )d x,ℓ j ( τ ) being the j th Lagrange polynomial of degree k − { τ i } . Moreover, given a positive integer s ≤ k , and considering thematrices defined in (20)–(21), we consider the matrix P s P Ts Ω ∈ R k × k with projects into the s -dimensional subspace spanned by the columns of P s .The class of Runge-Kutta methods we are interested in, is then defined by thetableau τ ... τ k A ≡ A ( P s P Ts Ω) ω . . . . . . ω k (29)11e note that the Butcher array A has rank which cannot exceed s , becauseit is defined by filtering A by the rank s matrix P s P Ts Ω. The following resultthen holds true, which clarifies the existing connections between classical Runge-Kutta collocation methods and HBVMs.
Theorem 2.
Provided that the quadrature formula defined by the weights { ω i } is exact for polynomials at least s − (i.e., the Runge-Kutta method definedby the tableau (29) satisfies the usual simplifying assumption B (2 s ) ), then thetableau (29) defines a HBVM ( k, s ) method based at the abscissae { τ i } . Proof Let us expand the basis { P ( τ ) , . . . , P s ( τ ) } along the Lagrange basis { ℓ j ( τ ) } , j = 1 , . . . , k , defined over the nodes τ i , i = 1 , . . . , k : P j ( τ ) = k X r =1 P j ( τ r ) ℓ r ( τ ) , j = 1 , . . . , s. It follows that, for i = 1 , . . . , k and j = 1 , . . . , s : Z τ i P j ( x )d x = k X r =1 P j ( τ r ) Z τ i ℓ r ( x )d x = k X r =1 P j ( τ r ) α ir , that is (see (20)–(21) and (28)), I s = AP s . Consequently, AP s P Ts Ω = I s P Ts Ω , so that one retrieves the tableau (22) which defines the method HBVM( k, s ). ✷ The resulting Runge-Kutta method (29) is then energy conserving if appliedto polynomial Hamiltonian systems (1), when the degree of H ( y ) is lower thanor equal to a quantity, say ν , depending on k and s . As an example, whena Gaussian distribution of the nodes { τ i } is considered, one obtains (23) and,moreover, HBVM( k, s ) is also related to the Gauss-Legendre method of order2 k , according to (29), whose Butcher array coincides with A , with this choiceof the nodes { τ i } . Remark 2.
It seems like the price paid to achieve such conservation propertyconsists in the lowering of the order of the new method with respect to the originalone (28) . Actually this is not true, because a fair comparison would be to relatemethod (22) – (29) to a collocation method constructed on s rather than on k stages, since the actual nonlinear system, deriving by the implementation ofHBVM ( k, s ) , turns out to have dimension s , as it has been shown in [3]. Further implications of the isospectral property of HBVMs, among which analternative proof for their order of convergence, may be found in [6]. A furtheralternative proof can be found in [7, 8].12 . Conclusions
In this paper, we have shown that the recently introduced class of energy-preserving methods { HBVM( k, s ) } , when recast as Runge-Kutta methods, havethe matrix of the corresponding Butcher tableau sharing the same nonzero eigen-values which, in turn, coincides with those of the matrix of the Butcher tableauof the Gauss method of order 2 s , for all k ≥ s such that B (2 s ) holds.Moreover, HBVM( k, s ) defined at the Gaussian nodes { τ , . . . , τ k } on theinterval [0 ,
1] are closely related to the Gauss method of order 2 k which, althoughsymplectic, is not in general energy-preserving. References [1] M. Abramovitz, I.A. Stegun.
Handbook of Mathematical Functions . Dover,1965.[2] L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian BVMs (HBVMs): afamily of “drift-free” methods for integrating polynomial Hamiltonian sys-tems. “Proceedings of ICNAAM 2009”,
AIP Conf. Proc. (2009) 715–718.[3] L. Brugnano, F. Iavernaro, D. Trigiante. Analisys of Hamiltonian Bound-ary Value Methods (HBVMs): a class of energy-preserving Runge-Kuttamethods for the numerical solution of polynomial Hamiltonian dynamicalsystems.
Preprint (2009). ( arXiv:0909.5659 )[4] L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian Boundary ValueMethods (Energy Preserving Discrete Line Integral Methods).
Jour.of Numer. Anal., Industr. and Appl. Math. ,1-2 (2010) 17–37.( arXiv:0910.3621 )[5] L.Brugnano, F.Iavernaro, D.Trigiante. The Hamiltonian BVMs (HBVMs)Homepage . arXiv:1002.2757 (URL: http://web.math.unifi.it/users/brugnano/HBVM/ )[6] L.Brugnano, F.Iavernaro, D.Trigiante. Isospectral Property of HamiltonianBoundary Value Methods (HBVMs) and their connections with Runge-Kutta collocation methods. Preprint , 2010. ( arXiv:1002.4394 )[7] L.Brugnano, F.Iavernaro, D.Trigiante. Numerical Solution of ODEs and theColumbus’ Egg: Three Simple Ideas for Three Difficult Problems.
Math-ematics in Engineering, Science and Aerospace , No. 4 (2010) 105–124.( arXiv:1008.4789 )[8] L.Brugnano, F.Iavernaro, D.Trigiante. A unifying framework for the deriva-tion and analysis of effective classes of one-step methods for ODEs. Preprint , 2010. ( arXiv:1009.3165 )139] L. Brugnano, D. Trigiante.
Solving Differential Problems by Multistep Ini-tial and Boundary Value Methods . Gordon and Breach, Amsterdam, 1998.[10] E. Celledoni, R.I. McLachlan, D. McLaren, B. Owren, G.R.W. Quispel,W.M. Wright. Energy preserving Runge-Kutta methods.
M2AN (2009)645–649.[11] P. Chartier, E. Faou, A. Murua. An algebraic approach to invariant preserv-ing integrators: The case of quadratic and Hamiltonian invariants. Numer.Math. (2006) 575–590.[12] K. Feng. On difference schemes and symplectic geometry.
Proceedings of the5-th Intern. Symposium on differential geometry & differential equations,August 1984, Beijing (1985) 42–58.[13] H. Goldstein, C.P. Poole, J.L. Safko.
Classical Mechanics , Addison Wesley,2001.[14] E. Hairer, C. Lubich, G. Wanner.
Geometric Numerical Integration.Structure-Preserving Algorithms for Ordinary Differential Equations, 2 nd ed. , Springer, Berlin, 2006.[15] E. Hairer, G. Wanner. Solving Ordinary Differential Equations II , Springer,Berlin, 1991.[16] F. Iavernaro, B. Pace. s -Stage Trapezoidal Methods for the Conservation ofHamiltonian Functions of Polynomial Type. AIP Conf. Proc. (2007)603–606.[17] F. Iavernaro, B. Pace. Conservative Block-Boundary Value Methods for theSolution of Polynomial Hamiltonian Systems.
AIP Conf. Proc. (2008)888–891.[18] F. Iavernaro, D. Trigiante. State-dependent symplecticity and area preserv-ing numerical methods.
J. Comput. Appl. Math. no. 2 (2007) 814–825.[19] F. Iavernaro, D. Trigiante. High-order symmetric schemes for the energyconservation of polynomial Hamiltonian problems.
J. Numer. Anal. Ind.Appl. Math. ,1-2 (2009) 87–101.[20] A. Iserles, A. Zanna. Preserving algebraic invariants with Runge-Kuttamethods. J. Comput. Appl. Math. (2000) 69–81.[21] B. Leimkuhler, S. Reich.
Simulating Hamiltonian Dynamics , CambridgeUniversity Press, Cambridge, 2004.[22] R.I. McLachlan, G.R.W. Quispel, N. Robidoux. Geometric integration usingdiscrete gradient.
Phil. Trans. R. Soc. Lond. A , (1999) 1021-1045.[23] G.R.W. Quispel, D.I. McLaren. A new class of energy-preserving numericalintegration methods. J. Phys. A: Math. Theor. (2008) 045206 (7pp).1424] R.D. Ruth. A canonical integration technique. IEEE Trans. Nuclear Science ,4 (1983) 2669–2671.[25] J.M. Sanz-Serna and M.P. Calvo. Numerical Hamiltonian Problems , Chap-man & Hall, London, 1994.[26] F. Iavernaro, D. Trigiante. Discrete Mathematics, Discrete Physics and Nu-merical methods.
Le Matematiche