1. Introduction

There is considerable literature on the numerical treatment of singular boundary value problems[1-4]. The numerical treatment of singular boundary value problems has always been far from trivial, because of the singularity at some points. Consider a class of nonlinear singular boundary value problems of the following form[4,9]

(1)

(2)

(3)

where α, β and are real numbers and we assume that f(x, y) є{ L2[0,1]×R} is continuous, exists and continuous, and. For the case m = 2, α = γand β = 1 the existence and uniqueness of the solution (1)-(3) has been given in[18].

This work is based on the Michaelis-Menten kinetics[5] for the steady state oxygendiffusion in spherical cells, in which

(4)

A similar problem arise in the study of the distribution of heat sources in the humanhead[6-8] in which

(5)

Point wise bounds and uniqueness results are given in[6] for this problem with f(x,y) of the form given by (4) and (5). Pandy and Singh[4] have used finite difference method (FD) based on uniform mesh, Kanth and Bhattacharya[9] usedcubic spline method of order O(h4) for solving (1)-(3) approximately. The objective of this paper is to use nonclassical pseudospectral method for approximation singular boundary value problem (1)-(3). Theoretical studies and numerical experiences have confirmed that for problems with smooth solution pseudospectral methods converge faster than other methods[10]. The idea of employing nonclassical weight functions first has been used by Shizgal[11] for solving the Boltzmann equations, Planck equations, and Shizgal and Heli Chen[12] used these basis for the solution eigenvalues and eigenfunction of Schrodinger equation.

2. Nonclassical Pseudospctral Discretization Method

The nonclassical pseudospectral methods[11,12,20] expand the function by using weighted interpolations of degree N of the form[11,12]

(7)

where x_j, j = 0, 1,…,N are a set of distinct collocation points on the interval[a,b], W(x) is a positive weight function, and L_j(x) are a set of Lagrange interpolating polynomials that satisfy L_j(xk) = δ_jk, i.e

2.1. Collocation Points

Consider the orthogonal polynomials p_n(x) with respect to some weight function w(x) on the interval [a,b], that is

The polynomials satisfy a three-term recurrence relation[13]

(8)

where the coefficients satisfy the inner product formulae [13]

(9)

(10)

(11)

where are obtained from the solution of the linear system of equations

(12)

The collocation points x_j, including the end points, are determined as the eigenvalues of J, and weights w_j obtained in terms of the first components v_1j of then ormalized eigenvector v as follows:

(13)

The integration of a function, f(x) may be approximated by the Gauss-Lobattoquadrature rule as

(14)

Where w_j are the weights associated with the collocation points x_j, j = 0, 1,…,N.

2.2. Computation of Differential Matrices

A number of algorithms have been proposed for generating pseudo spectral differential matrices numerically. In this section, the algorithm developed by Welfert[13] is employed and summarized as follows:

(1) Calculate the diagonal elements:

(a) Initialize using the equation

(b) Update recursively using the equation

(2) Calculate the off-diagonal elements:

(a) Initialize the off-diagonal elements using the equation

(b) Update using the following equation

For practical numbers of collocation points, the above algorithm is generally sufficient.

3. Discretization of Singular Boundary Value Problem

In this section, we solve the singular boundary value problem (1)-(3) by using nonclassicalpseudospectral method. From Eq. (1) we have that boundary valueproblem is singular at point x = 0. Now we collocate Eq. (1) in the points x_r, r = 1, 2,…,N -1 which at this points Eq. (1) is not singular. For this purpose we first substitute Eq. (7) for y(x) in Eqs. (1)-(3), i.e

(15)

(16)

(17)

Now collocate the equation (15) at x = x_r, r = 1, 2,…, N-1 and by using the differential matrices obtained in section 2, the above equations can be written as follows:

(18)

(19)

(20)

Now equations (18)-(20) are a system of nonlinear equations that can be solved by Newton iterative method for the resulting nonlinear systems.

4. Numerical Examples

In this section, we have used the method presented for different weighted functions w(x) and W(x) given in Table 1, and on two physical model examples: (i) oxygen diffusion; (ii) nonlinear heat conduction model of human head. These problems already been studied by Asaithambi and Goodman[19], Pandy and singh[4], Kanthand Bhattacharya[9]. The numerical results show that present method approximates the solution very well and computational time is less than other methods.

Table 1. Differentw(x) and W(x) on[0,1]

Table 2. Numerical results for N = 10 and different cases for Example 1

Table 3. Numerical results by other methods

Table 4. Numerical results for       and for N = 10.

Table 5. Numerical results for      and for N = 10

4.1. Example 1

Consider the oxygen diffusion corresponding Eqs. (1)-(4) with n = 0.76129, k = 0.03119, . Numerical results for cases given in Table 1, be given in Table 2, and numerical results with other methods be given in Table 3.

4.2. Example 2

The second example is of nonlinear heat conduction model of the human head, which correspond to Eqs. (1)-(3) and (5) with n = 1; k = 1. Numerical results by our method are given in Tables 4 and 6. We have performed calculations for the following two cases:

(i)

(ii)

Table 4, shows numerical results for case (i) and all weighted functions given in Table 1. Also in Table 5, numerical results for case (ii) together with results of cubic spline[9] of Kanth and Bhattacharya are given.

5. Conclusions

The nonclassical pseudo spectral method has been used to solve a class of singular boundary value problems arising in physiology. The main advantage of the presented approach is that it is possible to select arbitrary weight functions for the generation of the orthogonal polynomials. The cubic spline[9] and the method given in[4] have the disadvantage that the number of system is large and also the approximation is not good. The numerical results show that the proposed method is very accurate and needs less computational efforts respect to cubic spline[9] and the finite difference[4] methods.