Adaptive Stepsize Control for Extrapolation Semi-Implicit Multistep ODE Solvers

: Developing new and efﬁcient numerical integration techniques is of great importance in applied mathematics and computer science. Among the variety of available methods, multistep ODE solvers are broadly used in simulation software. Recently, semi-implicit integration proved to be an efﬁcient compromise between implicit and explicit ODE solvers, and multiple high-performance semi-implicit methods were proposed. However, the computational efﬁciency of any ODE solver can be signiﬁcantly increased through the introduction of an adaptive integration stepsize, but it requires the estimation of local truncation error. It is known that recently proposed extrapolation semi-implicit multistep methods (ESIMM) cannot operate with existing local truncation error (LTE) estimators, e.g., embedded methods approach, due to their speciﬁc right-hand side calculation algorithm. In this paper, we propose two different techniques for local truncation error estimation and study the performance of ESIMM methods with adaptive stepsize control. The ﬁrst considered approach is based on two parallel semi-implicit solutions with different commutation orders. The second estimator, called the “double extrapolation” method, is a modiﬁcation of the embedded method approach. The introduction of the double extrapolation LTE estimator allowed us to additionally increase the precision of the ESIMM solver. Using several known nonlinear systems, including stiff van der Pol oscillator, as the testbench, we explicitly show that ESIMM solvers can outperform both implicit and explicit linear multistep methods when implemented with an adaptive stepsize.


Introduction
Many real-world processes can be mathematically described by systems of ordinary differential equations. One of the common ways to implement the solution of ODEs in discrete computers is numerical integration. The increasing complexity and scale of simulated systems requires a corresponding increase in the efficiency of ODE solvers. Thus, the development of novel numerical integration methods is of great interest in applied mathematics. In the last decade, many highly efficient solvers were proposed, including Falkner-type block methods [1,2], symmetric multistep methods [3], Störmer-Cowell methods [4], and single-step composition schemes [5]. Great attention is usually paid to the convergence and stability of designed methods, especially in the case of partial differential equations solvers [6]. Long-term simulations are also of great interest in the field, and many highly efficient and stable multistep schemes were developed recently [7,8]. The latest studies show that semi-implicit methods provide a reasonable compromise between highly stable implicit solvers and computationally efficient explicit solvers [9]. Moreover, being implemented in composition [5] and splitting [10] schemes, semi-implicit methods possess some properties of geometric integrators [11], e.g., time-reversibility and symmetry [9]. It is also possible to construct highly efficient single-step ODE solvers using the Aitken-Neville extrapolation [12]. In paper [13], the authors proposed a novel extrapolation semiimplicit multistep method (ESIMM), which combines the strengths of multistep schemes with the benefits of the extrapolation solvers, using the semi-implicit method as a basic integrator. It has been shown that ESIMM methods outperform classical multistep solvers, such as Adams-Bashforth (AB), Adams-Moulton (AM) and backward differentiation formula (BDF) being implemented with a fixed integration stepsize. However, the most advanced approaches for solving ODEs usually involve adaptive stepsize techniques [1,4]. Conventional local truncation error (LTE) estimators, such as the embedded methods approach [14], appeared to be barely suitable for ESIMM due to its specific method of right-hand side (RHS) calculation. Each step of the ESIMM solver in the variable-stepsize version requires recalculating all extrapolation coefficients due to their explicit dependency on the stepsize values used in previous iterations. Therefore, a new efficient LTE estimator and step control algorithm are needed to perform the ESIMM with an adaptive stepsize efficiently. In this paper, we propose two new techniques of adaptive stepsize control for ESIMM solvers. We develop a formula for recalculating the method's coefficients for an arbitrary sequence of stepsizes and propose two algorithms for local truncation error estimation. The first LTE estimator is based on the difference between two parallel solutions with different commutation orders [15] and the second, more efficient approach uses the double-extrapolation technique. Simulating chaotic systems is of great interest in applied mathematics due to the complex behavior and high sensitivity of their finitedifference models to the truncation error accumulation. Thus, we have chosen Rössler [16], Dadras-Momeni [17], Nose-Hoover [18][19][20] and van der Pol nonlinear systems [21] as test ODEs. In several computational experiments, we have shown that the ESIMM methods with adaptive stepsizes outperform both explicit and implicit linear multistep solvers, while simulating dissipative and conservative chaotic systems, as well as stiff nonlinear equations. The rest of the paper is organized as follows: in Section 2, the basic architecture of ESIMM methods is recalled, and a formula to recalculate the coefficients of the method is given. Section 3 introduces the novel LTE estimators and step control technique for ESIMM methods and provides the experimental estimation of the computational efficiency of the proposed method using a test set of dynamical systems. Finally, in Section 4, some conclusions are given.

Materials and Methods
Let us revisit the general idea of ESIMM solvers. These multistep methods use an extrapolation of the values obtained on previous solution steps together with newly calculated values, combining the multistep and extrapolation approaches to increase the accuracy order of the resulting scheme.
Let H 1 be a stepsize used to calculate the solution for IVP of order p at the point x (n+1) from the point x (n) . Then, the sequence of local error terms e is as follows: In order to obtain the solution at the point x (n+1) by calculating the increment from the point x (n−1) we will use the step value H 2 = H 1 + H n , where H n is a step value which steps from point x (n−1) to x (n) (see Figure 1). Mathematics 2021, 9,950 2 of 18 single-step ODE solvers using the Aitken-Neville extrapolation [12]. In paper [13], the authors proposed a novel extrapolation semi-implicit multistep method (ESIMM), which combines the strengths of multistep schemes with the benefits of the extrapolation solvers, using the semi-implicit method as a basic integrator. It has been shown that ESIMM methods outperform classical multistep solvers, such as Adams-Bashforth (AB), Adams-Moulton (AM) and backward differentiation formula (BDF) being implemented with a fixed integration stepsize. However, the most advanced approaches for solving ODEs usually involve adaptive stepsize techniques [1,4]. Conventional local truncation error (LTE) estimators, such as the embedded methods approach [14], appeared to be barely suitable for ESIMM due to its specific method of right-hand side (RHS) calculation. Each step of the ESIMM solver in the variable-stepsize version requires recalculating all extrapolation coefficients due to their explicit dependency on the stepsize values used in previous iterations. Therefore, a new efficient LTE estimator and step control algorithm are needed to perform the ESIMM with an adaptive stepsize efficiently. In this paper, we propose two new techniques of adaptive stepsize control for ESIMM solvers. We develop a formula for recalculating the method's coefficients for an arbitrary sequence of stepsizes and propose two algorithms for local truncation error estimation. The first LTE estimator is based on the difference between two parallel solutions with different commutation orders [15] and the second, more efficient approach uses the double-extrapolation technique. Simulating chaotic systems is of great interest in applied mathematics due to the complex behavior and high sensitivity of their finite-difference models to the truncation error accumulation. Thus, we have chosen Rössler [16], Dadras-Momeni [17], Nose-Hoover [18][19][20] and van der Pol nonlinear systems [21] as test ODEs. In several computational experiments, we have shown that the ESIMM methods with adaptive stepsizes outperform both explicit and implicit linear multistep solvers, while simulating dissipative and conservative chaotic systems, as well as stiff nonlinear equations. The rest of the paper is organized as follows: in Section 2, the basic architecture of ESIMM methods is recalled, and a formula to recalculate the coefficients of the method is given. Section 3 introduces the novel LTE estimators and step control technique for ESIMM methods and provides the experimental estimation of the computational efficiency of the proposed method using a test set of dynamical systems. Finally, in Section 4, some conclusions are given.

Materials and Methods
Let us revisit the general idea of ESIMM solvers. These multistep methods use an extrapolation of the values obtained on previous solution steps together with newly calculated values, combining the multistep and extrapolation approaches to increase the accuracy order of the resulting scheme.
Let H1 be a stepsize used to calculate the solution for IVP of order p at the point x(n+1) from the point x(n). Then, the sequence of local error terms e is as follows: In order to obtain the solution at the point x(n+1) by calculating the increment from the point x(n−1) we will use the step value  Now the sequence of local error terms can be calculated as follows: Mathematics 2021, 9, 950 3 of 14 The further algorithm is described in [13] and generally follows the idea of the extrapolation methods [11,22] in multistep form.
The main idea of the ESIMM method is to use the solution with the stepsize H 1 (denote it as T 1 ) and the solution at the same point with the stepsize H 2 (denote it as T 2 ) in the following way: Coefficients k 1 and k 2 are chosen in accordance with general extrapolation principles, taking into account the ratio of truncation errors in adjacent stages so that the following conditions are satisfied: where H 1 and H 2 are the step values used to calculate the solutions at point x (n+1), x (n) and x (n−1) , respectively. Figure 2 illustrates the general scheme of the ESIMM method of order 3. Now the sequence of local error terms can be calculated as follows: The further algorithm is described in [13] and generally follows the idea of the extrapolation methods [11,22] in multistep form.
The main idea of the ESIMM method is to use the solution with the stepsize H1 (denote it as T1) and the solution at the same point with the stepsize H2 (denote it as T2) in the following way: Coefficients k1 and k2 are chosen in accordance with general extrapolation principles, taking into account the ratio of truncation errors in adjacent stages so that the following conditions are satisfied: where H1 and H2 are the step values used to calculate the solutions at point x(n+1), x(n) and x(n−1), respectively. Figure   By using matrix formalism, one can obtain coefficients for the first stage of extrapolation according to Equation (3). Because we have to eliminate the error term ep+1 in both (1) and (2), the coefficients are as follows: ; 0 After completing the first extrapolation step, the equations for two of our approximations for the point x(tn+1), which we will call T12 and T22, are as follows: By using matrix formalism, one can obtain coefficients for the first stage of extrapolation according to Equation (3). Because we have to eliminate the error term e p+1 in both (1) and (2), the coefficients are as follows: After completing the first extrapolation step, the equations for two of our approximations for the point x(t n+1 ), which we will call T 12 and T 22 , are as follows: The vector form can be used to rewrite those expressions as follows: Using a similar approach, one can obtain the final equations for k 12 and k 14 : Using this technique, one can obtain extrapolation coefficients for a "full" variation of the ESIMM method, which requires solving the matrix equation at every step. It is possible to significantly reduce the number of calculations using the so-called "short" version of the same method, previously described in [13]. We will use the "short" ESIMM method throughout this paper.
Let us modify a short ESIMM method to make it suitable for implementation with a variable stepsize by interpreting our step values as the following system of linear algebraic equations: Considering the Taylor expansion for k i* x(t n+1 ) through T i where T i represents the first components of the extrapolation table, we obtain coefficients for the short version of the ESIMM method of the desired order: Keeping in mind that that the sum of coefficients k i should equal 1, while the sum of all the other terms by any power of H equals zero, one can write the following system of algebraic equations: Mathematics 2021, 9, 950 5 of 14 This approach can be used for obtaining coefficients of the ESIMM method of an arbitrary accuracy order and proved to be efficient while using a symmetric basic method. In the current study, we use the semi-implicit self-adjoint CD method [9,22] of order 2 as a basic method for the ESIMM solver.
The important note here is that one needs to recalculate all extrapolation coefficients on each step of ESIMM while implementing it with variable stepsizes because these coefficients are directly influenced by the step values used in previous iterations. Thus, we recommend using the "short" version of the ESIMM method in practical applications because it is more computationally efficient, provides fewer overhead calculations, and maintains the same high accuracy as the "full" version.

Results
In this section, we propose two algorithms for estimating the local truncation error while controlling the stepsize of the ESIMM solver. Then, we analyze the performance of the extrapolation semi-implicit multistep (ESIMM) method with a variable stepsize in comparison with classical multistep methods, including Adams-Bashforth, Adams-Moulton and Backward Differentiation Formula. In our study, we use the Dormand-Prince 8 (DOPRI8) method [23] to obtain the reference solution. All experiments were conducted using NI LabVIEW 2020 environment with double floating-point data precision.

Rössler Attractor
This dissipative chaotic system was proposed by Otto E. Rössler [16] and can be described as follows: where a, b and c are system parameters. The basic semi-implicit CD method for the system is as follows: y n+0.5 = y n + H 2 (x n + ay n ); z n+0.5 = z n + H 2 (b + z n x n − cz n ); x n+1 = x n + H(−y n+0.5 − z n+0.5 ); We carried out experiments with the following parameter values: a = 0.2, b = 0.2, c = 5.7. The timing parameters and initial conditions are given in Table 1. To implement the adaptive stepsize control, one needs to estimate the local truncation error (LTE) of the method. The first LTE estimator for the ESIMM method considered in this study is based on the difference between two parallel solutions with different commutation orders, generally following the idea used in [15] to create new metrics for quantifying chaotic dynamics. The term "commutations order" basically means changing the order of the state variables calculation, as shown in Equation (14), where a second possible commutation of CD method for the Rössler system is given: (14) Figure 4 shows the scheme of the LTE evaluation using two different commutations of the basic CD method executed in parallel.  The local truncation error is then estimated as the difference between ESIMM solutions obtained with basic methods (13) and (14). The step control law is as follow [14]: where m is the accuracy order of the scheme, Tol is the desired tolerance value and d(xn) i a local error between two solutions. Figure 5 illustrates the behavior of the stepsize in the ESIMM scheme with the commutation-based LTE estimator being compared with the explicit Adams-Bashforth linear multistep method with a standard stepsize control algorithm [14]. The local truncation error is then estimated as the difference between ESIMM solutions obtained with basic methods (13) and (14). The step control law is as follows [14]: where m is the accuracy order of the scheme, Tol is the desired tolerance value and d(x n ) is a local error between two solutions. Figure 5 illustrates the behavior of the stepsize in the ESIMM scheme with the commutation-based LTE estimator being compared with the explicit Adams-Bashforth linear multistep method with a standard stepsize control algorithm [14].  Step behavior for Adams-Bashforth and ESIMM methods of order 4 while solving the Rössler system. A commutation-based LTE estimator is used for controlling the stepsize in the ESIMM scheme.
One can see from Figure 5 that the ESIMM solver demonstrates a sharper curve of the stepsize dynamics graph, but at the same time, it tends to use higher step values, resulting in fewer RHS function calculations. It follows from the fewer number of steps required for obtaining a solution with ESIMM. Although these results were already promising, we additionally developed an alternative step control approach to further increase the ESIMM solver efficiency.
The second proposed technique for local truncation error estimation derives from Step behavior for Adams-Bashforth and ESIMM methods of order 4 while solving the Rössler system. A commutation-based LTE estimator is used for controlling the stepsize in the ESIMM scheme.
One can see from Figure 5 that the ESIMM solver demonstrates a sharper curve of the stepsize dynamics graph, but at the same time, it tends to use higher step values, resulting in fewer RHS function calculations. It follows from the fewer number of steps required for obtaining a solution with ESIMM. Although these results were already promising, we additionally developed an alternative step control approach to further increase the ESIMM solver efficiency.
The second proposed technique for local truncation error estimation derives from the embedded method approach, which is broadly used in single-step extrapolation solvers and embedded Runge-Kutta methods [14]. To estimate the local truncation error, a difference between two basic methods with the same commutation orders but with different stepsizes was used. One solution was obtained using the step value H and the other contains a composition of two consequent steps with H/2. Taking into account that the basic method is of order 2, the local truncation error demonstrates a quadratic dependence from the chosen integration step. Both solutions are extrapolated to increase the overall accuracy as follows: where P H n+1 and P H/2 n+1 are two solutions found with step sizes of H and H/2, respectively. We called this approach "double extrapolation" because the ESIMM method already includes one extrapolation procedure in the computation algorithm. After acquiring two ESIMM solutions, one can use the difference between them as the estimation of LTE and calculate the new stepsize value using Equation (15) . Figure 6 illustrates the modification of the basic method in the ESIMM scheme, incorporating a double extrapolation estimator. Figure 7 illustrates the benefits of the double extrapolation approach in the form of performance plots obtained during the Rössler system simulation. The performance plots show the dependence between the achieved precision and the real program execution time needed to obtain the solution with the desired tolerance. The execution time was averaged from 100 simulation runs. The simulation parameters were set in accordance with Table 1.
As is shown in Figure 7, using the double-extrapolation technique results in a significant increase in the calculation speed on higher tolerance values and a slight improvement in the accuracy of the solution on low tolerances. Note that the stepsize curve in Figure 7b is smoother than that in Figure 7a. Figure 8 represents the performance plots for all methods under investigation of accuracy order 3 and 5 while solving the Rössler system.  As is shown in Figure 7, using the double-extrapolation technique results in a significant increase in the calculation speed on higher tolerance values and a slight improvement in the accuracy of the solution on low tolerances. Note that the stepsize curve in Figure 7b is smoother than that in Figure 7a. Figure 8 represents the performance plots for all methods under investigation of accuracy order 3 and 5 while solving the Rössler system. One can see from Figure 8 that ESIMM methods with adaptive stepsizes appear to be more computationally efficient than other multistep methods under investigation, namely, Adams-Bashforth, Adams-Moulton and Backward Differentiation Formula. The observed differences in performance with explicit methods, such as Adams-Bashforth One can see from Figure 8 that ESIMM methods with adaptive stepsizes appear to be more computationally efficient than other multistep methods under investigation, namely, Adams-Bashforth, Adams-Moulton and Backward Differentiation Formula. The observed differences in performance with explicit methods, such as Adams-Bashforth solver, become less noticeable with the increase in the accuracy order, but are still significant in terms of accuracy and calculation speed on lower tolerance values. This fact can be explained by the extra RHS calculations introduced to the ESIMM scheme by both proposed LTE estimation algorithms. However, the higher stability of ESIMM methods still provides some superiority over Adams-Bashforth solver on lower tolerance values. Let us consider some other test dynamical systems to verify the obtained results in a series of computational experiments.
We carried out experiments with the Dadras-Momeni system using the following parameter values: a = 3, b = 2.7, c = 4.7, d = 2, m = 9. Other simulation parameters are given in Table 2. Step Values (s.)

Maximum 1
One can see from Figure 9 that the efficiency of the ESIMM methods decreases with the increase in accuracy order due to the necessity of calculating additional RHS functions per order. Nevertheless, the ESIMM solver with the double-extrapolation LTE estimator provides the best computational efficiency among the investigated methods. It is of certain interest to estimate the performance of the proposed multistep methods while simulating a conservative nonlinear system.

Nose-Hoover Attractor
This conservative chaotic system was proposed by S. Nose and W.G. Hoover [18][19][20] and can be formulated as follows: 2 , dx ay dt dy x yz dt dz d y dt where a and d are system parameters. The basic semi-implicit algorithm for the system is as follows:

Nose-Hoover Attractor
This conservative chaotic system was proposed by S. Nose and W.G. Hoover [18][19][20] and can be formulated as follows: where a and d are system parameters. The basic semi-implicit algorithm for the system is as follows: Step Values (s.) From Figure 10, one can see that while simulating a conservative system, ESIMM methods maintain good performance even in the case of higher-order schemes. This fact can possibly be explained by the geometric properties of the basic CD method, which provides extra precision while solving conservative ODEs. However, the proposed step control algorithms do not possess energy-preserving properties, and this phenomenon still is to be investigated. Note that other linear multistep methods under investigation appear to be much less accurate compared with results obtained for dissipative systems, while ESIMM gives promising results in both cases.

Van der Pol System
Considering the possibility that the abovementioned differences in the methods' performance can originate from the varied stiffness of the previously studied systems, we considered a well-known stiff dynamical system as a test ODE. This nonlinear oscillator is a traditional testbench for numerical integration methods and was proposed by Balthasar van der Pol [21]:

Van der Pol System
Considering the possibility that the abovementioned differences in the methods' performance can originate from the varied stiffness of the previously studied systems, we considered a well-known stiff dynamical system as a test ODE. This nonlinear oscillator is a traditional testbench for numerical integration methods and was proposed by Balthasar van der Pol [21]: where m is a parameter that defines the stiffness of the system. The basic semi-implicit CD algorithm for the van der Pol oscillator is as follows: y n+0.5 = y n + H 2 (m(1 − x n x n )y n − x n ); x n+1 = x n + Hy n+0.5 ; Experiments were carried out with m = 55 which corresponds to the medium-stiff scenario. The rest of the simulation parameters are presented in Table 4. One can see from Figure 11 that the proposed adaptive ESIMM methods perform better than the other investigated multistep methods while simulating the stiff nonlinear system. The excellent numerical stability of the ESIMM scheme, analytically determined in [13], resulted in a better overall performance in the stiff system case. Moreover, it is possible to further increase the stability of the ESIMM solver by changing the complexity of the basic method and, therefore, achieving a perfect balance between the performance and stability of the solver. system. The excellent numerical stability of the ESIMM scheme, analytically determined in [13], resulted in a better overall performance in the stiff system case. Moreover, it is possible to further increase the stability of the ESIMM solver by changing the complexity of the basic method and, therefore, achieving a perfect balance between the performance and stability of the solver.

Conclusions and Discussion
This study considered two possible techniques to implement extrapolation semi-implicit multistep ODE solvers with adaptive stepsizes. The first method, which uses the idea of estimating the local truncation error as a difference between two parallel solutions with different string commutations, proven to be a reliable and computationally cheap step control algorithm but, in some cases, it failed to outperform explicit multistep methods. The second approach, named the double extrapolation estimator, is based on a more traditional technique of LTE evaluation, which is similar to the embedded Runge-Kutta methods or embedded LTE estimators in extrapolation solvers. A double extrapolation ESIMM solver assumes the parallel solution consisting of the composition of two basic methods performed with a half-step. This approach allowed us to introduce an extra extrapolation procedure per stage to the ESIMM method, thereby increasing its accuracy. The double extrapolation technique not only allowed estimating the local error more precisely but also additionally increased the overall performance of the scheme. The numerical experiments included the simulation of four non-linear ODEs: the dissipative Rössler and Dadras-Momeni systems, the conservative Nose-Hoover system and the classical stiff van der Pol stiff oscillator. ESIMM methods with the proposed adaptive step control techniques appeared to be the most efficient solvers among the tested integration schemes. Further studies will be devoted to the investigation of ESIMM methods with different basic methods, e.g., semi-implicit SED solver, semi-implicit midpoint methods and symmetrized explicit midpoint method. We will also consider the possibility of applying the proposed step control techniques to some other numerical methods, including PDE solvers.

Conclusions and Discussion
This study considered two possible techniques to implement extrapolation semiimplicit multistep ODE solvers with adaptive stepsizes. The first method, which uses the idea of estimating the local truncation error as a difference between two parallel solutions with different string commutations, proven to be a reliable and computationally cheap step control algorithm but, in some cases, it failed to outperform explicit multistep methods. The second approach, named the double extrapolation estimator, is based on a more traditional technique of LTE evaluation, which is similar to the embedded Runge-Kutta methods or embedded LTE estimators in extrapolation solvers. A double extrapolation ESIMM solver assumes the parallel solution consisting of the composition of two basic methods performed with a half-step. This approach allowed us to introduce an extra extrapolation procedure per stage to the ESIMM method, thereby increasing its accuracy. The double extrapolation technique not only allowed estimating the local error more precisely but also additionally increased the overall performance of the scheme. The numerical experiments included the simulation of four non-linear ODEs: the dissipative Rössler and Dadras-Momeni systems, the conservative Nose-Hoover system and the classical stiff van der Pol stiff oscillator. ESIMM methods with the proposed adaptive step control techniques appeared to be the most efficient solvers among the tested integration schemes. Further studies will be devoted to the investigation of ESIMM methods with different basic methods, e.g., semi-implicit SED solver, semi-implicit midpoint methods and symmetrized explicit midpoint method. We will also consider the possibility of applying the proposed step control techniques to some other numerical methods, including PDE solvers.