# A Numerical Algorithm for Solving Stiff Ordinary Differential Equations

Authors: S. A. M. Yatim; Z. B. Ibrahim; K. I. Othman; M. B. Suleiman
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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)y′=f(x,y)
with given initial values y(a)=y0 in the given interval x∈[a,b].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 yn+1,yn+2,yn+3,…,yn+k 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 Section2, 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 ofyn+1 and yn+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 (P3,  P4, and  P5) are distinguished by the number of back values contained in total blocks. The ratio distance between current (xn) and previous step (xn-1) is represented as r and q 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, q 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 (P3–P5).We find approximating polynomialsPk(x), by means of a  k-degree polynomial interpolating the values of y at given points, that are (xn-3,yn-3),(xn-2,yn-2),(xn-1,yn-1),…,(xn+2,yn+2):
(2)Pk=∑j=0ky(xn+1-j)·Lk,j(x),
where
(3)Lk,j(x)=∏i=0i≠jk(x-xn+1-i)(xn+1-j-xn+1-i)foreachj=0,1,…,k.Predictors for the first pointyn+1p and second point yn+2p were computed using back values as the interpolating points. The resulting Lagrange polynomial for each order was given as follows.For VSVO-BBDF of orderP3(P=3)(4)P(x)=P(xn+1+sh)=(r+1+s)ryn+(1+s)-ryn-1.
For VSVO-BBDF of order P4(P=4)(5)P(x)=P(xn+1+sh)=(2r+1+s)(r+1+s)2r2yn+(2r+1+s)(1+s)-r2yn-1+(r+1+s)(1+s)2r2yn-2.
For VSVO-BBDF of order P5(P=5)(6)P(x)=P(xn+1+sh)=(q+2r+1+s)(2r+1+s)(r+1+s)(2(q+2r))r2yn+(q+2r+1+s)(2r+1+s)(1+s)-(q+r)r2yn-1+(q+2r+1+s)(r+1+s)(1+s)(2qr2)yn-2+(2r+1+s)(r+1+s)(1+s)-q(-q-r)(-q-2r)yn-3.Substitutings=0 and s=1 gives the predictor for the first and second points, respectively. Therefore by letting r=1,q=1, r=2,q=2 and r=1,q=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 whenr=1, and r=1,q=1.
Order
y
n
y
n
-
1
y
n
-
2
y
n
-
3
P
=
3
h
f
n
+
1
2
-
1
h
f
n
+
2
3
-
2
P
=
4
h
f
n
+
1
3
-
3
1
h
f
n
+
2
6
-
8
3
P
=
5
h
f
n
+
1
4
-
6
4
-
1
h
f
n
+
2
10
-
20
15
-
4Table 2
Coefficients for the first and second point of predictor formulae for VSVO-BBDF whenr=2, and r=2,q=2.
Order
y
n
y
n
-
1
y
n
-
2
y
n
-
3
P
=
3
h
f
n
+
1
3
2
-
1
2
h
f
n
+
2
2
-
1
P
=
4
h
f
n
+
1
15
8
-
5
4
3
8
h
f
n
+
2
3
-
3
1
P
=
5
h
f
n
+
1
35
16
-
35
16
21
16
-
5
16
h
f
n
+
2
-
1
4
-
6
4Table 3
Coefficients for the first and second point of predictor formulae for VSVO-BBDF whenr=10/19, and r=1,q=10/19.
Order
y
n
y
n
-
1
y
n
-
2
y
n
-
3
P
=
3
h
f
n
+
1
29
10
-
19
10
h
f
n
+
2
24
5
-
19
5
P
=
4
h
f
n
+
1
1131
200
-
741
100
-
551
200
h
f
n
+
2
348
25
551
25
228
25
P
=
5
h
f
n
+
1
67
16
-
201
29
67
10
-
67
32
h
f
n
+
2
43
4
-
688
29
129
5
-
6859
580The interpolating polynomial of the functiony(x) using Lagrange polynomial in (2) gives the following corrector for the first point yn+1p and second point yn+2p. The resulting Lagrange polynomial for each order was given as follows.For VSVO-BBDF of orderP3(P=3)(7)P(x)=P(xn+1+sh)=(r+1+s)(s+1)(s)2r+4yn+2+(r+1+s)(s+1)(s-1)-1-ryn+1+(r+1+s)(s-1)(s)2ryn+(1+s)(s-1)(s)-r(-1-r)(-r-2)yn-1.
For VSVO-BBDF of order P4(P=4)(8)P(x)=P(xn+1+sh)=(2r+1+s)(r+1+s)(1+s)(s)2(2r+2)(r+2)yn+2+(2r+1+s)(r+1+s)(1+s)(s-1)-(2r+1)(r+1)yn+1+(2r+1+s)(r+1+s)(s)(s-1)4r2yn+(2r+1+s)(1+s)(s)(s-1)-r2(-r-1)(-r-2)yn-1+(r+1+s)(1+s)(s)(s-1)2r2(-2r-1)(-2r-2)yn-2.
For VSVO-BBDF of order P5(P=5)(9)P(x)=P(xn+1+sh)=((2(q+2r+2)(2r+2)(r+2))-1((q+2r+1+s)(2r+1+s)(r+1+s)(1+s)s)×(2(q+2r+2)(2r+2)(r+2))-1)yn+2+((2(q+2r+2)(2r+2)(r+2))-1((q+2r+1+s)(2r+1+s)(r+1+s)×(1+s)(s-1)(q+2r+1+s))×(-(q+2r+1)(2r+1)(r+1))-1)yn+1+(q+2r+1+s)(2r+1+s)(r+1+s)s(s-1)4(q+2r)r2yn+(q+2r+1+s)(2r+1+s)(1+s)s(s-1)-r2(q+r)(-r-1)(-r-2)yn-1+(q+2r+1+s)(r+1+s)(1+s)s(s-1)2qr2(-2r-1)(-2r-2)yn-2+((2(q+2r+2)(2r+2)(r+2))-1((2r+1+s)(r+1+s)(1+s)s(s-1))×((-q-2r-1)(-q-2r-2)-q(-q-r)(-q-2r)×(-q-2r-1)(-q-2r-2))-1)yn-3.Linear Multistep Method (LMM) given in [15] is given in the definition below.Definition 1.
The lineark-step method can be represented in standard form by an equation ∑j=0kαjyn+j=h∑j=0kβjfn+j, where yn+j≈y(xn+j) and fn+j≡f(xn+j,yn+j), coefficients αj,  βj are suitably chosen constants subject to conditions αk=1, |α0|+|β0|≠0, and k is defined as the order of the particular method employed. This method is said to be explicit if βk=0 and implicit otherwise.
Substitutings=0 and s=1 gives the corrector for the first and second points, respectively. Therefore by letting r=1,q=1, r=2,q=2 and r=1,q=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 whenr=1, and r=1,q=1.
Order
α
k
,
0
α
k
,
1
α
k
,
2
α
k
,
3
α
k
,
4
α
k
,
5
P
=
3
h
f
n
+
1
1
6
-
1
1
2
1
3
h
f
n
+
2
-
1
3
3
2
-
3
11
6
P
=
4
h
f
n
+
1
-
1
12
1
2
-
3
2
5
6
1
4
h
f
n
+
2
1
4
-
4
3
3
-
4
25
12
P
=
5
h
f
n
+
1
1
12
-
1
2
7
6
-
1
3
-
5
4
5
6
h
f
n
+
2
-
5
6
61
12
-
13
107
6
-
77
6
15
4Table 5
Coefficients for the first and second point of VSVO-BBDF whenr=2, and r=2,q=2.
Order
α
^
k
,
0
α
^
k
,
1
α
^
k
,
2
α
^
k
,
3
α
^
k
,
4
α
^
k
,
5
P
=
3
h
f
n
+
1
1
24
-
3
4
1
3
3
8
h
f
n
+
2
-
1
12
1
-
8
3
7
4
P
=
4
h
f
n
+
1
-
1
80
5
48
-
15
16
8
15
5
16
h
f
n
+
2
1
30
-
1
4
3
2
-
16
5
23
12
P
=
5
h
f
n
+
1
1
168
-
1
24
1
8
17
24
-
12
7
11
12
h
f
n
+
2
-
23
336
1
2
-
43
24
37
6
-
160
21
45
16Table 6
Coefficients for the first and second point of VSVO-BBDF whenr=10/19, and r=1,q=10/19.
Order
α
k
,
0
α
k
,
1
α
k
,
2
α
k
,
3
α
k
,
4
α
k
,
5
P
=
3
h
f
n
+
1
6859
13920
-
29
20
19
29
29
96
h
f
n
+
2
-
6859
6960
12
5
-
96
29
91
48
P
=
4
h
f
n
+
1
-
6854
15600
89167
46400
-
1131
400
1292
1131
13
64
h
f
n
+
2
13718
9425
-
6859
1200
174
25
-
64
13
3095
1392
P
=
5
h
f
n
+
1
2476099
8020704
-
7
8
124
87
-
47
96
-
80
67
851
1032
h
f
n
+
2
-
12380495
4010352
53
6
-
452
29
931
48
-
2692
201
165
43


## 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 L defined by
(10)L[z(x);h]=∑k=0j[αkz(x+kh)-hβkz′(x+kh)]
are said to be of order p if co=c1=⋯=cp=0, Cp+1≠0. The general form for the constant Cq is defined as
(11)Cq=∑k=0j[kqαk-1(q-1)!kq-1βk],q=2,3,…p+1.
Consequently, BBDF method can be represented in standard form by an equation ∑j=0kAjyn+j=h∑j=0kBjfn+j, where Aj and Bj are r by r matrices with elements al,m and bl,m for l,m=1,2,…,r. Since VSVO-BBDF for variable order (P) is a block method, we extend the Definition 2 in the form of
(12)L[z(x);h]=∑k=0j[Akz(x+kh)-hBkz′(x+kh)].
And the general form for the constant Cq is defined as
(13)Cq=∑k=0j[kqAk-1(q-1)!kq-1Bk],q=2,3,…p+1.Ak is equal to the coefficients of yk, where k=n=(P-2),…,n+1,n+2 and P=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(14)yn+1,n+2(i+1)-yn+1,n+2(i)=-[(I-A)-hB∂F∂y(yn+1,n+2(i))]-1×[(I-A)yn+1,n+2(i)-hBF(yn+1,n+2(i))-ξ].
The general form of VSVO-BBDF method is
(15)yn+1=α1hfn+1+θ1yn+2+ψ1,yn+2=α1hfn+2+θ1yn+1+ψ2,
with ψ1 and ψ2 are the back values. By setting,
(16)I=[1001],yn+1,n+2=[yn+1yn+2],B=[α100α2],Fn+1,n+2=[fn+1fn+2],ξn+1,n+2=[ψ1ψ2].
Equation (15) in matrix-vector form is equivalent to
(17)(I-A)yn+1,n+2=hBFn+1,n+2+ξn+1,n+2.
Equation (17) is simplified as
(18)f^n+1,n+2=(I-A)yn+1,n+2-hBFn+1,n+2-ξn+1,n+2=0.Newton iteration is performed to the systemf^n+1,n+2=0; by taking the analogous form of (14) where Jn+1,n+2=(∂F/∂Y)(Yn+1,n+2(i)) is the Jacobian matrix of F with respect to Y.  Equation (14) is separated to three different matrices denoted as
(19)E1,2(i+1)=yn+1,n+2(i+1)-yn+1,n+2(i),A^=(I-A)-hB∂F∂Y(yn+1,n+2(i)),B^=(I-A)yn+1,n+2(i)-hBF(yn+1,n+2(i))-ξn+1,n+2.Two-stage Newton iteration works to find the approximating solution to (1) with two simplified strategies based on evaluating the Jacobian (Jn+1,n+2) and LU factorization of A^ [16]. Two-stage Newton iteration is carried out as follows.Step 1.
Compute the values foren+1,n+2(i+1)=A^-1B^, where
(20)A^=[1-α1h∂Fn+1∂yn+1θ1θ21-α2h∂Fn+2∂yn+2],B^=[yn+1(i)+α1hFn+1(i)+ξ1θ1yn+2(i)θ2yn+1(i)yn+2(i)+α2hFn+2(i)+ξ2].Step 2.
Calculate the corrected value foryn+1,n+2(i+1) with the value en+1,n+2(i+1) from  (1).Step 3.
Solveen+1,n+2(i+1)=A^-1B^ for the second iteration.Step 4.
Letyn+1,n+2(i+1) be equal to the second iteration of en+1,n+2(i+1)=A^-1B^.
The error is defined as(21)error=|yapproximate-yactual|,
and the maximum error is defined as
(22)MAXE=(max1<i<TS(max1<i<n(error))),
where TS is the total steps and n 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ρ(z)=π(z,0) satisfy the condition z1=1≤|z2|,  z2≠1.Definition 4.
A method is said to be absolutely stable in a regionR for a given hλ if for that hλ all the roots rs of the stability polynomial π(r,hλ)=ρ(r)-hλσ(r)=0 satisfy |rs|<1,  s=1,2,…,k.Definition 5.
A method is said to be stiffly stable ifℜ1 and ℜ2 are contained in the absolute stability region, and it is accurate for all h∈ℜ2 when applied to the scalar test equation y′=λy; λ is a complex constant with Rλ<0, where ℜ1={hλ∣Rhλ<-t}, ℜ2={hλ∣-t≤Rhλ≤u,-v≤∐hλ<v}, and t,  u, and v are positive constants.The stability polynomial,R(t,h^) associated with the method of (15), is given by det(At2-Bt-C), while the absolute stability region of this method in the hλ plane is determined by solving det(At2-Bt-C=0). Consequently, to determine zero stable, we substitute h^=hλ=0 to R(t,h^) for each order of VSVO-BBDF method.Hence, the roots for the three different step size and order(P) selections obtained by using Maple are listed in Table 7.Table 7
Roots for different step sizes and orders.
Roots for different step size
P
=
3
P
=
4
P
=
5
r
=
1, and r=1,q=1
t
=
-
0.4347826087
e
-
1, and t=1
t
=
-
0.24414201370,
t
=
-
0.9671746307
e
-
1,
t
=
0.02079175991,
t
=
0.01195305716,
andt=1
andt=1
r
=
2, and r=2,q=2
t
=
-
0.01315789474, and t=1
t
=
-
0.052708171410,
t
=
-
0.06178960900,
t
=
0.003257621961,
t
=
-
0.001737052270,
andt=1
andt=1
r
=
10
19, and r=1,q=1019
t
=
-
0.1098846524, and t=-1
t
=
-
0.93455113330,
t
=
-
0.4375589404,
t
=
0.08581158625,
t
=
-
0.007529578675,
andt=1
andt=1The stability region was given by the set of points determined by the boundaryt=eiθ,  0≤θ≤2π. We obtained the stability region by finding the region for which |t|<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  P3,  P4, and  P5  of VSVO-BBDF method, respectively. The stability regions lie outside the closed region for each case.Stability regions for order (a)P3, (b)  P4, and (c)  P5.
(a)
(b)
(c)Based on Figure2, VSVO-BBDF method possesses the region absolute stability, which contains almost the whole of the half-plane Re(hλ)<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 P-1,  P,  and  P+1  are selected depending on the occurrence of every successful step. Consequently, the new step size hnew 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(23)LTEk=yn+2(P+1)-yn+2(P),P=3,4,5,
where yn+2(P+1) is the (P+1)th order method and yn+2(P) is the kth order method. By finding the LTEs, the maximum step size is defined as
(24)hP-1=hold×(TOLLTEP-1)1/P,hP=hold×(TOLLTEP)1/(P+1),hP+1=hold×(TOLLTEP+1)1/(P+2),
where hold is the step size from previous block and hmax 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,(r=1), (r=2) or (r=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,(r=1)(q=1),(r=2)(q=2), or (r=1)(q=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,(r=1)(r=2) or (r=10/19).The successful step is dependent on the conditionLTE<TOL. If this condition fails, the values of  yn+1,  yn+2 are rejected, and the current step is reiterated with step size selection (q=2). On the contrary, the step size increment for each successful step is defined ashnew=c×hmax  and ifhnew>1.9×hold  then hnew=1.9×hold,wherec is the safety factor, p is the order of the method while hold and hnew is the step size from previous and current block, respectively. In this paper, c 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 inMATLAB version 7 in Windows XP or later. The syntax for using ode15s and ode23s is similar for every test problem. For example,
(25)[X,Y]=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, [xo,xend], init is a column vector of initial conditions, [y0], 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,(26)y′=-100(y-x)+1,y(0)=1,0≤x≤10,
with solution
(27)y(x)=e-100x+x.Problem 2.
Consider,(28)y1′=-1002y1+1000y22,y1(0)=1,0≤x≤10,y2′=y1-y2(1+y2),y2(0)=1,
with solution
(29)y1=e-2x,y2=e-x.
See Kaps [18].Problem 3.
Consider,(30)y1′=-2y1+y2+2sin(x),y1(0)=2,0≤x≤10,y2′=998y1-999y2+999(cos(x)-sin(x)),y2(0)=3,
with solution
(31)y1(x)=2e-x+sin(x),y2(x)=2e-x+cos(x).
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 Table8, 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 Problem1.
TOL
MTD
TS
AVEE
MAXE
TIME
10−2
VSVO-BBDF
21
2.9370
e
-
005
2.8298
e
-
004
0.0103
ode15s
28
1.3000
e
-
003
8.4000
e
-
003
0.0313
ode23s
19
1.0000
e
-
003
4.5000
e
-
003
0.1406
10−4
VSVO-BBDF
48
1.0716
e
-
006
3.2212
e
-
006
0.0105
ode15s
60
3.0358
e
-
005
1.6621
e
-
004
0.0156
ode32s
42
1.1285
e
-
004
2.5683
e
-
004
0.0313
10−6
VSVO-BBDF
164
1.6733
e
-
008
3.1232
e
-
008
0.0115
ode15s
100
7.2564
e
-
007
2.7506
e
-
006
0.0313
ode23s
143
7.3558
e
-
006
1.2514
e
-
005
0.0469Efficiency curves for Problem1.
(a)
(b)Again by comparing the four factors mentioned earlier, we can see VSVO-BBDF in Table9 gives the minimum maximum error for every tolerance level except for TOL  10e-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 Problem2.
TOL
MTD
STs
AVEE
MAXE
TIME
10−2
VSVO-BBDF
22
7.1459
e
-
005
2.5736
e
-
004
0.0106
ode15s
29
9.0287
e
-
004
5.2000
e
-
003
0.0781
ode23s
25
3.3626
e
-
004
1.1000
e
-
003
0.0781
10−4
VSVO-BBDF
54
7.4173
e
-
006
3.7659
e
-
004
0.0111
ode15s
55
1.7139
e
-
005
8.5506
e
-
005
0.1250
ode32s
118
1.3783
e
-
005
6.9774
e
-
005
0.2031
10−6
VSVO-BBDF
194
6.3429
e
-
009
3.2882
e
-
008
0.0143
ode15s
197
2.1320
e
-
007
1.0790
e
-
006
0.2344
ode23s
773
3.3163
e
-
007
2.8081
e
-
006
0.4531Efficiency curves for Problem2.
(a)
(b)From Table10, 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 Problem3.
TOL
MTD
TSs
AVEE
MAXE
TIME
10−2
VSVO-BBDF
35
4.6584
e
-
005
3.0045
e
-
004
0.0111
ode15s
45
6.9000
e
-
003
1.4620
e
-
002
0.0156
ode23s
137
4.0745
e
-
004
3.0000
e
-
002
0.2031
10−4
VSVO-BBDF
84
2.5775
e
-
006
1.1002
e
-
005
0.3281
ode15s
93
7.6548
e
-
005
2.7591
e
-
004
0.0469
ode32s
1211
2.4915
e
-
005
6.3075
e
-
005
0.0156
10−6
VSVO-BBDF
380
2.4244
e
-
008
8.9627
e
-
008
0.0199
ode15s
186
2.1681
e
-
006
6.1936
e
-
006
0.0781
ode23s
2829
5.3155
e
-
007
1.5667
e
-
006
2.3281Efficiency curves for Problem3.
(a)
(b)


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



---

*Source: 989381-2013-02-14.xml*