Semi-Implicit Multistep Extrapolation ODE Solvers

: Multistep methods for the numerical solution of ordinary di ﬀ erential equations are an important class of applied mathematical techniques. This paper is motivated by recently reported advances in semi-implicit numerical integration methods, multistep and extrapolation solvers. Here we propose a novel type of multistep extrapolation method for solving ODEs based on the semi-implicit basic method of order 2. Considering several chaotic systems and van der Pol nonlinear oscillator as examples, we implemented a performance analysis of the proposed technique in comparison with well-known multistep methods: Adams–Bashforth, Adams–Moulton and the backward di ﬀ erentiation formula. We explicitly show that the multistep semi-implicit methods can outperform the classical linear multistep methods, providing more precision in the solutions for nonlinear di ﬀ erential equations. The analysis of stability regions reveals that the proposed methods are more stable than explicit linear multistep methods. The possible applications of the developed ODE solver are the long-term simulations of chaotic systems and processes, solving moderately sti ﬀ di ﬀ erential equations and advanced modeling systems.


Introduction
In recent years, tools for computer simulation have become increasingly important in various fields of science. Computational experiments allow one to obtain new knowledge about the properties and dynamics of the object under investigation without the significant economic and time costs and risks that attend a full-scale experiment. To study the continuous dynamical systems described by ordinary differential equations (ODE), numerical integration methods are traditionally used. While the scale and complexity of the systems under investigation increase, one needs more sophisticated and efficient integration techniques [1]. Therefore, the advanced numerical methods used to simulate dynamical systems should possess great stability, precision and low computational costs. In engineering applications, especially while solving high-dimensional initial value problems (IVP), multistep integration is widespread [2].
Being compared to their single-step counterparts, multistep methods employ one call to the right-hand-side function per each integration step. To obtain the numerical solution with the desired accuracy, such algorithms reuse values calculated at previous steps. Therefore, multistep algorithms significantly exceed single-step integration in computational speed. However, the well-known problem of such integration is a decrease in numerical stability with an increase in the accuracy order of the chosen scheme. This property restricts the application of high-order multistep methods for solving

Materials and Methods
Consider a basic symmetric integration method for solving IVP of order p. A sequence of local error terms after stepping from the point x n to x n+1 with a stepsize h reads x(t n+1 ) = x n+1 + h p+1 e p+1 + h p+2 e p+2 + . . .
Following the idea of extrapolation methods [4,12], let us combine the solution at the point n + 1 with the stepsize h (denote it as T 1 ) and the solution at the same point with the stepsize 2h (denote it as T 2 ) as follows: 3 of 18 The dataflow representation of this combination is depicted in Figure 1. Keeping in mind Aitken-Neville interpolation rule, let us select the coefficients k 1 and k 2 so as to satisfy the following conditions: k 1 + k 2 = 1 k 1 + 2 p+1 k 2 = 0 (4) The system of Equation (4) can be analytically solved as This allows eliminating the term by the error ep+1 and increasing the order by one.  A general form for increasing an order can be obtained if we consider how the operation (3) affects the error terms in (1) and (2). Introducing the matrix formalism, for the coefficients on the first stage one can obtain After the first stage of extrapolation, equations for two approximations of ( ) denoted as and , which are then multiplied by and respectively, are as follows: = + ( ℎ + (2ℎ) ) + ⋯ Figure 1. Graphical interpretation of the proposed multistep algorithm: solution T 1 obtained with stepsize h from the point x n and solution T 2 obtained with stepsize 2h from the point x n−1 are combined to obtain a higher-order solution. Basic integration methods are denoted as Φ h and Φ 2h with respect to the stepsize.
The system of Equation (4) can be analytically solved as For example, for the method of order 3 with p = 2 we have: This allows eliminating the term by the error e p+1 and increasing the order by one. Let us organize the computational process in a cascade manner as shown in Figure 2. The main term T ij represents a component in the extrapolation table on i-th stage on the j-th line. Indices by coefficients k m ij are the following: m denotes whether it is the coefficient by the first or the second term, with the sum giving the term T ij .
The system of Equation (4) can be analytically solved as − For example, for the method of order 3 with p = 2 we have: This allows eliminating the term by the error ep+1 and increasing the order by one.  A general form for increasing an order can be obtained if we consider how the operation (3) affects the error terms in (1) and (2). Introducing the matrix formalism, for the coefficients on the first stage one can obtain After the first stage of extrapolation, equations for two approximations of ( ) denoted as and , which are then multiplied by and respectively, are as follows: A general form for increasing an order can be obtained if we consider how the operation (3) affects the error terms in (1) and (2). Introducing the matrix formalism, for the coefficients on the first stage one can obtain After the first stage of extrapolation, equations for two approximations of x(t n+1 ) denoted as T 12 and T 22 , which are then multiplied by k 1 13 and k 2 13 respectively, are as follows: . . Let us rewrite them using a vector form: (5), the final equation for k 13 reads:

Similar to the Equation
Continuing such a sequence of operations, obtain an equation for k 14 : Continuing this process, we can derive the coefficients for the method of an arbitrary accuracy order. For the basic symmetric method of order 2, the coefficients for the 6-order multistep extrapolation method are given in Table 1. Let us call the method illustrated by Figure 2 and Table 1, a "full" ESIMM method. However, one can use a simplified ("short") version of the same method: In order to obtain coefficients k i for the "short" version of the method, consider the Taylor expansion for k i x(t n+1 ) through T i1 up to the desired order of the error term s + p: k 1 x(t n+1 ) = k 1 T 11 + k 1 h p+1 e p+1 + k 1 h p+2 e p+2 + . . . + k 1 h p+s e p+s , k 2 x(t n+1 ) = k 2 T 21 + k 2 (2h) p+1 e p+1 + k 1 (2h) p+2 e p+2 + . . . + k 1 (2h) p+s e p+s , k s x(t n+1 ) = k s T s1 + k s (sh) p+1 e p+1 + k s (sh) p+2 e p+2 + . . . + k s (sh) p+s e p+s One can see, that that the sum of k i equals 1, while the sum of any other terms by any power of h equals zero, which can be written as a system of linear algebraic equations: from where k i can be easily found. Table 2 represents the coefficients for the "short" method version. One can obtain coefficients for methods of higher order following the general logic of this approach. It should be noted that this technique is efficient only when a self-adjoint symmetric basic method is chosen. In our study, we used the semi-implicit CD [11,12] method of order 2, but the implicit midpoint method can be used as a basic method as well.

Results
In this section, we investigate the performance of two versions of the extrapolation semi-implicit multistep (ESIMM) method. Let us recall, that the first version is labeled as "full" and obtains a solution using a direct calculation scheme with the full set of coefficients per each solution step. The second version, named "short", uses pre-calculated coefficients values for each stage and sums the weighted results, which results in an increased calculation speed at the cost of a possibly insignificant decrease in accuracy compared to the "full" version due to the round-off error. All multistep methods under investigation were used with the Dormand-Prince 8 (DOPRI8) method as a starting solver and reference solution. All the experiments were performed in an NI LabVIEW 2018 environment with double floating-point precision.
We have chosen several well-known nonlinear dynamical systems as a test set for the numerical method investigation.

Rössler Attractor
The system was proposed by Otto E. Rossler [16], and can be described as follows: where a, b and c are system parameters. The semi-implicit basic method of order 2 for the system (8) can be formulated as [11][12][13]: We carried out experiments with the following parameter values: Other simulation parameters, including simulation time, stepsize values, initial conditions and number of measurements are given in Table 3. Table 3. Parameters used with the experiments with the Rossler system.

Order of Accuracy 3 4 5 6
Step values, sec. One can see from Figure 3 that the proposed semi-implicit extrapolation multistep method (ESIMM) possesses better performance than well-known linear multistep methods. The "short" version of ESIMM provides a great balance between computational cost and accuracy, while the "full" version is slightly more precise but requires more time to perform calculations. Classical explicit multistep methods such as Adams-Bashforth tend to have closer performance in higher-order schemes but are still significantly less accurate compared to the proposed method.

Sprott System Case A and Case E
These two systems were proposed by J.C. Sprott in his famous paper [17] and can be written as following initial value problems: Sprott Case A is a conservative chaotic system and exhibits chaotic behavior with parameter values a = 1, b = 1. The basic semi-implicit symmetric algorithm for the Sprott Case A system is as follows:

Sprott System Case A and Case E
These two systems were proposed by J.C. Sprott in his famous paper [17] and can be written as following initial value problems: Sprott Case A is a conservative chaotic system and exhibits chaotic behavior with parameter values a = 1, b = 1. The basic semi-implicit symmetric algorithm for the Sprott Case A system is as follows: The semi-implicit basic method for the Sprott Case E system can be written as: Parameter value d = 11 corresponds to the chaotic regime of oscillations. Figure 4 represents the performance analysis for all investigated methods including the proposed ESIMM algorithms.  Thus, the numerical efficiency of the ESIMM algorithm decreases with the increase in accuracy order due to the one extra calculation of RHS function per order, unlike classical multistep methods.
The simulation of the Sprott Case E system confirmed the preliminary results. One can see from Figure 5 that the performance gain of ESIMM methods over the Adams-Bashforth method decreases as the order of the scheme increases. However, the overall results are inspiring due to the known fact that the application of Adams-Bashforth method with orders of accuracy 6 and higher is rare due to the stability issues. Table 5 gives some information on the simulation parameters.   The rest of the parameters used to perform the numerical experiments with the system (9) are given in Table 4. One can see from Figure 4 that the performance of the ESIMM methods in the case of the conservative chaotic system Sprott Case A was unexpectedly lower despite the symmetry of the underlying basic semi-implicit method. Step values, sec. Thus, the numerical efficiency of the ESIMM algorithm decreases with the increase in accuracy order due to the one extra calculation of RHS function per order, unlike classical multistep methods.
The simulation of the Sprott Case E system confirmed the preliminary results. One can see from Figure 5 that the performance gain of ESIMM methods over the Adams-Bashforth method decreases as the order of the scheme increases. However, the overall results are inspiring due to the known fact that the application of Adams-Bashforth method with orders of accuracy 6 and higher is rare due to the stability issues. Table 5 gives some information on the simulation parameters. By analyzing results shown in Figures 4 and 5, one can conclude that the semi-implicit ESIMM methods perform well with both conservative (Sprott case A) and dissipative (Sprott case E) systems. They are significantly faster than all other listed methods and give us a rather accurate representation of systems behavior. However, implicit multistep methods such as Adams-Moulton can be more accurate than ESIMM in higher-order schemes, but the proposed method proposes the best trade-off between precision and simulation time, which is extremely important in long-term simulations. To  Step values, sec. By analyzing results shown in Figures 4 and 5, one can conclude that the semi-implicit ESIMM methods perform well with both conservative (Sprott case A) and dissipative (Sprott case E) systems. They are significantly faster than all other listed methods and give us a rather accurate representation of systems behavior. However, implicit multistep methods such as Adams-Moulton can be more accurate than ESIMM in higher-order schemes, but the proposed method proposes the best trade-off between precision and simulation time, which is extremely important in long-term simulations. To investigate the impact of stiffness on the investigated methods, we additionally considered one more system that is well known as a testbench for numerical integration methods.

Van der Pol System
This classical nonlinear oscillator was proposed by Balthasar van der Pol [18] and can be written as the following IVP: This system is important for our study because the value of parameter m determines the stiffness of the system. The algorithm of the basic symmetric CD method for the system (11) is as follows: The simulation parameters are given in Table 6. The parameter m initially was set to 1, but for higher-order schemes we reduced it to 0.5 to avoid the loss of stability of explicit methods. We varied the initial conditions during the experiments to ensure that they do not significantly affect the obtained results. The performance plots for the van der Pol system are shown in Figure 6.
The simulation parameters are given in Table 6. The parameter m initially was set to 1, but for higher-order schemes we reduced it to 0.5 to avoid the loss of stability of explicit methods. We varied the initial conditions during the experiments to ensure that they do not significantly affect the obtained results. The performance plots for the van der Pol system are shown in Figure 6. One can see from Figure 6 that implicit linear multistep methods are slightly more accurate than other considered algorithms. Explicit methods such as Adams-Bashforth are, on the other hand, as fast as semi-implicit ones but they lose their accuracy due to less numerical stability. Further increase in m leads to total stability loss for Adams-Bashforth methods. We illustrated this in the final experiment with this system, the results of which are shown in Figure 7. One can see the higher stability of semi-implicit multistep methods in comparison with their well-known explicit One can see from Figure 6 that implicit linear multistep methods are slightly more accurate than other considered algorithms. Explicit methods such as Adams-Bashforth are, on the other hand, as fast as semi-implicit ones but they lose their accuracy due to less numerical stability. Further increase in m leads to total stability loss for Adams-Bashforth methods. We illustrated this in the final experiment with this system, the results of which are shown in Figure 7. One can see the higher stability of semi-implicit multistep methods in comparison with their well-known explicit counterparts. However, the Adams-Moulton method appears to be the most precise among the investigated algorithms, but the "short" version of the ESIMM solver generally outperforms it due to lower computational costs. Figure 6. The performance plots of van der Pol system obtained for investigated multistep methods with the accuracy (a) order 3, (b) order 4, (c) order 5 and (d) order 6.
One can see from Figure 6 that implicit linear multistep methods are slightly more accurate than other considered algorithms. Explicit methods such as Adams-Bashforth are, on the other hand, as fast as semi-implicit ones but they lose their accuracy due to less numerical stability. Further increase in m leads to total stability loss for Adams-Bashforth methods. We illustrated this in the final experiment with this system, the results of which are shown in Figure 7. One can see the higher stability of semi-implicit multistep methods in comparison with their well-known explicit counterparts. However, the Adams-Moulton method appears to be the most precise among the investigated algorithms, but the "short" version of the ESIMM solver generally outperforms it due to lower computational costs.

Order Plots and Convergence Test
To experimentally evaluate the convergence of the methods under investigation, we applied the approach known as order plots. This type of analysis involves the definition of global truncation error for a numerical method of order N, which can be generally described as E = O (h) N , where O is a system-dependent constant and h is integration stepsize [19]. In our experiment, the truncation error was estimated numerically through the simulation of the Rossler nonlinear oscillator by all methods under investigation with stepsize h and h/2. Then, order plots can be constructed as E(h)/E (h/2) vs. stepsize h. An order plot for an "ideal" numerical method of accuracy order 4 should be a horizontal line equal to 16, or for order 5 equal to 32, etc. Thus, the investigated finite-difference scheme can be considered more stable if its order plot is situated higher on the Y-axis and converges faster to the "ideal" line while stepsize decreases. One should note that this method evaluates the convergence and order preservation only empirically and for a certain test problem. We used the 2-norm of maximum truncation errors of all state variables of the system to estimate E(h) and E(h/2). Figures 8 and 9 show the order plots for all investigated solvers of order 4 and 5, respectively. We omitted the plot for the "full" ESIMM method because it possesses minimal differences from the "short" version in this experiment. One can see from Figure 8 that the proposed method outperforms both explicit and implicit linear multistep methods on the Rossler problem, being closer to the ideal method of order 4 (E(h)/E(h/2) = 16). Figure 9 represents real order evaluation for multistep methods of order 5. The left-hand part of the order plot can be used to roughly evaluate the convergence empirically. One can see that the ESIMM method quickly converges to the "ideal" order 5 (E(h)/E(h/2) = 32) while stepsize decreases. Other methods under investigation show similar convergence, but the explicit Adams-Bashforth method loses stability at a certain stepsize value. Thus, we can conclude that the proposed semi-implicit extrapolation methods can be more stable than explicit linear multistep methods, and possess higher practical order of accuracy for a given set of problems. This confirms the theoretically predicted properties of such methods. However, the analytical investigation of numerical stability is still needed, including the plotting of stability regions. maximum truncation errors of all state variables of the system to estimate E(h) and E(h/2). Figures 8 and 9 show the order plots for all investigated solvers of order 4 and 5, respectively. We omitted the plot for the "full" ESIMM method because it possesses minimal differences from the "short" version in this experiment. One can see from Figure 8 that the proposed method outperforms both explicit and implicit linear multistep methods on the Rossler problem, being closer to the ideal method of order 4 (E(h)/E(h/2) = 16). Figure 9 represents real order evaluation for multistep methods of order 5. The left-hand part of the order plot can be used to roughly evaluate the convergence empirically. One can see that the ESIMM method quickly converges to the "ideal" order 5 (E(h)/E(h/2) = 32) while stepsize decreases. Other methods under investigation show similar convergence, but the explicit Adams-Bashforth method loses stability at a certain stepsize value. Thus, we can conclude that the proposed semi-implicit extrapolation methods can be more stable than explicit linear multistep methods, and possess higher practical order of accuracy for a given set of problems. This confirms the theoretically predicted properties of such methods. However, the analytical investigation of numerical stability is still needed, including the plotting of stability regions.

Stability regions
To analytically evaluate the stability of ESIMM methods we will use a technique based on the known approach with 2-dimensional test problems [11,20]. Let us recall that the semi-implicit methods exist only for the initial value problems of dimension 2 or greater, and one needs a special technique here. Consider a 2-dimensional autonomous test problem = , = .
Let eigenvalues of the matrix be = ± .
Then, the contents of the matrix can be determined by two free parameters, denote them ≥ 0 (eccentricity coefficient) and ≥ 0 (asymmetry coefficient): The shape of the stability region depends on , reaching its maximal size when = 1, which corresponds to the Jordan normal form of the matrix, and minimal size when = 0 or = ∞, which corresponds to the Frobenius normal form of the matrix. The stability region for the 2-dimensional problem is the area in the complex plane where the absolute values of complex conjugate eigenvalues of the matrix ( ℎ) are not greater than 1.

Stability Regions
To analytically evaluate the stability of ESIMM methods we will use a technique based on the known approach with 2-dimensional test problems [11,20]. Let us recall that the semi-implicit methods exist only for the initial value problems of dimension 2 or greater, and one needs a special technique here. Consider a 2-dimensional autonomous test problem Let eigenvalues of the matrix A be λ 12 = σ ± jω.
Then, the contents of the matrix can be determined by two free parameters, denote them r ≥ 0 (eccentricity coefficient) and k ≥ 0 (asymmetry coefficient): The shape of the stability region depends on k, reaching its maximal size when k = 1, which corresponds to the Jordan normal form of the matrix, and minimal size when k = 0 or k = ∞, which corresponds to the Frobenius normal form of the matrix. The stability region for the 2-dimensional problem is the area in the complex plane where the absolute values of complex conjugate eigenvalues of the matrix R(Ah) are not greater than 1.
Let stability region be R h (Ah) for stepsize h , R 2h (Ah) for stepsize 2h, and so on. Then, for the two-dimensional test problem, the s-stage multistep semi-implicit method can be written in terms of a shift operator z : x n+s = z s x n as follows: From the recurrence equation (12), implying that x 0 (0 0) , and denoting the identity matrix as I, one can obtain the characteristic equation of the method: Then, the stability region of the ESIMM method is a region where the maximal absolute value of the root z max of the characteristic polynomial det Iz s − k 1 R h (Ah)z s−1 − . . . − k s R sh (Ah) satisfies the condition |z max | ≤ 1.
Stability regions of "full" and "short" ESIMM methods are almost the same and are given in Figures 10 and 11. For comparison, we provide stability regions of explicit Adams-Bashforth and implicit Adams-Moulton linear multistep methods as Figures 12 and 13 Let stability region be ( ℎ) for stepsize ℎ , ( ℎ) for stepsize 2ℎ, and so on. Then, for the two-dimensional test problem, the -stage multistep semi-implicit method can be written in terms of a shift operator : = as follows: From the recurrence equation (12), implying that (0 0) , and denoting the identity matrix as , one can obtain the characteristic equation of the method: Then, the stability region of the ESIMM method is a region where the maximal absolute value of the root of the characteristic polynomial det( − ( ℎ) − ⋯ − ( ℎ)) satisfies the condition | | 1.
Stability regions of "full" and "short" ESIMM methods are almost the same and are given in Figure 10 and Figure 11. For comparison, we provide stability regions of explicit Adams-Bashforth and implicit Adams-Moulton linear multistep methods as Figure 12 and Figure 13, respectively. Let stability region be ( ℎ) for stepsize ℎ , ( ℎ) for stepsize 2ℎ, and so on. Then, for the two-dimensional test problem, the -stage multistep semi-implicit method can be written in terms of a shift operator : = as follows: From the recurrence equation (12), implying that (0 0) , and denoting the identity matrix as , one can obtain the characteristic equation of the method: Then, the stability region of the ESIMM method is a region where the maximal absolute value of the root of the characteristic polynomial det( − ( ℎ) − ⋯ − ( ℎ)) satisfies the condition | | 1.
Stability regions of "full" and "short" ESIMM methods are almost the same and are given in Figure 10 and Figure 11. For comparison, we provide stability regions of explicit Adams-Bashforth and implicit Adams-Moulton linear multistep methods as Figure 12 and Figure 13, respectively.  Observing the results of stability analysis, one can see that theoretically predicted properties of the ESIMM methods and the results of the order plot analysis were successfully confirmed. ESIMM methods possess numerical stability that is comparable with implicit Adams-Moulton methods and is far superior to the stability of explicit Adams-Bashforth methods. Despite the fact that the numerical stability of semi-implicit methods strongly depends on the symmetry of the investigated system, we can declare that even in a worth case (full asymmetry of the simulated system) their stability is better than the stability of explicit linear multistep methods which allows implementing these methods with larger stepsizes.

Conclusions and Discussion
In this paper, we proposed novel multistep extrapolation methods based on semi-implicit symmetric integration. We explicitly show that the proposed ESIMM integration technique in many cases of nonlinear systems possesses greater performance then well-known multistep ODE solvers. We investigated the order and convergence of the investigated methods experimentally through the order plots using Rossler chaotic system as a testbench. One of the important research questions was the analytical estimation of the numerical stability of ESIMM methods. The investigation of the numerical stability of semi-implicit methods is rather complicated due to their non-trivial call to the right-hand-side function. We used a technique based on second-order linear test problems [11,20] which allowed plotting stability regions for ESIMM and comparing them with stability regions of linear multistep methods.
We explicitly show that the proposed ESIMM solvers are more numerically stable than explicit Adams-Bashforth methods, and, in some cases, implicit Adams-Moulton methods. Moreover, the stepsize range where the order of ESIMM methods preserves, theoretically allows constructing Observing the results of stability analysis, one can see that theoretically predicted properties of the ESIMM methods and the results of the order plot analysis were successfully confirmed. ESIMM methods possess numerical stability that is comparable with implicit Adams-Moulton methods and is far superior to the stability of explicit Adams-Bashforth methods. Despite the fact that the numerical stability of semi-implicit methods strongly depends on the symmetry of the investigated system, we can declare that even in a worth case (full asymmetry of the simulated system) their stability is better than the stability of explicit linear multistep methods which allows implementing these methods with larger stepsizes.

Conclusions and Discussion
In this paper, we proposed novel multistep extrapolation methods based on semi-implicit symmetric integration. We explicitly show that the proposed ESIMM integration technique in many cases of nonlinear systems possesses greater performance then well-known multistep ODE solvers. We investigated the order and convergence of the investigated methods experimentally through the order plots using Rossler chaotic system as a testbench. One of the important research questions was the analytical estimation of the numerical stability of ESIMM methods. The investigation of the numerical stability of semi-implicit methods is rather complicated due to their non-trivial call to the right-hand-side function. We used a technique based on second-order linear test problems [11,20] which allowed plotting stability regions for ESIMM and comparing them with stability regions of linear multistep methods.
We explicitly show that the proposed ESIMM solvers are more numerically stable than explicit Adams-Bashforth methods, and, in some cases, implicit Adams-Moulton methods. Moreover, the stepsize range where the order of ESIMM methods preserves, theoretically allows constructing Observing the results of stability analysis, one can see that theoretically predicted properties of the ESIMM methods and the results of the order plot analysis were successfully confirmed. ESIMM methods possess numerical stability that is comparable with implicit Adams-Moulton methods and is far superior to the stability of explicit Adams-Bashforth methods. Despite the fact that the numerical stability of semi-implicit methods strongly depends on the symmetry of the investigated system, we can declare that even in a worth case (full asymmetry of the simulated system) their stability is better than the stability of explicit linear multistep methods which allows implementing these methods with larger stepsizes.

Conclusions and Discussion
In this paper, we proposed novel multistep extrapolation methods based on semi-implicit symmetric integration. We explicitly show that the proposed ESIMM integration technique in many cases of nonlinear systems possesses greater performance then well-known multistep ODE solvers. We investigated the order and convergence of the investigated methods experimentally through the order plots using Rossler chaotic system as a testbench. One of the important research questions was the analytical estimation of the numerical stability of ESIMM methods. The investigation of the numerical stability of semi-implicit methods is rather complicated due to their non-trivial call to the right-hand-side function. We used a technique based on second-order linear test problems [11,20] which allowed plotting stability regions for ESIMM and comparing them with stability regions of linear multistep methods.
We explicitly show that the proposed ESIMM solvers are more numerically stable than explicit Adams-Bashforth methods, and, in some cases, implicit Adams-Moulton methods. Moreover, the stepsize range where the order of ESIMM methods preserves, theoretically allows constructing highly efficient multistep extrapolation solvers with adaptive stepsize. Being based on numerical extrapolation, ESIMM methods possess the same capacity to construct built-in estimators of local truncation error as single-step extrapolation solvers. Besides, the ESIMM ODE solvers do not need the concatenation procedure for output data obtained during the starting method and main solution, which also simplifies and fastens the algorithm. The computational cost of the proposed methods was estimated experimentally and is visualized through the performance plots. The program implementation of the proposed methods is significantly simpler than for implicit multistep solvers. Speaking generally, the ESIMM methods are almost as simple as explicit multistep methods, being almost as precise as implicit multistep methods. One of the main drawbacks of the ESIMM method is the necessity to increase the number of basic method calculations by 1 per each increase in the accuracy order. Therefore, high-order solvers, based on the proposed technique, can be less efficient.
In our further studies we plan to implement ESIMM solvers with adaptive stepsize, which requires the development of special step control and coefficient recalculation techniques. ESIMM methods based on other symmetric integrators, e.g., implicit midpoint method or composition solvers [10], are also to be studied.