Research Article  Open Access
A Numerical Algorithm for Solving Stiff Ordinary Differential Equations
Abstract
An advanced method using block backward differentiation formula (BBDF) is introduced with efficient strategy in choosing the step size and order of the method. Variable step and variable order block backward differentiation formula (VSVOBBDF) approach is applied throughout the numerical computation. The stability regions of the VSVOBBDF method 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.
1. 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 wellknown 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–6]. The problems considered in this paper are for the numerical solution of the initial value problem, with given initial values 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 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 (VSVOBBDF).
In Section 2, we introduce the general formulation of VSVOBBDF of order 3 to 5. Order conditions are listed in Section 3 while in Section 4, we consider the implementation of VSVOBBDF. The analysis of stability regions is illustrated in Section 5. This is followed by the strategy in choosing the step size and order of VSVOBBDF in Section 6. Numerical results are presented in Section 7. The appendix describes the algorithm applied for VSVOBBDF method.
2. VSVOBBDF Method Formulation
Two values of and 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 (, , and ) are distinguished by the number of back values contained in total blocks. The ratio distance between current and previous step 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 : where
Predictors for the first point and second point were computed using back values as the interpolating points. The resulting Lagrange polynomial for each order was given as follows.
For VSVOBBDF of order For VSVOBBDF of order For VSVOBBDF of order
Substituting and gives the predictor for the first and second points, respectively. Therefore by letting , and . This produced the following coefficients (Tables 1, 2, and 3) for the first and second points of predictor formulae for VSVOBBDF method.



The interpolating polynomial of the function using Lagrange polynomial in (2) gives the following corrector for the first point and second point . The resulting Lagrange polynomial for each order was given as follows.
For VSVOBBDF of order For VSVOBBDF of order For VSVOBBDF of order
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 and implicit otherwise.
Substituting and gives the corrector for the first and second points, respectively. Therefore by letting , and . This produced the following coefficients (Tables 4, 5, and 6), as in Definition 1 for the first and second points of VSVOBBDF method.



3. Order Conditions for General VSVOBBDF
As similar to analysis for order of Linear Multistep Method (LMM) given in [15], we use the following to determine the order of VSVOBBDF method.
Definition 2. The LMM [15] and the associated difference operator defined by are said to be of order if , . The general form for the constant is defined as Consequently, BBDF method can be represented in standard form by an equation , where and are by matrices with elements and for . Since VSVOBBDF 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 and .
4. Implementation of VSVOBBDF Method
Throughout this section, we illustrate the effect of Newtontype scheme which in a general form of The general form of VSVOBBDF method is with and are the back values. By setting, Equation (15) in matrixvector form is equivalent to Equation (17) is simplified as
Newton iteration is performed to the system ; by taking the analogous form of (14) where is the Jacobian matrix of with respect to . Equation (14) is separated to three different matrices denoted as
Twostage Newton iteration works to find the approximating solution to (1) with two simplified strategies based on evaluating the Jacobian () and LU factorization of [16]. Twostage Newton iteration is carried out as follows.
Step 1. Compute the values for , where
Step 2. Calculate the corrected value for with the value from (1).
Step 3. Solve for the second iteration.
Step 4. Let 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.
5. Stability Conditions for VSVOBBDF Method
To begin this section, we provide some definitions for the stability conditions of VSVOBBDF method as in (15).
Definition 3. A method is said to be zero stable if the roots of the polynomial satisfy the condition , .
Definition 4. A method is said to be absolutely stable in a region for a given if for that all the roots of the stability polynomial satisfy , .
Definition 5. A method is said to be stiffly stable if and are contained in the absolute stability region, and it is accurate for all when applied to the scalar test equation ; is a complex constant with , where , , and , , and are positive constants.
The stability polynomial, associated with the method of (15), is given by , while the absolute stability region of this method in the plane is determined by solving . Consequently, to determine zero stable, we substitute to for each order of VSVOBBDF 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 , . We obtained the stability region by finding the region for which . 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 , , and of VSVOBBDF method, respectively. The stability regions lie outside the closed region for each case.
(a)
(b)
(c)
Based on Figure 2, VSVOBBDF method possesses the region absolute stability, which contains almost the whole of the halfplane .
6. Choosing Order and Step Size
To increase the efficiency in BBDF algorithm, VSVOBBDF 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 VSVOBBDF 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 , , and are selected depending on the occurrence of every successful step. Consequently, the new step size 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 is the th order method and is the th order method. By finding the LTEs, the maximum step size is defined as where is the step size from previous block and 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 1. From order 3 to order 4 and from order 4 to order 3, the step size is allowed to increase, decrease, or be maintained, that is, , or .
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, , or .
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, or .
The successful step is dependent on the condition . If this condition fails, the values of , are rejected, and the current step is reiterated with step size selection . On the contrary, the step size increment for each successful step is defined as and if then ,
where is the safety factor, is the order of the method while and 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.
7. Numerical Results
We carry out numerical experiments arising from problems in physics to compare the performance of VSVOBBDF 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, where odefun is a function stored that evaluates the right side of the differential equations, int is a vector of the integration interval, , init is a column vector of initial conditions, , 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.
Problem 1. Consider, with solution
Problem 2. Consider,
with solution
See Kaps [18].
Problem 3. Consider,
with solution
See Lambert [15].
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, VSVOBBDF 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 VSVOBBDF gives the lowest maximum error for every tolerance level.
