Diagonally Implicit Block Backward Differentiation Formulas for Solving Ordinary Differential Equations

Abstract

This paper focuses on the derivation of diagonally implicit two-point block backward differentiation formulas (DI2BBDF) for solving first-order initial value problem (IVP) with two fixed points. The method approximates the solution at two points simultaneously. The implementation and the stability of the proposed method are also discussed. A performance of the DI2BBDF is compared with the existing methods.

1. Introduction

Most real life problems that arise in various fields of study such as engineering or science are constructed as mathematical models before they are solved. These models often lead to differential equations. A differential equation can simply be defined as an equation that contains a derivative. In this paper, we are concerned with the numerical solution of IVP with two fixed points of the general form𝑦^′=𝑦−𝑣₁𝑦−𝑣₂𝑔𝑦, (1.1) given initial values 𝑦(𝑥_𝑛)=𝑦_𝑛, where 𝑣₁<𝑣₂∈ℝ and 𝑔(𝑦)≠𝟎 is a bounded real-valued function with continuous derivatives. See; [1] how the problem may be simplified and we can assume that the fixed points are 𝑦₁(𝑥)=𝟎 and 𝑦₂(𝑥)=𝟏. The aim of this paper is to introduce a new block method which is called diagonally implicit two point block backward differentiation formulas. Basically, block method is a method which is used to compute previous 𝑘 blocks and to calculate the current block where each block contain 𝑟 points. Ibrahim et al. [2] gave the general form for 𝑟 point 𝑘 block as ^𝑘_𝑗=0𝐴_𝑗𝑦_𝑛+𝑗=ℎ^𝑘_𝑗=0𝐵_𝑗𝑓_𝑛+𝑗, (1.2) where 𝐴_𝑗 and 𝐵_𝑗 are 𝑟 by 𝑟 matrices.

The rapid growth of the studies on the block methods for solving ODEs contribute to the competition in developing and deriving an accurate method for solving many types of ODEs. Among the earliest research on block methods was proposed by Shampine and Watts [3] with block implicit one step method, Cash [4] with modified extended backward differential formulas, Chu and Hamilton [5] with multi-block methods, Majid and Suleiman [6] with fully implicit block method for solving ODEs. The most recent study of block method is on computing block approximation 𝑦_𝑛+1,𝑦_𝑛+2…𝑦_𝑛+𝑘 known as block backward differentiation formulas (BBDF) presented by Ibrahim et al. [7]. Diagonally implicit for multistep method are discussed by a few researchers such as Majid and Suleiman [8] with 4-point diagonally implicit block method and Alexander [9] with diagonally implicit Runge-Kutta method for solving stiff ODEs. In this paper, we interested to compare the accuracy of the proposed method with the one-step predictor-corrector method (1SPCM) and P-C ADAMS method studied by Aguiar and Ramos [1].

2. Formulation of Diagonally Implicit BBDF

In Figure 1, diagonally implicit 2-point BBDF (DI2BBDF) will create two new equally spaced solution values simultaneously. Majid and Suleiman [8] stated that the method is called diagonally implicit because the coefficients of the upper triangular matrix entries are zero.

Figure 1 

Diagonally Implicit 2-point block method of constant step size.

We will derive the formula to compute the approximate values of 𝑦_𝑛+1 and 𝑦_𝑛+2 simultaneously with two previous back values 𝑥_𝑛−1 and 𝑥_𝑛. The formula was derived by using Lagrange interpolation polynomial which compute the solutions at 𝑥_𝑛+1 and 𝑥_𝑛+2 separately. Consider the polynomial 𝑃_𝑘(𝑥) of degree 𝑘 which interpolates the values𝑦_𝑛, 𝑦_𝑛−1,⋯, 𝑦_{𝑛−𝑘+1} of a function 𝑓 at interpolating points 𝑥_𝑛, 𝑥_𝑛−1,⋯, 𝑥_{𝑛−𝑘+1} in terms of Lagrange polynomial defined as follows: 𝑃_𝑘(𝑥)=^𝑘_𝑗=0𝐿_𝑘,𝑗(𝑥)𝑓𝑥_{𝑛+1−𝑗}, (2.1) where 𝐿_𝑘,𝑗(𝑥)=^𝑘_{𝑖=0𝑖≠𝑗}𝑥−𝑥_{𝑛+1−𝑖}𝑥_{𝑛+1−𝑗}−𝑥_{𝑛+1−𝑖} (2.2) for each 𝑗=0,1,…,𝑘.

For 𝑦_𝑛+1, we defined 𝑠=(𝑥−𝑥_𝑛+1)/ℎ and replace 𝑓(𝑥,𝑦) by polynomial (2.1) which interpolates only the values 𝑦_𝑛−1, 𝑦_𝑛, and 𝑦_𝑛+1 at the interpolating points 𝑥_𝑛−1,𝑥_𝑛,…,𝑥_𝑛+1. 𝑃𝑥_𝑛+1=𝑥+𝑠ℎ^𝑛+1+𝑠ℎ−𝑥_𝑛𝑥_𝑛+1+𝑠ℎ−𝑥_𝑛+1()𝑦−ℎ)(−2ℎ^𝑛−1+𝑥_𝑛+1+𝑠ℎ−𝑥_𝑛−1𝑥_𝑛+1+𝑠ℎ−𝑥_𝑛+1(ℎ)𝑦)(−ℎ^𝑛+𝑥_𝑛+1+𝑠ℎ−𝑥_𝑛−1𝑥_𝑛+1+𝑠ℎ−𝑥_𝑛(ℎ)𝑦2ℎ)(^𝑛+1. (2.3) Thus, differentiating the resulting polynomial once with respect to 𝑠 at the point 𝑥=𝑥_𝑛+1 and evaluating at 𝑠=0 gives the following: 𝑃^′(𝑥)=𝑃^′𝑥_𝑛+1=ℎ𝑓_𝑛+1=12𝑦_𝑛−1−2𝑦_𝑛+32𝑦_𝑛+1. (2.4)

Then, we interpolate the values 𝑦_𝑛−1, 𝑦_𝑛, 𝑦_𝑛+1, and 𝑦_𝑛+2 at the interpolating points 𝑥_𝑛−1, 𝑥_𝑛, 𝑥_𝑛+1,…,𝑥_𝑛+2 and differentiating the resulting polynomial once with respect to 𝑠 at the point 𝑥=𝑥_𝑛+2 and substituting 𝑠=0 yield 𝑃^′(𝑥)=𝑃^′𝑥_𝑛+2=ℎ𝑓_𝑛+21=−3𝑦_𝑛−1+32𝑦_𝑛−3𝑦_𝑛+1+116𝑦_𝑛+2. (2.5)

Therefore, the corrector formulae at 𝑥_𝑛+1 and 𝑥_𝑛+2 are given by 𝑦_𝑛+11=−3𝑦_𝑛−1+43𝑦_𝑛+23ℎ𝑓_𝑛+1,𝑦_𝑛+2=2𝑦11^𝑛−1−9𝑦11^𝑛+18𝑦11^𝑛+1+611ℎ𝑓_𝑛+2. (2.6)

3. Stability of Diagonally Implicit BBDF

Many numerical methods depend on speed of convergence, computational expense and accuracy for the solution of initial value problems (IVPs). In this section, our aim is to identify the properties of linear stability of DI2BBDF. According to Shampine and Watts [3], the most critical limitation of block method is due to the existence of stability problems. The basic definition of zero stable for block method described by Ibrahim et al. [10]. The block method is said to be zero-stable if the roots 𝑅_𝑗,𝑗=1(1)𝑘, of the first characteristic polynomial 𝜌(𝑅)=det[∑^𝑘_𝑖=0𝐴_𝑖𝑅^𝑘−𝑖]=0,𝐴₀=−𝐼, satisfy |𝑅_𝑗|≤1. If one of the roots is +1, we call this root the principal root of 𝜌(𝑅).

The linear stability properties are determined through application of the standard linear test problem 𝑦^′=𝜆𝑦,𝜆<0,𝜆complex. (3.1) Hence, we will construct the formulas that have been derived in the previous section, which is (2.6). Application of the standard linear test problem (3.1) to (2.6) gives the following 𝑦_𝑛+1−23𝜆ℎ𝑦_𝑛+11=−3𝑦_𝑛−1+43𝑦_𝑛.𝑦_𝑛+2−611𝜆ℎ𝑦_𝑛+2=2𝑦11^𝑛−1−9𝑦11^𝑛+18𝑦11^𝑛+1. (3.2) Setting ℎ=𝜆ℎ, we write (3.2) in the matrix coefficient specified as 2𝐴=1−3ℎ0061−ℎ11,𝑌_𝑚=𝑦_𝑛+1𝑦_𝑛+2−1,𝐵=3432−91111,𝑌_𝑚−1=𝑦_𝑛−1𝑦_𝑛. (3.3) The first characteristic polynomial of (3.2) is given by 𝜌(𝑡)=det[𝑡𝐴−𝐵]=𝑡²−40𝑡33²ℎ−34433𝑡+𝑡11²ℎ²−8𝑡111ℎ+3. (3.4)

By solving 𝜌(𝑡)=0 and ℎ=0, we can determine the zero stable and yield 𝑡=1 and 𝑡=1/33. Thus the diagonally implicit BBDF is zero stable. Since one of the roots is +1, we call this root as principal root.

4. Implementation of the Method

In this section, we were obtained the calculation of 𝑦_𝑛+1 and 𝑦_𝑛+2 using the Newton iteration. The DI2BBDF can be written in general form as 𝐹₁=𝑦_𝑛+1−23ℎ𝑓_𝑛+1−𝜇₁𝐹₂=𝑦_𝑛+2−18𝑦11^𝑛+1−611ℎ𝑓_𝑛+2−𝜇₂, (4.1) where 𝜇₁ and 𝜇₂ are the backvalues.

To specify the iteration, the following notation is introduced. 𝑦^(𝑖+1)_𝑛+1 will denote the (𝑖+1)th iterative value of 𝑦_𝑛+1, and let 𝑒^(𝑖+1)_𝑛+1=𝑦^(𝑖+1)_𝑛+1−𝑦^(𝑖)_𝑛+1 and 𝑒^(𝑖+1)_𝑛+2=𝑦^(𝑖+1)_𝑛+2−𝑦^(𝑖)_𝑛+2. Newton iteration takes the form as below: 𝑦^(𝑖+1)_𝑛+1=𝑦^(𝑖)_𝑛+1−𝐹₁𝑦^(𝑖)_𝑛+1𝐹^′₁𝑦^(𝑖)_𝑛+1⁻¹,𝑦^(𝑖+1)_𝑛+2=𝑦^(𝑖)_𝑛+2−𝐹₂𝑦^(𝑖)_𝑛+2𝐹^′₂𝑦^(𝑖)_𝑛+2⁻¹. (4.2)

Substituting𝐹₁𝐹^′₁, 𝐹₂ and 𝐹^′₂ take the form 𝑦^(𝑖+1)_𝑛+1−𝑦^(𝑖)_𝑛+1=−𝑦^(𝑖)_𝑛+1−23ℎ𝑓^(𝑖)_𝑛+1−𝜇₁21−3ℎ𝜕𝑓_𝑛+1𝜕𝑦_𝑛+1⁻¹,21−3ℎ𝜕𝑓_𝑛+1𝜕𝑦_𝑛+1𝑦^(𝑖+1)_𝑛+1−𝑦^(𝑖)_𝑛+1=−𝑦^(𝑖)_𝑛+1+23ℎ𝑓^(𝑖)_𝑛+1+𝜇₁. (4.3)

Similarly for the second point, 𝑦_𝑛+2 gives 𝑦^(𝑖+1)_𝑛+2−𝑦^(𝑖)_𝑛+2=−𝑦^(𝑖)_𝑛+2−18𝑦11^(𝑖)𝑛+1−611ℎ𝑓^(𝑖)_𝑛+2−𝜇₂61−ℎ11𝜕𝑓_𝑛+2𝜕𝑦_𝑛+2⁻¹,𝑦^(𝑖+1)_𝑛+2−𝑦^(𝑖)_𝑛+261−ℎ11𝜕𝑓_𝑛+2𝜕𝑦_𝑛+2=−𝑦^(𝑖)_𝑛+2+18𝑦11^(𝑖)𝑛+1+611ℎ𝑓^(𝑖)_𝑛+2+𝜇₂. (4.4)

Hence, equation (4.3) and (4.4) in matrix form 21−3ℎ𝜕𝑓_𝑛+1𝜕𝑦_𝑛+10186111−ℎ11𝜕𝑓_𝑛+2𝜕𝑦_𝑛+2𝑒^(𝑖+1)_𝑛+1𝑒^(𝑖+1)_𝑛+2=−1018𝑦11−1𝑦^{(𝑖)𝑛+1(𝑖)}_𝑛+22+ℎ3006𝑓11^(𝑖)𝑛+1𝑓^(𝑖)_𝑛+2+𝜇₁𝜇₂. (4.5)

Therefore, the approximation values of 𝑦_𝑛+1and 𝑦_𝑛+2 are 𝑦^(𝑖+1)_𝑛+1=𝑦^(𝑖)_𝑛+1+𝑒^(𝑖+1)_𝑛+1,𝑦^(𝑖+1)_𝑛+2=y^(𝑖)_𝑛+2+𝑒^(𝑖+1)_𝑛+2. (4.6)

5. Numerical Example

In this section, an autonomous problem with two fixed points considered for the purpose of validating the numerical results: 𝑦^′(𝑥)=𝑦(𝑥)𝑦(𝑥)−1𝑦(𝑥).−2 (5.1) The solution is 𝑦(0)=𝑦₀.

Exact solution for 0<𝑦₀<1 is 𝑦(𝑥)=𝑒^𝑥/2𝑒^𝑥/2𝑦²₀−𝑦₀𝑒^𝑥𝑦²₀−4𝑦₀+42𝑦₀.−1 (5.2) See Aguiar and Ramos [1].

We have tested the problem for 𝑦₀=0.1 and consider the step size, ℎ=2^−𝑗 where 𝑗=2,3,4,5,6,7,8 with interval [0,20]. The accuracy for all methods in terms of maximum error is shown in Table 1.

Table 1 

Numerical results for problem (5.1) for different step sizes, ℎ=2^−𝑗.

The notation used in the table takes the following meaning: ℎ: step size, MAX 𝐸: maximum error, DI2BBDF: diagonally implicit 2-point block backward differentiation formulas method, 1SPCM: one-step predictor-corrector method,  P-C ADAMS: explicit Euler as a predictor and trapezoidal rule as a corrector.

The maximum error is defined as follows: Maximum error = max_{1≤𝑖≤NS}(|𝑦_𝑖−𝑦(𝑥)_𝑖|) where NS is the number of steps.

In Figure 2, the numerical results obtained with the proposed method and existing methods are plotted. The exact and approximate solutions using ℎ=0.4 are also plotted in Figure 3.

Figure 2 

The graph of log10 (MAXE) against 𝑗.

Figure 3 

The graph of exact and approximate solutions against 𝑥∈[0,20] using ℎ=0.4.

6. Conclusion

In this work, we have applied DI2BBDF as a numerical method for solving ordinary initial value problem with two fixed points. By reducing the step size ℎ, we observe that the value of maximum errors for DI2BBDF gets smaller. Furthermore, by comparing our results with the solutions obtained by the existing methods, we can conclude that DI2BBDF yields less accurate approximations. However, for the case of high dimension problem, DI2BBDF still performs well although as known that block backward differentiation formula is more effective for solving stiff initial value problems. In Figure 3, we can observe that both exact and approximate solutions of DI2BBDF are still increase towards the stable fixed point 𝑦=1. For future research, we map out the strategy to improve the performance of this method by extending the order of the formula and modify the formula in a variable-step size for solving differential initial value problems with two fixed points.