CC BY
8
Discrete & Continuous Models 2019 27 (1) 33—41
. „ . http://journals.rudn.ru/miph
& Applied Computational Science
Research article
UDC 512, 519.6
DOI: 10.22363/2658-4670-2019-27-1-33-41
The symbolic problems associated with Runge—Kutta methods and their solving in Sage
Yu Ying12
1 Department of Applied Probability and Informatics Peoples' Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya str., Moscow, 117198, Russian Federation 2 Department of Algebra and Geometry
Kaili University 3, Kaiyuan Road, Kaili, 556011, China
(received: February 21, 2019; accepted: October 21, 2019)
Runge-Kutta schemes play a very important role in solving ordinary differential equations numerically. At first we want to present the Sage routine for calculation of Butcher matrix, we call it an rk package. We tested our Sage routine in several numerical experiments with standard and symplectic schemes and verified our result by corporation with results of the calculations made by hand.
Second, in Sage there are the excellent tools for investigation of algebraic sets, based on Grobner basis technique. As we all known, the choice of parameters in Runge-Kutta scheme is free. By the help of these tools we study the algebraic properties of the manifolds in affine space, coordinates of whose are Butcher coefficients in Runge-Kutta scheme. Results are given both for explicit Runge-Kutta scheme and implicit Runge-Kutta scheme by using our rk package. Examples are carried out to justify our results. All calculation are executed in the computer algebra system Sage.
Key words and phrases: Sympletic Runge-Kutta Scheme; Grobner basis, Sage, Sagemath
1. Introduction
The Runge-Kutta method is the most popular numerical method for solving of ordinary differential equations (ODE), however the development of this method indicates some symbolic problems.
Let's review some results on Runge-Kutta scheme [1]. For an autonomous system
X = F (x) (1)
Runge-Kutta scheme with s stages can be written as
ki = F ^x + dt JJ aijkjj , i = 1,2,..., s (2)
© YingYu., 2019
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/
and
s
x = x + dt^^biki. (3)
i= 1
We will write the coefficients «j and bi in Butcher matrix, for example, for s = 2
Cl «11 «12
C2 «21 «22
bi b2
Hereinafter ci = aij.
j=i
These coefficients are selected in such a way that the difference scheme approximates the ODE (1) with some order p. For s = p = 4 appropriate numerical values for the coefficients was found by Kutta in the XlX century. The usage of Runge-Kutta scheme for the solving of given ODEs means the numerical calculations in the floating-point arithmetic.
However the problem about finding the coefficients of Butcher matrix is pure algebraic, so now we can try to solve it with the help of computer algebra systems like Sage or Maple [2-4]. This is the first symbolic problem associated with Runge-Kutta method. In our paper we want to present our algorithms for symbolical calculation of the Butcher matrix and its realization in Sage.
It should be noted that conditions of approximation and other conditions don't define the coefficients of Butcher matrix unambiguously. From geometrical point of view the list of coefficients of Butcher matrix is a point in affine space (Butcher space) and appropriate points form a variety in this space. As in Sage there are some tools for a research of a set of solutions of systems of the algebraic equations, we can try to describe the varieties in Butcher spaces.
The Runge-Kutta scheme is explicit iff aij = 0, j ^ i, in this case the numerical calculation doesn't demand the solving of nonlinear algebraic equations. This is the most investigated case, for small s the Butcher varieties ware described by Butcher himself, now we know approximate values of aij for schemes with s = 9 [5-7].
The implicit Runge-Kutta scheme is interesting because they can save some symbolic properties of exact solution. By Cooper theorem symplectic Runge-Kutta scheme saves exactly all quadratic integrals of motion [8,9]. The symplecticity gives algebraic conditions for Butcher matrix, and we can try to investigate the varieties in Butcher space by the help of Sage.
In our numerical experiments these varieties were rational with singularity at point with aij = 0 at all values of i and j. Thus there are infinite set of Butcher matrices with rational coefficients.
Organization of calculations by the implicit Runge-Kutta scheme demands solving of algebraic equations at each step. Investigation of the algebraic system is the second symbolic problem associated with Runge-Kutta scheme which we can try to solve by Sage.
2. The problem of finding coefficients of Butcher matrix, approximative conditions
In this part, we try to set up a program for finding Runge-Kutta scheme when s is number of the stages and p is order of approximation. This scheme we can call rk(s,p). Butcher matrix is square matrix of degree s. Condition of approximation is algebraic equation with respect of elements of Butcher matrix. Firstly, we note that these equations don't depend on the number of the approximative differential equations, so we can only consider the first order equation X = f (x).
The Taylor series for the solution satisfying the zero initial condition
x(to) = 0 are
x(t) = D(f )lx=o (t - to) + D2(f)U + ...,
where
D = ff
dx
At t = ti we have
dt2
x(ti) = D(f )Uo dt + D (f )|x=o - g!" + ..., (4)
where dt = t1 — t0.
The first step of Runge-Kutta scheme with s stages can be written as
ki = f (^dt ^^ aij k^j , i = 1, 2,..., s (5)
and
s
X = dt^^biki. (6)
i=1
Let the power series for the solution of systems (5) be
hi = ki,o + ki,idt + ki,2dt2 +----
Manually or in computer algebra system we can find their coefficients by recursion. So at dt = 0 the equations (5) give us
ki,o = f (0).
To find the coefficients at dt, we differentiate the equations (5) with respect to dt and substitute dt = 0, then we have
ki,i = f'(0) ■ ^^ aijkj,o^ , i = 1, 2,..., s.
To find the coefficients of dt2, we differentiate the equations (5) twice and substitute dt = 0 again, etc. Substituting the power series for k in (6) we have the power series for X. By definition, the difference scheme has the p-th order approximation, if first p terms of the power series (4) for exact solution at t = ti and the power series for X coincide. So, these are the required conditions of approximation written as a system of algebraic equations with respect to Butcher coefficients.
In the area of symbolic calculation, we consider the Butcher coefficients a, b as symbolical variables. We also consider the function f as arbitrary function of x, thus we may manipulate with its derivatives
f> = f(0), fi = m f2 = f"(o),...
as independent symbolic variables. Comparing the coefficients at dt,... dtp in the way stated above, we get p algebraic equations, left sides of which belong to the polynomial ring
Q[a,b][fo ,...,fp].
By the arbitrariness of f, coefficients of all the monomials ... must be zero. This gives us a number of equations for the coefficients of Butcher matrix, which are the required equations. Thus the approximation conditions are algebraic equations, left sides of which belong to the ring Q[a,b] and form an ideal in Q[a, b]. For brevity we will say that the approximation conditions generate the ideal.
In our package rk.sage function rk_var(s) gives all the necessary symbolical variables for the coefficients a, b and ring
Kab = Q[aii,...,bi,... ]
for the scheme with s stages. As we know the condition of approximation with order p is the system of equations, the function rk_series(s,p) returns the left side of these equations. Let's review an example.
sage: var( 'x,t,dt') (x, t, dt)
sage: load("rk.sage") None
sage: rk_var(3) None sage: a [a00 a01 a02] [a10 all a12] [a20 a21 a22] sage: b (b0, bl, b2) sage: c
[a00 + a01 + a02, a10 + all + a12, a20 + a21 + a22]
So we have all Butcher coefficients. The construction sage: rk_series(3,3)
return the result which coincides with the statement in the book of E. Hairer et al. [1, §II 2, theorem 2.1].
Any partial solution a, b of this algebraic equation system gives scheme type rk(3,3). Also in our package it is possible to work only with explicit schemes.
sage: rk_var(3,implicit=false) None sage: a [ 0 0 0] [a10 0 0] [a20 a21 0] sage: b (b0, bl, b2) sage: c
[0, a10, a20 + a21] sage: rk_series(3,3,implicit=false) [-b0 - b1 - b2 + 1, -2*a10*b1 - 2*a20*b2 - 2*a21*b2 + 1, -6*a10*a21*b2 + 1,
-3*a10~2*b1 - 3*a20~2*b2 - 6*a20*a21*b2 - 3*a21~2*b2 + 1]
This means that the affine variety of 3-stage Runge-Kutta scheme, approximating the differential equation with order 3, has 6 Butcher coefficients satisfying the 4 equations. These equations give some algebraic sets in the affine 6d-space. To study this set under Sage, we consider the corresponding ideal:
sage: J=Kab*rk_series(3,3,implicit=false) sage: J.dimension() 2
sage: J.is_prime() True
Thus, the set of coefficients formulate an irreducible surface. Its projection in the 3-dimension space ai0, a2o, a2i is a cubic surface:
sage: J.elimination_ideaK[Kab(b0),Kab(b1),Kab(b2)]) Ideal
(3*a10~2*a21 - a10*a20 - 3*a10*a21 + a20~2 + 2*a20*a21 + a21~2) of Multivariate Polynomial Ring in a10, a20, a21, b0, b1, b2 ^ over Rational Field
Since the linear family of right lines passing through the point O intersects this surface at one moving point, so this surface is rational.
We try investigate rk(4,4) in Sage. Our program give system of equations very quickly but Sage can't calculate Grobner basis of the ideal, generated by this system. We try use CAS Maple (which have several commercial routines for calculation Grobner basis) but without success. It should be noted that this system was investigated by hand in the XIX century!
By using Sage we calculate the dimensions of the algebraic sets for explicit and implicit RK-schemes. Results are presented in the Tables 1 and 2. The analysis of RK-schemes with bigger p cannot be made in such a way that the Grobner basis can be calculated.
Table 1
Dimensions for implicit Runge-Kutta scheme in the case of 1 < s < 4; 1 < p < 4
s,p p =2 p =3
s = 2 4 2
s = 3 10 8
Table 2
Dimensions for explicit Runge-Kutta scheme on the case of 1 < s < 4; 1 < p < 4. The value — 1 means that the set of Runge-Kutta scheme with given s and p is empty
s,p p =2 p =3
s = 2 4 — 1
s = 3 4 2
3. Symplectic Runge-Kutta schemes
We describe the manifold in the affine space ab, coordinates of whose points can be used as coefficients in the Runge-Kutta scheme, satisfying both the symplecticity and p-order approximation [8-10]. We denote this manifold as S(s,p). We can find their equations in Sage like:
sage: rk_var(2) None
sage: rk_sympletic(3,2) [2*a00*b0 - b0~2, a01*b0 + a10*b1 - b0*b1, a01*b0 + a10*b1 - b0*b1, 2*a11*b1 - b1~2, -b0 - b1 + 1,
-2*a00*b0 - 2*a01*b0 - 2*a10*b1 - 2*a11*b1 + 1, -6*a00~2*b0 - 6*a00*a01*b0 - 6*a00*a10*b1 - 6*a01*a10*b0
- 6*a01*a10*b1 - 6*a01*a11*b0 - 6*a10*a11*b1 - 6*a11~2*b1 ^ + 1,
-3*a00~2*b0 - 6*a00*a01*b0 - 3*a01~2*b0 - 3*a10~2*b1 - 6*a10*a11*b1 - 3*a11~2*b1 + 1]
The tool of Sage allows you to describe this manifold. For example, for S(2, 3) we have
sage: eqs = rk_sympletic(3,2) sage: J=Kab*eqs sage: J.dimension() 1
sage: J.is_prime() True
So S(2, 3) is irreducible curve. We can find its projection on any of the coordinate planes. In the aooaio-plane its equation is
sage: J.elimination_ideal([Kab(b1),Kab(b0),Kab(a01),Kab(a11)]) Ideal
(12*a00~3 - 24*a00~2*a10 - 6*a00~2 + 12*a00*a10~2 + 12*a00*a10 ^ + a00 - 6*a10~2)
of Multivariate Polynomial Ring in a00, a01, a10, a11, b0, b1 ^ over Rational Field
We can calculate the genus of the curve of the 3rd order:
sage: elvars=[Kab(b1),Kab(b0),Kab(a00),Kab(a01)] sage: F=J.elimination_ideal(elvars).gens()[0] sage: KK=QQ[a10,a11] sage: Curve(KK(F)).genus() 0
If genus of the curve is equal to zero the curve is rational and in Sage we can find the rational parametrization for given curve with the help of the construction
sage: Curve(KK(F)).rational_parameterization()
By using the elimination ideal technique we makes it easy to find that b0, bi, a00, a01 are linearly of a10 and an, thus coefficients of Butcher matrix for S(2, 3) can be expressed as rational functions of one parameter t:
ci
Ci
3t2 - 12t +12 8t2 - 24t + 24
-t3 + 8t2 - 24t + 24 8t3 - 24t2 + 24t
3t3 - 8t2 + 24t - 24 t3 - 8t2 - 24t - 24 + -
8t3 - 24t2 + 24t
6t(t - 2)2
t2
8t2 - 24t + 24
1
1 - -
t2
4 t2-3t+3
4t2 — 12t + 12
For example at the value t = 3 ± v^ we have the scheme used by Nuan et al. [11, f. (24)]. At rational curve there are the infinite set of rational points thus there are infinite set of Runge-Kutta schemes of type S(2, 3) with rational coefficients. For example at t =1 we have
2
t
1/3 3/8 -1/24
1 7/8 1/8
3/4 1/4
Increase in number of stages is interfered by difficulties, associated with the algorithm of calculation of Grobner basis. We try investigating the case s = 4 in Maple, where there are several non-free algorithms of calculation of Grobner basis, but without success.
In our numerical experiments the variate S(s,p) with maximal value of p at fixed value of s has the dimension 1 and the genus 0 and thus is a rational curve in Butcher space. In practice it is important that there is an infinite number of sets of rational values for Butcher coefficients.
4. Conclusion
In this article we investigated one of the symbolic problems associated with Runge-Kutta method, namely the problem of calculation of Butcher matrix. We tested our Sage routine in several numerical experiments and verified our results by comparing them with the results of calculations made by hand.
We saw that computer algebra systems was ready to solve this problem and gave us some tools for investigation of the set of Butcher matrices. These tools indicate in our numerical experiments that Butcher coefficients were rational functions of one or two parameters and thus there is an infinite number of sets of rational values for Butcher coefficients. Perhaps, this property can be proved for Runge-Kutta schemes with any number of stages. This question requires further study. Also we can't calculate Grobner basis of the ideals generated at s > 3. Perhaps, a successful substitution can solve this issue.
Acknowledgments
The publication has been prepared with the support of the "RUDN University Program 5-100" and funded by RFBR according to the research projects No. 18-07-00567 and No. 18-51-18005.
References
1. E. Hairer, G. Wanner, S. P. N0rsett, Solving ordinary differential equations, 3rd Edition, Vol. 1, Springer, New York, 2008.
2. The Sage Developers, SageMath, the Sage Mathematics Software System (Version 7.4) (2016).
URL https://www.sagemath.org
3. W. Stein, D. Joyner, SAGE: system for algebra and geometry experimentation, ACM SIGSAM Bulletin 39 (2) (2005) 61-64.
4. A. Y. Golubkov, A. I. Zobnin, S. O. V., Computer algebra in the Sage system. Manual [Komp'yuternaya algebra v sisteme Sage. Uchebnoye posobiye], Bauman Moscow State Technical University, Moscow, 2013, in Russian.
5. S. I. Khashin, A symbolic-numeric approach to the solution of the Butcher equations, Canadian Applied Mathematics Quarterly 17 (3) (2009) 555-569.
6. S. I. Khashin, Butcher algebras for Butcher systems, Numerical Algorithms 63 (4) (2013) 679-689. doi:10.1007/s11075-012-9647-x.
7. J. H. Verner, Explicit Runge-Kutta pairs with lower stage-order, Numerical Algorithms 65 (3) (2014) 555-577. doi:10.1007/s11075-012-9647-x.
8. E. Hairer, G. Wanner, C. Lubich, Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations, 2nd Edition, Springer, New York, 2000.
9. J. M. Sanz-Serna, Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more, SIAM REVIEW 58 (1) (2016) 3-33. doi:10.1137/151002769.
10. M. N. Gevorkyan, J. V. Gladysheva, Symplectic integrators and the problem of wave propagation in layered media, Bulletin of Peoples' Friendship University of Russia. Series: Mathematics. Information Sciences. Physics (1) (2012) 50-60, in Russian.
11. Nuan Fang Xu, Zi-Chen Deng, Yan Wang, Kai Zhang, A symplectic Runge-Kutta method for the analysis of the tethered satellite system, Multidiscipline Modeling in Materials and Structures 13 (1) (2017) 26-35. doi:10.1108/MMMS-11-2016-0060.
For citation:
Yu. Ying, The symbolic problems associated with Runge-Kutta methods and their solving in Sage, Discrete and Continuous Models and Applied Computational Science 27 (1) (2019) 33-41. DOI: 10.22363/2658-4670-201927-1-33-41.
Information about the authors:
Yu Ying (China) — postgraduate student of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University); assistant professor of Department of Algebra and Geometry, Kaili University (e-mail: yingy6165@gmail.com, phone: +7(495)9550927, ORCID: https://orcid.org/0000-0002-4105-2566, ResearcherID: AAC-8344-2019, Scopus Author ID: 57208127921)
