- Email: [email protected]
A local mesh-less collocation method for solving a class of time-dependent fractional integral equations: 2D fractional evolution equation
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A local mesh-less collocation method for solving a class of time-dependent fractional integral equations: 2D fractional evolution equation Mehdi Radmanesh a, M.J. Ebadi b,∗ a b
Department of Mathematics, Graduate University of Advanced Technology of Kerman, Kerman, Iran Department of Mathematics, Chabahar Maritime University, Chabahar, Iran
a r t i c l e
i n f o
Keywords: Radial basis function Fractional evolution equation Local mesh-less collocation method Global RBF Finite difference
a b s t r a c t The local radial basis function (LRBF) method is an excellent tool for solving variable-order time-fractional evolution equations. This method reduces the expensive computational cost of the conventional global RBF collocation (GRBFC) approach. The present paper represents an efficient mesh-less LRBF collocation approach for solving the two-dimensional (2D) fractional evolution equation for the arbitrary fractional order in complex-shaped domains. Having been used an appropriate finite difference (FD) technique to discrete the time variable, the LRBF method is applied in order to solve the equation. Unlike the conventional GRBFC method, consideration of the local nodes in a subdomain surrounding the given local point is needed for the proposed LRBF method. This method can reduce the ill-conditioning of the global resultant RBF interpolation matrix while it is also well efficient for large scale problems. Some illustrative examples are given to demonstrate the efficiency and proficiency of the method.
1. Introduction In recent years, many studies have been carried out to solve fractional differential equations, which are valuable tools for modeling many phenomena in economics, physics, and engineering [1]. Some of their significant applications were introduced in the research works of Podlubny [2], Oldham and Spanier [3], Bagley and Trovik [4], and Metzler and Klafter [5]. Also, the authors in [6] have given several remarkable studies on the theory, analysis, and recent historical development of the fractional calculus. In this manuscript, the following 2D fractional evolution equation is considered [7–11]: 𝜕𝑤(𝑥, 𝑦, 𝑡) − 𝐼 𝛼 Δ𝑤(𝑥, 𝑦, 𝑡) = 𝑓 (𝑥, 𝑦, 𝑡), (𝑥, 𝑦) ∈ Σ ⊂ 𝑅2 , 𝜕𝑡 with the below initial and boundary conditions: ⋃ 𝑤(𝑥, 𝑦, 0) = 𝑤0 (𝑥, 𝑦), (𝑥, 𝑦) ∈ Σ = 𝜕Σ Σ,
(1)
(2)
and 𝑤(𝑥, 𝑦, 𝑡) = ℎ(𝑥, 𝑦, 𝑡),
(𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0,
(3)
where the functions f, h and w(x, y, 0) are known in their respective domains, Δ denotes the 2D Laplacian operator, Σ is a connected and bounded domain, and 𝜕Σ represents the boundary of Σ. I𝛼 denotes the Riemann–Liouville fractional integral operator of the order 𝛼, 𝛼 ≥ 0. ∗
Also, w(x, y, t) can be represented by w(X, t). Eq. (1) can be considered as the fractional diffusion-wave equation. We can find these equations in the wave propagation modeling including heat transfer in materials with memory, viscoelastic forces, and so on [12,13]. In recent years, ample research has been performed on the progression of numerical approaches to solve the fractional evolution equations. A detailed literature review of the proposed techniques for solving a partial integro-differential equations was presented in [14]. FD approaches for a partial integro-differential equation were analyzed in [15,16]. Spectral and spectral collocation methods were studied by Chang [17]. In [18], the authors solved 2D problem (1)–(3) by using H1-Galerkin mixed finite element techniques. Furthermore, the same problem also studied in [8] and a numerical approach on the base of alternating direction implicit (ADI) method was used to approximate the solution of the problem. Two different viewpoints of strong form meshless approaches were used to study the numerical solution of fractional evolution equations in [19] as follows: First, the authors apply a numerical method on the base of convolution sum and a time discretization method to estimate the appearing fractional integral operator and time derivative. Next, they analytically proved unconditional stability of the time discretization scheme. Then, they well utilized and formulated a fully Kansa’s mesh-free technique based on the Gaussian RBF to solve the governing problem. Recently, mesh-less methods have been of much interest due to their straight forward numerical implementation without
Corresponding author. E-mail addresses: [email protected] (M. Radmanesh), [email protected] (M.J. Ebadi).
https://doi.org/10.1016/j.enganabound.2020.01.017 Received 27 August 2019; Received in revised form 29 December 2019; Accepted 30 January 2020 0955-7997/© 2020 Elsevier Ltd. All rights reserved.
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
the need for meshing or re-meshing. There are several kinds of meshless methods including mesh-less approaches on the base of weak forms and mesh-less approaches on the base of collocation methods. The approaches based on weak forms can be mesh-less local radial point interpolation method [20,21], diffuse element method [22], the element-free Galerkin (EFG) method [23–25], and mesh-less local Petrov–Galerkin method [26,27]. The approaches based on collocation methods include the mesh-less collocation approach on the basis of the RBFs [28–30]. The RBF collocation technique is too simple to be implemented but the resulting matrices of the large scale problems are fully populated and thus ill-conditioned. Different adaptive methods have been recommended to tackle the ill-conditioning problem[10,11,31]. The proposed method in [32] is a mesh-less LRBF implying a local strategy, resulting in sparse and well-conditioned matrices. For further information about the LRBF, the interested readers are referred to [12,13,15,33]. In this paper, we implement a mesh-less LRBF technique for the 2D fractional evolution equation problem (1)–(3). On the base of discretization of the problem, the time variable is discrete by means of a suitable FD scheme. In this approach, to approximate the unknown solution u in an arbitrary point x, a local region x around this point has been considered and the collocation in this subdomain is applied. This approach is described in the next sections in details. The remaining of the paper is organized as follows: In Section 2, some definitions and properties of fractional integral equations and fractional calculus which will be needed in the next sections are presented. In Section 3, the LRBF method for discretization of fractional integral equations is described. The numerical results are summarized in Section 4. In the end, the conclusions are made in Section 5.
Definition 2.5. The convergence order (C1 -O) in time variable is defined as below [11] [ ] ∥ 𝐿∞ (2𝛿𝑡, ℎ) ∥ 𝐶1 − 𝑂 = log2 . ∥ 𝐿∞ (𝛿𝑡, ℎ) ∥
2. Basic definitions
𝑤(𝑋, 𝑡 + 𝛿𝑡) − 𝑤(𝑋, 𝑡) = 𝐼 𝛼 Δ𝑤(𝑋, 𝑡) + 𝑓 (𝑋, 𝑡). 𝛿𝑡
In this section some necessary definitions and properties of fractional integral equations and fractional calculus are presented (they are taken from [7–11]). Consider f(x) as a function on (a, b).
Using the trapezoidal rule for the integral part of Eq. (4) yields to
Remark 2.6. By definition, the computational value of C1 -O indicates whether the error changes in a time step are decreasing or increasing compared to the prior time step. If the value of C1 -O is positive, it shows a decreasing error by doubling the number of time steps and the larger this positive value, the error is more reduced compared to the prior step. Also, the negative value of C1 -O means when the number of time steps becomes twice, the error is increased compared to the preceding step. 3. Discretization by LRBF-MQ An FD plan is applied for the discretization of the time dimension to solve the time-dependent equations of type (1). Then, the LRBF method is used to estimate the solution of the equation at each point of the place. Let us define 𝑡𝑘 = 𝑘𝛿𝑡,
where 𝛿𝑡 = 𝑇 ∕𝐾 denotes the time-variable step size. So, using the below backward Euler finite difference formula on Eq. (1) 𝑤𝑡 (𝑋 , 𝑡) =
𝑥 ⎧ 1 (𝑥 − 𝑡)(𝛼−1) 𝑓 (𝑡)𝑑𝑡 ⎪ 𝐼 𝑓 (𝑥) = ⎨ Γ(𝛼) ∫0 ⎪ 𝑓 (𝑥) ⎩
𝛼 > 0, 𝑥 > 0
𝐼 𝛼 Δ𝑤𝑘 (𝑋, 𝑡) = =
(4)
(6)
∫0
(7)
(𝑡𝑘 − 𝑠)(𝛼−1) Δ𝑤(𝑋, 𝑠)𝑑𝑠 (8)
𝑘 −1 ∑
(𝑡𝑘 − 𝑠𝑖 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠𝑖 )],
𝑖=1
putting Eq. (8) into Eq. (7) leads to
where I is the fractional integral and Γ(.) is the Gamma function
𝑤(𝑋, 𝑡𝑘+1 ) = 𝑤(𝑋, 𝑡𝑘 ) +
Definition 2.2. The Riemann–Liouville fractional derivative of f(x) is defined as:
𝑅𝐷 𝛼 𝑓 (𝑥)
𝑡𝑘
𝑡𝑘 𝑘 [(𝑡 − 𝑠0 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠0 ) + (𝑡𝑘 − 𝑠𝑘 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠𝑘 ) 2𝑘 +2
𝛼 = 0,
𝑥 ⎧ 1 𝑑𝑚 (𝑥 − 𝑡)(𝑚−𝛼−1) 𝑓 (𝑡)𝑑𝑡 ⎪ 𝑚 ∫ Γ( 𝑚 − 𝛼) 𝑑𝑥 =⎨ 𝑎 ⎪ 𝑓 𝑚 (𝑥) ⎩
𝑤(𝑋 , 𝑡 + 𝛿𝑡) − 𝑤(𝑋, 𝑡) , 𝛿𝑡
leads to
Definition 2.1. The Riemann–Liouville (R–L) fractional integration operator I𝛼 of order 𝛼 ≥ 0 of a function f is defined as: 𝛼
𝑘 = 1, 2, … , 𝐾,
+2
𝑘 −1 ∑
𝛿𝑡𝑡𝑘 𝑘 [(𝑡 − 𝑡0 )𝛼−1 Δ𝑤𝑛 (𝑋, 𝑡0 ) 2𝑘
(𝑡𝑘 − 𝑠𝑖 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠𝑖 )] + 𝛿𝑡𝑓 𝑘 (𝑋).
(9)
𝑖=1
𝑚 − 1 ≤ 𝛼 < 𝑚, 𝑥 > 0
In this section, the LRBF technique is used to approximate w(X, tk ) in Eq. (9). A set of 𝑁 = 𝑁Σ + 𝑁𝜕Σ nodal points is considered where NΣ is in Σ and the remaining N𝜕Σ are placed on the boundary 𝜕Σ, which can be either regularly or irregularly distributed. There is a subdomain Σ𝑥𝑙 = {𝑥𝑠 }𝑛𝑠=1 (𝑛 < 𝑁), (𝑙 = 1, 2, … , 𝑁) for each point xl where n represents the number of nodes falling in the subdomain.
𝛼 = 𝑚, (5)
where 𝑚 = ⌈𝛼⌉ is the smallest integer such that 𝛼 < m and the standard derivatives of integer order.
𝑑𝑚 𝑑𝑥𝑚
represents
Remark 3.1. The nodal points si are used in the trapezoidal integral approximation corresponding to the points ti .
Definition 2.3. The root mean square(RMS) error is defined as follows ( )1∕2 ∑ (𝑤𝑒𝑥𝑎𝑐𝑡 (𝑋𝑖 , 𝑇 ) − 𝑤𝑎𝑝𝑝𝑟𝑜𝑥 (𝑋𝑖 , 𝑇 ))2 𝑅𝑀𝑆 = , 𝑁 𝑧 ∈𝑍
Definition 3.2. The number of nodes or each collocation point plus its 𝑛 − 1 neighbors used in discretization is called stencil.
𝑖
where wapprox and wexact are respectively the approximate and exact solutions, and N is the number of testing nodes.
The accuracy of the LRBF approach will be improved by the symmetric choice of 𝑛 − 1 neighbor nodes. So, the stencil of the inside and boundary nodes can be different from each other. Thus, 𝑤(𝑋𝑠 , 𝑡𝑘+1 ) in Eq. (9) is approximated by 𝑤̃ (𝑋𝑠 , 𝑡𝑘+1 ) with the following local interpolation form:
Definition 2.4. The infinity norm error ‖e‖∞ is defined as follows 𝐿∞ (𝑢) = ‖𝑤𝑒𝑥𝑎𝑐𝑡 − 𝑤𝑎𝑝𝑝𝑟𝑜𝑥 ‖∞ = max{|𝑤𝑒𝑥𝑎𝑐𝑡 (𝑋𝑖 , 𝑇 ) − 𝑤𝑎𝑝𝑝𝑟𝑜𝑥 (𝑋𝑖 , 𝑇 )| ∶ 𝑖 = 1, 2, … , 𝑁},
̃(𝑋𝑠 , 𝑡𝑘+1 ) = 𝑤(𝑋𝑠 , 𝑡𝑘+1 ) ≃ 𝑤
where wapprox and wexact are the same as Definition 2.3. 373
𝑛 ∑ 𝑗=1
𝜆𝑘𝑗 +1 𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 ),
(10)
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 1. Considered domains for the test problems. Table 1 Comparison of 𝐶1 − 𝑂𝑠, CPU run times, and corresponding absolute errors of the proposed method and the method of [11] for Ex. 1. 𝛼
𝛿t
Error[11]
L∞
𝐶1 − 𝑂[11]
𝐶1 − 𝑂
CPU(s)[11]
CPU(s)
0.25
1/32 1/64 1/128 1/256
8.0258 × 10−3 4.4150 × 10−3 2.4137 × 10−3 1.2957 × 10−3
2.8603 × 10−3 1.4108 × 10−3 1.0057 × 10−3 3.4907 × 10−4
– 0.8622 0.8712 0.8975
– 4.3416 0.4883 1.5266
20.59 64.77 237.58 964.76
16.94 60.98 140.93 281.81
0.50
1/32 1/64 1/128 1/256
1.3183 × 10−2 6.7454 × 10−3 3.4214 × 10−3 1.7248 × 10−3
3.2520 × 10−2 1.6092 × 10−3 8.0040 × 10−4 3.9913 × 10−4
– 0.9667 0.9793 0.9882
– 4.3369 1.0076 1.0039
22.59 65.42 220.73 1128.60
16.87 63.46 141.69 279.09
0.75
1/32 1/64 1/128 1/256
1.5980 × 10−2 8.2081 × 10−3 4.1640 × 10−3 2.0977 × 10−3
3.4713 × 10−3 1.7170 × 10−3 8.5781 × 10−3 4.2805 × 10−4
– 0.9611 0.9791 0.9892
– 1.0156 2.3208 4.3248
31.37 88.72 302.64 1328.20
16.33 71.10 139.99 278.30
0.95
1/32 1/64 1/128 1/256
1.3508 × 10−2 6.8090 × 10−3 3.4163 × 10−3 1.7106 × 10−3
3.5393 × 10−3 1.7593 × 10−3 8.7702 × 10−3 4.3784 × 10−4
– 0.9883 0.9950 0.9979
− 1.0085 1.3091 4.3241
30.64 87.86 304.28 1328.50
16.92 70.98 140.43 276.50
where
Eq. (10) in all nodal points 𝑋𝑖 ; 𝑖 = 1, 2, … , 𝑛, at each stencil. The below n × n matrix form can be obtained at each stencil:
is multiquadric(MQ). The parameter c > 0 is known as the “shape parameter” describing the respective RBFs width around their centers. Also, 𝜆𝑘𝑗 is an expansion vector including unknown coefficients at the (𝑛 + 1)th time layer, and ‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 represents the distance between Xs and 𝑋𝑗𝑠 . To obtain the values of the coefficients 𝜆𝑘𝑗 , we can first evaluate
̃𝑘+1 = Φ𝜆𝑘+1 , 𝑤
√ 𝜙𝑗 (𝑋) = ‖𝑋 − 𝑋𝑗 ‖2 + 𝑐 2 ,
(11)
̃𝑘+1 = [𝑤 ̃(𝑋1 , 𝑡𝑘+1 ), 𝑤 ̃(𝑋2 , 𝑡𝑘+1 ), … , 𝑤 ̃(𝑋𝑛 , 𝑡𝑘+1 )]𝑇 , Φ = 𝑤 where [𝜙(‖𝑋𝑖 − 𝑋𝑗𝑠 ‖2 )]1≤𝑖,𝑗≤𝑛 and 𝜆𝑘+1 = [𝜆𝑘1 +1 , 𝜆𝑘2 +1 , … , 𝜆𝑘𝑛 +1 ]𝑇 . If c ≠ 0 and all the collocation points are distinct, therefore, we can calculate 374
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 2. Approximate solutions and the corresponding absolute errors with 𝑁 = 324, 𝛼 = 0.35, 𝛿𝑡 = 0.005, and 𝑇 = 1 on [0, 1]2 with the regular (a) and irregular (b) distribution collocation points with 3 × 3 stencils for Ex. 1.
the vector of coefficients 𝜆𝑘+1 in Eq. (11) by 𝑘+1
𝜆
̃𝑘+1
=Φ 𝑤 −1
.
also, the second derivative of w(x) is given as below ̃(𝑋, 𝑡𝑘+1 ) 𝜕2 𝑤 = 𝜕𝑥2
(12)
̃(𝑥𝑠 , 𝑡𝑘+1 ) in terms of the certain We rewrite the approximated solution 𝑤 nodal values 𝑤𝑘+1 (𝑋𝑗𝑠 ) at the influence domain Σ𝑋𝑠 as follows ̃(𝑋𝑠 , 𝑡𝑘+1 ) = 𝑤
𝑛 ∑ 𝑗=1
=
= Φ𝑠 Φ
−1
(13)
Ψ𝑤𝑘+1 (𝑋) = Ψ𝑤𝑘 (𝑋) +
𝑤𝑘𝑠 +1
= Ψ𝑠 𝑤𝑘𝑠 +1 ,
+2
where Φ𝑠 = (𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 )) and Ψ𝑠 = Φ𝑠 Φ−1 = [𝜓1 , 𝜓1 , … , 𝜓𝑛 ]. The functions 𝜓𝑖 , 𝑖 = 1, 2, … , 𝑛 are called the shape functions for the local RBF interpolation. We can compute ΔΨ as follows ̃(𝑋, 𝑡𝑘+1 ) = Δ𝑤
𝑛 ∑ 𝑗=1
=
𝑤𝑘+1 (𝑋𝑙 ) = ℎ𝑘+1 (𝑋𝑙 ),
(14)
A𝑊 𝑘+1 = B𝑊 𝑘 +
ΔΨ𝑠 𝑤𝑘𝑠 +1 ,
=
𝜕(Φ𝑠 Φ−1 𝑤𝑘𝑠 +1 ) 𝜕𝑥
𝜕(Φ𝑠 ) −1 𝑘+1 Φ 𝑤𝑠 , 𝜕𝑥
(17)
𝑙 = 1, 2, … , 𝑁𝜕Σ .
(18)
Hence, there exist 𝑁Σ + 𝑁𝜕Σ equations which means that we have N equations with N unknowns. Therefore, the sparse system can be shown as the below matrix form
where ΔΦ𝑠 = [Δ𝜙(‖𝑋𝑠 − 𝑋1𝑠 ‖2 ), Δ𝜙(‖𝑋𝑠 − 𝑋2𝑠 ‖2 ), … , Δ𝜙(‖𝑋𝑠 − 𝑋𝑛𝑠 ‖2 )]. The numerical derivative can be then described as ̃(𝑋, 𝑡𝑘+1 ) 𝜕𝑤 = 𝜕𝑥
𝑛−1 ∑ (𝑡𝑘 − 𝑡𝑖 )𝛼−1 ΨΔ 𝑤𝑖 (𝑋)].
By using Eq. (3) for those nodes placed on the boundary 𝜕Σ, we have
= ΔΦ𝑠 Φ−1 𝑤𝑘𝑠 +1 =
𝛿𝑡𝑡𝑘 𝑘 [(𝑡 − 𝑡0 )𝛼−1 ΨΔ 𝑤0 (𝑋) 2𝑘
𝑖=1
𝜆𝑘𝑠𝑗+1 Δ𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 )
ΔΦ𝑠 𝜆𝑘𝑠 +1
(16)
𝜕𝑥 The first or higher derivatives with respect to y can be calculated in a similar manner. By substituting the approximation formulas Eqs. (13) and (14) into Eq. (9) we have
𝜆𝑘𝑠𝑗+1 𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 )
= Φ𝑠 𝜆𝑘𝑠 +1
𝜕 2 (Φ𝑠 Φ−1 𝑤𝑘𝑠 +1 ) 𝜕𝑥2 𝜕 2 (Φ𝑠 ) −1 𝑘+1 Φ 𝑤𝑠 . 2
𝑘
∑ 𝛿𝑡𝑡𝑘 C((𝑡𝑘 − 𝑡0 )𝛼−1 𝑊 0 + (𝑡𝑘 − 𝑡𝑖 )𝛼−1 𝑊 𝑖 ) + b, 𝑘 2 𝑖=1
where the matrices A, B, C, b and W𝑘 can be defined as: { Ψ𝑖𝑗 , 𝑖 = 1, 2, … , 𝑁Σ , A= 𝛿𝑖𝑗 , 𝑖 = 𝑁Σ + 1, … , 𝑁, 𝑗 = 1, 2, … , 𝑁,
(15) B=
375
{ Ψ𝑖𝑗 , 𝑖 = 1, 2, … , 𝑁Σ , 𝛿𝑖𝑗 , 𝑖 = 𝑁Σ + 1, … , 𝑁,
𝑗 = 1, 2, … , 𝑁,
(19)
(20)
(21)
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
C=
{ ΔΨ𝑖𝑗 , 𝑖 = 1, 2, … , 𝑁Σ , 0, 𝑖 = 𝑁Σ + 1, … , 𝑁, 𝑗 = 1, 2, … , 𝑁,
(22)
b=
{ 𝑘 𝛿𝑡𝑓 (𝑥𝑖 , 𝑦𝑗 ), ℎ𝑘 (𝑥𝑖 , 𝑦𝑖 ),
(23)
𝑖 = 1, 2, … , 𝑁Σ , 𝑖 = 𝑁Σ + 1, … , 𝑁,
W𝑘 = (𝑤𝑘 (𝑥1 , 𝑦1 ), 𝑤𝑘 (𝑥2 , 𝑦2 ), … , 𝑤𝑘 (𝑥𝑁 , 𝑦𝑁 )),
(24)
where 𝛿 ij denotes the Kronecker delta function which can be defined as follows: { 1, 𝑖 = 𝑗, 𝛿𝑖,𝑗 = (25) 0, 𝑖 ≠ 𝑗. Remark 3.3. The RBF numerical solution methods of Functional equations are interpolated by all the collocation points in the whole physical domain and boundary. Such methods of interest are called global approximation and have some difficulties. First, the resultant matrices are fully populated and thus ill-conditioned which leads to unstable computation. Second, the computation of a dense matrix equation is very expensive. Third, these RBF methods are not applicable for large-scale problems [28]. In recent years, several techniques have been developed to overcome the above-mentioned difficulties. The singular value decomposition (SVD) [34] performs well to regularize the ill-conditioning of a moderate-size RBF dense interpolation matrix [35–37]. Alternatively, one could also utilize the multi-grid approach [38], the greedy algorithm [39,40], and the extended precision arithmetic [41]. To overcome these difficulties, here we have introduced a new LRBF method. This method alleviates the ill-conditioning of the resultant matrix and the costly dense matrix of the RBF interpolation. Using strong form equation, collocation approach and the matrix operations require only inversion of matrices of small size resulting in modest computational costs. Unlike the traditional global RBF collocation method, our method is highly stable by dividing the collocation of the problem in the global domain into many local regions.
Fig. 3. Absolute errors against 𝛿𝑡={1∕10, 1∕20, 1∕30, 1∕40, 1∕50, 1∕60, 1∕70, 1∕80, 1∕90, 1∕120, 1∕320, 1∕640} on the domain Σ1 with 3 × 3 stencils for Ex. 1.
4. Numerical examples In the following examples, the numerical experiments of our proposed approach on the three test problems are considered. All examples are solved using MATLAB R2015b, on a 2.4 GHz Intel(R) Core(TM) i74510U laptop running Windows 10 Ultimate with 8.00 GB main memory. We use the root mean square(RMS) error or the infinity norm error ‖e‖∞ to measure the accuracies of the numerical results. Fig. 1 which is made by using MATLAB routine ’haltonset’ presents the considered domains in the test problems. Σ1 and Σ2 refer to the regular and irregular collocation nodes distribution respectively which cover the domain [0, 1]2 . The environment of the irregular domain Σ3 is elliptic as {(𝑥, 𝑦) ∶ 𝑥2 + 4𝑦2 ≤ 1}; the irregular domain circumference Σ4 is also defined by 𝑟 = cos(2𝜃). Example 1. Consider the following inhomogeneous fractional problem which was given and solved in [11] 𝜕𝑤(𝑥, 𝑦, 𝑡) − 𝐼 𝛼 Δ𝑤(𝑥, 𝑦, 𝑡) = 𝑓 (𝑥, 𝑦, 𝑡), (𝑥, 𝑦) ∈ Σ ⊂ 𝑅2 , 𝜕𝑡 with the following initial and boundary conditions 𝑤(𝑥, 𝑦, 0) = sin 𝑥 sin 𝑦,
(26)
(𝑥, 𝑦) ∈ Σ,
and 𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 −
𝑡𝛼+1 sin 2𝑥 sin 2𝑦, Γ(𝛼 + 2)
(𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0.
The inhomogeneous function f(x, y, t) is given as follows 𝑓 (𝑥, 𝑦, 𝑡) =
𝑡𝛼 8𝑡2𝛼+1 (2 sin 𝑥 sin 𝑦 − sin 2𝑥 sin 2𝑦) − (sin 2𝑥 sin 2𝑦). Γ(𝛼 + 1) Γ(2𝛼 + 2)
The exact analytical solution of this example is
Fig. 4. RMS errors against 𝛿𝑡 = {1∕10, 1∕20, 1∕30, 1∕40, 1∕50, 1∕60, 1∕70, 1∕80, 1∕90, 1∕120, 1∕320, 1∕640} on the domain Σ1 with 𝛼 = 0.95 for Ex. 1.
𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 − 376
𝑡𝛼+1 sin 2𝑥 sin 2𝑦. Γ(𝛼 + 2)
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 5. Approximation solutions and corresponding absolute errors at several terminal times, 𝛿𝑡 = 0.002 and 𝛼 = 0.5 on the irregular domain Σ4 for Ex.1.
Table 2 Comparison of the corresponding absolute errors, convergence orders and CPU run times taken of the proposed method and method of [11] for Ex.2. 𝛼
𝛿t
Error[11]
L∞
𝐶1 − 𝑂[11]
𝐶1 − 𝑂
CPU(s)[11]
CPU(s)
0.25
1/32 1/64 1/128 1/256
1.6691 × 10−2 9.4600 × 10−3 5.1907 × 10−3 2.7870 × 10−3
8.9068 × 10−3 2.9340 × 10−3 1.4739 × 10−3 7.2988 × 10−4
– 0.8192 0.8659 0.8972
– 1.6020 0.9932 1.0139
26.88 97.51 363.28 1274.30
20.46 74.35 171.24 302.17
0.50
1/32 1/64 1/128 1/256
1.7712 × 10−2 9.2508 × 10−3 4.7523 × 10−3 2.4166 × 10−3
1.0393 × 10−2 5.0290 × 10−3 2.0721 × 10−3 6.2254 × 10−4
– 0.9371 0.9610 0.9756
– 1.0473 1.2792 1.7349
36.57 105.55 288.79 1282.90
26.48 74.67 152.31 297.51
0.75
1/32 1/64 1/128 1/256
1.8961 × 10−2 9.7475 × 10−3 4.9431 × 10−3 2.4893 × 10−3
1.2586 × 10−2 6.1602 × 10−3 3.0463 × 10−3 4.5146 × 10−4
– 0.9599 0.9796 0.9897
– 1.0308 1.0159 2.7544
37.57 76.60 283.30 1292.70
26.01 34.29 142.36 294.49
0.95
1/32 1/64 1/128 1/256
1.4164 × 10−2 7.1532 × 10−3 3.5906 × 10−3 1.7983 × 10−3
1.3266 × 10−2 6.5297 × 10−3 3.2384 × 10−3 6.8636 × 10−4
– 0.9856 0.9944 0.9976
– 1.0226 1.0117 2.2382
37.34 107.80 344.43 1097.60
25.81 72.58 158.35 264.19
𝛿𝑡 = 0.01, 𝑁 = 324, 𝛼 = 0.35 and 𝑇 = 1 on the rectangular domain [0, 1]2 with the regular and irregular distributions of collocation points. From Fig. 2, one can find that, in most cases, the regular distribution of collocation points has better accuracy than the irregular one. Fig. 3 presents the absolute errors against 𝛿t with 𝑁 = 202 and 𝑇 = 1 on the domain [0, 1]2 for several values of 𝛼. From Fig. 3, it can be found that by increasing the number of time steps, the errors reduce. RMS error against 𝛿t with 𝑁 = 202 , 𝛼 = 0.95 and 𝑇 = 1 on the domain Σ1 with 3 × 3 and 5 × 5 stencils can be seen in Fig. 4. From this figure, we can conclude that the obtained results with 3 × 3 stencil are better than those obtained with
Here, the LRBF technique is applied to solve this problem for several values of N, 𝛼, and 𝛿t at the terminal time T with the regular and irregular distributions of collocation points on the domain [0, 1]2 and on the irregular domain Σ4 . Table 1 gives the errors comparison, convergence orders and CPU run times obtained from the present approach and the approach in [11] with 𝑁 = 202 . It can be clearly seen that the proposed method is more efficient than the proposed method in [11] from the errors comparison, convergence orders and CPU run times. Fig. 2 presents the approximate solutions and the corresponding absolute errors as acquired with 377
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 6. Approximate solution and corresponding absolute error with 𝑁 = 324, 𝛿𝑡 = 0.005, 𝛼 = 0.75 and 𝑇 = 1 on [0, 1]2 with the regular (a) and irregular (Halton quasi-random point set) (b) distribution collocation points with 3 × 3 stencils for Ex.2.
the 5 × 5. Fig. 5 presents the graphs of approximate solutions and corresponding absolute errors with the values 𝛼 = 0.5, 𝑁 = 202 and 𝛿𝑡 = 0.005 on the domain Σ4 at different terminal times. So, the proposed method is an accurate technique for irregular domains. Table 1 and Figs. 2–5 show that the generated results by our proposed method have an acceptable accuracy. Example 2. Consider the below inhomogeneous fractional problem which was given and solved in [11] 𝜕𝑤(𝑥, 𝑦, 𝑡) − 𝐼 𝛼 Δ𝑤(𝑥, 𝑦, 𝑡) = 𝑓 (𝑥, 𝑦, 𝑡), (𝑥, 𝑦) ∈ Σ ⊂ 𝑅, 𝜕𝑡 with the below initial and boundary conditions 𝑤(𝑥, 𝑦, 0) = sin 𝑥 sin 𝑦 + sin 𝑥 + sin 𝑦,
(27)
(𝑥, 𝑦) ∈ Σ,
and 𝑡𝛼+1 sin 2𝑥 sin 2𝑦 Γ(𝛼 + 2) + sin 𝑥 + sin 𝑦, (𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0,
𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 −
also, the inhomogeneous function f(x, y, t) is given as below 𝑓 (𝑥, 𝑦, 𝑡) = −
𝑡𝛼 (2 sin 𝑥 sin 𝑦 − sin 2𝑥 sin 2𝑦 + sin 𝑥 + sin 𝑦) Γ(𝛼 + 1) 8𝑡2𝛼+1 (sin 2𝑥 sin 2𝑦). Γ(2𝛼 + 2)
The exact analytical solution of the problem is 𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 −
𝑡𝛼+1 sin 2𝑥 sin 2𝑦 + sin 𝑥 + sin 𝑦. Γ(𝛼 + 2)
This example is solved with the LRBF method for several values of N, 𝛼 and 𝛿t at the terminal time T with the regular and irregular distributions of the collocation nodes on the famous domain [0, 1]2 and on
Fig. 7. Absolute errors against 𝛿𝑡={1∕10, 1∕20, 1∕30, 1∕40, 1∕50, 1∕60, 1∕70, 1∕80, 1∕90, 1∕120, 1∕320, 1∕640} on the domain Σ1 with 3 × 3 stencils for Ex.2. 378
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 8. Approximate solutions and corresponding absolute errors at several terminal times, 𝛿𝑡 = 0.002 and 𝛼 = 0.25 on the irregular domain Σ3 for Ex.2.
Table 3 Numerical results taken with 𝑁 = 122 , 𝛿𝑡 = 0.001 at various time levels on the rectangular domain [0, 1]2 for Ex. 3.
the irregular domain Σ3 . Table 2. presents the obtained absolute errors, convergence orders and CPU run times for 𝑁 = 202 of nodal points, 3 × 3 stencils, and several values of 𝛿t and 𝛼 of the proposed method in comparison with the proposed method in [11]. Fig. 6 presents the approximate solutions and the corresponding absolute errors as taken with 𝑁 = 324, 𝛿𝑡 = 0.01, 𝛼 = 0.35 and 𝑇 = 1 on the domain [0, 1]2 with the regular and irregular distribution of collocation points. From Fig. 6, we can conclude that in most cases, the regular distribution of collocation points has better accuracy than the irregular one. Fig. 7 shows the absolute errors against 𝛿t with 𝑁 = 202 and 𝑇 = 1 on [0, 1]2 for several values of 𝛼. Fig. 7 illustrates that by increasing the number of time steps, the errors are reduced. Fig. 8 presents the approximate solutions and the corresponding absolute errors as taken with 𝛿𝑡 = 0.002, 𝛼 = 0.25, and 𝑁 = 202 on the domain Σ4 at different terminal times. It can be clearly seen that the proposed approach is more accurate for the irregular domains. The findings in Table 2 and Figs. 6–8 show that the present approach results in an acceptable accuracy; the obtained results by the proposed method in this article are more efficient than those by [11].
5 × 5stencil L∞
Error [14]
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
4.0447 × 10−4 9.9651 × 10−5 2.9024 × 10−4 2.2550 × 10−4 5.5171 × 10−6 7.9265 × 10−5 2.9944 × 10−5 1.9320 × 10−5 2.3318 × 10−5 8.0032 × 10−6
4.0448 × 10−4 9.9551 × 10−5 2.9043 × 10−4 2.2555 × 10−4 6.3880 × 10−6 8.1349 × 10−5 2.9840 × 10−5 1.9573 × 10−5 2.3437 × 10−5 8.0929 × 10−6
4.3675 × 10−3 6.8088 × 10−3 5.2975 × 10−3 1.6229 × 10−3 2.2194 × 10−3 4.5764 × 10−3 4.9185 × 10−3 3.6973 × 10−3 1.8209 × 10−3 1.3248 × 10−4
5
3
𝑤(𝑥, 𝑦, 𝑡) = 𝑀(𝜋 2 𝑡 2 ) sin 𝜋𝑥 sin 𝜋𝑦, (𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0, where M(.) represents the following entire function 𝑀(𝑧) =
(28)
∞ ( )−1 ∑ 3 (−1)𝑛 Γ 𝑛 + 1 𝑧𝑛 . 2 𝑛=0
The numerical results of test problem 3 are shown in Fig. 9 and Table 3. The approximate solutions and the corresponding absolute error as taken with 𝑁 = 402 , 𝛿𝑡 = 0.005 and 3 × 3 stencils at the terminal time 𝑇 = 1 on the rectangular domain [0, 1]2 are presented in Fig. 9. Table 3 shows the maximum errors with 𝑁 = 122 number of nodal points, 𝛿𝑡 = 0.001, during the time levels up to 𝑇 = 1 on the rectangular domain [0, 1]2 with 3 × 3 and 5 × 5 stencils. As
with the below initial and boundary conditions 𝑤(𝑥, 𝑦, 0) = sin 𝜋𝑥 sin 𝜋𝑦, (𝑥, 𝑦) ∈ Σ, and 5
3 × 3stencil L∞
The exact analytical solution of the problem is
Example 3. Consider the below homogeneous equation in the specific case 𝛼 = 12 [8,14] 𝑡 1 𝜕𝑤(𝑥, 𝑦, 𝑡) 1 = (𝑡 − 𝑠) 2 Δ𝑤(𝑥, 𝑦, 𝑡)𝑑𝑠, ∫ 𝜕𝑡 2 0
t
3
𝑤(𝑥, 𝑦, 𝑡) = 𝑀(𝜋 2 𝑡 2 ) sin 𝜋𝑥 sin 𝜋𝑦, (𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0. 379
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 9. Approximate solution and corresponding absolute error with 𝑁 = 402 , 𝛿𝑡 = 0.005 and 𝑇 = 1 on the domain Σ1 for Ex.3.
shown in the present example, the LRBF method with the smaller stencil is faster than the LRBF method with the larger one. The runtime with the terminal time 𝑇 = 1, 𝑁 = 122 and 𝛿𝑡 = 0.001 for the 3 × 3 stencils is 620.21 s and for the 5 × 5 stencils is 2824.81 s. The obtained results in this paper using the LRBF method are more accurate than those in [14].
[8] Li L, Xu D. Alternating direction implicit-euler method for the two-dimensional fractional evolution equation. J Comput Phys 2013;236:157–68. [9] Khebchareon M, Pani AK, Fairweather G. Alternating direction implicit Galerkin methods for an evolution equation with a positive-type memory term. J Sci Comput 2015;65(3):1166–88. [10] Pani AK, Fairweather G, Fernandes RI. Adi orthogonal spline collocation methods for parabolic partial integro differential equations. IMA J Numer Anal 2009;41(1):248–76. [11] Chen H, Xu D, Peng Y. A second order BDF alternating direction implicit difference scheme for the two-dimensional fractional evolution equation. Appl Math Model 2017;41:54–67. [12] Fujita Y. Integrodifferential equation which interpolates the heat equation and the wave equation. Osaka J Math 1990;27:309–21. [13] McLean W, Thomée V. Numerical solution of an evolution equation with a positive– type memory term. J Austral Math Soc Ser B Appl Math 1993;35(1):23–70. [14] Shivanian E, Jafarabadi A. An improved spectral meshless radial point interpolation for a class of time-dependent fractional integral equations: 2d fractional evolution equation. J Comput Appl Math 2017;325:18–33. [15] Lopez-Marcos J. A difference scheme for a nonlinear partial integrodifferential equation. SIAM J Numer Anal 1990;27(01):20–31. [16] Tang T. A finite difference scheme for partial integro-differential equations with a weakly singular kernel. Appl Numer Math 1993;11(4):309–19. [17] Kim CH, Choi UJ. Spectral collocation methods for a partial integro-differential equation with a weakly singular kernel. J Austral Math Soc Ser B Appl Math 1998;39(3):408–30. [18] Pani AK, Fairweather G. H1-Galerkin mixed finite element methods for parabolic partial integro-differential equations. IMA J Numer Anal 2002;22(2):231–52. [19] Ghehsareh HR, Bateni SH, Zaghian A. A meshfree method based on the radial basis functions for solution of two-dimensional fractional evolution equation. Eng Anal Bound Elem 2015;61:52–60. [20] Shivanian E, Khodabandehlo HR. Application of meshless local radial point interpolation (MLRPI) on a one-dimensional inverse heat conduction problem. Ain Shams Eng J 2016;7(3):993–1000. [21] Shivanian E. Analysis of meshless local radial point interpolation (MLRPI) on a nonlinear partial integro-differential equation arising in population dynamics.. Eng Anal Bound Elem 2013;37(12):1693–702. [22] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: Diffuse approximation and diffuse elements. Comput Mech 1992;10(5):307–18. [23] Fili A, Naji A, Duan Y. Coupling three-field formulation and meshless mixed Galerkin methods using radial basis functions. J Comput Appl Math 2010;234(8):2456–68. [24] Peng M, Li D, Cheng Y. The complex variable element-free Galerkin (CVEFG) method for elasto-plasticity problems. Eng Struct 2011;33(1):127–35. [25] Fu-Nong B, Dong-Ming L, Jian-Fei W, Yu-Min C. An improved complex variable element-free Galerkin method for two-dimensional elasticity problems. Chin Phys B 2012;21(2):020–204. [26] Dehghan M, Mirzaei D. Meshless local petrovgalerkin (MLPG) method for the unsteady magnetohydrodynamic (MHD) flow through pipe with arbitrary wall conductivity. Appl Numer Math 2009;59(5):1043–58. [27] Shirzadi A, Takhtabnoos F. A local meshless method for cauchy problem of elliptic PDES in annulus domains. Inverse Probl Sci Eng 2016;24(5):729–43. [28] Chen W, Fu Z-J, Chen C-S. Recent advances in radial basis function collocation methods. Springer; 2014. [29] Hosseini VR, Shivanian E, Chen W. Local integration of 2-d fractional telegraph equation via local radial point interpolant approximation. Eur Phys J Plus 2015;130(2):33. [30] Shirzadi A, Takhtabnoos F. A local meshless collocation method for solving Landaulifschitzgilbert equation. Eng Anal Bound Elem 2015;61:104–13. [31] Lubich C. Discretized fractional calculus. Inverse Probl Sci Eng 1986;17(3):704—719. [32] Chen C, Thomée V, Wahlbin LB. Finite element approximation of a parabolic integro-differential equation with a weakly singular kernel. Math Comput 1992;198(58):587–602.
Remark 4.1. In Examples 1 and 2, it can be clearly seen that the values of C1 -O are positive with increasing the number of steps and do not obey a particular order or do not follow a certain pattern. It shows that the error decreases and a better fit is achieved by increasing the number of steps. In other words, the positive values of C1 -O are the evidence of error reduction even if they may not follow any particular pattern with increasing number of steps or nodes. 5. Conclusion In this paper, we implemented a mesh-less LRBF technique for the 2D fractional evolution equation. Based on the discretization of the problem, the time variable has become discrete by means of a suitable FD scheme. We have confirmed the accuracy of LRBF method for some numerical experiments with a curved boundary which have analytical solutions. The result analysis and numerical experiments showed that the LRBF method had a good performance for curved domains with the regular and irregular distributions of collocation nodes. In addition, a smaller time step resulted in a higher numerical precision. Furthermore, the numerical investigation confirmed that the LRBF technique could act appropriately for several boundary value problems; this was particularly the case with the Dirichlet boundary condition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] Uchaikin VV. Fractional derivatives for physicists and engineers. Springer; 2013. [2] Podlubny I. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, 198. Academic Press; 1999. [3] Oldham KB, Spanier J. The fractional calculus, 111. Academic Press; 1974. [4] Bagley RL, Torvik P. A theoretical basis for the application of fractional calculus to viscoelasticity. J Rheol 1983;27:201–10. [5] Metzler R, Klafter J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J Phys A: Math Gen 2004;37(31):161–208. [6] Machado JT, Kiryakova V, Mainardi F. Recent history of fractional calculus. Commun Nonlinear Sci Numer Simul 2011;16:1140–53. [7] Chen H, Xu D, Peng Y. An alternating direction implicit fractional trapezoidal rule type difference scheme for the two-dimensional fractional evolution equation. Int J Comput Math 2015;92(10):2178–97. 380
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
[33] Safinejad M, Moghaddam MM. A local meshless RBF method for solving fractional integro-differential equations with optimal shape parameters. Ital J Pure Appl Math 2019;41:382–98. [34] Hansen P. Regularization tools: a matlab package for analysis and solution of discrete illposed problems. Numer Algor 1994;6(1):1–35. [35] Ramachandran P. Method of fundamental solutions: singular value decomposition analysis. Commun Numer Methods Eng 2002;18(11):789–801. [36] Wei T, Hon Y, Ling L. Method of fundamental solutions with regularization techniques for cauchy problems of elliptic operators. Eng Anal Bound Elem 2007;31(4):373–85. [37] Emdadi A, Kansa E, Libre N, Rahimian M, Shekarchi M. Stable PDE solution methods for large multiquadric shape parameters. CMES Comput Model Eng Sci 2008;25(1):23–41 .
[38] Fasshauer G. Solving differential equations with radial basis functions: multilevel methods and smoothing. Adv Comput Math 1999;11:139–59. [39] Hon Y, Schaback R, Zhou X. An adaptive greedy algorithm for solving large RBF collocation problems. Numer Algor 2003;32(1):13–25. [40] Schaback R, Ling L. Stable and convergent unsymmetric meshless collocation methods. SIAM J Numer Anal 2008;46:1097–115. [41] Huang C, Lee C, Cheng A. Error estimate, optimal shape factor, and high precision computation of multiquadric collocation method. Eng Anal Boundary Elem 2007;31(7):614–23.
381
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A local mesh-less collocation method for solving a class of time-dependent fractional integral equations: 2D fractional evolution equation Mehdi Radmanesh a, M.J. Ebadi b,∗ a b
Department of Mathematics, Graduate University of Advanced Technology of Kerman, Kerman, Iran Department of Mathematics, Chabahar Maritime University, Chabahar, Iran
a r t i c l e
i n f o
Keywords: Radial basis function Fractional evolution equation Local mesh-less collocation method Global RBF Finite difference
a b s t r a c t The local radial basis function (LRBF) method is an excellent tool for solving variable-order time-fractional evolution equations. This method reduces the expensive computational cost of the conventional global RBF collocation (GRBFC) approach. The present paper represents an efficient mesh-less LRBF collocation approach for solving the two-dimensional (2D) fractional evolution equation for the arbitrary fractional order in complex-shaped domains. Having been used an appropriate finite difference (FD) technique to discrete the time variable, the LRBF method is applied in order to solve the equation. Unlike the conventional GRBFC method, consideration of the local nodes in a subdomain surrounding the given local point is needed for the proposed LRBF method. This method can reduce the ill-conditioning of the global resultant RBF interpolation matrix while it is also well efficient for large scale problems. Some illustrative examples are given to demonstrate the efficiency and proficiency of the method.
1. Introduction In recent years, many studies have been carried out to solve fractional differential equations, which are valuable tools for modeling many phenomena in economics, physics, and engineering [1]. Some of their significant applications were introduced in the research works of Podlubny [2], Oldham and Spanier [3], Bagley and Trovik [4], and Metzler and Klafter [5]. Also, the authors in [6] have given several remarkable studies on the theory, analysis, and recent historical development of the fractional calculus. In this manuscript, the following 2D fractional evolution equation is considered [7–11]: 𝜕𝑤(𝑥, 𝑦, 𝑡) − 𝐼 𝛼 Δ𝑤(𝑥, 𝑦, 𝑡) = 𝑓 (𝑥, 𝑦, 𝑡), (𝑥, 𝑦) ∈ Σ ⊂ 𝑅2 , 𝜕𝑡 with the below initial and boundary conditions: ⋃ 𝑤(𝑥, 𝑦, 0) = 𝑤0 (𝑥, 𝑦), (𝑥, 𝑦) ∈ Σ = 𝜕Σ Σ,
(1)
(2)
and 𝑤(𝑥, 𝑦, 𝑡) = ℎ(𝑥, 𝑦, 𝑡),
(𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0,
(3)
where the functions f, h and w(x, y, 0) are known in their respective domains, Δ denotes the 2D Laplacian operator, Σ is a connected and bounded domain, and 𝜕Σ represents the boundary of Σ. I𝛼 denotes the Riemann–Liouville fractional integral operator of the order 𝛼, 𝛼 ≥ 0. ∗
Also, w(x, y, t) can be represented by w(X, t). Eq. (1) can be considered as the fractional diffusion-wave equation. We can find these equations in the wave propagation modeling including heat transfer in materials with memory, viscoelastic forces, and so on [12,13]. In recent years, ample research has been performed on the progression of numerical approaches to solve the fractional evolution equations. A detailed literature review of the proposed techniques for solving a partial integro-differential equations was presented in [14]. FD approaches for a partial integro-differential equation were analyzed in [15,16]. Spectral and spectral collocation methods were studied by Chang [17]. In [18], the authors solved 2D problem (1)–(3) by using H1-Galerkin mixed finite element techniques. Furthermore, the same problem also studied in [8] and a numerical approach on the base of alternating direction implicit (ADI) method was used to approximate the solution of the problem. Two different viewpoints of strong form meshless approaches were used to study the numerical solution of fractional evolution equations in [19] as follows: First, the authors apply a numerical method on the base of convolution sum and a time discretization method to estimate the appearing fractional integral operator and time derivative. Next, they analytically proved unconditional stability of the time discretization scheme. Then, they well utilized and formulated a fully Kansa’s mesh-free technique based on the Gaussian RBF to solve the governing problem. Recently, mesh-less methods have been of much interest due to their straight forward numerical implementation without
Corresponding author. E-mail addresses: [email protected] (M. Radmanesh), [email protected] (M.J. Ebadi).
https://doi.org/10.1016/j.enganabound.2020.01.017 Received 27 August 2019; Received in revised form 29 December 2019; Accepted 30 January 2020 0955-7997/© 2020 Elsevier Ltd. All rights reserved.
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
the need for meshing or re-meshing. There are several kinds of meshless methods including mesh-less approaches on the base of weak forms and mesh-less approaches on the base of collocation methods. The approaches based on weak forms can be mesh-less local radial point interpolation method [20,21], diffuse element method [22], the element-free Galerkin (EFG) method [23–25], and mesh-less local Petrov–Galerkin method [26,27]. The approaches based on collocation methods include the mesh-less collocation approach on the basis of the RBFs [28–30]. The RBF collocation technique is too simple to be implemented but the resulting matrices of the large scale problems are fully populated and thus ill-conditioned. Different adaptive methods have been recommended to tackle the ill-conditioning problem[10,11,31]. The proposed method in [32] is a mesh-less LRBF implying a local strategy, resulting in sparse and well-conditioned matrices. For further information about the LRBF, the interested readers are referred to [12,13,15,33]. In this paper, we implement a mesh-less LRBF technique for the 2D fractional evolution equation problem (1)–(3). On the base of discretization of the problem, the time variable is discrete by means of a suitable FD scheme. In this approach, to approximate the unknown solution u in an arbitrary point x, a local region x around this point has been considered and the collocation in this subdomain is applied. This approach is described in the next sections in details. The remaining of the paper is organized as follows: In Section 2, some definitions and properties of fractional integral equations and fractional calculus which will be needed in the next sections are presented. In Section 3, the LRBF method for discretization of fractional integral equations is described. The numerical results are summarized in Section 4. In the end, the conclusions are made in Section 5.
Definition 2.5. The convergence order (C1 -O) in time variable is defined as below [11] [ ] ∥ 𝐿∞ (2𝛿𝑡, ℎ) ∥ 𝐶1 − 𝑂 = log2 . ∥ 𝐿∞ (𝛿𝑡, ℎ) ∥
2. Basic definitions
𝑤(𝑋, 𝑡 + 𝛿𝑡) − 𝑤(𝑋, 𝑡) = 𝐼 𝛼 Δ𝑤(𝑋, 𝑡) + 𝑓 (𝑋, 𝑡). 𝛿𝑡
In this section some necessary definitions and properties of fractional integral equations and fractional calculus are presented (they are taken from [7–11]). Consider f(x) as a function on (a, b).
Using the trapezoidal rule for the integral part of Eq. (4) yields to
Remark 2.6. By definition, the computational value of C1 -O indicates whether the error changes in a time step are decreasing or increasing compared to the prior time step. If the value of C1 -O is positive, it shows a decreasing error by doubling the number of time steps and the larger this positive value, the error is more reduced compared to the prior step. Also, the negative value of C1 -O means when the number of time steps becomes twice, the error is increased compared to the preceding step. 3. Discretization by LRBF-MQ An FD plan is applied for the discretization of the time dimension to solve the time-dependent equations of type (1). Then, the LRBF method is used to estimate the solution of the equation at each point of the place. Let us define 𝑡𝑘 = 𝑘𝛿𝑡,
where 𝛿𝑡 = 𝑇 ∕𝐾 denotes the time-variable step size. So, using the below backward Euler finite difference formula on Eq. (1) 𝑤𝑡 (𝑋 , 𝑡) =
𝑥 ⎧ 1 (𝑥 − 𝑡)(𝛼−1) 𝑓 (𝑡)𝑑𝑡 ⎪ 𝐼 𝑓 (𝑥) = ⎨ Γ(𝛼) ∫0 ⎪ 𝑓 (𝑥) ⎩
𝛼 > 0, 𝑥 > 0
𝐼 𝛼 Δ𝑤𝑘 (𝑋, 𝑡) = =
(4)
(6)
∫0
(7)
(𝑡𝑘 − 𝑠)(𝛼−1) Δ𝑤(𝑋, 𝑠)𝑑𝑠 (8)
𝑘 −1 ∑
(𝑡𝑘 − 𝑠𝑖 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠𝑖 )],
𝑖=1
putting Eq. (8) into Eq. (7) leads to
where I is the fractional integral and Γ(.) is the Gamma function
𝑤(𝑋, 𝑡𝑘+1 ) = 𝑤(𝑋, 𝑡𝑘 ) +
Definition 2.2. The Riemann–Liouville fractional derivative of f(x) is defined as:
𝑅𝐷 𝛼 𝑓 (𝑥)
𝑡𝑘
𝑡𝑘 𝑘 [(𝑡 − 𝑠0 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠0 ) + (𝑡𝑘 − 𝑠𝑘 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠𝑘 ) 2𝑘 +2
𝛼 = 0,
𝑥 ⎧ 1 𝑑𝑚 (𝑥 − 𝑡)(𝑚−𝛼−1) 𝑓 (𝑡)𝑑𝑡 ⎪ 𝑚 ∫ Γ( 𝑚 − 𝛼) 𝑑𝑥 =⎨ 𝑎 ⎪ 𝑓 𝑚 (𝑥) ⎩
𝑤(𝑋 , 𝑡 + 𝛿𝑡) − 𝑤(𝑋, 𝑡) , 𝛿𝑡
leads to
Definition 2.1. The Riemann–Liouville (R–L) fractional integration operator I𝛼 of order 𝛼 ≥ 0 of a function f is defined as: 𝛼
𝑘 = 1, 2, … , 𝐾,
+2
𝑘 −1 ∑
𝛿𝑡𝑡𝑘 𝑘 [(𝑡 − 𝑡0 )𝛼−1 Δ𝑤𝑛 (𝑋, 𝑡0 ) 2𝑘
(𝑡𝑘 − 𝑠𝑖 )𝛼−1 Δ𝑤𝑘 (𝑋, 𝑠𝑖 )] + 𝛿𝑡𝑓 𝑘 (𝑋).
(9)
𝑖=1
𝑚 − 1 ≤ 𝛼 < 𝑚, 𝑥 > 0
In this section, the LRBF technique is used to approximate w(X, tk ) in Eq. (9). A set of 𝑁 = 𝑁Σ + 𝑁𝜕Σ nodal points is considered where NΣ is in Σ and the remaining N𝜕Σ are placed on the boundary 𝜕Σ, which can be either regularly or irregularly distributed. There is a subdomain Σ𝑥𝑙 = {𝑥𝑠 }𝑛𝑠=1 (𝑛 < 𝑁), (𝑙 = 1, 2, … , 𝑁) for each point xl where n represents the number of nodes falling in the subdomain.
𝛼 = 𝑚, (5)
where 𝑚 = ⌈𝛼⌉ is the smallest integer such that 𝛼 < m and the standard derivatives of integer order.
𝑑𝑚 𝑑𝑥𝑚
represents
Remark 3.1. The nodal points si are used in the trapezoidal integral approximation corresponding to the points ti .
Definition 2.3. The root mean square(RMS) error is defined as follows ( )1∕2 ∑ (𝑤𝑒𝑥𝑎𝑐𝑡 (𝑋𝑖 , 𝑇 ) − 𝑤𝑎𝑝𝑝𝑟𝑜𝑥 (𝑋𝑖 , 𝑇 ))2 𝑅𝑀𝑆 = , 𝑁 𝑧 ∈𝑍
Definition 3.2. The number of nodes or each collocation point plus its 𝑛 − 1 neighbors used in discretization is called stencil.
𝑖
where wapprox and wexact are respectively the approximate and exact solutions, and N is the number of testing nodes.
The accuracy of the LRBF approach will be improved by the symmetric choice of 𝑛 − 1 neighbor nodes. So, the stencil of the inside and boundary nodes can be different from each other. Thus, 𝑤(𝑋𝑠 , 𝑡𝑘+1 ) in Eq. (9) is approximated by 𝑤̃ (𝑋𝑠 , 𝑡𝑘+1 ) with the following local interpolation form:
Definition 2.4. The infinity norm error ‖e‖∞ is defined as follows 𝐿∞ (𝑢) = ‖𝑤𝑒𝑥𝑎𝑐𝑡 − 𝑤𝑎𝑝𝑝𝑟𝑜𝑥 ‖∞ = max{|𝑤𝑒𝑥𝑎𝑐𝑡 (𝑋𝑖 , 𝑇 ) − 𝑤𝑎𝑝𝑝𝑟𝑜𝑥 (𝑋𝑖 , 𝑇 )| ∶ 𝑖 = 1, 2, … , 𝑁},
̃(𝑋𝑠 , 𝑡𝑘+1 ) = 𝑤(𝑋𝑠 , 𝑡𝑘+1 ) ≃ 𝑤
where wapprox and wexact are the same as Definition 2.3. 373
𝑛 ∑ 𝑗=1
𝜆𝑘𝑗 +1 𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 ),
(10)
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 1. Considered domains for the test problems. Table 1 Comparison of 𝐶1 − 𝑂𝑠, CPU run times, and corresponding absolute errors of the proposed method and the method of [11] for Ex. 1. 𝛼
𝛿t
Error[11]
L∞
𝐶1 − 𝑂[11]
𝐶1 − 𝑂
CPU(s)[11]
CPU(s)
0.25
1/32 1/64 1/128 1/256
8.0258 × 10−3 4.4150 × 10−3 2.4137 × 10−3 1.2957 × 10−3
2.8603 × 10−3 1.4108 × 10−3 1.0057 × 10−3 3.4907 × 10−4
– 0.8622 0.8712 0.8975
– 4.3416 0.4883 1.5266
20.59 64.77 237.58 964.76
16.94 60.98 140.93 281.81
0.50
1/32 1/64 1/128 1/256
1.3183 × 10−2 6.7454 × 10−3 3.4214 × 10−3 1.7248 × 10−3
3.2520 × 10−2 1.6092 × 10−3 8.0040 × 10−4 3.9913 × 10−4
– 0.9667 0.9793 0.9882
– 4.3369 1.0076 1.0039
22.59 65.42 220.73 1128.60
16.87 63.46 141.69 279.09
0.75
1/32 1/64 1/128 1/256
1.5980 × 10−2 8.2081 × 10−3 4.1640 × 10−3 2.0977 × 10−3
3.4713 × 10−3 1.7170 × 10−3 8.5781 × 10−3 4.2805 × 10−4
– 0.9611 0.9791 0.9892
– 1.0156 2.3208 4.3248
31.37 88.72 302.64 1328.20
16.33 71.10 139.99 278.30
0.95
1/32 1/64 1/128 1/256
1.3508 × 10−2 6.8090 × 10−3 3.4163 × 10−3 1.7106 × 10−3
3.5393 × 10−3 1.7593 × 10−3 8.7702 × 10−3 4.3784 × 10−4
– 0.9883 0.9950 0.9979
− 1.0085 1.3091 4.3241
30.64 87.86 304.28 1328.50
16.92 70.98 140.43 276.50
where
Eq. (10) in all nodal points 𝑋𝑖 ; 𝑖 = 1, 2, … , 𝑛, at each stencil. The below n × n matrix form can be obtained at each stencil:
is multiquadric(MQ). The parameter c > 0 is known as the “shape parameter” describing the respective RBFs width around their centers. Also, 𝜆𝑘𝑗 is an expansion vector including unknown coefficients at the (𝑛 + 1)th time layer, and ‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 represents the distance between Xs and 𝑋𝑗𝑠 . To obtain the values of the coefficients 𝜆𝑘𝑗 , we can first evaluate
̃𝑘+1 = Φ𝜆𝑘+1 , 𝑤
√ 𝜙𝑗 (𝑋) = ‖𝑋 − 𝑋𝑗 ‖2 + 𝑐 2 ,
(11)
̃𝑘+1 = [𝑤 ̃(𝑋1 , 𝑡𝑘+1 ), 𝑤 ̃(𝑋2 , 𝑡𝑘+1 ), … , 𝑤 ̃(𝑋𝑛 , 𝑡𝑘+1 )]𝑇 , Φ = 𝑤 where [𝜙(‖𝑋𝑖 − 𝑋𝑗𝑠 ‖2 )]1≤𝑖,𝑗≤𝑛 and 𝜆𝑘+1 = [𝜆𝑘1 +1 , 𝜆𝑘2 +1 , … , 𝜆𝑘𝑛 +1 ]𝑇 . If c ≠ 0 and all the collocation points are distinct, therefore, we can calculate 374
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 2. Approximate solutions and the corresponding absolute errors with 𝑁 = 324, 𝛼 = 0.35, 𝛿𝑡 = 0.005, and 𝑇 = 1 on [0, 1]2 with the regular (a) and irregular (b) distribution collocation points with 3 × 3 stencils for Ex. 1.
the vector of coefficients 𝜆𝑘+1 in Eq. (11) by 𝑘+1
𝜆
̃𝑘+1
=Φ 𝑤 −1
.
also, the second derivative of w(x) is given as below ̃(𝑋, 𝑡𝑘+1 ) 𝜕2 𝑤 = 𝜕𝑥2
(12)
̃(𝑥𝑠 , 𝑡𝑘+1 ) in terms of the certain We rewrite the approximated solution 𝑤 nodal values 𝑤𝑘+1 (𝑋𝑗𝑠 ) at the influence domain Σ𝑋𝑠 as follows ̃(𝑋𝑠 , 𝑡𝑘+1 ) = 𝑤
𝑛 ∑ 𝑗=1
=
= Φ𝑠 Φ
−1
(13)
Ψ𝑤𝑘+1 (𝑋) = Ψ𝑤𝑘 (𝑋) +
𝑤𝑘𝑠 +1
= Ψ𝑠 𝑤𝑘𝑠 +1 ,
+2
where Φ𝑠 = (𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 )) and Ψ𝑠 = Φ𝑠 Φ−1 = [𝜓1 , 𝜓1 , … , 𝜓𝑛 ]. The functions 𝜓𝑖 , 𝑖 = 1, 2, … , 𝑛 are called the shape functions for the local RBF interpolation. We can compute ΔΨ as follows ̃(𝑋, 𝑡𝑘+1 ) = Δ𝑤
𝑛 ∑ 𝑗=1
=
𝑤𝑘+1 (𝑋𝑙 ) = ℎ𝑘+1 (𝑋𝑙 ),
(14)
A𝑊 𝑘+1 = B𝑊 𝑘 +
ΔΨ𝑠 𝑤𝑘𝑠 +1 ,
=
𝜕(Φ𝑠 Φ−1 𝑤𝑘𝑠 +1 ) 𝜕𝑥
𝜕(Φ𝑠 ) −1 𝑘+1 Φ 𝑤𝑠 , 𝜕𝑥
(17)
𝑙 = 1, 2, … , 𝑁𝜕Σ .
(18)
Hence, there exist 𝑁Σ + 𝑁𝜕Σ equations which means that we have N equations with N unknowns. Therefore, the sparse system can be shown as the below matrix form
where ΔΦ𝑠 = [Δ𝜙(‖𝑋𝑠 − 𝑋1𝑠 ‖2 ), Δ𝜙(‖𝑋𝑠 − 𝑋2𝑠 ‖2 ), … , Δ𝜙(‖𝑋𝑠 − 𝑋𝑛𝑠 ‖2 )]. The numerical derivative can be then described as ̃(𝑋, 𝑡𝑘+1 ) 𝜕𝑤 = 𝜕𝑥
𝑛−1 ∑ (𝑡𝑘 − 𝑡𝑖 )𝛼−1 ΨΔ 𝑤𝑖 (𝑋)].
By using Eq. (3) for those nodes placed on the boundary 𝜕Σ, we have
= ΔΦ𝑠 Φ−1 𝑤𝑘𝑠 +1 =
𝛿𝑡𝑡𝑘 𝑘 [(𝑡 − 𝑡0 )𝛼−1 ΨΔ 𝑤0 (𝑋) 2𝑘
𝑖=1
𝜆𝑘𝑠𝑗+1 Δ𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 )
ΔΦ𝑠 𝜆𝑘𝑠 +1
(16)
𝜕𝑥 The first or higher derivatives with respect to y can be calculated in a similar manner. By substituting the approximation formulas Eqs. (13) and (14) into Eq. (9) we have
𝜆𝑘𝑠𝑗+1 𝜙(‖𝑋𝑠 − 𝑋𝑗𝑠 ‖2 )
= Φ𝑠 𝜆𝑘𝑠 +1
𝜕 2 (Φ𝑠 Φ−1 𝑤𝑘𝑠 +1 ) 𝜕𝑥2 𝜕 2 (Φ𝑠 ) −1 𝑘+1 Φ 𝑤𝑠 . 2
𝑘
∑ 𝛿𝑡𝑡𝑘 C((𝑡𝑘 − 𝑡0 )𝛼−1 𝑊 0 + (𝑡𝑘 − 𝑡𝑖 )𝛼−1 𝑊 𝑖 ) + b, 𝑘 2 𝑖=1
where the matrices A, B, C, b and W𝑘 can be defined as: { Ψ𝑖𝑗 , 𝑖 = 1, 2, … , 𝑁Σ , A= 𝛿𝑖𝑗 , 𝑖 = 𝑁Σ + 1, … , 𝑁, 𝑗 = 1, 2, … , 𝑁,
(15) B=
375
{ Ψ𝑖𝑗 , 𝑖 = 1, 2, … , 𝑁Σ , 𝛿𝑖𝑗 , 𝑖 = 𝑁Σ + 1, … , 𝑁,
𝑗 = 1, 2, … , 𝑁,
(19)
(20)
(21)
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
C=
{ ΔΨ𝑖𝑗 , 𝑖 = 1, 2, … , 𝑁Σ , 0, 𝑖 = 𝑁Σ + 1, … , 𝑁, 𝑗 = 1, 2, … , 𝑁,
(22)
b=
{ 𝑘 𝛿𝑡𝑓 (𝑥𝑖 , 𝑦𝑗 ), ℎ𝑘 (𝑥𝑖 , 𝑦𝑖 ),
(23)
𝑖 = 1, 2, … , 𝑁Σ , 𝑖 = 𝑁Σ + 1, … , 𝑁,
W𝑘 = (𝑤𝑘 (𝑥1 , 𝑦1 ), 𝑤𝑘 (𝑥2 , 𝑦2 ), … , 𝑤𝑘 (𝑥𝑁 , 𝑦𝑁 )),
(24)
where 𝛿 ij denotes the Kronecker delta function which can be defined as follows: { 1, 𝑖 = 𝑗, 𝛿𝑖,𝑗 = (25) 0, 𝑖 ≠ 𝑗. Remark 3.3. The RBF numerical solution methods of Functional equations are interpolated by all the collocation points in the whole physical domain and boundary. Such methods of interest are called global approximation and have some difficulties. First, the resultant matrices are fully populated and thus ill-conditioned which leads to unstable computation. Second, the computation of a dense matrix equation is very expensive. Third, these RBF methods are not applicable for large-scale problems [28]. In recent years, several techniques have been developed to overcome the above-mentioned difficulties. The singular value decomposition (SVD) [34] performs well to regularize the ill-conditioning of a moderate-size RBF dense interpolation matrix [35–37]. Alternatively, one could also utilize the multi-grid approach [38], the greedy algorithm [39,40], and the extended precision arithmetic [41]. To overcome these difficulties, here we have introduced a new LRBF method. This method alleviates the ill-conditioning of the resultant matrix and the costly dense matrix of the RBF interpolation. Using strong form equation, collocation approach and the matrix operations require only inversion of matrices of small size resulting in modest computational costs. Unlike the traditional global RBF collocation method, our method is highly stable by dividing the collocation of the problem in the global domain into many local regions.
Fig. 3. Absolute errors against 𝛿𝑡={1∕10, 1∕20, 1∕30, 1∕40, 1∕50, 1∕60, 1∕70, 1∕80, 1∕90, 1∕120, 1∕320, 1∕640} on the domain Σ1 with 3 × 3 stencils for Ex. 1.
4. Numerical examples In the following examples, the numerical experiments of our proposed approach on the three test problems are considered. All examples are solved using MATLAB R2015b, on a 2.4 GHz Intel(R) Core(TM) i74510U laptop running Windows 10 Ultimate with 8.00 GB main memory. We use the root mean square(RMS) error or the infinity norm error ‖e‖∞ to measure the accuracies of the numerical results. Fig. 1 which is made by using MATLAB routine ’haltonset’ presents the considered domains in the test problems. Σ1 and Σ2 refer to the regular and irregular collocation nodes distribution respectively which cover the domain [0, 1]2 . The environment of the irregular domain Σ3 is elliptic as {(𝑥, 𝑦) ∶ 𝑥2 + 4𝑦2 ≤ 1}; the irregular domain circumference Σ4 is also defined by 𝑟 = cos(2𝜃). Example 1. Consider the following inhomogeneous fractional problem which was given and solved in [11] 𝜕𝑤(𝑥, 𝑦, 𝑡) − 𝐼 𝛼 Δ𝑤(𝑥, 𝑦, 𝑡) = 𝑓 (𝑥, 𝑦, 𝑡), (𝑥, 𝑦) ∈ Σ ⊂ 𝑅2 , 𝜕𝑡 with the following initial and boundary conditions 𝑤(𝑥, 𝑦, 0) = sin 𝑥 sin 𝑦,
(26)
(𝑥, 𝑦) ∈ Σ,
and 𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 −
𝑡𝛼+1 sin 2𝑥 sin 2𝑦, Γ(𝛼 + 2)
(𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0.
The inhomogeneous function f(x, y, t) is given as follows 𝑓 (𝑥, 𝑦, 𝑡) =
𝑡𝛼 8𝑡2𝛼+1 (2 sin 𝑥 sin 𝑦 − sin 2𝑥 sin 2𝑦) − (sin 2𝑥 sin 2𝑦). Γ(𝛼 + 1) Γ(2𝛼 + 2)
The exact analytical solution of this example is
Fig. 4. RMS errors against 𝛿𝑡 = {1∕10, 1∕20, 1∕30, 1∕40, 1∕50, 1∕60, 1∕70, 1∕80, 1∕90, 1∕120, 1∕320, 1∕640} on the domain Σ1 with 𝛼 = 0.95 for Ex. 1.
𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 − 376
𝑡𝛼+1 sin 2𝑥 sin 2𝑦. Γ(𝛼 + 2)
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 5. Approximation solutions and corresponding absolute errors at several terminal times, 𝛿𝑡 = 0.002 and 𝛼 = 0.5 on the irregular domain Σ4 for Ex.1.
Table 2 Comparison of the corresponding absolute errors, convergence orders and CPU run times taken of the proposed method and method of [11] for Ex.2. 𝛼
𝛿t
Error[11]
L∞
𝐶1 − 𝑂[11]
𝐶1 − 𝑂
CPU(s)[11]
CPU(s)
0.25
1/32 1/64 1/128 1/256
1.6691 × 10−2 9.4600 × 10−3 5.1907 × 10−3 2.7870 × 10−3
8.9068 × 10−3 2.9340 × 10−3 1.4739 × 10−3 7.2988 × 10−4
– 0.8192 0.8659 0.8972
– 1.6020 0.9932 1.0139
26.88 97.51 363.28 1274.30
20.46 74.35 171.24 302.17
0.50
1/32 1/64 1/128 1/256
1.7712 × 10−2 9.2508 × 10−3 4.7523 × 10−3 2.4166 × 10−3
1.0393 × 10−2 5.0290 × 10−3 2.0721 × 10−3 6.2254 × 10−4
– 0.9371 0.9610 0.9756
– 1.0473 1.2792 1.7349
36.57 105.55 288.79 1282.90
26.48 74.67 152.31 297.51
0.75
1/32 1/64 1/128 1/256
1.8961 × 10−2 9.7475 × 10−3 4.9431 × 10−3 2.4893 × 10−3
1.2586 × 10−2 6.1602 × 10−3 3.0463 × 10−3 4.5146 × 10−4
– 0.9599 0.9796 0.9897
– 1.0308 1.0159 2.7544
37.57 76.60 283.30 1292.70
26.01 34.29 142.36 294.49
0.95
1/32 1/64 1/128 1/256
1.4164 × 10−2 7.1532 × 10−3 3.5906 × 10−3 1.7983 × 10−3
1.3266 × 10−2 6.5297 × 10−3 3.2384 × 10−3 6.8636 × 10−4
– 0.9856 0.9944 0.9976
– 1.0226 1.0117 2.2382
37.34 107.80 344.43 1097.60
25.81 72.58 158.35 264.19
𝛿𝑡 = 0.01, 𝑁 = 324, 𝛼 = 0.35 and 𝑇 = 1 on the rectangular domain [0, 1]2 with the regular and irregular distributions of collocation points. From Fig. 2, one can find that, in most cases, the regular distribution of collocation points has better accuracy than the irregular one. Fig. 3 presents the absolute errors against 𝛿t with 𝑁 = 202 and 𝑇 = 1 on the domain [0, 1]2 for several values of 𝛼. From Fig. 3, it can be found that by increasing the number of time steps, the errors reduce. RMS error against 𝛿t with 𝑁 = 202 , 𝛼 = 0.95 and 𝑇 = 1 on the domain Σ1 with 3 × 3 and 5 × 5 stencils can be seen in Fig. 4. From this figure, we can conclude that the obtained results with 3 × 3 stencil are better than those obtained with
Here, the LRBF technique is applied to solve this problem for several values of N, 𝛼, and 𝛿t at the terminal time T with the regular and irregular distributions of collocation points on the domain [0, 1]2 and on the irregular domain Σ4 . Table 1 gives the errors comparison, convergence orders and CPU run times obtained from the present approach and the approach in [11] with 𝑁 = 202 . It can be clearly seen that the proposed method is more efficient than the proposed method in [11] from the errors comparison, convergence orders and CPU run times. Fig. 2 presents the approximate solutions and the corresponding absolute errors as acquired with 377
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 6. Approximate solution and corresponding absolute error with 𝑁 = 324, 𝛿𝑡 = 0.005, 𝛼 = 0.75 and 𝑇 = 1 on [0, 1]2 with the regular (a) and irregular (Halton quasi-random point set) (b) distribution collocation points with 3 × 3 stencils for Ex.2.
the 5 × 5. Fig. 5 presents the graphs of approximate solutions and corresponding absolute errors with the values 𝛼 = 0.5, 𝑁 = 202 and 𝛿𝑡 = 0.005 on the domain Σ4 at different terminal times. So, the proposed method is an accurate technique for irregular domains. Table 1 and Figs. 2–5 show that the generated results by our proposed method have an acceptable accuracy. Example 2. Consider the below inhomogeneous fractional problem which was given and solved in [11] 𝜕𝑤(𝑥, 𝑦, 𝑡) − 𝐼 𝛼 Δ𝑤(𝑥, 𝑦, 𝑡) = 𝑓 (𝑥, 𝑦, 𝑡), (𝑥, 𝑦) ∈ Σ ⊂ 𝑅, 𝜕𝑡 with the below initial and boundary conditions 𝑤(𝑥, 𝑦, 0) = sin 𝑥 sin 𝑦 + sin 𝑥 + sin 𝑦,
(27)
(𝑥, 𝑦) ∈ Σ,
and 𝑡𝛼+1 sin 2𝑥 sin 2𝑦 Γ(𝛼 + 2) + sin 𝑥 + sin 𝑦, (𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0,
𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 −
also, the inhomogeneous function f(x, y, t) is given as below 𝑓 (𝑥, 𝑦, 𝑡) = −
𝑡𝛼 (2 sin 𝑥 sin 𝑦 − sin 2𝑥 sin 2𝑦 + sin 𝑥 + sin 𝑦) Γ(𝛼 + 1) 8𝑡2𝛼+1 (sin 2𝑥 sin 2𝑦). Γ(2𝛼 + 2)
The exact analytical solution of the problem is 𝑤(𝑥, 𝑦, 𝑡) = sin 𝑥 sin 𝑦 −
𝑡𝛼+1 sin 2𝑥 sin 2𝑦 + sin 𝑥 + sin 𝑦. Γ(𝛼 + 2)
This example is solved with the LRBF method for several values of N, 𝛼 and 𝛿t at the terminal time T with the regular and irregular distributions of the collocation nodes on the famous domain [0, 1]2 and on
Fig. 7. Absolute errors against 𝛿𝑡={1∕10, 1∕20, 1∕30, 1∕40, 1∕50, 1∕60, 1∕70, 1∕80, 1∕90, 1∕120, 1∕320, 1∕640} on the domain Σ1 with 3 × 3 stencils for Ex.2. 378
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 8. Approximate solutions and corresponding absolute errors at several terminal times, 𝛿𝑡 = 0.002 and 𝛼 = 0.25 on the irregular domain Σ3 for Ex.2.
Table 3 Numerical results taken with 𝑁 = 122 , 𝛿𝑡 = 0.001 at various time levels on the rectangular domain [0, 1]2 for Ex. 3.
the irregular domain Σ3 . Table 2. presents the obtained absolute errors, convergence orders and CPU run times for 𝑁 = 202 of nodal points, 3 × 3 stencils, and several values of 𝛿t and 𝛼 of the proposed method in comparison with the proposed method in [11]. Fig. 6 presents the approximate solutions and the corresponding absolute errors as taken with 𝑁 = 324, 𝛿𝑡 = 0.01, 𝛼 = 0.35 and 𝑇 = 1 on the domain [0, 1]2 with the regular and irregular distribution of collocation points. From Fig. 6, we can conclude that in most cases, the regular distribution of collocation points has better accuracy than the irregular one. Fig. 7 shows the absolute errors against 𝛿t with 𝑁 = 202 and 𝑇 = 1 on [0, 1]2 for several values of 𝛼. Fig. 7 illustrates that by increasing the number of time steps, the errors are reduced. Fig. 8 presents the approximate solutions and the corresponding absolute errors as taken with 𝛿𝑡 = 0.002, 𝛼 = 0.25, and 𝑁 = 202 on the domain Σ4 at different terminal times. It can be clearly seen that the proposed approach is more accurate for the irregular domains. The findings in Table 2 and Figs. 6–8 show that the present approach results in an acceptable accuracy; the obtained results by the proposed method in this article are more efficient than those by [11].
5 × 5stencil L∞
Error [14]
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
4.0447 × 10−4 9.9651 × 10−5 2.9024 × 10−4 2.2550 × 10−4 5.5171 × 10−6 7.9265 × 10−5 2.9944 × 10−5 1.9320 × 10−5 2.3318 × 10−5 8.0032 × 10−6
4.0448 × 10−4 9.9551 × 10−5 2.9043 × 10−4 2.2555 × 10−4 6.3880 × 10−6 8.1349 × 10−5 2.9840 × 10−5 1.9573 × 10−5 2.3437 × 10−5 8.0929 × 10−6
4.3675 × 10−3 6.8088 × 10−3 5.2975 × 10−3 1.6229 × 10−3 2.2194 × 10−3 4.5764 × 10−3 4.9185 × 10−3 3.6973 × 10−3 1.8209 × 10−3 1.3248 × 10−4
5
3
𝑤(𝑥, 𝑦, 𝑡) = 𝑀(𝜋 2 𝑡 2 ) sin 𝜋𝑥 sin 𝜋𝑦, (𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0, where M(.) represents the following entire function 𝑀(𝑧) =
(28)
∞ ( )−1 ∑ 3 (−1)𝑛 Γ 𝑛 + 1 𝑧𝑛 . 2 𝑛=0
The numerical results of test problem 3 are shown in Fig. 9 and Table 3. The approximate solutions and the corresponding absolute error as taken with 𝑁 = 402 , 𝛿𝑡 = 0.005 and 3 × 3 stencils at the terminal time 𝑇 = 1 on the rectangular domain [0, 1]2 are presented in Fig. 9. Table 3 shows the maximum errors with 𝑁 = 122 number of nodal points, 𝛿𝑡 = 0.001, during the time levels up to 𝑇 = 1 on the rectangular domain [0, 1]2 with 3 × 3 and 5 × 5 stencils. As
with the below initial and boundary conditions 𝑤(𝑥, 𝑦, 0) = sin 𝜋𝑥 sin 𝜋𝑦, (𝑥, 𝑦) ∈ Σ, and 5
3 × 3stencil L∞
The exact analytical solution of the problem is
Example 3. Consider the below homogeneous equation in the specific case 𝛼 = 12 [8,14] 𝑡 1 𝜕𝑤(𝑥, 𝑦, 𝑡) 1 = (𝑡 − 𝑠) 2 Δ𝑤(𝑥, 𝑦, 𝑡)𝑑𝑠, ∫ 𝜕𝑡 2 0
t
3
𝑤(𝑥, 𝑦, 𝑡) = 𝑀(𝜋 2 𝑡 2 ) sin 𝜋𝑥 sin 𝜋𝑦, (𝑥, 𝑦) ∈ 𝜕Σ, 𝑡 > 0. 379
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
Fig. 9. Approximate solution and corresponding absolute error with 𝑁 = 402 , 𝛿𝑡 = 0.005 and 𝑇 = 1 on the domain Σ1 for Ex.3.
shown in the present example, the LRBF method with the smaller stencil is faster than the LRBF method with the larger one. The runtime with the terminal time 𝑇 = 1, 𝑁 = 122 and 𝛿𝑡 = 0.001 for the 3 × 3 stencils is 620.21 s and for the 5 × 5 stencils is 2824.81 s. The obtained results in this paper using the LRBF method are more accurate than those in [14].
[8] Li L, Xu D. Alternating direction implicit-euler method for the two-dimensional fractional evolution equation. J Comput Phys 2013;236:157–68. [9] Khebchareon M, Pani AK, Fairweather G. Alternating direction implicit Galerkin methods for an evolution equation with a positive-type memory term. J Sci Comput 2015;65(3):1166–88. [10] Pani AK, Fairweather G, Fernandes RI. Adi orthogonal spline collocation methods for parabolic partial integro differential equations. IMA J Numer Anal 2009;41(1):248–76. [11] Chen H, Xu D, Peng Y. A second order BDF alternating direction implicit difference scheme for the two-dimensional fractional evolution equation. Appl Math Model 2017;41:54–67. [12] Fujita Y. Integrodifferential equation which interpolates the heat equation and the wave equation. Osaka J Math 1990;27:309–21. [13] McLean W, Thomée V. Numerical solution of an evolution equation with a positive– type memory term. J Austral Math Soc Ser B Appl Math 1993;35(1):23–70. [14] Shivanian E, Jafarabadi A. An improved spectral meshless radial point interpolation for a class of time-dependent fractional integral equations: 2d fractional evolution equation. J Comput Appl Math 2017;325:18–33. [15] Lopez-Marcos J. A difference scheme for a nonlinear partial integrodifferential equation. SIAM J Numer Anal 1990;27(01):20–31. [16] Tang T. A finite difference scheme for partial integro-differential equations with a weakly singular kernel. Appl Numer Math 1993;11(4):309–19. [17] Kim CH, Choi UJ. Spectral collocation methods for a partial integro-differential equation with a weakly singular kernel. J Austral Math Soc Ser B Appl Math 1998;39(3):408–30. [18] Pani AK, Fairweather G. H1-Galerkin mixed finite element methods for parabolic partial integro-differential equations. IMA J Numer Anal 2002;22(2):231–52. [19] Ghehsareh HR, Bateni SH, Zaghian A. A meshfree method based on the radial basis functions for solution of two-dimensional fractional evolution equation. Eng Anal Bound Elem 2015;61:52–60. [20] Shivanian E, Khodabandehlo HR. Application of meshless local radial point interpolation (MLRPI) on a one-dimensional inverse heat conduction problem. Ain Shams Eng J 2016;7(3):993–1000. [21] Shivanian E. Analysis of meshless local radial point interpolation (MLRPI) on a nonlinear partial integro-differential equation arising in population dynamics.. Eng Anal Bound Elem 2013;37(12):1693–702. [22] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: Diffuse approximation and diffuse elements. Comput Mech 1992;10(5):307–18. [23] Fili A, Naji A, Duan Y. Coupling three-field formulation and meshless mixed Galerkin methods using radial basis functions. J Comput Appl Math 2010;234(8):2456–68. [24] Peng M, Li D, Cheng Y. The complex variable element-free Galerkin (CVEFG) method for elasto-plasticity problems. Eng Struct 2011;33(1):127–35. [25] Fu-Nong B, Dong-Ming L, Jian-Fei W, Yu-Min C. An improved complex variable element-free Galerkin method for two-dimensional elasticity problems. Chin Phys B 2012;21(2):020–204. [26] Dehghan M, Mirzaei D. Meshless local petrovgalerkin (MLPG) method for the unsteady magnetohydrodynamic (MHD) flow through pipe with arbitrary wall conductivity. Appl Numer Math 2009;59(5):1043–58. [27] Shirzadi A, Takhtabnoos F. A local meshless method for cauchy problem of elliptic PDES in annulus domains. Inverse Probl Sci Eng 2016;24(5):729–43. [28] Chen W, Fu Z-J, Chen C-S. Recent advances in radial basis function collocation methods. Springer; 2014. [29] Hosseini VR, Shivanian E, Chen W. Local integration of 2-d fractional telegraph equation via local radial point interpolant approximation. Eur Phys J Plus 2015;130(2):33. [30] Shirzadi A, Takhtabnoos F. A local meshless collocation method for solving Landaulifschitzgilbert equation. Eng Anal Bound Elem 2015;61:104–13. [31] Lubich C. Discretized fractional calculus. Inverse Probl Sci Eng 1986;17(3):704—719. [32] Chen C, Thomée V, Wahlbin LB. Finite element approximation of a parabolic integro-differential equation with a weakly singular kernel. Math Comput 1992;198(58):587–602.
Remark 4.1. In Examples 1 and 2, it can be clearly seen that the values of C1 -O are positive with increasing the number of steps and do not obey a particular order or do not follow a certain pattern. It shows that the error decreases and a better fit is achieved by increasing the number of steps. In other words, the positive values of C1 -O are the evidence of error reduction even if they may not follow any particular pattern with increasing number of steps or nodes. 5. Conclusion In this paper, we implemented a mesh-less LRBF technique for the 2D fractional evolution equation. Based on the discretization of the problem, the time variable has become discrete by means of a suitable FD scheme. We have confirmed the accuracy of LRBF method for some numerical experiments with a curved boundary which have analytical solutions. The result analysis and numerical experiments showed that the LRBF method had a good performance for curved domains with the regular and irregular distributions of collocation nodes. In addition, a smaller time step resulted in a higher numerical precision. Furthermore, the numerical investigation confirmed that the LRBF technique could act appropriately for several boundary value problems; this was particularly the case with the Dirichlet boundary condition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] Uchaikin VV. Fractional derivatives for physicists and engineers. Springer; 2013. [2] Podlubny I. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, 198. Academic Press; 1999. [3] Oldham KB, Spanier J. The fractional calculus, 111. Academic Press; 1974. [4] Bagley RL, Torvik P. A theoretical basis for the application of fractional calculus to viscoelasticity. J Rheol 1983;27:201–10. [5] Metzler R, Klafter J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J Phys A: Math Gen 2004;37(31):161–208. [6] Machado JT, Kiryakova V, Mainardi F. Recent history of fractional calculus. Commun Nonlinear Sci Numer Simul 2011;16:1140–53. [7] Chen H, Xu D, Peng Y. An alternating direction implicit fractional trapezoidal rule type difference scheme for the two-dimensional fractional evolution equation. Int J Comput Math 2015;92(10):2178–97. 380
M. Radmanesh and M.J. Ebadi
Engineering Analysis with Boundary Elements 113 (2020) 372–381
[33] Safinejad M, Moghaddam MM. A local meshless RBF method for solving fractional integro-differential equations with optimal shape parameters. Ital J Pure Appl Math 2019;41:382–98. [34] Hansen P. Regularization tools: a matlab package for analysis and solution of discrete illposed problems. Numer Algor 1994;6(1):1–35. [35] Ramachandran P. Method of fundamental solutions: singular value decomposition analysis. Commun Numer Methods Eng 2002;18(11):789–801. [36] Wei T, Hon Y, Ling L. Method of fundamental solutions with regularization techniques for cauchy problems of elliptic operators. Eng Anal Bound Elem 2007;31(4):373–85. [37] Emdadi A, Kansa E, Libre N, Rahimian M, Shekarchi M. Stable PDE solution methods for large multiquadric shape parameters. CMES Comput Model Eng Sci 2008;25(1):23–41 .
[38] Fasshauer G. Solving differential equations with radial basis functions: multilevel methods and smoothing. Adv Comput Math 1999;11:139–59. [39] Hon Y, Schaback R, Zhou X. An adaptive greedy algorithm for solving large RBF collocation problems. Numer Algor 2003;32(1):13–25. [40] Schaback R, Ling L. Stable and convergent unsymmetric meshless collocation methods. SIAM J Numer Anal 2008;46:1097–115. [41] Huang C, Lee C, Cheng A. Error estimate, optimal shape factor, and high precision computation of multiquadric collocation method. Eng Anal Boundary Elem 2007;31(7):614–23.
381