Chebyshev symplectic methods based on continuous-stage Runge-Kutta methods

Chebyshev symplectic methods based on continuous-stage Runge-Kutta methods

Wensheng Tang tangws@lsec.cc.ac.cn College of Mathematics and Statistics,
Changsha University of Science and Technology,
Changsha 410114, China
Hunan Provincial Key Laboratory of Mathematical
Modeling and Analysis in Engineering,
Changsha 410114, China
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 newly-built theory of weighted continuous-stage Runge-Kutta methods. A few numerical experiments are well performed to verify the efficiency of our new methods.

keywords:
Hamiltonian systems; Symplectic methods; Continuous-stage Runge-Kutta 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 well-known 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 special-purpose 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, long-time near-conservation 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 Runge-Kutta (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 well-known -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 RK-type methods gained a new growth point tangs12ana (); Tangs12tfe (); Tang13tfe (); Tangs14cor (); Tanglx16cos (); Tangsc17dgm (); Tang18ano (); Tangz18spc (); Tang18csr ().

In this paper, we further develop symplectic RK-type methods within the newly-developed framework of continuous-stage 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 newly-developed theoretical results for constructing RK-type methods with general purpose. Section 3 is devoted to study the construction of Chebyshev RK-type methods with symplecticity-preserving property. Some numerical tests are given in Section 4. At last, we conclude this paper.

2 Theory of continuous-stage 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 one-step method given by

(2.2)

is called a continuous-stage Runge-Kutta (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 non-negative 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 () Suppose111We use the notation to stand for the one-variable 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 “Gauss-Christoffel 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 skew-symmetric, 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 methods222Then, 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 low-degree 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 right-hand side being represented by the basis . In view of the skew-symmetry 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 integrals333Of 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 Gauss-Christoffel-Chebyshev(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 Gauss-Christoffel-Chebyshev(I) quadrature rule we regain the well-known 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 symmetric444See Theorem 4.6 in Tang18csr (). and thus they possess an even order . By using the -point Gauss-Christoffel-Chebyshev(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 Gauss-Christoffel-Chebyshev(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.

Table 3.1: A family of one-parameter -stage -order symplectic RK methods, based on Chebyshev polynomials of the first kind.
Table 3.2: A -stage -order symplectic RK method, based on Chebyshev polynomials of the first kind.

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 Gauss-Christoffel-Chebyshev(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 Gauss-Christoffel-Chebyshev(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 Gauss-Christoffel-Chebyshev(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 Gauss-Christoffel-Chebyshev(II) quadrature rule we get a -stage -order symplectic RK method which is shown in Tab. 3.4.

Table 3.3: A family of one-parameter -stage -order symplectic RK methods, based on Chebyshev polynomials of the second kind.
Table 3.4: A -stage -order symplectic RK method, based on Chebyshev polynomials of the second kind.

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 well-known Gauss-Legendre 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 newly-derived 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 6-order Gauss-Legendre 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 newly-constructed Chebyshev methods are effective for solving Hamiltonian systems.

Figure 4.1: Solution error by three symplectic methods of order 4, step size .
Figure 4.2: Energy error by three symplectic methods of order 4, step size .
Figure 4.3: Solution error by three symplectic methods of order 6, step size .
Figure 4.4: Energy error by three symplectic methods of order 6, step size .

5 Conclusions

This paper intensively discusses the construction of Chebyshev symplectic RK-type methods with the help of the newly-built 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 Chebyshev-Gauss-Lobatto 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 non-symplectic 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 Near-to-the-Identity symplectic mappings with application to symplectic integration algorithms, J. Statist. Phys., 74 (1994), 1117-1143.
  • (3) J.C. Butcher, Implicit Runge-Kutta processes, Math. Comput. 18 (1964), 50–64.
  • (4) J.C. Butcher, An algebraic theory of integration methods, Math. Comp., 26 (1972), 79-106.
  • (5) J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods, John Wiley & Sons, 1987.
  • (6) M. Calvo, J.M. Franco, J.I. Montijano, L. Rández, Sixth-order symmetric and symplectic exponentially fitted Runge-Kutta 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 Springer-Verlag Berlin Heidelberg, 2010.
  • (9) D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM-CBMS, 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, Springer-Verlag, Berlin, 1993.
  • (12) E. Hairer, G. Wanner, Solving Ordiary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Edition, Springer Series in Computational Mathematics, 14, Springer-Verlag, Berlin, 1996.
  • (13) E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms For Ordinary Differential Equations, Second edition, Springer Series in Computational Mathematics, 31, Springer-Verlag, Berlin, 2006.
  • (14) E. Hairer, Energy-preserving variant of collocation methods, JNAIAM J. Numer. Anal. Indust. Appl. Math. 5 (2010), 73–84.
  • (15) F. Lasagni, Canonical Runge-Kutta 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. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian problems, Chapman & Hall, 1994.
  • (19) J. M. Sanz-Serna, Runge-Kutta methods for Hamiltonian systems, BIT 28 (1988), 877–883.
  • (20) J. M. Sanz-Serna, Symplectic Runge-Kutta 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 Runge-Kutta methods, J. Comput. Math., 11 (1993), 250–260.
  • (24) G. Sun, A simple way constructing symplectic Runge-Kutta 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 Runge-Kutta 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, continuous-stage Runge-Kutta methods and structure-preserving algorithms, PhD thesis, Chinese Academy of Sciences, 2013.
  • (30) W. Tang, Y. Sun, Construction of Runge-Kutta 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 continuous-stage Runge-Kutta-Nyström methods, preprint, 2015.
  • (32) W. Tang, G. Lang, X. Luo, Construction of symplectic (partitioned) Runge-Kutta 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, Symplecticity-preserving continuous-stage Runge-Kutta-Nyström methods, Appl. Math. Comput., 323 (2018), 204–219.
  • (35) W. Tang, A note on continuous-stage Runge-Kutta methods, submitted, 2018.
  • (36) W. Tang, Continuous-stage Runge-Kutta 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.
Comments 0
Request Comment
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
   
Add comment
Cancel
Loading ...
217668
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description