INTRODUCTION
Behaviour of many physical systems leads to a singularly perturbed differential
equation depending on a small physical parameter such that 0<ε<<1.
We consider the following class of singularly perturbed twopoint boundary
value problem arises in various fields of science such as fluid dynamics,
control theory, chemicalreactor theory aerodynamics theory, elasticity
and especially fluid motion,
with the boundary conditions
ε is the small parameter. Here q(x), p(x) and f(x) are assumed to
be sufficiently smooth, bounded and real functions as mentioned in (Kadalbajoo
and Patidar, 2003a; Khan et al., 2004). In many applied areas,
(1) possesses boundary layers, that is, regions of quick change in the
solution near the ends with widths o(1) as ε→0.
In recent years, much attention has been paid on finding the solutions
of these singularly perturbed equations. In literature, there are various
methods by which the solutions of these equations can be obtained. For
example, difference methods (Kellogg and Tsan, 1978; Berger et al.,
1981; Ilicasu and Schultz, 2004), finite element methods (Chen, 1997;
Chin and Krasny, 1983; Schatz and Wahlbin, 1983; Stynes and O`Riordan,
1986), differential transformation method (Chen and Liu, 1998), spline
methods (Flaherty and Mathon, 1980; Jain and Aziz, 1983; Sakai and Usmani,
1989; Stojanovic, 1996; Kadalbajoo and Patidar, 2003b), spectral methods
(Liu and Tang, 2001) and numerical integration method (Reddy and Reddy,
2002). In all these techniques, especially in finite difference techniques,
usually uniformity, convergence of the approximations, linearity, order
of the finite difference and some applications were successfully analyzed.
To the best knowledge of the author, the idea of the DQM has not been
implemented for the singular perturbation problems so far. The DQM is
an efficient discretization technique in solving initial and/or boundary
value problems accurately using a considerably small number of grid points.
Bellman et al. (1972) introduced the method in the early 1970s
and, since then, the technique has been successfully employed in finding
the solutions of many problems in applied and physical sciences (Shu,
1991; Shu and Richards, 1992; Yücel, 2006). The method has been projected
by its proponents as a potential alternative to the conventional numerical
solution techniques such as the finite difference and finite element methods.
In the DQ method, derivatives of a function with respect to a coordinate
direction are expressed as linear weighted sums of all the functional
values at all grid points along that direction. The weighting coefficients
in that weighting sum are determined using test functions. Among the many
kinds of test functions, the Lagrange interpolation polynomial is widely
employed since it has no limitation on the choice of the grid points.
This leads to PolynomialBased Differential Quadrature (PDQ) which is
suitable in most engineering problems. For problems with periodic behaviours,
polynomial approximation may not be the best choice for the true solution.
In contrast, Fourier series expansion can be the best approximation giving
the Fourier ExpansionBased Differential Quadrature (FDQ). The ease for
computation of weighting coefficients in explicit formulations (Shu, 2000)
for both cases is based on the analysis of function approximation and
linear vector space.
This study suggests the use of the DQM for solving the considered singularly
perturbed boundary value problems. The current study aims to demonstrate
that the DQM is capable of achieving high accuracy for the problem under
consideration.
MATERIALS AND METHODS
The DQM was presented for the first time by Bellman et al. (1972) as
mentioned before for solving differential equations. The method uses the basis
of the quadrature method in driving the derivatives of a function. It follows
that the partial derivative of a function with respect to a space variable can
be approximated by a weighted linear combination of function values at some
intermediate points in that variable.
In order to show the mathematical representation of the method, we consider
a single variable function u(x) on the domain a≤x≤b; then the n^{th}
order derivative of the function u(x) at an intermediate point (grid point)
x_{i} can be written as:
where, N is the number of grid points in the whole domain (a = x_{1}
< x_{2} < ... < x_{N} = b) and
are the weighting coefficients of the n^{th} derivative. As can
be seen from Eq. 3, two important factors control the
quality of the approximation resulting from the application of the DQM.
These are the values of weighting coefficients and the positions of the
discrete variables. Once the weighting coefficients are determined, the
bridge to link the derivatives in the governing differential equation
and the functional values at the mesh points is established. In other
words, with the weighting coefficients, one can easily use the functional
values to compute the derivatives. Note that for multidimensional problems
each derivative is approximated in the respective direction similarly.
In order to determine the weighting coefficients in Eq.
3, u(x) must be approximated by some test functions. To select a suitable
test function, one needs to satisfy the following conditions:
• 
Differentiability: The test function of the differential
equation must be differentiable at least up to the nth derivative
(here n is the highest order of the differential equation). 
• 
Smoothness: u(x) must be sufficiently smooth to be satisfied
the condition of the differentiability. 
PolynomialBased Differential Quadrature (PDQ): When the function
u(x) is approximated by a higher order polynomial, Shu (2000) and Shu
and Richards (1992) presented some explicit formulations to compute the
weighting coefficients within the scope of a higher order polynomial approximation
and a linear vector space. It is supposed that the solution of a onedimensional
differential equation is approximated by a (N1)^{th} degree polynomial:
If ,
are the base polynomials in V_{N} (Ndimensional linear vector
space), then u(x) can be expressed by
Here, the base polynomials w_{k}(x), k = 1,2,..., N, are chosen
as the Lagrange interpolated polynomials:
Where:
and x_{i}, i = 1,2,..., N, stand for the coordinates of grid
points which can be chosen arbitrarily. Substituting Eq.
5 into Eq. 3 and using Eq. 6 result
in the following weighting coefficients for the first and secondorder
derivatives:
It can be understood from the above equations that the weighting coefficients
of the secondorder derivative can be completely determined from those
of the firstorder derivative. Note that the explicit weighting coefficients
for the FDQ can be found in (Shu, 2000).
Choice of the grid point distributions: The selection of locations
of the sampling points plays an important role in the accuracy of the
solution of the differential equations. Using uniform grids can be considered
to be a convenient and easy selection method. Quite frequently, the DQM
delivers more accurate solutions using the socalled ChebyshevGaussLobatto
points. For a domain specified by a≤x≤b and discretized by a set
of unequally spaced points (nonuniform grid), then the coordinate of
any point i can be evaluated by:
Implementation of boundary conditions: In order to gain the accurate
numerical solution of differential equations, proper implementation of
the boundary conditions is also very important. Essential and natural
boundary conditions can be approximated by the DQM. Using the technique
in solving differential equations, the governing equations are actually
satisfied at each sampling point of the domain, so one has one equation
for each point, for each unknown. In the resulting system of algebraic
equations from the DQM, each boundary condition replaces the corresponding
field equation. This procedure is straightforward when there is one boundary
condition at each boundary and when we have distributed the sampling points
so that there is one point at each boundary.
Applications to singularly perturbed problems: For the approximate
solution of the boundary value problem Eq. 1 with the
boundary conditions given in Eq. 2 using the DQM, we
first discritize the interval [0,1] such that 0 = x_{1} < x_{2}
< ... < x_{N} = 1 where N is the number of grid points.
We denote y_{i} = y(x_{i}), f_{i} = f(x_{i}),
etc. Applications of the DQM, to discritize the derivatives in Eq.
1, lead to:
where,
and
are the weighting coefficients of the first and second order derivatives,
respectively. These weighting coefficients can be determined using the
explicit formulas given by Eq. 8 and 9.
Note that Eq. 11 should be applied at all interior
points to obtain a set of DQ algebraic equations. This will give us N2
equations with N unknowns. But we have two boundary conditions specified
at both ends. These boundary conditions can be used to eliminate two unknowns
(y_{l} and y_{N}) in Eq. 11.
If we have boundary conditions of the type given in Eq.
2, we obtain y_{l} = A, y_{N} = B. These can be substituted
into the discritized Eq. 11 to obtain a (N2)x(N2)
system of equations. This can be solved for the unknowns y_{2},
y_{3}, ..., y_{N1}.
Note that it is necessary to analyze the error resulting from the approximation
of a function and its derivatives. Shu (1991) has given a thorough error
analysis in his Ph.D thesis. Therefore it will not be discussed here in
this research.
We also note that the PDQ method is an extension of finite difference
methods and is actually the highest order finite difference scheme. Equation
4 can be applied to both interior points and boundary points and can
also be applied to a uniform mesh or a nonuniform mesh. As the highest
order finite difference scheme, the PDQ method is a global approximation
approach since it uses all the functional values in the whole computational
domain.
Numerical illustrations: To demonstrate the efficiency and accuracy
of the DQM, we have solved the following three linear problems whose exact
solutions are known as well as a nonlinear problem which is compared
to its possible solution in literature. The PDQ method is applied to the
problems for the nonuniform grid points distribution given in Eq.
10. All computations were carried out using doublelength arithmetic.
We have solved the following singularly perturbed twopoint boundary
value problems in the region 0≤x≤1 and for various values of ε.
Example 1: Consider the following problem:
with exact solution given by:
Table 1: 
Maximum absolute errors for each ε 

We solved this boundary value problem using the DQM with q(x) = 0, p(x)
= 1 and f(x) = 0 in Eq. 1. Applications of the DQM reduce
this boundary value problem to a (N2)x(N2) system of linear algebraic
equations for the unknowns y_{2}, y_{3}, ..., y_{N1}.
Although this system of equations can be solved by using several approaches,
we use here a FORTRAN IMSL routine called DLSARG which solves a real general
system of linear equations with iterative refinement.
In Table 1, the maximum absolute errors are given for various values
of N and ε. It can be observed from the table that the results obtained
are very accurate, even if the number of grid points is taken to be relatively
small, i.e., N = 8. As the number of grid points increases reasonably,
the errors decrease. This implies that the accuracy of the DQ results
can be improved by using a larger number of grid points.
Example 2: εy`` (x)+y`(x) = 1+2x, y(0) = 0, y(1) = 1 whose
exact solution is given by
This equation represents singularly perturbed linear equation and is
solved for q(x) = 1, p (x) = 0 and f(x) = 1+2x in Eq. 1.
The maximum absolute errors for various values of N and ε are given
in Table 2. As can be seen, using the reasonably small
number of grid points the method produces accurate results. Note that
decreasing in error is relatively slow for very small values of ε.
However, to overcome this drawback, number of grid points can be increased
at reasonable level.
In Table 3, we also compare the DQ results with the
exact solutions at selected points for different values of N and ε.
Note that the DQ results are given at uniform grids interpolated with
the use of a FORTRAN IMSL routine called DCSIEZ. This routine computes
the cubic spline interpolant with the notaknot condition and return
values of the interpolant at specified points. It can be seen that the
DQ results are seen to be in a very good agreement with the exact solutions.
To solve the system of linear algebraic equations for the unknowns here
and the following example, the same procedure with example 1 has been
used.
Table 2: 
Maximum absolute errors for each ε 

Table 3: 
Computational results for example 2 

Table 4: 
Computational results for example 3 

Example 3: In this example we consider the following homogeneous
singularly perturbed problem given by Bender and Orszag (1978) εy``
(x)+y`y(x) = 0, 0≤x≤1 with y(0) = y(1) = 1. The exact solution
is given by:
Where:
We solve this problem using the DQM for q(x) = 1 and p(x) = 1 with f(x)
= 0 in Eq. 1. The maximum errors in computed solutions
are given in Table 4 for various values of N and ε.
The agreement between the DQ solutions and the exact solutions is very
good. Similar observations can be made as in examples 1 and 2. The computational
results at grid points have been obtained with the use of same procedure
as in the previous example.
Table 5: 
Computational results for example 4 

The DQM can be used to solve nonlinear singular perturbation problems
as well as to linear singular perturbation problems. Even though we have
mostly intended to focus on linear singular perturbation problems, to
properly show the range of the DQ method we also include a nonlinear
singular perturbation problem as follows:
Example 4: Consider the following singularly perturbed problem
given by Bender and Orszag (1978):
with
y(0) = y(1) = 0.
Applications of the DQM reduce this problem to a system of nonlinear
equations. This system is solved here using a FORTRAN IMSL routine called
DNEQNF. This routine solves a system of nonlinear equations using a modified
Powell hybrid algorithm and a finitedifference approximation to the Jacobian.
We compare the DQ results with the Bender and Orszag`s uniformly valid
approximation
The comparison of the results can be seen in Table 5
for various values of N and ε. It can be noticed that the agreement
is considerably good.
CONCLUSIONS
In this study, the DQM has been applied to solve singularly perturbed
twopoint boundary value problems with a linear or nonlinear nature.
The applications presented here showed that the DQM has the capability
of solving singularly perturbed twopoint boundary value problems and
also is capable of producing highly accurate solutions with minimal computational
effort. The performance of the method for the considered problems was
measured by comparing with the exact solutions. It can be observed from
the tables that the DQM approximates the exact solution very well. This
shows the efficiency and accuracy of the method.
It has been observed that an increase in the number of grid points gives
rise to an increase in the accuracy of the DQM solution, as is the case
in most numerical techniques. However, using a considerably small number
of grid points in the DQM produces highly accurate results with the use
of nonuniform grids. This technique provides an alternative method to
the conventional ways of solving singular perturbation problems.
ACKNOWLEDGMENT
The author is grateful to Dr. U. Yucel (Pamukkale University, Faculty
of Science and Art, Department of Mathematics, Denizli, 20070, Turkey)
for his very useful comments on the paper.