- Email: [email protected]
Solving second order non-linear hyperbolic PDEs using generalized finite difference method (GFDM)
Journal of Computational and Applied Mathematics 363 (2020) 1–21
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Solving second order non-linear hyperbolic PDEs using generalized finite difference method (GFDM) ∗
F. Ureña a , , L. Gavete b , A. García a , J.J. Benito a , A.M. Vargas c a b c
UNED, ETSII, Madrid, Spain UPM, ETSIM, Madrid, Spain UCM, Madrid, Spain
article
info
Article history: Received 17 May 2018 Received in revised form 1 March 2019 Keywords: Meshless methods Generalized finite difference method Non-linear hyperbolic partial differential equations
a b s t r a c t 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 hyperbolic PDEs. This paper illustrates that the GFD explicit formulae developed to obtain the different derivatives of the PDEs are based on the existence of a positive definite matrix which is obtained using moving least squares approximation and Taylor series development. Consistency and stability are shown in this paper for semilinear and quasilinear hyperbolic equations. This paper shows the application of the GFDM for solving second order non-linear hyperbolic problems. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Modern numerical methods, in particular those for solving non-linear PDEs, have been developed in recent years using finite differences, finite elements, finite volumes or spectral methods. A review of numerical methods for non-linear partial differential equations is given by Polyanin [1] and Tadmor [2]. Nonlinear hyperbolic partial differential equations have been applied in different fields such as in hypoelastic solids [3], astrophysics [4], electromagnetic theory [5], propagation of heat waves [6] and other disciplines. In this paper we use a meshless method called generalized finite difference method (GFDM) for solving different hyperbolic non-linear PDEs. Researchers as Jensen [7], Liszka and Orkisz [8], Orkisz [9] and Perrone and Kao [10] have contributed to the development of the GFDM in different aspects of its applications. Benito, Gavete and Ureña [11–13] have studied the influence of several factors and developed the explicit formulae and h-adaptive method for the solution of the PDEs in 2D. GFDM was also used to solve both non-linear elliptic [14] and parabolic equations [15]. Various applications of the GFDM have been developed for solving different inverse problems by Fan et al. [16,17]. Hosseini [18,19] has applied the GFDM for natural frequency analysis of nanocomposite cylinders and for shock-induced two dimensional coupled non-Fickian diffusion-elasticity analysis. ∗ Corresponding author. E-mail addresses: [email protected] (F. Ureña), [email protected] (L. Gavete), [email protected] (A. García), [email protected] (J.J. Benito), [email protected] (A.M. Vargas). https://doi.org/10.1016/j.cam.2019.05.028 0377-0427/© 2019 Elsevier B.V. All rights reserved.
2
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
The GFDM has also been used to solve 3D partial differential equations by Izadian et al. [20], Gavete et al. [21] and Hua et al. [22]. The GFDM was also applied for different practical tasks as in [23] for solving two dimensional non-linear obstacle problems, in [24,25] for simulating two dimensional sloshing phenomena and for the propagation of nonlinear water waves, in [26] for the modelling of casting solidification, and in [27] for solving two-dimensional shallow water equations. Adaptive methods for solving elliptical problems with singularities using GFD have been proposed [28,29]. An adaptive Particle-in-Cloud method has been used for solving the Vlasov–Poisson equation [30] and a Lagrangian particle method for solving Euler equations [31]. Both methods use GFD to improve the approximation of the differential operators involved. This paper shows that the GFDM can be used to solve second order hyperbolic nonlinear PDEs with different irregular clouds of points in 2D. Hyperbolic equations have been solved using an explicit method and their consistency and stability have been studied taking into account the irregularity of the cloud of points. The numerical results show the high accuracy obtained. A brief outline of this paper is as follows. In Section 2 the explicit method and GFDM to solve non-linear hyperbolic partial differential equations are shown. In Section 3, consistency and stability are studied. In Section 4 the proposed method is applied to different test problems to show the high accuracy of the method. In the last section, we present some conclusions. 2. Explicit method and GFDM: application to non-linear hyperbolic partial differential equations Consider the following non-linear problem in the domain D = [0, T ] × Ω with Ω ⊂ R2
∂ 2U = LΩ [ U ] ∂t2 where LΩ [U ] is a non-linear operator, Γ is the boundary of the domain Ω , with the boundary condition UΓ = f (t)
(1)
(2)
and initial conditions U(x, y, 0) = g(x, y)
(3)
∂U (x, y, 0) = h(x, y) ∂t
(4)
where f , g and h are three known functions. To solve the problem described by Eqs. (1)–(4), using the explicit method, the time derivative is approximated by un0+1 − 2un0 + un0−1 ∂ 2 u(x0 , y0 , n∆t) = + Θ ((∆t)2 ) ∂t2 (∆t)2
(5)
and spatial derivatives are approximated using GFD [13], as described in the following subsection. 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 Ω , 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 [11] 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: Ui = U0 + hi
( ) 2 ∂ U0 1 2 ∂ 2 U0 ∂ 2 U0 ∂ U0 2 ∂ U0 + ki + hi + k + 2h k + ···, i i i ∂x ∂y 2 ∂ x2 ∂ y2 ∂ x∂ y
i = 1, . . . , s
(6)
We define: cTi = {hi , ki , DT = {
h2i k2i
, , hi ki }, 2 2 ∂ u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0
∂ u0 , , , , }. ∂ x ∂ y ∂ x2 ∂ y2 ∂ x∂ y
(7) (8)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
3
If in Eq. (6) the higher than second order terms are ignored, a second order approximation for the function Ui is obtained. We denote it by ui . It is then possible to define the function B(u) as s ∑ ∂ u0 ∂ u0 [(u0 − ui ) + hi + ki + ∂x ∂y
B(u) =
i=1
∂ 2 u0 2 2 ∂ 2 u0 ∂ 2 u0 + k2i + 2hi ki )] wi 2 2 ∂x ∂y ∂ x∂ y
1
+ (h2i 2
(9)
where wi = w (hi , ki ) are positive symmetrical weighting functions, having the property, as defined in [32], of decreasing monotonically in magnitude as the distance to the centre of the weighting function increases. See also [33]. 1 Some weighting functions as potential or exponential exp(−n(dist 2 )) can be used [11], where n ∈ N. If the norm dist n given by Eq. (9) 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 )
(10)
where h1 k1
⎛ ⎜
h2 k2
.. .
··· ··· .. .
hs ks
h2 k2
···
hs ks
.. .
A=⎜ ⎝
h1 k1
⎞⎛
ω12
ω22
⎟⎜ ⎟⎝ ⎠
.. .
··· ω
2 s
⎞⎛ ⎟⎜ ⎠⎜ ⎝
h1 h2
.. .
hs
k1 k2
.. .
··· ··· .. .
h1 k1 h2 k2 ⎟
ks
···
hs ks
⎞
.. .
⎟ ⎠
(11)
and s ∑
(−u0 + ui )hi wi2 ,
bT = {
i=1
s ∑
(−u0 + ui )ki wi2 ,
i=1
s ∑
( − u0 + ui )
i=1
k2i 2
wi2 ,
s ∑
( − u0 + ui )
i=1
h2i 2
wi2 ,
s ∑
(−u0 + ui )hi ki wi2 }.
i=1
The matrix A is a positive definite matrix and the approximation is of second order Θ (h2i , k2i ) [14]. The assumption stated on the nodes distribution is related with the criterion selected to choose the stars nodes, because in each one of the stars, the matrix A must be positive definite to obtain the approximate derivatives, hence the existence of ill-conditioned stars of nodes is restricted. If we define A−1 = Q Q T
(12)
D = Q Q T b.
(13)
then
Eq. (13) can be rewritten as D = −u0 Q Q T
s ∑
wi2 c i + Q Q T
s ∑
i=1
ui wi2 c i
(14)
i=1
or D = Q Q T W (u − u0 1)
(15)
where h1 w12
h2 w22
⎜ k w2 ⎜ 1 1 ⎜ 2 ⎜ h1 2 ⎜ W = ⎜ 2 w1 ⎜ 2 ⎜ k ⎜ 1 w2 ⎝ 2 1
k2 w22
⎛
h1 k1 w12
h22 2 k22
w22
··· ··· .. .
w22
.. .
2 h2 k2 ws2
···
hs ws2 ks ws2 h2s 2 k2s
ws2 w2
2 s hs ks ws2
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
and 1 = {1, 1, . . . , 1} ;
u = {u1 , u2 , . . . , us }T .
4
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Lemma 1. The truncation error of the partial derivatives is zero when (hi , ki ) → (0, 0). Proof. Taking into account Eq. (15) D = Q Q T W (U − U0 1) =
⎛⎛ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ T ⎜ QQ W ⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎝⎝
⎛
∂ ∂ + k1 )U0 + ∂x ∂y ∂ ∂ U0 + (h2 + k2 )U0 + ∂x ∂y U0 + (h1
U0 + (hs
h1
⎜ h2 QQTW ⎜ ⎝ .. . hs
k1 k2
∂ ∂x
.. .
··· ··· .. .
ks
···
∂ ∂ + k1 )2 )U0 + · · · 2! ∂x ∂y 1 ∂ ∂ (h2 + k2 )2 )U0 + · · · 2! ∂x ∂y .. . ∂ 1 ∂ ∂ + ks )U0 + (hs + ks )2 )U0 + · · · ∂y 2! ∂x ∂y ⎛ ∂U ⎞ 0 ⎜ ∂x ⎟ ⎜ ∂U ⎟ 0 ⎟ ⎜ ⎟ ⎞⎜ h1 k 1 ⎜ ∂y ⎟ ⎟ ⎜ 2 h2 k 2 ⎟ ⎜ ∂ U 0 ⎟ T ⎟ ⎜ .. ⎟ ⎠ ⎜ ∂ x2 ⎟ + Q Q W R . ⎟ ⎜ 2 ⎜ ∂ U0 ⎟ hs k s ⎟ ⎜ ⎜ ∂ y2 ⎟ ⎟ ⎜ ⎝ ∂ 2U ⎠ 1
(h1
⎞ ⎟ ⎛ ⎟ ⎟ ⎟ ⎜ ⎟−⎜ ⎟ ⎝ ⎟ ⎟ ⎠
⎞ U0 ⎟ ⎟ U0 ⎟⎟
⎞⎟
⎟ .. ⎟ ⎠⎟ = . ⎟ ⎟ U0 ⎠
(16)
0
∂ x∂ y where
⎛ ⎜ ⎜ ⎜ ⎜ R=⎜ ⎜ ⎜ ⎜ ⎝
∂ ∂ + k1 )3 )U0 + · · · 3! ∂x ∂y ∂ ∂ 1 (h2 + k2 )3 )U0 + · · · 3! ∂x ∂y .. . 1 ∂ ∂ (hs + ks )3 )U0 + · · · 3! ∂x ∂y 1
(h1
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
thus Q Q T W (U − U0 1) = Q Q T AD + Q Q T WR = A−1 AD + Q Q T WR = D + Q Q T WR
(17)
Hence truncation error for spatial partial derivatives is TE = Q Q T WR
(18)
and therefore lim
(hi ,ki )→(0,0)
TE → 0. □
(19)
Remark 1. The explicit expressions of TE may be seen in [13]. Thus, spatial derivatives using GFD [11–13,21] are denoted by
⎧ s s ∑ ∑ ⎪ ∂ u(x0 , y0 , n△t) ⎪ 2 2 ⎪ = −λ u + λ u + Θ (h , k ) , w ith λ = λi ⎪ 0 0 i i 0 i i ⎪ ∂x ⎪ ⎪ i = 1 i = 1 ⎪ ⎪ s s ⎪ ∑ ∑ ⎪ ∂ u(x0 , y0 , n△t) ⎪ 2 2 ⎪ , k ) , w ith µ = µi = −µ u + µ u + Θ (h ⎪ 0 0 i i 0 i i ⎨ ∂y i=1 i=1 s ∑ ⎪ ∂ 2 u(x0 , y0 , n△t) ∂ 2 u(x0 , y0 , n△t) ⎪ ⎪ + = − m u + mi ui + Θ (h2i , k2i ) ⎪ 0 0 ⎪ ∂ x2 ∂ y2 ⎪ ⎪ i=1 ⎪ ⎪ s ⎪ ∑ ⎪ ⎪ ⎪ mi ⎪ ⎩ with m0 = i=1
(20)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
5
3. Consistency and stability of the scheme for non-linear hyperbolic PDEs In this section consistency and stability of non-linear hyperbolic PDEs, using GFDM, are studied. We will do so by starting from the following definitions:
• A partial differential equation is semilinear if the coefficients of its highest derivatives are functions of the spaces 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 case LΩ [U ] =
∂ 2U ∂ 2U + + F (U). ∂ x2 ∂ y2
(21)
Quasilinear case LΩ [U ] = c(U)(
∂ 2U ∂ 2U + 2 ) + F (U). 2 ∂x ∂y
(22)
Theorem 1. Let us consider the quasilinear equation (Eq. (22)) where c(U) and F (U) are two differentiable functions. The non-linear explicit scheme u0n+1 − 2un0 + un0−1
∆t 2
= c(un0 )(−m0 un0 +
s ∑
mj unj ) + F (un0 )
(23)
j=1
is consistent. Proof. The truncation error for time partial derivatives is TE = −
∆t 2 ∂ 4 U + Θ1 ((∆t)4 ) 12 ∂ t 4
(24)
and for spatial partial derivatives, taking into account Eqs. (18) and (19), is TE = c(U )Q Q T WR
(25)
where
(
c(U ) =
0
0
c(U)
c(U)
0
)
.
The truncation error of the scheme is the sum of Eqs. (24) and (25), then TE = −
∆t 2 ∂ 4 U + Θ1 ((∆t)4 ) + c(U )Q Q T WR . 12 ∂ t 4
(26)
Taking into account Lemma 1, we obtain lim
(hi ,ki ,∆t)→(0,0,0)
TE → 0.
(27)
Thus, the non-linear explicit scheme for quasilinear and semilinear equations is consistent.
□
Notice that if in Eq. (22) we make c(U) = 1, the proof is valid for the semilinear case. 3.1. Semilinear hyperbolic PDEs Theorem 2. Let us consider the semilinear equation (Eq. (21)) where F (U) is a differentiable function. The non-linear explicit scheme un0+1 − 2un0 + un0−1
∆t 2 is conditionally stable.
= −m0 un0 +
s ∑ j=1
mj unj + F (un0 )
(28)
6
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Proof. Consider the non-linear scheme un0+1 − 2un0 + un0−1 (∆t)2
= −m0 un0 +
N ∑
mj unj + F (un0 )
(29)
j=1
and let the same expression for the exact solution be U0n+1 − 2U0n + U0n−1 (∆t)2
= −m0 U0n +
N ∑
mi Ujn + F (U0n ).
(30)
j=1
If we define the error as enj = unj − Ujn and en0 = un0 − U0n , using the Mean Value Theorem for F , it is obtained
∂F (χ j ) ∂U ξ ∈ [0, 1]. Subtracting Eqs. (29) and (30) together with Eq. (31) it is obtained
F (un0 ) − F (U0n ) = en0 where χj = U0n + ξ en0 ,
e0n+1 − 2en0 + e0n−1 = (∆t)2 (−m0 en0 +
N ∑
mj enj + en0
j=1
∂F (χj )). ∂U
(31)
(32)
Let us consider the von Neumann criterion for stability analysis: en0 = ξ n eiα⃗ ·r⃗0 ;
enj = ξ n eiα⃗ ·⃗rj ; r⃗j = r⃗0 + δ⃗rj
(33)
where: ξ is called the amplification factor, α ⃗ is the vector of the wave numbers, r⃗0 is the vector of coordinates of the central node of the star and r⃗j is the vector of coordinates of the other nodes of the star. Considering Eq. (20) and including Eq. (33) in Eq. (32), it is obtained
ξ n+1 eiα⃗ ·r⃗0 − 2ξ n eiα⃗ ·r⃗0 + ξ n−1 eiα⃗ ·r⃗0 = = (∆t)2 [−m0 ξ n eiα⃗ ·r⃗0 + ξ n eiα⃗ ·r⃗0
N ∑
mj eiα⃗ ·δ rj + ξ n eiα⃗ ·r⃗0 ⃗
j=1
∂F (χj )]. ∂U
(34)
Eq. (34) is simplified by cancelling the term ξ n−1 eiα⃗ ·r⃗0 , which leads to
ξ 2 − 2ξ + 1 = (∆t)2 [−
N ∑
mj (1 − eiα⃗ ·δ rj ) + ⃗
j=1
∂F (χj )]ξ ∂U
(35)
or rearranging
ξ 2 − 2bξ + 1 = 0
(36)
where b=1−
N (∆t)2 ∑
2
[
mj (1 − eiα⃗ ·δ rj ) − ⃗
j=1
∂F (χj )] ∂U
(37)
Eq. (36) has two complex solutions: ξ1 and ξ2 with
{
√ ξ1 = b + √b2 − 1 ξ2 = b − b2 − 1
(38)
Stability condition implies
|ξ1 | = |ξ2 | = 1 ⇔ |b| ≤ 1
(39)
Denoting ϕj = α ⃗ · δ⃗rj , it is possible to write:
⎧ ⎨ ξ1 = eiϕj ξ = e−iϕj ⎩ ξ2 + ξ = 2b = 2 cos ϕ 1 2 j ∑N and considering j=1 mj sin(ϕj ) = 0, then the stability condition is −1 ≤ 1 −
N ∆t 2 ∑
2
(
j=1
mj (1 − cos(ϕj )) −
∂F (χj )) ≤ 1. ∂U
(40)
(41)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
7
Fig. 1. Clouds 1, 2, 3.
Table 1 Table of the computing time steps ∆τ and optimal time steps ∆t for Example 1 and clouds of points (1, . . . , 9). Mesh
1
2
3
4
5
6
7
8
9
∆t ∆τ
0,0672 0,0625
0,0341 0,03125
0,0162 0,015625
0,0659 0,0625
0,0336 0,03125
0,0165 0,015625
0,1294 0,125
0,0631 0,0625
0,0322 0,03125
If we relate the terms of the inequality Eq. (41), it is obtained N ∆t 2 ∑
2
(
mj (1 − cos(ϕj )) −
j=1
∂F (χj )) ≤ 2. ∂U
(42)
The inequality of Eq. (42), taking into account Eq. (20), may be written as (∆t)2 ≤
4 2 |m0 | −
∂F (χj ) ∂U
.
(43)
Eq. (43) is the condition for stability of semilinear case. □ The minimum value obtained among all the stars is taken as the value ∆t in the stability condition. The denominator of the expression (43) must be positive, which requires us to demand that, since m0 is positive, either ∂F 1 ∂F F is a decreasing function or 2m0 > . However, considering that m0 (m0 ∼ 2 ) and that is bounded, as F is a ∂U ∂U hi ∂F differentiable function, one could always choose hi so that 2m0 > . ∂U
8
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 2. Clouds 4, 5, 6.
Fig. 3. Clouds 7, 8, 9.
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 4. Example 1 on clouds of points of Figs. 1–3.
Table 2 Table of errors of l2 and l∞ norms for Example 1 and clouds of points 1, 2 and 3. m0
l2
l∞
442 1629 8147
0,0003612 0,0001006 0,0000201
0,0006701 0,0001815 0,0000379
Table 3 Table of errors of l2 and l∞ norms for Example 1 and clouds of points 4, 5 and 6. m0
l2
l∞
119 502 1928
0,000689 0,000165 0,0000432
0,000892 0,000263 0,0000793
9
10
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 5. Exact, approximate solution, absolute values of the local errors of Example 1 on cloud of points 3 of Fig. 1. Table 4 Table of errors of l2 and l∞ norms for Example 1 and clouds of points 7, 8 and 9. m0
l2
l∞
460 1757 7346
0,000284 0,0000782 0,0000193
0,000502 0,000181 0,0000338
Table 5 Table of the computing time steps ∆τ and optimal time steps ∆t for Example 2 and clouds of points (1, . . . , 9). Mesh
1
2
3
4
5
6
7
8
9
∆t ∆τ
0,0275 0,025
0,0143 0,0125
0,0064 0,00625
0,0269 0,025
0,0138 0,0125
0,0067 0,00625
0,0529 0,050
0,0258 0,025
0,0131 0,0125
Table 6 Table of errors of l2 and l∞ norms for Example 2 and clouds of points 1, 2 and 3. m0
l2
l∞
442 1629 8147
0,000139 0,0000391 0,00000792
0,000639 0,000183 0,000038
3.2. Quasilinear hyperbolic PDEs Theorem 3. Let us consider the quasilinear equation (Eq. (22)) where c(U) and F (U) are two differentiable functions. The non-linear explicit scheme un0+1 − 2un0 + un0−1
∆t 2 is conditionally stable.
= c(un0 )(−m0 un0 +
s ∑ j=1
mj unj ) + F (un0 )
(44)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
11
Fig. 6. Example 2 on clouds of points of Figs. 1–3. Table 7 Table of errors of l2 and l∞ norms for Example 2 and clouds of points 4, 5 and 6. m0
l2
l∞
119 502 1928
0,001082 0,000263 0,0000691
0,001481 0,0003488 0,0000902
Table 8 Table of errors of l2 and l∞ norms for Example 2 and clouds of points 7, 8 and 9. m0
l2
l∞
460 1757 7346
0,000232 0,0000613 0,0000149
0,000356 0,0000923 0,0000235
Proof. We proceed as in the proof of Theorem 2. Let us consider the scheme for the exact solution U0n+1 − 2U0n + U0n−1
∆t 2
= c(U0n )(−m0 U0n +
s ∑ j=1
mj Ujn ) + F (U0n ).
(45)
12
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 7. Exact, approximate solution, absolute values of the local errors of Example 2 on cloud of points 6 of Fig. 2.
Table 9 Table of the computing time steps ∆τ and optimal time steps ∆t for Example 3 and clouds of points (1, . . . , 9). Mesh
1
2
3
4
5
6
7
8
9
∆t ∆τ
0,0672 0,0625
0,0341 0,03125
0,0162 0,015625
0,0659 0,0625
0,0336 0,03125
0,0165 0,015625
0,1294 0,125
0,0631 0,0625
0,0322 0,03125
Table 10 Table of errors of l2 and l∞ norms for Example 3 and clouds of points 1, 2 and 3. m0
l2
l∞
442 1629 8147
0,0000445 0,0000122 0,00000247
0,000112 0,000032 0,0000071
Table 11 Table of errors of l2 and l∞ norms for Example 3 and clouds of points 4, 5 and 6. m0
l2
l∞
119 502 1928
0,0000453 0,0000119 0,0000029
0,000124 0,0000342 0,0000083
Defining the error enj = unj − Ujn , it is obtained en0+1 − 2en0 + en0−1 (∆t)2
= c(un0 )(−m0 un0 +
s ∑ j=1
+ F (un0 ) − F (U0n ).
mj unj ) − c(U0n )(−m0 U0n +
s ∑
mj Ujn ) +
j=1
(46)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 8. Example 3 on clouds of points of Figs. 1–3.
Table 12 Table of errors of l2 and l∞ norms for Example 3 and clouds of points 7, 8 and 9. m0
l2
l∞
460 1757 7346
0,000148 0,000036 0,0000094
0,000379 0,0000981 0,0000233
Table 13 Table of errors (Eq. (56)) of l2 and l∞ norms for sine–Gordon equation. Time (s)
1
3
5
7
global error maximum error
5, 725 · 10−3 1, 596 · 10−2
9, 542 · 10−3 4, 131 · 10−2
1, 310 · 10−2 5, 154 · 10−2
1, 291 · 10−2 5, 867 · 10−2
13
14
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 9. Exact, approximate solution, absolute values of the local errors of Example 3 on cloud of points 9 of Fig. 3.
en0+1 − 2en0 + en0−1 (∆t)2
= [c(un0 )(−m0 un0 +
s ∑
mj unj ) − c(un0 )(−m0 U0n +
j=1
+ c(un0 )(−m0 U0n +
s ∑
mj Ujn )
j=1
s ∑
mj Ujn ) − c(U0n )(−m0 U0n +
j=1
s ∑
mj Ujn )] +
j=1
+ F (un0 ) − F (U0n )
(47)
and using the Mean Value Theorem
en0+1 − 2en0 + en0−1 = (∆t)2 [c(un0 )(−m0 en0 +
s ∑
mj enj ) +
j=1 s
∑ ∂c ∂F + en0 (χj )(−m0 U0n + mj Ujn ) + en0 (χj )] ∂U ∂U
(48)
j=1
Taking into account that enj = en0 eiϕj and rearranging Eq. (48), we obtain
e0n+1 = en0 [2 + (∆t)2 (c(un0 )(−m0 +
∂F + (χj ))] + en−1 . ∂U
s ∑ j=1
s
mj eiϕ ) +
∑ ∂c mj Ujn ) + (χj )(−m0 U0n + ∂U j=1
(49)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
15
Fig. 10. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 1.
If we substitute m0 =
∑s
j=1
mj (Eq. (20))
en0+1 = en0 [2 − (∆t)2 (c(un0 )
+
s ∑
∂F (χj ))] + en−1 . ∂U
s
mj (1 − eiϕ ) +
j=1
∑ ∂c (χj )(−m0 U0n + mj Ujn ) + ∂U j=1
(50)
Applying the von Neumann criterion for stability analysis, it holds
ξ 2 − [2 − 2R(∆t)2 ]ξ + 1 = 0
(51)
where N ∑
N
∑ ∂F ∂c (χj )(−m0 U0n + mj Ujn ) − (χj ). ∂U ∂U j=1 j=1 ∑s Operating as in previous subsections and considering j=1 mj sin(ϕj ) = 0 the following inequality is obtained 2R = c(un0 )
mj (1 − eiϕj ) −
−1 ≤ 1 − R∆t 2 ≤ 1
(52)
(53)
by relating the terms of the inequality Eq. (53), it is obtained 0 ≤ R∆t 2 ≤ 2.
(54)
The inequality of Eq. (54) may be written as (∆t)2 ≤
2
R Eq. (55) is the condition for stability of the quasilinear case. □ As before, the minimum value obtained among all the stars is taken as the value ∆t in the stability condition.
(55)
16
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 11. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 3.
Denominator of the expression Eq. (55) must be positive, which requires us to demand that R > 0. By the same ∂ 2U ∂ c ∂F ∂ 2U + , and are bounded as U , c , F are all discussion as in 3.1 and having in mind that the expressions 2 2 ∂x ∂y ∂U ∂U differentiable, one can choose hi so that R > 0. 4. Numerical results In this section we present the numerical results obtained by solving non-linear hyperbolic PDE’s with Dirichlet boundary conditions and the initial conditions, which correspond to non-linear problems in physics and engineering. In all the considered cases the global error has been computed according to l2 and l∞ norms
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨
NI ( )2 ∑ sol(i) − exac(i) √ i=1
l2 = ⎪ NI ⎪ ⎪ ⎪ ⎪ ⎩ l∞ = max|sol(i) − exac(i)|
(56)
where sol(i) is the GFD solution at node i, exac(i) is the exact value of the solution at node i and NI is the number of nodes of the domain Ω . Firstly three examples (1, 2 and 3) of irregular clouds of nodes are presented and for each one of them, three models are generated using uniform refinement, then nine irregular clouds of inner nodes are shown in Figs. 1–3. The optimum time step (∆t) as defined in Eqs. (43) and (55) depends, among others, on the value of parameter m0 defined in Eq. (20). This parameter depends on the number of nodes of the star (s), on the criterion to select these nodes 1 and on the weighting function used in Eq. (9). In this paper s = 8, the distance criterion and w (hi , ki ) = have been dist 4 used. The maximum value of m0 in each one of the nine clouds of points is as follows: cloud 1 (55 nodes, m0 = 442), cloud 2 (197 nodes, m0 = 1629), cloud 3 (743 nodes, m0 = 8147), cloud 4 (64 nodes, m0 = 119), cloud 5 (224 nodes, m0 = 502),
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
17
Fig. 12. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 5.
cloud 6 (832 nodes, m0 = 1928), cloud 7 (45 nodes, m0 = 460), cloud 8 (153 nodes, m0 = 1757) and cloud 9 (561 nodes, m0 = 7346). Clouds 1, 2 and 3 are shown in Fig. 1, clouds 4, 5 and 6 are shown in Fig. 2, and clouds 7, 8 and 9 are shown in Fig. 3. In all the cases uniform refinement has been used. Irregular clouds of nodes are generated and refined in such a way that the inter-particle distance is reduced approximately by a factor of two (see Figs. 1–3). The stability conditions given in Eqs. (43) and (55) are used to obtain the corresponding optimal time steps ∆t for each cloud of nodes. The total time interval used has been [0, 0.25]. This total interval is divided by a whole number to obtain the so called computing time step that is elected taking values very close but a little smaller than ∆t (see Tables 1, 5 and 9). 4.1. Example 1 Let us consider the equation
∂ 2U ∂ 2U ∂ 2U = + + U2 ∂t2 ∂ x2 ∂ y2
(57)
with the exact solution 1 x2 y2 (t + 2)2 −1 U(x, y, t) = − ( + − ) . 2 4 4 4
(58)
This example corresponds to the semilinear case, Eq. (21), whose consistency and stability have been studied in the previous section. Time increments ∆t have been computed by using the condition of stability. Recall the bound for the time increment given by Eq. (43):
∆t 2 ≤ min
4 2 |m0 | −
∂F ∂U
18
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 13. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 7.
which on the considered domain, since
−
4 ∂F = −2U = 2 ≤1 ∂U x + y2 + (t + 2)2
reduces to
∆t 2 ≤ min
4 2 |m0 | + 1
.
(59)
In Table 1 the optimum values of ∆t are given. In Fig. 4 the l2 and l∞ norms as a function of the values of parameter m0 are shown for the nine clouds of nodes in the case of example 1. The y-range of the plot has been shown in logarithmic scale. Tables of the convergence for l2 and l∞ errors are shown (see Tables 2–4). These Tables give precise information for numerical convergence. As it can be seen it appears an inverse relation between max. m0 and error l2 . In Fig. 5 the exact and approximated solution and the absolute values of the local errors of example 1 in cloud of points 3 are shown. 4.2. Example 2 Let us consider the equation
∂ 2U ∂ 2U ∂ 2U = (1 + U)( + ) ∂t2 ∂ x2 ∂ y2
(60)
with the exact solution: U(x, y, t) = t(1 + x3 − 3xy2 ).
(61)
This example corresponds to a quasilinear case, Eq. (22), whose consistency and stability have been studied in the previous section.
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
19
∆t has been calculated by using the condition of stability: In this example of quasilinear PDE the bound for the time increment is given by (∆t)2 ≤ min
2 R
where 2R = c(U0k )
s ∑
s
⃗j ) − mj (1 − eiα⃗ h
∑ ∂c k mj Ujk ). (U0 )(−m0 U0k + ∂U
(62)
j=1
j=1
This expression can be approximated by 2R ≃ 2m0 c(U) −
∂ ln(c(U)) ∂ 2 U . ∂U ∂t2
Applying this to the example and using that
∆t 2 ≤ min
2 m0 (1 + U)
≤
1 3m0
(63)
∂ 2U = 0 and c(U) = 1 + U it is obtained: ∂t2
.
(64)
In Table 2 the optimum values of ∆t are given. In Fig. 6 the l2 and l∞ norms as a function of the values of parameter m0 are shown for the nine clouds of nodes in the case of example 2. The y-range of the plot has been shown in logarithmic scale. Tables of the convergence for l2 and l∞ errors are shown (see Tables 6–8). These Tables give precise information for numerical convergence. As it can be seen it appears an inverse relation between max. m0 and error l2 . In Fig. 7 the exact and approximated solution and the absolute values of the local errors of example 2 in cloud of points 6 are shown. 4.3. Example 3 Let us consider the equation
∂ 2U ∂ 2U ∂ 2U = + 2 − 2U −1 2 2 ∂t ∂x ∂y
(65)
with the exact solution U(x, y, t) =
√
(x + 1)2 + (y + 1)2 − t 2 .
(66)
This example corresponds to a semilinear case, Eq. (21), whose consistency and stability have been studied in the previous section. ∆t has been calculated by using the condition of stability: As in Example 1, it is a semilinear PDE whose time bound is given by 4
∆t ≤
2 |m0 | −
∂F ∂U
.
In this case, having in mind that
−
∂F −2 −2 = 2 = ≤2 ∂U U (x + 1)2 + (y + 1)2 − t 2
it is obtained (∆t)2 ≤
4 2 |m0 | + 2
.
(67)
In Table 3 the optimum values of ∆t are given. In Fig. 8 the l2 and l∞ norms as a function of the values of parameter m0 are shown for the nine clouds of nodes in the case of example 3. The y-range of the plot has been shown in logarithmic scale. Tables of the convergence for l2 and l∞ errors are shown (see Tables 10–12). These Tables give precise information for numerical convergence. As it can be seen it appears an inverse relation between max. m0 and error l2 . In Fig. 9 the exact and approximate solution and the absolute values of the local errors of example 3 in cloud of points 9 are shown.
20
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
4.4. Sine–Gordon equation The two-dimensional sine–Gordon equation is given by
∂ 2U ∂U ∂ 2U ∂ 2U + β = + + φ (x, y) sin(U) ∂t2 ∂t ∂ x2 ∂ y2 in some continuous domain with appropriate boundary and initial conditions. The parameter β > 0 is a dissipative term and the function φ (x, y) is called the Josephson current density. The sine–Gordon (SG) equation is one of the most important dynamical models in nonlinear science. The sine–Gordon equation possesses soliton wave type solutions. A soliton is a certain kind of wave which travels with permanent speed and form. Solitons are used in communications, because they can transfer waves over large distances with no errors, and have various physical applications such as in electronics, optical fibres signals, superconductivity and so on. The sine–Gordon equation is a semilinear wave equation. Let us consider the sine–Gordon equation with β = 0 and φ (x, y) = −1
∂ 2U ∂ 2U ∂ 2U = + − sin(U) ∂t2 ∂ x2 ∂ y2
(68)
with the exact solution, see [34], U(x, y, t) = 4 tan−1 (ex+y−t ).
(69)
The computational domain [−7, 7] × [−7, 7] is discretized into 3249 nodes and ∆t = 0.25 s. The results obtained for t = 1, 3, 5, 7 s respectively are given in Table 13 and Figs. 10–13. 5. Conclusions In this paper an explicit method using GFDM for solving different non-linear hyperbolic PDEs, has been considered. Consistency has been proven, for semilinear and quasilinear equations, and conditions of stability have been obtained and implemented in different numerical examples. All these illustrate the viability of the application of GFDM for solving second order hyperbolic non-linear PDEs in 2D. The efficiency of the developed method is clearly shown. The numerical experiments show that the explicit-GFDM can be used for solving very different second order non-linear hyperbolic problems. Acknowledgements The authors acknowledge the support of the Escuela Técnica Superior de Ingenieros Industriales (UNED) of Spain, project 2019-IFC02 and of the Universidad Politécnica de Madrid (UPM) (Research groups 2019). 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, Bull. Amer. Math. Soc. 42 (4) (2012) 507–554. [3] S.-T.J. Yu, L. Yang, R.L. Lowe, S.E. Bechtel, Numerical simulation of linear and nonlinear waves in hypoelastic solids by the cese method, Wave Motion 47 (2010) 168–182. [4] S. Bonazzola, E. Gourgoulhon, J.A. Marck, Spectral methods in general relativistic astrophysics, J. Comput. Appl. Math. 109 (1999) 433–473. [5] F. Bloom, Systems of nonlinear hyperbolic equations associated with problems of classical electromagnetic theory, Comput. Math. Appl. 11 (1985) 261–279. [6] H. Qiu, Y. Zhang, Decay of the 3D quasilinear hyperbolic equations with nonlinear damping, Adv. Math. Phys. (2017) 2708483, 13. [7] P.S. Jensen, Finite difference technique for variable grids, Comput. Struct. 2 (1972) 17–29. [8] T. Liszka, J. Orkisz, The finite difference method at arbitrary irregular grids and its application in applied mechanics, Comput. Struct. 11 (1980) 83–95. [9] J. Orkisz, Finite difference method (Part, III), in: M. Kleiber (Ed.), Handbook of Computational Solid Mechanics, Spriger-Verlag, Berlin, 1998. [10] N. Perrone, R. Kao, A general finite difference method for arbitrary meshes, Comput. Struct. 5 (1975) 45–58. [11] J.J. Benito, F. Ureña, L. Gavete, Influence of several factors in the generalized finite difference method, Appl. Math. Model. 25 (2001) 1039–1053. [12] F. Ureña, J.J. Benito, L. Gavete, Application of the generalized finite difference method to solve the advection-diffusion equation, J. Comput. Appl. Math. 235 (2011) 1849–1855. [13] J.J. Benito, F. Ureña, L. Gavete, Solving parabolic and hyperbolic equations by the generalized finite difference method, J. Comput. Appl. Math. 209 (2007) 208–233. [14] L. Gavete, F. Ureña, J.J. Benito, A. Garcia, M. Ureña, E. Salete, Solving second order non-linear elliptic partial differential equations using generalized finite difference method, J. Comput. Appl. Math. 318 (2017) 378–387. [15] F. Ureña, L. Gavete, A. Garcia, J.J. Benito, A.M. Vargas, Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM), J. Comput. Appl. Math. http://dx.doi.org/10.1016/j.cam2018.02.016. [16] C.M. Fan, Y.K. Huang, P.W. Li, C.L. Chiu, Application of the generalized finite-difference method to inverse biharmonic boundary value problems, Numer. Heat Transfer B 65 (2) (2014) 129–154. [17] C.M. Fan, P.W. Li, W. Yeih, Generalized finite-difference method for solving two-dimensional inverse Cauchy problems, Inverse Probl. Sci. Eng. 23 (5) (2015) 737–759.
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
21
[18] S.M. Hosseini, Application of a hybrid mesh-free method based on generalized finite difference (GFD) method for natural frequency analysis of functionally graded nanocomposite cylinders reinforced by carbon nanotubes, Comput. Model. Eng. Sci. 95 (2013) 1–29. [19] S.M. Hosseini, Shock-induced two dimensional coupled non-Fickian diffusion-elasticity analysis using meshless generalized finite difference (GFD) method, Eng. Anal. Bound. Elem. 61 (2015) 232–240. [20] J. Izadian, N. Ranjbar, M. Jalili, The generalized finite difference method for solving elliptic equation on irregular mesh, World Appl. Sci. J. (ISSN: 1818-4952) 21 (2013) 95–100 (Special Issue of Applied Math.). [21] L. Gavete, J.J. Benito, F. Ureña, Generalized finite differences for solving 3D elliptic and parabolic equations, Appl. Math. Model. 40 (2016) 955–965. [22] Q. Hua, Y. Gu, W. Qu, W. Chen, C. Zang, A meshless generalized finite-difference method for inverse Cauchy problems associated with three-dimensional inhomogeneous Helmholtz-type equations, Eng. Anal. Bound. Elem. 82 (2017) 162–171. [23] H.F. Chan, C.M. Fan, C.W. Kuo, Generalized finite difference method for solving two-dimensional non-linear obstacle problems, Eng. Anal. Bound. Elem. 37 (2013) 1189–1196. [24] T. Zhang, Y.F. Ren, H.F. Chan, P.W. Li, Simulation of two-dimensional sloshing phenomenon by generalized finite difference method, Eng. Anal. Bound. Elem. 63 (2016) 82–91. [25] T. Zhang, Y.F. Ren, Z.Q. Yang, C.M. Fan, P.W. Li, Application of generalized finite-difference method to propagation of nonlinear water waves in numerical wave plume, Ocean Eng. 123 (2016) 278–290. [26] B. Mochnacki, E. Majchrzak, Numerical modelling of casting solidification using generalized finite difference method, Mater. Sci. Forum 638–642 (2010) 2676–2681. [27] P.W. Li, C.M. Fan, Generalized finite-difference method for two-dimensional shallow water equations, Eng. Anal. Bound. Elem. 80 (2017) 58–71. [28] L. Gavete, M.L. Gavete, F. Ureña, J.J. Benito, An approach to refinement of irregular clouds using generalized finite differences, Math. Probl. Eng. 2015 (2015) http://dx.doi.org/10.1155/2015/283757, 9 pages. [29] L. Gavete, F. Ureña, J.J. Benito, M. Ureña, M.L. Gavete, Solving elliptical equations in 3D by means of an adaptive refinement in generalized finite differences, Math. Probl. Eng. (2018) http://dx.doi.org/10.1155/2018/9678473, 14 pages. [30] X. Wang, R. Samulyak, X. Jiao, K. Yu, AP-cloud: Adaptive particle-in-cloud method for optimal solutions to Vlasov–Poisson equation, J. Comput. Phys. 316 (2016) 682–699. [31] R. Samulyak, X. Wang, H. Ch. Chen, Lagrangian particle method for compressible fluid dynamics, J. Comput. Phys. 362 (2018) 1–19. [32] P. Lancaster, K. Salkauskas, Curve and Surface Fitting, Ed. Academic Press, 1986. [33] D. Levin, The approximation power of moving least squares, Math. Comp. 67 (224) (1998) 1517–1531. [34] H.S. Shukla, Mohammad Tamsir, V.K. Srivastava, Srivastava Numerical simulation of two dimensional sine-Gordon solitons using modied cubic B-spline differential quadrature method, AIP Adv. 5 (2015) 017121, http://dx.doi.org/10.1063/1-4906256.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Solving second order non-linear hyperbolic PDEs using generalized finite difference method (GFDM) ∗
F. Ureña a , , L. Gavete b , A. García a , J.J. Benito a , A.M. Vargas c a b c
UNED, ETSII, Madrid, Spain UPM, ETSIM, Madrid, Spain UCM, Madrid, Spain
article
info
Article history: Received 17 May 2018 Received in revised form 1 March 2019 Keywords: Meshless methods Generalized finite difference method Non-linear hyperbolic partial differential equations
a b s t r a c t 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 hyperbolic PDEs. This paper illustrates that the GFD explicit formulae developed to obtain the different derivatives of the PDEs are based on the existence of a positive definite matrix which is obtained using moving least squares approximation and Taylor series development. Consistency and stability are shown in this paper for semilinear and quasilinear hyperbolic equations. This paper shows the application of the GFDM for solving second order non-linear hyperbolic problems. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Modern numerical methods, in particular those for solving non-linear PDEs, have been developed in recent years using finite differences, finite elements, finite volumes or spectral methods. A review of numerical methods for non-linear partial differential equations is given by Polyanin [1] and Tadmor [2]. Nonlinear hyperbolic partial differential equations have been applied in different fields such as in hypoelastic solids [3], astrophysics [4], electromagnetic theory [5], propagation of heat waves [6] and other disciplines. In this paper we use a meshless method called generalized finite difference method (GFDM) for solving different hyperbolic non-linear PDEs. Researchers as Jensen [7], Liszka and Orkisz [8], Orkisz [9] and Perrone and Kao [10] have contributed to the development of the GFDM in different aspects of its applications. Benito, Gavete and Ureña [11–13] have studied the influence of several factors and developed the explicit formulae and h-adaptive method for the solution of the PDEs in 2D. GFDM was also used to solve both non-linear elliptic [14] and parabolic equations [15]. Various applications of the GFDM have been developed for solving different inverse problems by Fan et al. [16,17]. Hosseini [18,19] has applied the GFDM for natural frequency analysis of nanocomposite cylinders and for shock-induced two dimensional coupled non-Fickian diffusion-elasticity analysis. ∗ Corresponding author. E-mail addresses: [email protected] (F. Ureña), [email protected] (L. Gavete), [email protected] (A. García), [email protected] (J.J. Benito), [email protected] (A.M. Vargas). https://doi.org/10.1016/j.cam.2019.05.028 0377-0427/© 2019 Elsevier B.V. All rights reserved.
2
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
The GFDM has also been used to solve 3D partial differential equations by Izadian et al. [20], Gavete et al. [21] and Hua et al. [22]. The GFDM was also applied for different practical tasks as in [23] for solving two dimensional non-linear obstacle problems, in [24,25] for simulating two dimensional sloshing phenomena and for the propagation of nonlinear water waves, in [26] for the modelling of casting solidification, and in [27] for solving two-dimensional shallow water equations. Adaptive methods for solving elliptical problems with singularities using GFD have been proposed [28,29]. An adaptive Particle-in-Cloud method has been used for solving the Vlasov–Poisson equation [30] and a Lagrangian particle method for solving Euler equations [31]. Both methods use GFD to improve the approximation of the differential operators involved. This paper shows that the GFDM can be used to solve second order hyperbolic nonlinear PDEs with different irregular clouds of points in 2D. Hyperbolic equations have been solved using an explicit method and their consistency and stability have been studied taking into account the irregularity of the cloud of points. The numerical results show the high accuracy obtained. A brief outline of this paper is as follows. In Section 2 the explicit method and GFDM to solve non-linear hyperbolic partial differential equations are shown. In Section 3, consistency and stability are studied. In Section 4 the proposed method is applied to different test problems to show the high accuracy of the method. In the last section, we present some conclusions. 2. Explicit method and GFDM: application to non-linear hyperbolic partial differential equations Consider the following non-linear problem in the domain D = [0, T ] × Ω with Ω ⊂ R2
∂ 2U = LΩ [ U ] ∂t2 where LΩ [U ] is a non-linear operator, Γ is the boundary of the domain Ω , with the boundary condition UΓ = f (t)
(1)
(2)
and initial conditions U(x, y, 0) = g(x, y)
(3)
∂U (x, y, 0) = h(x, y) ∂t
(4)
where f , g and h are three known functions. To solve the problem described by Eqs. (1)–(4), using the explicit method, the time derivative is approximated by un0+1 − 2un0 + un0−1 ∂ 2 u(x0 , y0 , n∆t) = + Θ ((∆t)2 ) ∂t2 (∆t)2
(5)
and spatial derivatives are approximated using GFD [13], as described in the following subsection. 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 Ω , 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 [11] 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: Ui = U0 + hi
( ) 2 ∂ U0 1 2 ∂ 2 U0 ∂ 2 U0 ∂ U0 2 ∂ U0 + ki + hi + k + 2h k + ···, i i i ∂x ∂y 2 ∂ x2 ∂ y2 ∂ x∂ y
i = 1, . . . , s
(6)
We define: cTi = {hi , ki , DT = {
h2i k2i
, , hi ki }, 2 2 ∂ u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0
∂ u0 , , , , }. ∂ x ∂ y ∂ x2 ∂ y2 ∂ x∂ y
(7) (8)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
3
If in Eq. (6) the higher than second order terms are ignored, a second order approximation for the function Ui is obtained. We denote it by ui . It is then possible to define the function B(u) as s ∑ ∂ u0 ∂ u0 [(u0 − ui ) + hi + ki + ∂x ∂y
B(u) =
i=1
∂ 2 u0 2 2 ∂ 2 u0 ∂ 2 u0 + k2i + 2hi ki )] wi 2 2 ∂x ∂y ∂ x∂ y
1
+ (h2i 2
(9)
where wi = w (hi , ki ) are positive symmetrical weighting functions, having the property, as defined in [32], of decreasing monotonically in magnitude as the distance to the centre of the weighting function increases. See also [33]. 1 Some weighting functions as potential or exponential exp(−n(dist 2 )) can be used [11], where n ∈ N. If the norm dist n given by Eq. (9) 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 )
(10)
where h1 k1
⎛ ⎜
h2 k2
.. .
··· ··· .. .
hs ks
h2 k2
···
hs ks
.. .
A=⎜ ⎝
h1 k1
⎞⎛
ω12
ω22
⎟⎜ ⎟⎝ ⎠
.. .
··· ω
2 s
⎞⎛ ⎟⎜ ⎠⎜ ⎝
h1 h2
.. .
hs
k1 k2
.. .
··· ··· .. .
h1 k1 h2 k2 ⎟
ks
···
hs ks
⎞
.. .
⎟ ⎠
(11)
and s ∑
(−u0 + ui )hi wi2 ,
bT = {
i=1
s ∑
(−u0 + ui )ki wi2 ,
i=1
s ∑
( − u0 + ui )
i=1
k2i 2
wi2 ,
s ∑
( − u0 + ui )
i=1
h2i 2
wi2 ,
s ∑
(−u0 + ui )hi ki wi2 }.
i=1
The matrix A is a positive definite matrix and the approximation is of second order Θ (h2i , k2i ) [14]. The assumption stated on the nodes distribution is related with the criterion selected to choose the stars nodes, because in each one of the stars, the matrix A must be positive definite to obtain the approximate derivatives, hence the existence of ill-conditioned stars of nodes is restricted. If we define A−1 = Q Q T
(12)
D = Q Q T b.
(13)
then
Eq. (13) can be rewritten as D = −u0 Q Q T
s ∑
wi2 c i + Q Q T
s ∑
i=1
ui wi2 c i
(14)
i=1
or D = Q Q T W (u − u0 1)
(15)
where h1 w12
h2 w22
⎜ k w2 ⎜ 1 1 ⎜ 2 ⎜ h1 2 ⎜ W = ⎜ 2 w1 ⎜ 2 ⎜ k ⎜ 1 w2 ⎝ 2 1
k2 w22
⎛
h1 k1 w12
h22 2 k22
w22
··· ··· .. .
w22
.. .
2 h2 k2 ws2
···
hs ws2 ks ws2 h2s 2 k2s
ws2 w2
2 s hs ks ws2
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
and 1 = {1, 1, . . . , 1} ;
u = {u1 , u2 , . . . , us }T .
4
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Lemma 1. The truncation error of the partial derivatives is zero when (hi , ki ) → (0, 0). Proof. Taking into account Eq. (15) D = Q Q T W (U − U0 1) =
⎛⎛ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ T ⎜ QQ W ⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎝⎝
⎛
∂ ∂ + k1 )U0 + ∂x ∂y ∂ ∂ U0 + (h2 + k2 )U0 + ∂x ∂y U0 + (h1
U0 + (hs
h1
⎜ h2 QQTW ⎜ ⎝ .. . hs
k1 k2
∂ ∂x
.. .
··· ··· .. .
ks
···
∂ ∂ + k1 )2 )U0 + · · · 2! ∂x ∂y 1 ∂ ∂ (h2 + k2 )2 )U0 + · · · 2! ∂x ∂y .. . ∂ 1 ∂ ∂ + ks )U0 + (hs + ks )2 )U0 + · · · ∂y 2! ∂x ∂y ⎛ ∂U ⎞ 0 ⎜ ∂x ⎟ ⎜ ∂U ⎟ 0 ⎟ ⎜ ⎟ ⎞⎜ h1 k 1 ⎜ ∂y ⎟ ⎟ ⎜ 2 h2 k 2 ⎟ ⎜ ∂ U 0 ⎟ T ⎟ ⎜ .. ⎟ ⎠ ⎜ ∂ x2 ⎟ + Q Q W R . ⎟ ⎜ 2 ⎜ ∂ U0 ⎟ hs k s ⎟ ⎜ ⎜ ∂ y2 ⎟ ⎟ ⎜ ⎝ ∂ 2U ⎠ 1
(h1
⎞ ⎟ ⎛ ⎟ ⎟ ⎟ ⎜ ⎟−⎜ ⎟ ⎝ ⎟ ⎟ ⎠
⎞ U0 ⎟ ⎟ U0 ⎟⎟
⎞⎟
⎟ .. ⎟ ⎠⎟ = . ⎟ ⎟ U0 ⎠
(16)
0
∂ x∂ y where
⎛ ⎜ ⎜ ⎜ ⎜ R=⎜ ⎜ ⎜ ⎜ ⎝
∂ ∂ + k1 )3 )U0 + · · · 3! ∂x ∂y ∂ ∂ 1 (h2 + k2 )3 )U0 + · · · 3! ∂x ∂y .. . 1 ∂ ∂ (hs + ks )3 )U0 + · · · 3! ∂x ∂y 1
(h1
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
thus Q Q T W (U − U0 1) = Q Q T AD + Q Q T WR = A−1 AD + Q Q T WR = D + Q Q T WR
(17)
Hence truncation error for spatial partial derivatives is TE = Q Q T WR
(18)
and therefore lim
(hi ,ki )→(0,0)
TE → 0. □
(19)
Remark 1. The explicit expressions of TE may be seen in [13]. Thus, spatial derivatives using GFD [11–13,21] are denoted by
⎧ s s ∑ ∑ ⎪ ∂ u(x0 , y0 , n△t) ⎪ 2 2 ⎪ = −λ u + λ u + Θ (h , k ) , w ith λ = λi ⎪ 0 0 i i 0 i i ⎪ ∂x ⎪ ⎪ i = 1 i = 1 ⎪ ⎪ s s ⎪ ∑ ∑ ⎪ ∂ u(x0 , y0 , n△t) ⎪ 2 2 ⎪ , k ) , w ith µ = µi = −µ u + µ u + Θ (h ⎪ 0 0 i i 0 i i ⎨ ∂y i=1 i=1 s ∑ ⎪ ∂ 2 u(x0 , y0 , n△t) ∂ 2 u(x0 , y0 , n△t) ⎪ ⎪ + = − m u + mi ui + Θ (h2i , k2i ) ⎪ 0 0 ⎪ ∂ x2 ∂ y2 ⎪ ⎪ i=1 ⎪ ⎪ s ⎪ ∑ ⎪ ⎪ ⎪ mi ⎪ ⎩ with m0 = i=1
(20)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
5
3. Consistency and stability of the scheme for non-linear hyperbolic PDEs In this section consistency and stability of non-linear hyperbolic PDEs, using GFDM, are studied. We will do so by starting from the following definitions:
• A partial differential equation is semilinear if the coefficients of its highest derivatives are functions of the spaces 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 case LΩ [U ] =
∂ 2U ∂ 2U + + F (U). ∂ x2 ∂ y2
(21)
Quasilinear case LΩ [U ] = c(U)(
∂ 2U ∂ 2U + 2 ) + F (U). 2 ∂x ∂y
(22)
Theorem 1. Let us consider the quasilinear equation (Eq. (22)) where c(U) and F (U) are two differentiable functions. The non-linear explicit scheme u0n+1 − 2un0 + un0−1
∆t 2
= c(un0 )(−m0 un0 +
s ∑
mj unj ) + F (un0 )
(23)
j=1
is consistent. Proof. The truncation error for time partial derivatives is TE = −
∆t 2 ∂ 4 U + Θ1 ((∆t)4 ) 12 ∂ t 4
(24)
and for spatial partial derivatives, taking into account Eqs. (18) and (19), is TE = c(U )Q Q T WR
(25)
where
(
c(U ) =
0
0
c(U)
c(U)
0
)
.
The truncation error of the scheme is the sum of Eqs. (24) and (25), then TE = −
∆t 2 ∂ 4 U + Θ1 ((∆t)4 ) + c(U )Q Q T WR . 12 ∂ t 4
(26)
Taking into account Lemma 1, we obtain lim
(hi ,ki ,∆t)→(0,0,0)
TE → 0.
(27)
Thus, the non-linear explicit scheme for quasilinear and semilinear equations is consistent.
□
Notice that if in Eq. (22) we make c(U) = 1, the proof is valid for the semilinear case. 3.1. Semilinear hyperbolic PDEs Theorem 2. Let us consider the semilinear equation (Eq. (21)) where F (U) is a differentiable function. The non-linear explicit scheme un0+1 − 2un0 + un0−1
∆t 2 is conditionally stable.
= −m0 un0 +
s ∑ j=1
mj unj + F (un0 )
(28)
6
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Proof. Consider the non-linear scheme un0+1 − 2un0 + un0−1 (∆t)2
= −m0 un0 +
N ∑
mj unj + F (un0 )
(29)
j=1
and let the same expression for the exact solution be U0n+1 − 2U0n + U0n−1 (∆t)2
= −m0 U0n +
N ∑
mi Ujn + F (U0n ).
(30)
j=1
If we define the error as enj = unj − Ujn and en0 = un0 − U0n , using the Mean Value Theorem for F , it is obtained
∂F (χ j ) ∂U ξ ∈ [0, 1]. Subtracting Eqs. (29) and (30) together with Eq. (31) it is obtained
F (un0 ) − F (U0n ) = en0 where χj = U0n + ξ en0 ,
e0n+1 − 2en0 + e0n−1 = (∆t)2 (−m0 en0 +
N ∑
mj enj + en0
j=1
∂F (χj )). ∂U
(31)
(32)
Let us consider the von Neumann criterion for stability analysis: en0 = ξ n eiα⃗ ·r⃗0 ;
enj = ξ n eiα⃗ ·⃗rj ; r⃗j = r⃗0 + δ⃗rj
(33)
where: ξ is called the amplification factor, α ⃗ is the vector of the wave numbers, r⃗0 is the vector of coordinates of the central node of the star and r⃗j is the vector of coordinates of the other nodes of the star. Considering Eq. (20) and including Eq. (33) in Eq. (32), it is obtained
ξ n+1 eiα⃗ ·r⃗0 − 2ξ n eiα⃗ ·r⃗0 + ξ n−1 eiα⃗ ·r⃗0 = = (∆t)2 [−m0 ξ n eiα⃗ ·r⃗0 + ξ n eiα⃗ ·r⃗0
N ∑
mj eiα⃗ ·δ rj + ξ n eiα⃗ ·r⃗0 ⃗
j=1
∂F (χj )]. ∂U
(34)
Eq. (34) is simplified by cancelling the term ξ n−1 eiα⃗ ·r⃗0 , which leads to
ξ 2 − 2ξ + 1 = (∆t)2 [−
N ∑
mj (1 − eiα⃗ ·δ rj ) + ⃗
j=1
∂F (χj )]ξ ∂U
(35)
or rearranging
ξ 2 − 2bξ + 1 = 0
(36)
where b=1−
N (∆t)2 ∑
2
[
mj (1 − eiα⃗ ·δ rj ) − ⃗
j=1
∂F (χj )] ∂U
(37)
Eq. (36) has two complex solutions: ξ1 and ξ2 with
{
√ ξ1 = b + √b2 − 1 ξ2 = b − b2 − 1
(38)
Stability condition implies
|ξ1 | = |ξ2 | = 1 ⇔ |b| ≤ 1
(39)
Denoting ϕj = α ⃗ · δ⃗rj , it is possible to write:
⎧ ⎨ ξ1 = eiϕj ξ = e−iϕj ⎩ ξ2 + ξ = 2b = 2 cos ϕ 1 2 j ∑N and considering j=1 mj sin(ϕj ) = 0, then the stability condition is −1 ≤ 1 −
N ∆t 2 ∑
2
(
j=1
mj (1 − cos(ϕj )) −
∂F (χj )) ≤ 1. ∂U
(40)
(41)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
7
Fig. 1. Clouds 1, 2, 3.
Table 1 Table of the computing time steps ∆τ and optimal time steps ∆t for Example 1 and clouds of points (1, . . . , 9). Mesh
1
2
3
4
5
6
7
8
9
∆t ∆τ
0,0672 0,0625
0,0341 0,03125
0,0162 0,015625
0,0659 0,0625
0,0336 0,03125
0,0165 0,015625
0,1294 0,125
0,0631 0,0625
0,0322 0,03125
If we relate the terms of the inequality Eq. (41), it is obtained N ∆t 2 ∑
2
(
mj (1 − cos(ϕj )) −
j=1
∂F (χj )) ≤ 2. ∂U
(42)
The inequality of Eq. (42), taking into account Eq. (20), may be written as (∆t)2 ≤
4 2 |m0 | −
∂F (χj ) ∂U
.
(43)
Eq. (43) is the condition for stability of semilinear case. □ The minimum value obtained among all the stars is taken as the value ∆t in the stability condition. The denominator of the expression (43) must be positive, which requires us to demand that, since m0 is positive, either ∂F 1 ∂F F is a decreasing function or 2m0 > . However, considering that m0 (m0 ∼ 2 ) and that is bounded, as F is a ∂U ∂U hi ∂F differentiable function, one could always choose hi so that 2m0 > . ∂U
8
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 2. Clouds 4, 5, 6.
Fig. 3. Clouds 7, 8, 9.
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 4. Example 1 on clouds of points of Figs. 1–3.
Table 2 Table of errors of l2 and l∞ norms for Example 1 and clouds of points 1, 2 and 3. m0
l2
l∞
442 1629 8147
0,0003612 0,0001006 0,0000201
0,0006701 0,0001815 0,0000379
Table 3 Table of errors of l2 and l∞ norms for Example 1 and clouds of points 4, 5 and 6. m0
l2
l∞
119 502 1928
0,000689 0,000165 0,0000432
0,000892 0,000263 0,0000793
9
10
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 5. Exact, approximate solution, absolute values of the local errors of Example 1 on cloud of points 3 of Fig. 1. Table 4 Table of errors of l2 and l∞ norms for Example 1 and clouds of points 7, 8 and 9. m0
l2
l∞
460 1757 7346
0,000284 0,0000782 0,0000193
0,000502 0,000181 0,0000338
Table 5 Table of the computing time steps ∆τ and optimal time steps ∆t for Example 2 and clouds of points (1, . . . , 9). Mesh
1
2
3
4
5
6
7
8
9
∆t ∆τ
0,0275 0,025
0,0143 0,0125
0,0064 0,00625
0,0269 0,025
0,0138 0,0125
0,0067 0,00625
0,0529 0,050
0,0258 0,025
0,0131 0,0125
Table 6 Table of errors of l2 and l∞ norms for Example 2 and clouds of points 1, 2 and 3. m0
l2
l∞
442 1629 8147
0,000139 0,0000391 0,00000792
0,000639 0,000183 0,000038
3.2. Quasilinear hyperbolic PDEs Theorem 3. Let us consider the quasilinear equation (Eq. (22)) where c(U) and F (U) are two differentiable functions. The non-linear explicit scheme un0+1 − 2un0 + un0−1
∆t 2 is conditionally stable.
= c(un0 )(−m0 un0 +
s ∑ j=1
mj unj ) + F (un0 )
(44)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
11
Fig. 6. Example 2 on clouds of points of Figs. 1–3. Table 7 Table of errors of l2 and l∞ norms for Example 2 and clouds of points 4, 5 and 6. m0
l2
l∞
119 502 1928
0,001082 0,000263 0,0000691
0,001481 0,0003488 0,0000902
Table 8 Table of errors of l2 and l∞ norms for Example 2 and clouds of points 7, 8 and 9. m0
l2
l∞
460 1757 7346
0,000232 0,0000613 0,0000149
0,000356 0,0000923 0,0000235
Proof. We proceed as in the proof of Theorem 2. Let us consider the scheme for the exact solution U0n+1 − 2U0n + U0n−1
∆t 2
= c(U0n )(−m0 U0n +
s ∑ j=1
mj Ujn ) + F (U0n ).
(45)
12
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 7. Exact, approximate solution, absolute values of the local errors of Example 2 on cloud of points 6 of Fig. 2.
Table 9 Table of the computing time steps ∆τ and optimal time steps ∆t for Example 3 and clouds of points (1, . . . , 9). Mesh
1
2
3
4
5
6
7
8
9
∆t ∆τ
0,0672 0,0625
0,0341 0,03125
0,0162 0,015625
0,0659 0,0625
0,0336 0,03125
0,0165 0,015625
0,1294 0,125
0,0631 0,0625
0,0322 0,03125
Table 10 Table of errors of l2 and l∞ norms for Example 3 and clouds of points 1, 2 and 3. m0
l2
l∞
442 1629 8147
0,0000445 0,0000122 0,00000247
0,000112 0,000032 0,0000071
Table 11 Table of errors of l2 and l∞ norms for Example 3 and clouds of points 4, 5 and 6. m0
l2
l∞
119 502 1928
0,0000453 0,0000119 0,0000029
0,000124 0,0000342 0,0000083
Defining the error enj = unj − Ujn , it is obtained en0+1 − 2en0 + en0−1 (∆t)2
= c(un0 )(−m0 un0 +
s ∑ j=1
+ F (un0 ) − F (U0n ).
mj unj ) − c(U0n )(−m0 U0n +
s ∑
mj Ujn ) +
j=1
(46)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 8. Example 3 on clouds of points of Figs. 1–3.
Table 12 Table of errors of l2 and l∞ norms for Example 3 and clouds of points 7, 8 and 9. m0
l2
l∞
460 1757 7346
0,000148 0,000036 0,0000094
0,000379 0,0000981 0,0000233
Table 13 Table of errors (Eq. (56)) of l2 and l∞ norms for sine–Gordon equation. Time (s)
1
3
5
7
global error maximum error
5, 725 · 10−3 1, 596 · 10−2
9, 542 · 10−3 4, 131 · 10−2
1, 310 · 10−2 5, 154 · 10−2
1, 291 · 10−2 5, 867 · 10−2
13
14
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 9. Exact, approximate solution, absolute values of the local errors of Example 3 on cloud of points 9 of Fig. 3.
en0+1 − 2en0 + en0−1 (∆t)2
= [c(un0 )(−m0 un0 +
s ∑
mj unj ) − c(un0 )(−m0 U0n +
j=1
+ c(un0 )(−m0 U0n +
s ∑
mj Ujn )
j=1
s ∑
mj Ujn ) − c(U0n )(−m0 U0n +
j=1
s ∑
mj Ujn )] +
j=1
+ F (un0 ) − F (U0n )
(47)
and using the Mean Value Theorem
en0+1 − 2en0 + en0−1 = (∆t)2 [c(un0 )(−m0 en0 +
s ∑
mj enj ) +
j=1 s
∑ ∂c ∂F + en0 (χj )(−m0 U0n + mj Ujn ) + en0 (χj )] ∂U ∂U
(48)
j=1
Taking into account that enj = en0 eiϕj and rearranging Eq. (48), we obtain
e0n+1 = en0 [2 + (∆t)2 (c(un0 )(−m0 +
∂F + (χj ))] + en−1 . ∂U
s ∑ j=1
s
mj eiϕ ) +
∑ ∂c mj Ujn ) + (χj )(−m0 U0n + ∂U j=1
(49)
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
15
Fig. 10. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 1.
If we substitute m0 =
∑s
j=1
mj (Eq. (20))
en0+1 = en0 [2 − (∆t)2 (c(un0 )
+
s ∑
∂F (χj ))] + en−1 . ∂U
s
mj (1 − eiϕ ) +
j=1
∑ ∂c (χj )(−m0 U0n + mj Ujn ) + ∂U j=1
(50)
Applying the von Neumann criterion for stability analysis, it holds
ξ 2 − [2 − 2R(∆t)2 ]ξ + 1 = 0
(51)
where N ∑
N
∑ ∂F ∂c (χj )(−m0 U0n + mj Ujn ) − (χj ). ∂U ∂U j=1 j=1 ∑s Operating as in previous subsections and considering j=1 mj sin(ϕj ) = 0 the following inequality is obtained 2R = c(un0 )
mj (1 − eiϕj ) −
−1 ≤ 1 − R∆t 2 ≤ 1
(52)
(53)
by relating the terms of the inequality Eq. (53), it is obtained 0 ≤ R∆t 2 ≤ 2.
(54)
The inequality of Eq. (54) may be written as (∆t)2 ≤
2
R Eq. (55) is the condition for stability of the quasilinear case. □ As before, the minimum value obtained among all the stars is taken as the value ∆t in the stability condition.
(55)
16
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 11. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 3.
Denominator of the expression Eq. (55) must be positive, which requires us to demand that R > 0. By the same ∂ 2U ∂ c ∂F ∂ 2U + , and are bounded as U , c , F are all discussion as in 3.1 and having in mind that the expressions 2 2 ∂x ∂y ∂U ∂U differentiable, one can choose hi so that R > 0. 4. Numerical results In this section we present the numerical results obtained by solving non-linear hyperbolic PDE’s with Dirichlet boundary conditions and the initial conditions, which correspond to non-linear problems in physics and engineering. In all the considered cases the global error has been computed according to l2 and l∞ norms
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨
NI ( )2 ∑ sol(i) − exac(i) √ i=1
l2 = ⎪ NI ⎪ ⎪ ⎪ ⎪ ⎩ l∞ = max|sol(i) − exac(i)|
(56)
where sol(i) is the GFD solution at node i, exac(i) is the exact value of the solution at node i and NI is the number of nodes of the domain Ω . Firstly three examples (1, 2 and 3) of irregular clouds of nodes are presented and for each one of them, three models are generated using uniform refinement, then nine irregular clouds of inner nodes are shown in Figs. 1–3. The optimum time step (∆t) as defined in Eqs. (43) and (55) depends, among others, on the value of parameter m0 defined in Eq. (20). This parameter depends on the number of nodes of the star (s), on the criterion to select these nodes 1 and on the weighting function used in Eq. (9). In this paper s = 8, the distance criterion and w (hi , ki ) = have been dist 4 used. The maximum value of m0 in each one of the nine clouds of points is as follows: cloud 1 (55 nodes, m0 = 442), cloud 2 (197 nodes, m0 = 1629), cloud 3 (743 nodes, m0 = 8147), cloud 4 (64 nodes, m0 = 119), cloud 5 (224 nodes, m0 = 502),
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
17
Fig. 12. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 5.
cloud 6 (832 nodes, m0 = 1928), cloud 7 (45 nodes, m0 = 460), cloud 8 (153 nodes, m0 = 1757) and cloud 9 (561 nodes, m0 = 7346). Clouds 1, 2 and 3 are shown in Fig. 1, clouds 4, 5 and 6 are shown in Fig. 2, and clouds 7, 8 and 9 are shown in Fig. 3. In all the cases uniform refinement has been used. Irregular clouds of nodes are generated and refined in such a way that the inter-particle distance is reduced approximately by a factor of two (see Figs. 1–3). The stability conditions given in Eqs. (43) and (55) are used to obtain the corresponding optimal time steps ∆t for each cloud of nodes. The total time interval used has been [0, 0.25]. This total interval is divided by a whole number to obtain the so called computing time step that is elected taking values very close but a little smaller than ∆t (see Tables 1, 5 and 9). 4.1. Example 1 Let us consider the equation
∂ 2U ∂ 2U ∂ 2U = + + U2 ∂t2 ∂ x2 ∂ y2
(57)
with the exact solution 1 x2 y2 (t + 2)2 −1 U(x, y, t) = − ( + − ) . 2 4 4 4
(58)
This example corresponds to the semilinear case, Eq. (21), whose consistency and stability have been studied in the previous section. Time increments ∆t have been computed by using the condition of stability. Recall the bound for the time increment given by Eq. (43):
∆t 2 ≤ min
4 2 |m0 | −
∂F ∂U
18
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
Fig. 13. Exact, approximate solution and absolute values of the local errors of the sine–Gordon equation for t = 7.
which on the considered domain, since
−
4 ∂F = −2U = 2 ≤1 ∂U x + y2 + (t + 2)2
reduces to
∆t 2 ≤ min
4 2 |m0 | + 1
.
(59)
In Table 1 the optimum values of ∆t are given. In Fig. 4 the l2 and l∞ norms as a function of the values of parameter m0 are shown for the nine clouds of nodes in the case of example 1. The y-range of the plot has been shown in logarithmic scale. Tables of the convergence for l2 and l∞ errors are shown (see Tables 2–4). These Tables give precise information for numerical convergence. As it can be seen it appears an inverse relation between max. m0 and error l2 . In Fig. 5 the exact and approximated solution and the absolute values of the local errors of example 1 in cloud of points 3 are shown. 4.2. Example 2 Let us consider the equation
∂ 2U ∂ 2U ∂ 2U = (1 + U)( + ) ∂t2 ∂ x2 ∂ y2
(60)
with the exact solution: U(x, y, t) = t(1 + x3 − 3xy2 ).
(61)
This example corresponds to a quasilinear case, Eq. (22), whose consistency and stability have been studied in the previous section.
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
19
∆t has been calculated by using the condition of stability: In this example of quasilinear PDE the bound for the time increment is given by (∆t)2 ≤ min
2 R
where 2R = c(U0k )
s ∑
s
⃗j ) − mj (1 − eiα⃗ h
∑ ∂c k mj Ujk ). (U0 )(−m0 U0k + ∂U
(62)
j=1
j=1
This expression can be approximated by 2R ≃ 2m0 c(U) −
∂ ln(c(U)) ∂ 2 U . ∂U ∂t2
Applying this to the example and using that
∆t 2 ≤ min
2 m0 (1 + U)
≤
1 3m0
(63)
∂ 2U = 0 and c(U) = 1 + U it is obtained: ∂t2
.
(64)
In Table 2 the optimum values of ∆t are given. In Fig. 6 the l2 and l∞ norms as a function of the values of parameter m0 are shown for the nine clouds of nodes in the case of example 2. The y-range of the plot has been shown in logarithmic scale. Tables of the convergence for l2 and l∞ errors are shown (see Tables 6–8). These Tables give precise information for numerical convergence. As it can be seen it appears an inverse relation between max. m0 and error l2 . In Fig. 7 the exact and approximated solution and the absolute values of the local errors of example 2 in cloud of points 6 are shown. 4.3. Example 3 Let us consider the equation
∂ 2U ∂ 2U ∂ 2U = + 2 − 2U −1 2 2 ∂t ∂x ∂y
(65)
with the exact solution U(x, y, t) =
√
(x + 1)2 + (y + 1)2 − t 2 .
(66)
This example corresponds to a semilinear case, Eq. (21), whose consistency and stability have been studied in the previous section. ∆t has been calculated by using the condition of stability: As in Example 1, it is a semilinear PDE whose time bound is given by 4
∆t ≤
2 |m0 | −
∂F ∂U
.
In this case, having in mind that
−
∂F −2 −2 = 2 = ≤2 ∂U U (x + 1)2 + (y + 1)2 − t 2
it is obtained (∆t)2 ≤
4 2 |m0 | + 2
.
(67)
In Table 3 the optimum values of ∆t are given. In Fig. 8 the l2 and l∞ norms as a function of the values of parameter m0 are shown for the nine clouds of nodes in the case of example 3. The y-range of the plot has been shown in logarithmic scale. Tables of the convergence for l2 and l∞ errors are shown (see Tables 10–12). These Tables give precise information for numerical convergence. As it can be seen it appears an inverse relation between max. m0 and error l2 . In Fig. 9 the exact and approximate solution and the absolute values of the local errors of example 3 in cloud of points 9 are shown.
20
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
4.4. Sine–Gordon equation The two-dimensional sine–Gordon equation is given by
∂ 2U ∂U ∂ 2U ∂ 2U + β = + + φ (x, y) sin(U) ∂t2 ∂t ∂ x2 ∂ y2 in some continuous domain with appropriate boundary and initial conditions. The parameter β > 0 is a dissipative term and the function φ (x, y) is called the Josephson current density. The sine–Gordon (SG) equation is one of the most important dynamical models in nonlinear science. The sine–Gordon equation possesses soliton wave type solutions. A soliton is a certain kind of wave which travels with permanent speed and form. Solitons are used in communications, because they can transfer waves over large distances with no errors, and have various physical applications such as in electronics, optical fibres signals, superconductivity and so on. The sine–Gordon equation is a semilinear wave equation. Let us consider the sine–Gordon equation with β = 0 and φ (x, y) = −1
∂ 2U ∂ 2U ∂ 2U = + − sin(U) ∂t2 ∂ x2 ∂ y2
(68)
with the exact solution, see [34], U(x, y, t) = 4 tan−1 (ex+y−t ).
(69)
The computational domain [−7, 7] × [−7, 7] is discretized into 3249 nodes and ∆t = 0.25 s. The results obtained for t = 1, 3, 5, 7 s respectively are given in Table 13 and Figs. 10–13. 5. Conclusions In this paper an explicit method using GFDM for solving different non-linear hyperbolic PDEs, has been considered. Consistency has been proven, for semilinear and quasilinear equations, and conditions of stability have been obtained and implemented in different numerical examples. All these illustrate the viability of the application of GFDM for solving second order hyperbolic non-linear PDEs in 2D. The efficiency of the developed method is clearly shown. The numerical experiments show that the explicit-GFDM can be used for solving very different second order non-linear hyperbolic problems. Acknowledgements The authors acknowledge the support of the Escuela Técnica Superior de Ingenieros Industriales (UNED) of Spain, project 2019-IFC02 and of the Universidad Politécnica de Madrid (UPM) (Research groups 2019). 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, Bull. Amer. Math. Soc. 42 (4) (2012) 507–554. [3] S.-T.J. Yu, L. Yang, R.L. Lowe, S.E. Bechtel, Numerical simulation of linear and nonlinear waves in hypoelastic solids by the cese method, Wave Motion 47 (2010) 168–182. [4] S. Bonazzola, E. Gourgoulhon, J.A. Marck, Spectral methods in general relativistic astrophysics, J. Comput. Appl. Math. 109 (1999) 433–473. [5] F. Bloom, Systems of nonlinear hyperbolic equations associated with problems of classical electromagnetic theory, Comput. Math. Appl. 11 (1985) 261–279. [6] H. Qiu, Y. Zhang, Decay of the 3D quasilinear hyperbolic equations with nonlinear damping, Adv. Math. Phys. (2017) 2708483, 13. [7] P.S. Jensen, Finite difference technique for variable grids, Comput. Struct. 2 (1972) 17–29. [8] T. Liszka, J. Orkisz, The finite difference method at arbitrary irregular grids and its application in applied mechanics, Comput. Struct. 11 (1980) 83–95. [9] J. Orkisz, Finite difference method (Part, III), in: M. Kleiber (Ed.), Handbook of Computational Solid Mechanics, Spriger-Verlag, Berlin, 1998. [10] N. Perrone, R. Kao, A general finite difference method for arbitrary meshes, Comput. Struct. 5 (1975) 45–58. [11] J.J. Benito, F. Ureña, L. Gavete, Influence of several factors in the generalized finite difference method, Appl. Math. Model. 25 (2001) 1039–1053. [12] F. Ureña, J.J. Benito, L. Gavete, Application of the generalized finite difference method to solve the advection-diffusion equation, J. Comput. Appl. Math. 235 (2011) 1849–1855. [13] J.J. Benito, F. Ureña, L. Gavete, Solving parabolic and hyperbolic equations by the generalized finite difference method, J. Comput. Appl. Math. 209 (2007) 208–233. [14] L. Gavete, F. Ureña, J.J. Benito, A. Garcia, M. Ureña, E. Salete, Solving second order non-linear elliptic partial differential equations using generalized finite difference method, J. Comput. Appl. Math. 318 (2017) 378–387. [15] F. Ureña, L. Gavete, A. Garcia, J.J. Benito, A.M. Vargas, Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM), J. Comput. Appl. Math. http://dx.doi.org/10.1016/j.cam2018.02.016. [16] C.M. Fan, Y.K. Huang, P.W. Li, C.L. Chiu, Application of the generalized finite-difference method to inverse biharmonic boundary value problems, Numer. Heat Transfer B 65 (2) (2014) 129–154. [17] C.M. Fan, P.W. Li, W. Yeih, Generalized finite-difference method for solving two-dimensional inverse Cauchy problems, Inverse Probl. Sci. Eng. 23 (5) (2015) 737–759.
F. Ureña, L. Gavete, A. García et al. / Journal of Computational and Applied Mathematics 363 (2020) 1–21
21
[18] S.M. Hosseini, Application of a hybrid mesh-free method based on generalized finite difference (GFD) method for natural frequency analysis of functionally graded nanocomposite cylinders reinforced by carbon nanotubes, Comput. Model. Eng. Sci. 95 (2013) 1–29. [19] S.M. Hosseini, Shock-induced two dimensional coupled non-Fickian diffusion-elasticity analysis using meshless generalized finite difference (GFD) method, Eng. Anal. Bound. Elem. 61 (2015) 232–240. [20] J. Izadian, N. Ranjbar, M. Jalili, The generalized finite difference method for solving elliptic equation on irregular mesh, World Appl. Sci. J. (ISSN: 1818-4952) 21 (2013) 95–100 (Special Issue of Applied Math.). [21] L. Gavete, J.J. Benito, F. Ureña, Generalized finite differences for solving 3D elliptic and parabolic equations, Appl. Math. Model. 40 (2016) 955–965. [22] Q. Hua, Y. Gu, W. Qu, W. Chen, C. Zang, A meshless generalized finite-difference method for inverse Cauchy problems associated with three-dimensional inhomogeneous Helmholtz-type equations, Eng. Anal. Bound. Elem. 82 (2017) 162–171. [23] H.F. Chan, C.M. Fan, C.W. Kuo, Generalized finite difference method for solving two-dimensional non-linear obstacle problems, Eng. Anal. Bound. Elem. 37 (2013) 1189–1196. [24] T. Zhang, Y.F. Ren, H.F. Chan, P.W. Li, Simulation of two-dimensional sloshing phenomenon by generalized finite difference method, Eng. Anal. Bound. Elem. 63 (2016) 82–91. [25] T. Zhang, Y.F. Ren, Z.Q. Yang, C.M. Fan, P.W. Li, Application of generalized finite-difference method to propagation of nonlinear water waves in numerical wave plume, Ocean Eng. 123 (2016) 278–290. [26] B. Mochnacki, E. Majchrzak, Numerical modelling of casting solidification using generalized finite difference method, Mater. Sci. Forum 638–642 (2010) 2676–2681. [27] P.W. Li, C.M. Fan, Generalized finite-difference method for two-dimensional shallow water equations, Eng. Anal. Bound. Elem. 80 (2017) 58–71. [28] L. Gavete, M.L. Gavete, F. Ureña, J.J. Benito, An approach to refinement of irregular clouds using generalized finite differences, Math. Probl. Eng. 2015 (2015) http://dx.doi.org/10.1155/2015/283757, 9 pages. [29] L. Gavete, F. Ureña, J.J. Benito, M. Ureña, M.L. Gavete, Solving elliptical equations in 3D by means of an adaptive refinement in generalized finite differences, Math. Probl. Eng. (2018) http://dx.doi.org/10.1155/2018/9678473, 14 pages. [30] X. Wang, R. Samulyak, X. Jiao, K. Yu, AP-cloud: Adaptive particle-in-cloud method for optimal solutions to Vlasov–Poisson equation, J. Comput. Phys. 316 (2016) 682–699. [31] R. Samulyak, X. Wang, H. Ch. Chen, Lagrangian particle method for compressible fluid dynamics, J. Comput. Phys. 362 (2018) 1–19. [32] P. Lancaster, K. Salkauskas, Curve and Surface Fitting, Ed. Academic Press, 1986. [33] D. Levin, The approximation power of moving least squares, Math. Comp. 67 (224) (1998) 1517–1531. [34] H.S. Shukla, Mohammad Tamsir, V.K. Srivastava, Srivastava Numerical simulation of two dimensional sine-Gordon solitons using modied cubic B-spline differential quadrature method, AIP Adv. 5 (2015) 017121, http://dx.doi.org/10.1063/1-4906256.