Research Article
Differential Quadrature Method for Singularly Perturbed Two-Point Boundary Value Problems
Department of Mathematics, Faculty of Science and Art, Pamukkale University, Denizli, 20070, Turkey
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 two-point boundary value problem arises in various fields of science such as fluid dynamics, control theory, chemical-reactor theory aerodynamics theory, elasticity and especially fluid motion,
(1) |
(2) |
ε 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 Polynomial-Based 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 Expansion-Based 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.
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 nth order derivative of the function u(x) at an intermediate point (grid point) xi can be written as:
(3) |
where, N is the number of grid points in the whole domain (a = x1 < x2 < ... < xN = b) and are the weighting coefficients of the nth 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. |
Polynomial-Based 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 one-dimensional differential equation is approximated by a (N-1)th degree polynomial:
(4) |
If , are the base polynomials in VN (N-dimensional linear vector space), then u(x) can be expressed by
(5) |
Here, the base polynomials wk(x), k = 1,2,..., N, are chosen as the Lagrange interpolated polynomials:
(6) |
(7) |
and xi, 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 second-order derivatives:
(8) |
(9) |
It can be understood from the above equations that the weighting coefficients of the second-order derivative can be completely determined from those of the first-order 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 so-called Chebyshev-Gauss-Lobatto points. For a domain specified by a≤x≤b and discretized by a set of unequally spaced points (non-uniform grid), then the coordinate of any point i can be evaluated by:
(10) |
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 = x1 < x2 < ... < xN = 1 where N is the number of grid points. We denote yi = y(xi), fi = f(xi), etc. Applications of the DQM, to discritize the derivatives in Eq. 1, lead to:
(11) |
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 N-2 equations with N unknowns. But we have two boundary conditions specified at both ends. These boundary conditions can be used to eliminate two unknowns (yl and yN) in Eq. 11.
If we have boundary conditions of the type given in Eq. 2, we obtain yl = A, yN = B. These can be substituted into the discritized Eq. 11 to obtain a (N-2)x(N-2) system of equations. This can be solved for the unknowns y2, y3, ..., yN-1.
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 non-uniform 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 non-linear problem which is compared to its possible solution in literature. The PDQ method is applied to the problems for the non-uniform grid points distribution given in Eq. 10. All computations were carried out using double-length arithmetic.
We have solved the following singularly perturbed two-point 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 (N-2)x(N-2) system of linear algebraic equations for the unknowns y2, y3, ..., yN-1. 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 not-a-knot 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 non-linear 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 non-linear 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 non-linear equations. This system is solved here using a FORTRAN IMSL routine called DNEQNF. This routine solves a system of non-linear equations using a modified Powell hybrid algorithm and a finite-difference 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.
In this study, the DQM has been applied to solve singularly perturbed two-point boundary value problems with a linear or non-linear nature. The applications presented here showed that the DQM has the capability of solving singularly perturbed two-point 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 non-uniform grids. This technique provides an alternative method to the conventional ways of solving singular perturbation problems.
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.