A new one-step method with three intermediate points in a variable step-size mode for stiff differential systems

This work introduces a new one-step method with three intermediate points for solving stiff differential systems. These types of problems appear in different disciplines and, in particular, in problems derived from chemical reactions. In fact, the term “stiff”’ was coined by Curtiss and Hirschfelder in an article on problems of chemical kinetics (Hirschfelder, Proc Natl Acad Sci USA 38:235–243, 1952). The techniques of interpolation and collocation are used in the construction of the scheme. We consider a suitable polynomial to approximate the theoretical solution of the problem under consideration. The basic properties of the new scheme are analyzed. An embedded strategy is adopted to formulate the proposed scheme in a variable stepsize mode to get better performance. Finally, some models of initial-value problems, including ordinary and time-dependent partial differential equations, are solved numerically to assess the performance and efficiency of the proposed technique, with applications to real-world problems.


Introduction
In this article, we are concerned with constructing a strategy to obtain an approximate solution to the following first-order IVP: (1) where t 0 and t N denote the initial and final points of the integration interval, and the function w(t, u) fulfils the Lipchitz's condition that guarantees the existence and uniqueness theorem given in [2].
Many physical phenomena in which a quantity varies in time and at the same time relies upon its past values could be modelled by (1). For example, equation (1) has been used to model many real-world application problems in different disciplines, including physics, electrical and mechanical engineering, economics, mathematical chemistry, or mathematical biology.
Numerical treatments of a type of problem (1) have received a lot of consideration because most of the differential equations arising from the modelling of physical phenomena do not have known theoretical solutions. Thus, it is of practical interest to consider numerical solutions for such model problems.
We remark that to choose the most efficient and stable numerical technique for computing (1), it is preferable to have qualitative information about the solution. The primary information is the stiffness of the problem. The explicit Runge-Kutta methods (ERKMs) are generally unsuitable for the solution of stiff problems because ERKMs may fail to provide the solution or give wrong results. The instability of ERKM methods motivated researchers in this field to consider implicit Runge-Kutta methods (IRKMs). Several IRKMs have been proposed and applied by different scholars to solve stiff IVPs of the form in (1). Among the various IRKMs you can consult [3][4][5][6][7][8][9][10].
Another numerical strategy for integrating (1) is the implicit multistep-methods (IMMs); IMMs utilize the estimation of the approximation at more than one mesh point to calculate approximate solutions at the next point. The IMMs need another strategy to start them, leading to an increase in the total number of function evaluations and CPU time. Numerous outstanding researchers in Numerical Analysis have used different IMMs for solving (1). Some interesting monographs that have employed IMMs to provide approximate solutions to stiff and non-stiff problems are those in [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28].
Several approximate techniques have been proposed recently for the numerical treatment of (1). Among those techniques, we can mention a higher-order embedded method proposed by Kim et al. [29], the Chebyshev collocation technique by Piao et al. [30], the highly stable implicit-explicit schemes presented by Izzo and Jackiewicz [31], the general linear methods developed by Abdi and Jackiewicz [32], the new hybrid symmetric technique by Kovalnogov et al. [33], a fitted finite difference approach by Chen and Simos [34], or the Numerov-type methods proposed by Medvedeva et al. [35].
This article presents and applies a suitable error estimation and stepsize controller for a newly developed one-step A-stable method (OSASM) for solving (1). Although the development of the method and its theoretical analysis has been carried out for a scalar problem, we must point out that it can also be applied for solving initial-value differential systems, as shown in Sect. 5.

Mathematical formulation of the OSASM
Consider a partition t 0 < t 1 < · · · < t N with h = t j+1 − t j and suppose that the theoretical solution of the IVP in (1) on a generic subinterval [t n , t n+1 ], is approximated by the following polynomial whose first derivative can be approximated by where a j ∈ R are unknown coefficients that would be obtained imposing some collocation and interpolation conditions at specified nodes. The OSASM is formulated by considering three off-step points t n+c 1 = t n + c 1 h, t n+c 2 = t n + c 2 h, and t n+c 3 Equation (2) is evaluated at t n , and the one in (3) is also evaluated at t n , t n+c 1 , t n+c 2 , t n+c 3 , t n+1 to get a system of equations with six unknowns a n , n = 0(1)5, where u n+ j and w n+ j denote respectively approximations of u(t n+ j ) and u (t n+ j ) = w(t n+ j , u(t n+ j )).
After getting the values of a n , n = 0(1)5, and substituting them into the Eq. (2), and replacing the variable t = t n + xh, then, the polynomial in (2) may be written as where α 0 (x) = 1, and the coefficients {β i (x)} i=0,c 1 ,c 2 ,c 3 ,1 depend on c 1 , c 2 , c 3 . Taking x = 1 in (4), which corresponds to the point t n+1 = x n + h, after inserting , we obtain the following formula, which is the main scheme of the proposed OSASM: We note that the values of the c i are not taken arbitrarily, but after optimizing the local truncation error (LTE) of the formula in (4). In order to provide numerical solutions to problem (1) simultaneously at the points t n+c 1 , t n+c 2 , t n+c 3 , t n+1 , we also take x = c 1 , x = c 2 , x = c 3 . Doing this, we obtain a total of four formulas that form the OSASM. The remaining formulas are

Theoretical analysis of the OSASM technique
This section gives the theoretical analysis of the OSASM technique, concerning the order of accuracy, convergence and linear stability analysis.

Order of accuracy, consistency, and convergence of the OSASM
We rewrite the scheme in (5)-(6) as whereR,S denote constant matrices containing the coefficients of the Formulas (5)-(6), and We define the operator γ related to the OSASM in (5)-(6) as where j and j are respectively vector columns ofR andS, and I = {0, c 1 , c 2 , c 3 , 1}. Expanding u (t n + jh) and u (t n + jh) in Taylor series about t n , we get and q = 0, 1, 2, . . . Considering the Formulas (5)-(6), we get θ 0 = 0, θ 1 = 0, . . . , θ 5 = 0 and which is the vector of LTEs of the OSASM formulas. The method in (5)-(6) may be reformulated as a Runge-Kutta method, and thus the theory for these methods can be applied here. The vector of stage values is given by c = (0, c 1 , c 2 , c 3 , 1) and the coefficients of the Butcher array can be obtained readily from the Formulas (5)-(6). It is a routine task to check that the order conditions are verified up to order 6, which ultimately is the order of the method. It is worth remarking that the scheme is zero-stable since the OSASM is a onestep method. Since the order is greater than one, then OSASM is consistent with the differential problems under consideration. It is also important to note that the OSASM is convergent because OSASM is zero-stable and consistent.

Linear stability analysis of the OSASM
The concept of linear stability analysis (LSA) plays an essential role in the solution of IVPs in ODEs. This LSA is frequently referred to as absolute stability, and is crucial when solving stiff differential equations. Owing to this, we study the linear stability of the newly developed technique by applying the method in (7) to the following standard test equation where η stands for a complex parameter and the above equation has a bounded solution that will be zero as t tends to infinity. Applying the scheme in (5)-(6) to (11) to have where η is a complex parameter and S(z) denotes the following stability matrix Fig. 1 Absolute stability region in the complex z-plane for the OSASM using (11) with The obtained eigenvalues for S(z) are 0, 0, 0, 5z 4 +198z 3 +2076z 2 +10,080z+20,160 5z 4 −198z 3 +2076z 2 −10,080z+20,160 . It is observed that the numerical solution relies upon the reported eigenvalues, and the stability properties of the new scheme will be characterized by the spectral radius, ρ(S(z)). The region of absolute stability S for the proposed method is represented as S = {z ∈ C : |ρ(S(z))| < 1} and its domain is displayed in Fig. 1. The proposed OSASM is A-stable because the left half complex plane is contained in S.

Error estimation and control strategy of the proposed method
In order to obtain an efficient implementation of any numerical technique and solve challenging problems like the stiff differential equations, estimating and controlling some measure of errors by varying the step size is required. The proposed method given in (5)-(6) may be formulated and computed in variable step-size mode (VSS) by using a lower order technique to estimate the local error (LE) at the final point on each iteration. To achieve a robust estimation of the LTE, we adopt a similar approach as in the book by Ascher and Petzold [27]. The following multi-step formula of order p = 4 and local truncation error (LT E) = − 17h 5 u (5) (t n ) 120,960 + O h 6 has been used to estimate the LE at the end-point. This error estimate provides the basis for determining the step-size for the next step. The estimation of the LE is obtained through the quotient where we have taken · as the maximum norm, and ATOL and RTOL are predefined absolute and relative tolerances, respectively. If E ST ≤ ATOL, then we accept the results and take the next step as h new = 2 × h old . If E ST > ATOL, then we reject the obtained results by decreasing the step-size and repeat the calculations with the following new step where 0 < μ < 1 is a suitable adjustment factor whose purpose is to avoid failed steps during implementation.

Implementation details of the OSASM
The proposed A-stable technique is implemented in a step by step mode. We rewrite the system in (5)-(6) as W (u) = 0 where the unknowns arẽ Since the OSASM is an implicit scheme, we use modified Newton's method (MNM) to solve the obtained system of non-linear equations. The MNM is given bỹ where J represents the Jacobian matrix of W. The starting values to use MNM for solving the systems on each iteration are taken as u n+ j = u n + ( jh)w n , n = 0, 1, 2, . . . , N − 1, j = c 1 , c 2 , c 3 , 1.
The developed method can also be applied using a componentwise implementation to solve systems of first-order IVPs with m equations of the form:

Numerical experiments
Here, we illustrate the numerical performance of the proposed OSASM on some real-life model first-order IVPs. We note that the effectiveness and efficiency of the method in (5)-(6) and the strategy in (14) will rely on the VSS approach given in (15). To measure the accuracy of the proposed OSASM with other existing numerical strategies, the following error formulas have been considered where u(t j ) is the exact solution, and u j is the computed result at each point t j of the discrete grid. The following methods have been considered for comparisons: • OSASM: The one-step A-stable method presented in (5)-(6).
• CCM: The Chebyshev Collocation Method of order eight in [30].
• NDSolve : The built-in package solver in Mathematica.

Test Problem 1
In the first problem, we solve a scalar stiff problem [36], given by The true solution of this problem is u(t) = cos(2π t).
In Table 1, we can observe the data for this problem. Not only the proposed method is the one that provides the best maximum absolute error (MAXAE), but it presents the lowest number of steps (NS), being the most efficient of the methods considered. Figure 2 shows the absolute errors along the integration interval and the exact and the discrete solutions for the proposed method for ATOL = RTOL = 10 −6 , h ini = 10 −2 .

Test Problem 2
Consider a non-linear stiff first-order system of ODEs [37] given by u 1 (t) = −10,004u 1 (t) + 10,000u 2 (t) 4 , u 1 (0) = 1, The exact solution to this system is The results in Table 2 show again the good performance of the proposed method.

Test Problem 3
We also consider the non-linear well-known Robertson's chemical problem which has been studied earlier by Mazzia and Iavernaro [38] To compare the maximum absolute errors, the following reference solution at the final point t N = 40 given in [38] has been used, The data in Table 3 indicates the outstanding performance of the proposed method for this problem, where the total number of function evaluations (NFEVAL) is considerably reduced compared with the methods used for comparison.

Test Problem 4
We consider the Chapman model problem that describes how ozone is generated and destroyed with conditions useful for the lower to the middle stratosphere. In this case, where k 1 , k 2 , k 3 , k 4 stand for the reaction rates, M stands for an additional molecule required to carry off excess energy, and hv represents the absorption of light, a photochemical process. The first reaction produces ozone, but the second and third reactions produce ozone destruction. The amount of molecular oxygen, O 2 , is kept constant for this problem (a reasonable assumption in the atmosphere). The first two reaction rates are constant (only functions of temperature), whereas the subsequent two, which are determined by light absorption, change during the day. The resulting system has the following form if we set u 1 , u 2 , u 3 to be the concentrations of O, O 2 , O 3 , respectively:   Figure 3 shows the numerical results from the OSASM with ATOL = RTOL = 10 −6 , h ini = 10 −3 . This figure shows that the solutions provided by the proposed method are in good agreement with the numerical results in [40].

Test Problem 5
We also consider the Belousov-Zhabotinskii chemical reaction Oregonator problem originated from chemistry, given as follows: [37] The numerical results in Table 4 show again the good efficiency of the OSASM method.

Test Problem 6
Consider the following one-dimensional stiff Burgers' equation model problem given in [30] whose exact solution is unknown. Burgers' equation is a mathematical model of the motion of a viscous compressible gas, where u is the speed of the gas, v is the kinematic viscosity, x is the spatial coordinate, and t is the time. Test Problem 4 is solved by discretizing the first-order and second-order spatial derivatives, and then applying the proposed method, following a similar procedure as in the method of lines. The u x and u x x in Problem 4 are discretized by the classical where δx = (x M+1 − x 0 ) M + 1 , and M is the internal number of spatial nodes. The data in Table 5 indicates the good performance of the proposed method for the Burgers' problem. The approximate solutions provided by NDSolve and the OSASM method for Problem 6 are displayed in Fig. 4.

Conclusion
This article introduced a VSS formulation of a new OSASM. The proposed OSASM is sixth-order convergent with the favorable characteristic of A-stability. The proposed OSASM has been successfully used to integrate problem (1), and the obtained numerical results confirmed that the OSASM is a suitable solver for solving stiff differential systems. The procedure of controlling the step size using a reasonable error estimation and an adaptive strategy to get better accuracy is also described. The efficiency of the OSASM, when compared with some existing techniques, is remarkable.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.