An Accurate Block Solver for Stiff Initial Value Problems

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.

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) 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 𝐴-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 𝐴-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 𝑟. 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]: ⁴_𝑗=0𝛼_{𝑗,𝑖,𝑟}𝑦_{𝑛+𝑗−2}=ℎ𝛽_{𝑘,𝑖,𝑟}𝑓_𝑛+𝑘,𝑘=𝑖=1,2, (2) where ℎ = the variable step size used; 𝑥_𝑛=𝑎+𝑛ℎ; 𝑦_𝑛 = computed solution; 𝑓_𝑛=𝑓(𝑥_𝑛,𝑦_𝑛), 𝑛≥0; 𝛽_0,𝑖,𝑟=𝛽_1,𝑖,𝑟=⋯=𝛽_{𝑘−1,𝑖,𝑟}=0; 𝑘=𝑖=1 represents the formula for the first point; 𝑘=𝑖=2 represents the formula for the second point; and 𝑟 = the value used in the step size changing strategy.

In this paper, a nonzero coefficient 𝛽_{𝑘−1,𝑖,𝑟}=𝜌𝛽_{𝑘,𝑖,𝑟} is introduced in (2) to have a formula of the following form: ⁴_𝑗=0𝛼_{𝑗,𝑖,𝑟}𝑦_{𝑛+𝑗−2}=ℎ𝛽_{𝑘,𝑖,𝑟}𝑓_𝑛+𝑘−𝜌𝑓_{𝑛+𝑘−1},𝑘=𝑖=1,2, (3) where 𝛽_{𝑘−1,𝑖,𝑟}=𝜌𝛽_{𝑘,𝑖,𝑟}.

Define the linear difference operator 𝐿_𝑖 associated with the method (3) by 𝐿_𝑖𝑦𝑥_𝑛,ℎ=𝛼_0,𝑖𝑦𝑥_𝑛−2𝑟ℎ+𝛼_1,𝑖𝑦𝑥_𝑛−𝑟ℎ+𝛼_2,𝑖𝑦𝑥_𝑛+𝛼_3,𝑖𝑦𝑥_𝑛+ℎ+𝛼_4,𝑖𝑦𝑥_𝑛+2ℎ−ℎ𝛽_𝑘,𝑖𝑓𝑥_𝑛+𝑘ℎ−𝜌𝑓𝑥_𝑛+(𝑘−1)ℎ=0, (4) where 𝑘=𝑖=1 stands for the first point and 𝑘=𝑖=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 𝛼_{𝑗,𝑖′𝑠} and 𝛽_{𝑘,𝑖′𝑠}. 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𝑦_𝑛+1=−1−𝑟−2𝜌𝑟4𝑟²−2−3𝑟+4𝜌𝑟²𝑦_𝑛−2+1+4𝑟+4𝜌𝑟+4𝑟²+8𝜌𝑟²𝑟²(2+𝑟)−2−3𝑟+4𝜌𝑟²𝑦_𝑛−1+−1−6𝑟−6𝜌𝑟−13𝑟²−12𝜌𝑟²−12𝑟³+6𝜌𝑟³−4𝑟⁴+12𝜌𝑟⁴4𝑟²−2−3𝑟+4𝜌𝑟²𝑦_𝑛+1+5𝑟+8𝑟²+2𝜌𝑟²+4𝑟³+4𝜌𝑟³4(2+𝑟)−2−3𝑟+4𝜌𝑟²𝑦_𝑛+2+1+3𝑟+2𝑟²−2−3𝑟+4𝜌𝑟²𝜌ℎ𝑓_𝑛−1+3𝑟+2𝑟²−2−3𝑟+4𝜌𝑟²ℎ𝑓_𝑛+1,𝑦_𝑛+2=8+2𝜌+8𝑟+3𝜌𝑟+2𝑟²+𝜌𝑟²𝑟²(1+2𝑟)−20+𝜌−24𝑟+3𝜌𝑟−6𝑟²+2𝜌𝑟²𝑦_𝑛−2−44+𝜌+4𝑟+2𝜌𝑟𝑟²−20+𝜌−24𝑟+3𝜌𝑟−6𝑟²+2𝜌𝑟²𝑦_𝑛−1+8+2𝜌+24𝑟+9𝜌𝑟+26𝑟²+14𝜌𝑟²+12𝑟³+9𝜌𝑟³+2𝑟⁴+2𝜌𝑟⁴𝑟²−20+𝜌−24𝑟+3𝜌𝑟−6𝑟²+2𝜌𝑟²𝑦_𝑛−416+4𝜌+32𝑟+8𝜌𝑟+20𝑟²+3𝜌𝑟²+4𝑟³(1+2𝑟)−20+𝜌−24𝑟+3𝜌𝑟−6𝑟²+2𝜌𝑟²𝑦_𝑛+1+42+3𝑟+𝑟²−20+𝜌−24𝑟+3𝜌𝑟−6𝑟²+2𝜌𝑟²ℎ𝜌𝑓_𝑛+1−42+3𝑟+𝑟²−20+𝜌−24𝑟+3𝜌𝑟−6𝑟²+2𝜌𝑟²ℎ𝑓_𝑛+2. (5)For optimal performance, we choose different parameters from those used in [15]. We choose 𝜌=3/8 and 𝑟=1,2 and 10/17 where the values of 𝑟 correspond to the step size control strategy of maintaining constant, halving, and increasing the step size by a multiple of 1.7.

When 𝑟=1 we have 𝑦_𝑛+1=11𝑦56^𝑛−2−97𝑦_𝑛−1+187𝑦_𝑛−27𝑦56^𝑛+2−914ℎ𝑓_𝑛+127ℎ𝑓_𝑛+1,𝑦_𝑛+2=−27𝑦191^𝑛−2+146𝑦191^𝑛−1−342𝑦191^𝑛+414𝑦191^𝑛+1−36191ℎ𝑓_𝑛+1+96191ℎ𝑓_𝑛+2. (6) When 𝑟=2 we have 𝑦_𝑛+1=9𝑦64^𝑛−2−54𝑦_𝑛−1+315𝑦64^𝑛−45𝑦16^𝑛+2−4516ℎ𝑓_𝑛+152ℎ𝑓_𝑛+1,𝑦_𝑛+2=−73𝑦3455^𝑛−2+111𝑦691^𝑛−1−711𝑦691^𝑛+6528𝑦3455^𝑛+1−144691ℎ𝑓_𝑛+1+384691ℎ𝑓_𝑛+2. (7) And when 𝑟=10/17 we have 𝑦_𝑛+1=338997𝑦750400^𝑛−2−2363153𝑦1031800^𝑛−1+2310687𝑦750400^𝑛−19869𝑦82544^𝑛+2−29977504ℎ𝑓_𝑛+999938ℎ𝑓_𝑛+1,𝑦_𝑛+2=−8484751𝑦14926355^𝑛−2+191607𝑦80683^𝑛−1−1307097𝑦403415^𝑛+7265280𝑦2985271^𝑛+1−1425680683ℎ𝑓_𝑛+1[phone]ℎ𝑓_𝑛+2. (8) 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 𝑟=1, the formulae in (6) can be represented as: 10𝑦01𝑦^𝑛+1_𝑛+2=011560−27𝑦191^𝑛−3𝑦_𝑛−2+−97187146−191342𝑦191^𝑛−1𝑦_𝑛+0−27564140𝑦191^𝑛+1𝑦_𝑛+29+ℎ0−𝑓1400𝑓^𝑛−1_𝑛+ℎ1270−3619196𝑓191^𝑛+1𝑓_𝑛+2. (9) Collecting like terms we have 127−564141𝑦191^𝑛+1𝑦_𝑛+2=011560−27𝑦191^𝑛−3𝑦_𝑛−2+−97187146−191342𝑦191^𝑛−1𝑦_𝑛9+ℎ0−𝑓1400𝑓^𝑛−1_𝑛+ℎ1270−3619196𝑓191^𝑛+1𝑓_𝑛+2. (10)

Let ℎ𝜆=ℎ. Applying the scalar test equation, 𝑦^′=𝜆𝑦𝜆<0, (11) to (10) gives 1−127ℎ27−56414+19136191ℎ1−96191ℎ𝑦_𝑛+1𝑦_𝑛+2=−97187−914ℎ146−191342𝑦191^𝑛−1𝑦_𝑛+011560−27𝑦191^𝑛−3𝑦_𝑛−2. (12)

Equation (12) is equivalent to 𝐴𝑌_𝑚=𝐵𝑌_𝑚−1+𝐶𝑌_𝑚−2,𝑛=2𝑚, (13) where 𝐴=1−127ℎ27−56414+19136191ℎ1−96191ℎ,−9𝐵=7187−914ℎ146−191342,0191𝐶=11560−27191,𝑌_𝑚=𝑦_𝑛+1𝑦_𝑛+2=𝑦_2𝑚+1𝑦_2𝑚+2,𝑌_𝑚−1=𝑦_𝑛−1𝑦_𝑛=𝑦_2𝑚−1𝑦_2𝑚=𝑦_{2(𝑚−1)+1}𝑦_{2(𝑚−1)+2},𝑌_𝑚−2=𝑦_𝑛−3𝑦_𝑛−2=𝑦_{2(𝑚−2)+1}𝑦_{2(𝑚−2)+2}. (14) The stability polynomial 𝑅(𝑡,ℎ) is obtained by evaluating the determinant of (𝐴𝑡²−𝐵𝑡−𝐶) to obtain 𝑅𝑡,ℎ=169𝑡+5348279𝑡²+5348765ℎ𝑡²−267411385𝑡³−5348351ℎ𝑡³−191162ℎ²𝑡³[phone]𝑡⁴−53486171ℎ𝑡⁴[phone]ℎ²𝑡⁴1337=0. (15) To show that the method is zero stable, we let ℎ=0 in (15) and the following equation for the first characteristics polynomial is obtained as 169𝑡+5348279𝑡²−534811385𝑡³[phone]𝑡⁴5348=0. (16) Solving (16) for 𝑡, we have 𝑡=0,𝑡=1,𝑡=−0.105502,𝑡=0.146463. (17) The values of 𝑡 above indicate that the method is zero stable since no magnitude of the root is greater than one and the root 𝑡=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 𝐴-stable since the region covers the entire negative half plane.

Figure 1 

Stability region of the 2-point VSSBBDF when 𝑟=1.

The same procedure is applied to obtain the stability polynomial of the method (7) (corresponding to 𝑟=2) as 𝑅𝑡,ℎ=169𝑡+4422411099𝑡²[phone]ℎ𝑡²−1105672625𝑡³−110565727ℎ𝑡³−2764405ℎ²𝑡³+6914363𝑡⁴−69111943ℎ𝑡⁴[phone]ℎ²𝑡⁴691=0. (18) When ℎ=0 in (18), the first characteristic polynomial of the method (7) becomes 169𝑡+4422411099𝑡²−4422472625𝑡³[phone]𝑡⁴691=0 (19) and solving (19) for 𝑡 gives the following roots: 𝑡=0,𝑡=1,𝑡=−0.0116404,𝑡=0.051994. (20)

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 𝐴-stable.

Figure 2 

Stability region of the 2-point VSSBBDF when 𝑟=2.

Similarly, the stability polynomial for the method (8) (corresponding to 𝑟=10/17) is given by 𝑅𝑡,ℎ=69347235737𝑡−302722616000128013466997𝑡²+3027226160006400808703ℎ𝑡²−1513613080021069763237𝑡³−151361308001629222471ℎ𝑡³−5405761002670327ℎ²𝑡³+3784032760007687𝑡⁴−37840327119480103ℎ𝑡⁴+7568065418988992ℎ²𝑡⁴37840327=0 (21) and when ℎ=0 in (21), the first characteristics polynomial is given by 69347235737𝑡−302722616000128013466997𝑡²−30272261600021069763237𝑡³+1513613080060007687𝑡⁴37840327=0. (22) When (22) is solved for 𝑡, the following values are obtained: 𝑡=0,𝑡=1,𝑡=−0.446055,𝑡=0.32385. (23) Hence, the method (8) is also zero stable.

The stability plot for (21) is given in Figure 3. The method (8) is therefore almost 𝐴-stable.

Figure 3 

Stability region of the 2-point VSSBBDF when 𝑟=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: 𝑦_𝑛+1=𝜃₁𝑦_𝑛+2+𝛼₁ℎ𝑓_𝑛+𝛼₂ℎ𝑓_𝑛+1+Ψ₁,𝑦_𝑛+2=𝜃₂𝑦_𝑛+1+𝛼₃ℎ𝑓_𝑛+1+𝛼₄ℎ𝑓_𝑛+2+Ψ₂, (24) which is equivalent to the following matrix equation: 10𝑦01𝑦^𝑛+1_𝑛+2=0𝜃₁𝜃₂0𝑦_𝑛+1𝑦_𝑛+2+ℎ0𝛼₁𝑓00𝑓^𝑛−1_𝑛𝛼+ℎ²0𝛼₃𝛼₄𝑓_𝑛+1𝑓_𝑛+2+Ψ₁Ψ₂. (25) Equation (25) can be represented as (𝐼−𝐴)𝑌=ℎ𝐵₁𝐹₁+𝐵₂𝐹₂+𝜁, (26) where 𝐼=1001,𝐴=0𝜃₁𝜃₂0,𝐵₁=0𝛼₁,𝐵00₂=𝛼₂0𝛼₃𝛼₄𝑦,𝑌=^𝑛+1𝑦_𝑛+2,𝐹₁=𝑓_𝑛−1𝑓_𝑛,𝐹₂=𝑓_𝑛+1𝑓_𝑛+2Ψ,𝜁=¹Ψ₂. (27)

Let 𝐹=(𝐼−𝐴)𝑌−ℎ𝐵₁𝐹₁+𝐵₂𝐹₂−𝜁=0. (28) Newton’s iteration takes the form 𝑌^(𝑖+1)_{𝑛+1,𝑛+2}−𝑌^(𝑖)_{𝑛+1,𝑛+2}=−𝐹^′_𝑗𝑦^(𝑖)_{𝑛+1,𝑛+2}⁻¹𝐹_𝑗𝑦^(𝑖)_{𝑛+1,𝑛+2}. (29) Equation (32) is equivalent to 𝑌^(𝑖+1)_{𝑛+1,𝑛+2}−𝑌^(𝑖)_{𝑛+1,𝑛+2}=−(𝐼−𝐴)−ℎ𝐵₁𝛿𝐹₁𝑌𝛿𝑌^{(𝑖)𝑛+1,𝑛+2}−ℎ𝐵₂𝛿𝐹₂𝑌𝛿𝑌^{(𝑖)𝑛+1,𝑛+2−1}×(𝐼−𝐴)𝑌^(𝑖)_{𝑛+1,𝑛+2}−ℎ𝐵₁𝐹₁−ℎ𝐵₂𝐹₂−𝜁. (30) We show the accuracy of the method by finding the maximum error in the following way. Let 𝑦_𝑖 and 𝑦(𝑥_𝑖) be the computed and the exact solutions of (1), respectively. The absolute error is defined by error_𝑖𝑡=𝑦_𝑖𝑡−𝑦𝑥_𝑖𝑡. (31)

The maximum error is given by: MAXE=max_{1≤𝑡≤𝑇}max_{1≤𝑖≤𝑁}error_𝑖𝑡, (32) where 𝑇 is the total number of steps and 𝑁 is the number of equations.

Let 𝑌^(𝑖+1)_𝑛+1 denote the (𝑖+1)th iterate and 𝐸^(𝑖+1)_1,2=𝑌^(𝑖+1)_{𝑛+1,𝑛+2}−𝑌^(𝑖)_{𝑛+1,𝑛+2}. (33) Then (33) can be written as 𝐸^(𝑖+1)_1,2=𝐴⁻¹𝐵 (34) which is equivalent to 𝐴𝐸^(𝑖+1)_1,2=𝐵, (35) where 𝐴=(𝐼−𝐴)−ℎ𝐵₁𝛿𝐹₁𝑌𝛿𝑌^{(𝑖)𝑛+1,𝑛+2}−ℎ𝐵₂𝛿𝐹₂𝑌𝛿𝑌^{(𝑖)𝑛+1,𝑛+2−1},𝐵=−(𝐼−𝐴)𝑌^(𝑖)_{𝑛+1,𝑛+2}−ℎ𝐵₁𝐹₁−ℎ𝐵₂𝐹₂−𝜁. (36) LU decomposition method is used to solve the system (35). For the different values of 𝑟, 𝐴=1−𝛼₂ℎ𝛿𝑓_𝑛+1𝛿𝑦_𝑛+1−𝜃₁−𝜃₂−𝛼₃ℎ𝛿𝑓_𝑛+1𝛿𝑦_𝑛+11−𝛼₄𝑓ℎ𝛿^𝑛+2𝛿𝑦_𝑛+2,𝐵=−𝑦^𝑖_𝑛+1+𝜃₁𝑦^𝑖_𝑛+2+𝛼₁ℎ𝑓^𝑖_𝑛+𝛼₂ℎ𝑓^𝑖_𝑛+1+Ψ−𝑦^𝑖_𝑛+2+𝜃₂𝑦^𝑖_𝑛+1+𝛼₃ℎ𝑓^𝑖_𝑛+1+𝛼₄ℎ𝑓^𝑖_𝑛+2.+Ψ (37) When 𝑟=1𝐴=1−127ℎ𝛿𝑓_𝑛+1𝛿𝑦_𝑛+127−56414+19136ℎ191𝛿𝑓_𝑛+1𝛿𝑦_𝑛+11−96𝑓191ℎ𝛿^𝑛+2𝛿𝑦_𝑛+2,𝐵=−𝑦^𝑖_𝑛+1−27𝑦56^𝑖𝑛+2−914ℎ𝑓^𝑖_𝑛+127ℎ𝑓^𝑖_𝑛+1+11𝑦56^𝑖𝑛−2−97𝑦_𝑛−1+187𝑦_𝑛−𝑦^𝑖_𝑛+2+414𝑦191^𝑖𝑛+1−36191ℎ𝑓^𝑖_𝑛+1+96191ℎ𝑓^𝑖_𝑛+2−27𝑦191^𝑛−2+146𝑦191^𝑛−1−342𝑦191^𝑛. (38) When 𝑟=2𝐴=1−152ℎ𝛿𝑓_𝑛+1𝛿𝑦_𝑛+145−166528+3455144ℎ691𝛿𝑓_𝑛+1𝛿𝑦_𝑛+11−384𝑓691ℎ𝛿^𝑛+2𝛿𝑦_𝑛+2,𝐵=−𝑦^𝑖_𝑛+1−45𝑦16^𝑖𝑛+2−4516ℎ𝑓^𝑖_𝑛+152ℎ𝑓^𝑖_𝑛+1+9𝑦64^𝑛−2−54𝑦_𝑛−1+315𝑦64^𝑛−𝑦^𝑖_𝑛+2+6528𝑦3455^𝑖𝑛+1−144691ℎ𝑓^𝑖_𝑛+1+384691ℎ𝑓^𝑖_𝑛+2−73𝑦3455^𝑛−2+111𝑦691^𝑛−1−711𝑦691^𝑛 (39) and when 𝑟=10/17𝐴=^{1−999ℎ938𝛿𝑓𝑛+1𝛿𝑦𝑛+119869}₋⁸²⁵⁴⁴_{7265280+298527114256ℎ80683𝛿𝑓𝑛+1𝛿𝑦𝑛+11−38016𝑓80683ℎ𝛿𝑛+2𝛿𝑦𝑛+2},𝐵=^{−𝑦𝑖𝑛+1−19869𝑦82544𝑖𝑛+2−29977504ℎ𝑓𝑖𝑛+999938ℎ𝑓𝑖𝑛+1+338997𝑦750400𝑛−2−2363153𝑦1031800𝑛−1+2310687𝑦750400𝑛}_{−𝑦𝑖𝑛+2+7265280𝑦2985271𝑖𝑛+1−1425680683ℎ𝑓𝑖𝑛+1+3801680683ℎ𝑓𝑖𝑛+2−8484751𝑦14926355𝑛−2+191607𝑦80683𝑛−1−1307097𝑦403415𝑛}. (40)

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: LTE=𝑦^𝑘+1_𝑛+2−𝑦^𝑘_𝑛+2,𝑘=3, (41) where 𝑦^𝑘+1_𝑛+2 refers to the formula of order 𝑘+1 and 𝑦^𝑘_𝑛+2 refers to the formula of order 𝑘.

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 𝑟=1) and a new step size is computed as _newℎ=𝑐×_oldℎ×TOLLTE^1/𝑝, (42) where 𝑝 is the order of the method and 𝑐 is a safety factor. Again we choose a safety factor (𝑐=0.3) different from that in [15].

If _newℎ>(1.7×_oldℎ), then the step size ℎ becomes ℎ=1.7×_oldℎ (equivalent to the formula when 𝑟=10/17). Otherwise if LTE>TOL, the step size is halved (equivalent to the formula when 𝑟=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: 𝑦^′=−20𝑦+24,𝑦(0)=0,0≤𝑥≤10. (43) The exact solution is as follows: 𝑦(𝑥)=65−65𝑒^−20𝑥. (44)(2) The cosine problem taken from [8, 15, 19] is as follows: 𝑦^′()−1=−2𝜋sin2𝜋𝑥𝜖(),𝑦(0)𝑦−cos2𝜋𝑥=1,0≤𝑥≤10. (45) The problem becomes increasingly stiff as 𝜖→0. We choose 𝜖=10⁻³. The exact solution is as follows: 𝑦(𝑥)=cos(2𝜋𝑥). (46)(3) We consider the system of nonlinear stiff IVP taken from [18, 20] as follows: 𝑦^′₁=−𝜖⁻¹+2𝑦₁+𝜖⁻¹𝑦²₂,𝑦₁(0)𝑦=1,0≤𝑥≤20_′2=𝑦₁−𝑦₂1+𝑦₂,𝑦₂(0)=1. (47) The solution, which is independent of 𝜖, is 𝑦₁(𝑥)=𝑒^−2𝑥,𝑦₂(𝑥)=𝑒^−𝑥. (48) The eigenvalues are 𝜆₁=−1 and 𝜆₂=−(𝜖⁻¹+2). This allows the degree of stiffness to be regulated. We choose 𝜖=10⁻⁵. (4) This physical problem is taken from [21] as follows: 𝑦^′₁=−21𝑦₁+19𝑦₂−20𝑦₃,𝑦₁(0)𝑦=1,0≤𝑥≤10,_′2=19𝑦₁−21𝑦₂+20𝑦₃,𝑦₂(0)𝑦=0,_′3=40𝑦₁−40𝑦₂−40𝑦₃,𝑦₃(0)=−1. (49) The exact solution is as follows: 𝑦₁(𝑥)=0.5𝑒^−2𝑥+𝑒^−40𝑥(,𝑦cos40𝑥+sin40𝑥)₂(𝑥)=0.5𝑒^−2𝑥−𝑒^−40𝑥(,𝑦cos40𝑥+sin40𝑥)₃(𝑥)=2𝑒^−40𝑥−121cos40𝑥+2sin40𝑥. (50) The eigenvalues are −2, −40+40𝑖, and −40−40𝑖.

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⁻², 10⁻⁴, and 10⁻⁶. 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 Tables 1, 2, 3, and 4: TOL: Tolerance value used Method: Name of the method IST: Number of successful steps FS: Number of failed steps TS: Total steps MAXE: Maximum error.

Table 1 

Numerical results for problem (1) (the circuit problem).

Table 2 

Numerical results for problem (2) (the cosine problem).

Table 3 

Numerical results for the nonlinear problem (3).

Table 4 

Numerical results for problem (4).

From Tables 1, 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 of Log10MAXE versus TOL for problem (1).

Figure 5 

Graph of Log10MAXE versus TOL for problem (2).

Figure 6 

Graph of Log10MAXE versus TOL for problem (3).

Figure 7 

Graph of Log10MAXE 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, 𝐴-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 𝐴-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.

Acknowledgments

The authors thank the anonymous reviewers for their comments which improved the paper and the Institute for Mathematical Research, Universiti Putra Malaysia, for funding the research.