Self-consistent boundary conditions with blurred derivatives

Engineering Analysis with Boundary Elements 29 (2005) 326–333 www.elsevier.com/locate/enganabound

Self-consistent boundary conditions with blurred derivatives Enrique Pardo* Facultad de Ingenierı´a, Universidad de Mar del Plata, Juan B. Justo 4302, Mar del Plata 7600, Argentina Received 13 February 2004; accepted 26 August 2004 Available online 10 March 2005

Abstract It is shown in this paper that self-consistent boundary conditions for numerical methods based on blurred derivatives can be derived from a suitable change of variables of the fundamental blurred approximation of the differential equation, followed by application of Leibnitz theorem for differentiation of an integral. The simplest scheme obtained in this way resembles the weak Local Petrov-Galerkin approximation, although interpretation of the operators appearing in the final equations is quite different—as is the derivation itself. Subsequent transformation leads to integral equations similar to the starting point for boundary integral methods of solution. In this way, a number of well-known computational methods are shown to be derivable from adequate manipulation of the blurred derivative technique. However, other approximations, which are not derivable with standard methods can also be obtained, hinting at a greater generality of blurred derivatives. q 2005 Published by Elsevier Ltd. Keywords: Meshless methods; Local weak forms; Blurred derivatives

1. Introduction Even though a large number of so-called meshless methods have been devised in the last decade [1–17], development of an effective mesh-free numerical solution method still remains a challenge. Actually, this has been identified by researchers as one of the eight critical research tasks facing the field of computational mechanics [18]. In fact, most meshless methods are not clearly advantageous over the robust, versatile and well-understood finite element method. Methods based on global principles, such as the Galerkin approximation over he whole domain under consideration or variational principles, still require some kind of background mesh to perform integrations, so that they are not truly meshless. In contrast, methods based on strong formulations, such as point collocation, are usually less accurate than the former. Hence, it is reasonable to infer that an effective meshless method should meet two requirements: it must have a local character—in order to avoid global operations on the whole domain—and must be

* Tel.: C54 223 481 6600; fax: C54 223 481 0046. E-mail address: [email protected] 0955-7997/$ - see front matter q 2005 Published by Elsevier Ltd. doi:10.1016/j.enganabound.2004.08.009

devised to avoid high order differentiation of the approximating functions. One method that satisfies both conditions is the Local Petrov-Galerkin Method (LPGM) [16,17]. It starts from the usual weighted residual statement but, unlike finite elements, enforces the integral identity on each weight function’s support separately thus eliminating ‘global’ operations. This in turn requires the weight functions to be zero over their boundary in order to eliminate boundary terms. For boundary nodes, only the intersection of the global boundary with the weight function contributes to the boundary integrals of the weak principle. A completely unrelated method which also has a strictly local character is the functional integral method (FIM). Although this method can be introduced as an ad-hoc procedure [5–7], it is most conveniently derived from the formalism of blurred derivatives [15]. The main drawback with such method is that boundary conditions must be enforced through some artifice, not emerging from the formulation itself. In this paper it is shown, however, that boundary conditions also emerge in FIM in a natural manner, through suitable transformations applied to the fundamental blurred approximation. In doing so, unexpected relationships between FIM, LPGM and other methods also emerge.

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

2. Blurred derivatives Briefly, the technique of blurred derivates [15] starts by rewriting the value of any continuous function at a point as: ð
x
f ðxÞ Z lim P0 ðx
K x; 3Þf ðxÞd (1) 3/0

U0

327

The concept of blurred derivative can be used to obtain a number of numerical schemes. The simplest one is to replace actual derivatives in the differential equations by blurred derivatives. As an example, let us consider the onedimensional Poisson equation: v2 uðxÞ C gðxÞ Z 0 vx2

(7)

The function P0 in Eq. (1) can be any function that yields the Dirac delta function in the limit 3/0, and U0 is its support. Hence, it must satisfy the condition: Ð 0 PðxKx;3Þ dx
Z 1 (2)


Using the blurred derivative kernel of order 2 we rewrite Eq. (7) as: Ð Ð 2
dx
C P0ðxKx;3Þ
dx
Z 0 uðxÞ gðxÞ (8) PðxKx;3Þ



The simplest and most versatile choice is a Gaussian:

If the Gaussian family of kernels is used Eq. (8) takes the specific form:

2

0


KðxKxÞ 2

e

3

P ðx
K x; 3Þ Z pffiffiffi p3

(3)

From Eq. (1) it follows that the derivatives of f(x) can be computed as: ðN n 0 dn f ðxÞ d P ðx
K x; 3Þ
x
Z lim f ðxÞd 3/0 KN dxn dxn ðN
x
Z lim Pn ðx
K x; 3Þf ðxÞd (4)

2
 

KðxKxÞ 1 e 32 x
K x 2 p ffiffiffi ffi
dx
K 2 uðxÞ 4 2 3 p3 KN 3 2
ðN KðxKxÞ e 32 pffiffiffiffi gðxÞ
dx
Z 0 C p3 KN

ðN

To obtain a discrete scheme the exact continuous field ~ u(x) is replaced by an approximation uðxÞ in terms of N unknowns (i.e. nodal values):

3/0 KN

n

The functions P are called Blurred derivative kernels. After integration they transform a given function, f(x), into a new one, f31 , which in the limit 3/0 yields the derivative of the former one: df ðxÞ hf 1 ðxÞ Z lim f31 ðxÞ 3/0 dx

(5)

If the Gaussian (3) is chosen as the zero order kernel, each derivative kernel turns out to be just the product of P0 by the corresponding Hermite polynomial Hn(x) apart from a factor 3Kn. Details on this respect can be found elsewhere [15]. The main drawback of Gaussian kernels is that their support is infinite, so that they extend beyond the domain where the differential equation to be approximated is defined. However, it has been shown elsewhere [5] that the optimum choice of parameter 3 is such that the Gaussian becomes negligible beyond the closest neighbor to a node, so that they give excellent results in practice. But it is possible to choose a number of kernels of finite support. The simplest, of course, is a truncated Gaussian. Another very convenient choice are functions of the form: 8  n ðx
K xÞ2 < 1 G½1=ðnC3Þ if jx
K xj! 3 P0ðxKx;3Þ Z 3 G½1=2 G½nC1 1 K 32
: 0 otherwise (6) The exponent ‘n’ in this formula can be chosen so that derivatives at the support’s boundary are zero up to the desired order.

(9)

~ Z uðxÞ

N X

4i u i

(10)

iZ1

Replacing (10) in (9) and evaluating the later at each node leads to a system of N equations which are solved for the unknowns. A more detailed description of discretization will be given in Section 6. This scheme is easily generalized to several space dimensions. As an example, the bidimensional Gaussian kernel of order zero is:
KðrKrÞ 2

0

P ðr
K r; 3Þ Z

e

2

3

p32

;


yÞ
t r Z ðx; yÞt r
Z ðx;

(11)

Integrations as in formulas (1) to (4) are now carried over both space variables. Higher order kernels are computed by just taking the corresponding partial derivatives of P0. Hence, the kernel corresponding to the gradient is in this case the vector:

1 0 1
P ðr K r; 3Þ Z (12) P ðr
K r; 3Þ2ðr
K rÞ 3 And the kernel for the laplacian is: 2 2 P0 ðr
K r; 3Þ½ðr
K rÞ2 K 1 P2 ðr
K r; 3Þ Z 3

(13)

In this way, the gradient and laplacian of any continuous function f(r) are computed as: Ð
V2 f ðrÞ ~ f ðrÞ Z lim P1 ðr
K r; 3Þf ðrÞd
U; V 3/0

Ð
dU
Z lim3/0 P2 ðr
K r; 3Þf ðrÞ

(14)

328

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

r of a fixed shape so that it will be written as U0(r) to make the dependence on r explicit. Hence, in order to evaluate the partial derivative of f(r) theorem (15) must be used yielding:

3. Self-consistent boundary conditions with blurred derivatives It should be noticed that the scheme described above is only applicable to interior points since boundary conditions are not taken into account. Within this scheme Dirichlet conditions can be enforced in the usual manner (i.e. penalty methods, Lagrange multipliers, etc.), while Neumann conditions must be enforced by some ad-hoc procedure, such as the use of special shape functions for boundary nodes. In contrast, weak formulation naturally include Neumann conditions through suitable boundary integrals that arise from integration by parts. It will be shown in this section that boundary conditions also arise in a natural manner from a suitable manipulation of the blurred derivative technique, yielding a self-consistent scheme. For this, we first notice that if blurred derivatives as defined in Eq. (3) are applied to points on the boundary, only the portion of P0 that lies entirely inside the domain should be considered. This in turn requires that the kernel should be properly scaled to ensure condition (2). For such points differentiation of the blurred approximation (1) includes a boundary integral term which stems from the non-zero value of P0 at the boundary. This can be formalized as follows. Let F be a function of two sets of spatial variables, r
and r, and let the domain of integration depend on r. Differentiation of an integral of such function satisfies: ð ð v
rÞdV
Z
rÞdV
C
rÞni dS
Fðr; vxi Fðr; Fðr; vxi

%

UðrÞ

UðrÞ

vUðrÞ


V
P0ðrKr;3Þ f ðrÞd


U0 ðrÞ

ð Z

i
V
C PxðrKr;3Þ f ðrÞd


0
i dS
f ðrÞn
ðrKr;3Þ

(16)

where Pxi is the i-th partial derivative kernel. Notice that theorem (15) is unnecessary when the integration domain is fixed shape not dependent on the differentiation variable, as in this case the boundary integral vanishes identically. This is a crucial point for the derivation that will be given below, so it is worth emphasizing it. The fundamental approximation of blurred derivatives, Eq. (1), can be expressed in two different ways in the following manner: ð


V
Z P0ðrKr;3Þ f ðrÞd


U0 ðrÞ

ð


P0ðr;3Þ
f ðr
C rÞdV

(17)

U0

In the left hand side of Eq. (17) the kernel P0 is translated to point r so that its differentiation requires the use of theorem (15) yielding Eq. (16). However, the partial derivative of the right hand side of identity (17) does not contain a boundary integral because the domain of integration, U0, is now fixed since the translation is applied to the function f(r). It is simply given by: ð v xi


P0ðr;3Þ
f ðr
C rÞdV

ð Z

U0


P0ðr;3Þ
vxi f ðr
C rÞdV

(18)

U0

An important point to consider is that kernels of finite support can be chosen so that the boundary term in Eq. (16) vanishes for any interior point as long as the support U0 does not intersect the boundary of the global domain U as depicted in Fig. 1. This is the case for instance of the onedimensional polynomial kernel (6) and its multi-dimensional analogs. However, for a point on the boundary of the global domain there will be a contribution to the last term in Eq. (16) corresponding to the intersection of U0 with the boundary of the global domain. This region is indicated as G0 in Fig. 1. We can now apply the preceding ideas to Poisson equation, starting from a blurred approximation of it: V2 uðrÞ C gðrÞ Z 00 ð

Fig. 1. Schematic domain showing kernel’s support for interior and boundary nodes.

%P vU0 ðrÞ

U0 ðrÞ

(15) where ni is the i-th component of the outward normal. Eq. (15) is a generalization to several space dimensions of Leibnitz theorem for differentiation of an integral. It can be used to evaluate blurred derivatives as follows. Consider the blurred approximation—Eq. (1)—of a function f(r). The integral must be carried over the support of the kernel, U0, as depicted in Fig. 1. This support is the translation to point

ð

vxi f3 ðrÞ Z vxi

C U0

ð

2
P0ðr;3Þ
V uðr
C rÞdV

U0


P0ðr;3Þ
gðr
C rÞdV Z 0

(19)

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

The integral of the Laplacian in the Eq. (19) can be transformed into: ð ð 0 2

Pðr;3Þ P0ðr;3Þ
Vx uðr
C rÞdV Z
vxi vxi uðr
C rÞdV U0

ð Z vx i

U0


V
Z vxi P0ðr;3Þ
vxi uðr
C rÞd

ð


V
P0ðrKr;3Þ vx
i uðrÞd


U0 ðrÞ

U0

(20) Notice that in the last equality of (20) we have made use of identity (17), so that the partial derivative of the function u inside the last integral had to be changed from xi to x
i . In doing so, the domain of integration now depends on the differentiation variable x. Hence, theorem (15) can be applied to obtain: ð ð 0

VZ
V
vx i PðrKr;3Þ vx
i uðrÞd vxi P0ðrKr;3Þ vx
i uðrÞd

U0 ðrÞ

C

U0 ðrÞ

%

(21)


i dS
P0ðrKr;3Þ vx
i uðrÞn


vU0 ðrÞ

So that the approximation of Poisson equation reads: ð
V
C
i dS
vxi P0ðrKr;3Þ vx
i uðrÞd P0ðrKr;3Þ vx
i uðrÞn



%

U0 ðrÞ

ð

C

vU0 ðrÞ


V
P0ðrKr;3Þ gðrÞd


4. Relation with other methods

Z0

The first term of Eq. (22) contains differentiations over two different sets of variables: variable ‘r’ that is used to shift the origin of the kernel, and the variable with an overbar which is the integration variable. But it is possible to switch the differentiation variable in the first factor of the leftmost integral of (22) from r
to r. In doing so, the first term changes sign and a single ‘type’ of differentiation— namely, with respect to barred variables—applies everywhere. Using compact notation the equation reads: ð
V
C
S
K VP0ðrKr;3Þ VuðrÞd P0ðrKr;3Þ qðrÞd



%

ð

C

vU0 ðrÞ

(23)


V
Z 0 P0ðrKr;3Þ gðrÞd


V2 uðrÞ C gðrÞ Z 0

u Z u0 on Gu

vu hq Z q0 on Gq vn

(25)

(26)

The standard procedure is to replace the differential statement (25) by a weighted residual statement: ½V2 uðrÞ C gðrÞwðrÞdv Z 0

(27)

U

where, as usual, we have called vu Zq vn

Except for the use of the kernel P0 and the corresponding integration over its support, Eq. (23) is a familiar one in computational mechanics. Indeed, a very similar integral identity is the cornerstone of several numerical methods based on weak formulations. To fix ideas let us consider the solution of Poisson equation on a domain U:

ð

U0 ðrÞ

Vu:n^ Z

equation on a domain U, N nodes with coordinates {r1,r2,.,rN} are placed in the domain and Eq. (23) is formulated for each one (i.e. taking rZr1,r2,.,rN in (23)). After approximating the continuous field u(r) in terms of N nodal unknowns {u1,u2,.,uN}, this leads to a set of N equations that are solved for the N unknowns. Clearly, the set of supports U0 do not define a partition of the global domain U. Consider first kernels of compact support such as the polynomials given by Eq. (6) and its multidimensional analogs. Their supports are just a set of circles (or spheres in 3D) for internal nodes. As stated above derivatives of P0 up to the desired order can be made to vanish over the border of its support by a proper choice. Hence, the boundary integral in (23) will vanish for all internal points as long as the radius of the support is small enough so that the circle (sphere) lies inside the global domain U. For nodes on the boundary the portion of the circle (sphere) centered at the node that falls outside U must be eliminated as depicted in Fig. 1, and the kernel re-scaled to ensure condition (2). As a consequence, for each node that lies on the boundary of U there will be a non-zero contribution to the boundary integral in (23) stemming from the intersection of the kernel’s support with the global boundary. This region is identified as G0 in Fig. 1. The case of Gaussian kernels will be discussed in Section 6.

(22)

U0 ðrÞ

U0 ðrÞ

329

(24)

The volume integrals in Eq. (23) are carried over the support of the kernel P0, U0, and the surface integral over its boundary. Hence, to numerically solve Poisson

Eq. (27) is interpreted in terms of the fundamental lemma of the calculus of variations, which states that (27) is satisfied for all functions w(r) in U—the weight functions— if and only if the factor in brackets is identically zero. This ensures equivalence between (27) and (25). Integration by parts through application of Green’s theorem leads to

330

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

the well known weak statement: ð ð ½Vt uVw K gwdV K qwdS Z 0

(28)

Gq

U

We now apply identity (17) and theorem (15) to the last member of (29) to obtain: ð ð 0


V
vx
i P0ðrKr;3Þ uðrÞd vxi vx
i Pðr;3Þ
uðr C rÞdV Z vxi
U0 ðrÞ

U0

where the weight functions are assumed to be zero over the Dirichlet portion of the boundary. The well known Galerkin method corresponds to the special case where the weight functions w(r) are the same as the approximating—or ‘shape’—functions. Several computational methods are based on (28) including the standard Finite Element method and a number of so-called meshless methods. It should be noticed, however, that methods based on (28) cannot be fully meshless, because of need to integrate over a prescribed domain U requires a partition of it to perform the integration. In contrast, integrations in (23) are performed over arbitrary, user defined geometrical shapes—the supports of kernels P0— which are not subject to this restriction, thus allowing for a truly meshless methodology. Actually, there is another method that uses this same methodology: the local Petrov-Galerkin method (LPGM). In this method, the set of integral identities (28) are solved using special weight functions of local support that vanish over their boundary for all internal nodes. In this way, the surface integral in (28) is identically zero for all w(r) centered at internal nodes. For nodes located at the boundary of U, only the intersection of the weight function’s support with the global boundary contributes to the surface integral in (28). Hence, we see that the basic equations behind LPGM and Eq. (23) of blurred derivatives are essentially equivalent. The only difference lies on the rationale behind those equations, and the corresponding selection of approximating functions and selection of the zero order kernel/weight function. Moreover, it has recently been shown that LPGM can be considered as the general basis for several recently published variations of meshless methods. Consequently, the same applies to blurred derivatives. At this point, it is worth wondering if the formal equivalence of (23) and LPGM is merely casual or, rather, one of those methods is endowed with more generality. In the following section we show that some specific forms derived from blurred derivatives hint at the latter. But before doing so it is very illustrative to show yet another equivalence. Let us return to Eq. (23), and recall that it was derived from suitable changes of variables from the fundamental blurred approximation (19). Let us apply the same transformation once again to the first term, which originated in the laplacian. We have ð ð 0


V Z vx
i P0ðr;3Þ vx
i PðrKr;3Þ vx
i uðrÞd

vxi uðr
C rÞdV U0 ðrÞ

U0

ð Z v xi U0


vx
i P0ðr;3Þ
uðr
C rÞdV

(29)

ð


V
vxi vx
i P0ðrKr;3Þ uðrÞd


Z U0 ðrÞ

C

%vP x
i

0
i dS
uðrÞn
ðrKr;3Þ

(30)

vU0 ðrÞ

Finally, switching to derivatives on the barred variables alone as was done before we have: ð ð
V
C
V
P2ðrKr;3Þ uðrÞd P0ðrKr;3Þ gðrÞd

U0 ðrÞ

C

U0 ðrÞ

%P

%P

0
nd ^ S
K qðrÞ
ðrKr;3Þ

vU0 ðrÞ

1
S
^ rÞd nuð
ðrKr;3Þ

Z0

vU0 ðrÞ

(31) 2

1

Where P and P are the laplacian and gradient of P0, respectively. If the latter is Gaussian they are given by formulas (14) above. The salient feature of Eq. (31) is that both boundary conditions—Neumann and Dirichlet—are present. This may be useful if the field u(r) is discretized in such a way that the unknowns of the discrete problem are not values of the function at prescribed location. This is the case, for instance, of moving least squares approximations. Also, notice that Eq. (31) is frequently used as a starting point for boundary integral methods of solution. This class of methods are substantiated by just choosing a Green function of the differential problem as the zero order kernel P0 in (31). But the equation can also be implemented numerically in many other ways, including local boundary equations through proper choice of P0.

5. The functional integral formulation A very peculiar property of the blurred derivative technique is that it allows to derive functional integral formulations of evolutionary problems. As an example, we apply it to the diffusion equation: vt uðr; tÞ Z V2 uðr; tÞ C gðxÞ uðr; 0Þ Z vðrÞ u Z u0 on Gu

(32) vu hq Z q0 on Gq vn

(33)

The kernel P0 for this problem depends on both space and time variables, and integration must be carried out for the whole set of variables. But this kernel now contains two parameters: 3 for the space variables and a new parameter d

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

for time. An adequate bidimensional Gaussian kernel is: h i 2
ðt
KtÞ2 exp K ðrKrÞ K 2 2 3 d pffiffiffiffi P0ðrKr;3; (34)
t
Kt;dÞ Z p32 pd The transformations leading to (31) can be applied to the right hand side of (32), while in the left hand side we apply the time partial derivative kernel. Performing a linear approximation of the time variation of u(r,t) between two ‘time slices’ t and (tCd), integrating in time, and solving for u(r,tCd) leads to: ð
tÞdV
uðr; t C dÞ Z ðP0ðrKr;3Þ C dP2ðrKr;3Þ Þuðr;

U0 ðrÞ

ð

C

0
tÞnd ^ S
qðr;
ðrKr;3Þ

vU0 ðrÞ

U0 ðrÞ

K

%P


V
C P0ðrKr;3Þ gðrÞd


%P

1
tÞdS
^ r; nuð
ðrKr;3Þ

vU0 ðrÞ

(35) pffi Finally, the ‘space’ parameter 3 is chosen as 3Z 2 t and a special first zero order kernel P0 is used in (35) so that it is in transformed into: ð ð uðr;t CdÞ Z P0ð
rKr;3Þ uð
r; tÞdV
C P0ð
rKr;3Þ gð
rÞdV
U0 ðrÞ

C

%P

vU0 ðrÞ

U0 ðrÞ 0
^ SK r; tÞnd ð
rKr;3Þ qð


%P

1 ^ r;tÞdS
ð
rKr;3Þ nuð


vU0 ðrÞ

(36) 0

where P is simply given by (11). Details in the derivation are given elsewhere. Eq. (36) is an updating formula for the field u(r,t). However, it is important to notice that it is an approximate statement. An exact statement can be obtained by repeated application of it in such a way that the time step d tends to zero as the number of nested integrals, M, tends to infinity, keeping the product MdZDtZconstant. Such limit is known as a functional integral. However, as regards numerical implementation the approximate Eq. (36) can be used with confidence as is provides excellent results provided the time step is adequate. Also, this formulation can be applied to stationary problems by just equating the fields at both time slices: u(r,tCd)Zu(r,t). This amounts to enforcing a no-evolution condition. Details on the procedure and numerical examples are given elsewhere.

6. Discretization The main advantage of the present formulation is that it yields quite accurate results when using very simple approximations of the continuous field. The simplest one

331

is to fill the domain with N nodes—not necessarily evenly distributed—and define a set of M neighbors for each node in order to fit a polynomial around it. The set of these MC1 nodes is the local cloud. For instance, second degree approximation of the field in two dimensions is:
Z a1 C a2 x
C a3 y
C a4 x
2 C a5 x
y
C a6 y
2 ~ rÞ uð

(37)

The coefficients ai in Eq. (37) are linear combinations of the nodal values of the field, uj, on all points of the local
yÞ
in (37) can be the actual position cloud. The variables ðx; in the global frame of reference but this is highly inconvenient since it leads to large errors for points that are far from the origin. Eq. (36) indicates that for each node located at ri the kernels must be translated to that point and integrations performed around it. Hence, consistent with the local character of the approximation it is convenient to shift
yÞ
the origin to ri and measure the integration variables ðx; from it. Notice that this approximation is intrinsically meshless as it does not require a partition of the domain. Also, this approximation is valid only within a rather fuzzy region around the point but is assumed to be accurate enough inside the kernel’s support. The global continuous approximation of the solution in terms of the discrete set of nodes then becomes a post-processing problem completely separated from the numerical solution stage. This is a major difference with most methods—like FEM, BEM, EFG, etc.—where the same continuous mapping serves for both the solution and post-processing stages. A further advantage of this simple polynomial approximation is that integrations in (36) can be performed analytically even for boundary nodes. If more involved approximations are used it becomes necessary to resort to numerical quadrature. A number of schemes for integration on circles and spheres are available for this. Gaussian kernels have infinite support so that the preceding is not strictly applicable. However, they can be used for actual computations provided the mean deviation of the Gaussian is so small that P0 becomes negligible beyond the closest node. In particular, notice that all integrals involving Gaussian kernels can be computed numerically using Hermite quadrature formula [21] namely: ðN n X 2 eKx f ðxÞdx Z wi f ðxi Þ C Rn (38) KN

iZ1

where xi is the ith zero of Hermite polynomial Hn(x), and the weights and remainder are given by: pffiffiffiffi 2n n! p wi Z 2 ; Rn n ½HnK1 ðxi Þ2 pffiffiffiffi n! p 2n f ðxÞ; KN! x!N Z n (39) 2 ½2n! From the above formula for the remainder it is clear that nZ3 integrates exactly a fifth degree polynomial and thus approximation (37). In this case the most distant quadrature

332

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

point is situated a distance dz1.223 from the field point which can be considered as the radius of the ‘disc’ of integration. Hence, if parameter 3 (as defined in Eq. (3)) is sufficiently small the above reasoning for kernels of compact support is applicable. Actually this is the case as the optimum parameter 3 has been shown to be around onehalf the mean nodal distance [6].

7. Relation with kernel-based particle methods A number of meshless particle methods are based on kernel approximations of a function as defined in Eq. (1). These include Smooth Particle Hydrodynamics [12], Reproducing Kernel Particle Methods [19] and Reproducing Kernel Element Method [20]. Hence, it might then seem that the methods presented in this paper are closely related to particle methods based on kernel approximations. However, the relationship between this two sets of methods lies entirely in the point of departure: the way in which these equations are used and the final discrete equations are completely different. Briefly, in particle methods based on kernel representations Eq. (1) is used to approximate the function whose computational estimate is sought. This is fed into the governing differential equations or variational principles yielding a set of discrete equations. On the contrary, in the present paper the kernel representation (1) and its derivatives are used to transform the governing differential equations as described in this section. Approximation of the unknown function is a separate problem that will be tackled in Section 9 on numerical implementations. In particular, it will be shown that the present formulation allows to obtain very high precision when using a simple polynomial approximation of the unknown field.

Fig. 2. Unit square used to solve diffusion equation with initial temperature TZ1. Only a quarter is modeled due to symmetry.

Dirichlet conditions with TZ0. The other two have Neumann conditions qZ0. The square domain was discretized using a regular array of 21!21 nodes. A Gaussian kernel was used, and the field T(x,y,t) was approximated using simple polynomial fitting around each node. Only the first neighbors to each node where considered in the approximation. This procedure is described in detail and discussed in [6–7] and [15]. Comparison of numerical and analytical solutions along the line yZ0.2 at different times is shown in Fig. 3. A time step of 0.005 s. was used. As can be seen the agreement is excellent. Moreover, the mean nodal error after 450 time steps is EZ4!10K4. This indicates that very accurate solutions can be obtained with the simplest approximation procedure for the unknown field—ordinary polynomial fitting. 1.20

8. Numerical example 1.00

ð40Þ The numerical solution was obtained by modeling one quarter of the square due to symmetry. This is depicted in Fig. 2. The two boundaries xZ1 and yZ1 have fixed

0.80

Temperature

We illustrate the use of Eq. (36) through the solution of the diffusion equation on a unit square with initially constant temperature TZ1. The boundaries are kept at a fixed temperature TZ0. The analytical solution to this problem is: "

N 16 X ðK1Þn ð2n C 1Þ2 p2 t Tðx; y; tÞ Z 2 exp K 4 p nZ0 ð2n C 1Þ "   X N ð2n C 1Þpx ðK1Þm ! cos 2 ð2m C 1Þ mZ0   2 2

ð2m C 1Þ p t ð2m C 1Þpy cos !exp K 4 2

0.60

0.40

0.20

0.00 0.00

0.20

0.40

0.60

0.80

1.00

Fig. 3. Temperatures along the line YZ0.2 for 100, 500 and 800 time steps. Continuous lines are analytical solutions.

E. Pardo / Engineering Analysis with Boundary Elements 29 (2005) 326–333

9. Conclusions It was shown in this paper that blurred approximations of differential equations can be transformed in a number of ways by suitable change of variables and application of Leibnitz theorem for differentiation of an integral. The primary purpose for doing so was to derive self-consistent boundary condition approximations for the Functional Integral Methods described elsewhere. However, we have also shown that some well-known numerical methods also emerge as special cases, evidencing deep relationships among them. This may be an indication that the blurred derivative technique offers an appropriate framework for deriving efficient meshless methods.

References [1] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307–18. [2] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng 1996;139:3–47. [3] Duarte CA, Oden JT. H-p adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. [4] Sukumar N, Moran B, Belytschko T. The natural element method in solid mechanics. Int J Numer Methods Eng 1998;43:839–87. [5] Pardo E. Meshless method for linear elastostatics based on a path integral formulation. Int J Numer Methods Eng 2000;47:1463–80. [6] Pardo E. Convergence and accuracy of the path integral approach for elastostatics. Comput Methods Appl Mech Eng 2002;191(20): 2219–47. [7] Pardo E. Functional integral formulation of classical wave equations. J Sound Vibration 2003;261:819–37.

333

[8] Liszka T, Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980; 11:83–95. [9] Liszka T. An interpolation method for an irregular net of nodes. Int J Numer Methods Eng 1984;20:1599–612. [10] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [11] Liu WK, Jun S, Zhang YF. Reproducing Kernel particle methods. Int J Numer Methods Eng 1995;20:1081–106. [12] Monaghan JJ. An introduction to SPH. Comput Phys Commun 1988; 48:89–96. [13] Aluru NR. A point collocation method based on reproducing kernel approximations. Int J Numer Methods Eng 2000;47: 1083–121. [14] Lu YY, Belytschko T, Gu L. A new implementation of the elementfree Galerkin method. Int J Numer Methods Eng 1994. Lu YY, Belytschko T, Gu L. A new implementation of the element-free Galerkin method. Comput Methods Appl Mech Eng 1994;113: 397–414. [15] Pardo E. Blurred derivatives and meshless methods. Int J Numer Methods Eng 2003;56:295–324. [16] Atluri SN, Zhu T. A new meshless local Petrov-Galerkin approach in computational mechanics. Comput Mech 1998;22:117–27. [17] Atluri SN, Shen SP. The mesh local Petrov-Galerkin (MLPG) method: a simple & less-costly alternative to the finite element and boundary element methods, Computer Modeling in Engineering and Sciences, (2002) 3: 11–51. [18] Bathe KJ, editor. Computational fluid and solid mechanics 2003. Proceedings of the second MIT conference on computational fluid and solid mechanics, Cambridge, MA, USA. [19] Liu WK, Jun S, Li S, Adee J, Belytschko T. Reproducing kernel particle methods for structural dynamics. Int J Numer Methods Eng 1995;38:1655–79. [20] Liu WK, Han W, Lu H, Li S, Cao J. Reproducing kernel element method part I: theoretical formulation. Comput Methods Appl Mech Eng 2004;193:933–51. [21] Abramowitz M, Stegun I. Handbook of mathematical functions. New York: Dover; 1965.