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 (VSVO-BBDF) approach is applied throughout the numerical computation. The stability regions of the VSVO-BBDF 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 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–6]. The problems considered in this paper are for the numerical solution of the initial value problem, 𝑦^′=𝑓𝑥,𝑦 (1) 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 𝑦_𝑛+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 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.

2. 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.

Figure 1 

VSVO-BBDF method of order (𝑃3–𝑃5).

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): 𝑃_𝑘=^𝑘_𝑗=0𝑦𝑥_{𝑛+1−𝑗}⋅𝐿_𝑘,𝑗(𝑥), (2) where 𝐿_𝑘,𝑗(𝑥)=^𝑘_{𝑖=0𝑖≠𝑗}𝑥−𝑥_{𝑛+1−𝑖}𝑥_{𝑛+1−𝑗}−𝑥_{𝑛+1−𝑖}foreach𝑗=0,1,…,𝑘. (3)

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)𝑃(𝑥)=𝑃𝑥_𝑛+1=()+𝑠ℎ𝑟+1+𝑠𝑟𝑦_𝑛+(1+𝑠)𝑦−𝑟^𝑛−1. (4) For VSVO-BBDF of order 𝑃4(𝑃=4)𝑃(𝑥)=𝑃𝑥_𝑛+1=()+𝑠ℎ2𝑟+1+𝑠)(𝑟+1+𝑠2𝑟²𝑦_𝑛+(2𝑟+1+𝑠)(1+𝑠)−𝑟²𝑦_𝑛−1+(𝑟+1+𝑠)(1+𝑠)2𝑟²𝑦_𝑛−2. (5) For VSVO-BBDF of order 𝑃5(𝑃=5)𝑃(𝑥)=𝑃𝑥_𝑛+1=)+𝑠ℎ𝑞+2𝑟+1+𝑠(2𝑟+1+𝑠)(𝑟+1+𝑠2𝑞+2𝑟𝑟²𝑦_𝑛+𝑞+2𝑟+1+𝑠(2𝑟+1+𝑠)(1+𝑠)−𝑞+𝑟𝑟²𝑦_𝑛−1+𝑞+2𝑟+1+𝑠(𝑟+1+𝑠)(1+𝑠)2𝑞𝑟²𝑦_𝑛−2+(2𝑟+1+𝑠)(𝑟+1+𝑠)(1+𝑠)𝑦−𝑞−𝑞−𝑟−𝑞−2𝑟^𝑛−3. (6)

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.

Table 1 

Coefficients for the first and second point of predictor formulae for VSVO-BBDF when 𝑟=1, and 𝑟=1,𝑞=1.

Table 2 

Coefficients for the first and second point of predictor formulae for VSVO-BBDF when 𝑟=2, and 𝑟=2,𝑞=2.

Table 3 

Coefficients for the first and second point of predictor formulae for VSVO-BBDF when 𝑟=10/19, and 𝑟=1,𝑞=10/19.

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=(𝑠)+𝑠ℎ𝑟+1+𝑠)(𝑠+1)(𝑦2𝑟+4^𝑛+2+(𝑟+1+𝑠)(𝑠+1)(𝑠−1)𝑦−1−𝑟^𝑛+1+(𝑟+1+𝑠)(𝑠−1)(𝑠)𝑦2𝑟^𝑛+(1+𝑠)(𝑠−1)(𝑠)()𝑦−𝑟−1−𝑟)(−𝑟−2^𝑛−1. (7) For VSVO-BBDF of order 𝑃4(𝑃=4)𝑃(𝑥)=𝑃𝑥_𝑛+1=(𝑠)+𝑠ℎ2𝑟+1+𝑠)(𝑟+1+𝑠)(1+𝑠)(2()𝑦2𝑟+2)(𝑟+2^𝑛+2+(2𝑟+1+𝑠)(𝑟+1+𝑠)(1+𝑠)(𝑠−1)−()𝑦2𝑟+1)(𝑟+1^𝑛+1+(2𝑟+1+𝑠)(𝑟+1+𝑠)(𝑠)(𝑠−1)4𝑟²𝑦_𝑛+(2𝑟+1+𝑠)(1+𝑠)(𝑠)(𝑠−1)−𝑟²()𝑦−𝑟−1)(−𝑟−2^𝑛−1+(𝑟+1+𝑠)(1+𝑠)(𝑠)(𝑠−1)2𝑟²()𝑦−2𝑟−1)(−2𝑟−2^𝑛−2. (8) For VSVO-BBDF of order 𝑃5(𝑃=5)𝑃(𝑥)=𝑃𝑥_𝑛+1=()𝑠×2()+𝑠ℎ𝑞+2𝑟+1+𝑠2𝑟+1+𝑠)(𝑟+1+𝑠)(1+𝑠𝑞+2𝑟+22𝑟+2)(𝑟+2⁻¹𝑦_𝑛+2+)×()×−()𝑞+2𝑟+1+𝑠(2𝑟+1+𝑠)(𝑟+1+𝑠1+𝑠)(𝑠−1𝑞+2𝑟+12𝑟+1)(𝑟+1⁻¹𝑦_𝑛+1+𝑞+2𝑟+1+𝑠(2𝑟+1+𝑠)(𝑟+1+𝑠)𝑠(𝑠−1)4𝑞+2𝑟𝑟²𝑦_𝑛+𝑞+2𝑟+1+𝑠(2𝑟+1+𝑠)(1+𝑠)𝑠(𝑠−1)−𝑟²()𝑦𝑞+𝑟−𝑟−1)(−𝑟−2^𝑛−1+𝑞+2𝑟+1+𝑠(𝑟+1+𝑠)(1+𝑠)𝑠(𝑠−1)2𝑞𝑟²()𝑦−2𝑟−1)(−2𝑟−2^𝑛−2+)𝑠(××((2𝑟+1+𝑠)(𝑟+1+𝑠)(1+𝑠𝑠−1))−𝑞−𝑞−𝑟−𝑞−2𝑟−𝑞−2𝑟−1−𝑞−2𝑟−2⁻¹𝑦_𝑛−3. (9)

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 ∑^𝑘_𝑗=0𝛼_𝑗𝑦_𝑛+𝑗=ℎ∑^𝑘_𝑗=0𝛽_𝑗𝑓_𝑛+𝑗, where 𝑦_𝑛+𝑗≈𝑦(𝑥_𝑛+𝑗) and 𝑓_𝑛+𝑗≡𝑓(𝑥_𝑛+𝑗,𝑦_𝑛+𝑗), coefficients 𝛼_𝑗,  𝛽_𝑗 are suitably chosen constants subject to conditions 𝛼_𝑘=1, |𝛼₀|+|𝛽₀|≠0, and 𝑘 is defined as the order of the particular method employed. This method is said to be explicit if 𝛽_𝑘=0 and implicit otherwise.
Substituting 𝑠=0 and 𝑠=1 gives the corrector for the first and second points, respectively. Therefore by letting 𝑟=1,𝑞=1, 𝑟=2,𝑞=2 and 𝑟=1,𝑞=10/19. This produced the following coefficients (Tables 4, 5, and 6), as in Definition 1 for the first and second points of VSVO-BBDF method.

Table 4 

Coefficients for the first and second point of VSVO-BBDF when 𝑟=1, and 𝑟=1,𝑞=1.

Table 5 

Coefficients for the first and second point of VSVO-BBDF when 𝑟=2, and 𝑟=2,𝑞=2.

Table 6 

Coefficients for the first and second point of VSVO-BBDF when 𝑟=10/19, and 𝑟=1,𝑞=10/19.

3. 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 𝐿[𝑧(𝑥);ℎ]=^𝑗_𝑘=0𝛼_𝑘𝑧(𝑥+𝑘ℎ)−ℎ𝛽_𝑘𝑧^′(𝑥+𝑘ℎ) (10) are said to be of order 𝑝 if 𝑐_𝑜=𝑐₁=⋯=𝑐_𝑝=0, 𝐶_𝑝+1≠0. The general form for the constant 𝐶_𝑞 is defined as 𝐶_𝑞=^𝑗_𝑘=0𝑘^𝑞𝛼_𝑘−1!𝑘𝑞−1^𝑞−1𝛽_𝑘,𝑞=2,3,…𝑝+1. (11) 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 𝐿[𝑧(𝑥);ℎ]=^𝑗_𝑘=0𝐴_𝑘𝑧(𝑥+𝑘ℎ)−ℎ𝐵_𝑘𝑧^′(𝑥+𝑘ℎ). (12) And the general form for the constant 𝐶_𝑞 is defined as 𝐶_𝑞=^𝑗_𝑘=0𝑘^𝑞𝐴_𝑘−1!𝑘𝑞−1^𝑞−1𝐵_𝑘,𝑞=2,3,…𝑝+1. (13)𝐴_𝑘 is equal to the coefficients of 𝑦_𝑘, where 𝑘=𝑛=(𝑃−2),…,𝑛+1,𝑛+2 and 𝑃=3,4,5.

4. Implementation of VSVO-BBDF Method

Throughout this section, we illustrate the effect of Newton-type scheme which in a general form of 𝑦^(𝑖+1)_{𝑛+1,𝑛+2}−𝑦^(𝑖)_{𝑛+1,𝑛+2}=−(𝐼−𝐴)−ℎ𝐵𝜕𝐹𝑦𝜕𝑦^{(𝑖)𝑛+1,𝑛+2−1}×(𝐼−𝐴)𝑦^(𝑖)_{𝑛+1,𝑛+2}−ℎ𝐵𝐹𝑦^(𝑖)_{𝑛+1,𝑛+2}−𝜉. (14) The general form of VSVO-BBDF method is 𝑦_𝑛+1=𝛼₁ℎ𝑓_𝑛+1+𝜃₁𝑦_𝑛+2+𝜓₁,𝑦_𝑛+2=𝛼₁ℎ𝑓_𝑛+2+𝜃₁𝑦_𝑛+1+𝜓₂, (15) with 𝜓₁ and 𝜓₂ are the back values. By setting, 𝐼=1001,𝑦_{𝑛+1,𝑛+2}=𝑦_𝑛+1𝑦_𝑛+2𝛼,𝐵=¹00𝛼₂,𝐹_{𝑛+1,𝑛+2}=𝑓_𝑛+1𝑓_𝑛+2,𝜉_{𝑛+1,𝑛+2}=𝜓₁𝜓₂. (16) Equation (15) in matrix-vector form is equivalent to (𝐼−𝐴)𝑦_{𝑛+1,𝑛+2}=ℎ𝐵𝐹_{𝑛+1,𝑛+2}+𝜉_{𝑛+1,𝑛+2}. (17) Equation (17) is simplified as 𝑓_{𝑛+1,𝑛+2}=(𝐼−𝐴)𝑦_{𝑛+1,𝑛+2}−ℎ𝐵𝐹_{𝑛+1,𝑛+2}−𝜉_{𝑛+1,𝑛+2}=0. (18)

Newton iteration is performed to the system 𝑓_{𝑛+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 𝐸^(𝑖+1)_1,2=𝑦^(𝑖+1)_{𝑛+1,𝑛+2}−𝑦^(𝑖)_{𝑛+1,𝑛+2},𝐴=(𝐼−𝐴)−ℎ𝐵𝜕𝐹𝑦𝜕𝑌^{(𝑖)𝑛+1,𝑛+2},𝐵=(𝐼−𝐴)𝑦^(𝑖)_{𝑛+1,𝑛+2}−ℎ𝐵𝐹𝑦^(𝑖)_{𝑛+1,𝑛+2}−𝜉_{𝑛+1,𝑛+2}. (19)

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 1. Compute the values for 𝑒^(𝑖+1)_{𝑛+1,𝑛+2}=𝐴⁻¹𝐵, where 𝐴=1−𝛼₁ℎ𝜕𝐹_𝑛+1𝜕𝑦_𝑛+1𝜃₁𝜃₂1−𝛼₂ℎ𝜕𝐹_𝑛+2𝜕𝑦_𝑛+2,𝑦𝐵=^(𝑖)𝑛+1+𝛼₁ℎ𝐹^(𝑖)_𝑛+1+𝜉₁𝜃₁𝑦^(𝑖)_𝑛+2𝜃₂𝑦^(𝑖)_𝑛+1𝑦^(𝑖)_𝑛+2+𝛼₂ℎ𝐹^(𝑖)_𝑛+2+𝜉₂. (20)

Step 2. Calculate the corrected value for 𝑦^(𝑖+1)_{𝑛+1,𝑛+2} with the value 𝑒^(𝑖+1)_{𝑛+1,𝑛+2} from  (1).

Step 3. Solve 𝑒^(𝑖+1)_{𝑛+1,𝑛+2}=𝐴⁻¹𝐵 for the second iteration.

Step 4. Let 𝑦^(𝑖+1)_{𝑛+1,𝑛+2} be equal to the second iteration of 𝑒^(𝑖+1)_{𝑛+1,𝑛+2}=𝐴⁻¹𝐵.
The error is defined as error=𝑦_approximate−𝑦_actual, (21) and the maximum error is defined as MAXE=max_1<𝑖<TSmax_1<𝑖<𝑛(error), (22) 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 VSVO-BBDF Method

To begin this section, we provide some definitions for the stability conditions of VSVO-BBDF method as in (15).

Definition 3. A method is said to be zero stable if the roots of the polynomial 𝜌(𝑧)=𝜋(𝑧,0) satisfy the condition 𝑧₁=1≤|𝑧₂|,  𝑧₂≠1.

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 𝜋(𝑟,ℎ𝜆)=𝜌(𝑟)−ℎ𝜆𝜎(𝑟)=0 satisfy |𝑟_𝑠|<1,  𝑠=1,2,…,𝑘.

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 𝑅𝜆<0, where ℜ₁={ℎ𝜆∣𝑅ℎ𝜆<−𝑡}, ℜ₂={ℎ𝜆∣−𝑡≤𝑅ℎ𝜆≤𝑢,−𝑣≤∐ℎ𝜆<𝑣}, and 𝑡,  𝑢, and 𝑣 are positive constants.

The stability polynomial, 𝑅(𝑡,ℎ) associated with the method of (15), is given by det(𝐴𝑡²−𝐵𝑡−𝐶), while the absolute stability region of this method in the ℎ𝜆 plane is determined by solving det(𝐴𝑡²−𝐵𝑡−𝐶=0). Consequently, to determine zero stable, we substitute ℎ=ℎ𝜆=0 to 𝑅(𝑡,ℎ) for each order of VSVO-BBDF method.

Hence, the roots for the three different step size and order (𝑃) selections obtained by using Maple are listed in Table 7.

Table 7 

Roots for different step sizes and orders.

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.

(a)

(b)

(c)

Figure 2 

Stability regions for order (a)  𝑃3, (b)  𝑃4, and (c)  𝑃5.

Based on Figure 2, VSVO-BBDF method possesses the region absolute stability, which contains almost the whole of the half-plane Re(ℎ𝜆)<0.

6. 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 LTE_𝑘=𝑦^(𝑃+1)_𝑛+2−𝑦^(𝑃)_𝑛+2,𝑃=3,4,5, (23) 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 ℎ_𝑃−1=ℎ_old×TOLLTE_𝑃−1^1/𝑃,ℎ_𝑃=ℎ_old×TOLLTE_𝑃^1/(𝑃+1),ℎ_𝑃+1=ℎ_old×TOLLTE_𝑃+1^1/(𝑃+2), (24) 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 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, (𝑟=1), (𝑟=2) or (𝑟=10/19).

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, (𝑟=1)(𝑟=2) or (𝑟=10/19).

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.

7. 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, (25) 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, [𝑦₀], 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⁻², (b) 10⁻⁴, and (c) 10⁻⁶.

The test problems and solutions are listed below.

Problem 1. Consider, 𝑦^′=−100𝑦−𝑥+1,𝑦(0)=1,0≤𝑥≤10, (26) with solution 𝑦(𝑥)=𝑒^−100𝑥+𝑥. (27)

Problem 2. Consider, 𝑦^′₁=−1002𝑦₁+1000𝑦₂²,𝑦₁(0)𝑦=1,0≤𝑥≤10,_′2=𝑦₁−𝑦₂1+𝑦₂,𝑦₂(0)=1, (28) with solution 𝑦₁=𝑒^−2𝑥,𝑦₂=𝑒^−𝑥. (29)
See Kaps [18].

Problem 3. Consider, 𝑦^′₁=−2𝑦₁+𝑦₂+2sin(𝑥),𝑦₁(0)𝑦=2,0≤𝑥≤10,_′2=998𝑦₁−999𝑦₂+999(cos(𝑥)−sin(𝑥)),𝑦₂(0)=3, (30) with solution 𝑦₁(𝑥)=2𝑒^−𝑥(𝑥),𝑦+sin₂(𝑥)=2𝑒^−𝑥+cos(𝑥). (31)
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, 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.

Table 8 

Numerical results for Problem 1.

(a)

(b)

Figure 3 

Efficiency curves for Problem 1.

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.

Table 9 

Numerical results for Problem 2.

(a)

(b)

Figure 4 

Efficiency curves for Problem 2.

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.

Table 10 

Numerical results for Problem 3.

(a)

(b)

Figure 5 

Efficiency curves for Problem 3.

8. 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.

Step 1. Let order = 3.

Step 2. Calculate initial step size (ℎ) and 𝜀.

Step 3. Calculate initial array 𝑦₁,…𝑦_𝑖,  𝑖=order−1.

Step 4. Calculate the predictor values for 𝑦_𝑖+1,  𝑦_𝑖+2,  𝑖=order−1.

Step 5. Calculate 𝐹_𝑖 for 𝑖=1,…,order+1.

Step 6. Calculate the corrector values for 𝑦_𝑖+1,  𝑦_𝑖+2,  𝑖=order−1.

Step 7. Calculate 𝐹_𝑖 for 𝑖=1,…,order+1.

Implementation of VSVO-BBDF.

Step 8. Calculate Jacobian (JACBN).

Step 9. Perform LU-Decomposition to calculate the increment 𝐸.
Consider, 𝐸^(𝑖+1)_1,2=𝑦^(𝑖+1)_{𝑛+1,𝑛+2}−𝑦^(𝑖+1)_{𝑛+1,𝑛+2},𝑛=order−1. (A.1) Test for convergence.

Step 10. If 𝐸>𝜀, convergence is false. Else go to Step 11.

Step 11. Calculate LTE_𝑖,  𝑖=order.

Change order.

Step 12. Let order equal 4 and order equal 5.

Step 13. Repeat from Step 2 until Step 11 for each order correspondingly.

Step size and order control.

Step 14. Calculate ℎ_max for choosing the new step size.
Consider, ℎ_𝑘=ℎ_old×TOLLTE_2,𝑘^1/(𝑘+1). (A.2)
The final step size after a successful step (LTE<𝜀) is given by ℎ_new=𝑐×ℎ_old×TOLLTE^1/𝑝. (A.3)
The final step size after a failure step (LTE>𝜀) is given by ℎ_new≥12×ℎ_old. (A.4)

Step 15. If ℎ_new≥1.9×ℎ_old then ℎ_new=1.9×ℎ_old.
Else ℎ_new=ℎ_old.

Step 16. Update values: 𝑦_𝑖=𝑦_𝑖+2,𝑥_𝑖=𝑥_𝑖+2for𝑖=1,…,order−1. (A.5)

Step 17. Repeat Step 1 to Step 16 until  ℎ  = end point.

The Abbreviations Used in Figures  3–5 and Tables  8–10

TSs:	The total number of steps taken
TOL:	The initial value for the local error estimate
MAXE:	The maximum error
AVEE:	The average error
MTD:	The method used
TIME:	The total execution time (seconds).