Exponentially Fitted Element-Free Galerkin Approach for Nonlinear Singularly Perturbed Problems

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.*e 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 L∞ errors have been presented. Comparison of the EFG solutions with those available in the literature depicts the superiority of the proposed scheme.


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. e 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 reactiondiffusion 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. erefore, 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.
e authors [11] devised a modified reproducing kernel method (RKM) to solve singularly perturbed delay boundary value problems later in 2015. e 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 by where ϵ is a small positive parameter known as singular perturbation parameter with 0 < ϵ ≪ 1, and we suppose that f(v ′ , v, x) is continuously differentiable function. ese 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 elementfree Galerkin (EFG) method, one of the most popular and versatile mesh-free methods, for approximating nonlinear SPP solutions. e EFG technique is based on the theory of moving least squares (MLS) approximation and standard weight functions to construct the shape functions. e 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. ese 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. e 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. e paper is organized as follows. e 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.

Moving Least Squares Approximation
e field function v(x) is approximated by v h (x) using MLS approximants in the domain Ω. ese 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. e domain of influence of any nodal point x I is defined as a local neighborhood of x I or the domain that the particular node x I exerts an influence upon. e MLS's approximation of the unknown function v h (x) is defined as follows: where c(x) is a vector of monomial basis functions, m is the number of components in the basis vector c, and y(x) is an unknown vector of coefficients of the basis functions. e choice of basis functions c(x) is essential in the construction of moving least squares based shape functions. e following are the examples of the basic functions: (i) Linear basis: (ii) Quadratic basis: e field function's approximate value at any particular node x I in the support domain of x is defined by We define the weighted residual functional as follows to determine the unknown coefficient vector y(x) in the support domain of the point of interest x: where n denotes the number of nodes x I in the support domain of the point of interest x. e nodes in the support domain of any point x are used to support or approximate where α s signifies the dimensionless size of the support domain and the characteristic length d c denotes the nodal spacing for the point x. w(x − x I ) in expression (7) is some weight function which is generally chosen as monotonically decreasing as |x − x I | increases. In the literature, numerous weight functions have been introduced and employed by researchers, as mentioned in Section 3. In MLS approximation, the condition m > n generally leads to the ill-conditioned system of equations. On simplification, equation (7) can be rewritten as where and W is the diagonal matrix defined as e minimization of J with respect to y(x) leads to is gives where erefore, from (3), the approximated field function v h can be computed as where ϕ I (x), I � 1, 2, 3, . . . , n, are the shape functions defined by e graphs of typical shape functions for different nodes have been plotted in Figure 1. We can observe that the shape functions are bell-shaped. e derivatives of the MLS shape functions can be obtained by differentiating (16) with respect to the spatial coordinates.

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. e weight functions must be continuous and positive in the corresponding support domains. e weight functions should satisfy the following properties: represents the Dirac delta function and h is measure of the size of the support domain Ω x . e weight functions are usually based on a normalized distance s � (|x − x j |/d s ), where d s denotes the radius of the influence domain. In [15][16][17][18][19][20][21][22][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.
ough 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 d s � α s * d c , where the scaling parameter d c satisfies 0 < d c < 1.5, whereas exponential weight functions provide good results for the scaling parameter 1 < d c < 2.0. 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. e authors constructed cubic spline weight functions to impose C 2 continuity requirements while solving linear elastostatic problems. ese are some of the most frequent weight functions which are widely adopted by many researchers while employing the elementfree 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 ese weight functions are also widely used by researchers as they provide C 2 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. e 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 advectiondiffusion [29]. e impact of cubic and quartic spline weight functions on computational precision has been discussed in both the articles. e proposed method has higher computational accuracy while employing quartic spline weight functions. e 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 of exponential type is one of the most commonly used weight functions. ese 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]. e 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. e 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. e graphical results indicate that error estimation of conical formula increases by increasing the value of scaling parameter d c . Zhuang et al. [35] highlighted issues and sources of errors for meshless methods. e 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. e 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. e 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 ough only a few researchers have utilized the cosine weight functions in the element-free Galerkin formulation, these weight functions provide C 1 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 by respectively. Elliptical weight functions have been referenced by Singh in [15]. e author also considered exponential and cosine weight functions to carry out error and steady-state analysis for 3D heat transfer problems under the EFG framework. e shapes of different weight functions at any point x in the domain have been presented in Figure 2.

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 v (k+1) stands for the approximate solution at (k + 1)th iteration level. For example, the nonlinear problem 3 in Section 6 can be linearized by applying the quasi-linearization process (24) as below. e considered nonlinear differential equation in Example 3, given by is linearized as Similarly, the nonlinear singularly perturbed differential equation in Example 4,

Journal of Mathematics
is linearized as

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 as where (k) denotes the iteration level. We assume v (k�0) 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 e above differential equation is of second order in (v (k+1) − v (k) ). Next, with the help of Green's function, equation (33) can be converted into an integral equation as where Green's function G(x, r) is defined as  Journal of Mathematics Now, by using mean value theorem, we obtain where Next, using equation (37) in (34), we get the following expression: Using maximum norm over the spatial domain and using equations (36) and (39) in (38) Simplifying the above inequality, we obtain where Q � L/(8ϵ − 2R). us, the quasi-linearization process converges quadratically by choosing appropriate initial iterative approximation v (k�a) . erefore, 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).

Element-Free Galerkin Formulation
In this section, element-free Galerkin formulation has been discussed in detail.

Node Generation.
e 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 x � a, the distribution of nodes is carried out by dividing the spatial domain Ω � [a, b] into two subintervals [a, a + δ] and [a + δ, b], where the transition parameter δ is chosen as where M is a constant chosen appropriately and (N + 1) denotes the total number of nodes in the spatial domain. In each of these subintervals, the nodes are generated as follows.
We define Journal of Mathematics where the mesh spacings h i ′ s are given by For 32 and 64 number of nodes and ϵ � 2 − 5 and ϵ � 2 − 10 , the node distribution is given in Figure 3. It can be easily observed that nodes condense near x � a as ϵ becomes smaller. A similar procedure can be carried out for boundary layers appearing at x � 1.

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).

Weak Formulation.
Using the above discussed moving least squares methodology, the element-free Galerkin (EFG) weak formulation for the linearized problem (24) is given by where ϕ i s are the test functions generated by moving least square approach as discussed in Section 2. e 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 v(x) and simplifying equation (45), we obtain 8

Journal of Mathematics
On solving the above nodal weak formulations and assembling, the resulted algebraic linear equations can be written in the matrix form as where λ � n l�1 N l λ l and ρ � n l�1 R l ρ l are Lagrange multipliers, and N I and R I are shape functions for I th node on the essential boundary. Journal of Mathematics

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.

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. (7) Solve the global system using some efficient solver to get the nodal parameter values v i .
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. e pointwise errors and maximum absolute errors have been presented for the considered problems. ese errors have been estimated by e accuracy of the proposed method significantly depends on the choice of the size (radius) of the support domain. erefore, 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 2.8 ≤ α c ≤ 3.4, the EFGM provides better numerical approximation for SPP. Computations have been performed using MATLABR2016a software on an HP machine with Intel Core i3 processor running at 1.70 GHz speed with 4 GB RAM.
Example 1. Consider the following variable coefficients' singularly perturbed problem [39]: where f(x) is so chosen to satisfy the exact solution: v e proposed element-free Galerkin scheme has been employed to solve the abovementioned problem. e 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 α � 0.3 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. e 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. e L ∞ 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. is clearly depicts that the proposed scheme converges numerically. e 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.
Example 2. For the second example, the following variable coefficients convection-diffusion singularly perturbed problem [39] has been considered: Problem 2 has boundary layers at x � 0 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 ϵ. e maximum absolute errors obtained by the proposed method for different values of ϵ and N 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    Figure 5, the comparison of the EFG and the exact solutions and the L ∞ errors has been plotted. e figure shows that very sharp boundary layers appear in the solution at x � 0, and the proposed EFG scheme is efficient enough to capture these layers efficiently even for very small values of ϵ. e L ∞ error plot also depicts the stability of the proposed scheme.
Example 3. Consider the following nonlinear singularly perturbed problem [40]: e exact solution of the above problem is not available; however, a close approximation of the exact solution of the problem is considered to be v(x) � − ln 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. e L ∞ errors have been presented in Table 11 for distinct values of ϵ and the different number of nodes. e comparison of the EFG solution with that given by (58), depicting the sharp boundary layer, has been displayed in Figure 6(a). e pointwise absolute error plot in Figure 6(b) clearly demonstrates the compatibility and robustness of the proposed scheme.
Example 4. For the fourth example, we have considered the following nonlinear singularly perturbed problem [41,42]: with uniformly valid approximation:   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.

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

16
Journal of Mathematics 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 α s ∈ (2.8, 3.2), the proposed method provides better approximate results. e convergence of the proposed method has been shown numerically for different problems. e 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 e datasets supporting the results presented in the study have been generated using the MATLAB code.

Conflicts of Interest
e authors declare no conflicts of interest in this study.