Received: 26 July 2019
DOI: 10.1002/mma.6352
S P EC I A L I S S U E PAPER
A new approach to solving multiorder time-fractional
advection–diffusion–reaction equations using BEM and
Chebyshev matrix
Moein Khalighi1,2 Mohammad Amirianmatlob1,3 Alaeddin Malek1
1Department of Applied Mathematics,
Tarbiat Modares University, Tehran, Iran
2Department of Future Technologies,
University of Turku, Turku, Finland
3Department of Mathematics and
Statistics, Dalhousie University, Halifax,
Canada
Correspondence
Moein Khalighi, Department of Applied
Mathematics, Tarbiat Modares University,
Tehran, Iran.
Email: moein.khalighi@utu.fi
Communicated by: Y. Zhou
Present Address
Moein Khalighi, Department of Future
Technologies, Faculty of Sciences and
Engineering, Yliopistonmaki, FI-20014,
University of Turku, Finland.
In this paper, the boundary element method is combined with Chebyshev oper-
ational matrix technique to solve two-dimensional multiorder time-fractional
partial differential equations: nonlinear and linear in respect to spatial and
temporal variables, respectively. Fractional derivatives are estimated by Caputo
sense. Boundary element method is used to convert the main problem into a
system of a multiorder fractional ordinary differential equation. Then, the pro-
duced system is approximated by Chebyshev operational matrix technique, and
its conditionnumber is analyzed.Accuracy and efficiency of the proposedhybrid
scheme are demonstrated by solving three different types of two-dimensional
time-fractional convection–diffusion equations numerically. The convergent
rates are calculated for different meshing within the boundary element tech-
nique. Numerical results are given by graphs and tables for solutions and
different type of error norms.
KEYWORDS
boundary element method, Chebyshev operational matrix, multiorder fractional differential
equations
MSC CLASSIFICATION
65N35; 65N38; 35R11
1 INTRODUCTION
Fractional order differential operators are the representative of nonlocal phenomena, whereas many integer-order differ-
ential operators are mostly applied to examine local phenomena1; therefore, fractional calculus can be useful to describe
many of real-world problems that cannot be covered in the classic mathematical literature.2,3 Because the next state of
many systems depends on its current and historical states, there is a great demand to improve topical methods for the
real life problems.4,5 These problems happen in bioengineering,6 solid mechanics,7 anomalous transport,8 continuum
and statistical mechanics,9 economics,10 relaxation electrochemistry,11 diffusion procedures,12 complex networks,13,14
and optimal control problems.15,16 Fractional diffusion equations are largely used in describing abnormal convection
phenomenon of liquid in medium. Models of convection–diffusion quantities play significant roles in many practi-
cal applications,4,5 especially those involving fluid flow and heat transfer, such as thermal pollution in river system,
leaching of salts in soils for computational simulations, oil reservoir simulations, transport of mass and energy, and
global weather production. Numerical methods for convection–diffusion equations described by derivatives with inte-
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the
original work is properly cited.
© 2020 The Authors. Mathematical Methods in the Applied Sciences published by John Wiley & Sons Ltd.
Math Meth Appl Sci. 2020;1–21. wileyonlinelibrary.com/journal/mma 1
2 KHALIGHI ET AL.
ger order have been studied extensively.17 Due to the mathematical complexity, analytical solutions are very few and
are restricted to the solution of simple fractional ordinary differential equations (ODEs).18 Several numerical techniques
for solving fractional partial differential equations (PDEs) have been reported, such as variational iteration,19 Ado-
mian decomposition,20 operational matrix of B-spline functions,21 operational matrix of Jacobi polynomials,12 Jacobi
collocation,22 operational matrix of Chebyshev polynomials,23,24 Legendre collocation,25 pseudo-spectral,26 and opera-
tional matrix of Laguerre polynomials,27 Pade approximation, and two-sided Laplace transformations.28 Besides finite
elements and finite differences,29 spectral methods are one of the three main methodologies for solving various types of
fractional differential equations.30-34 The main idea of spectral methods is to express the solution of differential equations
as a sum of basis functions and then to choose the coefficients in order to minimize the error between the numerical and
exact solutions as much as possible.35 Therefore, high accuracy and ease of implementing are two of the main features
that have encouraged many researchers to apply such methods to solve various types of fractional integral and differen-
tial equation. In this article, shifted chebyshev polynomials are used,36 indeed it would be promising to apply spectral tau
methods with other orthogonal polynomials, for example, Legendre37 or Jacobi.38
The boundary element method (BEM) always requires a fundamental solution to the original differential equation in
order to avoid domain integrals in the formulation of the boundary integral equation. In addition, the nonhomogeneous
and nonlinear terms are incorporated in the formulation bymeans of domain integrals. The basic idea of BEM is the trans-
formation of the original differential equation into an equivalent integral equation only on the boundary, which has been
widely applied inmany areas of engineering such as fluid mechanics,39 magnetohydrodynamic,40 and electrodynamics.41
In this paper, the BEM is developed for the numerical solution of time-fractional partial differential equations (TFPDEs)
for nonhomogeneous bodies, which converts the main problem into a system of fractional ODE with initial conditions,
described by an equation having a known fundamental solution. The proposedmethod introduces an additional unknown
domain function, which represents the fictitious source function for the equivalent problem. This function is determined
from a supplementary domain integral equation, which is converted to a boundary integral equation using a meshless
technique based on global approximation by a radial basis function (RBF) series.
The Delaunay graph mapping method can be viewed as a fast interpolation scheme. Despite its efficiency, the mesh
quality for large deformation may deteriorate near the boundary, in particular, if the deformation involves large rotation,
whichmay even lead to an invalid Delaunay graph. Furthermore, the RBFmethod can generally better preserve themesh
quality near the boundary, but the computational cost is much higher, as the mesh size increases. In order to develop
methods that are more efficient and because of their flexibility and simplicity, the Delaunay graph-based RBF method
(DG-RBF) were proposed42 to overcome the difficulty of meshing and remeshing the entire structure. Thus, the pure
boundary character of the method is maintained, because elements discretization and the integrations are limited only
to the boundary. To obtain the fictitious source, we use the Chebyshev spectral method based on operational matrix.
The primary aim of this method is to propose a suitable way to approximate linear multiorder fractional ODEs with
constant coefficients using a shifted Chebyshev Tau method that guarantees an exponential convergence speed.43 Once
the fictitious source is established, the solution of the original problem can be calculated from the integral representation
of the solution in the substituted problem.
The outline of this paper is as follows. In Section 2, we introduce the multiorder time-fractional convection diffusion
equation (TFCDE) for a class of the TFPDEs as a mathematical modeling of natural phenomena, and some basic prelimi-
naries are also given. Section 3 is devoted to applying the BEM for converting themain problem into a systemofmultiorder
fractional ODE with initial conditions. In Section 4, the Chebyshev operational matrix (COM) of fractional derivative is
obtained by applying the spectral methods to solve the generated multiorder fractional ODE. In order to demonstrate the
efficiency and accuracy of the proposed method, along with the analysis of the condition number of COM, some numer-
ical experiments are presented in Section 6 using the definitions and lemmas of Section 5. Eventually, we conclude the
paper with some remarks and add the appendix to make more convenient understanding of the proposed algorithm.
2 PROBLEM STATEMENT
Assume that we are given the following initial boundary value problem for the multiorder time-fractional PDE in the
two-dimensional domain Ω with boundary 𝛤 ,
k∑
𝑗=0
𝛾𝑗(x)D
𝛼𝑗
c u = A(x)uxx + 2B(x)ux𝑦 + C(x)u𝑦𝑦
+ D(x)ux + E(x)u𝑦 + F(x)u + g(x, t),
(1)
KHALIGHI ET AL. 3
where A(x), B(x), C(x), D(x), E(x), F(x), and 𝛾 j(x) for j = 0, 1, 2, … , k and g(x, t) are specified functions their physical
meaning depends on that of the field function u(x, t), and
0 < 𝛼0 < 𝛼1 < … < 𝛼k−1 < 𝛼k, m − 1 < 𝛼k ⩽ m,
x = (x, 𝑦) ∈ Ω ∪ Γ, t > 0,
subject to the boundary conditions
B(u) = h(x, t), x ∈ Γ, (2)
and the initial conditions
Dicu(x, 0) = di(x), i = 0, 1, … ,m − 1 (3)
in whichm is an integer number and D𝛼𝑗c u is the Caputo fractional time derivative of order 𝛼j. The Caputo derivative,12 is
employed because initial conditions having direct physical meaning can be prescribed. This derivative is defined as
D𝛼c u(x, 𝑦, t) =
⎧⎪⎪⎨⎪⎪⎩
1
Γ(m−𝛼)
t
∫
0
u(m)(x,𝑦,𝜏)d𝜏
(t−𝜏)1+𝛼−m , m − 1 < 𝛼 < m,
dm
dtm u(x, 𝑦, t), 𝛼 = m ∈ N.
(4)
B(·) is a linear operator with respect to spatial variables x, y of order one. h(x, t) and di(x) (i = 0, … ,m− 1) are specified
functions in Equations (2) and (3), respectively. It seems that we could be able to recover the multiterm of classical diffu-
sion equation for 𝛼k = 1, 𝛾 j(x) = 0 for j = 0, … , k − 1, 𝛾k(x) ≠ 0 and the classical wave equation in presence of viscous
damping for 𝛼k = 2, 𝛼k−1 = 1, 𝛾k(x) ≠ 0, 𝛾k−1(x) ≠ 0 and 𝛾 j(x) = 0 for j = 0, … , k − 2.
3 IMPLEMENTATION OF BOUNDARY ELEMENT METHOD
Taking advantage of the following boundary element, the initial boundary value of the Equations (1–3) is reformed into
an ODE problem.
Let u(x, t) be the sought solution to the problem (1–3) and assume that u is twice continuously differentiable inΩ. After
applying Laplace operator on u, we have39
∇2u(x, t) = 𝔅(x, t), (5)
where𝔅(x, t) known as an unknown fictitious source function. That is the solution of Equation (1) could be established
by solving Equation (5) under the boundary condition (2), if convenient𝔅(x, t) is first established. This is accomplished
by adhering to the following procedure.
We write the solution of Equation (5) in the integral form. Thus, for un, as the normal derivative of u, we have44
𝜀u(x, t) = ∫
Ω
u∗𝔅dΩ − ∫
Γ
(u∗un − u∗nu)dΓ, x ∈ Ω ∪ Γ, (6)
where u∗ = ln r∕2𝜋 is the fundamental solution to Equation (5), r is the distance between any two points, and also u∗n
stands for its normal derivative on the boundary. 𝜀 is the free term coefficient taking the values 𝜀 = 1 if x ∈ Ω, 𝜀 = 𝜃∕2𝜋
if x ∈ 𝛤 , otherwise 𝜀 = 0; 𝜃 is the interior angle between the tangents of boundary at point x. 𝜀 = 1∕2 for points where
the boundary is smooth. After applying Equation (6), to boundary points by means of Greens second identity,45 we yield
the boundary integral equation
𝜃
2𝜋 u(x, t) =
M∑
𝑗=1
b𝑗
⎡⎢⎢⎣12û𝑗 + ∫Γ
(
u∗(ûn)𝑗 − u∗nû𝑗
)
dΓ
⎤⎥⎥⎦ − ∫Γ (u∗un − u∗nu)dΓ, x ∈ Γ, (7)
where û𝑗 is a certain solution of the equation
∇2û𝑗 = 𝑓𝑗, 𝑗 = 1, 2, … ,M, (8)
4 KHALIGHI ET AL.
AlsoM is the number of interior points inside Ω. Here bj are the coefficients that must be determined to satisfy
𝔅(x, t) =
M∑
𝑗=1
b𝑗𝑓𝑗 ,
where fj = fj(r), r = ||x − x𝑗|| is a set of radial basis approximating functions; xj are collocation points in Ω. The radial
basis function method is used to map the nodes rather than that based on surface or volume ratios.42 The algorithm is set
out in the following procedure. At first, we generate the Delaunay graph by using all the boundary nodes of the original
mesh, and then we locate the mesh points in the graph, after that, wemove the Delaunay graph according to the specified
geometric motion/deformation, and the final step is mapping the mesh points in the new graph according to the RBF
matrix andDelaunay triangle. Procedures before the last step are exactly the same as the original Delaunay graphmapping
method;46 hence, the details of these steps are skipped in this paper. The difference is in the last step, where the radial
basis function interpolation is used to calculate the displacement of the internal mesh nodes from the given displacement
of the Delaunay triangle nodes on the boundary, while the original Delaunay mapping method uses surface or volume
ratios to calculate the displacement of inner nodes. Equation (7) is solved numerically by using the BEM. The boundary
integrals in Equation (7) are approximated usingN boundary nodal points. Here, 𝜀 = 1∕2, as we ace the smooth boundary.
The domain integral can be evaluated when the fictitious source is estimated by a radial basis function series, and
subsequently, it is reformed to a boundary line integral using the GreenâA˘Z´s reciprocal identity.45 For the sake of simpli-
fication, we use multiquadric radial basis function in practice.M internal nodals are here used to define Delaunay linear
triangular elements in Ω. Therefore, after discretization and application of the boundary integral Equation (7) at the N
boundary nodal points, we have
Hu = Gun + Āb, (9)
where H and G are N × N known as coefficient matrices originating from the integration of the kernel functions on the
boundary elements and Ā is an N ×M coefficient matrix originating from the integration of the kernel function on the
domain elements; u and un are vectors containing the nodal values of the boundary displacements and their normal
derivatives. Also, b is the vector of the nodal values of the fictitious source at theM internal nodal points.
For a second order differential equation, the boundary condition is a correlation of 𝛿1(x)u + 𝛿2(x)un = h(x, t); after
applying it at the N boundary nodal points yields
𝛿1u + 𝛿2un = h(t), (10)
where 𝜹1 and 𝜹2 are N × N known diagonal matrices and h(t) =[h
(
xBP1 , t
)
, … , h(xBPN , t)]T is a known boundary vector,
where xBP𝑗 are N boundary nodal points. Equations (9) and (10) can be combined and solved for u and un. This yields
u =
[
(𝛿1+ 𝛿2G−1H)−1(𝛿2G−1 „A)
]
b + (𝛿1 + 𝛿2G−1H)−1h(t),
un =
[
(𝛿1H−1G + 𝛿2)−1(−𝛿1H−1 „A)
]
b + (𝛿1H−1G + 𝛿2)−1h(t).
(11)
Further, differentiating Equation (6) for points inside the domain (𝜀 = 1) with respect to x and y, using the same
discretization and collocating at theM internal nodal points, we have the following expression for the spatial derivatives
ûpq = Ĥpqu + Ĝpqun + Âpqb, p, q = 0, x, 𝑦 (12)
where the ûpq is vector of values for u and its derivatives at the M internal nodal points; Ĥpq and Ĝpq are M × N known
coefficientmatrices originating from the integration of the kernel functions on the boundary elements and Âpq is anM×M
coefficient matrix originating from the integration of the kernel functions on the domain elements.
Eliminating u and un from Equation (12) using Equations (11) yields
ûpq = Upqb + cpq, p, q = 0, x, 𝑦 (13)
where
Upq = Ĥpq(𝛿1 + 𝛿2G−1H)−1
(
𝛿2G−1 „A
)
+ Ĝpq
(
𝛿1H−1G + 𝛿2
)−1 (−𝛿1H−1 „A) + Âpq,
cpq =
[
Ĥpq(𝛿1 + 𝛿2G−1H)−1 + Ĝpq(𝛿1H−1G + 𝛿2)−1
]
h(t).
(14)
KHALIGHI ET AL. 5
The final step of the method is to apply Equation (1) at theM internal nodal points. This gives
k∑
𝑗=0
𝛾𝑗D
𝛼𝑗
c û = Aûxx + 2Bûx𝑦 + Cû𝑦𝑦 +Dûx + Eû𝑦 + Fû + g(t), (15)
where û = û00 and 𝜸j, A, B, C, D, E, and F are M × M known diagonal matrices including the nodal values of the
corresponding functions 𝛾 j(x)A(x), B(x), C(x),D(x), E(x), and F(x), respectively, and g(t) =
[
g
(
xIP1 , t
)
, … , g
(
xIPM , t
)]T is
a known internal vector, where xIP𝑗 areM internal nodal points. Substituting the corresponding terms from Equation (13)
into Equation (15) yields
k∑
𝑗=0
S𝑗D
𝛼𝑗
c b(t) = Nb(t) + f(t), (16)
where
S𝑗 = 𝛾𝑗U,
N = AUxx + 2BUx𝑦 + CU𝑦𝑦 +DUx + EU𝑦 + FU,
f(t) = Acxx + 2Bcx𝑦 + Cc𝑦𝑦 +Dcx + Ec𝑦 + Fc + g(t) −
k∑
𝑗=0
𝛾𝑗D
𝛼𝑗
c (c),
(17)
in which U = U00 and c = c00 for j = 0, … , k. Now, from Equation (13), we can write the initial conditions (3) for
i = 0, 1, … ,m − 1 in the form
b(i)(0) = U−1
[
di − c(i)(0)
]
, (18)
where di(t) =
[
di
(
xIP1 , t
)
, … , di
(
xIPM , t
)]T .
The above proposed procedure reduces the problem of multiorder two-dimensional time-fractional PDE (1–3) to a
simpler system of multiterm fractional ODE (16) with initial condition (18). The existence, uniqueness, and continuous
dependence of the system Equations (16–18) when Sk = 1 can be rigorously discussed (see, e.g., Diethelm and Neville's
paper47). In the next section, we show the implementation of Chebyshev operational matrix, as a spectral technique36 for
fractional calculus, to solve the system of initial value problem (16–18).
4 COM METHOD FOR SYSTEM OF MULTIORDER FRACTIONAL ODES
The Chebyshev polynomials Ti(t) are defined on the interval (−1, 1).43 Thus, by changing variable t → 2tL − 1, the shifted
Chebyshev polynomialsTL,i(t) of degree i on the interval t ∈ (0,L), with an orthogonality relation can be introduced by36,48
TL,i(t) = i
i∑
𝑗=0
(−1)i−𝑗 (i + 𝑗 − 1)!2
2𝑗
(i − 𝑗)!(2𝑗)!L𝑗 t
𝑗 ,
where TL,i(0) = (−1)i and TL,i(L) = 1. In this form, TL,i(t)may be generated by the following recurrence formula:
TL,i+1(t) = 2(2t∕L − 1)TL,i(t) − TL,i−1(t), i = 1, 2, ... (19)
where TL,0(t) = 1 and TL,1(t) = 2t∕L − 1. Therefore, a given function f ∈ L2[0, 1]may be approximated by K + 1 terms of
shifted Chebyshev polynomials as
𝑓 (t) ≃ 𝑓K(t) =
K∑
i=0
ciTL,i(t),
where the coefficients ci are described by weight functions wL(t) = 1√Lt−t2 as c𝑗 = 1hi ∫ L0 𝑓 (t)TL,i(t)wL(t)dt; in which hi = 𝜋
for i = 0, otherwise hi = 𝜋2 . If we set
Φ(t) =
[
TL,0(t),TL,1(t), … ,TL,K(t)
]T (20)
and suppose 𝜐 > 0 and the ceiling function ⌈𝜐⌉ denotes the smallest integer greater than or equal to 𝜐, then
D𝜐cΦ(t) ≃ 𝔇(𝜐)Φ(t), (21)
6 KHALIGHI ET AL.
where𝔇(𝜐) is the (K + 1) × (K + 1) COM of derivatives of order 𝜐 in the Caputo sense and is defined by36,48
𝔇(𝜐) =
⎛⎜⎜⎜⎜⎜⎜⎜⎝
0 0 0 … 0
⋮ ⋮ ⋮ … ⋮
0 0 0 … 0
S𝜐 (⌈𝜐⌉ , 0) S𝜐 (⌈𝜐⌉ , 1) S𝜐 (⌈𝜐⌉ , 2) … S𝜐 (⌈𝜐⌉ ,K)
⋮ ⋮ ⋮ … ⋮
S𝜐 (i, 0) S𝜐 (i, 1) S𝜐 (i, 2) … S𝜐 (i,K)
⋮ ⋮ ⋮ … ⋮
S𝜐 (K, 0) S𝜐 (K, 1) S𝜐 (K, 2) … S𝜐 (K,K)
⎞⎟⎟⎟⎟⎟⎟⎟⎠
, (22)
where
S𝜐(i, 𝑗) =
i∑
𝜆=⌈𝜐⌉
(−1)i−𝜆2i(i + 𝜆 − 1)!Γ
(
𝜆 − 𝜐 + 12
)
𝜌𝑗L𝜐Γ
(
𝜆 + 12
)
(i − 𝜆)!Γ(𝜆 − 𝜐 − 𝑗 + 1)Γ(𝜆 + 𝑗 − 𝜐 + 1)
,
where 𝜌0 = 2, 𝜌𝜆 = 1, 𝜆 ≥ 1. Note that in𝔇(𝜐), the first ⌈𝜐⌉ rows are all zero. In order to solve Equation (16) with initial
conditions (18), we approximate b(t) and f(t) in terms of shifted Chebyshev polynomials as
b(t) =
⎡⎢⎢⎢⎢⎢⎢⎣
K∑
i=0
𝜓1i TL,i(t)
⋮
K∑
i=0
𝜓Mi TL,i(t)
⎤⎥⎥⎥⎥⎥⎥⎦
=
⎡⎢⎢⎣
Ψ1
T
Φ(t)
⋮
ΨM
T
Φ(t)
⎤⎥⎥⎦ = ΨΦ(t), (23)
f (t) =
⎡⎢⎢⎢⎢⎢⎢⎣
K∑
i=0
𝓁1i TL,i(t)
⋮
K∑
i=0
𝓁Mi TL,i(t)
⎤⎥⎥⎥⎥⎥⎥⎦
=
⎡⎢⎢⎣
𝜔1
T
Φ(t)
⋮
𝜔M
T
Φ(t)
⎤⎥⎥⎦ = 𝜔Φ(t), (24)
where for j = 1, … ,M
Ψ𝑗
T
=
[
𝜓
𝑗
0 , … , 𝜓
𝑗
K
]
,
𝜔𝑗
T
=
[
𝓁𝑗0, … ,𝓁
𝑗
K
]
,
and Ψ, 𝝎 areM × (K + 1)matrices that are defined as
Ψ =
⎡⎢⎢⎣
Ψ1
T
⋮
ΨM
T
⎤⎥⎥⎦ , 𝜔 =
⎡⎢⎢⎣
𝜔1
T
⋮
𝜔M
T
⎤⎥⎥⎦ .
For j = 0, … , k, Equations (21) and (23) can be used to write
D𝛼𝑗c b(t) ≃ ΨD
𝛼𝑗
c Φ(t) ≃ Ψ𝔇(𝛼𝑗 )Φ(t). (25)
Employing Equations (23), (24), and (25), then the residual for Equation (16) can be written as
R(t) =
( k∑
𝑗=0
S𝑗Ψ𝔇(𝛼𝑗 ) −NΨ − 𝜔
)
Φ(t). (26)
R(t) is a M vector with respect to t. If Rj(t) be the jth component of R(t), then in a typical Tau method,43 we generate
M(K −m + 1) linear equations withM(K + 1) unknown coefficients of Ψ by applying
⟨
Ri(t),TL,𝑗(t)
⟩
= ∫
L
0
Ri(t)TL,𝑗(t)dt = 0, i = 1, … ,M, 𝑗 = 0, 1, … ,K −m. (27)
KHALIGHI ET AL. 7
Also, by substituting Equation (23) into Equation (18), and with the fact that
𝔇(n) = (𝔇(1))n, n ∈ N,
we get a system ofMm linear equations withM(K + 1) unknown coefficients for Ψ as following
b(i)(0) = Ψ𝔇(i)Φ(0) = U−1
[
di − c(i)(0)
]
, i = 0, 1, … ,m − 1. (28)
Equations (27) and (28) can be rewritten in the matrix form
AΨ = B, (29)
where A is anM(K+ 1) ×M(K+ 1) coefficient matrix. The system of algebraic Equations (29) can be easily solved for the
unknown vector Ψ. Consequently, b(t) given in Equation (23) can be calculated, which gives a solution of Equation (16)
with the initial conditions (18). Once the vector b(t) of the values of the fictitious source at the M internal nodal points
has been established, then the solution of Equation (1) and its derivatives can be computed from Equation (13). For the
points x = (x, y) that do not coincide with the prespecified internal nodal points, the solution could be drawn from the
discretized counterpart of Equation (6)with 𝜀 = 1 using the same boundary andnewdomain discretization.Note that here
the matrices Ĥpq, Ĝpq, and Âpq corresponding to previous internal nodes plus the additional points must be recomputed.
5 ERROR ESTIMATION
The convergence of the proposed method is shown by employing the following error norms, maximum error (L∞), max-
imum relative error (MRE) to assess the accuracy of the method in multiscale problems and root mean square (RMS) to
globally examine the method efficiency,
L∞ = max1⩽i⩽M
|||uiex − uiapp||| , (30)
MRE = max
1⩽i⩽M
|||||
uiex − uiapp
uiex
||||| , (31)
RMS =
√
1
M
M∑
i=1
(
uiex − uiapp
)2
, (32)
where uiex and uiapp denote the ith components of the exact and approximated solutions, respectively, andM denotes the
number of internal points. It is not convenient to certainly determine what is the convergence rate of the proposed hybrid
method; for example, for the number of nodal points N andM, and the size of the shifted Chebyshev polynomials K, the
convergence order of the method would be O(f(rp, 𝜏q)), where f is a function of the convergence rate of BEM with the
order p and the convergence rate of COM with the order q. The accuracy of the method depends on several factors, the
convergence speed for BEM, the domain and boundary discretization, the shape parameters of radial basis functions, the
orders of the derivatives, and condition number of COMs, as such the analysis of truncation errors for methods solving a
two-dimensional multiterm time-fractional differential equation is not straightforward. Nonetheless, information about
matrixA from the algebraic system (29), particularly its condition number, will be useful. The condition number is defined
by43
Cond(A) = max{|𝜆| ∶ det(A − 𝜆I) = 0}min{|𝜆| ∶ det(A − 𝜆I) = 0} ,
such that amatrix with a large condition number is so-called ill conditioned, whereas thematrix is namedwell conditioned
if its condition number is of a moderate size. We also suggest two P orders in the following lemmas to examine the
rate of convergence for BEM and COM distinctly. The first one is directly tested by the exact solution and the effect of
domain discretization, whereas the second one is addressed by comparing a sequence of numerical solutions of the ODE
system (16) with different degree sizes of COMs that have been offered exponential rates of convergence accuracy for
smooth problems in simple geometries.43
Lemma5.1. Let the vectorU be the exact solution of the initial boundary value problem (1–3) andu1,u2 the approximate
solutionswithN1,M1 andN2,M2 of nodal points, respectively. Then, the computational order of the BEMmethodproposed
8 KHALIGHI ET AL.
in Section 3 can be calculated with Pr order≃ log
(
Er1
/
Er2
)
∕log (r1 ∕ r2) in which Er1 and Er2 are corresponded RMS
errors (32) with the relative boundary mesh size r1 = 1∕N1 and r2 = 1∕N2, respectively.
Proof. When the leading terms in the spatial-discretization error are proportional to r1p and r2p, and ‖.‖RMS denoting
the root mean square norm (32),
‖U − u1‖RMS = Er1 = c1rp1 ≃ crp1 , ‖U − u2‖RMS = Er2 = c2rp2 ≃ crp2 .
Hence,
Er1
Er2
≃
rp1
rp2
,
then taking logarithm from both sides yields
p ≃
log
(
Er1
Er2
)
log
(
r1
r2
) .
Lemma 5.2. Let the vector bex be the exact solution and b1,b2, and b3 be the approximate solutions of the multiterm
fractional ODE (16) with the initial condition (18) at the sameM fictitious source points using K1, K2, and K3 the numbers
of shifted Chebyshev polynomials, respectively. With considering this proportion
K1
K2
= K2K3
, (33)
the temporal convergence order for the COM method presented in Section 4 is estimated using P𝜏 order≃
log
(‖b2−b1‖b‖b3−b2‖ b
)/
log
(
𝜏1
𝜏2
)
in which the norm ‖.‖b define as ‖bs − bt‖b =
√√√√ 1
M
M∑
i=1
(
bis − bit
)2, where bis and bit denote the
ith components, and 𝜏1 = 1 ∕ K1, 𝜏2 = 1 ∕ K2, and 𝜏3 = 1 ∕ K3.
Proof. When the leading terms in the error of COM are proportional to 𝜏1p, 𝜏2p, and 𝜏3p,
‖bex − b1‖b = c′1𝜏p1 ≃ c′𝜏p1 , ‖bex − b2‖b = c′2𝜏p2 ≃ c′𝜏p2 , ‖bex − b3‖b = c′3𝜏p3 ≃ c′𝜏p3 ,
thus
bex ≃
‖‖‖𝜏p2b1 − 𝜏p1b2‖‖‖b
𝜏
p
2 − 𝜏
p
1
, bex ≃
‖‖‖𝜏p3b2 − 𝜏p2b3‖‖‖b
𝜏
p
3 − 𝜏
p
2
,
according to Equation (33), we have ‖b2 − b1‖b‖b3 − b2‖b ≃
(
𝜏1
𝜏2
)p
.
Hence,
p ≃
log
(‖b2−b1‖b‖b3−b2‖b)
log
(
𝜏1
𝜏2
) .
In the following section, the numerical errors are computed based on assumptions described in Lemmas 5.1 and 5.2.
6 NUMERICAL RESULTS AND DISCUSSION
On the basis of the described procedure, some problems are solved to illustrate the efficiency and the accuracy of the pro-
posedmethod. In the first example, a simple two-dimensional fractional heat-like equation is considered for two different
conditions. In the second example, a nonlinear two-dimensional fractional wave-like equation is tested. In the third and
fourth problems, two linear TFCDEs are solved to test the impact of external force (g) and the final time on the convergence
KHALIGHI ET AL. 9
FIGURE 1 The location of nodal points and relative distances for constant element discretization [Colour figure can be viewed at
wileyonlinelibrary.com]
rate of the method. In the fifth and sixth test problems, multiorder time-fractional diffusion-wave equations in bounded
homogeneous anisotropic plane bodies are solved. The condition number of system 29 is examined for each example.
Because the size of the matrix A depends on the number of internal points,M, and the degree of COM, K, the condition
number of A can be compared versus K and the length of distance between nodal points. However, most domains are not
discretized uniformly. In this regard, suppose r denotes the mean length of all the distances between the internal points
and their adjacent points (e.g., see xIP in Figure 1). Thus, the numerical results show that the condition number behaves
as CondA ≃ r−2(K + 1)2 for Example 6.1 and CondA ≃ r−2(K + 1)3 for other examples.
Example 6.1. Consider the following two-dimensional time-fractional heat-like equation:
D𝛼c u = uxx + u𝑦𝑦, 0 < 𝛼 ⩽ 1, t > 0
subjected to different initial conditions with different domains49,50:
(I) u(x, 𝑦, 0) = sin x sin 𝑦,0 < x, 𝑦 < 2𝜋, (II)
u(x, 𝑦, 0) = cos
(
𝜋
2 x
)
cos
(
𝜋
2 𝑦
)
.
0 < x, 𝑦 < 1,
Here, boundary conditions satisfy the exact solutions:
(I) u(x, 𝑦, t) = E𝛼 (−2t𝛼) sin x sin 𝑦,
(II) u(x, 𝑦, t) = E𝛼
(
−12𝜋
2t𝛼
)
cos
(
𝜋
2 x
)
cos
(
𝜋
2 𝑦
)
,
(34)
where the following one-parameter Mittag-Leffler function is defined as
E𝛼(z) =
∞∑
k=0
zk
Γ(𝛼k + 1) , 𝛼 > 0.
Figure 2 demonstrates L∞ errors andMRE versus the degree K (right) for Example 6.1 (case I) when N = 160 and 𝛼 = 0.5
at t = 1.5. The convergence rate of COM is estimated that P𝜏 order> 5. Furthermore, the condition numbers of A, on
the Figure 2 (left), are shown versus the polynomial degree, K, and the mean length of discrete elements, r. It can be
numerically deduced that the condition number behaves as Cond(A) ≃ r−2(K + 1)2; for example, when K = 8, then
Cond(A) ≃ 76.4× r−2, and when r = 0.3, then Cond(A) ≃ 10.43×K2. In Table 1, numerical results are compared with the
exact solutions (34) for Example 6.1, case I, for fixed K = 12, t = 1.5, with the differential orders 𝛼 = 0.5 and 𝛼 = 0.75, and
different number of nodal points; the convergence rate of BEM is algebraic (Pr order> 4.4) when the number of nodal
points is increased from N = 40 to N = 80 and it is quadratic (Pr order> 2) when N = 80 goes to N = 160. Apart from the
value of 𝛼, it can be inferred that the computation cost of the second discretization for moderateN andM is more effective
10 KHALIGHI ET AL.
FIGURE 2 The condition number of A versus the polynomial degree K and the mean length of discrete elements r (left) and a comparison
of L∞ errors andMRE versus K (right) at t = 1.5 when N = 160 and 𝛼 = 0.5 for Example 6.1 case I; P𝜏 order= 5.13 is estimated for K = 4, 8, 16
[Colour figure can be viewed at wileyonlinelibrary.com]
FIGURE 3 The condition number of
A versus the polynomial degree K and
the mean length of discrete elements r
(left) and the relative absolute errors
(right) obtained for Example 6.1 case II
when N = 200, K = 12, 𝛼 = 0.5, and
t = 0.5 [Colour figure can be viewed at
wileyonlinelibrary.com]
TABLE 1 The error norms and
the estimated order of convergence
Pr for the vector solutionU
according to the Lemma 5.1, in
Example 6.1 case I when K = 12
and t = 1.5
N 𝛼 = 0.5 𝛼 = 0.75
MRE RMS Pr order MRE RMS Pr order
40 5.48066 × 10−3 2.38289 × 10−4 − 8.21377 × 10−3 3.57120 × 10−4 −
80 2.52922 × 10−4 1.09966 × 10−5 4.4376 3.72158 × 10−4 1.61807 × 10−5 4.4641
160 6.10497 × 10−5 2.65433 × 10−6 2.0506 8.43489 × 10−5 3.66734 × 10−6 2.1415
than the third one. For case (II), a similar behavior of Cond(A) versus K and r is shown in Figure 3 (left). The relative
absolute error (right) withN = 200, and K = 12 for the final time t = 1 is exhibited. Intuitively, the relative absolute errors
are approached to 10−5, which it could be expected for K = 12 and N = 200 based on the information from Table 2. This
table shows the estimated convergence for two terms of shifted Chebyshev series for the final time t = 0.5; in general,
there is an improvement for errors when the degree K increases, but no relationship between Pr order and degree K is
observed. It may also be concluded N = 64 is more computational cost-effective in this case.
Example 6.2. Consider the two-dimensional time-fractional wave-like equation51:
D𝛼c u =
1
12
(
x2uxx + 𝑦2u𝑦𝑦
)
, 0 < x, 𝑦 < 1, 1 < 𝛼 ⩽ 2, t > 0,
subjected to boundary conditions
u(0, 𝑦, t) = 0, u(1, 𝑦, t) = 4 cosh t,
u(x, 0, t) = 0, u(x, 1, t) = 4 sinh t,
and the initial condition
u(x, 0) = x4,
u′(x, 0) = 𝑦4.
KHALIGHI ET AL. 11
N K = 12 K = 18
L∞ RMS Pr order L∞ RMS Pr order
8 1.9970 × 10−3 1.8930 × 10−3 − 1.3222 × 10−3 1.0553 × 10−3 −
16 1.4673 × 10−3 7.1121 × 10−4 1.4123 1.0571 × 10−3 5.1238 × 10−4 1.0424
32 6.0501 × 10−4 1.2220 × 10−4 2.5410 3.5453 × 10−4 7.1607 × 10−5 2.8390
64 9.3468 × 10−5 9.9095 × 10−6 3.6243 6.9059 × 10−5 7.3216 × 10−6 3.2899
128 3.8455 × 10−5 1.8098 × 10−6 2.4530 2.4404 × 10−5 1.1486 × 10−6 2.6723
TABLE 2 A comparison of the error
norms and the estimated order of
convergence Pr for the vector solutionU
according to the Lemma 5.1, in Example
6.1 case II for two fixed K
K L∞ MRE RMS ‖bi − bi−1‖b P𝜏 order Cond(A)
8 7.15174 × 10−3 2.25751 × 10−2 2.91855 × 10−3 − − 4.54024 × 105
16 2.82029 × 10−5 8.43598 × 10−5 1.06609 × 10−5 5.5885 × 10−6 − 2.96189 × 106
32 1.07633 × 10−5 2.08835 × 10−5 6.63201 × 10−6 5.0329 × 10−7 3.4730 2.14632 × 107
64 8.17112 × 10−6 5.69019 × 10−6 9.87053 × 10−7 3.2380 × 10−8 3.9582 1.71381 × 108
TABLE 3 The error norms for
the vector solutionU and the
convergence orders of the vector
solution b in Example 6.2 for
𝛼 = 2, N = 100 and t = 0.5
The exact solution for 𝛼 = 2 is found to be
u(x, 𝑦, t) = x4 cosh t + 𝑦4 sinh t. (35)
This problem is solved for differentKwithN = 100 at t = 0.5 and integer order 𝛼 = 2 to compare with the exact solution
(35). The results from Table 3 offer an error improvement by increasing the number of degree K of COMs. The condition
numbers of the system (29) illustrate such behavior Cond(A) ≃ r−2(K + 1)3 (see Table 3). Due to the fact that the domain
and boundary nodal points are fixed here, the numerical solutions of U are directly affected by the numerical solutions
of b; in other words, affected by the accuracy of COM. However, the exact solution of the generated ODE system is not
clear, and the convergence order of COMs is estimated by norm ‖.‖b for three distinct degrees with a same proportion.
Interestingly, a comparison between the two columns RMS and ‖.‖b of Table 3 suggests direct relationships but with
different speed between the approximation solution ofU and b. In addition, by considering the scale of the solutions, and
a comprehensive assessment of the absolute errors for distinct degrees K, it could be concluded that the Chebyshev Tau
method converges with an oscillating manner around the exact solution of fractional ODE system (16).
Example 6.3. Consider the following TFCDE52:
D𝛼c u = uxx + u𝑦𝑦 − 5ux − 5u𝑦 + g(x, 𝑦, t),
1 < 𝛼 ⩽ 2, (x, 𝑦) ∈ Ω,
with the boundary condition and initial conditions
u(x, 𝑦, t) = 0, (x, 𝑦) ∈ Γ, 0 < t ⩽ 1,
u(x, 𝑦, 0) = 0, 𝜓(x, 𝑦) = 0, (x, 𝑦) ∈ Ω,
where Ω = [0, 1]2. Then we have the following exact solution:
u(x, 𝑦, t) = 212t2+𝛼x3(1 − x)3𝑦3(1 − 𝑦)3.
It is easy to check
g(x, 𝑦, t) = 211Γ(𝛼 + 3)x3(1 − x)3𝑦3(1 − 𝑦)3t2
− 212(6x − 51x2 + 120x3 − 150x3 − 150x4 + 30x5)𝑦3(1 − 𝑦)3t2+𝛼
− 212x3(1 − x)3(6𝑦 − 51𝑦2 + 120𝑦3 − 105𝑦4 + 30𝑦5)t2+𝛼.
This problem is challenging and sensitive because of the large numbers included in the function g. However, it could be
compensated by multiplying to power functions of decimal numbers, and considering the final time t = 1. Figure 4 (left
plan) shows the estimated error ranged around 10−4, for 𝛼 = 1.5 and final time t = 1, with the degree K = 10, and (right
plan) demonstrates the plot of the error versus the number of boundary nodes, N, with K = 10 for three different values
𝛼, illustrating that the smoothness roughly occurred after N = 135. The behavior of the condition number matrix A is
estimated as Cond(A) ≃ r−2(K + 1)3. Table 4 gives Cond(A), the RMS error and the convergence rates are obtained by
12 KHALIGHI ET AL.
FIGURE 4 Graphs of the absolute error with 𝛼 = 1.5 and K = 10 (left) and a comparison of errors for different values of 𝛼 (right) at finite
time t = 1 for Example 6.3. [Colour figure can be viewed at wileyonlinelibrary.com]
TABLE 4 The RMS error of the vector
solutionU and the convergence rate of
COM and spatial-discretization base on
Lemmas 5.1 and 5.2 for Example 6.3 at
time t = 1
K=16 𝜶=1.9 𝜶=1.1
N RMS Cond(A) Pr order RMS Cond(A) Pr order
40 8.8704 × 10−3 4.86686 × 105 − 7.3330 × 10−3 4.89036 × 105 −
80 2.0840 × 10−3 1.95421 × 106 2.0897 1.6591 × 10−3 1.96064 × 106 2.1439
160 4.8292 × 10−4 7.78881 × 106 2.1095 3.5514 × 10−4 7.80051 × 106 2.2240
320 8.1966 × 10−5 3.13058 × 107 2.5587 7.2452 × 10−5 3.13134 × 107 2.2933
N=200 𝜶=1.9 𝜶=1.6
K RMS ‖bi − bi−1‖b P𝜏 order RMS ‖bi − bi−1‖b P𝜏 order
10 1.3010 × 10−4 − − 7.5314 × 10−5 − −
20 4.2187 × 10−5 1.3840 × 10−6 − 1.6938 × 10−5 1.4839 × 10−6 −
40 1.2781 × 10−5 3.7002 × 10−7 1.9032 3.8515 × 10−6 2.8012 × 10−7 2.4053
80 3.7745 × 10−6 9.8207 × 10−8 1.9137 8.8808 × 10−7 4.8922 × 10−8 2.5175
solving Example 6.3 for different values of 𝛼. It indicates better RMS errors for 𝛼 near to 1 than 2, that is not true for Pr
orders.
Example 6.4. Consider the linear TFCDE53:
D𝛼c u = uxx + u𝑦𝑦 − 0.1ux − 0.1u𝑦 + g(x, 𝑦, t),
1 < 𝛼 ⩽ 2, (x, 𝑦) ∈ Ω,
with the initial conditions
u(x, 𝑦, 0) = 0, u′(x, 𝑦, 0) = 0, (x, 𝑦) ∈ Ω,
KHALIGHI ET AL. 13
FIGURE 5 Graphs of approximation solution and the absolute error obtained for Example 6.4 for 𝛼 = 1.05 and 𝛼 = 1.95 with K = 10 and
N = 255 for final time t = 2 [Colour figure can be viewed at wileyonlinelibrary.com]
RMS − error
N 𝛼 = 1.1 𝛼 = 1.4 𝛼 = 1.7 𝛼 = 1.9
10 1.5030 × 10−1 1.0628 × 10−1 1.2802 × 10−1 1.1669 × 10−1
20 3.5223 × 10−2 3.8318 × 10−2 5.0480 × 10−2 4.7404 × 10−2
40 1.6236 × 10−2 1.1808 × 10−2 1.4409 × 10−2 1.1033 × 10−2
80 2.4654 × 10−3 4.9065 × 10−3 2.0146 × 10−3 4.5839 × 10−3
160 1.2483 × 10−3 9.5059 × 10−4 9.5252 × 10−4 1.1706 × 10−3
320 3.7848 × 10−4 2.0804 × 10−4 1.7544 × 10−4 2.1400 × 10−4
TABLE 5 The RMS error for the vector solutionU, in
Example 6.4 for K = 8 and t = 2
and the Dirichlet boundary conditions. The exact solution of the current test problem is
u(x, 𝑦, t) = t3+𝛼 sin
(
𝜋
6 x
)
sin
(7𝜋
4 x
)
sin
(3𝜋
4 𝑦
)
sin
(5𝜋
4 𝑦
)
,
whereΩ is the computational domain as shown in Figure 5 (left plan). The approximate solution and its relative absolute
error are shown in Figure 5 for 𝛼 = 1.05 and 𝛼 = 1.95 with K = 10 and N = 255 for the final time t = 2. Although
N = 255 is considered to depict clear data points in the figure, N = 150 would be sufficient to achieve the semiequivalent
errors. Importantly, by considering the results of the previous example 6.3 and Figure 5 (right plan), it may convey that
the method has a better performance for the less values of 𝛼. In contrast, Table 5 refuses this idea; there are irrelevant
outcomes versus the values of 𝛼, although the table shows a reliable numerical convergence.
Example 6.5. The multiorder time-fractional diffusion-wave equation44
14 KHALIGHI ET AL.
FIGURE 6 The geometry of the plane
body discretization with N = 210,
M = 132 with DG-RBF technique and
contour plots of absolute errors when
K = 16 for Example 6.5 [Colour figure
can be viewed at wileyonlinelibrary.com]
TABLE 6 The condition number of A, the
RMS error, and the convergent order for the
vector solution b according to the Lemma
5.2 in Example 6.5 for K = 16 and t = 4
RMS − error
N u ux uxy Pr order r ≃ Cond(A)
50 6.27668 × 10−2 6.69339 × 10−2 6.37693 × 10−2 − 2.4 8.5723 × 102
100 6.74090 × 10−3 5.97502 × 10−3 6.03444 × 10−3 3.40157 1.2 3.4289 × 103
120 3.70614 × 10−3 3.60496 × 10−3 3.86762 × 10−3 2.43988 0.1 4.9376 × 103
200 7.71852 × 10−4 7.85751 × 10−4 7.03289 × 10−4 3.33700 0.6 1.3715 × 104
FIGURE 7 The condition number of matrix A versus the polynomial degree K (left) and versus the mean length of discrete element edges
r (right) for Example 6.6 when t = 5
KHALIGHI ET AL. 15
𝛾1D1.7c u + 𝛾0D0.8c u = Auxx + 2Bux𝑦 + Cu𝑦𝑦 + g(x, t), x(x, 𝑦) ∈ Ω, t > 0, (36)
in the plane inhomogeneous anisotropic body that is shown in Figure 6 has been solved, subject to boundary conditions
u(x, t) = 0, x(x, 𝑦) ∈ Γ,
and the initial condition
u(x, 0) = 0,
u′(x, 0) = U(x, 𝑦),
where A = (𝑦
2−x2+50)
50 , B =
2x𝑦
50 , C =
(x2−𝑦2+50)
50 , 𝛾1 = 5e
(−0.1(|x|+|𝑦|)), 𝛾0 = 0.4(x2 + 𝑦2)1∕2. The external source g is g(x, t) =
U(𝛾1D1.7c T + 𝛾0D0.8c T) − T(AUxx + 2BUx𝑦 + CU𝑦𝑦), where U(x, 𝑦) = a2b2 −
((
x
a
)2
+
(
𝑦
b
)2)(( x
b
)2
+
(
𝑦
a
)2)
and T(t) =
t − t
3
6 +
t5
200 . The boundary of the domain is defined by the curve:
Γ = (ab)
1∕2((
cos 𝜃a
)2
+
(
sin 𝜃b
)2)(1∕4)((
cos 𝜃b
)2
+
(
sin 𝜃a
)2)(1∕4) ,
where 0 ⩽ 𝜃 ⩽ 2𝜋, a = 3, b = 1.3. The problem admits an exact solution uexact = T(t)U(x, y). This problem is solved
using BEM and COM for variousM and N when K = 16 at t = 4. The value of u, ux, and uxy are compared with the exact
solution. Numerical results are given in Table 6 showing the efficiency of the proposed method by Pr order > 2, and the
condition number of matrix A with the manner as r−2 × (K + 1)3. In Figure 6, the contour plots illustrate the absolute
errors distributions of the approximations of u, ux, and uxy on the plane, including L∞,MRE,RMS, at t = 4 for specific
nodal points and Delaunay triangulation when N = 210,M = 132, K = 16.
Example 6.6. Consider the large terms time-fractional diffusion equation
6∑
𝑗=0
𝛾𝑗(x)D
𝛼𝑗
c u = D(x)ux + E(x)u𝑦 + F(x)u + g(x, t), x(x, 𝑦) ∈ Ω, t > 0, (37)
in a ‘C-shape’ made by the elimination of a circle with radius r2 = 3 and null origin from the inside of a circle with radius
r1 = 5 and the same origin, and extracting the space between the lines y = −1 and y = 1 from the right side of the outcome
(see Figure A1) with the Dirichlet boundary conditions, and the initial condition
16 KHALIGHI ET AL.
u(x, 0) = 0,
u′(x, 0) = U(x),
where U(x) = cos(x)esin(y) and the external force g is g(x, t) = UT − T(DUx + EU𝑦 + FU) such that T(t) = t
3
6 − t, and
with the order of the derivatives 𝛼6 = 53 , 𝛼5 =
7
5 , 𝛼4 =
4
3 , 𝛼3 =
6
5 , 𝛼2 =
3
2 , 𝛼1 =
2
3 , and 𝛼0 =
1
2 , and their coefficients
𝛾6 = Γ
(
1
3
)
, 𝛾5 = 4𝜋, 𝛾4 = 2𝜋, 𝛾3 = 4𝜋, 𝛾2 =
√
𝜋, 𝛾1 =
Γ
(
1
3
)
3 , 𝛾0 =
√
𝜋
2 , D =
r22−x
2
r21−𝑦2
, E = x
2−𝑦2
r1−r2
, and F =
√
x2 + 𝑦2 + r3, we
can set
T(t) = t
5
3
⎛⎜⎜⎜⎝
2𝜋
Γ
(
8
3
) + Γ
(
1
3
)
3Γ
(
11
3
) t⎞⎟⎟⎟⎠ + t
3
2
(20 − 8t
15
)
+ t
⎛⎜⎜⎜⎝
9
4 t
1
3 + 4𝜋
Γ
(
13
5
) t 35 + 4𝜋t 45
Γ
(
14
5
)⎞⎟⎟⎟⎠
−
(
t
−2
3 + t
1
3 + t
1
2 + t
−1
2 +
√
3Γ
(1
3
)
t
−1
3 +
(√
10 + 2
√
5
)
Γ
(2
5
)
t
−2
5 +
(√
10 − 2
√
5
)
Γ
(1
5
)
t
−1
5
)
to find the exact solution as uexact = T(t)U(x). This problem is solved with N = 216, M = 741, and K = 16, for t = 5. In
Figure 7, the distribution of the absolute errors on the domain,L∞,MRE, andRMS of u, uxx, and uyy are illustrated. Figure 6
exhibits the behavior of the condition numbermatrixA as r−2×(K+1)3; for example, when r = 1,Cond(A) ≃ 0.92×(K+1)3,
and when K = 10, Cond(A) ≃ 1202.46 × r−2. Among these six exams, an interesting point can be concluded that the
error distribution into a plan depends not only on the positions of the nodal points but also on the final time solving
the problem; with an increasingly asymmetric discretization, and longer computing time, the error distribution becomes
more ‘random’.
7 CONCLUSION
Here, we have proposed a hybrid algorithm to solve two-dimensional multiorder time-fractional partial differential
equations. Their general form is given in Equations (1–3). The method consists of the boundary element method com-
bined with spectral Chebyshev operational matrix. The BEM is used to transfer the corresponding time-fractional PDE
into a system ofODEs, whereas COM is used to solve the system efficiently. Thismethod is applied to the two-dimensional
fractional heat-like, wave-like, and diffusion-wave equations, which shows that the errors of the approximate solution
decay exponentially. When the exact solution exists, comparison is made with ‖‖Uex − uapp‖‖, and the convergence rate
is calculated using Lemma 5.1. When the exact solution is not in hand, the order of convergence is estimated by three
approximate solutions with various degrees of Chebyshev polynomials in a same grid-point based on Lemma 5.2. By
applying the assumptions of the Lemmas, the numerical results show the efficiency and convergence rate for the proposed
hybrid method. Notwithstanding, it is not easy to emphasize a unique conclusion for the accuracy of the method on the
ground that given the vast range of architectures, spectral methods, boundary element methods, fractional calculus, and
meshing used with in such a hybrid-technique framework. In general, for multiorder two-dimensional time-fractional
PDE (1–3), the current method calculate U for the test exams, with this range of convergence rate: 2 < Pr order< 4.5 for
the moderate values of N andM. And for multiterm fractional ODE (16) with initial condition (18), the current method
works well to calculate solution b with a range of the convergence rate around 1 < P𝜏 order< 5.5. Moreover, The condi-
tion number of matrix A from linear system (29) behaves like CondA ≃ r−2(K + 1)2 for the problems with ⌈𝛼⌉ = 1, and
CondA ≃ r−2(K + 1)3 for the problems with ⌈𝛼⌉ = 2. For the future direction, the authors believe that establishing new
methods to examine long-term effects of memory in complex systems modeled by fractional calculus is highly required
as fractional calculus is a proper mathematical tool for describing memory,14 although the proposed COM technique is
not an appropriate scheme for long-term problems. It is instead efficient for problems with multiterm orders.
ACKNOWLEDGEMENT
The authors are very grateful to the referees for carefully reading the paper and for their comments and suggestions which
have improved the paper.
KHALIGHI ET AL. 17
8 CONFLICT OF INTEREST
The authors declare no conflict of interest.
REFERENCES
1. Teodoro GS, Machado JAT, De Oliveira EC. A review of definitions of fractional derivatives and other operators. J Comput Phys.
2019;388:195-208.
2. Metzler R, Klafter J. The random walk's guide to anomalous diffusion: a fractional dynamics approach. Phys Reports. 2000;339(1):
1-77.
3. MatlobMA, Jamali Y. The concepts and applications of fractional order differential calculus in modeling of viscoelastic systems: a primer.
Crit Rev Biomed Eng. 2019;47(4).
4. Sun HG, Zhang Y, Baleanu D, Chen W, Chen YQ. A new collection of real world applications of fractional calculus in science and
engineering. Commun Nonlinear Sci Numer Simul. 2018;64:213-231.
5. Almeida R, Bastos NR, Monteiro MTT. Modeling some real phenomena by fractional differential equations. Math Methods Appl Sci.
2016;39(16):4846-4855.
6. Magin RL. Fractional calculus in bioengineering Begell House Redding; 2006.
7. Rossikhin YA, Shitikova MV. Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of
solids. Appl Mech Rev. 1997;50(1):15-67.
8. 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):R161.
9. Mainardi F. Fractional Calculus: Springer; 1997;291-348.
10. Baillie RT. Long memory processes and fractional integration in econometrics. J Econ. 1996;73(1):5-59.
11. Oldham KB. Fractional differential equations in electrochemistry. Adv Eng Softw. 2010;41(1):9-12.
12. Podlubny I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of
Their Solution and Some of Their Applications Elsevier; 1998.
13. Safdari H, Kamali MZ, Shirazi A, Khalighi M, Jafari G, Ausloos M, Fractional dynamics of network growth constrained by aging node
interactions. PLOS one. 2016;11(5):e0154983.
14. Saeedian M, Khalighi M, Azimi-Tafreshi N, Jafari GR, Ausloos M. Memory effects on epidemic evolution: the
susceptible-infected-recovered epidemic model. Phys Rev E. 2017;95(2):022409.
15. Zeid SS, Kamyad AV, Effati S, Rakhshan SA, Hosseinpour S. Numerical solutions for solving a class of fractional optimal control problems
via fixed-point approach. SeMA J. 2017;74(4):585-603.
16. Hosseinpour S, Nazemi A. A collocation method via block-pulse functions for solving delay fractional optimal control problems. IMA J
Math Control Info. 2016;34(4):1215-1237.
17. Heidari H, Malek A. Null boundary controllability for hyperdiffusion equation. Internat J Appl Math. 2009;22:615-626.
18. Kilbas AAA, Srivastava HM, Trujillo JJ. Theory and Applications of Fractional Differential Equations: Elsevier Science Limited; 2006.
19. DehghanM, Yousefi SA, Lotfi A. The use of He's variational iteration method for solving the telegraph and fractional telegraph equations.
Int J Numer Methods Biomed Eng. 2011;27(2):219-231.
20. Ford NJ, Connolly JA. Systems-based decomposition schemes for the approximate solution of multi-term fractional differential equations.
J Comput Appl Math. 2009;229(2):382-391.
21. Lakestani M, Dehghan M, Irandoust-Pakchin S. The construction of operational matrix of fractional derivatives using B-spline functions.
Commun Nonlinear Sci Numer Sim. 2012;17(3):1149-1162.
22. Lin Y, Xu Ch. Finite difference/spectral approximations for the time-fractional diffusion equation. J Comput Phys. 2007;225(2):
1533-1552.
23. Taher AHS, Malek A, Momeni-Masuleh SH. Chebyshev differentiation matrices for efficient computation of the eigenvalues of
fourth-order Sturm–Liouville problems. Appl Math Model. 2013;37(7):4634-4642.
24. Bhrawy AH, Alofi AS. The operational matrix of fractional integration for shifted Chebyshev polynomials. Appl Math Lett.
2013;26(1):25-31.
25. Mokhtary P. Reconstruction of exponentially rate of convergence to Legendre collocation solution of a class of fractional
integro-differential equations. J Comput Appl Math. 2015;279:145-158.
26. Esmaeili S, Shamsi M. A pseudo-spectral scheme for the approximate solution of a family of fractional differential equations. Commun
Nonlinear Sci Numer Sim. 2011;16(9):3646-3654.
27. Abdelkawy MA, Taha TM. An operational matrix of fractional derivatives of Laguerre polynomials. Walailak J Sci Technol(WJST).
2014;11(12):1041-1055.
28. Hosseinpour S, Nazemi A, Tohidi E. A new approach for solving a class of delay fractional partial differential equations.Mediterranean J
Math. 2018;15(6):218.
29. Li Ch, Zeng F. Finite difference methods for fractional differential equations. Int J Bifurcation Chaos. 2012;22(04):1230014.
30. Bhrawy AH, Zaky MA. Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. Nonlinear Dyn.
2015;80(1-2):101-116.
18 KHALIGHI ET AL.
31. Bhrawy AH, Zaky MA, Van Gorder RA. A space-time Legendre spectral tau method for the two-sided space-time Caputo fractional
diffusion-wave equation. Numer Algorithms. 2016;71(1):151-180.
32. Zaky MA, Hendy AS, Macías-Díaz JE. Semi-implicit Galerkin–Legendre spectral schemes for nonlinear time-space fractional
diffusion–reaction equations with smooth and nonsmooth solutions. J Sci Comput. 2020;82(1):1-27.
33. Zaky MA, Machado JT. Multi-dimensional spectral tau methods for distributed-order fractional diffusion equations. Comput Math Appli.
2020;79(2):476-488.
34. Hafez RM, Zaky MA. High-order continuous Galerkin methods for multi-dimensional advection–reaction–diffusion problems. Eng
Comput:1-17.
35. Zaky MA. A Legendre spectral quadrature tau method for the multi-term time-fractional diffusion equations. Comput Appl Math.
2018;37(3):3525-3538.
36. Bhrawy AH, Taha TM, Machado JAT. A review of operational matrices and spectral techniques for fractional calculus. Nonlinear Dyn.
2015;81(3):1023-1052.
37. Zaky MA. An improved tau method for the multi-dimensional fractional Rayleigh–Stokes problem for a heated generalized second grade
fluid. Comput Math Appli. 2018;75(7):2243-2258.
38. Bhrawy AH, Zaky MA. A method based on the Jacobi tau approximation for solving multi-term time–space fractional partial differential
equations. J Comput Phys. 2015;281:876-895.
39. Katsikadelis JT. The Boundary Element Method for Engineers and Scientists: Theory and Applications Academic Press; 2016.
40. Morovati V, Malek A. Solving inhomogeneous magnetohydrodynamic flow equations in an infinite region using boundary element
method. Eng Anal Bound Elem. 2015;58:202-221.
41. Brebbia CA, Walker S. Boundary Element Techniques in Engineering. Elsevier; 2016.
42. Wang Y, Qin N, Zhao N. Delaunay graph and radial basis function for fast quality mesh deformation. J Comput Phys. 2015;294:149-172.
43. Canuto C, Hussaini MY, Quarteroni A, et al. Spectral Methods in Fluid Dynamics: Springer Science & Business Media; 2012.
44. Katsikadelis JT. The BEM for numerical solution of partial fractional differential equations. Comput Math Appli. 2011;62(3):891-901.
45. Katsikadelis JT. The BEM for nonhomogeneous bodies. Archive Appl Mech. 2005;74(11-12):780-789.
46. Liu X, Qin N, Xia H. Fast dynamic grid deformation based on Delaunay graph mapping. J Comput Phys. 2006;211(2):405-423.
47. Diethelm K, Ford NJ. Multi-order fractional differential equations and their numerical solution. Appl Math Comput. 2004;154(3):
621-640.
48. Doha EH, Bhrawy AH, Ezz-Eldien SS. A Chebyshev spectral method based on operational matrix for initial and boundary value problems
of fractional order. Comput Math Appli. 2011;62(5):2364-2373.
49. Momani S. Analytical approximate solution for fractional heat-like and wave-like equations with variable coefficients using the decom-
position method. Appl Math Comput. 2005;165(2):459-472.
50. Shirzadi A, Ling L, Abbasbandy S. Meshless simulations of the two-dimensional fractional-time convection–diffusion–reaction equations.
Eng Anal Bound Elem. 2012;36(11):1522-1527.
51. Shou DH, method HEJH. Beyond Adomian The variational iteration method for solving heat-like and wave-like equations with variable
coefficients. Phys Lett A. 2008;372(3):233-237.
52. Wang Z, Vong S. Comput Math Appl. 2014;68(3):185-196.
53. Dehghan M, Safarpoor M. The dual reciprocity boundary elements method for the linear and nonlinear two-dimensional time-fractional
partial differential equations.Math Method Appl Sci. 2016;39(14):3979-3995.
How to cite this article: Khalighi M, Amirianmatlob M, Malek A. A new approach to solving multiorder
time-fractional advection–diffusion–reaction equations using BEM and Chebyshev matrix. Math Meth Appl Sci.
2020;1–21. https://doi.org/10.1002/mma.6352
APPENDIX A: ALGORITHM IN A NUTSHELL
For convenience, a brief description of the proposed algorithm is given in this section.
At the beginning, determine N boundary points and divideM internal nodal points into groups by the Delaunay graph,
then use the RBF method to interpolate the nodal points to its new position. Hence, there is no need to optimize shape
parameters for the domain discretization. After discretization, consider ∫k indicating integration on k element on the
boundary (see Figure A1 and Table A1) and set i, k = 1, … ,N, j = 1, … ,M for the boundary points xBPi , the domain
points xIP𝑗 , and the algorithm implementation. Notice that different 𝜃 must be considered for computing the boundary
integrals at the corner points, particularly in inhomogeneous shapes, in comparing to smooth boundary.39
KHALIGHI ET AL. 19
FIGURE A1 The geometry of the plane body
discretization with N = 216,M = 741 with
DG-RBF technique and contour plots of absolute
errors when K = 16 for Example 6.6 [Colour
figure can be viewed at wileyonlinelibrary.com]
20 KHALIGHI ET AL.
TABLE A1 Notations and abbreviations used in this paper
Notation Description
ODE Ordinary differential equation
PDE Partial differential equation
BEM Boundary element method
TFPDE Time-fractional PDE
TFCDE Time-fractional convection–diffusion equation
COM Chebyshev operational matrix
RBF Radial basis function
DG-RBF Delaunay graph-based RBF method
Ω Two-dimensional domain
𝛤 Two-dimensional boundary
D𝛼𝑗c Caputo fractional time derivative of order 𝛼j
u(x, t) Unknown field function of spatial x ∈ Ω ∪ 𝛤 and time t
un Normal derivative of u
A,B,C, Given coefficient functions of x and j = 0, 1, … , k
D,E,F, 𝛾 j
g Given independent function of x and t
B Linear operator with respect to x of order one
h(x, t) Given function in the boundary condition; x ∈ 𝛤
di(x) Given function in the initial condition; i = 0, 1, … ,m − 1
𝔅(x, t) Unknown fictitious source function
u*, u∗n Fundamental solution of Equation (5), and normal derivative of u* on the boundary
û, ûn Particular solution of Equation (9), and normal derivative of û
𝜀 Free term coefficient; 𝜀 = 1 if x ∈ Ω, 𝜀 = 𝜃∕2𝜋 if x ∈ 𝛤 , else 𝜀 = 0, see details in book39
𝜃 Interior angle between the tangents of boundary at point x, see details in book39
r Distance between two points (or mean of all distances of internal points)
M Number of interior points after discretization
xIP𝑗 M internal nodal points; j = 1, … ,M
N Number of boundary nodal points after discretization
xBP𝑗 N boundary nodal points; j = 1, … ,N
fj(r) Radial basis approximating functions, j = 1, 2, … ,M
Continued
KHALIGHI ET AL. 21
TABLE A1 Continues
Notation Description
H, G N × N known coefficient matrices from the integration of the kernel functions on the boundary
„A N ×M known coefficient matrix from the integration of the kernel function on the domain
u,un Unknown vectors of the nodal values of the boundary displacements and their normal derivatives
b(t) Vector of the nodal values of the fictitious source at theM internal nodal points
𝛿1, 𝛿2 N × N known diagonal matrices
h(t) Known boundary vector of h(xBP𝑗 , t), j = 1, … ,N
ûpq Vector of values for u and its derivatives at theM internal nodal points; û00 = û
Ĥpq, Ĝpq M × N known coefficient matrices from the integration of the kernel functions on the boundary
Âpq M ×M known coefficient matrix from the integration of the kernel functions on the domain
A,B,C, M ×M known diagonal matrices including the nodal values of the corresponding
functions A(x), B(x), C(x), D(x), E(x), F(x), 𝛾 j(x)
D,E,F, 𝛾 j
g(t) Known internal vector of g(xIP𝑗 , t), 𝑗 = 1, … ,M
TL,i(t) Shifted Chebyshev polynomials of degree i on the interval t ∈ (0,L)
𝔇(𝜐) (K + 1) × (K + 1) Chebychev operational matrix of derivatives of order 𝜐 in the sense of Caputo
Ψ M × (K + 1) unknown matrix
𝝎 M × (K + 1) known matrix
R(t) Residual vector of Equation (16) with lengthM
L∞ Maximum error
MRE Maximum relative error
RMS Root mean square
Pr order Convergence order of BEM approximated by the Lemma 5.1
P𝜏 order Convergence order of COM approximated by the Lemma 5.2
Cond(A) Condition number of matrix A of the algebraic system Equations 29