Sage Journals: Discover world-class research

Abstract

Stability region of two-point variable step–block backward differentiation formulae method is analysed in this article where we took into consideration on the increment of step size to a factor 1.8. The stability graph for three distinct step size ratios are plotted using Maple software and presented in different graphs. The stability properties are also become part of the analysis that have been studied in this article.

Keywords

Backward differentiation formulae block backward differentiation formulae stiff ordinary differential equations stability region

Introduction

Apparently, the study on numerical methods for the solution of initial value problem (IVP) of Ordinary Differential Equations (ODEs) especially stiff problem has become a famous topic and advanced in study as many researchers have produced new findings in this particular research areas. The existing numerical methods used for solving ODE are designed to find the approximations to the solutions of ODE where it provides an alternative way to some complex problems in real-world nowadays. Apart from that, the numerical methods are classified as single-step methods or multi-step methods. In addition, the methods can be divided into explicit methods or implicit methods.

Since numerical method has stability limitation on the step size, there are only few methods that can solve stiff problems.¹ The ideas of solving first-order stiff problem using block backward differentiation formulae (BBDF) method has been introduced by Ibrahim et al.² Currently, BBDF method becomes frequently used for solving stiff ODEs. Furthermore, there are variety of solver, which are based on BBDF method that are available to solve stiff ODE developed in literature.^3–10 Basically, for method to be of practical importance, it must have a region of absolute stability to ensure that the method will be able to solve at least for the mildly stiff problems.¹¹ In this article, we are interested to analyse the stability of VS-BBDF method with constant step size, half the step size and increment the step size by a factor 1.8. Besides, the variable step–block backward differentiation formula (VS-BBDF) method with increment of step size to a factor 1.6 and 1.9 have been studied by Ibrahim et al.² and Yatim et al.⁷

The formulae for VS-BBDF method of variable step size approach with increment of the step size to a factor 1.8 is as follows (Table 1).

Table 1.

The formulae for two-point VS-BBDF method.

Step size ratio	Computed points	Coefficients of the points
r = 1	$y n + 1$	$\frac{6}{5} hf n + 1 - \frac{3}{10} y n + 2 + \frac{1}{10} y n - 2 - \frac{3}{5} y n - 1 + \frac{9}{5} y n$
	$y n + 2$	$\frac{12}{25} hf n + 2 + \frac{48}{25} y n + 1 - \frac{3}{25} y n - 2 + \frac{16}{25} y n - 1 - \frac{36}{25} y n$
r = 2	$y n + 1$	$\frac{15}{8} hf n + 1 - \frac{75}{128} y n + 2 + \frac{3}{128} y n - 2 - \frac{25}{128} y n - 1 + \frac{225}{128} y n$
	$y n + 2$	$\frac{12}{23} hf n + 2 + \frac{192}{115} y n + 1 - \frac{2}{115} y n - 2 + \frac{3}{23} y n - 1 - \frac{18}{23} y n$
r = 5/9	$y n + 1$	$\frac{266}{297} hf n + 1 - \frac{2527}{13662} y n + 2 + \frac{189}{550} y n - 2 - \frac{9747}{6325} y n - 1 + \frac{17689}{7425} y n$
	$y n + 2$	$\frac{644}{1425} hf n + 2 + \frac{59248}{27075} y n + 1 - \frac{128547}{225625} y n - 2 + \frac{27216}{11875} y n - 1 - \frac{103684}{35625} y n$

The general form of the formulae is

[1001] [y n + 1 y n + 2] = [0 α 1 α 2 0] [y n + 1 y n + 2] + h [β 1 00 β 2] [f n + 1 f n + 2] + [μ 1 τ 1 μ 2 τ 2] [y n - 1 y n] + [0 ψ 1 0 ψ 2] [y n - 3 y n - 2]

(1)

where

ψ 1

and

ψ 2

are the value of previous points.

Stability and its properties

Our goal is to determine the stability properties for VS-BBDF method when the step size is constant, halved and increased to a factor 1.8. There are some definition on stability of linear multistep method (LMM) from Lambert¹² that will be defined later in this section.

The general LMM is

\sum_{j = 0}^{k} α j y n + j = h \sum_{j = 0}^{k} β j f n + j

(2)

Therefore, the general form of a block LMM is

\sum_{j = 0}^{k} A j Y n + j = h \sum_{j = 0}^{k} B j F n + j

(3)

where A_j and B_j are r by r matrices with elements

α im, β im

for

i \cdot m = 0, 1, \dots, r

Then, apply (3) to VS-BBDF method and yield the general form of VS-BBDF method as below

\sum_{j = 0}^{k} α ij y n + 2 - j = h \sum_{j = 0}^{k} β ij f n + 2 - j

(4)

or the general form of VS-BBDF method in matrix-vector as shown in (1).

Definition 2.1

The LMM (2) is said to be absolute stable in a region R for a given $h λ$ , all the roots r_s of the stability polynomial $π (r; h λ) : = ρ (r) - h λ σ (r) = 0$ , satisfy $| r s | < 1$ , where $s = 1, 2, \dots, k$ .

Definition 2.2

The LMM (2) is said to be A-stable if its region of absolute stability contains the whole of the left-hand half-plane $Re (h λ) < 0$ .

Definition 2.3

The LMM (2) is said to be zero stable of no root of the first characteristic polynomial $p (r)$ has modulus greater than one, and if every root with unit modulus is simple.

Stability analysis of VS-BBDF method

Apply (1) to the test equation, $f = y' = λ y$ and then obtain

[1001] [y n + 1 y n + 2] = [0 α 1 α 2 0] [y n + 1 y n + 2] + h [β 1 00 β 2] [λ y n + 1 λ y n + 2] + [μ 1 τ 1 μ 2 τ 2] [y n - 1 y n] + [0 ψ 1 0 ψ 2] [y n - 3 y n - 2]

(5)

or let

H = h λ

and yield

[1 - β 1 H - α 1 - α 2 1 - β 2 H] [y n + 1 y n + 2] = [μ 1 τ 1 μ 2 τ 2] [y n - 1 y n] + [0 ψ 1 0 ψ 2] [y n - 3 y n - 2]

(6)

In simpler way,

A 0 Y j = A 1 Y j - 1 + A 2 Y j - 2

(7)

where n = 2j.

The stability polynomial of the method is

R (t; H) = det (\sum_{j = 0}^{k} A j t^{k - j}) = 0

(8)

Hence (8) is equivalent to

R (t; H) = det | A 0 t^{2} - A 1 t - A 2 | = 0

(9)

For r = 1,

A 0 = [1 - \frac{6}{5} H \frac{3}{10} - \frac{48}{25} 1 - \frac{12}{25} H], A 1 = [- \frac{3}{5} \frac{9}{5} \frac{16}{25} - \frac{36}{25}], A 2 = [0 \frac{1}{10} 0 - \frac{3}{25}]

For r = 2,

A 0 = [1 - \frac{15}{8} H \frac{75}{128} - \frac{192}{115} 1 - \frac{12}{23} H], A 1 = [- \frac{25}{128} \frac{225}{128} \frac{3}{23} - \frac{18}{23}], A 2 = [0 \frac{3}{280} 0 - \frac{2}{115}]

For r = 5/9,

A 0 = [1 - \frac{266}{297} H \frac{2527}{13662} - \frac{59248}{27075} 1 - \frac{644}{1425} H], A 1 = [- \frac{9747}{6325} \frac{17689}{7425} \frac{27216}{11875} - \frac{103684}{35625}], A 2 = [0 \frac{189}{550} 0 - \frac{128547}{225625}]

Stability polynomial of all step size ratios are as follows,

r = 1,

R (t; H) = \frac{197}{125} t^{4} - \frac{153}{125} t^{3} - \frac{9}{25} t^{2} + \frac{1}{125} t - \frac{42}{25} t^{4} H + \frac{72}{125} t^{4} H^{2} - \frac{252}{125} t^{3} H - \frac{18}{125} t^{2} H

(10)

r = 2,

R (t; H) = \frac{91}{46} t^{4} - \frac{173}{93} t^{3} - \frac{289}{2944} t^{2} + \frac{1}{2944} t - \frac{441}{184} t^{4} H + \frac{45}{46} t^{4} H^{2} - \frac{1155}{736} t^{3} H - \frac{3}{92} t^{2} H

(11)

r = 5/9,

R (t; H) = \frac{31291}{22275} t^{4} - \frac{3575389}{10580625} t^{3} - \frac{755829}{653125} t^{2} + \frac{59049}{653125} t - \frac{190106}{141075} t^{4} H + \frac{9016}{22275} t^{4} H^{2} - \frac{1839404}{556875} t^{3} H - \frac{66654}{130625} t^{2} H

(12)

The absolute stability region for each step size ratios in the $h λ$ plane is investigated by evaluating the stability polynomial. The stability region is the region enclosed by the set of points determined by replacing $t = e^{i θ} = \sin θ + i \cos θ, 0 < θ < 2 π$ in stability polynomial. By the definition mentioned earlier, the stability region is determined by finding the region for which $| t | < 1$ . Furthermore, the region of stability of VS-BBDF method is plotted using MAPLE software.

To obtain zero stability, let $H = h λ = 0$ and replace H into all stability polynomials in (10), (11), and (12). Thus, we yield

r = 1,

R (t; 0) = \frac{197}{125} t^{4} - \frac{153}{125} t^{3} - \frac{9}{25} t^{2} + \frac{1}{125} t

(13)

r = 2,

R (t; 0) = \frac{91}{46} t^{4} - \frac{173}{93} t^{3} - \frac{289}{2944} t^{2} + \frac{1}{2944} t

(14)

r = 5/9,

Solve the updated equations will yield the roots for the stability polynomial. Thus, the stability polynomial and roots for all distinct step size ratios of VS-BBDF method are listed in Table 2.

Table 2.

Lists of stability polynomials and roots for the three step size ratios.

Step size ratio, r	Stability polynomial, $R (t; H)$	Roots
r = 1	$\frac{197}{125} t^{4} - \frac{153}{125} t^{3} - \frac{9}{25} t^{2} + \frac{1}{125} t - \frac{42}{25} t^{4} H + \frac{72}{125} t^{4} H^{2} - \frac{252}{125} t^{3} H - \frac{18}{125} t^{2} H$	1, 0.0207917599, −0.2441420137
r = 2	$\frac{91}{46} t^{4} - \frac{173}{93} t^{3} - \frac{289}{2944} t^{2} + \frac{1}{2944} t - \frac{441}{184} t^{4} H + \frac{45}{46} t^{4} H^{2} - \frac{1155}{736} t^{3} H - \frac{3}{92} t^{2} H$	1, 0.00325762197, −0.05270817143
r = 5/9	$\frac{31291}{22275} t^{4} - \frac{3575389}{10580625} t^{3} - \frac{755829}{653125} t^{2} + \frac{59049}{653125} t - \frac{190106}{141075} t^{4} H + \frac{9016}{22275} t^{4} H^{2} - \frac{1839404}{556875} t^{3} H - \frac{66654}{130625} t^{2} H$	1, 0.0769489074 −0.8363962010

Numerical examples

The test problems of stiff ODEs are solved using the VS-BBDF method and we compared the numerical results with ode15s and ode23s in term of maximum errors and number of total steps.

The notations are used in the tables and figures.

VS-BBDF	Variable step–block backward differentiation formulae method
TSs	Total steps taken during the computation of approximate solution
TOL	Tolerance limit
MAXE	Maximum error

Test Problem 1:

$y' = - 1000 y + 3000 - 2000 e^{- x}$

Interval: $[0, 20]$

Exact solution: $y (x) = 3 - 0.998 e^{- 1000 x} - 2.002 e^{- x}$

Initial conditions: $y (0) = 0$

Eigenvalue: $λ = - 100$

Source: Voss and Abbas¹³

Test Problem 2:

$y' 1 (x) = 998 y 1 + 1998 y 2$

$y' 2 (x) = - 999 y 1 - 1999 y 2$

Interval: $[0, 20]$

Exact solution: $y 1 (x) = 2 e^{- x} - e^{- 1000 x} y 2 (x) = - e^{- x} + e^{- 1000 x}$

Initial conditions: $y 1 (0) = 1, y 2 (0) = 0$

Eigenvalues: $λ = - 1, - 1000$

Source: Gear⁴

Test Problem 3:

$y' 1 = 1195 y 1 - 1995 y 2 y' 2 = 1197 y 1 - 1997 y 2$

Interval: $[0, 20]$

Exact solution: $y 1 (x) = 10 e^{- 2 x} - 8 e^{- 800 x} y 2 (x) = 6 e^{- 2 x} - 8 e^{- 800 x}$

Initial conditions: $y 1 (0) = 2, y 2 (0) = - 2$

Eigenvalues: $λ = - 2, - 800$

Source: Gerald and Wheatley¹⁴

Results for Test Problem 1

Table 3.

Numerical results for Test Problem 1.

Problem	Method	TOL	TSs	MAXE
'1	VS-BBDF	10⁻²	29	1.1090e-004
		10⁻⁴	56	1.5807e-006
		10⁻⁶	135	9.9454e-007
	ode15s	10⁻²	40	8.4000e-003
		10⁻⁴	93	1.6634e-004
		10⁻⁶	165	3.0953e-006
	ode23s	10⁻²	36	4.5000e-003
		10⁻⁴	181	2.5500e-004
		10⁻⁶	1193	1.0911e-005

Results for Test Problem 2

Table 4.

Numerical results for Test Problem 2.

Problem	Method	TOL	TSs	MAXE
2	VS-BBDF	10⁻²	30	1.1291e-004
		10⁻⁴	61	1.0537e-007
		10⁻⁶	152	1.1345e-008
	ode15s	10⁻²	37	1.7600e-002
		10⁻⁴	89	1.8659e-004
		10⁻⁶	167	3.9569e-006
	ode23s	10⁻²	22	7.3100e-003
		10⁻⁴	67	3.6837e-004
		10⁻⁶	287	1.7039e-005

Results for Test Problem 3

Table 5.

Numerical results for Test Problem 3.

Problem	Method	TOL	TSs	MAXE
3	VS-BBDF	10⁻²	36	1.1792e-004
		10⁻⁴	76	9.2833e-007
		10⁻⁶	196	1.0238e-008
	ode15s	10⁻²	52	1.0940e-001
		10⁻⁴	67	1.0940e-001
		10⁻⁶	85	1.0940e-001
	ode23s	10⁻²	30	4.4400e-002
		10⁻⁴	37	4.4400e-002
		10⁻⁶	43	4.4400e-002

Conclusions

From Figures 1 –3, the stability regions for all the step size ratios lie outside the closed region. Since all the roots for the step size ratios have modulus less than or equal to 1, thus the proposed method satisfy the zero stability condition. The figures also show that the absolute stability region for step size ratios 1 and 5/9 are almost A-stable while the stability region for step size ratio 2 is A-stable since the absolute stability region covers the entire left half-plane of the complex plane, $Re (h λ) < 0$ . Tables 3 –5 present the numerical results for the solver in terms of maximum errors. From Figures 4 to 9, it is proven that the maximum error obtained are within the tolerance given and VS-BBDF method converges faster as compared with ode15s and ode23s. Therefore, we can conclude that the proposed method is suitable for solving the stiff ODE problems.

Figure 1.

Stability region of VS-BBDF method when r = 1.

Figure 2.

Stability region of VS-BBDF method when r = 2.

Figure 3.

Stability region of VS-BBDF method when r = 5/9.

Figure 4.

Efficiency curves VS-BBDF and Matlab’s ODE solvers for Test Problem 1.

Figure 5.

Total steps curves of VS-BBDF and Matlab’s ODE solvers for Test Problem 1.

Figure 6.

Efficiency curves VS-BBDF and Matlab’s ODE solvers for Test Problem 2.

Figure 7.

Total steps curves of VS-BBDF and Matlab’s ODE solvers for Test Problem 2.

Figure 8.

Efficiency curves VS-BBDF and Matlab’s ODE solvers for Test Problem 3.

Figure 9.

Total steps curves of VS-BBDF and Matlab’s ODE solvers for Test Problem 3.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This paper is supported by Universiti Sains Malaysia under Short Term Grant USM, 304/PJJAUH/6313177.