Chebyshev symplectic methods based on continuousstage RungeKutta methods
Abstract
We develop Chebyshev symplectic methods based on Chebyshev orthogonal polynomials of the first and second kind separately in this paper. Such type of symplectic methods can be conveniently constructed with the newlybuilt theory of weighted continuousstage RungeKutta methods. A few numerical experiments are well performed to verify the efficiency of our new methods.
keywords:
Hamiltonian systems; Symplectic methods; Continuousstage RungeKutta methods; Jacobi polynomials; Legendre polynomials; Chebyshev polynomials.1 Introduction
Geometric numerical integration of Hamiltonian systems has been a central topic in numerical solution of differential equations since the late 1980s Feng95kfc (); Fengq10sga (); hairerlw06gni (); Leimkuhlerr04shd (); ruth83aci (); sanzc94nhp (). The wellknown Hamiltonian systems can be written in a compact form, i.e.,
(1.1) 
where is a standard structure matrix, is the Hamiltonian function. Symplecticity (Poincaré 1899) has been discovered to be a characteristic property of Hamiltonian systems (see hairerlw06gni (), page 185), and thus it is suggested to construct numerical methods that share this geometric property. Such type of specialpurpose methods were naturally named to be “symplectic” Feng95kfc (); hairerlw06gni (); Leimkuhlerr04shd (); ruth83aci (); sanzc94nhp (), which states that the discrete numerical flow induced by the algorithms is a symplectic transformation, i.e., satisfying
where represents the Jacobian matrix of the numerical flow . Usually, especially for those (near)integrable systems, symplectic methods can produce many excellent numerical behaviors including linear error growth, longtime nearconservation of first integrals, and existence of invariant tori hairerlw06gni (); Shang99kam (). Moreover, by backward error analysis the numerical flow of the symplectic methods lies in the trajectories of the interpolating Hamiltonian systems, which implies that it exactly preserves a modified Hamiltonian Benetting94oth (); Tang94feo ().
There exists a particularly important class of symplectic methods called “symplectic RungeKutta (RK) methods”, which were discovered independently by three authors in 1988 sanz88rkm (); suris88otc (); lasagni88crk (). Afterwards, symplectic RK methods were fully explored in the context of classic RK methods (see, for example, sanz92srk (); sun93coh (); sun00asw ()), and the wellknown transformation technique proposed by Hairer & Wanner hairerw96sod () was frequently used. More recently, however, RK methods have been creatively extended to RK methods with continuous stage butcher72ato (); butcher87tna (); hairer10epv (), and thus symplectic RKtype methods gained a new growth point tangs12ana (); Tangs12tfe (); Tang13tfe (); Tangs14cor (); Tanglx16cos (); Tangsc17dgm (); Tang18ano (); Tangz18spc (); Tang18csr ().
In this paper, we further develop symplectic RKtype methods within the newlydeveloped framework of continuousstage RK methods. By using Chebyshev polynomials, it enables us to get rich production of Chebyshev symplectic methods. It should be recognized that transformation is closely related to Legendre polynomials while our approach can be applied to any other weighted orthogonal polynomials (although Chebyshev polynomials are mainly involved in the construction of our methods in this paper). On account of this, our approach to construct symplectic methods is rather different from the transformation technique previously used.
This paper will be organized as follows. In the next section, we give a brief revisit of some newlydeveloped theoretical results for constructing RKtype methods with general purpose. Section 3 is devoted to study the construction of Chebyshev RKtype methods with symplecticitypreserving property. Some numerical tests are given in Section 4. At last, we conclude this paper.
2 Theory of continuousstage RK methods
We are concerned with the following initial value problem
(2.1) 
with being sufficiently differentiable.
Definition 2.1.
hairer10epv () Let be a function of two variables , , and , be functions of . The onestep method given by
(2.2) 
is called a continuousstage RungeKutta (csRK) method, where Here, we always assume
(2.3) 
and often use a triple to represent such a method. In this paper, we will hold on the following assumption almost everywhere as previously done in hairer10epv (); Tangs14cor (); Tang18ano (); Tang18csr ()
(2.4) 
We introduce the following simplifying assumptions proposed by Hairer in hairer10epv ()
(2.5) 
The following result is useful for analyzing the order of csRK methods, which is a counterpart of the classic result by Butcher in 1964 butcher64ipa ().
Theorem 2.1.
hairer10epv (); Tangs14cor () If the coefficients of method (2.2) satisfy , and , then the method is of order at least .
Lemma 2.0.
Tang18csr () Under the assumption (2.4), the simplifying assumptions and are equivalent to
(2.6)  
(2.7)  
(2.8) 
where stands for the degree of polynomial function .
The concept of weight function is rather important for our discussions later, which can be found in almost every textbook of numerical analysis (see, for example, sulim03ait ()).
Definition 2.2.
A nonnegative function is called a weight function on , if it satisfies the following two conditions:

The th moment exists;

For , .
It is known that for a given weight function , there exists a sequence of orthogonal polynomials in the weighted function space (Hilbert space) Szeg85op ()
with respect to the inner product
In what follows, we denote the orthogonal polynomials by and assume they have been normalized in , i.e.,
It is well to be reminded that these polynomials make up a complete orthogonal set in the Hilbert space and the degree polynomial has exactly real simple zeros in the open interval .
Assume and have the following decompositions
where is a weight function defined on , and then the csRK method (2.2) can be written as
(2.9) 
Theorem 2.3.
Tang18csr () Suppose^{1}^{1}1We use the notation to stand for the onevariable function in terms of , and can be understood likewise. , then, under the assumption (2.4) we have

holds has the following form in terms of the normalized orthogonal polynomials in :
(2.10) where are any real parameters;

holds has the following form in terms of the normalized orthogonal polynomials in :
(2.11) where are any real functions;

holds has the following form in terms of the normalized orthogonal polynomials in :
(2.12) where are any real functions.
For simplicity and practical application, we have to truncate the series (2.10) and (2.11) suitably according to our needs. Consequently, only the polynomial case of and needs to be considered. Besides, generally it is impossible to exactly compute the integrals of a csRK scheme (except that is a polynomial vector field), thus we have to approximate them with an point weighted interpolatory quadrature formula
(2.13) 
where
Here, we remark that for the simplest case , we define .
Thus, by applying the quadrature rule (2.13) to the weighted csRK method (2.9), it leads up to a traditional stage RK method
(2.14) 
where . After that, we can use the following result to determine the order of the resulting RK methods.
Theorem 2.4.
Tang18csr () Assume the underlying quadrature formula (2.13) is of order , and is of degree with respect to and of degree with respect to , and is of degree . If all the simplifying assumptions , and in (2.5) are fulfilled, then the standard RK method (2.14) is at least of order
where , and .
Proof.
Please refer to Tang18csr () for the details of proof. ∎
Next, we introduce the following optimal quadrature technique named “GaussChristoffel type” for practical use, though other suboptimal quadrature rules can also be considered Abramowitzs65hom (); sulim03ait ().
Theorem 2.5.
If are chosen as the distinct zeros of the normalized orthogonal polynomial of degree in , then the interpolatory quadrature formula (2.13) is exact for polynomials of degree , i.e., of the optimal order . If , then it has the following error estimate
for some , where is the leading coefficient of .
3 Construction of Chebyshev symplectic methods
It is known that Chebyshev polynomials as a special class of Jacobi polynomials are frequently used in various fields especially in the study of spectral methods (see Gottliebo77nao (); Guo98sma () and references therein). Particularly, zeros of Chebyshev polynomials of the first kind are often used in polynomial interpolation because the resulting interpolation polynomial minimizes the effect of Runge’s phenomenon. But unfortunately, so far as we know, there are few Chebyshev symplectic methods available in the scientific literature except for two methods given in Tang18csr (). On account of this, we are interested in such subject and try to develop these methods based on the previous work of Tang18csr ().
The construction of symplectic methods is mainly dependent on the following results (please refer to Tang18ano (); Tang18csr () for more information).
Theorem 3.6.
Tang18ano () If the coefficients of a csRK method (2.2) satisfy
(3.1) 
then it is symplectic. In addition, the RK scheme with coefficients (derived by using quadrature formula, c.f., (2.14)) based on the underlying symplectic csRK method with coefficients satisfying (3.1) is always symplectic.
Theorem 3.7.
Tang18csr () Under the assumption (2.4), for a symplectic csRK method with coefficients satisfying (3.1), we have the following statements:

and , where ;

and , where .
Theorem 3.8.
Tang18csr () Suppose that , then symplectic condition (3.6) is equivalent to the fact that has the following form in terms of the orthogonal polynomials in
(3.2) 
where is skewsymmetric, i.e., .
By virtue of these theorems and the relevant results given in the
previous section, we can introduce the following procedure
for constructing symplectic csRK methods^{2}^{2}2Then, symplectic
RK methods can be obtained easily by using any quadrature rule, as
revealed by Theorem
3.6. Tang18csr ():
Step 1. Make an ansatz for which satisfies
with according to (2.10), and a
finite number of could be kept as parameters;
Step 2. Suppose is in the form
(according to Theorem 3.8)
where are kept as parameters with a finite number, and then substitute into (see (2.7), usually we let ):
for the sake of settling ;
Step 3. Write down and
(satisfy and automatically),
which results in a symplectic method of order at least
with
by Theorem 2.1 and
3.7.
However, the procedure above only provides a general framework for establishing symplectic methods. For simplicity and practical use, it needs to be more refined or particularized. Actually, in view of Theorem 2.4 and 3.7, it is suggested to design Butcher coefficients with lowdegree and , and is better to take as . Besides, for the sake of conveniently computing those integrals of in the second step, the following ansatz may be advisable (with given by (2.4) and let and )
(3.3) 
where . Because of the index restricted by in the second formula of (3.3), we can use to arrive at (please c.f. (2.6))
Therefore, implies that
(3.4) 
Finally, it needs to settle by transposing, comparing or merging similar items of (3.4) after the polynomial on righthand side being represented by the basis . In view of the skewsymmetry of , if we let , then actually the degrees of freedom of these parameters is , by noticing that
When (number of equations), i.e., , we can appropriately reduce the degrees of freedom of these parameters by imposing some of them to be zero in pairs, if needed.
3.1 Chebyshev symplectic methods of the first kind
Firstly, let us consider using the following shifted normalized Chebyshev polynomials of the first kind denoted by , i.e.,
It is known that these Chebyshev polynomials have the following properties:
(3.5) 
Notice that the properties given in (3.5) are helpful for computing the integrals^{3}^{3}3Of course, we can use some symbolic computing tool or softwares (e.g., Mathematica, Maple, Maxima etc.) to treat these integrals alternatively. of (3.4), hence we can conveniently construct Chebyshev symplectic methods of the first kind. Next, we give some examples and the following shifted GaussChristoffelChebyshev(I) quadrature rule will be used Abramowitzs65hom ()
(3.6) 
where
with being the zeros of Chebyshev polynomial .
Example 3.1.
With the orthogonal polynomials in (3.3) replaced by , we consider the following three cases separately,
 (i)

Let , we have only one degree of freedom. After some elementary calculations, it gives a unique solution
which results in a symplectic method of order . By using the point GaussChristoffelChebyshev(I) quadrature rule we regain the wellknown implicit midpoint rule;
 (ii)

Let , it will lead to
If we let be a free parameter, then we get a family of parameter symplectic csRK methods of order . Actually, it is easy to verify that the resulting methods are also symmetric^{4}^{4}4See Theorem 4.6 in Tang18csr (). and thus they possess an even order . By using the point GaussChristoffelChebyshev(I) quadrature rule we get a family of stage order symplectic RK methods which are shown in Tab. 3.1, with . We find that this class of methods is exactly the same one as shown in Tang18csr ().
 (iii)

If we take , then it gives a unique solution
The resulting symplectic csRK method is symmetric and of order . By using the point GaussChristoffelChebyshev(I) quadrature rule we get a stage order symplectic RK method which is shown numerically (the exact Butcher tableau is too complicated to be exhibited) in Tab. 3.2. It is tested that such method satisfies the classic symplectic condition (i.e., stability matrix sanz88rkm ()) and order conditions (from order to order ) up to the machine error.
3.2 Chebyshev symplectic methods of the second kind
Secondly, let us consider the shifted normalized Chebyshev polynomials of the second kind denoted by , i.e.,
The following properties can be easily verified (define )
(3.7) 
where the last formula is deduced from
Applying the properties given in (3.7) to the integrals of (3.4), it produces Chebyshev symplectic methods of the second kind. In our examples below, the following shifted GaussChristoffelChebyshev(II) quadrature rule will be used Abramowitzs65hom ()
(3.8) 
where
with being the zeros of as well as the inner extrema on of .
Example 3.2.
With the orthogonal polynomials in (3.3) replaced by , we consider the following three cases separately,
 (i)

Let , we have only one degree of freedom. After some elementary calculations, it gives a unique solution
which results in a symplectic csRK method of order . By using the point GaussChristoffelChebyshev(II) quadrature rule it gives the implicit midpoint rule;
 (ii)

Let , after some elementary calculations, it gives
If we regard as a free parameter, then we get a family of parameter symplectic and symmetric csRK methods of order . By using the point GaussChristoffelChebyshev(II) quadrature rule we get a family of stage order symplectic RK methods which are shown in Tab. 3.3, with .
 (iii)

Alternatively, if we take , then it gives a unique solution
The resulting symplectic csRK method is symmetric and of order . By using the point GaussChristoffelChebyshev(II) quadrature rule we get a stage order symplectic RK method which is shown in Tab. 3.4.
4 Numerical tests
We consider the perturbed Kepler’s problem given by the Hamiltonian function Calvofmr09sos ()
with the initial value condition . The exact solution is
In our numerical experiments, we take and use the step size . The Chebyshev symplectic methods of order 4 given in Tab. 3.1 (with , denoted by “Chebyshev I order 4”) and Tab. 3.3 (with , denoted by “Chebyshev II order 4”) are tested comparing with the wellknown GaussLegendre RK method of order 4 (denoted by “Gauss order 4”). It is observed from Fig. 4.2 and Fig. 4.2 that our Chebyshev symplectic methods of order 4 share very similar numerical behaviors with the classic method “Gauss order 4”, although the latter has a little bit better result in the aspects of growth of solution error and conservation of energy. As is expected, we have a bounded error in energy conservation and a linear growth of solution error, which coincides well with the common view in general symplectic integration Fengq10sga (); hairerlw06gni (). Besides, the newlyderived Chebyshev symplectic methods of order 6, denoted by “Chebyshev I order 6” and “Chebyshev II order 6” respectively (see Tab. 3.2 and 3.4), are also tested comparing with the 6order GaussLegendre RK method, the numerical results of which are shown in Fig. 4.4 and Fig. 4.4. It is seen that these symplectic methods almost exhibit the same numerical results. These numerical tests are well conformed with our expects and the associated theoretical results. Therefore, the newlyconstructed Chebyshev methods are effective for solving Hamiltonian systems.
5 Conclusions
This paper intensively discusses the construction of Chebyshev symplectic RKtype methods with the help of the newlybuilt theory for csRK methods. We present a new family of symplectic RK methods in use of the Chebyshev polynomials of the first and second kind separately. Although these methods are developed mainly in terms of Chebyshev polynomials, they essentially can be directly extended to other types of orthogonal polynomials. In addition, we notice that ChebyshevGaussLobatto collocation methods have been considered in zhangk11coa () for solving Hamiltonian systems, stating that these spectral collocation methods can preserve both energy and symplectic structure up to the machine error in each time step. But their methods are nonsymplectic after all, it can not guarantee the correct qualitative behaviors for a rather long term. In contrast to this, by using the interpolatory quadrature rules with Chebyshev abscissae, we have constructed the Chebyshev methods which are exactly symplectic.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (11401055), China Scholarship Council and Scientific Research Fund of Hunan Provincial Education Department (15C0028).
References
 (1) M. Abramowitz, I.A. Stegun, Eds., Handbook of Mathematical Functions, Dover, New York, 1965.
 (2) G. Benettin, A. Giorgilli, On the Hamiltonian interpolation of NeartotheIdentity symplectic mappings with application to symplectic integration algorithms, J. Statist. Phys., 74 (1994), 11171143.
 (3) J.C. Butcher, Implicit RungeKutta processes, Math. Comput. 18 (1964), 50–64.
 (4) J.C. Butcher, An algebraic theory of integration methods, Math. Comp., 26 (1972), 79106.
 (5) J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations: RungeKutta and General Linear Methods, John Wiley & Sons, 1987.
 (6) M. Calvo, J.M. Franco, J.I. Montijano, L. Rández, Sixthorder symmetric and symplectic exponentially fitted RungeKutta methods of the Gauss type, J. Comput. Appl. Math., 223 (2009), 387–398.
 (7) K. Feng, K. Feng’s Collection of Works, Vol. 2, Beijing: National Defence Industry Press, 1995.
 (8) K. Feng, Q. Mengzhao, Symplectic Geometric Algorithms for Hamiltonian Systems, Zhejiang Publishing United Group, Zhejiang Science and Technology Publishing House, Hangzhou and SpringerVerlag Berlin Heidelberg, 2010.
 (9) D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAMCBMS, Philadelphia, 1977.
 (10) B. Guo, Spectral Methods and their Applications, World Scietific, Singapore, 1998.
 (11) E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordiary Differential Equations I: Nonstiff Problems, Springer Series in Computational Mathematics, 8, SpringerVerlag, Berlin, 1993.
 (12) E. Hairer, G. Wanner, Solving Ordiary Differential Equations II: Stiff and DifferentialAlgebraic Problems, Second Edition, Springer Series in Computational Mathematics, 14, SpringerVerlag, Berlin, 1996.
 (13) E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: StructurePreserving Algorithms For Ordinary Differential Equations, Second edition, Springer Series in Computational Mathematics, 31, SpringerVerlag, Berlin, 2006.
 (14) E. Hairer, Energypreserving variant of collocation methods, JNAIAM J. Numer. Anal. Indust. Appl. Math. 5 (2010), 73–84.
 (15) F. Lasagni, Canonical RungeKutta methods, ZAMP 39 (1988), 952–953.
 (16) B. Leimkuhler, S. Reich, Simulating Hamiltonian dynamics, Cambridge University Press, Cambridge, 2004.
 (17) R. Ruth, A canonical integration technique, IEEE Trans. Nucl. Sci., 30 (1983), 2669–2671.
 (18) J.M. SanzSerna, M.P. Calvo, Numerical Hamiltonian problems, Chapman & Hall, 1994.
 (19) J. M. SanzSerna, RungeKutta methods for Hamiltonian systems, BIT 28 (1988), 877–883.
 (20) J. M. SanzSerna, Symplectic RungeKutta and related methods: recent results, Physica. D., 60 (1992), 293–302.
 (21) Z. Shang, KAM theorem of symplectic algorithms for Hamiltonian systems, Numer. Math., 83 (1999), 477–496.
 (22) E. Süli, D. F. Mayers, An Introduction to Numerical Analysis, Cambridge University Press, 2003.
 (23) G. Sun, Construction of high order symplectic RungeKutta methods, J. Comput. Math., 11 (1993), 250–260.
 (24) G. Sun, A simple way constructing symplectic RungeKutta methods, J. Comput. Math., 18 (2000), 61–68.
 (25) Y. B. Suris, On the conservation of the symplectic structure in the numerical solution of Hamiltonian systems (in Russian), In: Numerical Solution of Ordinary Differential Equations, ed. S.S. Filippov, Keldysh Institute of Applied Mathematics, USSR Academy of Sciences, Moscow, 1988, 148–160.
 (26) G. Szegö, Orthogonal Polynomials, vol. 23, Amer. Math. Soc., 1985.
 (27) W. Tang, Y. Sun, A new approach to construct RungeKutta type methods and geometric numerical integrators, AIP. Conf. Proc. 1479 (2012), 1291–1294.
 (28) W. Tang, Y. Sun, Time finite element methods: A unified framework for numerical discretizations of ODEs, Appl. Math. Comput. 219 (2012), 2158–2179.
 (29) W. Tang, Time finite element methods, continuousstage RungeKutta methods and structurepreserving algorithms, PhD thesis, Chinese Academy of Sciences, 2013.
 (30) W. Tang, Y. Sun, Construction of RungeKutta type methods for solving ordinary differential equations, Appl. Math. Comput., 234 (2014), 179–191.
 (31) W. Tang, Y. Sun, J. Zhang, High order symplectic integrators based on continuousstage RungeKuttaNyström methods, preprint, 2015.
 (32) W. Tang, G. Lang, X. Luo, Construction of symplectic (partitioned) RungeKutta methods with continuous stage, Appl. Math. Comput., 286 (2016), 279–287.
 (33) W. Tang, Y. Sun, W. Cai, Discontinuous Galerkin methods for Hamiltonian ODEs and PDEs, J. Comput. Phys., 330 (2017), 340–364.
 (34) W. Tang, J. Zhang, Symplecticitypreserving continuousstage RungeKuttaNyström methods, Appl. Math. Comput., 323 (2018), 204–219.
 (35) W. Tang, A note on continuousstage RungeKutta methods, submitted, 2018.
 (36) W. Tang, Continuousstage RungeKutta methods based on weighted orthogonal polynomials, preprint, 2018.
 (37) Y. Tang, Formal energy of a symplectic scheme for Hamiltonian systems and its applications (I), Computers Math. Applic., 27 (1994), 31–39.
 (38) N. Kanyamee, Z. Zhang, Comparison of a spectral collocation method and symplectic methods for Hamiltonian systems, INT. J. Numer. Anal. Mod., 8 (2011), 86–104.