Error Estimates for a Multidimensional Meshfree Galerkin Method with Diffuse Derivatives and Stabilization

A meshfree method with diffuse derivatives and a penalty stabilization is developed. An error analysis for the approximation of the solution of a general elliptic differential equation, in several dimensions, with Neumann boundary conditions is provided. Theoretical and numerical results show that the approximation error and the convergence rate are better than the diffuse element method.


Introduction
Numerical methods based on moving least square (MLS) approximations and Galerkin formulations form a popular class of meshfree schemes.However, the high computational expense in the evaluation of the shape functions and their derivatives are drawbacks to the Galerkin approach.For this purpose Belytschko et al. [1] and Breitkopf et al. [2] have introduced efficient computational approaches for the evaluation of the MLS shape functions and their derivatives.
An alternative for the computation of derivatives, the diffuse derivative, was used by Nayroles in [3] in the DEM.In the diffuse derivative approximation, only the derivatives of the polynomial basis need to be included in computing the gradients of the local field variables.Belytshko et al. [4], [5] argued that diffuse derivatives are not attractive in Galerkin methods because they degrade the accuracy due to their lack of integrability.However, recently, the diffuse derivative has been used in a class of novel meshfree methods (Huerta et al [6]) for Stokes problems.Because of their simplicity, diffuse derivatives, unlike the full derivatives, retain the same subspace structure as their defining functions.This special feature allowed Huerta et al [6] to circumvent the complicated

|54
Ingeniería y Ciencia incompressibility constraint and define a class of divergence free meshfree approximation functions.Beyond fluid mechanics, we think this new approach could be used to enhance the common mixed method approach.
There are very few papers on meshfree methods with a complete error analysis [7]; among them we note [8] on the RKPM, [9] on the EFG and [10] on the MPCM.The last paper provides a complete mathematical analysis of the MPCM with diffuse derivatives, applied to a Poisson problem with Dirichlet boundary conditions.However it was not until our previous work in [11] that a paper with a complete error analysis of a Galerkin meshfree method with diffuse derivatives (DEM) was done.
In this paper, we introduce an extension of our previous work in [11] where we introduced a Galerkin MLS scheme developed entirely with diffuse derivatives, that we called Stabilized Diffuse Galerkin Method (SDGM).That previous work was completely applied to a 1D problem.We used a novel stabilization procedure and, unlike any of the previous work on diffuse derivative schemes for differential equations, provided a full error analysis of the new approach as well as example computations.The new scheme, when applied to self-adjoint elliptic problems, leaded to fully invertible symmetric positive matrices and has rates of convergence that improve as polynomial degree m is increased.Now in this paper we provided the necesary modificatios to extend these results to a general elliptic differential equation in any dimensions.As we will show in sections 3, 4 and 5, additional consideretions have to be made.
The use of stabilization (or penalty) terms is not uncommon; they have been used in finite element methods [12], [13] and have been proved to be quite effective and, indeed, they are often introduced intuitively.For instance, Beissel and Belytschko [14] introduce a penalty term to stabilize nodal integration in the EFG.We will show that our SDGM is accurate and its rate of convergence increases as the approximating polynomial degree increases and the width of the support domain (R) decreases.
This paper is organized as follows.In section 2 we introduce the model probles.In section 3 we review some aspects of MLS approximations.In section 4 we describe the diffuse derivatives and prove an approximation theorem.In section 5, we introduce the SDGM and prove our main error estimate.Section 6 gives numerical results and provides a comparison with the theoretical results.

Model Problem
For a bounded domain Ω with C 1 -boundary ∂Ω we consider problems of the form: where and ν is a unit normal vector to ∂Ω.
Here c(x) ≥ c 1 > 0 for all x ∈ Ω and A(x) = (a ij (x)) n×n is assumed to be uniformly elliptic in Ω, i.e. there exists θ > 0 such that for all x ∈ Ω and all η ∈ R n , After multiplying this differential equation by function v and using Green's formulas we find that: ) dx and

|56
Ingeniería y Ciencia And this problem has a unique solution by the Lax-Milgram theorem.
We have used Neumann or flux boundary conditions here for simplicity; they are natural and do not need to be imposed explicitly.We expect that Dirichlet conditions could be implemented in this framework using, for instance, the Nitsche approach (See, perhaps, [6]).But, we do not pursue that here in order to focus on this new diffuse derivative approach.
We, further, expect that this approach could be extended to problems with higher order derivatives where the payoff in computational cost from the use of diffuse derivatives instead of full derivatives would be more pronounced.
3 Preliminaries on moving least squares (MLS)

The moving least squares method
Let Λ = {x 1 , x 2 , ..., x N } be a set of N distinct points inside and in the boundary of Ω ⊂ R n which is an open and bounded set with Lipschitz boundary ∂Ω and u 1 , u 2 , ...u N be the values of an unknown scalar function u(x) at the points in Λ (i.e.u i = u(x i ), 1 ≤ i ≤ N ).Also let R > 0 (usually called dilation parameter) and consider a positive even weight function W (x) with compact support in B 1 (0) and Let m << N and p(z) = {p 0 (z), p 1 (z), p 2 (z), ..., p s (z)} T be a basis of the subspace of polynomials of degree less or equal than m (denoted P m ) in R n , placed in multi-index ordering.Note that s+1 = (n+m)!/(n!m!).For each x ∈ Ω consider where a(x) = {a 0 (x), a 1 (x), a 2 (x), ..., a s (x)} T is chosen such that it ing.cienc., vol.9, no.17, pp.53-76, enero-junio.2013.
minimizes the functional for a fixed x.
For each fixed x, P u (x, y) is a polynomial in y of degree less than or equal to m that represents the best local approximation of u(x) for y in a small neighborhood of width 2R of x.Then, since the weight function W usually favors the points closer to x it is natural to define the following approximation of u(x): A short calculation ( [15]) shows that where Thus, This expression can be written in the standard interpolation form where φ i (x), 1 ≤ i ≤ N, are called the shape functions and are given by and u R (x) corresponds to the moving least squares approximation of the function u at the point x ∈ Ω.Note also that from the minimization of functional J(a(x)) and, in computing the minimum via partial derivatives, the following orthogonality relationship holds: where q = q(z) is a polynomial of degree less than or equal to m.

Some error estimates for MLS approximations
We use the notation W m,p (Ω) to denote the Sobolev space consisting of functions with m derivatives in , and H m (Ω) for the special case where p = 2.We use the following notation for the norms and seminorms: We will also denote ∥v∥ := ∥v∥ L 2 (Ω) .
The moving least squares method (MLS) has been used for the numerical solution of ordinary and partial differential equations in many papers, but very few of these include error estimates.The papers [9] and [8] provide such analysis and require the following properties: (i) For each x ∈ Ω there exists at least s + 1 distinct points from Λ in BR 2 (x), where Lagrange interpolation is possible.
(ii) There exists (iii) There exists c # such that for all x ∈ Ω, (iv) For any x ∈ Ω there exists a constant c L such that the Lagrange basis functions associated with the set the points in property (i) are bounded by c L in B 2R (x).
Note, in particular, that c 0 , c # , c L and c 1 are independent of R and N .Here (iii) implies that card(Λ(x)) ≤ c # .Throughout the rest of this paper, C will denote a positive constant that is independent of R (and thus N as well).
Remark 3.1.In [8], Han and Meng stated that if conditions (ii) and (iv) are satisfied, then the family of particle distributions, Λ, is called regular, and it is enough to guarantee that there exists a constant L 0 such that With these assumptions the following theorem was established (see proof in [9], [16], [8], [17], [15]):

The diffuse derivative
Computation of the derivatives of a MLS function involves differentiation of a(x) which, in turn, involves differentiation of the M −1 and B matrices.
On the other hand, the derivative of polynomials in P m is trivial and can be evaluated a priori.The concept of diffuse derivative, proposed in [3], exploits the fact that, for a MLS function u, P u (x, y) is a good approximation of u = u(y) near the point x.Thus, in the onedimensional case: where ) .
With obvious modifications for the multidimensional case, using multi-index notation.Below we indicate why u ′ (x) ≈ δu(x).It can be shown (see [18]): where Ω is a bounded open set in R n with Lipschitz boundary, and suppose (where D β v and δ β v represent the full and diffuse derivatives, of order |β| of v in multi-index notation).

Now consider a function
where v i , 1 ≤ i ≤ N, are constants and φ i (x) are the MLS shape functions.Note that V can also been written as with The analysis in [9] provides arguments to show that derivatives of the a(x) functions are small, which is crucial in our theory.The next lemma shows how the diffuse derivative is controlled by differences with the local MLS functional P V (x, y).This lemma, somewhat of an inverse estimate, provides the key idea behind our stabilization. where Proof.Our proof follows very closely the proof of lemma 2.2 in [9] and Lemma 3 in [11].
Let also, e j = jth unit vector in R n .A short calculation shows that So, ] .
Notice that for each x, Q(x, y) is a polynomial in y of degree less or equal than m.Also, if z 1 , z 2 , ..., z s+1 are the points in B R/2 (x) given by property (i), then using property (ii), By the orthogonality property (7) the second term on the right hand side above is zero.Thus ing.cienc., vol.9, no.17, pp.53-76, enero-junio.2013.
By property (iv) and the Meanvalue theorem, we can guarantee that, for h small enough, ∃ξ l such that and using again the orthogonality condition (in x + he j ) and property (v), we have where On the other hand, we can use the Lagrange's polynomial basis functions l 0 , l 1 , . . ., l s , each of degree m at the s + 1 distinct points in B R/2 (x) to write So, using the Cauchy-Schwarz inequality, Using this result and Cauchy-Schwarz inequality in equation (12) we have and thus, Using the equation ( 14), choosing z = y + he j and property (iv) in ( 13) we obtain, for all y ∈ B R (x) ∩ Ω, Finally, using this result in the equation ( 11) yields: Remark 4.1.In the one dimensional case we presented in [15] and [11], we had that if V = u R (i.e.v i = u i = u(x i )) then P u (x, •) is an interpolating polynomial of the function u at least on m + 1 distinct ing.cienc., vol.9, no.17, pp.53-76, enero-junio.2013.
points in Λ(x); however, in the multidimensional case we are not aware of a similar result or whether this is true or not.Therefore, let us introduce for any x ∈ Ω the polynomial I u (x, •) in P m that interpolates u at the points that satisfy the property (i).We can prove the following theorem (see [9] or [15]): Therefore, in particular taking y = x we have

A Galerkin Approximation Scheme
The stabilization of our Galerkin scheme will involve where U and V are MLS functions in V R as defined in (9).So, by lemma 4.1 we have Let us also notice that from theorem 4.1 we have Thus, by standard arguments from polynomial interpolation

|66
Ingeniería y Ciencia Hence, For MLS functions u and w, the bilinear form for our Galerkin scheme is and A is a matrix whose ij entry is a ij So, our stabilized diffuse Galerkin method (SDGM) is as follows: SDGM: Here, we require γ > 0 and guidelines for it will be provided in the next theorem.
from which we can identify R −2γ P(V, V ), as a penalty or stabilization term.
We can now state and prove our main theorem: Theorem 5.1.Let u ∈ C m+1 (Ω) be the exact solution of the equation ( 1), u R be its MLS approximation and U be the solution given by the numerical scheme SDGM, then if γ = m/2 + 1, there exists a constant C independent of R such that, where Proof.The proof of this theorem is similar to the proof of theorem 1 in [11].

Numerical results
In this section we present some numerical results on the convergence orders of the SDGM.The numerical results confirm the theoretical predictions.
Solutions are reported for the following numerical methods: 1. Diffuse element method (DEM).

The stabilized diffuse Galerkin method (SDGM).
In all cases a background mesh of subintervals on cells was used for numerical integration.Within each integration cell, there was a set of Gauss-Legendre quadrature points.We kept the number of cells large enough so that numerical integration did not affect the convergence rates.The weight function was chosen to be the cubic spline:

|70
Ingeniería y Ciencia Figures 1 and 2 show numerical results comparing EFG, DEM and SDGM for different dimensions (m) of the polynomial basis used in the MLS approximation.We report errors for both the numerical approximation of the solution of the equation ( 27) and the approximation of its first derivative with respect to x using diffuse derivatives, in the L 2 norm.It is important to notice that similar results can be obtained for the full and diffuse derivatives with respect to y, and therefore for the gradient and diffuse gradient of u.The convergence rates are summarized in tables 1 and 2. It can be seen that the numerical results suggest that for the SDGM while we only get about ∥u − U ∥ ≤ CR 1.5 and |δ x u − δ U ∥ ≤ CR 0.5 for

|72
Ingeniería y Ciencia the DEM.The numerical convergence rates are slightly better than our theoretical results on ||δ x u − δ x U || and ||u − U ||, as is often true.We also observe that, in general, the errors obtained using SDGM are smaller than those using the standard DEM.The convergence rates for the DEM in the L 2 norm are about 1.5 independent of the value of m and 1 for its derivative (similar observations were made in [19] and [11], although a small convergence rate for the DEM was obtained in the 1D case).Our SDGM performs better as m increases.Also, as expected the standard EFG gives the best convergence rates, but unfortunately it is very expensive from the computational point of view, which makes diffuse derivatives more attractive.It is important to notice that all these results agree with those found for the 1D case in [11].Remark 6.1.A proof of the enhanced L 2 convergence observed in the equation ( 28) remains an open question.The fact that our penalty term is not exactly zero at the true solution and integration by parts with diffuse derivatives is not possible, seem to make the standard duality argument impossible.

Conclusions
A modification to the traditional DEM, the stabilized diffuse Galerkin method (SDGM) has been proposed, in which a stabilization term is introduced to improve the overall accuracy and stability.The new scheme, like DEM, does not require the evaluation of full derivatives.This method is shown to give better results than DEM and converges to the true solution as the dilation parameter (R) goes to zero, or the order of the polynomial basis is increased, as demonstrated numerically and proved theoretically.The procedure described in this paper can be applied to more general multidimensional problems (see [15]).
We see SDGM as enhancing the viability of the diffuse derivative approach.Again, as suggested in Huerta et al [6], we think the versatility of the diffuse derivative could be helpful in fluid flows or mixed method computations.

Figure 2 :
Figure 2: Numerical solution error in L 2 norm.The continuous line with diamonds is for the SDGM, the dashed line with squares for the DEM and the pointed line with circles for the EFG.

Table 1 :
Convergence Rates for Solution Approximation in L 2

Table 2 :
Convergence Rates for Diffuse Derivative Approximation (with respect to x) in L 2