Generalized finite differences for solving 3D elliptic and parabolic equations

Accepted Manuscript

Generalized Finite Differences for Solving 3D Elliptic and Parabolic Equations L. Gavete, J.J. Benito, F. Urena ˜ PII: DOI: Reference:

S0307-904X(15)00400-X 10.1016/j.apm.2015.07.003 APM 10635

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

11 March 2014 15 April 2015 3 July 2015

Please cite this article as: L. Gavete, J.J. Benito, F. Urena, Generalized Finite Differences ˜ for Solving 3D Elliptic and Parabolic Equations, Applied Mathematical Modelling (2015), doi: 10.1016/j.apm.2015.07.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • The GFDM can be applied for solving problems defined over irregular clods of points.

CR IP T

• The study of the influence of the main parameters involved in the approximation is performed. • The treatment of the Neumann boundary conditions.

AC

CE

PT

ED

M

• A new limit for stability.

AN US

• The measure of the irregularity the cloud of points

1

ACCEPTED MANUSCRIPT

Generalized Finite Differences for Solving 3D Elliptic and Parabolic Equations

CR IP T

L.Gavetea , J.J.Benitob , F.Ure˜ nac,∗ a

UPM, ETSIM, Madrid, Spain UNED, ETSII, Madrid, Spain c UCLM, IMACI, Ciudad Real, Spain b

AN US

Abstract

ED

M

The Generalized finite difference method (GFDM) is a meshfree method that can be applied for solving problems defined over irregular clouds of points. The GFDM uses the Taylor series development and the moving least squares approximation to obtain explicit formulae for the partial derivatives. In this paper, this meshfree method is used for solving elliptic and parabolic partial differential equations in 3-D. The influence of the main parameters involved in the approximation and the treatment of the Neumann boundary condition are shown. Parabolic equations have been solved using an explicit method and the criterion for stability has been improved taking into account the irregularity of the cloud of points. The numerical results show the high accuracy obtained.

PT

Keywords: meshfree method, generalized finite difference method, moving least squares, elliptic, parabolic.

CE

1. Introduction

AC

GFDM is included in the so-called Meshless Methods or Meshfree Methods(MM) [1]. The GFDM has been developed for many years. Jensen [2] introduces fully arbitrary mesh, Perrone and Kao [3] suggested the use of additional nodes in the star. Liszka [4,5] and Orkisz [5,6] introduce the weight ∗

Corresponding author Email addresses: [email protected] (L.Gavete), [email protected] (J.J.Benito), [email protected] (F.Ure˜ na)

Preprint submitted to Applied Mathematical Modeling

July 13, 2015

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

function and a star with eight nodes, using moving least squares approximation. Benito et al.[7] show explicit formulae for GFDM in the bi-dimensional case. An h-adaptive method in GFDM for 2-D and 3-D are described in Benito et al.[8] and Ure˜ na et al.[9] respectively. Gavete et al.[10] show improvements of GFDM and comparison with other meshless method. Ure˜ na et al.[11] show the solution of the advection-diffusion equation using GFDM. Gavete et al. [12] describes how to solve dynamic analysis of different problems and Ure˜ na et al. [13] show the application to seismic wave propagation in 2D. Gavete et al. [14] describes the application to electrical conductivity of a tissue. GFDM was also discussed and applied for different practical tasks as in Chang, Fan and Kuo [15] for solving two-dimensional non linear obstacle problems, or Mochnacki, Majchrzak [16] for numerical modelling of casting solidification. Paper [9] presented an h-adaptive algorithm in GFDM and a posterior error indicator to estimate the error solving second order partial differential equations. In [17] the extension of the GFDM to the explicit solution of parabolic and hyperbolic equations was introduced. This paper describes how the GFDM can be applied for solving elliptic and parabolic partial differential equations in 3-D. The differences with previous papers [9,17], are related with the study of the influence of key parameters involved (weighting function, shape of the domain, criterium of nodes selection), the definition and application of the index of irregularity of a cloud of nodes, the treatment of Neumann boundary conditions and the improvement of the stability limit. The influence of the selection of the nodes of the star and the weight functions, etc, is analyzed in three different irregular domains cube, cylinder and sphere. For each criterium of nodes selection, the accuracy of the results depends on the weighting function and the index of irregularity of the clouds of nodes considered in each one of the domains. Also the treatment of the Neumann boundary condition has been included. Parabolic equations have been solved using an explicit finite difference formulae and the new limit of stability, improved as compared with [17]. The numerical results show the accuracy obtained over irregular clouds of points. The outline of the paper follows. In Section 2 the explicit formulae in the GFD method are given. In Section 3 the procedure to approximate the truncation error is described and in Section 4 an explicit scheme is considered including the new stability limit. The key parameters of the method are described in Section 5, and in section 6 the performance of the method is analyzed for different clouds of points. We close with some conclusions in 3

ACCEPTED MANUSCRIPT

Section 7. 2. Explicit formulae in Generalized Finite Differences

∂ 2U ∂ 2U ∂U ∂ 2U ∂ 2U ∂ 2U ∂ 2U + a + a + a + + a + a + a 5 6 7 2 3 4 ∂x2 ∂y 2 ∂z 2 ∂x∂y ∂x∂z ∂y∂z ∂x ∂U ∂U a8 + a9 = f1 (x, y, z) in (x, y, z) ∈ Ω ⊂ R3 (1) ∂y ∂z

with the boundary condition: b1

AN US

a1

CR IP T

Let us assume a problem governed by the following second-order partial differential equation in 3-D of the unknown function U(x,y,z)

∂U ∂U ∂U + b2 + b3 + b4 U = f2 (x, y, z) in (x, y, z) ∈ Γ ⊂ R3 ∂x ∂y ∂z

(2)

M

where Γ is the boundary of Ω, f1 and f2 are two known functions, ai (i = 1, , 9) and bj (j = 1, , 4) are constants. To obtain the explicit difference formulae of the partial derivatives [9,17] the norm B(u) is defined as

li2 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 + h k + h l + k l )w(hi , ki , li )]2 (3) i i i i i i 2 2 ∂z ∂x∂y ∂x∂z ∂y∂z

PT

+

ED

N X ∂u0 ∂u0 h2i ∂ 2 u0 ki2 ∂ 2 u0 ∂u0 + ki + li + + B(u) = [u0 − ui + hi ∂x ∂y ∂z 2 ∂x2 2 ∂y 2 i=1

AC

CE

where U0 is the value of the function at the node of the coordinates (x0 , y0 , z0 ) and Ui is the function values at the nodes i, with i = 1, ..., N the coordinates (xi , yi , zi ), hi = xi − x0 , ki = yi − y0 , li = zi − z0 and w(hi , ki , li ) is the weighting function. Then, by minimizing the norm B(u) the following linear equation system is obtained A · Du = b

4

(4)

ACCEPTED MANUSCRIPT

The matrix A (of 9x9) and the vectors Du and b are given by: hi k i w

2

i=1 N X

hi l i w

N X h3 i

2

i=1 2

i=1 2 2 ki w

i=1

N X

ki li w

2

i=1

2

N X 2 2 li w

N X hi ki2

2

i=1

N X h2 i ki

i=1

i=1

w

w

2

N X h2 i li 2 w 2

i=1

N X h4 i

i=1 4

w

2

2

w

N X hi l 2 i

2

w

N X ki l2 i

i=1

N X ki2 li 2 w 2

N 2 X h2 i ki

i=1

4

i=1 4

2

i=1 2

w

w

2

2

w

2

N X ki2 l2 i

w

i=1

Y

M

hi k i w

N X

2 2 hi ki w

2

i=1

5

PT

                                          

N X

hi ki li w

N X

4

4

hi ki li w

2

2

2

N X h3 i ki

w

N X hi ki3

w

2

2

i=1

2

2

hi li w

2

i=1

2 2

2

N X h3 i li

2

i=1

N X hi ki l2 i N X

N X

2

i=1

i=1

w

2

i=1

2

N X hi ki2 li

2

i=1

w

hi ki w

2

N X hi l3 i

2

i=1

2

i=1

w

N X

2

w

w

2

2

i=1 N X

2 2

hi li w

i=1

2

2

hi ki li w



N X

2

i=1

∂u0 ∂u0 ∂u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 ∂x ∂y ∂z ∂x2 ∂y 2 ∂z 2 ∂x∂y ∂x∂z ∂y∂z                                            

2

hi l i w

i=1

M

 N X    (−u0 + ui )hi w2     i=1    N  X    (−u0 + ui )ki w2     i=1    N  X    (−u0 + ui )li w2     i=1   N  X  h2    (−u0 + ui ) i w2   2   i=1   N  X k2 b= (−u0 + ui ) i w2 2   i=1   N  2 X  l   (−u0 + ui ) i w2    2  i=1    N  X    (−u0 + ui )hi ki w2     i=1    X N     (−u0 + ui )hi li w2     i=1   N  X    (−u0 + ui )ki li w2  

N X

2

i=1

w

i=1

N X

i=1 2

N 2 X h2 i li

i=1 4



2

i=1

N X l4 i

S

w

N X l3 2 i w 2

i=1

N X ki4

CE AC

2

N X ki3

i=1

ED

Du =

N X

CR IP T

N X

AN US

A

N X 2 2  hi w  i=1                        =                          

M



2

2 hi ki li w    i=1   N X 2 2   ki li w   i=1   N  X 2 2  ki li w    i=1   N 2 X hi ki li 2  w   2  i=1   N 3 X ki li 2   w   2  i=1   N 3 X ki li 2   w   2  i=1   N  X 2 2 hi ki li w    i=1   N X 2 2 hi ki li w    i=1   N X 2 2 2  k l w  i i

i=1

T

(5)

(6)

(7)

ACCEPTED MANUSCRIPT

ED

M

AN US

CR IP T

On solving the system of Eqs.4, the following explicit difference formulae are obtained  N X    ∂U0 = −m1 u0 +  m1i ui  ∂x   i=1    N  X  ∂U0   = −m u + m2i ui 2 0  ∂y    i=1   N  X   ∂U0  = −m u + m3i ui 3 0    ∂z  i=1   N  X  2  ∂ U0  = −m u + m4i ui  4 0 ∂x2    i=1   N  X 2 ∂ U0 (8) = −m5 u0 + m5i ui + Θ(h2i , ki2 , li2 ) ∂y 2   i=1   N  X   ∂ 2 U0   = −m6 u0 + m6i ui  ∂z 2   i=1    N  X  2U ∂  0  m7i ui  ∂x∂y = −m7 u0 +    i=1   N  X   ∂ 2 U0  = −m8 u0 + m8i ui  ∂x∂z    i=1   N  X   ∂ 2 U0  m9i ui   ∂y∂z = −m9 u0 + i=1

aj (−mj u0 +

CE

9 X

PT

Inserting Eqs.8 into Eq.1 leads to

j=1

N X i=1

mji ui ) = f0 ⇔ −u0

AC

hence

−λ0 u0 +

N X

λi ui = f0 ,

i=1

9 X

aj mj +

j=1

with

N X 9 X ( aj mji )ui = f0 (9) i=1 j=1

N X

λi = λ0

(10)

i=1

On using the same procedure for each one of the stars leads to obtain of a system of equations. On solving it and using Eqs.4 in the interior nodes of the domain, the partial derivatives can be calculated. 6

ACCEPTED MANUSCRIPT

3. Truncation error in GFD

CR IP T

To initiate the calculation of the truncation error for space derivatives, Taylors expansion including higher order derivatives is used in function B ∗ (u) N X ∂u0 ∂u0 ∂u0 1 ∂u0 ∂u0 ∂u0 2) +ki +li + (hi +ki +li ) + B (u) = [(u0 −ui +hi ∂x ∂y ∂z 2 ∂x ∂y ∂z i=1 ∗

1 ∂u0 ∂u0 ∂u0 3) 1 ∂u0 ∂u0 ∂u0 4) (hi + ki + li ) + (hi + ki + li ) +)w]2 (11) 3! ∂x ∂y ∂z 4! ∂x ∂y ∂z

A · Du =

(

N X

Ξhi

i=1

N X

Ξki

i=1

N X

Ξli

i=1

N X

Ξ

i=1

AN US

If the function Eq.11 is minimized with respect partial derivatives to the second order, the following linear equation is defined N N N N N X X h2i X ki2 X li2 X Ξ Ξ Ξhi ki Ξhi li Ξki li 2 i=1 2 i=1 2 i=1 i=1 i=1

where

)T

(12)

where

(

N X

Υhi

i=1

N X

Υki

i=1

N X i=1

Υli

N X i=1

Υ

N N N N N X X h2i X ki2 X li2 X Υ Υ Υhi ki Υhi li Υki li 2 i=1 2 i=1 2 i=1 i=1 i=1

)T (14)

PT

T E = A−1

ED

M

∂u0 ∂u0 ∂u0 3) 1 ∂u0 ∂u0 ∂u0 4) 1 Ξ = (−u0 +ui − (hi +ki +li ) − (hi +ki +li ) +· · · )w2 3! ∂x ∂y ∂z 4! ∂x ∂y ∂z (13) the matrix A and Du are defined by Eqs.5 and 6 respectively. Then, the truncation error is

AC

CE

1 ∂u0 ∂u0 ∂u0 3) 1 ∂u0 ∂u0 ∂u0 4) Υ = (− (hi + ki + li ) − (hi + ki + li ) + · · · )w2 3! ∂x ∂y ∂z 4! ∂x ∂y ∂z (15) By considering bounded derivatives in Eq. 14 lim

(hi ,ki ,li )→(0,0,0)

TE → 0

(16)

Then, the truncation error condition given in Eq. 14 shows the consistency of the approximation for the Eq. 1 with constant coefficients.

7

ACCEPTED MANUSCRIPT

4. Parabolic Equations.Convergence The equation considered, is:

CR IP T

∂ 2U ∂ 2U ∂U ∂ 2U ∂ 2U ∂ 2U = a1 2 + a2 2 + a3 2 + a4 + a5 ∂t ∂x ∂y ∂z ∂x∂y ∂x∂z 2 ∂ U ∂U ∂U ∂U + a6 + a7 + a8 + a9 t > 0, (x, y, z) ∈ Ω ⊂ R3 (17) ∂y∂z ∂x ∂y ∂z with the boundary condition:

∂U ∂U ∂U + b2 + b3 + b4 U = g1 (t) in (x, y, z) ∈ Γ ⊂ R3 ∂x ∂y ∂z

AN US

b1

and initial condition

U (x, y, z, 0) = g2 (x, y, z)

(18)

(19)

where g1 and g2 are two known functions.

un+1 − un0 ∂U (x0 , y0 , z0 , n4t) = 0 + Θ(4t) ∂t 4t

(20)

M

Substituting Eq.(10) and Eq.(20) into Eq.(17) leads to N

PT

ED

X un+1 − un0 0 = −λ0 un0 + λi uni ⇔ un+1 = 0 4t i=1 (1 − λ0 4t)un0 + 4t

N X

λi uni + Θ(4t, h2i , ki2 , li2 ) (21)

i=1

AC

CE

The Lax theorem [18] states that for a consistent finite difference method for a well-posed linear initial problem, the method is convergent if and only if it is stable. In this paper an improvement of the limit of stability of [17] is obtained as follows, using Von Neumann stability analysis. ξ = (1 − λ0 4t) + 4t 1 − 4t

N X

λj eiν

T 4r

j

=

j=1

N X j=1

λj (1 − cos(ν T 4rj )) + i4t 8

N X j=1

λj sin(ν T 4rj ) (22)

ACCEPTED MANUSCRIPT

Then kξk = 1 − 24t

N X j=1

T

λj (1 − cos(ν 4rj )) + (4t (4t

N X j=1

with Eq.10 and kξk2 ≤ 1 N X j=1

j=1

λj sin(ν T 4rj ))2 +

λj (1 − cos(ν T 4rj )))2 (23)

λj (1 − cos(ν T 4rj )) ⇒ 4t ≤

AN US

(2λ0 4t)2 + (2λ0 4t)2 ≤ 24t

N X

CR IP T

2

4 5|λ0 |

(24)

Taking into account that λ0 will be different in each star of irregular domain, in order to ensure stability it is necessary to select 4t as 4 5|λ0 |max

(25)

M

4t ≤

This formula (Eq.25) improves previously calculated stability limit in [17].

ED

5. Essential Factors for obtaining the Explicit Formulae in GFD From the above it may be seen that the coefficients mi and mij of the explicit formulae in GFD (Eq.(8)) depends on the following factors:

CE

PT

1. N 2. (hi , ki , li ) 3. w(hi , ki , li )

AC

For problems in 3-D, the minimum number of nodes of the star (N ) is nine, however much more nodes are normally necessary to assure the quality of the star equations. The main drawback of the GFDM is the possibility of obtaining ill-conditioned stars of nodes. However, the quality of numeral results can be evaluated by obtaining the residual at each point, which must be very small (near zero) and also with a uniform distribution over the entire domain. It is also possible to increase N to obtain correct residual values over the entire domain. The relative coordinates (hi , ki , li ) are defined by the selection of the nodes 9

ACCEPTED MANUSCRIPT

AN US

CR IP T

in the star, and have a great influence on the results. Jensen [2] only considered the distance criterion. This criterion may produce distorted stars with unevenly distributed nodes around the central node which is, subsequently, reflected by more imprecise results [9]. A viable alternative would be to correct this last method in terms of distance, assuring that none of the nodes are set at greater distances to those indicated, and thereby correcting the effect of node density irregularity. A second method, denominated the eight octant criterion, consists of the selection of the two or three nearest nodes per octant, the said octant being established in the same manner as that indicated above by Cartesian axes around the central node of the star. The said method corrects the problems of the distance criterion. The equation Eq.(10) shows the influence of the weighting function giving more influence to the nodes closest to the central node. In this paper we have used two different weighting functions: Potential and Exponential. 6. Numerical Results

1 ; (dist)3

ED

wA =

M

In this section some elliptic and parabolic partial differential equations are solved in different domains and the results are compared with analytical solutions. The weighting functions used have been wB = e−2.5(dist)

2

(26)

CE

PT

where (dist) is the euclidean distance from each node of the star to the central node. The % global error has been evaluated using the formula v u NT uX u (sol(i) − exac(i))2 u t i=1

AC

% Global error =

NT |exacmax |

× 100

(27)

where sol(i) is the GFDM solution at the node i, exac(i) is the exact value of the solution at the node i, exacmax is the maximum value of the exact values in the cloud of nodes considered and N T is the total number of nodes of the domain considered. Three different irregular clouds of nodes (cube, cylinder and sphere) have been considered corresponding to the following domains(see Fig. 1, 2 and 3). 10

ACCEPTED MANUSCRIPT

CR IP T

1. Cube: Ω1 = {(x, y, z) ∈ R3 /0 ≤ x, y, z ≤ 1} To generate the nodes a initial regular cloud is used and then the interior nodes are moved just a little from the initial location using aleatory numbers. 2. Cylinder: Ω2 = {(x, y, z) ∈ R3 /0 ≤ x2 + y 2 ≤ 22 , −2 ≤ z ≤ 2} In this case the nodes are generated using cylindrical coordinates with

((r = 0.5, 2, step = 0.5), (θ = 0◦ , 337.5◦ , step = 22.5◦ ), (z = −2, 2, step = 0.5))

3. Sphere: Ω3 = {(x, y, z) ∈ R3 /0 ≤ x2 + y 2 + z 2 ≤ 52 } Spherical coordinates are used to generate the nodes with

AN US

((ρ = 1, 5, step = 1), (ϕ = 0◦ , 157.5◦ , step = 22.5◦ ), (θ = 0◦ , 337.5◦ , step = 22.5◦ )) The number of nodes corresponding to the domains are Ω1 = 729, Ω2 = 576 and Ω3 = 720 nodes. 6.1. Elliptic equation Let us consider firstly the elliptic equation:

M

∂ 2U ∂ 2U ∂ 2U + + = 0 in (x, y, z) ∈ Ωi (i = 1, 2, 3) ∂x2 ∂y 2 ∂z 2

(28)

ED

with Dirichlet boundary conditions in Γ. The exact solution is:

CE

PT

U (x, y, z) = 4x2 − 2y 2 − 2z 2

AC

Clouds of nodes Cube Cylinder Sphere

(29)

Table 1: % Global Error

wA and distance 0.29×10−13 0.33×10−11 0.18×10−2

wA and octant 0.15×10−13 0.91×10−11 0.36×10−7

wB and distance 0.17×10−13 0.20×10−13 0.57×10−3

wB and octant 0.16×10−13 0.14×10−13 0.12×10−8

24 nodes have been selected for each one of the stars using both criteria distance and octant (3 nodes by octant). The results obtained for the % global error for two different weighting functions Eq.29 are given in Table 1. 11

Fig.2 Cube

AC

CE

PT

ED

M

Fig.1: Cylinder

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig.3: Sphere

12

ACCEPTED MANUSCRIPT

CR IP T

These results as shown in Table 1 are very accurate. The accuracy of the results depends of the weighting function and the index of irregularity of the cloud of nodes considered in each one of the domains. We are going to define the index of irregularity of a cloud of nodes (IIC). Let us consider in each node (i) of a star (o) the distance to the central node of the star q (30) di = h2i + ki2 + li2

then we calculate the medium value of the distances to the central node of the star as nns X di i=1

(31) nns being nns the number of nodes of the star (o), excluding the central node. We can calculate the medium radius of the cloud of nodes as

AN US

rm0 =

rm0

i=1

M

rmc =

nni X

nni

(32)

PT

ED

being nni the number of interior nodes of the domain, then the index of irregularity of a cloud of nodes (IIC) can be calculated as v u nni uX u (rm0 − rmc)2 u t i=1 IIC = (33) nni

AC

CE

In the following Table 2 the different values of (IIC) obtained for the different cases considered are given As it is shown in Table 2 the (IIC) obtained is greater for the spherical domain, and accordingly the results of the errors obtained previously in Table 1 are greater. To solve problems with Neumann boundary conditions it is important to select an accurate procedure. In this paper the selected procedure is as follows: For each one of Neumann node ξi , two new nodes are added one ηi , is located in the interior of the domain and the other one µi is located out of the boundary. Both nodes must be in the same line in the direction normal to the boundary and with the same distance h to the corresponding Neumann 13

ACCEPTED MANUSCRIPT

IIC distance octant Cube 0.034 0.07 Cylinder 0.084 0.106 Sphere 0.201 0.214

CR IP T

Table 2: Irregularity Index (IIC)

u(ηi ) − u(µi ) = 2h

AN US

node. Then the Neumann nodes are considered as interior nodes and the new equations needed to solve the system of equations are obtained using the classical central finite difference approximation of first derivative ∂U (ξi ) , ∂n

i = 1, · · · , number of N eumann nodes (34)

M

being n the normal external to the domain. Let us consider now the case of Neumann condition for Eq.31, being the exact solution U (x, y, z) = ex sin y + ey sin z + ez sin x (35)

AC

CE

PT

ED

Let us consider an irregular cloud of nodes, in the domain Ω = {x ∈ R3 |0 < x < 1, 0 < y < 1, 0 < z < 1} with the boundary Γ = {x ∈ R3 |x = 0, x = 1, y = 0, y = 1, z = 0, z = 1} with Neumann boundary conditions in Γ0 =}x ∈ R3 |0 < x < 1, 0 < y < 1, z = 0} and Dirichlet boundary conditions in Γ − Γ0 . In the following Table 3 the % error obtained using 18 nodes stars with the distance criterium are given considering the Neumann-Dirichlet boundary condition and also for the case of considering only Dirchlet boundary condition. Potential weighting function has been used. Table 3: % Error

% Error

Dirichlet 4.19 × 10−4

Neumann-Dirichlet 2.325 × 10−3

Note that in Eqs. 1 and 2, the coefficients a and b are constants, however the 14

ACCEPTED MANUSCRIPT

CR IP T

method can be used also in the case of problems with variable coefficients (a = a(x, y, z) and b = b(x, y, z)), see [14] where the GFDM is used to simulate the electrical conductivity of different tissues including the case when the direction of the fibres varies throughout the tissue that corresponds to the case of a poisson equation with variable coefficients. 6.2. Parabolic equation Let us consider secondly equation: ∂ 2U ∂U ∂ 2U ∂ 2U = + + ∂t ∂x2 ∂y 2 ∂z 2

t > 0 in (x, y, z) ∈ Ωi (i = 1, 2, 3)

(36)

AN US

with Dirichlet boundary conditions and the initial condition U (x, y, z, 0) = sin π the exact solution is:

x+y+z 3

in Γ

(37)

x+y+z (38) 3 The global error obtained for the three domains previously defined is given in Table 4. Only the case of three nodes by octant with the weighting function 1 wA = (dist) 3 has been considered. π2 t 3

sin π

ED

M

U (x, y, z, t) = e−

Table 4: % Global Error (after 20 time steps)

CE

PT

Clouds of nodes % global error Cube 0.14×10−3 Cylinder 0.20×10−3 Sphere 0.35×10−3

AC

The influence on the global error by using different values of the time increment is given in Figures 4, 5 and 6. As it is shown in Figures 4, 5 and 6 the global error decreases when the time step decreases. The maximum local error is evaluated for each time step using the formula M aximum local error = maxi |sol(i) − exac(i)|

(39)

As it is shown in figures 7, 8 and 9 the Maximum local error increases with the number of time steps until it remain stable. 15

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Fig.4: % global error (irregular cube) in the last time step calculated versus 4t

Fig.5: % global error (irregular cylinder) in the last time step calculated versus 4t

16

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Fig.6: % global error (irregular sphere) in the last time step calculated versus 4t

Fig.7: Maximum local error (irregular cube) in the last time step calculated versus n, for 4t = 0.0001

17

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Fig.8: Maximum local error (irregular cylinder) in the last time step calculated versus n, for 4t = 0.0001

Fig.9: Maximum local error (irregular sphere) in the last time step calculated versus n, for 4t = 0.0001

18

ACCEPTED MANUSCRIPT

7. Conclusions

AN US

CR IP T

In this paper we have shown that the meshfree GFDM is an interesting way of solving elliptic and parabolic partial differential equations in 3-D problems defined over irregular clouds of points. The truncation errors have been defined. The treatment of the Neumann boundary condition has been developed. Different examples of irregular clouds of nodes have been solved. The index of irregularity of a cloud of nodes (IIC) has been defined and it is related with the global error obtained. The stability limit for solving parabolic problems using an explicit method has been improved. The numerical results obtained show the high accuracy obtained. In the case of parabolic equations the global error decreases when the time step decreases, and the maximum local error tends to remain bounded and stable when the number of time steps is increased. 8. Ackonwledgement

M

The authors acknowledge the support of the Escuela T´ ecnica Superior de Ingenieros Industriales (UNED) of Spain, project 2014-IFC02. References

ED

[1] G. R. Liu, Meshfree Methods, CRC press (2003).

PT

[2] P. S. Jensen, Finite difference technique for variable grids, Computer& Structures 2 (1972) 17–29.

CE

[3] N. Perrone, R. Kao, A general finite difference method for arbitrary meshes, Computer& Structures 5 (1975) 45–58.

AC

[4] T. Liszka, An interpolation method for an irregular net of nodes, International Journal for Numerical Methods in Engineering 20 (1984) 1599–1612.

[5] T. Liszka, J. Orkisz, The Finite Difference Method at Arbitrary Irregular Grids and its Application in Applied Mechanics, Computer & Structures 11 (1980) 83–95. [6] J. Orkisz, Finite Difference Method (Part, III) in handbook of Computational Solid Mechanics, M. Kleiber (Ed.), Spriger-Verlag, Berlin 1998. 19

ACCEPTED MANUSCRIPT

[7] J. J. Benito, F. Ure˜ na, L. Gavete, Influence several factors in the generalized finite difference method, Applied Mathematical Modeling 25 (2001) 1039–1053.

CR IP T

[8] J. J. Benito, F. Ure˜ na, L. Gavete, R. Alvarez, An h-adaptive method in the generalized finite difference, Comput. Methods Appl. Mech. Eng. 192 (2003) 735–759.

AN US

[9] F. Ure˜ na, J. J. Benito, R. Alvarez, L. Gavete, Computational error approximation and h-adaptive algorithm for the 3-D generalized finite difference method, International Journal for Computational Methods in Engineering Science and Mechanics 6 (1) (2005) 31–39. [10] L. Gavete, M. L. Gavete, J. J. Benito, Improvements of generalized finite difference method and comparison with other meshless method, Applied Mathematical Modeling 27 (2003) 831–847.

M

[11] F. Ure˜ na, J. J. Benito, L. Gavete, Application of the generalized finite difference method to solve the advection-diffusion equation, Journal of Computational and Applied Mathematics 235 (2011) 1849–1855

ED

[12] L. Gavete, F. Ure˜ na, J. J. Benito, E. Salete, M. L. Gavete , A note on the dynamic analysis using the generalized finite difference method, Journal of Computational and Applied Mathematics 252 (2013) 132–147.

CE

PT

[13] F. Ure˜ na, J. J. Benito, E. Salete, L. Gavete, A note on the application of the generalized finite difference method to seismic wave propagation in 2D, Journal of Computational and Applied Mathematics 236 (2012) 3016–3025.

AC

[14] M.L, Gavete, F. Vicente, L. Gavete, F.Ure˜ na, J. J. Benito, Solving anisotropic elliptic and parabolic equations by a meshless method; simulation of the electrical conductivity of a tissue, International Journal of Computer Mathematics 89 (13-14) (2012) 1914–1926.

[15] H. F. Chan, C. M. Fan, C. W. Kuo, Generalized finite difference method for solving two-dimensional non-linear obstacle problems, Engineering Analysis with Boundary Elements 37 (2013) 1189–1196.

20

ACCEPTED MANUSCRIPT

[16] B. Mochnacki, E. Majchrzak, Numerical modelling of casting solidification using generalized finite difference method, Materials Science Forum 638-642 (2010) 2676–2681.

CR IP T

[17] J. J. Benito, F. Ure˜ na, L. Gavete, Solving parabolic and hyperbolic equations by the generalized finite difference method, Journal of Computational and Applied Mathematics 209 (2007) 208–233.

AC

CE

PT

ED

M

AN US

[18] A. R. Mitchell, D. F. Griffiths,The Finite Difference Method in Partial Differential Equations, John Wiley & Sons (1980)

21