Semi-Implicit and Semi-Explicit Adams-Bashforth-Moulton Methods

: Multistep integration methods are widespread in the simulation of high-dimensional dynamical systems due to their low computational costs. However, the stability of these methods decreases with the increase of the accuracy order, so there is a known room for improvement. One of the possible ways to increase stability is implicit integration, but it consequently leads to sufﬁcient growth in computational costs. Recently, the development of semi-implicit techniques achieved great success in the construction of highly efﬁcient single-step ordinary differential equations (ODE) solvers. Thus, the development of multistep semi-implicit integration methods is of interest. In this paper, we propose the simple solution to increase the numerical efﬁciency of Adams-Bashforth-Moulton predictor-corrector methods using semi-implicit integration. We present a general description of the proposed methods and explicitly show the superiority of ODE solvers based on semi-implicit predictor-corrector methods over their explicit and implicit counterparts. To validate this, performance plots are given for simulation of the van der Pol oscillator and the Rossler chaotic system with ﬁxed and variable stepsize. The obtained results can be applied in the development of advanced simulation


Introduction
Multistep integration is among the prospective approaches for solving ordinary differential equations (ODE). Due to a single call to the right-hand side function per integration step, such algorithms can be efficient in real-time and long-term simulations. However, the well-known problem of such methods is a decrease in numerical stability with an increase in the accuracy order of the applied scheme. This property restricts the application of high-order multistep methods for solving stiff ODE systems, excluding implicit methods, for example, backward differentiation formulas [1][2][3]. Another promising class of numerical integration methods is semi-explicit and semi-implicit integrators [4][5][6]. Initially developed among so-called symplectic methods [7] for solving Hamiltonian problems they were later extended on a broader class of dynamical systems [8][9][10][11][12][13]. Currently, many symplectic Runge-Kutta methods were developed [7] as well as the exactly symplectic methods of higher order [8]. These algorithms were initially designed to preserve the geometric properties of the simulated continuous system in the discrete model. However, later it was discovered that even in the case of non-Hamiltonian systems, extrapolation [14] and composition [15] ODE solvers based on symmetric and semi-implicit methods possess lower solution error and computational costs than their explicit and implicit counterparts [16][17][18][19]. The hypothesis of this study is that introducing semi-implicit integration can improve the performance of multistep ODE solvers because the semi-implicit Euler method possesses better precision and stability than an explicit Euler integrator. Therefore, it is of interest to develop new multistep methods that involve semi-implicit and semi-explicit integration techniques with increased stability and precision.

Semi-Explicit and Semi-Implicit Multistep Methods
In this study, we propose and investigate semi-explicit and semi-implicit modifications of the integration scheme called the Adams-Bashforth-Moulton formula (ABM) that is one of the well-known predictor-corrector multistep methods. These methods possess relatively good stability and convergence properties [20]. The proposed integration techniques exist for ODE systems of order two and higher, degenerating to traditional ABM method for a system of order one. Let us consider the following two-dimensional initial value problem (IVP) where f and g are known right-hand side functions. The system (1) is to be solved with the initial condition x(0) = x 0 and y(0) = y 0 . Let us denote two values of state variables obtained by the explicit Adams-Bashforth method as x p n+1 and y p n+1 . Then, for the two-dimensional system (1) the corrector method is as follows where h is the integration step, b i are coefficients of the implicit Adams-Moulton method and k is the number of stages.

Semi-Explicit Adams-Bashforth-Moulton Method
The proposed semi-explicit technique is based on the idea of the symplectic Euler method [7]. The semi-explicit technique supposes using of already calculated values of state variables to more accurately approximate remaining values at the current integration step. We apply this approach to multistep integration as follows. Replacing x p n+1 in (3) by x n+1 we get We will call formulae (4) the semi-explicit Adams-Bashforth-Moulton method. One can see, that if the one-dimensional system is under consideration, then the proposed technique reduces to the conventional Adams-Bashforth-Moulton method.

Semi-Implicit Adams-Bashforth-Moulton Method
In our previous studies [14,17,18], the generalized semi-implicit modification of the Euler method for non-Hamiltonian systems, called the D-method, was considered. Each equation of the finite-difference scheme obtained using the D-method contains one implicit calculation of the variable that corresponds to the equation being solved. Thus, if we write a matrix of size N × N (N is the number of state variables) and mark the implicit calculation of 1, and 0 is explicit, then in the D-method ones will fill the main diagonal of the matrix. Moreover, each obtained value of the state variable is used in the same step to calculate the remaining values as in the case of the semi-explicit Euler method. Using this approach in the multistep interpretation, for system (1) we obtain the semi-implicit Adams-Bashforth-Moulton scheme where x n+1 in (5a) and y n+1 in (5b) can be calculated using the Newton method or even the simple iterations method, which converges to the implicit solution in 2-4 iterations [16]. In simple cases the analytical solution of algebraic equations with respect to the values of state variables at (n + 1) point is possible. This simple method is applicable to the chosen test systems and we used it in our study. Both considered algorithms were applied to the two-dimensional system for easy understanding. However, the proposed approaches can be used to solve ODE systems of arbitrary order. Moreover, changing the order of lines in the calculation algorithm for both considered methods, one can obtain a family of semi-implicit and semi-implicit multistep algorithms with similar properties. Such methods can be used to create new tools for dynamical systems analysis, numerical stability control, local truncation error estimation and more accurate simulation of orbits using the interval approach [21].
It is expected that, as in the case of single-step ODE solvers [14,17,18], the introduction of semi-explicit and semi-implicit calculations will increase the accuracy of the obtained numerical solution and the computational costs will not change in comparison with the original Adams-Bashforth-Moulton method. Let us consider some results of experimental studies of the proposed techniques in comparison with well-known multistep ODE solvers with fixed and adaptive stepsize.

Experimental Results
In the experimental part of our study, we compared semi-explicit (SEABM) and semi-implicit (SIABM) techniques with the Adams-Bashforth (AB), Adams-Moulton (AM), Adams-Bashforth-Moulton (ABM) methods and the Backward Differentiation Formula (BDF). We examined four-stage methods with the fixed integration step and fifth-stage methods with variable stepsize. Recurrence formulas were applied for the step control and the method for estimating the local error described in Reference [1]. As a reference, we used the numerical solution obtained by the Dormand-Prince method of order 8. We studied two sample ODE systems.

Van der Pol System
In our first experiment, the well-known benchmark nonlinear model-the van der Pol oscillator-was considered. In the IVP form its equations are as follows where µ is the stiffness parameter. We simulated system (6) with µ = 1. We estimated time costs needed for obtaining numerical solutions by ODE solvers under investigation as well as the maximal global truncation error. We simulated system (6) T = 50 s from the initial conditions x 0 = 0.1, y 0 = 0 using fixed and variable integration step. In experiments with the fixed step, we simulated system (6) by the semi-explicit method using the following finite-difference scheme where f and g are the values of the right-hand side functions of system (6) at a certain point in time.
For simulation using semi-implicit integration, we applied the same predictor scheme and the corrector algorithm was The obtained performance plots for system (6) are shown in Figures 1 and 2. One can see, that in the case of fixed step integration (Figure 1) solvers based on the ABM method and its SE and SI modifications are the most efficient. However, the proposed SEABM method demonstrates the best results.  Figure 2 represents the results of performance analysis for studied solvers with the variable stepsize. It is known that for multistep integration methods step-size control supposes more computational costs than for single-step solvers. In our study, we implemented methods with an uneven grid, for which the recursive calculation of the method coefficients is required when the integration step value is changed [1]. As one can see, this leads to a reduction of the difference in the computational costs between explicit and implicit algorithms (Figure 2). The proposed SIABM and SEABM solvers possess a higher speed of computation with similar accuracy being compared with the known ABM technique.

Rossler Chaotic System
As the second sample system, the Rossler oscillator [22] was chosen. This system is described by the following ODEs where a, b, c are free parameters. We investigated system (7) with a = 0.2, b = 0.2, c = 0.7, which correspond to the case of chaotic oscillations [22]. System (7) was simulated with variable and fixed integration step for time T = 50 s from initial conditions x 0 = 0.1, y 0 = 0, z 0 = −0.1 by all compared methods. Figures 3 and 4 represent the results of this experiment. As one can see, in the case of the three-dimensional system, superiority in the performance of fixed-step solvers based on the proposed methods over the conventional ABM scheme is more noticeable. The accuracy of the solutions calculated by semi-explicit and semi-implicit techniques is comparable with the results obtained for the implicit AM method. All of the predictor-corrector methods appear to be the most efficient among the considered integration schemes.  Step-size control reinforces the difference in the performance of ODE solvers based on predictor-corrector schemes in comparison with the explicit AB method (Figure 4). Moreover, the explicit AB algorithm shows a significant decrease in accuracy when the system with chaotic behavior is simulated. The most effective method is the semi-explicit technique, as was in the case of the van der Pol oscillator.  Thus, we can conclude that semi-implicit and semi-implicit multistep integration for both considered systems are more efficient when implemented with the variable step.

Discussion
The computational costs of the proposed semi-implicit solver can increase significantly if Newton's iterations will be required for approximation of diagonally implicit variable. Figure 5 shows the simulation results for the Rossler system using proposed four-stage ODE solvers with various approximation techniques for diagonal implicitness. We compared the SEAM algorithm with methods which involve Newton's iterations (SIABMN) and simple iterations (SIABMSI) as well as the SIABM method with the analytical solution of diagonally implicit algebraic equations. To clarify the comparison and increase the understandability we provide listings of the used algorithms in Appendix A, excluding code for the (SIABMN) method due to its massiveness. As one can see from Figure 5, the semi-implicit ABM modification with the Newton method is the most computationally expensive solver. This method should be applied only if it is impossible to resolve a system of diagonally implicit algebraic equations. The examples of such systems can be the simple model of spiking neurons [23] and the Hodgkin-Huxley system [24] which are broadly used in computational neuroscience. An alternative to Newton's method is the simple iterations algorithm (SIABMSI in Figure 5), which is less computationally expensive and is proven to be convergent for arbitrary right-hand side function [16].  One can note that the semi-explicit integration is free from this limitation but potentially possesses less stability and convergence, especially for high-dimensional systems.

Conclusions
In this study we proposed the simple technique to enhance the performance of predictor-corrector multistep ODE solvers. We described two versions of the modified Adams-Bashforth-Moulton method with semi-implicit and semi-explicit integration. The proposed numerical integration methods possess high precision and low computational costs. In the experimental part of the study, we have explicitly shown that the performance of semi-implicit and semi-explicit multistep integration methods is superior to their explicit and implicit counterparts. To illustrate this, we simulated two nonlinear systems: the van der Pol oscillator and the chaotic Rossler system. In both cases, the computational efficiency of semi-explicit and semi-implicit Adams-Bashforth-Moulton methods was highest among the tested algorithms. Implementation with adaptive stepsize confirmed this property of proposed methods. In our further research, we will compare these algorithms with high-order single-step solvers, including semi-explicit and semi-implicit extrapolation and composition methods, as well as the effects arising in the long-term simulation of nonlinear systems. #Predictor Listing A2. Single integration step for the semi-implicit Adams-Bashforth-Moulton method. In this case the analytic solution of the diagonally implicit algebraic equations is derived.