This numerical study presents the diagonal block method of order four for solving the second-order boundary value problems (BVPs) with Robin boundary conditions at two-point concurrently using constant step size. The solution is obtained directly without reducing to a system of first-order differential equations using a combination of predictor-corrector mode via shooting technique. The shooting method was adapted with the Newton divided difference interpolation approach as the strategy of seeking for the new initial estimate. Five numerical examples are included to examine and illustrate the practical usefulness of the proposed method. Numerical tested problem is also highlighted on the diffusion of heat generated application that imposed the Robin boundary conditions. The present findings revealed that the proposed method gives an efficient performance in terms of accuracy, total function calls, and execution time as compared with the existing method.

1. Introduction

This study is focusing on the numerical approach for solving second-order boundary value problems (BVPs) associated with Robin boundary conditions. Generally, this type of BVPs is given as follows:withwhere , and are all constants and , and all nonzero. In the case of Robin type, both functional value and derivative of the solutions are given in (2). If only functional value is given, then condition in (2) is known as Dirichlet type; otherwise BVPs will be subject to Neumann condition when only derivative values exist. Robin boundary conditions arise in several branches of applications such as in electromagnetic problem and heat transfer problem where these Robin type conditions are called impedance boundary conditions and the convective boundary conditions are, respectively, as explained in [1]. The Bernoulli polynomial together with Galerkin approximation in solving linear and nonlinear Robin boundary condition problems had been studied by Islam and Shirin [2]. The eminent scholars including Duan et al. [3] and Rach et al. [4] had derived the Adomian decomposition methods for solving BVPs imposing this condition where the approach involved analytic stage and numeric simulations. Meanwhile, Bhatta and Sastri [5] and Lang and Xu [6] had transformed respective BVPs into discretization form and solved them using symmetric global (continuous) spline and Quintic B-spline collocation method, respectively.

To optimize the computational cost, the development of proposed algorithm must at least satisfy the facility to generate solutions at several points simultaneously as suggested in Fatunla [7]. This implementation has been shown in the literature discussed by Phang et al. [8], Majid et al. [9], and Omar and Adeyeye [10]. They obtained the approximate solution of (1) at two points concurrently. Therefore, the advantages of outcome from their discussion have motivated us in this study.

Diagonal block method for solving differential equations was widely studied in the previous literatures. These include the discussion on solving first-order differential equations in [11, 12]. Solving second-order ordinary differential equations using diagonal block method has been discussed in [13]. The formulation in [12, 13] is based on backward differentiation approach. In [11], the author has derived the diagonal block method using Lagrange interpolation polynomial for solving first-order ordinary differential equations. In this research, we have extended the derivation in [11] to obtain the formulation of direct integration for solving second-order differential equations. Then the new formulae obtained have been implemented to solve the boundary value problems. The investigation will be focusing on the Newton divided difference interpolation as an iterative technique in seeking as well as updating the missing initial guesses while performing the shooting technique.

The organization of this paper is as follows. Section 2 presents the derivation of the two-point diagonal block method. Section 3 elaborates on the analysis of the method including the order, consistency, and stability. Section 4 presents the idea of the procedure of shooting method together with Newton divided difference interpolation as an iterative part. For validation and a clear overview, five tested problems will be discussed in Section 5. Finally, Section 6 concludes the finding from this study.

2. Derivation of the Diagonal Block Method

The interval of is divided into a series of blocks so that each block contains two subintervals with equal distance step size, , as depicted in Figure 1. Both numerical solutions of and in these subintervals are simultaneously generated using the appropriate number of back values. The approximate solution of at the point will be computed using three back values at the points, , and While the approximation of at the point required an additional one back value, The derivation formulae of and will be obtained by integrate equation (1) once and twice over the intervals and , respectively, as follows.

First PointSecond PointThe function in these equations will be approximated using Lagrange interpolating polynomial, By following the standard mathematical process, replace the in (4) and (6) with Lagrange interpolation polynomial that interpolates the set of pointsNow, introducing the variable substitution and into interpolating polynomial, hence, evaluate these integrals using MAPLE with the limit of the integration from to This will generate the corrector formula for the first point. Similar procedure applied to the in (8) and (10) with Lagrange interpolation polynomial that interpolates the set of pointsNext, introduce the variable substitutions and into interpolating polynomial. Again, evaluate these integrals using MAPLE with the limit of the integration from to This will generate the corrector formula for the second point. This two-point one block method is the combination of mode where P is the predictor formula, C is the corrector, and E is the evaluation of the function, The predictor formulae were derived similarly as the corrector formulae but with the order being one less. In this study, we called the proposed predictor-corrector formula as 2PDD4 method. The 2PDD4 formulae were given as follows.

PredictorCorrectorOne-step method will be used at the beginning of the proposed block method in order to get the starting initial points since multistep method needs more than one previous points before generating the remaining values over the interval.

3. Analysis of the Method

In this section, the order, consistency, stability, and convergence of the proposed two-point block method are discussed.

3.1. Order of the Method

The main proposed method of (13)–(16) is classified as a member of the Linear Multistep Method (LMM) which generally can be represented asThe local truncation error (LTE) associated with LMM in (17) is defined in the form of linear difference operator asAssume that is sufficiently differentiable, so that expanding the terms in (18) using Taylor’s series about the point will give the following expression:Substituting (19)–(21) into (18) results inSimplifying this givesand now we obtained the expressionwhere

Definition 1. According to Fatunla [7] and Lambert [14], the method in (17) is said to be of order with error constant ifThis concept was used to calculate the order and error constant of the proposed formula as stated in (13)–(16). Now, rewrite the corrector formula in matrix difference form aswith

By choosing the calculation givesThus, the 2PDD4 satisfied order four with the error constant,

3.2. Consistency of the Method

Definition 2. The linear multistep method is said to be consistent if it has order according to [14].
Since the order of the proposed method is , the method is consistent.

3.3. Stability of the Method

Definition 3. The linear multistep method is zero stable provided that the root of the first characteristics polynomial specified as satisfies and for those roots with , the multiplicity must not exceed two as in Lambert [14] and See et al. [15].
Rewrite (15) to (16) in matrix form as follows:From (30), the first characteristic polynomial , whereAccording to Definition 3 and (31), the diagonal two-point one block method is zero stable.

The test equation used in order to obtain the stability polynomial of the two-point block follows the idea from [11] as follows:The stability polynomial obtained is as follows:where and

The boundary of the absolute stability region in plane is determined by substituting in the stability polynomial with and for Figure 2 illustrates the regions of the absolute stability for the block method.

The stability region shows that the proposed numerical method will be able to produce reasonable results with the given values of the time steps.

3.4. Convergence of the Method

Definition 4. The linear multistep method is convergent if and only if it is consistent and zero stable [14].

Since the consistency and zero stability of the method have been achieved, in conclusion, the proposed two-point block method is convergent.

4. Implementation of the Method

In this study, shooting technique will be adapted throughout the solution process where the calculation starts with deciding a set of initial guesses. The behaviour of shooting method is ‘hit or miss’ the target. Due to that, we attempt to obtain the result as close as possible to the required solution with less number of guessing values. Therefore, choosing the right initial guesses and implementing the best method for refining the guessing value are two important elements in this shooting procedure. Initial values of and need to be guessed for the case of Robin type since none of them are given in (2). Starting with the initial guess, for , then the initial guess of is given explicitly from the first boundary condition as follows:where and Compute the approximate solutions using the formula in (13) to (16). We obtained the first stopping condition aswhere , with If this is sufficiently close to the condition in (35), then the BVPs are solved. Otherwise, update the new set of guessing values and the process continues as described in the following algorithm.

Algorithm of 2PDD4

Step 1. Set TOL and

Step 2. Set and calculate the approximate values, and , using the formula in (13) to (16) with where and depend on the test of convergence at each iteration.

Step 3. If , then repeat Step 2. If , then go to Step 4.

Step 4. If the following stopping condition is fulfilled: , then go to Step 6.
Else, continue to Step 5.

Step 5. Generate the new guessing values, and for based on the previous guess using Newton divided difference interpolation formula. Repeat Step 2.

Step 6. Complete.

The first two estimates were decided to be and based on the consideration in Roberts [16]. The calculation for the corrector part involved the test for convergence. The formulae for the calculation of error and convergence test are defined as follows.


Convergence Testwhere r is the number of iterations and is the component of the approximation. result in absolute error test, result in mixed error test, and result in relative error test. All the calculations were done using the C programing code.

5. Numerical Results

In this section, we have applied the algorithm of 2PDD4 to five numerical tested problems to illustrate its accuracy and efficiency. Problem 1 used mixed error test and Problem 2 applied the relative error test whereas Problems 3, 4, and 5 used absolute error test throughout the calculation for obtaining the required result. The following notations are used in the following result.MAXE: Maximum absolute error as stated with AVE: Average absolute error: Step sizeTOL: ToleranceTS: Total steps at last iterationFCN: Total function callITN: Total iteration of guessTIME: Time computation in seconds2PDD4: Direct two-point diagonal block method of order four proposed in this study2PDAM4: Direct two-step Adams Moulton method of order four as in Phang et al. [17]DAM4: Direct Adams Moulton method of order four as in Majid et al. [18]BP5: Bernoulli polynomials of degree five proposed by Islam and Shirin [2].

Problem 1. One has linear second-order differential equationwith and
Exact solution:
Source: Islam and Shirin [2].

Problem 2. One has nonlinear second-order differential equationwith and
Exact solution:
Source: Duan et al. [3].

Problem 3. One has nonlinear second-order differential equationwith and
Exact solution:
Source: Duan et al. [3].

Problem 4. One has nonlinear second-order differential equationwith and
Exact solution: Source: Bhatta and Sastri [5].

Problem 5. One has nonlinear second-order differential equationwith and
Exact solution:
Source: Lang and Xu [6].
The general equation of Problem 5 given as follows:arises in applications involving the diffusion of heat generated by positive temperature-dependent sources. According to discussion in Agarwal and O’Regan [19], p.295, if , the diffusion of heat arises in the two different analyses as Joule losses either in electrically conducting solids or in frictional heating. For the first analysis, represents the square of the constant current and the temperature-dependent resistance. Meanwhile, for the second part, represents the square of the constant shear stress and the temperature-dependent fluidity. In addition, if and , then problem in (43) has a unique solution.

Data in Tables 1–4 summarized the result obtained for all tested problems at and with In addition, the values ending with in Tables 1–4 denote the maximum error for those particular results. In order to investigate the competency of this proposed method, 2PDD4 has been compared with 2PDAM4 and DAM4. These two existing methods have the same order as 2PDD4 and will solve the BVP problems directly. The same shooting procedure has been implemented in the algorithm as explained in Section 4.

In Table 1, the maximum error of 2PDD4 is comparable with 2PDAM4 and DAM4 but better than BP5. The total function calls and execution times for 2PDD4 are less compared to 2PDAM4 and DAM4. As can be seen in Table 2, the total function calls and number of initial guessing for 2PDD4 are greater than DAM4 at but comparable in terms of computation time. Besides that, 2PDD4 achieved an excellent accuracy result compared to 2PDAM4. In detail, 2PDD4 has less three- or two-decimal accuracy than 2PDAM4 for solving Problem 2 with and , respectively. Method 2PDD4 needs less total function call and obtained smallest maximum error compared to other methods as displayed in Table 3. Table 4 has shown that 2PDD4 required only one initial guessing as compared to 2PDAM4 and DAM4 for solving Problem 4 at Due to this efficient iterative performance, 2PDD4 gives a superiority result in terms of execution times and total function calls.

Overall, a similar number of total steps are observed from all the data generated by 2PDD4 and 2PDAM4 as presented in Tables 1–4. These results justify that both methods have potential to reduce the total step approximately by half when computing the results at two-point simultaneously.

Finally, Table 5 describes the comparison of maximum error obtained in Lang and Xu [6] using the Quintic b-spline collocation method of order four with 2PDD4 method. The first three results demonstrate that the maximum error for both methods is comparable. However, at and , the results showed that 2PDD4 are able to achieve smaller maximum error than the method in [6]. Figures 3, 4, 5, and 6 display clear information about the comparison of the maximum errors versus the total function call as well as time for the numerical results summarized in Tables 1–4. Comparison result for Problem 5 is as depicted in Figure 7.

6. Conclusion

This study has presented a diagonal two-point block method formula of order four to solve directly the linear and nonlinear second-order Robin boundary conditions. Overview from the results showed a significant finding that the proposed diagonal block method manages to give faster execution time and less number of total function calls in all tested problems. The proposed method is able to achieve good accuracy when solving the problems. Future work should be devoted to solving third- and fourth-order Robin boundary value problems.

Data Availability

All data generated or analysed from the five tested problems in the article during this study are included in this published article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The authors gratefully acknowledge the financial support received from Putra Grant (Project Code: GP-IPS/2018/9625100), Universiti Putra Malaysia (UPM), and SLAB Scholarship sponsored by the Ministry of Higher Education (MOHE), Malaysia.