# An Accurate Block Solver for Stiff Initial Value Problems

Authors: H. Musa; M. B. Suleiman; F. Ismail; N. Senu; Z. B. Ibrahim
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

New implicit block formulae that compute solution of stiff initial value problems at two points simultaneously are derived and implemented in a variable step size mode. The strategy for changing the step size for
optimum performance involves halving, increasing by a multiple of 1.7, or maintaining the current step size. The stability analysis of the methods indicates their suitability for solving stiff problems. Numerical results are given and compared with some existing backward differentiation formula algorithms. The results indicate an improvement in terms of
accuracy.


---

## Body

## 1. Introduction

Stiff problems occur in many fields of engineering science and they represent coupled physical systems having components varying with very different time scales [1]. Some physical examples include chemical kinetics where reactions go at very different speeds and in circuit simulation where components respond at widely different time scales. A considerable research has been done on numerical methods for solving stiff initial value problems (IVPs) of the following form:
(1)y′=f(x,y)y(x0)=y0,x∈[a,b].
The most successful ones are implicit in nature. Examples include single step methods like the diagonally implicit Runge Kutta (DIRK) methods, single diagonally implicit Runge-Kutta (SDIRK), and singly implicit Runge-Kutta methods [2–4]. Also well known are methods based on Richardson extrapolation and deferred correction schemes [5–8]. Others are multistep methods like the backward differentiation formula (BDF) (see, e.g., [4]). Initially, the BDF methods are known to be A-stable only up to order 2 [9]. However, Cash in [10] proved that the barrier in [9] can be circumvented by adding extra future point to produce higher order A-stable methods. Methods of this type are known as extended backward differentiation formula [10–12]. Some block methods for stiff problems are also known to be A-stable and have order greater than two (see, e.g., [13–18]).One of the motivation of this paper is to develop a new block method based on the block backward differentiation formula [13, 14, 17]. The method is similar in fashion with that in [15] but differs greatly by the choice of the constant ρ. Furthermore, our choice of ρ=3/8 gives rise to different formulae with better stability regions than in [15], even with the choice of the same step ratio r. In the following sections of the paper, we will present the derivation of the method, its stability analysis, and implementation. Standard test problems will also be presented and a comparison of their numerical results with some known stiff solvers will be given.


## 2. Formulation of the Method

Consider a special case of the following BBDF given in [14, 17]:
(2)∑j=04‍αj,i,ryn+j-2=hβk,i,rfn+k,k=i=1,2,
where h = the variable step size used;  xn=a+nh; yn = computed solution; fn=f(xn,yn), n≥0; β0,i,r=β1,i,r=⋯=βk-1,i,r=0; k=i=1 represents the formula for the first point; k=i=2 represents the formula for the second point; and r = the value used in the step size changing strategy.In this paper, a nonzero coefficientβk-1,i,r=ρβk,i,r is introduced in (2) to have a formula of the following form:
(3)∑j=04‍αj,i,ryn+j-2=hβk,i,r(fn+k-ρfn+k-1),k=i=1,2,
where βk-1,i,r=ρβk,i,r.Define the linear difference operatorLi associated with the method (3) by
(4)Li(y(xn),h)=α0,iy(xn-2rh)+α1,iy(xn-rh)+α2,iy(xn)+α3,iy(xn+h)+α4,iy(xn+2h)-hβk,i(f(xn+kh)-ρf(xn+(k-1)h))=0,
where k=i=1 stands for the first point and k=i=2 stands for the second point. Expanding (4) using Taylor’s series and collecting like terms give a set of equations to be solved for αj,i′s and βk,i′s. For the derivation of the first point, α3,1 is normalized to 1 and for the second point, α4,2 is normalized to 1; the following 2-point implicit block method is obtained as(5)yn+1=-1-r-2ρr4r2(-2-3r+4ρr2)yn-2+1+4r+4ρr+4r2+8ρr2r2(2+r)(-2-3r+4ρr2)yn-1+-1-6r-6ρr-13r2-12ρr2-12r3+6ρr3-4r4+12ρr44r2(-2-3r+4ρr2)yn+1+5r+8r2+2ρr2+4r3+4ρr34(2+r)(-2-3r+4ρr2)yn+2+1+3r+2r2-2-3r+4ρr2ρhfn-1+3r+2r2-2-3r+4ρr2hfn+1,yn+2=8+2ρ+8r+3ρr+2r2+ρr2r2(1+2r)(-20+ρ-24r+3ρr-6r2+2ρr2)yn-2-4(4+ρ+4r+2ρr)r2(-20+ρ-24r+3ρr-6r2+2ρr2)yn-1+8+2ρ+24r+9ρr+26r2+14ρr2+12r3+9ρr3+2r4+2ρr4r2(-20+ρ-24r+3ρr-6r2+2ρr2)yn-4(16+4ρ+32r+8ρr+20r2+3ρr2+4r3)(1+2r)(-20+ρ-24r+3ρr-6r2+2ρr2)yn+1+4(2+3r+r2)-20+ρ-24r+3ρr-6r2+2ρr2hρfn+1-4(2+3r+r2)-20+ρ-24r+3ρr-6r2+2ρr2hfn+2.For optimal performance, we choose different parameters from those used in [15]. We choose ρ=3/8 and r=1,2 and 10/17 where the values of r correspond to the step size control strategy of maintaining constant, halving, and increasing the step size by a multiple of 1.7.Whenr=1 we have
(6)yn+1=1156yn-2-97yn-1+187yn-2756yn+2-914hfn+127hfn+1,yn+2=-27191yn-2+146191yn-1-342191yn+414191yn+1-36191hfn+1+96191hfn+2.
When r=2 we have
(7)yn+1=964yn-2-54yn-1+31564yn-4516yn+2-4516hfn+152hfn+1,yn+2=-733455yn-2+111691yn-1-711691yn+65283455yn+1-144691hfn+1+384691hfn+2.
And when r=10/17 we have
(8)yn+1=338997750400yn-2-23631531031800yn-1+2310687750400yn-1986982544yn+2-29977504hfn+999938hfn+1,yn+2=-848475114926355yn-2+19160780683yn-1-1307097403415yn+72652802985271yn+1-1425680683hfn+1+3801680683hfn+2.
Formulae (6), (7), and (8) compute two solution values concurrently at each step of the integration thereby having the advantage of reducing the total number of steps when a conventional nonblock method is used. Also, the formulae obtained are entirely different from those in [15]. It should be noted that when ρ=0 in (5), formula (3) reduces to (2) and we say that formula (2) is contained in (3). Hence, (3) is a superclass of (2) and contains (2) as a subclass. We will refer to our method as variable step size superclass of block backward differentiation formula (VSSBBDF). In the following section, the stability analysis of the method will be presented.


## 3. Stability of the Method

This section investigates the stability regions of the methods (6), (7), and (8). When r=1, the formulae in (6) can be represented as:
(9)(1001)(yn+1yn+2)=(011560-27191)(yn-3yn-2)+(-97187146191-342191)(yn-1yn)+(0-27564141910)(yn+1yn+2)+h(0-91400)(fn-1fn)+h(1270-3619196191)(fn+1fn+2).
Collecting like terms we have
(10)(12756-4141911)(yn+1yn+2)=(011560-27191)(yn-3yn-2)+(-97187146191-342191)(yn-1yn)+h(0-91400)(fn-1fn)+h(1270-3619196191)(fn+1fn+2).Lethλ=h-. Applying the scalar test equation,
(11)y′=λyλ<0,
to (10) gives
(12)(1-127h-2756-414191+36191h-1-96191h-)(yn+1yn+2)=(-97187-914h-146191-342191)(yn-1yn)+(011560-27191)(yn-3yn-2).Equation (12) is equivalent to
(13)AYm=BYm-1+CYm-2,n=2m,
where
(14)A=(1-127h-2756-414191+36191h-1-96191h-),B=(-97187-914h-146191-342191),C=(011560-27191),Ym=(yn+1yn+2)=(y2m+1y2m+2),Ym-1=(yn-1yn)=(y2m-1y2m)=(y2(m-1)+1y2(m-1)+2),Ym-2=(yn-3yn-2)=(y2(m-2)+1y2(m-2)+2).
The stability polynomial R(t,h-) is obtained by evaluating the determinant of (At2-Bt-C) to obtain
(15)R(t,h-)=169t5348+279t25348+765h-t22674-11385t35348-351h-t3191-162h-2t31337+10937t45348-6171h-t42674+1152h-2t41337=0.
To show that the method is zero stable, we let h-=0 in (15) and the following equation for the first characteristics polynomial is obtained as
(16)169t5348+279t25348-11385t35348+10937t45348=0.
Solving (16) for t, we have
(17)t=0,t=1,t=-0.105502,t=0.146463.
The values of t above indicate that the method is zero stable since no magnitude of the root is greater than one and the root t=1 is simple.The plot of the stability polynomial (15) is given in Figure 1. The region of stability is the region outside the circular shape. It indicates that the method (6) is A-stable since the region covers the entire negative half plane.Figure 1
Stability region of the 2-point VSSBBDF whenr=1.The same procedure is applied to obtain the stability polynomial of the method (7) (corresponding to r=2) as
(18)R(t,h-)=169t44224+11099t244224+3567h-t211056-72625t311056-5727h-t32764-405h-2t3691+4363t4691-11943h-t41382+2880h-2t4691=0.
When h-=0 in (18), the first characteristic polynomial of the method (7) becomes
(19)169t44224+11099t244224-72625t311056+4363t4691=0
and solving (19) for t gives the following roots:
(20)t=0,t=1,t=-0.0116404,t=0.051994.Again, it shows that the method (7) is zero stable.The stability plot of the polynomial (18) is given in Figure 2 and the region outside the circle indicates that the method is also A-stable.Figure 2
Stability region of the 2-point VSSBBDF whenr=2.Similarly, the stability polynomial for the method (8) (corresponding to r=10/17) is given by
(21)R(t,h-)=69347235737t302722616000-128013466997t2302722616000+6400808703h-t215136130800-21069763237t315136130800-1629222471h-t3540576100-2670327h-2t337840327+60007687t437840327-119480103h-t475680654+18988992h-2t437840327=0
and when h-=0 in (21), the first characteristics polynomial is given by
(22)69347235737t302722616000-128013466997t2302722616000-21069763237t315136130800+60007687t437840327=0.
When (22) is solved for t, the following values are obtained:
(23)t=0,t=1,t=-0.446055,t=0.32385.
Hence, the method (8) is also zero stable.The stability plot for (21) is given in Figure 3. The method (8) is therefore almost A-stable.Figure 3
Stability region of the 2-point VSSBBDF whenr=10/17.A-stability is a desirable property for numerical methods of solving stiff problems. From the forgoing analysis, the methods developed are suitable for solving stiff initial value problems. In the following section, the procedure involved in the implementation of the method will be explained.


## 4. Implementation of the Method

The implementation involves Newton’s iteration. To simplify the code, methods (6), (7), and (8) are written in the following form:
(24)yn+1=θ1yn+2+α1hfn+α2hfn+1+Ψ1,yn+2=θ2yn+1+α3hfn+1+α4hfn+2+Ψ2,
which is equivalent to the following matrix equation:
(25)(1001)(yn+1yn+2)=(0θ1θ20)(yn+1yn+2)+h(0α100)(fn-1fn)+h(α20α3α4)(fn+1fn+2)+(Ψ1Ψ2).
Equation (25) can be represented as
(26)(I-A)Y=h(B1F1+B2F2)+ζ,
where
(27)I=(1001),A=(0θ1θ20),B1=(0α100),B2=(α20α3α4),Y=(yn+1yn+2),F1=(fn-1fn),F2=(fn+1fn+2),ζ=(Ψ1Ψ2).Let(28)F=(I-A)Y-h(B1F1+B2F2)-ζ=0.
Newton’s iteration takes the form
(29)Yn+1,n+2(i+1)-Yn+1,n+2(i)=-[Fj′(yn+1,n+2(i))]-1[Fj(yn+1,n+2(i))].
Equation (32) is equivalent to
(30)Yn+1,n+2(i+1)-Yn+1,n+2(i)=-[(I-A)-hB1δF1δY(Yn+1,n+2(i))-hB2δF2δY(Yn+1,n+2(i))]-1×[(I-A)(Yn+1,n+2(i))-hB1F1-hB2F2-ζ].
We show the accuracy of the method by finding the maximum error in the following way. Let yi and y(xi) be the computed and the exact solutions of (1), respectively. The absolute error is defined by
(31)(errori)t=|(yi)t-(y(xi))t|.The maximum error is given by:(32)MAXE=max︸1≤t≤T(max︸1≤i≤N(errori)t),
where T is the total number of steps and N is the number of equations.LetYn+1(i+1) denote the (i+1)th iterate and
(33)E1,2(i+1)=Yn+1,n+2(i+1)-Yn+1,n+2(i).
Then (33) can be written as
(34)E1,2(i+1)=A¯-1B-
which is equivalent to
(35)A-E1,2(i+1)=B-,
where
(36)A-=[(I-A)-hB1δF1δY(Yn+1,n+2(i))-hB2δF2δY(Yn+1,n+2(i))]-1,B-=-[(I-A)(Yn+1,n+2(i))-hB1F1-hB2F2-ζ].
LU decomposition method is used to solve the system (35). For the different values of r,
(37)A-=(1-α2hδfn+1δyn+1-θ1-θ2-α3hδfn+1δyn+11-α4hδfn+2δyn+2),B-=(-yn+1i+θ1yn+2i+α1hfni+α2hfn+1i+Ψ-yn+2i+θ2yn+1i+α3hfn+1i+α4hfn+2i+Ψ).
When r=1(38)A-=(1-127hδfn+1δyn+12756-414191+36191hδfn+1δyn+11-96191hδfn+2δyn+2),B-=(-yn+1i-2756yn+2i-914hfni+127hfn+1i+1156yn-2i-97yn-1+187yn-yn+2i+414191yn+1i-36191hfn+1i+96191hfn+2i-27191yn-2+146191yn-1-342191yn).
When r=2(39)A-=(1-152hδfn+1δyn+14516-65283455+144691hδfn+1δyn+11-384691hδfn+2δyn+2),B-=(-yn+1i-4516yn+2i-4516hfni+152hfn+1i+964yn-2-54yn-1+31564yn-yn+2i+65283455yn+1i-144691hfn+1i+384691hfn+2i-733455yn-2+111691yn-1-711691yn)
and when r=10/17(40)A-=(1-999938hδfn+1δyn+11986982544-72652802985271+1425680683hδfn+1δyn+11-3801680683hδfn+2δyn+2),B-=(-yn+1i-1986982544yn+2i-29977504hfni+999938hfn+1i+338997750400yn-2-23631531031800yn-1+2310687750400yn-yn+2i+72652802985271yn+1i-1425680683hfn+1i+3801680683hfn+2i-848475114926355yn-2+19160780683yn-1-1307097403415yn).


### 4.1. Step Size Selection

The step size increase in this paper differs with that in [15] in the sense that an increase is made by a multiple of 1.7 of the previous step size. This is done in order to optimize total number of steps and computation time. Three techniques are employed in controlling the step size. At the beginning of the integration, an initial step size is determined and the local truncation error (LTE) is computed as follows:
(41)LTE=|yn+2k+1-yn+2k|,k=3,
where yn+2k+1 refers to the formula of order k+1 and yn+2k refers to the formula of order k.If the LTE is less than or equal to a specified error tolerance (TOL) (i.e.,LTE≤TOL), the step size (from the previous block) is maintained as constant (equivalent to the formula when r=1) and a new step size is computed as
(42)newh=c×oldh×(TOLLTE)1/p,
where p is the order of the method and c is a safety factor. Again we choose a safety factor (c=0.3) different from that in [15].Ifnewh>(1.7×oldh), then the step size h becomes h=1.7×oldh (equivalent to the formula when r=10/17). Otherwise if LTE>TOL, the step size is halved (equivalent to the formula when r=2).


## 4.1. Step Size Selection

The step size increase in this paper differs with that in [15] in the sense that an increase is made by a multiple of 1.7 of the previous step size. This is done in order to optimize total number of steps and computation time. Three techniques are employed in controlling the step size. At the beginning of the integration, an initial step size is determined and the local truncation error (LTE) is computed as follows:
(41)LTE=|yn+2k+1-yn+2k|,k=3,
where yn+2k+1 refers to the formula of order k+1 and yn+2k refers to the formula of order k.If the LTE is less than or equal to a specified error tolerance (TOL) (i.e.,LTE≤TOL), the step size (from the previous block) is maintained as constant (equivalent to the formula when r=1) and a new step size is computed as
(42)newh=c×oldh×(TOLLTE)1/p,
where p is the order of the method and c is a safety factor. Again we choose a safety factor (c=0.3) different from that in [15].Ifnewh>(1.7×oldh), then the step size h becomes h=1.7×oldh (equivalent to the formula when r=10/17). Otherwise if LTE>TOL, the step size is halved (equivalent to the formula when r=2).


## 5. Test Problems

The performance of the method is tested on the following stiff problems that are commonly found in engineering and physical sciences.(1)
The circuit problem taken from [17] is as follows:
(43)y′=-20y+24,y(0)=0,0≤x≤10.The exact solution is as follows:(44)y(x)=65-65e-20x.(2)
The cosine problem taken from [8, 15, 19] is as follows:
(45)y′=-2πsin(2πx)-1ϵ(y-cos(2πx)),y(0)=1,0≤x≤10.The problem becomes increasingly stiff asϵ→0. We choose ϵ=10-3.The exact solution is as follows:(46)y(x)=cos(2πx).(3)
We consider the system of nonlinear stiff IVP taken from [18, 20] as follows:
(47)y1′=-(ϵ-1+2)y1+ϵ-1y22,y1(0)=1,0≤x≤20y2′=y1-y2(1+y2),y2(0)=1.The solution, which is independent ofϵ, is
(48)y1(x)=e-2x,y2(x)=e-x.The eigenvalues areλ1=-1 and λ2=-(ϵ-1+2). This allows the degree of stiffness to be regulated. We choose ϵ=10-5.(4)
This physical problem is taken from [21] as follows:
(49)y1′=-21y1+19y2-20y3,y1(0)=1,0≤x≤10,y2′=19y1-21y2+20y3,y2(0)=0,y3′=40y1-40y2-40y3,y3(0)=-1.The exact solution is as follows:(50)y1(x)=0.5[e-2x+e-40x(cos40x+sin40x)],y2(x)=0.5[e-2x-e-40x(cos40x+sin40x)],y3(x)=2e-40x[-12cos40x+12sin40x].The eigenvalues are-2, -40+40i, and -40-40i.


## 6. Numerical Results

In the following are tables of results indicating the maximum error and the total steps taken to solve each of the problems given in the previous section. Each problem is solved at a prescribed tolerance, that is,10-2, 10-4, and 10-6. The results are compared with some known successful BDF based algorithms like the nonblock BDF method (Gear’s method), the VSBBDF [14] method, and the NVSBBDF [15] method. Also given are the graphs of the performance profile for the methods.The following notations are used in presenting the results in Tables1, 2, 3, and 4:TOL: Tolerance value usedMethod: Name of the methodIST: Number of successful stepsFS: Number of failed stepsTS: Total stepsMAXE: Maximum error.Table 1
Numerical results for problem (1) (the circuit problem).
TOL
Method
IST
FS
TS
MAXE
1
0
-
2
NBDF
23
0
23
2.74641
e
-
003
VSBBDF
63
0
63
2.24782
e
-
005
NVSBBDF
43
0
43
2.76983
e
-
005
VSSBBDF
44
0
44
2.34920
e
-
006
1
0
-
4
NBDF
121
0
121
3.62021
e
-
005
VSBBDF
77
0
77
2.24754
e
-
006
NVSBBDF
62
0
62
3.06417
e
-
007
VSSBBDF
70
0
70
2.73764
e
-
008
1
0
-
6
NBDF
365
0
365
3.82424
e
-
007
VSBBDF
107
0
107
2.24755
e
-
007
NVSBBDF
101
0
101
3.85316
e
-
009
VSSBBDF
142
0
142
3.45127
e
-
010Table 2
Numerical results for problem (2) (the cosine problem).
TOL
Method
IST
FS
TS
MAXE
1
0
-
2
NBDF
126
38
164
8.46687
e
-
005
VSBBDF
72
4
76
5.24918
e
-
005
NVSBBDF
65
1
66
5.16437
e
-
005
VSSBBDF
72
0
72
1.23612
e
-
005
1
0
-
4
NBDF
249
70
319
3.42565
e
-
006
VSBBDF
159
7
166
2.12255
e
-
006
NVSBBDF
116
0
116
1.54385
e
-
006
VSSBBDF
149
0
149
6.74134
e
-
007
1
0
-
6
NBDF
474
100
574
1.33439
e
-
007
VSBBDF
457
19
476
4.32393
e
-
008
NVSBBDF
398
4
402
3.11839
e
-
008
VSSBBDF
472
0
472
3.19856
e
-
009Table 3
Numerical results for the nonlinear problem (3).
TOL
Method
IST
FS
TS
MAXE
1
0
-
2
NBDF
65
4
69
9.38499
e
-
002
VSBBDF
45
0
45
6.48394
e
-
004
NVSBBDF
48
1
49
3.53673
e
-
004
VSSBBDF
49
0
49
3.77477
e
-
006
1
0
-
4
NBDF
107
7
114
1.02574
e
-
004
VSBBDF
63
0
63
6.48780
e
-
005
NVSBBDF
70
4
74
1.59269
e
-
005
VSSBBDF
80
0
80
1.17305
e
-
008
1
0
-
6
NBDF
134
10
144
1.03536
e
-
005
VSBBDF
99
0
99
6.48812
e
-
006
NVSBBDF
110
8
118
8.11971
e
-
007
VSSBBDF
169
0
169
2.94571
e
-
010Table 4
Numerical results for problem (4).
TOL
Method
IST
FS
TS
MAXE
1
0
-
2
NBDF
48
8
56
4.92270
e
-
002
VSBBDF
46
0
46
4.09523
e
-
003
NVSBBDF
48
0
48
5.56533
e
-
005
VSSBBDF
50
0
50
5.27731
e
-
006
1
0
-
4
NBDF
100
15
115
5.96300
e
-
004
VSBBDF
66
0
66
4.10262
e
-
004
NVSBBDF
73
0
73
6.47609
e
-
007
VSSBBDF
96
0
96
7.42231
e
-
008
1
0
-
6
NBDF
167
19
186
1.04720
e
-
005
VSBBDF
117
1
118
2.31132
e
-
006
NVSBBDF
142
0
142
6.61553
e
-
009
VSSBBDF
224
0
224
7.93169
e
-
010From Tables1, 2, 3, and 4, the maximum error for the VSSBBDF method is less when compared to the NBDF, VSBBDF, and the NVSBBDF methods. This is also depicted clearly in the graphs shown in Figures 4, 5, 6, and 7. To compare between the methods, we can therefore conclude that the present scheme is more accurate. This is achieved by an optimal choice of the constant ρ that gives better stability region and ensures accuracy at the same time. Furthermore, the total number of steps to complete the integration is competitive for all the problems tested.Figure 4
Graph ofLog10MAXE versus TOL for problem (1).Figure 5
Graph ofLog10MAXE versus TOL for problem (2).Figure 6
Graph ofLog10MAXE versus TOL for problem (3).Figure 7
Graph ofLog10MAXE versus TOL for problem (4).


## 7. Conclusion

A new variable step size algorithm based on block backward differentiation formula is developed for the integration of stiff problems occurring in the fields of engineering and physical sciences. The step size changing algorithm involves a strategy of increase by a multiple of 1.7, halving, or maintaining the current step size. By choosingρ=3/8, A-stable methods are obtained corresponding to the step size ratio of 1 and 2. This is an added advantage to the method because in most BDF algorithms, the method is not A-stable, but almost A-stable when the step ratio is 1. In addition, choosing ρ=0 reduced the algorithm developed to the BBDF, thus making it its superclass. The results obtained show the suitability of the method for solving stiff physical problems and indicate better accuracy when compared with some existing BDF algorithms. Thus, the algorithm developed can be an alternative solver for initial value problems arising in engineering and physical sciences.



---

*Source: 567451-2013-08-28.xml*