The symbolic problems associated with Runge-Kutta methods and their solving in Sage

Full Text

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 ( s \ �ki = F� �x + dt aij�kj j=1 , i = 1, 2, . . . , s (2) © Ying Yu., 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 bi�ki. (3) i=1 We will write the coefficients aij and bi in Butcher matrix, for example, for s = 2 c1 a11 a12 c2 a21 a22 b1 b2 s Hereinafter ci = aij . j=1 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 XIX 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 geomet- rical 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. 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(t0) = 0 are where 2 x(t) = D(f )|x=0 (t - t0) + D (f ) x=0 (t - t0)2 2! + . . . , At t = t1 we have d D = f . dx 2 dt2 (f ) x(t1) = D(f )|x=0 dt + D where dt = t1 - t0. x=0 + . . . , (4) 2! The first step of Runge-Kutta scheme with s stages can be written as and ki = f ( s \ dt aijkj j=1 s , i = 1, 2, . . . , s (5) xˆ = dt biki. (6) i=1 Let the power series for the solution of systems (5) be ki = ki,0 + ki,1dt + 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,0 = 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,1 = f t(0) · ( s \ aijkj,0 j=1 , 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 ki 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 = t1 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 f0 = f (0), f1 = f t(0), f2 = f tt(0), . . . 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][f0, . . . , fp]. 0 By the arbitrariness of f , coefficients of all the monomials fm0 . . . 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 sym- bolical variables for the coefficients a, b and ring Kab = Q[a11, . . . , b1, . . . ] 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 a11 a12] [a20 a21 a22] sage: b (b0, b1, b2) sage: c [a00 + a01 + a02, a10 + a11 + 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, b1, 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, ap- proximating 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 a10, a20, a21 is a cubic surface: sage: J.elimination_ideal([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 c→ 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 Gröbner basis of the ideal, generated by this system. We try use CAS Maple (which have several commercial routines for calculation Gröbner 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 Gröbner 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 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 c→ + 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 a00a10-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 c→ + a00 - 6*a10^2) of Multivariate Polynomial Ring in a00, a01, a10, a11, b0, b1 c→ 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, b1, a00, a01 are linearly of a10 and a11, thus coefficients of Butcher matrix for S(2, 3) can be expressed as rational functions of one parameter t: c1 c1 - 3t2 - 12t + 12 8t2 - 24t + 24 t3 + 8t2 - 24t + 24 + 3t3 - 8t2 + 24t - 24 8t3 - 24t2 + 24t t2 t3 - 8t2 - 24t - 24 6t(t - 2)2 8t3 - 24t2 + 24t 8t2 - 24t + 24 1 t2 t2 1 - 4 t2-3t+3 4t2 - 12t + 12 For example at the value t = 3 ± √3 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 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 Gröbner basis. We try investigating the case s = 4 in Maple, where there are several non-free algorithms of calculation of Gröbner 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. 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 Gröbner basis of the ideals generated at s > 3. Perhaps, a successful substitution can solve this issue.

About the authors

Yu Ying

Peoples’ Friendship University of Russia (RUDN University); Kaili University

Author for correspondence.
Email: yingy6165@gmail.com

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

6, Miklukho-Maklaya str., Moscow, 117198, Russian Federation; 3, Kaiyuan Road, Kaili, 556011, China

References