An Embedded Fourth Order Method for Solving Structurally Partitioned Systems of Ordinary Diﬀerential Equations

An embedded one-step numerical method for structurally partitioned systems of ordinary diﬀerential equations (ODEs) is considered. Two-parametric families of methods of order four with automatic step-control are constructed for systems of ODEs of ﬁrst and second order. The methods have fewer stages than classic Runge–Kutta methods.

The groups of equations (3) and (4) are structurally identical. Equations in them are ordered in such way that in every equation the right-hand side depends only on the previous functions of the same group. The group (2) contains all the equations that cannot be placed into (3) or (4), so it can be called general group. Sometimes it can be missing, as well as the group (3).
In works [6,9,10] a generalisation of explicit Runge-Kutta methods (RKs) for structurally partitioned systems (SPS) of ODEs (2)-(4) was suggested. The effectiveness of the methods presented there is provided by the fact, that the number of stages m for the general group (3) is equal to that for classic RKs (m = Q, q ≤ 4), but for the groups (3) and (4), taking in account their special structure, it can be reduced while keeping the method's order. Moreover, when the general group is absent the number of necessary stages can be reduced even more significantly. Thus, a fifth order method for a system with only groups (3) and (4) was constructed with four stages [9] (while classic RKs have at least six).
The considered generalisation of RKs is based on usage of the system's special structure in the algorithm. It is natural to apply the idea to embedded RKs. Dormand-Prince type methods [2] are widely used nowadays. Combining their approach with the First Same As Last (FSAL) idea -the last stage of a step is the first of the next step [3,4] -fast explicit one-step methods with stable step-size control were constructed. For systems with groups (3) and (4) an embedded Dormand-Prince type method of order five with automatic step-size control was constructed with five stages [11].

The integration method
In the paper an explicit embedded fourth order method with automatic stepsize control for systems y i = f i (x, y 1 , . . . , y i−1 , y l+1 , . . . , y n ), i = 1, . . . , l, y j = f j (x, y 1 , . . . , y j−1 ), j = l + 1, . . . , n, n ≥ 3 is constructed. Such systems appear, for example, in problems of optimal control and stabilisation [7,8], medical modelling [13], celestial mechanics or high-energy physics. For instance, large-scale oscillations of systems with rotational symmetry [6] are described (after some reducing) with equations Here x is a dimensionless moment of inertia about the rotational axis; z is a dimensionless moment of inertia about the equatorial plane; K is a kinetic energy of motions parallel to the equatorial plane; and L is a kinetic energy of vertical motions. Change of variables ξ 1 = x , ξ 2 = z , ξ 3 = x, ξ 4 = z, ξ 5 = K, ξ 6 = L reduces the system to the form (5), which is an SPS without the general group. The approximated solution of (5) is found as where k sv ≡ k sv (h) are calculated in the strict order k 11 , k 21 , . . . , k n1 , k 12 , k 22 , . . . , The method will be referred as structural.
For fast and effective practical implementation of numerical methods for ODEs one needs a step-size control algorithm. Here the embedded Dormand-Prince type method of order four, corresponding to the form (7)-(8), is constructed on base of three-stage fourth order methods from [6]. In Dormand-Prince methods two approximations of different orders are found simultaneously; the higher order result is used in further computations, and the lower order result is used only for the step-size control (technically, it is not the error being estimated, but the last term equal in the Taylor series of the solution and the approximation).
In [6] a number of important relations between the method's parameters was obtained. They combine the parameters a spvw , b sp and c sp into quite few groups and can be referred as group equalities: Here and further the indices i, j,î, of the parameters c iv , c jv , b iv , b jv , a ivîη , a ivjw , a jviw , a jvw correspond to the group attribute, not to a certain number. Namely, the index at first place connects parameters to the system (5) structural groups -i to the first and j to the second; and the third index shows which argument of f i or f j in (8) the parameter corresponds -i andî up to l-th function, j and from l + 1 to n.
The use of restrictions (9) significantly simplifies the design of order conditions and methods construction and realisation.
We use the slightly modified Butcher tableau to present the embedded structural method (see table 1).
The control term estimation is done with the same k iv and k jv values as the result but with d iv and d jv parameters The parameters of an m = (m i , m j )-stage embedded method (7)-(8) of order p and estimator order q must provide that y s (x + h) − z s = O(h p+1 ) and y s (x + h) −z s = O(h q+1 ) for all s = 1, . . . , n.
We consider methods where p > q. The difference z s −z s is used for the step-size control.
The method's (7)-(8), (10) abbreviation RKSp(q)(m i , m j )F or RKSp(q)mF tells about its type -RK (Runge-Kutta); its generalisation to systems with special structure (5) -S; its main order p and the estimator order (q); number of stages m if m i = m j = m or (m i , m j ); and about the FSAL idea used -F.
We construct fourth order methods with third order estimators. FSAL idea is also used.
New methods will be compared to the classic four-stages Runge-Kutta method ("The" Runge-Kutta method) with five-stages third-order estimator, named RK4(3)T in [2]; and to the same Runge-Kutta method but with an embedded second-order estimator (named in Russian mathematical tradition "Egorov control term" after Prof. Vsevolod A. Egorov [1]) -RK4(2).
The constructed method has better number of stages to order ratio than RK4(3)T and RK4 (2).
Often the implementation of Runge-Kutta methods is done so, that the step-size is changed every step. If so, k j4 values of RKS4(3)4F cannot be used as k j1 of the next step. It reduces the effectiveness of the method, but still it excels both RK4(3)T and RK4 (2).
There is an option to control the step-size on base of error estimation of only first l functions (which can work fine when two groups correspond to two different physical values or are scaled differently [1]). The fourth stage for the second group of (5) can be omitted then. Such method RKS4(4)(4,2)F (Table 3) in fact takes only three stages per step, while RK4(3)T and RK4(2) both take four.