A Numerical Algorithm for Solving Stiff Ordinary Differential Equations

An advancedmethod using block backward differentiation formula (BBDF) is introducedwith efficient strategy in choosing the step size and order of the method. Variable step and variable order block backward differentiation formula (VSVO-BBDF) approach is applied throughout the numerical computation.The stability regions of the VSVO-BBDFmethod are investigated and presented in distinct graphs.The improved performances in terms of accuracy and computation time are presented in the numerical results with different sets of test problems. Comparisons are made between the proposed method and MATLAB’s suite of ordinary differential equations (ODEs) solvers, namely, ode15s and ode23s.


Introduction
Many studies on solving the equations of stiff ordinary differential equations (ODEs) have been done by researchers or mathematicians specifically.With the numbers of numerical methods that currently exist in the literature, extensive research has been done to unveil the comparison between their rate of convergence, number of computations, accuracy, and capability to solve certain type of test problems [1].The well-known numerical methods that are used widely are from the class of BDFs or commonly understood as Gear's Method [2].However, many other methods that have evolved to this date are for solving stiff ODEs which arise in many fields of the applied sciences [3][4][5][6].The problems considered in this paper are for the numerical solution of the initial value problem, with given initial values () =  0 in the given interval  ∈ [, ].
The existing numerical methods have often been compared to one another to find the best approximation and the best method.MATLAB ODE suite is one of the most preferred solvers to be used for comparison purposes [7,8].Stiff ODE solvers that are available in MATLAB ODE suite are ode15s and ode23s which are based on the numerical differentiation formulas [8,9] and the modified Rosenbrock formula of order 2, respectively [8].
An improved method has been identified in [10], which was expanded from the method incorporated with BDF proposed by Gear.Since then, the study on producing block approximations  +1 ,  +2 ,  +3 , . . .,  + also known as Block Backward Differentiation Formulae (BBDF) [11] has attracted much attention.BBDF method by using variable step size in [12] demonstrates the competency of computing concurrent solution values at different points.Consequently, study in [13,14] is an extension of a previous study in a way that the accuracy is improved by increasing the order of the method up to order 5. Thus, these studies lead to enhancing the existing method to become Variable Order Variable Step Block Backward Differentiation Formula of order 3 until 5 (VSVO-BBDF).
In Section 2, we introduce the general formulation of VSVO-BBDF of order 3 to 5. Order conditions are listed in Section 3 while in Section 4, we consider the implementation of VSVO-BBDF.The analysis of stability regions is illustrated qh rh h P 5 Figure 1: VSVO-BBDF method of order (3-5).
in Section 5.This is followed by the strategy in choosing the step size and order of VSVO-BBDF in Section 6. Numerical results are presented in Section 7. The appendix describes the algorithm applied for VSVO-BBDF method.

VSVO-BBDF Method Formulation
Two values of  +1 and  +2 were computed simultaneously in block by using earlier blocks with each block containing a maximum of two points (Figure 1).The orders of the method (3, 4, and 5) are distinguished by the number of back values contained in total blocks.The ratio distance between current (  ) and previous step ( −1 ) is represented as  and  in Figure 1.In this paper, the step size is given selection to decrease to half of the previous steps or increase up to a factor of 1.9.For simplicity,  is assigned as 1, 2, and 10/19 for the case of constant, halving and increasing the step size, respectively.The zero stability is achieved for each of these cases and explained in the next section.We find approximating polynomials   (), by means of a -degree polynomial interpolating the values of  at given points, that are ( −3 ,  −3 ), ( −2 ,  −2 ), ( −1 ,  −1 ), . . ., ( +2 ,  +2 ): where for each  = 0, 1, . . ., . ( Predictors for the first point   +1 and second point   +2 were computed using back values as the interpolating points.The resulting Lagrange polynomial for each order was given as follows. For VSVO-BBDF of order 3 ( = 3) For VSVO-BBDF of order 4 ( = 4) ( For VSVO-BBDF of order 5 ( = 5) Substituting  = 0 and  = 1 gives the predictor for the first and second points, respectively.Therefore by letting  = 1,  = 1,  = 2,  = 2 and  = 1,  = 10/19.This produced the following coefficients (Tables 1, 2, and 3) for the first and second points of predictor formulae for VSVO-BBDF method.
The interpolating polynomial of the function () using Lagrange polynomial in (2) gives the following corrector for the first point   +1 and second point   +2 .The resulting Lagrange polynomial for each order was given as follows.
For VSVO-BBDF of order 3 ( = 3)  () =  ( +1 + ℎ) For VSVO-BBDF of order 4 ( = 4) For VSVO-BBDF of order 5 ( = 5) Linear Multistep Method (LMM) given in [15] is given in the definition below.Definition 1.The linear -step method can be represented in standard form by an equation where  + ≈ ( + ) and  + ≡ ( + ,  + ), coefficients   ,   are suitably chosen constants subject to conditions , and  is defined as the order of the particular method employed.This method is said to be explicit if   = 0 and implicit otherwise.

Order Conditions for General VSVO-BBDF
As similar to analysis for order of Linear Multistep Method (LMM) given in [15], we use the following to determine the order of VSVO-BBDF method.
Definition 2. The LMM [15] and the associated difference operator  defined by are said to be of order  if   =  1 = ⋅ ⋅ ⋅ =   = 0,  +1 ̸ = 0.The general form for the constant   is defined as Consequently, BBDF method can be represented in standard form by an equation ∑  =0    + = ℎ ∑  =0    + , where   and   are  by  matrices with elements  , and  , for ,  = 1, 2, . . ., .Since VSVO-BBDF for variable order () is a block method, we extend the Definition 2 in the form of And the general form for the constant   is defined as is equal to the coefficients of   , where  =  = ( − 2), . . .,  + 1,  + 2 and  = 3, 4, 5.

Implementation of VSVO-BBDF Method
Throughout this section, we illustrate the effect of Newtontype scheme which in a general form of The general form of VSVO-BBDF method is with  1 and  2 are the back values.By setting, Equation ( 15) in matrix-vector form is equivalent to Equation ( 17) is simplified as Newton iteration is performed to the system f+1,+2 = 0; by taking the analogous form of (14) where  +1,+2 = (/)( ()  +1,+2 ) is the Jacobian matrix of  with respect to .Equation ( 14) is separated to three different matrices denoted as Two-stage Newton iteration works to find the approximating solution to (1) with two simplified strategies based on evaluating the Jacobian ( +1,+2 ) and LU factorization of Â [16].Two-stage Newton iteration is carried out as follows.
Step 4. Let  (+1) +1,+2 be equal to the second iteration of The error is defined as and the maximum error is defined as where TS is the total steps and  is the number of equations.Consequently, the rest of the method is carried out as the appendix.The stability polynomial, (, ĥ) associated with the method of (15), is given by det( 2 −−), while the absolute stability region of this method in the ℎ plane is determined by solving det( 2 −  −  = 0).Consequently, to determine zero stable, we substitute ĥ = ℎ = 0 to (, ĥ) for each order of VSVO-BBDF method.

Stability Conditions for VSVO-BBDF Method
Hence, the roots for the three different step size and order () selections obtained by using Maple are listed in Table 7.
The stability region was given by the set of points determined by the boundary  =   , 0 ≤  ≤ 2.We obtained the stability region by finding the region for which || < 1.Since all of the roots in Table 7 have modulus less than or equal to 1, the method ( 15) is said to be zero stable.Figure 2 shows the stability for orders 3, 4, and 5 of VSVO-BBDF method, respectively.The stability regions lie outside the closed region for each case.
Based on Figure 2, VSVO-BBDF method possesses the region absolute stability, which contains almost the whole of the half-plane Re(ℎ) < 0.

Choosing Order and Step Size
To increase the efficiency in BBDF algorithm, VSVO-BBDF algorithm is designed to have the capacity to vary not only the step size, but also the order of the method employed.The importance of choosing the step size is to achieve reduction in computation time and number of iterations.
Meanwhile changing the order of the method is designed for finding the best approximation.The essential components of VSVO-BBDF algorithm are the local truncation error (LTE), strategies for deciding when to change step size and order, and techniques for changing the step size and order.Strategies proposed in [17] are applied in this study for choosing the step size and order.The strategy is to estimate the maximum step size for the following step.Methods of order  − 1, , and  + 1 are selected depending on the occurrence of every successful step.Consequently, the new step size ℎ new is obtained from which order produces the maximum step size.
The user initially will have to provide an error tolerance limit, TOL, on any given step and obtain the LTE for each iteration.The LTE is obtained from where  (+1) +2 is the ( + 1)th order method and  () +2 is the th order method.By finding the LTEs, the maximum step size is defined as where ℎ old is the step size from previous block and ℎ max is obtained from the maximum step size given in the above equations.
There are 3 cases of the possibilities in choosing the step size.
Case 2. From order 3 to order 5 and from order 4 to order 5, the step size is allowed to decrease, or be maintained, that is, ( = 1) ( = 1), ( = 2) ( = 2), or ( = 1) ( = 10/19).Case 3. From order 5 to order 3 and from order 5 to order 4, the step size is allowed to decrease, or be maintained, that is, The successful step is dependent on the condition LTE < TOL.If this condition fails, the values of  +1 ,  +2 are rejected, and the current step is reiterated with step size selection ( = 2).On the contrary, the step size increment for each successful step is defined as ℎ new =  × ℎ max and if ℎ new > 1.9 × ℎ old then ℎ new = 1.9 × ℎ old , where  is the safety factor,  is the order of the method while ℎ old and ℎ new is the step size from previous and current block, respectively.In this paper,  is set to be 0.8 so as to make sure the rejected step is being reduced.

Numerical Results
We carry out numerical experiments arising from problems in physics to compare the performance of VSVO-BBDF method with stiff ODE solvers in MATLAB version 7 in Windows XP or later.The syntax for using ode15s and ode23s is similar for every test problem.For example, [, ] = ode15s ("odefun", int, init, options) , where odefun is a function stored that evaluates the right side of the differential equations, int is a vector of the integration interval, [  ,  end ], init is a column vector of initial conditions, [ 0 ], and options is an adjustable format of optional parameters to change the default integration properties.We created the above option structure to vary the relative tolerance (RelTol) and absolute tolerance (AbsTol) according to specific error tolerance.These test problems are performed under different conditions of error tolerances-(a) 10 −2 , (b) 10 −4 , and (c) 10 −6 .The test problems and solutions are listed below.
This paper considers the comparison of four different factors, namely, number of steps taken, average error, maximum error, and computation time.From Table 8, among the three methods tested, our method, VSVO-BBDF method, requires the shortest execution time, smallest maximum error, and average error for each given tolerance level.From Figure 3, we can see more clearly that VSVO-BBDF gives the lowest maximum error for every tolerance level.
Again by comparing the four factors mentioned earlier, we can see VSVO-BBDF in Table 9 gives the minimum maximum error for every tolerance level except for TOL 10− 4.However, our method prevails in terms of average error for each given tolerance level.VSVO-BBDF once again requires the shortest execution time for each given tolerance level.Figure 4 illustrates the efficiency of VSVO-BBDF as compared to ode15s and ode23s.
From Table 10, VSVO-BBDF method gives the shortest execution time for every given tolerance level.While in terms of average error and maximum error, VSVO-BBDF method once again gives the best result as compared to ode15s and ode23s.Figure 5 clarifies the efficiency of the proposed method based on its total steps and tolerance level.

Conclusion
The objective is met when VSVO-BBDF method outperformed ode15s and ode23s in terms of average error as well as maximum error.In most of the cases, VSVO-BBDF has successfully managed to reduce the number of total steps taken.As for the computation time wise, it gave lesser values  for all cases.Therefore, we can conclude that VSVO-BBDF can serve as an alternative solver for solving stiff ordinary differential equations which arise in engineering and applied sciences.

Appendix
The following notation is used in the algorithm of VSVO-BBDF: ℎ: step size, : tolerance, TS: total steps, Order: order of the method, JCBN: jacobian, E: increment.
The algorithm for VSVO-BBDF code is given as follows.

Table 8 :
Numerical results for Problem 1.

Table 9 :
Numerical results for Problem 2.

Table 10 :
Numerical results for Problem 3.