Abstract
As it is well recognized that conventional numerical schemes are inefficient in approximating the solutions of the singularly perturbed problems (SPP) in the boundary layer region, in the present work, an effort has been made to propose a robust and efficient numerical approach known as element-free Galerkin (EFG) technique to capture these solutions with a high precision of accuracy. Since a lot of weight functions exist in the literature which plays a crucial role in the moving least square (MLS) approximations for generating the shape functions and hence affect the accuracy of the numerical solution, in the present work, due emphasis has been given to propose a robust weight function for the element-free Galerkin scheme for SPP. The key feature of nonrequirement of elements or node connectivity of the EFG method has also been utilized by proposing a way to generate nonuniformly distributed nodes. In order to verify the computational consistency and robustness of the proposed scheme, a variety of linear and nonlinear numerical examples have been considered and errors have been presented. Comparison of the EFG solutions with those available in the literature depicts the superiority of the proposed scheme.
1. Introduction
Singularly perturbed problems (SPP) and their numerical solutions have attained great attention from researchers from the past few decades. Pearson [1] is credited with being the first researcher to solve the SPP numerically using traditional three-point finite-difference techniques on nonuniform meshes. The problem of nonsmoothness in the solution was firstly overcome by utilizing layer adaptive mesh methodologies. Bakhvalov [2] was the first who proposed special type of meshes, known as Bakhvalov meshes, to capture the boundary layers found in reaction-diffusion problems’ solutions, and later on, these meshes were used and modified by Gartland [3] and others for convection-diffusion problems. Later on, Shishkin [4] introduced another special type of meshes, called Shishkin meshes, to generate parameter-uniform numerical schemes under the finite-difference framework. Because of the simple structure of the Shishkin meshes, these meshes were utilized by many other researchers also for generating parameter-uniform numerical schemes for singularly perturbed problems. With the passage of time, other numerical schemes, called fitted operator schemes, were also developed by mathematicians and other researchers for solving singularly perturbed differential equations. For instance, Roos et al. [5] proposed an exponentially fitted finite element method for solving singularly perturbed problems.
Recently, considering the benefits and advantages of mesh-free methods over other numerical schemes, researchers have paid considerable attention to the developments of mesh-free methods. To name a few, the diffuse element method by Nayroles et al. [6], the EFG method by Belytschko et al. [7], reproducing kernel particle method (RKPM) by Liu et al. [8], the local Petrov–Galerkin method by Atluri and Zhu [9] etc., are some very famous mesh-free methods. However, for singularly perturbed problems, rarely, any mesh-free method has been proposed. Therefore, the main motive of the present work is to propose a suitable mesh-free method for SPP.
Geng and Qian [10] employed the reproducing kernel method for obtaining numerical solutions of singularly perturbed turning point problems having dual boundary layers. The authors [11] devised a modified reproducing kernel method (RKM) to solve singularly perturbed delay boundary value problems later in 2015. The suggested approach was based on the reproducing kernel theory, and the method’s error estimates were also discussed. Nadjafi and Ghassabzade [12] proposed another MLS-based meshless technique to deal with singularly perturbed delay differential equations (DDEs). Another meshless approach, namely, the radial basis collocation method with the coordinate stretching procedure was presented by Akhavan Ghassabzade et al. [13] for the mathematical treatment of singularly perturbed DDEs.
In this study, we have investigated the following nonlinear singularly perturbed problem:under essential boundary conditions defined bywhere is a small positive parameter known as singular perturbation parameter with , and we suppose that is continuously differentiable function. These types of problems appear in almost all the fields of engineering and science, e.g., plasma dynamics, biochemistry, electric power, aerospace systems, geophysics, fluid mechanics, and heat transport problems [14].
In the present work, we have proposed the element-free Galerkin (EFG) method, one of the most popular and versatile mesh-free methods, for approximating nonlinear SPP solutions. The EFG technique is based on the theory of moving least squares (MLS) approximation and standard weight functions to construct the shape functions. The method employs background cells to compute numerical integrations and is based on the global weak form. Since the Kronecker delta condition is not satisfied by the MLS shape function approximants, as a result, the imposition of essential boundary conditions is not straightforward. To overcome this problem, Lagrange multiplier approach [7] has been utilized. Looking on the bright side, the EFG technique is one among the very efficient techniques for approximating the solutions of nonlinear differential equations. Due to the absence of element connectivity, adaptive refinement of the node particles can be achieved very efficiently and easily. These features make the EFG method more adaptable than the finite element method (FEM). On the contrary, the EFG technique still has a few downsides also in comparison to FEM. The EFG method is computationally more expensive and has a more complex implementation algorithm than FEM. Since the method has not been applied to SPP so far, therefore, after analyzing various weight functions in the present work, efficient weight functions have been proposed for SPP.
The paper is organized as follows.
The MLS approximation and construction of shape functions have been discussed in Section 2. Different weight functions and their properties have been presented in Section 3. In Section 4, the quasi-linearization technique and its convergence have been presented. Section 5 describes the element-free Galerkin (EFG) technique. Section 6 contains the algorithm of the proposed method and numerical results, followed by a conclusion and bibliography.
2. Moving Least Squares Approximation
The field function is approximated by using MLS approximants in the domain . These approximants contain three components: basis functions which are usually polynomials, weight functions that are related with nodes present in the influence domain of a particular node, and the set of coefficients that depend on the position of nodes. The domain of influence of any nodal point is defined as a local neighborhood of or the domain that the particular node exerts an influence upon. The MLS’s approximation of the unknown function is defined as follows:where is a vector of monomial basis functions, is the number of components in the basis vector , and is an unknown vector of coefficients of the basis functions.
The choice of basis functions is essential in the construction of moving least squares based shape functions. The following are the examples of the basic functions:(i)Linear basis:(ii)Quadratic basis:
The field function’s approximate value at any particular node in the support domain of is defined by
We define the weighted residual functional as follows to determine the unknown coefficient vector in the support domain of the point of interest :where denotes the number of nodes in the support domain of the point of interest . The nodes in the support domain of any point are used to support or approximate the function value at . The dimension of the support domain of the point is defined bywhere signifies the dimensionless size of the support domain and the characteristic length denotes the nodal spacing for the point . in expression (7) is some weight function which is generally chosen as monotonically decreasing as increases. In the literature, numerous weight functions have been introduced and employed by researchers, as mentioned in Section 3. In MLS approximation, the condition generally leads to the ill-conditioned system of equations.
On simplification, equation (7) can be rewritten aswhereand is the diagonal matrix defined as
The minimization of with respect to leads to
This giveswhere
Therefore, from (3), the approximated field function can be computed aswhere , are the shape functions defined by
The graphs of typical shape functions for different nodes have been plotted in Figure 1. We can observe that the shape functions are bell-shaped. The derivatives of the MLS shape functions can be obtained by differentiating (16) with respect to the spatial coordinates.
3. Choice of Weight Functions
Weight function corresponding to any node provides weights for the residuals at different nodes within the support domain of the node of interest, and it ensures that nodes should enter and leave the support domain smoothly so that the shape functions satisfy the compatibility condition and the approximation is globally continuous. The weight functions must be continuous and positive in the corresponding support domains. The weight functions should satisfy the following properties:(i) in the support domain of the node (ii) outside the support domain (iii) is a monotonically decreasing function(iv) as
represents the Dirac delta function and is measure of the size of the support domain . The weight functions are usually based on a normalized distance , where denotes the radius of the influence domain. In [15–23], there are variety of weight functions for the MLS approach. Some of the weight functions, which are widely used in literature, are defined below. Though these weight functions are utilized by many researchers, here we have referenced only a few of them.(i)Quadratic spline: the weight function of a quadratic spline is defined as follows: Singh [15] compared different weight functions (exponential, elliptical, cosine, and quadratic) under the EFG method for heat transfer problems and found that elliptical, quadratic, and cosine weight functions provide acceptable results for the support domain with radius , where the scaling parameter satisfies , whereas exponential weight functions provide good results for the scaling parameter . In overall comparison, the author found that the exponential weight functions provide much better results in comparison to the other weight functions (i.e., elliptical, cosine, and quadratic) for heat transfer problems.(ii)Cubic spline: the weight function of a cubic spline is defined by Cubic spline weight functions were firstly proposed by Dolbow and Belytschko [16] in the EFG method. The authors constructed cubic spline weight functions to impose continuity requirements while solving linear elastostatic problems. These are some of the most frequent weight functions which are widely adopted by many researchers while employing the element-free Galerkin methodology, e.g., Park and Leap [17] utilized the EFG method with cubic spline weight functions for solving groundwater flow problems. Yang and He [18] utilized these weight functions in the EFG method for heat transfer problems. Juan et al. [19] used cubic spline weight functions to solve the topology optimization problem of the continuum structures. Zhang et al. [24] also applied the EFG method with cubic spline weight functions for two-dimensional elastodynamics problems. Rajesh Sharma [25] also successfully applied these weight functions along with the EFG method to study the heat transfer fluid flow problems over a stretching sheet.(iii)Quartic spline: the quartic spline weight function is given by These weight functions are also widely used by researchers as they provide continuity. Shibata and Murakami [26] presented a stability procedure for soil-water coupled problems using the EFG method with quartic spline weight function. He et al. [27] established the hybrid EFG method and the scaled boundary method known as the EFG scaled boundary method. The authors employed quartic spline weight function in EFG formulation to carry out numerical results. Jaberzadeh et al. [28] also considered quartic weight function in EFG methodology to study thermal buckling of trapezoidal plates. Cheng et al. demonstrated improved the complex variable EFG method for advection-diffusion [29]. The impact of cubic and quartic spline weight functions on computational precision has been discussed in both the articles. The proposed method has higher computational accuracy while employing quartic spline weight functions. The quartic spline weight functions in the EFG technique have been adopted by various other researchers also for solving different types of problems such as nonlinear shallow water equations [20] and static and buckling analysis of the thin plate [30].(iv)Gaussian exponential spline weight function: the Gaussian exponential spline weight function is defined as Here, the parameter is generally chosen from the interval . The Gaussian weighting function of exponential type is one of the most commonly used weight functions. These weight functions have been initially utilized by Belytschko et al. to solve elasticity and heat conduction problems [7], crack propagation problems [31], and fracture problems [23]. Following these, many researchers employed Gaussian exponential weight functions for solving various types of other problems such as fractional cable equation [32], quasi-three-dimensional free vibration analysis of multilayered composite plates [33], and bimaterial problems [34]. The performance of Gaussian weight functions has been analyzed by comparison with other weight functions in [35, 36].(v)Conical weight function: conical weight function is defined as Liu et al. [37] performed EFG formulation with conical weight function to obtain the numerical solutions of 3D electromagnetic problems. The authors compared the performance of different weight functions (cubic spline, rational formula, conical formula, negative exponential, and quartic spline) by correlating the EFG method results with the exact one. The graphical results indicate that error estimation of conical formula increases by increasing the value of scaling parameter . Zhuang et al. [35] highlighted issues and sources of errors for meshless methods. The control of errors has been shown by considering three types of weight functions, i.e., a rational function, an exponential function, and a conical function in the element-free Galerkin methodology. Graphical results of curve fitting problems show that the conical weight function exhibits lower error in comparison to the exponential weight function. The steady heat transfer analysis for orthotropic structure has been carried out by Zhang et al. [21] based on the proposed EFG technique. Different weight functions such as cubic spline, parabola function, conical function, triangle weight function, and exponential weight function have been considered. The numerical results suggest preferring cubic spline and parabola weight functions in comparison to the conical weight function for the heat transfer problems.(vi)Cosine weight function: the cosine weight function is defined as Though only a few researchers have utilized the cosine weight functions in the element-free Galerkin formulation, these weight functions provide continuous approximations. Bouillard and Suleau [22] successfully applied these weight functions for numerical assessment of the pollution effect governed by the Helmholtz model problem. After that, Singh [15] also mentioned cosine weight function along with some other weight functions for EFG formulation to solve three-dimensional heat transfer problems.(vii)Hyperbolic and elliptical weight functions: hyperbolic and elliptic weight functions are given byrespectively. Elliptical weight functions have been referenced by Singh in [15]. The author also considered exponential and cosine weight functions to carry out error and steady-state analysis for heat transfer problems under the EFG framework.
The shapes of different weight functions at any point in the domain have been presented in Figure 2.
4. Quasi-Linearization Technique
Since the problem under consideration is nonlinear in nature, and quasi-linearization is a standard technique that is widely used for approximating the nonlinear problems by their linearized version for obtaining the approximate solutions; in this section, we will discuss the quasi-linearization process and its convergence for the problem under consideration, i.e., equation (1).
Using the quasi-linearization process by Bellman and Kalaba [38], equation (1) can be linearized as follows:where stands for the approximate solution at th iteration level. For example, the nonlinear problem 3 in Section 6 can be linearized by applying the quasi-linearization process (24) as below.
The considered nonlinear differential equation in Example 3, given byis linearized as
Similarly, the nonlinear singularly perturbed differential equation in Example 4,is linearized as
4.1. Convergence of Quasi-Linearization Process
In order to prove the convergence of the quasi-linearization process, we consider
By using quasi-linearization process, the following recurrence relation provides a sequence of linear differential equations aswhere denotes the iteration level. We assume as the initial guess which also incorporates the boundary conditions. Rewriting the above relation at the previous iteration step, we obtain
Subtracting equation (32) from equation (30), we obtain
The above differential equation is of second order in . Next, with the help of Green’s function, equation (33) can be converted into an integral equation aswhere Green’s function is defined as
Now, by using mean value theorem, we obtainwhere .
Next, using equation (37) in (34), we get the following expression:
Suppose
Using maximum norm over the spatial domain and using equations (36) and (39) in (38), we obtain
Simplifying the above inequality, we obtainwhere .
Thus, the quasi-linearization process converges quadratically by choosing appropriate initial iterative approximation . Therefore, in order to find the solution of the nonlinear problem (1), it is sufficient to find the solution of its approximate linearized problem (24). Next, we will discuss the proposed numerical strategy to approximate the solution of the linearized problem (24).
5. Element-Free Galerkin Formulation
In this section, element-free Galerkin formulation has been discussed in detail.
5.1. Node Generation
The node generation has been done in a way so that more points can be put in the boundary layer region. As discussed by Roos et al. [5], when the singular perturbation parameter tends to 0, boundary layers appear in the solution, and we need to adopt some special technique to capture these boundary layers. In the present work, we have utilized Shishkin’s approach to generate more nodes in the boundary layer region and to capture these layers nicely. For boundary layers at , the distribution of nodes is carried out by dividing the spatial domain into two subintervals and , where the transition parameter is chosen aswhere is a constant chosen appropriately and denotes the total number of nodes in the spatial domain. In each of these subintervals, the nodes are generated as follows.
We definewhere the mesh spacings are given by
For 32 and 64 number of nodes and and , the node distribution is given in Figure 3. It can be easily observed that nodes condense near as becomes smaller. A similar procedure can be carried out for boundary layers appearing at .
Remark 1. Note that here Shishkin’s idea is used only to generate nodes and not the mesh as the present method does not require element’s connectivity.
Next, we will discuss the element-free Galerkin weak formulation of (24).
5.2. Weak Formulation
Using the above discussed moving least squares methodology, the element-free Galerkin (EFG) weak formulation for the linearized problem (24) is given bywhere s are the test functions generated by moving least square approach as discussed in Section 2. The last terms in equation (45) are introduced due to the method of Lagrange multipliers to enforce the essential boundary conditions (2). Substituting the EFG approximation of field function and simplifying equation (45), we obtain
On solving the above nodal weak formulations and assembling, the resulted algebraic linear equations can be written in the matrix form aswhere
and are Lagrange multipliers, and and are shape functions for node on the essential boundary.
6. Numerical Results and Discussion
In the present section, before proceeding to numerical examples, firstly, we present below the algorithm of the proposed scheme to help the interested researchers.
6.1. Algorithm for the EFG Scheme
We mention below the step-by-step algorithm for the proposed scheme:(1)Enter the given parameter values and material properties.(2)Set up the adaptive nodal coordinates in the domain of interest.(3)Set up the Gauss quadrature background points (to carry out numerical integrations), weights, and Jacobian.(4)Define the support domain of nodal points.(5)Set parameters for weight functions and calculate the weight functions for each node.(6)Define loop over Gauss integration points and do(a)Construct shape functions and their derivatives at each node.(b)Assemble the contributions from each node to matrix and the source vector based on the weak formulation as defined in the last section.(c)Incorporate the essential boundary conditions using the Lagrange multiplier method ( and matrices defined above) or using the penalty method.(7)Solve the global system using some efficient solver to get the nodal parameter values .
Next, we will consider some linear and nonlinear singularly perturbed problems to test the efficiency and robustness of the proposed scheme. In order to demonstrate the potentials of the present approach, the numerical results obtained using the proposed numerical methodology have been compared with those cited in the literature. The pointwise errors and maximum absolute errors have been presented for the considered problems. These errors have been estimated by
The accuracy of the proposed method significantly depends on the choice of the size (radius) of the support domain. Therefore, a lot of experiments have been carried out to find the range of the size of the support domain for singularly perturbed problems. It has been found that, for , the EFGM provides better numerical approximation for SPP. Computations have been performed using software on an HP machine with Intel Core processor running at speed with RAM.
Example 1. Consider the following variable coefficients’ singularly perturbed problem [39]:under boundary conditionswhere is so chosen to satisfy the exact solution:
The proposed element-free Galerkin scheme has been employed to solve the abovementioned problem. The problem has been solved by considering different weight functions and by using both the linear and the quadratic basis functions. Comparisons have been made for different weight functions in both the cases, and the results have been tabulated in Tables 1 and 2 for 128 number of nodes and different values of singular perturbation parameter . One can easily observe that the present scheme provides more accurate results by utilizing Gaussian exponential spline and cubic spline weight functions with shape parameter for both linear and quadratic basis functions. Also, as expected, the numerical results are more accurate in the case of quadratic basis functions in comparison to the linear basis functions.
The problem has been solved again for different values of and the different numbers of nodes and by considering quadratic basis and Gaussian exponential spline weight functions. The errors have been presented in Table 3. It can clearly be seen that, for each , the maximum pointwise absolute errors are decreasing as the number of nodes are increasing. This clearly depicts that the proposed scheme converges numerically. The results obtained from the proposed scheme have also been compared with those cited in the literature in Table 4. It can be noticed that the proposed scheme provides comparative better accurate results. Figure 4 also depicts that the EFG scheme is very much efficient in capturing the sharp boundary layers.
(a)
(b)
Example 2. For the second example, the following variable coefficients convection-diffusion singularly perturbed problem [39] has been considered:under boundary conditionsIts exact solution is given by
Problem 2 has boundary layers at for small values of the singular perturbation parameter . For this problem also, it can be easily noticed in both the cases, i.e., linear basis functions (Table 5) and quadratic basis functions (Table 6), that the Gaussian exponential spline weight functions provide better results in comparison to the other weight functions for all values of the singular perturbation parameter . The maximum absolute errors obtained by the proposed method for different values of and have been presented in Table 7. From all these tables, one can easily conclude the numerical convergence of the proposed scheme. Comparison of the EFGM results in Table 8 with those cited in the literature for different values of also proves the efficiency of the proposed EFG scheme. In plots (a) and (b) of Figure 5, the comparison of the EFG and the exact solutions and the errors has been plotted. The figure shows that very sharp boundary layers appear in the solution at , and the proposed EFG scheme is efficient enough to capture these layers efficiently even for very small values of . The error plot also depicts the stability of the proposed scheme.
(a)
(b)
Example 3. Consider the following nonlinear singularly perturbed problem [40]:under boundary conditions,The exact solution of the above problem is not available; however, a close approximation of the exact solution of the problem is considered to be
Alike the above test problems, comparisons of the EFG solutions obtained using different weight functions for both the linear and the quadratic basis functions have been presented in Tables 9 and 10, respectively. It can clearly be observed that, for the nonlinear singularly perturbed problem, Gaussian exponential spline weight functions perform much better in comparison to the other weight functions. The errors have been presented in Table 11 for distinct values of and the different number of nodes. The comparison of the EFG solution with that given by (58), depicting the sharp boundary layer, has been displayed in Figure 6(a). The pointwise absolute error plot in Figure 6(b) clearly demonstrates the compatibility and robustness of the proposed scheme.
(a)
(b)
Example 4. For the fourth example, we have considered the following nonlinear singularly perturbed problem [41, 42]:with uniformly valid approximation:where . From Table 12, one can analyze the validity and accuracy of the considered EFG scheme. The pointwise EFG solution of the considered problem and the solutions obtained by different methods [41, 42] have also been tabulated in this table. From the table, it is clear that the EFG results are in good agreement with those cited in the literature. The comparison of the EFG solution and the approximated solution given by equation (60) has been made in Figure 7. It can be easily seen that the EFG solution approximates the close approximation of the exact one up to a good precision of accuracy.
7. Conclusion
Meshless methods have rarely been implemented for solving singularly perturbed problems in literature. In particular, the element-free Galerkin method has not been utilized for solving SPP to the best of the authors’ knowledge. In the present work, an effort has been made to propose the appropriate weight functions in order to implement the EFG method for SPP. Since boundary layers appear in the solutions of SPP, non-equidistributed nodes have been generated which condense in the boundary layer region and thus utilize the benefit of the proposed scheme. As the choice of weight functions affects the resulting approximation in the EFG method, a variety of weight functions have been presented and tested in the EFG formulation for solving both the linear and nonlinear singularly perturbed problems. It has been concluded that the exponential weight functions are most reliable and produce the most accurate results for SPP in comparison to other weight functions in the case of both the linear and quadratic basis functions. Moving least square approximations and the Lagrange multiplier method have been utilized for generating the shape functions and to implement the boundary conditions respectively. For nonlinear problems, the quasi-linearization process has been utilized and its convergence has been discussed. Since the radius of the support domain also affects the EFG solution a lot, after a lot of experiments, it has been concluded that, for the range , the proposed method provides better approximate results. The convergence of the proposed method has been shown numerically for different problems. The EFG results have been compared with those cited in the literature and are found in good agreement for both the linear and nonlinear singularly perturbed problems.
Data Availability
The datasets supporting the results presented in the study have been generated using the MATLAB code.
Conflicts of Interest
The authors declare no conflicts of interest in this study.
Acknowledgments
The authors take this opportunity to acknowledge the use of Lab facilities provided by the School of Mathematics, Thapar Institute of Engineering and Technology, Patiala, due to the DST-FIST grant SR/FST/MS-1/2017/13.