Comments on the use of block methods for solving singular boundary value problems

Singular boundary-value problems appear frequently on the modellization of many physical phenomena as in catalytic diffusion reactions, chemical kinetics, thermal-explosion theory, or electro hydrodynamics, among others. The singular Lane-Endem equation is a typical kind of equation modelling some of those problems. Unfortunately, just in few occasions the exact solutions can be obtained. In this situation the block methods have been used largely for approximating different kind of differential problems. We propose its use for solving singular boundary value problems. The proposed strategy consist in a block method combined with an appropriate set of formulas which are developed at the first subinterval to circumvent the singularity at the left end of the integration interval. Some examples are presented to validate the efficiency of the proposed strategy.


Introduction
It seems that block methods were firstly proposed by Milne for solving initial-value problems in differential equations (see [1]). Since then, dozens of such block methods have appeared in literature. Notably scholars such as L. F. Shampine, S. O. Fatunla, S. N. Jator, Z. A. Majid, L. Brugnano, T. E. Simos, J. R. Cash, etc., have worked on block methods. Roughly speaking, a k-step block method is a set of k multistep formulas that produces at the same time k approximate values of the solution of the differential problem at k grid points. This may vary according to the order and the kind of equation considered.
It is well-known that for a second order differential problem y ′′ = f (x, y, y ′ ), the boundary conditions may be of different types: • Dirichlet boundary conditions, y(a) = y a , y(b) = y b • Neumann boundary conditions, y ′ (a) = y ′ a , y ′ (b) = y ′ b • Robin boundary conditions, g 1 (y(a), y ′ (a)) = v a , g 2 (y(b), y ′ (b)) = v b . * e-mail: higra@usal.es Let us see how to develop a two-step block method that uses third derivatives for solving a boundary value problem (BVP) with any of the above boundary conditions. Suppose that the function f satisfies the conditions to assure the existence and uniqueness of a solution. We consider a generic two-step interval [x n , x n+2 ] where x n+ j = x n + jh, and assume that the solution y(x) can be approximated on this interval by a polynomial p(x), that is, where the c j are unknown coefficients that will be determined. To determine those coefficients we impose that the following conditions hold where the y n+ j , y ′ n+ j , f n+ j = f (x n+ j , y n+ j , y ′ n+ j ), g n+ j = g(x n+ j , y n+ j , y ′ n+ j ) are respectively approximations for y(x n+ j ), y ′ (x n+ j ), y ′′ (x n+ j ) and y ′′′ (x n+ j ) with A similar approach was used in [2], resulting in the following formulas These formulas must be considered along the grid points in order to form a global method that provides approximate solutions of a considered BVP. That is, if we consider the grid points we take the formulas in (4)-(5) for n = 0, 1, 2, . . . , N − 2, making a total of 2(N − 1) formulas. The number of unknowns is 2N + 2 and we have two boundary conditions, so we need two additional formulas that can be obtained at any set of grid points. The global method consists in a system of 2N + 2 equations with 2N + 2 unknowns: y 0 , y 1 , . . . , y N , y ′ 0 , y ′ 1 , . . . , y ′ N . The solution of the system by some iterative scheme provides the approximate values of the BVP.

Singular boundary value problems of the form y
Now assume that our problem is of singular type. Those problems are characterized by the presence of a singularity in the coefficients of the differential equation. There are many models that use this type of problem in he applied sciences, for example in thermal explosions, in oxygen diffusion, in a non-linear heat conduction model of the human head, in human physiology, etc. (see [3,5,6]).
There are many approaches in the literature for solving these problems, as the variational iteration method, the homotopy perturbation method, different methods based on spline approximations, methods that use Padé or Taylor approximations, methods using different types of wavelets, the Adomian decomposition method, etc. (see [8][9][10][11]13]).
Let us consider the following singular BVP [3]: whose exact solution is If we try to solve this problem with the global method in the previous section we will run into troubles, since the value f 0 cannot be obtained.
To overcome this difficulty we propose to develop a set of multistep formulas to be applied at the first subinterval, [x 0 , x 1 ], with the purpose of specifically avoid the use of f 0 . A possibility to develop these formulas has been used in [7]. We consider three different intermediate points x r = x 0 + rh, x s = x 0 + sh and x t = x 0 + th with 0 < r < s < t < 1, on [x 0 , x 1 ], and suppose that the approximation to the true solution on this interval can be done by a polynomial, that is, Now, evaluating q(x) and its first derivative at x n , and its second derivative at the points x r , x s , x t , x 1 , we obtain a system of equations with the unknowns a n , n = 0(1)5, given by where similarly as before, the notations f j stand respectively for approximations of y ′′ (x j ). After solving this system we obtain the coefficients a n in terms of the parameters r, s, t. If we impose that the two first terms of the local truncation error for y 1 and the principal term of the local truncation error for y ′ 1 are null, we obtain the values r = 0.08858795951270394739554614376945, Using these values, the set of formulas for the first subinterval are readily obtained as

Main result
The formulas in (4) are used only for the first subinterval [x 0 , x 1 ], while for the remaining interval [x 1 , x N ] we consider the formulas in (8) for n = 1, 3, . . . , N − 2 with N odd. All these formulas give a global method for solving a singular BVP. It can be shown that this method is convergent, as stated in the following theorem: Theorem 3.1 Let y(t) be the exact solution of a SBVP y ′′ = f (x, y, y ′ ) subject to the boundary conditions y(a) = y a , y(b) = y b , and {y j } N j=0 the discrete solution provided by the proposed global method describe above. Then the proposed method is sixth order convergent, that is, for sufficiently small h, there exists a constant K (independent of h) such that We note that in the above theorem the boundary conditions considered are of Dirichlet type, but any other conditions may be considered. In the numerical examples in the following section we have considered mixed boundary conditions.

Example 1
We consider the singular BVP in (6) to asses the performance of the proposed method. Note that the solution presents at x = 0 a removable singularity. We have solved the problem considering a total number of steps M = 11, 21, 41. The obtained results are collected in Table 1 where we have included a row with the maximum absolute errors at the grid points,

Max Err(2h) Max Err(h)
) .  Figure 1 shows the exact solution and the discrete one obtained for M = 11.
We have solve the problem for r = 1/4, considering a total number of steps M = 11, 21, 41. The obtained results are collected in Table 2 where we have included as before the maximum absolute errors on the grid points and a numerical estimation of the order of convergence, given by ROC.

Future research
We have presented a strategy for solving singular boundary value problems using block methods. The numerical results presented show the good performance of the proposed strategy for a particular block method. Nevertheless, there are different tasks to be addressed in future works: • Analyze how the stepsize of the first subinterval can be chosen in order to get a better accuracy • Analyze which is the best option for joining the different formulas • Study how to choose the starting points in the Newtons method • Extend the approach to other problems, for example with singular points at both ends • Simplify the formulation of the method, for example reformulating it in such a way that the number of function evaluations is minimized (see [4]) • Reduce the number of equations in case of special problems where the first derivative is not present.