Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM)

Accepted Manuscript Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM) F. Ureña, L. Gavete, A. García, J.J. Benito, A.M. Vargas

PII: DOI: Reference:

S0377-0427(18)30087-6 https://doi.org/10.1016/j.cam.2018.02.016 CAM 11522

To appear in:

Journal of Computational and Applied Mathematics

Received date : 12 September 2017 Revised date : 9 January 2018 Please cite this article as: F. Ureña, L. Gavete, A. García, J.J. Benito, A.M. Vargas, Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM), Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.02.016 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.

Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM) F.Ure˜ naa,∗, L.Gaveteb , A.Garc´ıaa , J.J.Benitoa , A.M.Vargasc a b

UNED, ETSII, Madrid, Spain UPM, ETSIM, Madrid, Spain c UCM, Madrid, Spain

Abstract The generalized finite difference method (GFDM) has been proved to be a good meshless method to solve several linear partial differential equations (PDEs): wave propagation, advection-diffusion, plates, beams, etc. The GFDM allows us to use irregular clouds of nodes that can be of interest for modelling non-linear parabolic PDEs. This paper illustres that the GFD explicit formulae developed to obtain the different derivatives of the pde’s are based in the existence of a positive definite matrix that it is obtained using moving least squares approximation and Taylor series development. Criteria for convergence of fully explicit method using GFDM for different non linear parabolic pdes are given. This paper shows the application of the GFDM to solving different non-linear problems including applications to heat transfer, acoustics and problems of mass transfer. Keywords: meshless methods, generalized finite difference method, non-linear parabolic partial differential equations. 1. Introduction Modern numerical methods, in particular those for solving non-linear PDEs, have been developed in recent years using finite differences, finite ∗

Corresponding author Email addresses: [email protected] (F.Ure˜ na), [email protected] (L.Gavete), [email protected] (A.Garc´ıa), [email protected] (J.J.Benito), [email protected] (A.M.Vargas) Preprint submitted to Journal Computational and Applied MathematicsFebruary 14, 2018

elements, finite volume or spectral methods. A review of numerical methods for non-linear partial differential equations is given by Polyanin [1] and Tadmor [2]. In this paper we use a meshless method called generalizad finite difference method (GFDM) for solving different parabolic non-linear pdes. Researchers as Jensen [3], Liszka and Orkisz [4], Orkisz [5] and Perrone and Kao [6] have contributed to develop the GFDM in different aspects of its applications. Benito, Gavete and Ure˜ na [7, 8, 9, 10, 11] have developed the explicit formulae and h-adaptive method for the solution of the pdes in 2-D. GFDM was also applied for different practical tasks as in Chan, Fan and Kuo [12] for solving two dimensional non-linear obstacle problems, Zhang, Ren, Fan and Li for simulating two dimensional sloshing phenomenon [13] or Mochnacki and Majchrak [14] for numerical modeling of casting solidification. This paper shows that the GFDM can be applied for solving parabolic nonlinear pde’s with different irregular clouds of points in 2D. In this paper, this meshless method is used for solving non-linear parabolic partial differential equations in 2-D. Parabolic equations have been solved using an explicit method and the convergence has been studied taking into account the irregularity of the cloud of points. The numerical results show the high accuracy obtained. The paper is organized as follows. In section 2 the explicit method and GFDM to solve non-linear parabolic partial differential are showed. In section 3, convergence is studied. Sections 4 exposes the results obtained when solving different non-linear problems. Finally, in section 5, some conclusions are obtained. 2. Explicit method and GFDM: application to non-linear parabolic partial differential equations Consider the following non-linear problem in the domain D = [0, T ] × Ω with Ω ⊂ R2 ∂U = LΩ [U ] (1) ∂t where LΩ [U ] is a non-linear operator, Γ is the boundary of the domain Ω, with boundary condition: UΓ = f (t) (2)

2

and initial condition U (x, y, 0) = g(x, y)

(3)

where f and g are two known functions. To solve the problem described by Eqs.(1), (2) and (3), using the explicit method, time derivative is approximated by un+1 − un0 ∂u(x0 , y0 , n4t) 0 = + Θ(4t) ∂t 4t

(4)

and spatial derivatives are approximated using GFD [7, 8, 9, 10], as described in the following section. 2.1. Explicit formulae Let Ω ⊂ R2 be a domain and M = {x1 , . . . , xN } ⊂ Ω a discretization of the domain Ω with N points. Each point of the discretization M is denoted as a node. For each one of the nodes of the domain, where the value of U is unknown, Es -star is defined as a set of selected points Es = {x0 ; x1 , . . . , xs } ⊂ M with the central node x0 ∈ M and xi (i = 1, . . . , s) ∈ M is a set of points located in the neighbourhood of x0 . In order to select the points different criteria as four quadrants or distance can be used. Consider an Es -star with the central node x0 , where (x0 , y0 ) are the coordinates of the central node, (xi , yi ) are the coordinates of the ith node in the Es -star, and hi = xi − x0 and ki = yi − y0 . If U0 = U (x0 ) is the value of the function at the central node of the star and Ui = U (xi ) are the function values at the rest of the nodes, with i = 1, ..., s, then, according to the Taylor series expansion it is known that:   ∂U0 1 2 ∂ 2 U0 2 ∂ 2 U0 ∂ 2 U0 ∂U0 +ki + h +ki +2hi ki +..., i = 1, ..., s Ui = U0 +hi ∂x ∂y 2 i ∂x2 ∂y 2 ∂x∂y (5) We define: h2 k 2 ci T = {hi , ki , i , i , hi ki } (6) 2 2 ∂u0 ∂u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 DT = { , , , , } (7) ∂x ∂y ∂x2 ∂y 2 ∂x∂y 3

If in Eq.(5) the higher than second order terms are ignored, a second order approximation for the Ui function is obtained. This is indicated as ui . It is then possible to define the function B(u) as: B(u) =

s X ∂u0 ∂u0 + ki + [(u0 − ui ) + hi ∂x ∂y i=1

2 ∂ 2 u0 1 ∂ 2 u0 2 2 2 ∂ u0 )] wi (8) + (h2i + k + 2h k i i i 2 ∂x2 ∂y 2 ∂x∂y

where wi = w(hi , ki ) are positive symmetrical weighting functions, having the property as defined in Lankaster and Salkauskas [19] that are functions decreasing monotonically in magnitude as the distance to the center of the weighting function increases. See also Levin [20]. 1 Some weighting functions as potential or exponential exp(−n(dist2 )) n dist can be used, where n ∈ N. If the norm given by Eq.(8) is minimized with respect to the partial derivatives, the following system of linear equations is obtained (for each node of the discretization where u is unknown) A(hi , ki , wi )D = b(hi , ki , wi , u0 , ui ) where    A= 

h1 k1 .. .

h2 k2 .. .

··· ··· .. .

h1 k1 h2 k2 · · ·

hs ks .. . hs ks

    

ω12

ω22 ···

ωs2

    

(9)

h1 k1 · · · h2 k2 · · · .. .. .. . . . hs ks · · ·

h1 k1 h2 k2 .. .

hs ks (10) The matrix A is positive definite matrix and the approximation is of second order Θ(h2i , ki2 ) [15]. Assumption stated on the nodes distribution is related with the criterium selected to choose the stars nodes, because in each one of the stars, the Matrix A must be positive definite to obtain the approximated derivatives, hence the existence of ill-conditioned stars of nodes is restricted.

4

    

Thus, spatial derivatives using GFD [7, 8, 9, 10] are denoted by  s s X X ∂u(x0 , y0 , n4t)  2 2   = −λ0 u0 + λi ui + Θ(hi , ki ), with λ0 = λi   ∂x   i=1 i=1  s s  X X  ∂u(x0 , y0 , n4t)  2 2  = −µ u + µ u + Θ(h , k ), with µ = µi  0 0 i i 0 i i  ∂y i=1

s

i=1

X  ∂ 2 u(x0 , y0 , n4t) ∂ 2 u(x0 , y0 , n4t)   + = −m u + mi ui + Θ(h2i , ki2 ) 0 0  2 2  ∂x ∂y   i=1  s  X     mi  with m0 = i=1

(11)

3. Convergence of the scheme for non-linear parabolic pde’s In this section convergence of non-linear parabolic pde’s, using GFDM, is studied. We will do so by introducing the following definitions: • A partial differential equation is semilinear if the coefficients of its highest derivatives are functions of the space variables only. • A partial differential equation is quasi-linear if it is linear in its highest derivatives. In this paper, the following cases are considered for LΩ [U ] of Eq.(1), for their importance in applied science and technology: Semilinear ∂ 2U ∂ 2U + + F (U ) (12) LΩ [U ] = ∂x2 ∂y 2 Quasilinear (Case a) ∂ 2U ∂ 2U + ) ∂x2 ∂y 2

(13)

∂ 2U ∂ 2U + ) + F (U ) ∂x2 ∂y 2

(14)

LΩ [U ] = K(U )( Quasilinear (Case b) LΩ [U ] = K(U )(

The method which has been used consists in finding an error bound directly from the discretized expressions and, therefore, the concept of stability is not 5

used. However this concept appears when we consider that the error must be bounded. In the following paragraphs four different cases are considered. 3.1. Semilinear parabolic PDEs: explicit case The following theorem holds for explicit case Theorem 1. Let us consider the semilinear equation (Eq.(12)) where F (U ) is a differentiable function. The non-linear explicit scheme is convergent if  ∂F (U ) PN    + i=1 |mi | − m0 < 0 ∂U (15) 2  ∆t <  P N  |m | + m − ∂F (U ) i=1

i

0

∂U

Proof. By considering the non-linear scheme together with the explicit expressions Eqs.(4) and (11), the following expression is obtained N

X un+1 − un0 0 = −m0 un0 + mi uni + F (un0 ) ∆t i=1

(16)

and let the same expression for the exact solution be N

X U0n+1 − U0n = −m0 U0n + mi Uin + F (U0n ) ∆t i=1

(17)

If we define the error as eni = uni − Uin and en0 = un0 − U0n , using the mean value theorem for F , it is obtained F (un0 ) − F (U0n ) = en0

∂F (χj ) ∂U

(18)

where χj = U0n + ξen0 , ξ ∈ [0, 1] Subtracting Eqs.(16) and (17) and introducing Eq.(18) it is obtained N

X − en0 en+1 ∂F 0 = −m0 en0 + mi eni + en0 (χj ) ∆t ∂U i=1

(19)

N

en+1 0

X ∂F = (1 − m0 ∆t + ∆t (χj ))en0 + ∆t mi eni ∂U i=1 6

(20)

N

|en+1 | ≤ |(1 − m0 ∆t + ∆t 0

X ∂F (χj ))||en0 | + ∆t| mi eni | ∂U i=1

(21)

Let en = max |eni |, taking into account Eqs.(4) and (1) it is obtained i=0,..,N

N

e

n+1

X ∂F ≤ e (|1 − m0 ∆t + ∆t |+ ∆t|mi |) + Θ1 (∆t(h2i + ki2 + ∆t)) (22) ∂U i=1 n

N

X ∂F α = |1 − m0 ∆t + ∆t |+ ∆t|mi | ∂U i=1

(23)

e1 = αe0 + C(∆t(h2i + ki2 + ∆t))

(24)

by considering n = 1 and including a constant C

as e0 = 0 by the initial condition Eq.(3) e1 = C(∆t(h2i + ki2 + ∆t))

(25)

and as in Eq.(24) e2 = αe1 + C(∆t(h2i + ki2 + ∆t)) = C(∆t(h2i + ki2 + ∆t))(1 + α) en+1 = C(∆t(h2i + ki2 + ∆t))(1 + α + α2 + · · · + αn )

(26) (27)

1 + α + α2 + · · · + αn + · · · is a geometric series, where the condition of convergence is |α| < 1. Hence, in order for the error not to diverge as time increment tends to 0, it must be: N

X ∂F |α| < 1 ⇔ |1 − m0 ∆t + ∆t |+ ∆t|mi | < 1 ∂U i=1

(28)

Thus, N

X ∂F |1 − m0 ∆t + ∆t |<1− ∆t|mi | ∂U i=1

−1 +

N X i=1

(29)

N

X ∂F <1− ∆t|mi | ∆t|mi | < 1 − m0 ∆t + ∆t ∂U i=1 7

(30)

by comparing first and second terms of inequality Eq.(30) and second and third terms the following two conditions can be obtained.   P P ∂F ∂F <1− N ∆t|mi | + N 1 − m0 ∆t + ∆t ∂U ∆t( ∂U i=1 i=1 |mi | − m0 ) < 0 PN PN ∂F ⇔ ∂F −1 + i=1 ∆t|mi | < 1 − m0 ∆t + ∆t ∂U ∆t( i=1 |mi | + m0 − ∂U )<2 (31) As ∆t > 0, the convergence conditions are  ∂F PN  + i=1 |mi | − m0 < 0  ∂U (32) 2   ∆t < PN ∂F i=1 |mi | + m0 − ∂U ∂F From first inequality of Eq.(31), taking account Eq.(11) thus ∂U < 0 is obtained. The second inequality of Eq.(32) gives us the value ∆t for convergence of each one of the stars of the domain. Then the minimum value obtained among all the stars is taken as value ∆t for convergence condition.

Remark 1. The method is stable because small perturbations of initial data lead to small variations in the solution at a later time. Let us consider a perturbation  in the initial instant e0 =  e1 = α + C(∆t(h2i + ki2 + ∆t)) e2 = αe1 + C(∆t(h2i + ki2 + ∆t)) = α2  + C(∆t(h2i + ki2 + ∆t))(1 + α) en+1 = αn+1  + C(∆t(h2i + ki2 + ∆t))(1 + α + α2 + · · · + αn )

As |α| < 1 → lim αn  → 0 n→∞

3.2. Semilinear parabolic PDEs: implicit case The following theorem holds for implicit case Theorem 2. Let us consider the semilinear equation (Eq.(12)) where F (U ) is a differentiable function. The non-linear implicit scheme is convergent

8

Proof. Let us consider Eq.(12) and the implicit scheme, the following expression is obtained N

X U0n+1 − U0n = −m0 U0n+1 + mi Uin+1 + F (U0n+1 ) ∆t i=1

(33)

Introducing the error as in the explicit case and expressing the difference of the values of the function F (U ) using the mean value theorem, it is obtained as in Eq.(20) ∂F en+1 (1+m0 ∆t−∆t 0

∂U

(χj )) =

en0 +∆t

N X

mi en+1 +Θ2 (∆t(h2i +ki2 +∆t)) (34) i

i=1

Let en = max |eni |, it is obtained i=0,..,N

N

X ∂F (χj )| ≤ en + en+1 ∆t |mi | + C1 (∆t(h2i + ki2 + ∆t)) ∂U i=1 (35) where C1 is a constant. Denoting

en+1 |1 + m0 ∆t − ∆t

N

X ∂F | − ∆t |mi | β = |1 + m0 ∆t − ∆t ∂U i=1 en+1 ≤ Thus e1 ≤

en C1 (∆t(h2i + ki2 + ∆t)) + β β

e0 C1 (∆t(h2i + ki2 + ∆t)) + β β

(36)

(37)

(38)

as e0 = 0 by initial condition Eq.(3), then as in Eq.(27) en ≤ C1 (∆t(h2i + ki2 + ∆t))[

1 1 1 + 2 + ··· + n] β β β

(39)

Taking into account Eq.(32) N

N

X ∂F X ∂F + |mi | − m0 < 0 ⇔ m0 − > |mi | > 0 ∂U ∂U i=1 i=1 9

(40)

hence β of Eq.(36) N N X ∂F ∂F X β = |1+m0 ∆t−∆t |−∆t |mi | = 1+∆t[m0 − − |mi |] > 1 (41) ∂U ∂U i=1 i=1

1 1 + 2 + · · · is convergent and taking β β into account Eq.(39) the error is bounded and, thus, the convergence of the implicit scheme is obtained for any time increment. then β > 1 and the geometric serie

3.3. Quasilinear (case a) parabolic PDEs: explicit case The following theorem holds for explicit case Theorem 3. Let us consider the quasilinear (case a) equation (Eq.(13)) where K(U ) is a differentiable function. The non-linear explicit scheme is convergent if  P  −R + m0 K(un0 ) > |K(un0 )| N i=1 |mi | 1 (42)  ∆t < m0 K(un0 ) − R where R =

P ∂K n (−m0 U0n + N i=1 mi Ui ) ∂U

Proof. Let us consider Eq.(13), where K(U ) is a differentiable function. Using fully explicit approximations N

X un+1 − un0 0 = K(un0 )(−m0 un0 + mi uni ) ∆t i=1

(43)

and similarly for the exact solution N

X U0n+1 − U0n = K(U0n )(−m0 U0n + mi Uin ) ∆t i=1

(44)

Defining the error en = un − U n , it is obtained N N X X e0n+1 − en0 n n n n n = K(u0 )(−m0 u0 + mi ui ) − K(U0 )(−m0 U0 + mi Uin ) (45) ∆t i=1 i=1

10

N X

en+1 = en0 + ∆t[K(un0 )(−m0 un0 + 0

i=1

+ K(un0 )(−m0 U0n +

N X

mi uni ) − K(un0 )(−m0 U0n +

mi Uin ) − K(U0n )(−m0 U0n +

i=1

N X

N X

mi Uin )

i=1

mi Uin )] (46)

i=1

and using the mean value theorem en+1 0

=

en0 + ∆t[K(un0 )(−m0 en0 +

N X

∂K (χj )(−m0 U0n + mi eni ) + en0 ∂U

i=1

N X

mi Uin )]

i=1

(47)

and rearranging en+1 0

=

∂K en0 [1+∆t (−m0 U0n + ∂U

If, again, we make e = max

i=0,..,N

shall obtain the relation: e

mi Uin )−m0 ∆tK(un0 )]+∆tK(un0 )

i=1

n

n+1

N X

|eni |

N X

eni mi

i=1

(48) in every node of the cloud of points, we

N N X X ∂K n n n n ≤ e (|1+∆t (−m0 U0 + mi Ui )−m0 ∆tK(u0 )|+∆t|K(u0 )| |mi |) ∂U i=1 i=1 n

+ Θ3 (∆t(h2i + ki2 + ∆t)) (49)

Denoting: N

X ∂K (−m0 U0n + mi Uin ) R= ∂U i=1 γ = |1 + R∆t −

Thus,

m0 ∆tK(un0 )|

+

∆t|K(un0 )|

(50) N X i=1

|mi |

(51)

en+1 ≤ en γ + Θ3 (∆t(h2i + ki2 + ∆t))

(52)

e1 ≤ e0 γ + C2 (∆t(h2i + ki2 + ∆t)) = C2 (∆t(h2i + ki2 + ∆t))

(53)

where e0 = 0 by initial condition Eq.(3) en ≤ C2 (∆t(h2i + ki2 + ∆t))(1 + γ + · · · ) 11

(54)

Taking into account that 1 + γ + γ 2 + · · · is a geometric series, then the condition of convergence is |γ| < 1. |1 + R∆t − m0 ∆tK(un0 )| + ∆t|K(un0 )| −1 + ∆t|K(un0 )|

N X i=1

N X i=1

|mi | < 1

|mi | < 1 + R∆t − m0 ∆tK(un0 ) < 1 − ∆t|K(un0 )|

by relating second and third terms and first and second ones  P 1 + R∆t − m0 ∆tK(un0 ) < 1 − ∆t|K(un0 )| N i=1 |mi | PN n −1 + ∆t|K(u0 )| i=1 |mi | < 1 + R∆t − m0 ∆tK(un0 )

then the following two inequalities are obtained  P n ∆t(R + |K(un0 )| N i=1 |mi | − m0 K(u0 )) < 0 P N ∆t(|K(un0 )| i=1 |mi | + K(un0 )m0 − R) < 2 and operating in Eq.(58)  PN  −R + m0 K(un0 ) > |K(un0 )| i=1 |mi | 2 P  ∆t < N |K(un0 )| i=1 |mi | + m0 K(un0 ) − R

(55) N X i=1

|mi | (56)

(57)

(58)

(59)

by taking first inequality of the Eq.(59) and substituting it in the second one  P  −R + m0 K(un0 ) > |K(un0 )| N i=1 |mi | 1 (60)  ∆t < n m0 K(u0 ) − R As it is shown in the second inequality of Eq.(58) the values of the coefficients (mo , mi ) are important to determine ∆t for each one of the stars. For each cloud of nodes the value ∆t, convergence condition, is the minimum value obtained considering all the stars of the domain. Remark 2. The method is stable because small perturbations of initial data lead to small variations in the solution at a later time. 12

Let us consider a perturbation  in the initial instant e0 =  e1 = γ + C(∆t(h2i + ki2 + ∆t)) e2 = γe1 + C(∆t(h2i + ki2 + ∆t)) = γ 2  + C(∆t(h2i + ki2 + ∆t))(1 + γ) en+1 = γ n+1  + C(∆t(h2i + ki2 + ∆t))(1 + γ + γ 2 + · · · + γ n )

As |γ| < 1 → lim γ n  → 0 n→∞

3.4. Quasilinear (case b) parabolic PDEs: explicit case The following theorem holds for explicit case Theorem 4. Let us consider the quasilinear (case b) equation (Eq.(13)) where F (U ) and K(U ) are differentiable functions. The non-linear explicit scheme is convergent if  PN ∂F ∂F  n n  −R − |m | ⇔ −R − + m K(u ) > |K(u )| >0  i 0 0 0 i=1  ∂U ∂U 1 (61) ∆t <   ∂F   m0 K(un0 ) − R − ∂U

where R =

P ∂K n (−m0 U0n + N i=1 mi Ui ) ∂U

Proof. Let us consider now Eq.(14), where F (U ) and K(U ) are differentiable functions, and operating as in previous subsections Eqs.(22), (49) and (50), the following inequality is obtained N

e

n+1

n

≤ e (|1 + R∆t −

m0 K(U0n )∆t

X ∂F + ∆t | + ∆t|K(U0n )| |mi |) (62) ∂U i=1

Denoting N

τ = |1 + R∆t −

m0 K(U0n )∆t

X ∂F + ∆t | + ∆t|K(U0n )| |mi | ∂U i=1

(63)

and operating as in previous subsections en ≤ C3 (∆t(h2i + ki2 + ∆t))(1 + τ + · · · ) 13

(64)

Taking into account that 1 + τ + τ 2 + · · · is a geometric series, then the condition of convergence is |τ | < 1. N

|1 + R∆t − −1+∆t|K(un0 )|

m0 K(U0n )∆t

N X i=1

|mi | <

X ∂F + ∆t | + ∆t|K(U0n )| |mi | < 1 ∂U i=1

∂F 1+R∆t−m0 K(U0n )∆t+∆t

∂U

<

(65)

1−∆t|K(un0 )|

N X i=1

(66) and relating second and third terms and first and second ones, the following two conditions are obtained  PN ∂F ∂F  n n    −R − ∂U + m0 K(u0 ) > |K(u0 )| i=1 |mi | ⇔ −R − ∂U > 0 1 (67) ∆t <   ∂F   m0 K(un0 ) − R − ∂U Then ∆t is bounded as it is shown in Eq.(67). For each cloud of nodes the value ∆t, convergence condition, is the minimum value obtained of all the stars of the domain. Remark 3. The method is stable because small perturbations of initial data lead to small variations in the solution at a later time. Let us consider a perturbation  in the initial instant e0 =  e1 = τ  + C(∆t(h2i + ki2 + ∆t)) e2 = τ e1 + C(∆t(h2i + ki2 + ∆t)) = τ 2  + C(∆t(h2i + ki2 + ∆t))(1 + τ ) en+1 = τ n+1  + C(∆t(h2i + ki2 + ∆t))(1 + τ + τ 2 + · · · + τ n )

As |τ | < 1 → lim τ n  → 0 n→∞

14

|mi |

4. Numerical results In this section it will be shown the numerical results obtained by solving five non-linear parabolic PDEs with Dirichlet boundary conditions the initial conditions (which are obtained by substitution in the exact solution), which correspond to non-linear problems in physics and engineering: heat equation, diffusion of liquids in porous medium, advection-diffusion, etc. In fig. 1 a variety of clouds has been used to show the possibilities of the method. These clouds of points are quasi-uniform except in the neighbourhood of some curved boundaries, however cloud 4 of fig. 1 is irregular. Also nodes distribution, based in convex quadtree cells as proposed in [11], can be recommended. Firstly we consider 3 examples and apply them in 6 different discretizations (see figure 1), obtaining the solution for t = 0.25 sec. To do it we divide the total time into a different number of steps, in particular, 50, 100, 250, 500 and 1000 steps, which is equivalent to use as time increments 0.005, 0.0025, 0.001, 0.0005 and 0.00025 sec. respectively. In all the cases considered the global error has been calculated according to the norm v u NI  2 uX u sol(i) − exac(i) u t i=1 NI |exacmax |

%Global Error =

× 100

(68)

where sol(i) is the GFD solution in node i, exac(i) is the exact value of the solution at node i, exacmax is the maximum value of the exact values in the cloud of nodes considered and N I is the number of nodes of the domain Ω. Notice that in the numerical results the time step used, ∆t, has been selected according to the convergence condition defined in Eqs.(15), (42) and (61). However these values depend on the stars of nodes, and hence they are different for the six clouds of points considering in fig.1 Then as consequence only the converging cases have been considered in tables [1-5]. Hence the non convergent cases for such ∆t are left blank in tables. Example 1 Let us consider the equation ∂U ∂ 2U ∂ 2U = + − U ln(U ) ∂t ∂x2 ∂y 2 15

(69)

Figure 1: Clouds of points 1, 2, 3, 4 , 5 and 6

with exact solution: −t

x2

xy

y2

1

U (x, y, t) = ee e 5 + 5 + 20 + 2

(70)

This example corresponds to a semilinear case, Eq.(12), whose convergence has been studied in previous section. Time increments, ∆t, 16

calculated using the convergence condition (Eq. 15) for each cloud of nodes are, respectively: 0.0022, 0.00125, 0.001, 0.0078, 0.00125 and 0.0022. Thus, the numerical results obtained Table 1 and Figures. 2 and 3 are according to the convergence condition previously developed. mesh 4t 0.005 1 2 3 4 3.8457 · 10−4 5 6 -

0.0025 1.6820 · 10−4 -

0.001 2.7582 · 10−5 2.7802 · 10−3 5.5248 · 10−5 1.2141 · 10−4 1.6033 · 10−5 4.1119 · 10−5

0.0005 8.2719 · 10−6 2.0688 · 10−5 2.7699 · 10−5 4.4682 · 10−1 1.0019 · 10−5 2.3447 · 10−5

0.00025 2.4258 · 10−5 6.8425 · 10−6 2.5353 · 10−5 1.5638 · 10−4 9.2320 · 10−6 5.3133 · 10−5

Table 1: Table of global errors Eq.(68) for Example 1

Figure 2: Exact, approximated solution and error of the example 1 on cloud of points 1

17

Figure 3: Exact, approximated solution and error of the example 1 on cloud of points 4

18

Example 2 Let us consider the equation ∂U ∂ 2U ∂ 2U = U( 2 + ) ∂t ∂x ∂y 2

(71)

with the exact solution: U (x, y, t) =

1 x2 + y 2 + 1 ( + xy + x + y) 1−t 4

(72)

This case corresponds to a heat transmission problem where the conductivity K depends on the temperature U . This example corresponds to a quasilinear (case a), Eq.(13), whose convergence has been studied in previous section. Time increments, ∆t, calculated using the convergence condition (Eq. 42) for each cloud of nodes being, respectively: 0.00088, 0.00035, 0.00067, 0.0038, 0.00076 and 0.00088. Thus, the numerical results obtained Table 2 and Figures 4 and 5 are according to the convergence condition, previously developed. mesh 4t 0.005 0.0025 1 2 3 4 1.5380 · 10−4 5 6 -

0.001 6.1535 · 10−5 -

0.0005 2.1974 · 10−5 9.8553 · 10−2 3.0915 · 10−5 9.8553 · 10−2 3.8040 · 10−5

Table 2: Table of global errors Eq.(68) for Example 2

19

0.00025 1.1149 · 10−5 1.3456 · 10−5 1.2612 · 10−5 1.5726 · 10−5 1.2612 · 10−5 1.9095 · 10−5

Figure 4: Exact, approximated solution and error of the example 2 on cloud of points 2

20

Figure 5: Exact, approximated solution and error of the example 2 on cloud of points 4

21

Example 3 Let us consider the equation ∂U ∂ 2U ∂ 2U = U( 2 + 2 ) + U ∂t ∂ x ∂ y

(73)

with the exact solution: U (x, y, t) = et (1 + ex sin(y))

(74)

This example corresponds to a quasilinear (case b), Eq.(14), whose convergence has been studied in previous section. Time increments, ∆t, calculated using the convergence condition (Eq. 61) for each cloud of nodes are, respectively: 0.0017, 0.00065, 0.00062, 0.0031, 0.00072 and 0.0017. Thus, the numerical results obtained Table 3 and Figures 6 and 7 are according to the convergence condition, previously developed. mesh ∆t 1 2 3 4 5 6

0.005 0.0025 2.1002 · 10−4 -

0.001 2.2895 · 10−1 1.9940 · 10−4 5.2962 · 10−1

0.0005 1.3241 · 10−5 1.497 · 10−1 4.9862 · 10−5 1.9655 · 10−4 3.8149 · 10−1 2.838 · 10−5

Table 3: Table of global errors Eq.(68) for Example 3

22

0.00025 1.0168 · 10−5 1.3643 · 10−5 4.7147 · 10−5 1.9527 · 10−4 2.1619 · 10−5 2.2815 · 10−5

Figure 6: Exact, approximated solution and error of the example 3 on cloud of points 3

23

Figure 7: Exact, approximated solution and error of the example 3 on cloud of points 4

24

Secondly we consider now another three non-linear partial differential equations corresponding to physical problems: Non-linear heat equation, Burgers’ equation and porous medium equation. 4.1. Non-linear heat equation Let us consider the equation ∂ 2U ∂U 2 ∂U ∂ 2U ∂U 2 = U( 2 + ) +( ) )+( 2 ∂t ∂x ∂y ∂x ∂y

(75)

with the exact solution U (x, y, t) = −

2 1 y2 3 + x(1 + t)− 3 + (1 + t) 3 6(1 + t) 2

(76)

Figure 8: Non-linear heat equation, exact and approximated solutions of case 4.1 on clouds of points 4

25

Figure 9: Non-linear heat equation, exact and approximated solutions of case 4.1 on clouds of points 5

mesh 4t 0.005 0.0025 1 2 3 4 4.8093 · 10−6 5 6 -

0.001 3.3081 · 10−6 -

0.0005 2.9100 · 10−6 6.6274 · 10−2 3.4911 · 10−6 2.9993 · 10−6 2.9163 · 10−6 2.9699 · 10−6

0.00025 2.8894 · 10−6 2.9067 · 10−6 3.3823 · 10−6 2.9016 · 10−6 2.9205 · 10−6 2.8860 · 10−6

Table 4: Table of global errors Eq.(68) for non-linear heat equation (case 4.1)

. This case corresponds to the complete heat equation where the conductivity depends on the temperature. The results obtained are shown in Table 4 and Figures 8 and 9.

26

4.2. Burgers’ equation Burgers’ equation is a fundamental partial differential equation occurring in various areas of applied mathematics, such as fluid mechanics, nonlinear acoustics, gas dynamics, traffic flow. This equation corresponds a semilinear case, in the 2-D case is defined as ∂U ∂U ∂U ∂ 2U ∂ 2U + U( + )= + ∂t ∂x ∂y ∂x2 ∂y 2

(77)

We consider the case correspondents to the following exact solution: U (x, y, t) =

1 1−x−y ( ) 1−t 2

(78)

The results obtained are given in Table 5 and Figures 10 and 11.

Figure 10: Burgers’ equation, exact and approximated solutions on clouds of points 4

27

Figure 11: Burgers’ equation, exact and approximated solutions on clouds of points 6

mesh 4t 0.005 1 2 3 4 2.5466 · 10−5 5 6 -

0.0025 1.2433 · 10−5 -

0.001 7.1157 · 10−6 1.7489 · 10−5 5.5798 · 10−6 4.9985 · 10−6 4.0343 · 10−6 7.6230 · 10−6

0.0005 4.2539 · 10−6 4.2360 · 10−6 3.6631 · 10−6 3.0502 · 10−6 3.1601 · 10−6 4.4709 · 10−6

Table 5: Table of global errors Eq.(68) for Burgers’ equation

28

0.00025 3.1617 · 10−6 3.1731 · 10−6 2.9980 · 10−6 2.5682 · 10−6 2.8694 · 10−6 3.2399 · 10−6

4.3. Porous medium equation In this case, we consider the two-dimensional porous medium equation (PME) (quasilinear case b). In this test which is taken from [16, 17, 18], ∂U ∂ ∂ = (U 2 ) + 2 (U 2 ) 2 ∂t ∂x ∂y with the initial data u(x, y, 0) given by two bumps  −1  ) , (x − 2)2 + (y + 2)2 < 6 exp(    6 − (x − 2)2 − (y + 2)2 −1 ) , (x + 2)2 + (y − 2)2 < 6 exp(  2 − (y − 2)2  6 − (x + 2)   0 , o.w.

(79)

(80)

The computational domain [−8, 8] × [−8, 8] is discretized into 81 × 81 nodes and 4t = 0.01. Fig.12 shows the numerical approximation at time t = 0, t = 0.5, t = 1.5, t = 3 and t = 4.5 that are the ones used in [16, 17, 18]. Figure 12 compares quite well with the figures given in [16, 17, 18]. As we can see from Fig.12, the symmetries of the initial condition are preserved. Moreover in Table 6 the maximum values obtained at the different time values including the point location are given. Point T=0.0 T=0.5 T=1.5 T=3.0 T=4.5

(2,-2) 0.8465 0.6872 0.4821 0.3623 0.3048

(-2,2) 0.8465 0.6872 0.4821 0.3623 0.3048

Table 6: Points with maximum value for each time value

29

Figure 12: From top left to bottom right, we show the numerical solution at times T = 0.0, 0.5, 1.5, 3.0, 4.5.

30

5. Conclusions In this paper explicit and implicit methods using GFDM for solving different non-linear parabolic PDEs, have been considered. Convergence has been studied, for semilinear and quasilinear equations, and limits of convergence have been developed and implemented in different examples. The examples provided illustrates the viability of the application of GFDM for solving parabolic non-linear PDEs in 2D. The efficiency of the developed methods is clearly shown. The accuracy of the GFDM has been tested in different non-linear PDEs, including different cases related with acoustics, heat transfer, mass transfer, heat extinction, combustion. Numerical results for several non-linear problems validate the use of GFDM to solve this type of practical problems. 6. Ackonwledgements The authors acknowledge the support of the Escuela T´ ecnica Superior de Ingenieros Industriales (UNED) of Spain, project 2017-IFC02, and of the Universidad Polit´ ecnica de Madrid (UPM) (Research groups 2017). 7. References [1] A. D. Polyanin, V. F. Zaitsev, Handbook of nonlinear partial differential equations, Chapman & Hall/CRC, ISBN 1-58488-355-3. [2] A. Tadmor, A review of numerical methods for non-linear partial differential equations, Bulletin of the American Mathematical Society 2012; 42(4):507-554. [3] P. S. Jensen, Finite difference technique for variable grids, Computer& Structures 1972; 2: 17–29. [4] T. Liszka, J. Orkisz, The Finite Difference Method at Arbitrary Irregular Grids and its Application in Applied Mechanics, Computer & Structures 1980; 11: 83–95. [5] J. Orkisz, Finite Difference Method (Part, III) in handbook of Computational Solid Mechanics, M. Kleiber (Ed.), Spriger-Verlag, Berlin 1998.

31

[6] N. Perrone, R. Kao, A general finite difference method for arbitrary meshes, Computer& Structures 1975; 5: 45–58. [7] J. J. Benito, F. Ure˜ na, L. Gavete, Influence of several factors in the generalized finite difference method, Applied Mathematical Modeling 2001; 25: 1039–1053. [8] 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 2011; 235: 1849–1855 [9] 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 2007; 209: 208–233. [10] L. Gavete, J. J. Benito, F. Ure˜ na, Generalized finite differences for solving 3D elliptic and parabolic eqwuations, Applied Mathematical Modelling 2016; 40: 955–965. [11] L. Gavete, M. L. Gavete, F. Ure˜ na, J. J. Benito, An Approach to Refinement of Irregular Clouds of Points Using Generalized Finite Differences, Mathematical Problems in engineering, DOI:10.1155/2015/283757, Vol. 2015 (2015), 9 pages. [12] 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 2013; 37: 1189–1196. [13] T. Zhang, Y. F. Ren, C. M. Fan, P. W. Li, Simulation of two-dimensional sloshing phenomenon by generalized finite difference method, Engineering Analysis with Boundary Elements 2016; 63: 82–91 [14] B. Mochnacki, E. Majchrzak, Numerical modelling of casting solidification using generalized finite difference method, Materials Science Forum 2010; 638-642: 2676–2681. [15] L. Gavete, F. Ure˜ na, J. J. Benito, A. Garcia, M. Ure˜ na, E. Salete, Solving second order non-linear elliptic partial differential equations using generalized finite difference method, Journal of Computational and Applied Mathematics 2017; 318: 378–387. 32

[16] F. Cavalli, G. Naldi, G. Puppo, M. Semplice, High order relaxation schemes for nonlinear degenerate diffusion problems, SIAM J. Numer. Anal. 2007; 45: 2098–2119. [17] Y. Liu, C.-W. Shu, M. Zhang, High order finite difference WENO schemes for nonlinear degenerate parabolic equations, SIAM J. Sci. Comput. 2011; 33: 939–965. [18] R. Abedian, H. Adibi, M. Dehghan, A high-order weighted essentially non-oscillatory (WENO) finite difference scheme for nonlinear degenerate parabolic equations, Computer Physics Communications 2013; 184:1874-1888 [19] P. Lancaster, K. Salkauskas, Curve and surface fitting, Ed. Acasemic Press, 1986. [20] D. Levin, The approximation power of moving least squares, Math. Comp. !998; 67,224:1517-1531.

33