On Reusing the Stages of a Rejected Runge-Kutta Step

: Runge-Kutta (RK) pairs are amongst the most popular methods for numerically solving Initial Value Problems. While using an RK pair, we may experience rejection of some steps through the interval of integration. Traditionally, all of the evaluations are then dropped, and we proceed with a completely new round of computations. In this work, we propose avoiding this waste and continuing by reusing the rejected RK stages. We focus especially on an RK pair of orders six and ﬁve. After step rejection, we reuse all the previously evaluated stages and only add three new stages. We proceed by evaluating the output using a smaller step. By this technique, we manage to signiﬁcantly reduce the cost in a set of problems that are known to pose difﬁculties in RK algorithms with changing step sizes.


Introduction
A fundamental issue in mathematical analysis, the Initial Value Problem (IVP) has applications in a variety of disciplines, including physics, engineering, and economics. The IVP has the following form with f : R × R m → R m continuously differentiable. The IVP's purpose is to determine the solution to a differential equation given the starting condition, which is the value of the dependent variable at a specified initial time. The most frequent approach for solving the IVP is to use numerical methods, which entail estimating the solution of the differential equation using discrete time steps.
The Runge-Kutta (RK) method is one of the most prominent and commonly used numerical methods for solving the IVP. The RK technique was developed in the late 1800s by German mathematicians Carl Runge and Martin Kutta. It is a class of numerical methods that estimate the value of the solution at the end of a time step by taking a weighted average of function values at different points in the time step. The algebraic order of the RK method relates to the number of function evaluations for each time step, which impacts the technique's accuracy.
For solving the IVP, the RK approach offers significant benefits over other numerical methods. It is a basic and easy-to-implement approach that is understandable to nonmathematicians. It is also a flexible approach that may be used to solve a broad variety of differential equations, including stiff problems that are difficult to solve with other numerical methods. Furthermore, the RK approach is a self-starting technique, which implies that no extra equations or conditions must be solved.
However, the RK method's accuracy is dependent on the sequence of the process, which increases the technique's computing cost and complexity. To overcome this issue, the Runge-Kutta pair, a version of the RK technique, was devised. A Runge-Kutta pair is a pair of RK methods using the same body of function evaluations but distinct weight coefficients that are used to estimate the solution and error at each time step. The estimated error is utilised to adjust the step size in the following time step, leading to a more accurate and efficient technique than the single RK method.
The general form of an s-stage Runge-Kutta pair of orders p and p − 1 can be written as follows (remark that clearly s > p for p > 4): where y n is the approximate solution at time x n , h n is the time step size, k ni are the intermediate function values, and b i ,b i , c i , and a ij are the RK coefficients. We assume that the coefficients form the column vectors b T ,b T c ∈ R s and the matrix A ∈ R s×s . The error estimate can be calculated as follows: By comparing the error estimate to a specified tolerance level TOL and modifying the time-step size accordingly, the step size may be adjusted. This enables a more efficient and precise technique for resolving the IVP. For achieving this, we use the formula In case e n+1 > TOL, the step is rejected, i.e., we do not move forward from x n to x n+1 . Then, we again use (3) but replace the left-hand side with h n . This is a classical procedure: see [1] (p. 167).
There have been several studies on the performance and accuracy of Runge-Kutta pairs in solving the IVP. Among the first such studies, Fehlberg compared the performance of different RK pairs of order 5 and 6 on a set of test problems and found that a particular RK pair, known as the Fehlberg pair, was the most efficient and accurate pair for a wide range of problems [2]. These pairs by Fehlberg are quadrature defective. Thus, for quadratic problems, the lower and higher order results coincide, and error control fails.
In another study, Prince and Dormand introduced a new RK pair of order 6 and 5 that avoided this defect and showed improved performance compared to other RK pairs [3].
In recent years, there have been developments in using RK pairs in combination with adaptive time-stepping algorithms to solve the IVP. One such algorithm is the Dormand-Prince method, which uses an embedded RK pair of order 5 and 4 with adaptive time-step size control to improve the accuracy and efficiency of the method [4]. Another algorithm is the Tsitouras method, which uses an embedded RK pair of order 5 and 6 with adaptive time-step size control; it has been shown to be more efficient and accurate than other RK pairs on certain types of problems [5].
In conclusion, the Runge-Kutta pair is a powerful and versatile numerical method for solving the IVP. It provides a more accurate and efficient method than using a single RK method by using a pair of methods of the same order with different coefficients to estimate the value of the solution and the error at each time step. There have been several studies on the performance and accuracy of RK pairs, with some pairs designed specifically for stiff equations and others optimized for efficiency and accuracy on a wide range of problems. Adaptive time-stepping algorithms, such as the Dormand-Prince and Tsitouras pairs, have been developed to further improve the accuracy and efficiency of the RK pair method.

The New Step-Size Control Scheme
After using the devised Equations (2) and (3) for changing the step size, we may experience step rejection. Then, traditionally, all the stages are completely wasted and we proceed by evaluating a new smaller step. This unjustified loss will be discussed here. We intend to suggest extension of existing pairs in order to handle these stages of a rejected step again and reduce the cost.
Thus, after an s-stage RK pair is applied and the error estimation fails to stay below n+1 , we reduce the current step to have a new length, say τh n , τ < 1. We then extend the pair by adding some more stages. After having at hand the s stages of the underlying pair, given by (1), we additionally compute and combine all theses stages in a new pair of formulas i.e., the new weights b * ,b * are used in order to approximate y * n+1 ≈ y(x n + τh n ). The error estimate now becomes and the slightly modified formula is used for advancing the integration. In consequence, we spend only (s − s) stages after a step rejection instead of the s-stages lost with the standard step-size-control algorithm. Obviously, we clearly demand (s − s) < s.
For safety reasons, we only apply the new algorithm in case for some λ > 1. Otherwise, we continue with the standard policy (2)-(3). We assume that the new extended set of coefficients now forms the column vectors b * T ,b * T c ∈ R˜s and the new extended matrix A ∈ R˜s ×s .

Evaluation of a Mew-Extended Runge-Kutta Pair
In the present work, we concentrate on a pair of orders six and five. We choose the DLMP6(5) pair that appeared in [6]. This pair seems to be among the best choices for moderate tolerances (i.e., 10 −4 > TOL > 10 −9 , see [5]). For this case, s = 9, but since a 9j = b j , j = 1, 2, · · · , 9, the last stage of each step can be used as the first stage of the next step. This is called FSAL (First Stage As Last), and then DLMP6(5) only actually spends 8 stages per step. The same number of stages is wasted after a step rejection using (2) and (3).
We set τ = 0.8 ands = 12. The new extended pair (named DLMP6(5)ext) shares an enhanced set of coefficients, as shown in the previous subsection. When a step is rejected, we evaluate these three new stages, k n,10 , k n,11 and k n,12 , and combine them with the new set of weights b * ,b * to form the new approximations for the point x n + τh n .
For solving the above equations, we make some simplifying assumptions. Namely and b · (A + C − τIs) = 0s.
In the above, we use the symbolization v (i+) with the elements of some vector v but with its first (i − 1) elements dropped. We set C = diag(c), and by Is and 0s, we mean, respectively, the identity and zero matrix in R˜s ×s .
Then, there is a severe reduction in the number of the order conditions to be solved. We are able to analytically solve these remaining order conditions with respect to the new coefficients. Indeed, when the above assumptions hold, only the order conditions have to be solved. The task is becoming a lot easier. We even manage to solve the seventh-order equations of condition for the higher-order formula. The resulting pair is shown in Table 3. These coefficients can be retrieved in Mathematica [7] format from http://users.uoa.gr/~tsitourasc/dlmp65ext.m (accessed on 3 June 2023).
An interpolant of the sixth-order at the same cost can be constructed for this type of pair [8]. Then, we may freely choose the magnitude of τ and vary the lengths of the new step after rejection. By fixing τ = 0.8 and achieving locally seventh-order accuracy after rejection, we obtain better results overall.

Numerical Results
We have chosen the detest problems D4, D5, and E2 as the set to perform our tests on. Integrating these problems by any 6(5) pair, we observe that almost 25% of the computational cost is wasted upon rejected steps. The problems chosen are given below.
This information is perhaps more significant because it ensures that we do not simply decrease the number of rejected steps with a very conservative step-size algorithm, which, in reverse, decreases efficiency. Table 4. Results of DLMP6(5) over the problem D4 using the standard step-size-change algorithm.
The corresponding results are given in Tables 8 and 9.
Interpreting the results, we observe about 28%-29% improvement in the efficiency for the problems and accuracies chosen. This can be found by dividing the corresponding values in the last column of Tables 4-11: e.g., for problem D4 and TOL = 10 −4 , we observe that the standard algorithm ((2) and (3)) is about 213. 6 150.3 = 42.1% more expensive than the new algorithm. In all 24 runs, we experience better efficiency with the new algorithm. Another interesting issue is that the error is decreasing, since after step rejection, we apply locally a seventh-order approximation.
The new algorithm proposed in the paper could be applied to fractional-order differential equations after some proper modifications [13,14].

Conclusions
A new step-control algorithm was presented in this article that, after the rejection of a Runge-Kutta step, avoids the unnecessary waste of the already-made computations. Numerical results on problems with a high percentage of rejected steps show the clear advantages of the new step-size policy.

Conflicts of Interest:
The authors declare no conflict of interest.