A Numerical Method for Singular Boundary-Value Problems

This note is concerned with an iterative method for the solution of singular boundary value problems. It can be considered as a predictor-corrector method. Sufficient conditions for the convergence of the method are introduced. A number of numerical examples are used to study the applicability of the method.


Introduction
In this note we consider a numerical method for singular linear boundary value problems.Such problems arise very naturally in various applications including gas dynamics, chemical reactions, and structural mechanics.Traditional methods fail to produce good approximations for such equations.As a result, a number of investigators have considered various non-classical methods, including Chebyshev polynomials, B-splines, and cubic splines [1]- [3].Recent results also include methods based on reproducing kernel space [4] [5], and Sinc collocation method [6].
The purpose of this note is to develop an iterative method for singular and singularly perturbed boundary value problems.The method is explicit in nature, and can be considered to be an iterative predictor-corrector method.Section 2 introduces the method in details.Section 3 provides sufficient conditions for the convergence of the method.Section 4 uses a number of examples to investigate the applicability of the method, and compares the results to exact solutions.
( ) ( ) ( ), 0 With Dirichlet-type boundary condition where, the functions ( ) P x , ( ) ( ) . The functions ( ) P x and ( ) Q x can vanish at the boundary points.Assume that the domain is divided into e n equal intervals x ∆ , which leads to 1 e n n = + nodes with 1 0 x = , ( )  .The boundary conditions provide the val- ues for 1 y and n y .Consider a finite difference approximation of the above equation given by where ( ) ( ) The boundary conditions provide the values for 1 y and n y .It is then possible to formulate a two-step iterative formulation.First, using the given boundary conditions, an initial value can be assigned to the unknown function, i.e. j y , 1, , j n =  .Starting from the left ( ) 0 x = , and marching to the right, the first step is to solve Equation (3) for j y (here, we are naming it j z and refer to it as an intermediate variable) according to ( ) for 2, , 1 j n = −  Note that, marching from the left to right, the above equation is explicit for j z .The second step is to march from the right to the left according to ( ) z is an intermediate variable for which ( ) ( ) Marching from right to left, the above equations can be solved for ˆj y explicitly.The above two steps can be repeated af- ter setting ĵ j y y = , for 2, 1 j n = − .In the next section, we provide sufficient condition for convergence of the above iteration.

Convergence of the Method
For simplicity, consider a second-order singular differential equation given by The same analysis can be performed for a non-zero ( ) P x .The above two-step iterative method simplifies to ( ) ( ) In a matrix form, the first step can be written as The above equations can be written in a more compact form according to 2 1 1 where, the bold lower case letters indicate vectors corresponding to the above terms.After eliminating the intermediate term z in the above equations one can arrive at the equation given by ( )( ) For convergence of the above iteration, it is sufficient that all eigenvalues of the coefficient matrix, i.e., be inside the unit circle.The form of the coefficient matrices are such that it is possible to obtain explicit expressions for the inverses.For 1 A , we have ( ) ( ) Note that matrix 1 A is a lower triangular matrix, and its inverse is also a lower triangular matrix.It is also possible to explicitly compute the triple product ( ) The matrix 2 A is an upper triangular and its inverse which is an upper triangular is given by Writing the coefficient matrix in Equation (13) as a product of two matrices according to The first matrix on the right hand side is upper triangular, and the second matrix is lower triangular.The diagonal entries are the eigenvalues and using the spectral radius, it is sufficient to have 1 . This condition can be satisfied by choosing a small x ∆ for bounded values of ( )

Numerical Examples
Example 1: Consider a singular boundary value problem [4] given by ( ) ( ) ( ) ( ) The exact solution is given by ( ) 2 3 y x x x = − .Using a second-order accurate finite-difference approximation, the first step of the algorithm which is marching from the left ( ) to the right ( ) The second step is to march back from 1 x = to 0 x = using ( ) ( ) Both steps involve the explicit calculations only.After dividing the domain into equal intervals, Table 1 presents the results at selected points inside the domain and compares their values to the exact solution.The iterations are continued until the results do not change.In all the cases, the error is within the order of the approximation of the finite difference scheme.The relative error in the table is computed according to Example 2: Consider another singular boundary value problem [7] given by ( ) ( ) ( ) The exact solution is given by ( ) Dividing the domain into equal intervals and using the above method, Table 2 compares the computed solution with the exact solution at different points in the domain.The computed values are after 90,000 iterations.But the calculations are explicit and pose little burden.
Example 3: Consider a singularly perturbed boundary-value problem [7] given by for which the exact solution is given by ( ) ( ) This problem is singularly perturbed, and large gradients exist close to the boundaries.The present method can still be used to obtain an accurate solution.Table 3 presents the convergence of the solution close to the left ( ) boundary as a function of the mesh size, and Table 4 presents the convergence close to the right boundary ( ) The above results are for finite difference approximation with equal intervals.The accuracy of the results can be improved by using a smaller mesh sizes close to the boundaries.However, the calculations are explicit and using equal intervals does not pose a computational burden.

Conclusions
In this note we presented a numerical method for obtaining the solution of linear singular boundary value prob-   lems.The method can be considered as an iterative predictor-corrector method.Sufficient conditions for the convergence of the iteration was also presented.The proposed method is fully explicit and requires little computational time.It can also be applied to singularly perturbed boundary value problems.
Three numerical examples were used to study the applicability of the method.
− .Similarly, writing the second step in a matrix form leads to

Table 1 .
Computed value and the relative error at different values of x.

Table 2 .
Computed value and the relative error at different values of x for Example 2.

Table 3 .
Computed values and comparison to the exact solution close to x = 0.

Table 4 .
Computed values and comparison to the exact solution close to x = 1.