AN EXPONENTIALLY FITTED INTEGRATION SCHEME FOR A CLASS OF QUASILINEAR SINGULAR PERTURBATION PROBLEMS

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract. A class of quasilinear singularly perturbed two point boundary value problems in ordinary differential equations exhibiting boundary layer at the left end of the underlying interval is considered and a new exponentially fitted integration scheme on a uniform mesh is devised for its numerical solution. The devised scheme is obtained by introducing a suitable constant fitting factor in a new three term recurrence relationship derived by the application of outer region solution obtained by asymptotic expansion procedure and a combination of the exact and approximate rules of integration with finite difference approximations of the first derivative. Value of the fitting factor is determined using the theory of singular perturbations and used to take care of rapid changes in the solution. Resulting tridiagonal algebraic system of equations is solved by the Thomas algorithm. Convergence of the scheme is analyzed. Three numerical example problems are solved and computational results are tabulated to show the accuracy and efficiency of the method. It is easily observed that the derived scheme is able to produce accurate results with minimal computational effort for all the values of the mesh size h when the perturbation parameter ε tends to zero. Both the theoretical and numerical analysis of the method reveals that the method is able to produce uniformly convergent results with quadratic convergence rate.


INTRODUCTION
The problem of finding solution of a differential equation in which highest order derivative term is multiplied with a small positive parameter ε is known as singular perturbation problems. The investigation of numerous physical problems in science and engineering leads to the boundary value problems for singularly perturbed differential equations whose solutions show the capricious behaviour when the values of the ε tends to zero. That is, there exist sharp boundary and/or interior layer(s) in which solution varies most quickly while away from the layers(s) solution behaves regularly and varies slowly. These sharp/thin layers are usually referred to as shock/edge layers in solid mechanics, transition points in quantum mechanics, skin layers in electrical applications, and the boundary layers in fluid mechanics. Singularly perturbed problems are most ubiquitous due to such layer behaviour of the solutions. The existence and capriciousness of such problems can easily be observed in the modelling of various modern complicated processes, such as the fluid flow at high Reynolds numbers, water quality problems in river networks, convective heat transport problem with large Peclet numbers, drift diffusion equation of semiconductor device modelling, magneto-hydrodynamics duct problems at high Hartman numbers, electromagnetic field problem in moving media, chemical reactor theory, and financial modelling of option pricing. The development of a suitable numerical or analytical method for the solutions of singularly perturbed problems(SPPs) has been a challenging and interesting task for the applied scientists and engineers due to the failure of classical standard methods in providing accurate reliable results in case when ε tends to zero. Existing standard numerical methods produces an oscillatory or unsatisfactory results. The numerical treatment of a singularly perturbed quasilinear problem exhibiting layers is more difficult than the ordinary SPPs due to the appearance of the non linear term in the equation. A variety of analytical and numerical techniques have been suggested for finding the solution of such type of singularly perturbed problems. For the more detailed information about SPPs with/without quasilinear term and its treatments, we may refer to the high level monographs [3-9, 16, 19-22, 27-28] and the articles [1, 10-14, 17-18, 23-26, 29-30].
In the recent past, the authors in the articles [13] presented an initial value technique for solving quasilinear singularly perturbed two point boundary value problems. Technique is based on division of the given region into two inner and outer region; the asymptotic expansion; and reduction of the the original problem in to an asymptotically equivalent first order initial value problem with an exponentially fitted finite difference scheme. Authors in [10] described an approximate method in which they divided the given region into two inner and outer region and using the deviating argument they solved first order problems over each of the regions. In this paper, we proposed a simple and computationally efficient integration scheme for the solution of the quasilinear singularly perturbed problem [13] exhibiting a layer at the left end point of the underlying domain. Method is based on asymptotic expansion and the usual rules of the exact and approximate integration with finite difference approximations of the first derivative.
The main benefit of this proposed scheme is that it neither depend on very fine mesh [26] nor any other parameter like deviating arguments, transition parameter [10,13]. Further there is no need of the reduction of the original problem into an asymptotically equivalent first order initial value problem by approximating the zeroth order term by outer solution obtained by asymptotic expansion [13].
The outline of this paper is as follows: Problem statement with the proposed method of solution is described in section: 2. The convergence analysis is given in section: 3. Numerical illustration with computational results for the considered three example problems are presented in section: 4. The conclusions are presented in section: 5. The paper ends with the references.

STATEMENT OF THE PROBLEM
Consider a class of singularly perturbed quasilinear two point boundary value problems in ordinary differential equations with layer behaviour of the form: under the boundary conditions (2) x(0) = α and x(1) = β where ε is a small positive perturbation parameter (0 < ε ≤ 1); α, β are given finite constants; and the functions f (t), g(t, x), r(t) are defined in such a way that the boundary value problem (1) with (2) is guaranteed to have a unique solution (see [13,15]

DESCRIPTION OF THE METHOD FOR LEFT-END BOUNDARY LAYER PROBLEMS
To describe the method, we consider the division of the interval [0, 1] into N equal parts/ subintervals by means of the mesh/grid points The points t 0 and t N are called the boundary grid points, while the points t 1 ,t 2 ,t 2 , ....,t N−1 are called interior grid points. Now, we write the equation (1) in the form: Replacing x(t) by x out (t) in the right hand side expression in equation (3) we get: where x out (t) is the asymptotic expansion for x(t) in the form [19][20][21][22]: On substituting x(t) from equation (5) into equation (1) we get with the equation (6).
x 0 (t)+x 1 (t)(ε)+....] using Taylor's series expansion procedure and equating the coefficients of like power of ε in the equation (6) and equation (7), we get successively and hence x(t) given by equation (5) is obtained.

Integrating the equation (4) in the interval
Now using the non symmetric first derivative finite difference approximations: in the equation (9) and simplifying we get By introducing a constant fitting factor into equation (10), we obtain the proposed exponentially fitted three term integration scheme: where value of the fitting factor σ (ρ) is to be determined in such a way that the solution of (11) converges uniformly to the solution of the equations (1) with (2). The procedure for determining the value of the fitting component introduced in the equation (10) is described in the following subsection.

Determination of the fitting factor for Left-End Boundary Layer Problems
In this section, we will determine the value of the fitting factor introduced in the equation (10) using the theory of singular perturbations.
Taking limits as h → 0 in equation (11), we obtain The equation (12) at t = t i = ih becomes: Now, we have the solutions of equations (1) with (2) in the following form (cf. [19], pp.22-26): Using Taylor's expansions for f (t) about the point t = 0 upto their first terms only, the equation (14) becomes: Furthermore, considering the equation (16) at the point t = t i = ih, (i = 0, 1, 2, ..., N) and taking the limit as h tends to zero, we obtain Using the equation (17) (13) and then simplifying, we get the value of the fitting factor as Finally, by making use of equation (11) and σ (ρ) given by equation (18), we get the following exponentially fitted three-term integration scheme of the form: Here, the relation (19) represents a system of (N − 1) equations with (N + 1) unknowns x 0 to is an increasing function of t with f (t) > 0 in [0, 1], the inequalities: , and E i < G i are satisfied for the coefficient matrix of the resulting tridiagonal system of equations. Hence, the condition system (19) together with the boundary conditions equation (2) is irreducibly diagonally dominant and non-singular.We have utilized the Thomas Algorithm portrayed in [10,23] for solving the resulting tridiagonal system of equations.

CONVERGENCE ANALYSIS
In this section, we discuss the convergence analysis of the method in similar way as given in [2]. For this purpose we write the tri-diagonal system (19) in matrix-vector form, we get we also have where X = X 0 , X 1 , X 2 , ..., X N t denotes the actual solution and τ(h) = (τ 1 (h), τ 2 (h), ..., τ N (h)) t is the local truncation error.
From the equations (20) and (22), we have Thus, we obtained the error equation is where E = X − X = (e 0 , e 1 , e 2 , ..., e N ) t . Let x i be the sum of elements of i th row of D, then we Since 0 < ε << 1 for sufficiently small D the matrix is irreducible and monotone Hence D −1 exists and D −1 ≥ 0. Hence, from the error equation (24) we have Let u k,i be the (k, i) th element of D −1 since u k,i ≥ 0, from the theory of matrices we have, for some i 0 lies between 1 and N-1.
Therefore, from the equations (21), (25) and (26), we obtain and therefore The above results can be summarized as follows: Theorem 1. Let x(t) be the exact solution of singularly perturbed second-order quasilinear differential equation (1) with the boundary conditions (2) and let x i be the numerical solution obtained by the scheme (19) with the boundary conditions equation (2). Then, for sufficiently small h, relation (19) with equation (2) provides the solution with the quadratic convergence rate on uniform mesh.

NUMERICAL ILLUSTRATION
To examine the effectiveness of the derived scheme, we carried out the numerical analyses on the three standard test model for different values of N = 1/h and perturbation parameter ε, and presented the computational results in term of maximum absolute errors (E ε ) and numerical rate of convergence (R N ) in the Tables: 1, 2 and 3 respectively. One can easily observe from the tables that the method is capable of solving the problems accurately with quadratic rates of convergence.
It should be noted that the maximum absolute errors (E ε ) and numerical rate of convergence (R N ) for all the three example problems are computed by the double mesh principle [5]: 2i and the formula [5]: , respectively. Example 1: Consider the following non-linear singular perturbation problem from [13]: with boundary conditions x(0) = 0 and x(1) = 0.
For this example, outer region solution is x out = log 2 1+t The exact solution of this example problem-1 is given by The computational results (Maximum absolute errors and numerical rates of convergence) for example problem-1 are introduced in Table-1, for different estimations of ε and N. Example 2: Consider the following homogeneous singular perturbation problem from [13]: where outer region solution is x out = 1 2−t . The exact solution for this example problem-2 is given by The computational results (Maximum absolute errors and numerical rates of convergence) for example problem-2 are introduced in Table-2, for different estimations of ε and N. Example 3: Consider the following singular perturbation problem from [13]: with outer region solution x out (t) = 1 1+t . The exact solution for this example problem-3 is given by The computational results (Maximum absolute errors and numerical rates of convergence) for example problem-3 are introduced in Table-3, for different estimations of ε and N.

CONCLUSION
We have developed and described an exponentially fitted tridiagonal scheme for solving singularly perturbed quasilinear boundary value problems with boundary layer at left end point of the domain. To examine the performance of the proposed integration scheme, we performed the numerical experiments on three standard test example problems for different values of N and ε, and tabulated the computational results in terms of the maximum absolute errors and numerical rates of convergence. We carried out the convergence analysis for the proposed scheme.
The numerical investigations performed approves the theoretical rate of convergence obtained.
Further, we effectively see from these Tables: 1, 2, and 3 that the introduced fitted plan is fit for approximating the solution well with insignificant computational exertion, when perturbation parameter ε tends to zero for any fixed value of the mesh size h = 1 N . As stated in the introduction section, the proposed scheme neither depend on very fine mesh [26] nor any other parameter like deviating arguments and the transition parameter [10,13]. There is no need of the reduction of the original problem into an asymptotically equivalent first order initial value problem by approximating the zeroth order term by outer solution obtained by asymptotic expansion [13]. Finally, the present method may considered as one of the best choices among the methods developed for solving the considered quasilinear problem (1) with (2) with a small amount of computational time and problems preparation.

CONFLICT OF INTERESTS
The author(s) declare that there is no conflict of interests.