A Family of Boundary Value Methods for Systems of Second-Order Boundary Value Problems

A family of boundary value methods (BVMs) with continuous coefficients is derived and used to obtain methods which are applied via the block unification approach. The methods obtained from these continuous BVMs are weighted the same and are used to simultaneously generate approximations to the exact solution of systems of second-order boundary value problems (BVPs) on the entire interval of integration. The convergence of the methods is analyzed. Numerical experiments were performed to show efficiency and accuracy advantages.


Introduction
In what follows, we consider the general system of secondorder boundary value problems: where  : R × R 2 → R  are continuous functions, ,   , and   ∈ R  , and  is the dimension of the system.These second-order boundary value problems are encountered in several areas of engineering and applied sciences such as celestial mechanics, circuit theory, astrophysics, chemical kinetics, and biology.Most of these problems cannot be solved analytically, thus the need for a numerical approach.In practice, (1) is solved by the multiple shooting technique and the finite difference methods.The construction and implementation of higher order methods for the latter approach are difficult while the former approach suffers from numerical instability if the BVP is stiff [1][2][3] and singularly perturbed.
In the past few decades, the boundary value methods (BVMs) have been used to solve first-order initial and boundary value problems [4][5][6][7][8].Their stability and convergence properties have been fully discussed in [5].These BVMs are also used to solve higher order initial and boundary value problems by first reducing the higher order differential equations into an equivalent first-order system.This approach increases the computational costs and time and also does not utilize additional information associated with specific differential equations such as the oscillatory nature of some solutions [9,10].
Lambert and Watson [11] have derived symmetric schemes for periodic initial value problems of the special second-order   = (, ).Brugnano and Trigiante [4][5][6] have also derived BVMs for the first-order initial and boundary value problems.Amodio and Iavernaro [12] used BVMs to solve the special second-order problem   = (, ).Biala, Biala and Jator, Jator and Li [13][14][15] applied the BVMs to solve the general second-order problem   = (, ,   ) and Aceto et al. [16] constructed symmetric linear multistep methods (LMMs) which were used as BVMs for the special second-order problem   = (, ).In this paper, we have derived a class of BVMs and given a general framework via the block unification approach on how to use the BVMs on systems of BVPs for the general second-order differential equations (ODEs).

International Journal of Differential Equations
The boundary value technique simultaneously generates approximate solution ( 1 ,  2 , . . .,  −1 )  to the exact solution (( 1 ), ( 2 ), . . ., ( −1 ))  of (1) on the entire interval of integration.The BVMs can only be successfully implemented if used together with appropriate additional methods [5].In this regard, we have proposed methods which are obtained from the same continuous scheme and are derived via the interpolation and collocation approach [15,[17][18][19].
The paper is organised as follows.In Section 2, we derive a continuous approximation () of the exact solution ().Section 3 gives the specification of the methods.The convergence of the methods is discussed in Section 4. The use and implementation of the methods on ODEs and partial differential equations (PDEs) are detailed in Section 5. Numerical tests and concluding remarks are given in Sections 6 and 7, respectively.

Derivation of Methods
In this section, we shall use the interpolation and collocation approach [17] to construct a 2]-step continuous LMM (CLMM) which will be used to produce the main and additional formulas for solving (1).
Our starting point is to construct the CLMM which has the form where  0 (),  ] (), and   () are continuous coefficients and ] is chosen to be half the step number so that each formula, derived from (2), satisfies the root condition.The main and additional methods are then obtained by evaluating (2) at  + ( = 1(1)2],  ̸ = ]) to obtain the formulas of the form obtained from the first derivative of (2).Next, we discuss the construction of (2) in the theorem that follows.
Using Crammer's rule, the elements of  are determined and given as where   is obtained by replacing the th column of  by .
We rewrite (12) using the newly found elements of  as in (6); that is,

Specification of Methods
In this section, we specify the family of methods by evaluating the CLMM (2) at  + ,  = , ] − 1, ] + 1, . . ., 2], which is also used to obtain the derivative formula given by which is effectively applied by imposing that to produce derivative formulas of the form (4). Orders 4,6,and 8.For ] = 1, the BVM of order 4 is given as follows (where we have denoted a BVM with  step number as BVM):

BVM of
with the derivative formulas For ] = 2, we obtain the BVM of order 6 given as follows: with the derivative formulas For ] = 3, we obtain the BVM of order 8 given as follows: with the derivatives
We introduce the matrices  and  such that systems (24) and 4 are given by  −  () +  = 0, (26) and the exact form of the system is where 22 is an identity matrix so that ‖ −1 22 ‖ = 1.Thus, to obtain an estimate for ‖ −1 ‖, it suffices to show the existence of the inverse of  11 .Now, we define where  11 = diag( 11 ) so that Proof.
Theorem 4. Let  be an approximation of the solution vector  for the system obtained on a partition Proof.Subtracting (27) from (26), we obtain Under the conditions of Lemma 3,  −1 exists and is nonnegative.Therefore, International Journal of Differential Equations provided ℎ 2 ‖ −1 ‖‖‖‖‖ < 1.Hence,

Use of Methods
In this section, we discuss the use of methods in ( 16) and ( 17) for  = 0(2])( − 2]), where  is a multiple of 2].
We emphasize that the methods in ( 16) and ( 17) are all main methods since they are weighted the same and their use leads to a single matrix equation which can be solved for the unknowns.For example, for BVM6, we make use of each of the methods above in steps of 6; that is,  = 0, 6, . . .,  − 6.This results in a system of 2 equations in 2 unknowns which can be easily solved for the unknowns.Below is an algorithm for the use of the methods.The methods are implemented as BVMs by efficiently using the following steps.
Step 1. Use the methods in ( 16) and ( 17 Step 2. The unified block given by the system Step 1 results in a system of 2 equations in 2 unknowns which can be easily solved. Step 3. The values of the solution and the first derivatives of (1) are generated by the sequence {  }, {   },  = 0, . . ., , obtained as the solution in Step 2.
We note that all computations were carried out in Mathematica 10.0 enhanced by the feature FindRoot[ ].

Numerical Examples
In this section, we consider seven numerical examples.Examples 1 to 5 were solved using the BVMs  = 4,  = 6, and  = 8 (derived in this paper) of orders 6, 8, and 10, respectively.Also, these examples were solved using the Extended Trapezoidal Methods of the second kind (ETRs) and the Top Order Methods (TOMs) given in [5] of orders 6 and 10, respectively.Comparisons are made between the BVM  = 4 and the ETRs [5] as well as between the BVM  = 8 and the TOMs [5] by obtaining the maximum errors  in the interval of integration.We also compared our methods with the Sinc-Collocation method [20].Examples 6 and 7 were solved using the BVMs of order 6.We note that the number of function evaluations (NFEs) involved in implementing the BVMs is  × 2] in the entire range of integration.The code was based on Newton's method which uses the feature FindRoot[ ] or NSolve[ ] for linear problems in Mathematica.The efficiency curves show the plot of the logarithm of  against the number of function evaluations for each method.
Example 1.We consider the linear system of second-order boundary value problems given in [20] where This problem was solved using the ETRs and BVM of order 6 as well as the TOMs and BVM of order 10.The maximum Euclidean norm of the absolute errors in  1 and  2 was obtained in the entire interval of integration.In Table 1, we compared the Sinc-Collocation method [20] with the BVM of order 8. Table 2 shows the comparison between the ETRs, BVM4, TOMs, and BVM8.While the results of these methods are of approximate accuracy, we emphasize that the TOMs and ETRs use 20 function evaluations per step while the BVM4 and BVM8 use 8 and 16 function evaluations for this system.Hence, the BVMs are quite accurate and efficient.We also calculated the Rate of Convergence (ROC) using the formula log 2 ( 2ℎ / ℎ ), where  ℎ is the error obtained using step size ℎ.The ROC of the BVM4 and ETRs shows that these methods are consistent with the theoretical order (order 6) behavior of the methods.We omit the ROC of the TOMs and BVM8 because their errors are mainly due to round-off errors rather than to truncation errors.Figure 1 also shows the efficiency curves of these methods.

Example 2. Consider the nonlinear BVP given in [22]
where   2 () =  3 (1 − ) The maximum Euclidean norm of the absolute errors in  1 and  2 was obtained in the range of integration.Table 3 shows the comparison between the ETRs, BVM4, TOMs, and BVM8.While the results of these methods are of approximate accuracy, we emphasize that the TOMs and ETRs use 20 function evaluations per step while the BVM4 and BVM8 use 8 and 16 function evaluations for this system.We also calculated the ROC of the BVM4 and ETRs which shows that  these methods are consistent with the theoretical order (order 6) behavior of the methods.We do not calculate the ROC of the TOMs and BVM8 because their errors are mainly due to round-off errors rather than to truncation errors.Figure 2 shows the efficiency curves of these methods.
Example 3. We consider the nonlinear BVP with mixed boundary conditions given in [23] This problem was chosen to demonstrate the performance of the BVMs on nonlinear BVPs with mixed boundary  conditions.The maximum absolute errors were obtained in the range of integration.Table 4 shows the comparison between the ETRs, BVM4, TOMs, and BVM8. Figure 3 shows the efficiency curves of these methods.
Example 4. Consider the second-order BVP given in [24] (bvpT17) In order to assess the efficiency of our methods, we solve the boundary layer problem given in [24] (bvpT17).The maximum absolute errors were obtained in the range of integration.Tables 5 and 6 show the comparison between the ETRs, BVM4, TOMs, and BVM8 with  = 1 and 0.1, respectively.Figure 4  (45) Also, the efficiency of the scheme is shown by solving the problem given in [24] (bvpT20).The maximum absolute errors were obtained in the range of integration.Tables 7 and  8 show the comparison between the ETRs, BVM4, TOMs, and BVM8 with  = 1 and 0.1, respectively.This example shows the performance of the BVMs on the Poisson equation.In order to solve the equation using the BVMs, we carry out the semidiscretization of the spatial variable  using the second-order finite difference method International Journal of Differential Equations to obtain the following second-order system in the second variable : where Δ = (−)/,   = +Δ,  = 0, 1, . . ., , u = [ 1 (), . . .,   ()]  , g = [ 1 (), . . .,   ()]  ,   () ≈ (  , ), and   () ≈ (  , ) which can be written in the form subject to the boundary conditions u( 0 ) = u( /2 ), u(  ) = u  , where f(, u) = Au + g and A is an  − 1 ×  − 1 matrix arising from the semidiscretized system and g is a vector of constants.Table 9 shows the comparison between the BVM and the method in [25].Figure 6 shows the plot of the exact, approximate, and error function of the problem.

Conclusions
This paper is concerned with the solution of systems of second-order boundary value problems.This has been achieved by the construction and implementation of a family of BVMs.The methods are applied as a block unification method to obtain the solution on the entire interval of integration.We established the convergence of the methods.
We have also shown that the methods are competitive with existing methods cited in the literature.
In the future, we would like to develop a variable step size version of the BVMs with an automatic error estimation.

Example 5 .
shows the plot of the solution for values of  = 1, 0.1, 0.01, 0.001 and the solution has a boundary layer at  = 0. Consider the second-order BVP given in[24]

Figure 6 :
Figure 6: Plot of solution for the Poisson equation.

Figure 7 :
Figure 7: Plot of solution for the Sine-Gordon equation.

Table 1 :
Maximum errors for Example 1.

Table 2 :
Maximum Errors for different stepsizes for Example 1.

Table 3 :
Maximum Errors for different stepsizes for Example 2.

Table 4 :
Maximum errors for different step sizes for Example 3.

Table 5 :
Maximum Errors for different stepsizes for Example 4 for  = 1.

Table 6 :
Maximum Errors for different stepsizes for Example 4 for  = 0.1.

Table 7 :
Maximum Errors for different stepsizes for Example 5 for  = 1.

Table 8 :
Maximum Errors for different stepsizes for Example 5 for  = 0.1.

Table 9 :
Maximum error for the Poisson equation on  = 1.

Table 10 :
Maximum errors for the Sine-Gordon equation.