Diagonal Block Method for Solving Two-Point Boundary Value Problems with Robin Boundary Conditions

This numerical study presents the diagonal block method of order four for solving the second-order boundary value problems (BVPs)withRobin boundary conditions at two-point concurrently using constant step size.The solution is obtaineddirectlywithout reducing to a system of first-order differential equations using a combination of predictor-corrector mode via shooting technique. The shootingmethod 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.


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: () =  (, ,   ) for  ≤  ≤  (1) with where , ,  1 ,  2 ,  3 ,  4 , , and  are all constants and  1 ,  2 ,  3 , and  4 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 Bspline 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.

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  +1 and  +2 in these subintervals are simultaneously generated using the appropriate number of back values.The approximate solution of  +1 at the point  +1 will be computed using three back values at the points,   ,  −1 , and  −2 .While the approximation of  +2 at the point  +2 required an additional one back value,  +1 .The derivation formulae of  +1 and  +2 will be obtained by integrate equation (1) once and twice over the intervals [  ,  +1 ] and [  ,  +2 ], respectively, as follows.
Second Point The 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 points Now, introducing the variable substitution  =  +1 + ℎ and  = ℎ into interpolating polynomial, hence, evaluate these integrals using MAPLE with the limit of the integration from −1 to 0. 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 points Next, introduce the variable substitutions  =  +2 + ℎ and  = ℎ into interpolating polynomial.Again, evaluate these integrals using MAPLE with the limit of the integration from −2 to 0. 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. Predictor Corrector One-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.

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 as The local truncation error (LTE) associated with LMM in (17) is defined in the form of linear difference operator as Assume 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 in Simplifying this gives and now we obtained the expression where Definition 1.According to Fatunla [7] and Lambert [14], the method in ( 17) is said to be of order  with error constant This 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 as with By choosing  = 7 the calculation gives Thus, the 2PDD4 satisfied order four with the error constant,

Consistency of the Method
Definition 2. The linear multistep method is said to be consistent if it has order  ≥ 1 according to [14].
Since the order of the proposed method is  = 4 ≥ 1, the method is consistent.satisfies |  | ≤ 1 and for those roots with |  | = 1, the multiplicity must not exceed two as in Lambert [14] and See et al. [15].Rewrite (15) to (16) in matrix form as follows:
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 1 = ℎ and 2 = ℎ 2 .The boundary of the absolute stability region in 1 − 2 plane is determined by substituting  in the stability polynomial with 1, −1 and   for 0 ≤  ≤ 2. 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.

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.

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,  1 for (), then the initial guess of   () is given explicitly from the first boundary condition as follows: where  1 = / 1 and  1 =  2 / 1 .Compute the approximate solutions using the formula in ( 13) to (16).We obtained the first stopping condition as where ℎ((),   ()) =   () +  2 (), with  2 =  4 / 3 .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.
The first two estimates were decided to be  1 = 0 and  2 = 1 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.

Error
where r is the number of iterations and ()  is the  th component of the approximation. = 1, = 0 result in absolute error test,  = 1,  = 1 result in mixed error test, and  = 0,  = 1 result in relative error test.All the calculations were done using the C programing code.

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.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].
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  = 1, 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 temperaturedependent fluidity.In addition, if  = 0 and  > 0, then problem in (43) has a unique solution.
Data in Tables 1-4 summarized the result obtained for all tested problems at ℎ = 0.1 and ℎ = 0.01 with  = 10 −6 .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 ℎ = 0.1 but comparable in terms of computation time.Besides that, 2PDD4 achieved an excellent accuracy result compared to 2PDAM4.In detail, 2PDD4 has 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 ℎ = 0.025 and ℎ = 0.0125, the results showed that 2PDD4 are able to achieve smaller maximum error than the method in [6].  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.

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.

Figure 2 :
Figure 2: Stability region of the two-point diagonal block method.

MAXE:
Maximum absolute error as stated with * AVE: Average absolute error ℎ: Step size TOL: Tolerance TS: Total steps at last iteration FCN: Total function call ITN: Total iteration of guess TIME: Time computation in seconds 2PDD4: Direct two-point diagonal block method of order four proposed in this study 2PDAM4: Direct two-step Adams Moulton method of order four as in Phang et al.[17]

Figure 5 :
Figure 5: Graph of numerical result for solving Problem 3.

Figure 6 :
Figure 6: Graph of numerical result for solving Problem 4.

Table 1 :
Comparison of the numerical result for solving Problem 1.

Table 2 :
Comparison of the numerical result for solving Problem 2.

Table 3 :
Comparison of the numerical result for solving Problem 3.

Table 4 :
Comparison of the numerical result for solving Problem 4.

Table 3 .
Table 4 has shown that 2PDD4 required only one initial guessing as compared to 2PDAM4 and DAM4 for solving Problem 4 at ℎ = 0.1.Due to this efficient iterative performance, 2PDD4 gives a superiority result in terms of execution times and total function calls.