Computational and Applied Mathematics (2020) 39:259

https://doi.org/10.1007/s40314-020-01314-4

The novel operational matrices based on 2D-Genocchi

polynomials: solving a general class of variable-order
fractional partial integro-differential equations

Haniye Dehestani1 · Yadollah Ordokhani1 · Mohsen Razzaghi2

Received: 5 April 2020 / Revised: 19 August 2020 / Accepted: 21 August 2020

© SBMAC - Sociedade Brasileira de Matemática Aplicada e Computacional 2020

Abstract
The main purpose of this study was to introduce an efficient approximate approach for
solving a general class of variable-order fractional partial integro-differential equations (VO-
FPIDEs). First, the Genocchi polynomials properties and the pseudo-operational matrix of
the VO-fractional derivative and fractional integration are presented. Then, to approximate
an integral part of the problems, we obtain the dual pseudo-operational matrix of fractional
order with a new technique. The pseudo-operational matrices of fractional order and Genocchi
collocation method are applied to reduce the VO-FPIDEs to a system of algebraic equations.
The error estimation indicates that the approximate solution converges to the exact solution.
Also, we discuss in detail on the upper bound of error for the pseudo-operational matrix of
fractional integration. In addition, to illustrate the efficiency and applicability of the approach,
we present several numerical examples.

Keywords Genocchi polynomials · Fractional pseudo-operational matrix · Variable-order

fractional partial integro-differential equations · Mixed Riemann–Liouville integral ·
Variable-order fractional derivative

Mathematics Subject Classification 26A33 · 33F05 · 35R09

Communicated by Agnieszka Malinowska.

B Yadollah Ordokhani
[email protected]
Haniye Dehestani
[email protected]
Mohsen Razzaghi
[email protected]
1 Department of Mathematics, Faculty of Mathematical Sciences, Alzahra University, Tehran, Iran
2 Department of Mathematics and Statistics, Mississippi State University, Starkville, MS 39762, USA

0123456789().: V,-vol 123

259 Page 2 of 32 H. Dehestani et al.

1 Introduction

Partial integro-differential equations play an important role in various fields of engineer-

ing and science such as plasma physics (Meleshko et al. 2010), electromagnetic theory
(Bloom 1981), nuclear reactor dynamics (Pachpatte 1983; Pao et al. 1979; Pao 1974), geo-
physics (Grasselli et al. 1990), reaction diffusion problems (Engler 1983), poro-viscoelastic
media (Habetler and Schiffman 1970). Many researchers are investigating various forms
of partial integro-differential equations and constructed numerous analytical and numeri-
cal methods. Express some of these methods: Spline collocation methods (Greenwell-Yanik
and Fairweather 1986; Yan and Fairweather 1992), discontinuous Galerkin method (Larsson
et al. 1998), finite element methods (Yanik and Fairweather 1988; Ma 2007), radial basis
function-based differential quadrature method (Shu et al. 2004), Homotopy perturbation
method (Babolian and Dastani 2011), differential transform method (Jang 2009), and there
are other methods that readers can refer to (Aghazadeh and Khajehnasiri 2013; Fakhar-Izadi
and Dehghan 2011; Avazzadeh et al. 2012; Tari et al. 2009; Khajehnasiri 2016).
Fractional calculus has application in the wide area of fluid mechanics, biology, control
theory of dynamical systems, probability and statistics, viscoelasticity, polymer modeling,
finance, physics and engineering (Gorenflo and Francesco 1997; Magin et al. 2013; Beyer
and Kempfle 1995; Torvik and Bagley 1984; Kulish and Lage 2002).
Recently, the concepts of the fractional calculus of variable order have received con-
siderable attention, which is created from constant-order fractional calculus. Variable-order
fractional derivatives are employed in several physical branches, for more information about
this subject refer to (Coimbra 2003; Ingman and Suzdalnitsky 2004; Ramirez and Coimbra
2007; Ostalczyk and Rybicki 2008; Sun et al. 2011; Diaz and Coimbra 2009). Various defini-
tions for the VO operators exist such as the Riemann–Liouville, Marchaud, Coimbra, Caputo
and Grunwald (Soon et al. 2005; Samko 1995).
The VO fractional problems have been investigated in recent years. For example, Li
and Wu (2015) introduced a numerical technique for variable fractional functional bound-
ary value problems. Heydari (2016) employed Chebyshev wavelets for the variable-order
time fractional mobile–immobile advection–dispersion model. Moghaddam and Machado
(2017) considered spline finite difference scheme for solving nonlinear time variable order
fractional partial differential equations. Yaghoobi et al. (2017) discussed an efficient cubic
Spline approximation for VO fractional differential equations with time delay. For more
information on this topic, see (Babaei et al. 2020; Ganji et al. 2020; Ganji and Jafari 2020).
For the first time, Genocchi numbers and polynomials were introduced between 1817
and 1889 by Angelo Genocchi. Recently, many researchers have applied them in various
fields of Mathematics (Araci et al. 2013; Isah and Phang 2016). Also, researchers have
worked on the Genocchi polynomials for the solutions of variable-order time-fractional
partial differential equations and various kinds of fractional partial differential equations
(Dehestani et al. 2019b, c). Furthermore, (Sadeghi et al. 2020) applied operational matrix
for Atangana–Baleanu derivative based on Genocchi polynomials for solving fractional dif-
ferential equations. The fractional Genocchi functions have been applied to the numerical
solution of variable-order fractional partial integro-differential equations (Dehestani et al.
2020a).

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 3 of 32 259

Consider the following general class of nonlinear variable-order fractional partial integro-
differential equations:

γ (x,t) ν(x,t)

2
∂ i u(x, t) 
2
∂ j u(x, t)
Dx u(x, t) + Dt u(x, t) + μi (x, t) + ϑ j (x, t)
∂xi ∂t j
i=1 j=1
+u(x, t) = g(x, t) + λIθr (k(x, t, ξ, η)u(ξ, η)), (x, t) ∈ Ω = [0, 1] × [0, 1],
(1)

0 < γ (x, t), ν(x, t) ≤ 1, r = (r1 , r2 ) ≥ 0, θ = (0, 0),

with the supplementary conditions

1
∂i u 
1
∂ ju
(x, 0) + (x, 0) = αk (x), k = 0, 1, 0 ≤ x ≤ 1,
∂x i ∂t j
i=0 j=0


1
∂i u 
1
∂ ju
(0, t) + (0, t) = βl (t), l = 0, 1, 0 ≤ t ≤ 1, (2)
∂x i ∂t j
i=0 j=0

where u(x, t) is an unknown function. The known functions g(x, t), k(x, t, ξ, η), μi (x, t),
i = 1, 2, ϑ j (x, t), j = 1, 2, αk (x), l = 0, 1 and βl (t), l = 0, 1 are defined on interval
γ (x,t) ν(x,t)
Ω = [0, 1] × [0, 1] and λ is a real or complex constant. And also Dx and Dt are
the variable-order fractional derivative operator with respect to space and time, respectively,
and Iθr is a mixed Riemann–Liouville integral operator.
In this work, we introduce new matrices, called the pseudo-operational matrices, which
are more accurate than the operational matrices. Then, by applying these matrices and
approximating the various terms in the underlying VO-FPIDEs by linear combinations of the
Genocchi polynomials, we convert the problem to a system of algebraic equations. Finally,
to demonstrate the applicability and validity of the novel technique we examine several
examples.
The present article is organized as follows: In Sect. 2, we briefly describe the definitions
and properties of the fractional calculus. Section 3 is devoted to the properties of Genoc-
chi polynomial and function approximation. In Sects. 4 and 5, we obtain the Genocchi
polynomial’s pseudo-operational matrices of variable-fractional order derivative, pseudo-
operational matrices of fractional order integration and dual pseudo-operational matrix of
fractional order for the solution of VO-FPIDEs. In Sect. 6, we present the error vector of the
pseudo-operational matrix of fractional integration. The numerical technique is presented in
Sect. 7. In Sect. 8, we report the numerical experiments to demonstrate the accuracy of the
proposed method. Section 9 concludes this article with a brief summary.

2 Preliminaries

In this section, we focus on some definitions and properties of the fractional calculus theory
(Dehestani et al. 2019a; Chen et al. 2014; Vityuk and Golushkov 2004; Najafalizadeh and
Ezzati 2016).

Definition 1 Let the fractional order be a continuous function of two variables γ : [0, κ]2 →
(q − 1, q) and u ∈ C q ([0, κ]2 ). The variable-order fractional derivative of u(x, t) in the

123
259 Page 4 of 32 H. Dehestani et al.

Caputo sense is defined as follows:

 t q−γ (x,t)−1 ∂ q u(x,s) ds,
Γ (q−γ (x,t)) 0 (t − s) q − 1 < γ (x, t) < q,
1
γ (x,t) ∂s q
D
0 t u(x, t) = ∂ q u(x,t)
∂t q , γ (x, t) = q ∈ z+ .
It has a following useful property:

Γ (β+1) β−γ (x,t) ,
γ (x,t) β Γ (β−γ (x,t)+1) t q ≤ β ∈ z+ ,
0 Dt t =
0, otherwise.

Definition 2 Let r = (r1 , r2 ) ∈ (0, ∞) × (0, ∞), θ = (0, 0) and u ∈ L 1 ([0, a] × [0, b]).
The left-sided mixed Riemann–Liouville integral of order r of u is defined by
 x t
1
(Iθr u)(x, t) = (x − ξ )r1 −1 (t − η)r2 −1 u(ξ, η)dηdξ.
Γ (r1 )Γ (r2 ) 0 0
It has a following useful property:
• (Iθθ u)(x, t) = u(x, t).
x t
• (Iθδ u)(x, t) = 0 0 u(ξ, η)dξ dη, ∀(x, t) ∈ [0, a] × [0, b], δ = (1, 1).
• If u ∈ L 1 ([0, a] × [0, b]), then Iθr exists for all r1 , r2 ∈ (0, ∞). Moreover, Iθr u ∈
C([0, a] × [0, b]) provided u ∈ C([0, a] × [0, b]) and
(Iθr u)(x, 0) = (Iθr u)(0, t) = 0.
• Let λ, ω ∈ (−1, ∞) and r1 , r2 ∈ (0, ∞). Then,
Γ (1 + λ)Γ (1 + ω)
Iθr x λ t ω = x λ+r1 t ω+r2 .
Γ (1 + λ + r1 )Γ (1 + λ + r2 )

3 Properties of Genocchi polynomials

3.1 Genocchi polynomials

The Genocchi polynomial G m (x) are defined by the exponential generating functions as
(Isah and Phang 2018; Dehestani et al. 2019a) follows:
∞

2te xt tm
= G m (x) , (|t| < π).
e +1
t m!
m=0

The Genocchi polynomials of degree m is defined as follows:

m  
m
G m (x) = gm−k x k , m = 1, 2, . . . , M, (3)
k
k=0

where gk = 2Bk − 2k+1 Bk , is called the Genocchi number and Bm is the Bernoulli number
(Dehestani et al. 2019a). We introduce some important properties of Genocchi polynomials:
• G m (1) + G m (0) = 0, m > 1.
m (x)
• dGdx = mG
⎧ m−1 (x), m ≥ 1.
⎨ 0,  m ≤ k,
• d Gdxmk(x) =
k
m k, m ∈ N ∪ 0.
⎩ k! G m−k (x), m > k,
k

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 5 of 32 259

1 2(−1)m m!n!
• G m (x)G n (x)dx = (m+n)! gm+n , m, n ≥ 1.
0b G m+1 (b)−G m+1 (a)
• G m (x)dx = .
ax m+1
(x)−gm+1
• 0 G m (t)dt = G m+1m+1 .

3.2 Function approximation

A function f defined over [0, 1] may be expressed in terms of Genocchi polynomials as

follows:

M
f (x) cm G m (x) = C T G(x);
m=1

the vector coefficient C is given by

C = D −1
f (x), G(x), D =
G(x), G(x),
where
C = [c1 , c2 , . . . , c M ]T , G(x) = [G 1 (x), G 2 (x), . . . , G M (x)]T ,
in which
., . denotes the inner product.
Theorem 1 Assume  Mthat f is a sufficiently smooth function approximated by the truncated
Genocchi series m=1 cm G m (x); then the coefficients cm for m = 1, 2, . . . , M, can be
calculated from the following relation (Isah and Phang 2018):
1
cm =( f (m−1) (1) + f (m−1) (0)). (4)
2m!
Assume that {G 11 (x, t), . . . , G 1N (x, t), . . . , G M1 (x, t), . . . , G M N (x, t)} is the set of 2D-
Genocchi polynomials. Let f (x, t) be an arbitrary element defined on Ω and
YG = span{G 11 (x, t), . . . , G 1N (x, t), . . . , G M1 (x, t), . . . , G M N (x, t)},
where
G mn (x, t) = G m (x)G n (t), m = 1, 2, . . . , M, n = 1, 2, . . . , N .
Since YG is a finite dimensional vector space, f (x, t) has the unique best approximation out
of SG as

M 
N
f (x, t) f mn G mn (x, t) = F T G(x, t), (5)
m=1 n=1
where
G(x, t) = [G 11 (x, t), G 12 (x, t), . . . , G 1N (x, t),
. . . , G M1 (x, t), G M2 (x, t), . . . , G M N (x, t)]T , (6)
and
F = [ f 11 , f 12 , . . . , f 1N , f 21 , f 22 , . . . , f 2N , . . . , f M1 , f M2 , . . . , f M N ]T .
Each element of F is represented as follows:

f , G m , G n 
f mn = , m = 1, 2, . . . , M, n = 1, 2, . . . , N ,

G m , G m 
G n , G n 

123
259 Page 6 of 32 H. Dehestani et al.

in which
., . denotes the inner product.

Theorem 2 Suppose that f (x, t) is a sufficiently smooth function in Ω and let


M 
N
f (x, t) = f mn G mn (x, t)
m=1 n=1

be its expansion in terms of 2D-Genocchi polynomials as described in Eq. (5). Then, the
coefficients f mn can be calculated as follows (Dehestani et al. 2020b; Loh and Phang 2018):
 m+n−2 
1 ∂ f (1, 1) ∂ m+n−2 f (1, 0) ∂ m+n−2 f (0, 0) ∂ m+n−2 f (0, 1)
f mn = + + + ,
4m!n! ∂ x m−1 ∂t n−1 ∂ x m−1 ∂t n−1 ∂ x m−1 ∂t n−1 ∂ x m−1 ∂t n−1
(7)

m = 1, 2, . . . , M, n = 1, 2, . . . , N .

3.3 Convergence

In this part, we will investigate the error analysis. We suppose that f (x, t) is a sufficiently
smooth function on Ω and p M N (x, t) is the interpolating polynomial to f at points (xi , t j );
then (Nemati and Ordokhani 2013)

∂ M+1 f (ξ, t)  ∂ N +1 f (x, η) 

M N
f (x, t) − p M N (x, t) = (x − x m ) + (t − tn )
∂ x M+1 (M + 1)! ∂t N +1 (N + 1)!
m=1 n=1

∂ M+N +2 f (ξ̄ , η̄) 

M 
N
− (x − xm ) (t − tn ),
∂ x M+1 ∂t N +1 (M+ 1)!(N + 1)!
m=1 n=1

where ξ, ξ̄ , η, η̄ ∈ [0, 1]. Therefore, due to the result obtained in Nemati and Ordokhani
(2013), we have
   
1 1 M+1 1 1 N +1
| f (x, t) − p M N (x, t)| ≤ C1 + C2
4 M 4 N
  M+1   N +1
1 1 1
+ C3 . (8)
16 M N
where
M+1 N +1
∂ f (x, t) ∂ f (x, t)
C1 = sup , C2 = sup ,
(x,t)∈Ω ∂ x M+1 (x,t)∈Ω ∂t N +1
M+N +2
∂ f (x, t)
C3 = sup .
(x,t)∈Ω ∂ x M+1 ∂t N +1

Theorem 3 Let f M N (x, t) be an approximation of f (x, t), which expand by the Genocchi
polynomials


M 
N
f M N (x, t) = f mn G mn (x, t) = F T G(x, t),
m=1 n=1

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 7 of 32 259

where G(x, t), m = 1, 2, . . . , M, n = 1, 2, . . . , N , are defined in Eq. (6). Then, we have

   
1 1 M+1 1 1 N +1
 f (x, t) − f M N (x, t)2 ≤ C1 + C2
4 M 4 N
  M+1   N +1
1 1 1
+ C3 .
16 M N

Proof By considering Eq. (8) and definition of the best approximation

 f − f M N 2 ≤  f − p M N 2 ,

we have

 f − f M N 2
 1  1  21
≤ | f (x, t) − p M N (x, t)| dxdt
2
0 0
 1
1
1 1 1 1 1 1 1 1
2 2
≤ ( ) M+1 C1 + ( ) N +1 C2 + ( ) M+1 ( ) N +1 C3 dxdt
0 0 4 M 4 N 16 M N
      M+1   N +1
1 1 M+1 1 1 N +1 1 1 1
= C1 + C2 + C3 .
4 M 4 N 16 M N

Hence, the result is achieved easily. 

4 Operational and pseudo-operational matrices of derivative

In this section, we derive the pseudo-operational matrix of derivative with integer and variable
fractional order for Genocchi polynomials. For this aim, using the properties of the Caputo
fractional derivative and Genocchi polynomials.

4.1 Integer operational matrices of derivative

The operational matrix of the derivative for the vector of 2D-Genocchi polynomials G(x, t)
with respect to x is obtained as follows:

Dx G(x, t) L 1 G(x, t),

Dx2 G(x, t) L 1 Dx G(x, t) L 21 G(x, t),
..
.
Dxk G(x, t) L k1 G(x, t), (9)

123
259 Page 8 of 32 H. Dehestani et al.

where
⎡ ⎤
O O O ··· O O O
⎢ 2I O O · · · O O O⎥
⎢ ⎥
⎢ O 3I O · · · O O O⎥
⎢ ⎥
⎢ O O 4I · · · O O O⎥
L1 = ⎢ ⎥,
⎢ .. .. .. . . .. .. .. ⎥
⎢ . . . . . . . ⎥
⎢ ⎥
⎣ O O O · · · (M − 1)I O O⎦
O O O ··· O MI O

I and O are the identity and zero matrix of dimensions N , respectively. Also, the operational
matrix of the derivative for the vector of 2D-Genocchi polynomials G(x, t) with respect to
t can be obtained as follows:

Dt G(x, t) L 2 G(x, t),

Dt2 G(x, t) L 2 Dt G(x, t) L 22 G(x, t),
..
.
Dtk G(x, t) L k2 G(x, t), (10)

where
⎡ ⎤
l O ··· O
⎢ O l ··· O ⎥
⎢ ⎥
L2 = ⎢ . . . ,
⎣ .. .. . . ... ⎥
⎦
O O ··· l

and Isah and Phang (2019)

⎡ ⎤
0 0 0 ··· 0 0 0
⎢2 0 0 ··· 0 0 0⎥
⎢ ⎥
⎢0 3 0 ··· 0 0 0⎥
⎢ ⎥
⎢ ··· 0⎥
l = ⎢0 0 4 0 0 ⎥.
⎢ .. .. .. .. .. .. .. ⎥
⎢. . . . . . .⎥
⎢ ⎥
⎣0 0 0 ··· N − 1 0 0⎦
0 0 0 ··· 0 N 0

4.2 Variable-order fractional pseudo-operational matrices of derivative

We obtain the pseudo-operational matrices of VO fractional derivative of order γ (x, t), ν(x, t)
with respect to two variables x and t for 2D-Genocchi polynomials.

Theorem 4 The VO fractional derivative of order q − 1 < γ (x, t) ≤ q in the Caputo sense
for Genocchi polynomials with respect to x is given by

γ (x,t)

M
q,γ (x,t)
Dx G m (x) x q−γ (x,t) sm,i,k G i (x). (11)
i=1

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 9 of 32 259

Proof By applying properties of the Caputo fractional derivative and Genocchi polynomials
defined in Eq. (3), we have
m  
γ (x,t) m γ (x,t) k
Dx G m (x) = gm−k Dx (x )
k
k=0
m  
m Γ (k + 1)
= gm−k x k−γ (x,t)
k Γ (k + 1 − γ (x, t))
k=q
m  
 m Γ (k + 1)
= x q−γ (x,t) gm−k x k−q
k Γ (k + 1 − γ (x, t))
k=q

m
γ (x,t)
= x q−γ (x,t) bm,k x k−q , (12)
k=q

where
 
γ (x,t) m Γ (k + 1)
bm,k = gm−k .
k Γ (k + 1 − γ (x, t))

By using Genocchi polynomials, we expand x k−q as follows:


M
q
x k−q ck,i G i (x).
i=1

As a result,
⎛ ⎞
γ (x,t)

m
γ (x,t)

M 
M m
q,γ (x,t)
Dx G m (x) x q−γ (x,t) bm,k ck,i G i (x) = x q−γ (x,t)
q ⎝ dm,k,i ⎠ G i (x)
k=q i=1 i=1 k=q


M
q,γ (x,t)
= x q−γ (x,t) sm,i G i (x), (13)
i=1

where

q,γ (x,t)

m
q,γ (x,t) q,γ (x,t) γ (x,t) q
sm,i = dm,k,i , dm,k,i = bm,k ck,i .
k=q

Therefore, we can rewrite the above equation as follows:

 
γ (x,t) q,γ (x,t) q,γ (x,t) q,γ (x,t)
Dx G m (x) x q−γ (x,t) sm,1 sm,2 · · · sm,M G(x).

According to the above theorem, the pseudo-operational matrix of the variable-order frac-
tional derivative of order γ (x, t) with respect to x is obtained as follows:
γ (x,t) γ (x,t)
Dx G(x) x q−γ (x,t) S1 G(x), q − 1 < γ (x, t) ≤ q, q ∈ z+ . (14)

Similarly, the pseudo-operational matrix of the variable-order fractional derivative of order

ν(x, t) with respect to t is gained as follows:
ν(x,t) ν(x,t)
Dt G(t) t p−ν(x,t) S2 G(t), p − 1 < ν(x, t) ≤ p, p ∈ z+ . (15)

123
259 Page 10 of 32 H. Dehestani et al.

Now, we calculate the pseudo-operational matrix of variable-order fractional derivative for

2D-Genocchi polynomials with respect to x as follows:
γ (x,t) γ (x,t)
Dx G(x, t) x q−γ (x,t) W1 G(x, t), (16)

in which each element of this matrix is computed as follows:

⎛ ⎞
γ (x,t) γ (x,t)

M m
(x,t)
Dx G mn (x, t) = Dx (G m (x)G n (t)) x q−γ (x,t) ⎝ d
q,γ ⎠ G i (x)G n (t)
m,k,i
i=1 k=q
⎛ ⎞

M 
m
q,γ (x,t)
= x q−γ (x,t) ⎝ dm,k,i ⎠ G in (x, t)
i=1 k=q


M
q,γ (x,t)
= x q−γ (x,t) sm,i G in (x, t).
i=1

In other words,
⎡ ⎤
O ··· O O O
⎢ s q,γ (x,t) × I ··· O O O⎥
⎢ 2,1 ⎥
⎢ q,γ (x,t) ⎥
W1
γ (x,t)
=⎢ O s3,2 × I ··· O O ⎥.
⎢ ⎥
⎢ .. .. .. .. .. ⎥
⎣ . . . . . ⎦
q,γ (x,t)
O O · · · s M,M−1 × I O

where O and I are the zero and identity matrix with dimension N , respectively. Similarly,
the pseudo-operational matrix of variable-order fractional derivative for 2D-Genocchi poly-
nomials with respect to t is achieved as follows:
ν(x,t) ν(x,t)
Dt G(x, t) t p−ν(x,t) W2 G(x, t), (17)

in which each element is

ν(x,t) ν(x,t)

N
p,ν(x,t)
Dt G mn (x, t) = Dt (G m (x)G n (t)) t p−ν(x,t) sm, j G m (x)G j (t)
j=1


N
p,ν(x,t)
= t p−ν(x,t) sm, j G m j (x, t). (18)
j=1

Thus, the pseudo-operational matrix can be written as follows:

⎡ ν(x,t)
⎤
S2 O ··· O
⎢ ν(x,t) ⎥
ν(x,t) ⎢ O S2 ··· O ⎥
W2 =⎢
⎢ .. .. .. .. ⎥ ⎥.
⎣ . . . . ⎦
ν(x,t)
O O · · · S2

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 11 of 32 259

5 Pseudo-operational matrices of integration

The purpose of the present section was to introduce the pseudo-operational matrices of inte-
gration and dual pseudo-operational matrix of the fractional-order for Genocchi polynomials.

5.1 Pseudo-operational matrices of fractional integration

The main objective of this section was to introduce the pseudo-operational matrix of the
fractional-order integration for 2D-Genocchi polynomials by using Riemann-Liouville inte-
gral properties.

Theorem 5 The fractional integration of order r1 > 0 for Genocchi polynomials is given by

M
r1
I xr1 G m (x) x r1 ξm,k, j G j (x). (19)
j=1

Proof In order to obtain the desired result, we use Riemann–Liouville integral properties and
Eq. (3), so that
 m    m  
 m  m
I x G m (x) = I x
r1 r1
gm−k x =
k
gm−k I xr1 (x k )
k k
k=0 k=0
 m   m
m Γ (k + 1) r1
= gm−k x k+r1 = x r1 ρm,k xk, m = 1, 2, . . . , M,
k Γ (k + r1 + 1)
k=0 k=0
(20)
where
 
r1 m Γ (k + 1)
ρm,k = gm−k .
k Γ (k + r1 + 1)
Suppose that the approximation of x k by Genocchi polynomials is as follows:

M
r1
xk ωk, j G j (x). (21)
j=1

By replacing the above equation in Eq. (20), we get

 m 

m
r1

M
r1

M  r1 r1
I xr1 G m (x) x r1
ρm,k ωk, j G j (x) =x r1
ρm,k ωk, j G j (x)
k=0 j=1 j=1 k=0


M
r1 r1

m
r1 r1
= x r1 ξm,k, j G j (x), ξm,k, j = ρm,k ωk, j. (22)
j=1 k=0

In other words, from the above theorem the following can be written:
 m 
 r1 r1  m
r1 r1 
m
r1 r1
I x G m (x) x
r1 r1 ρm,k ωk,1 ρm,k ωk,2 · · · ρm,k ωk,M G(x),
k=0 k=0 k=0
m = 1, 2, . . . , M.
(23)

123
259 Page 12 of 32 H. Dehestani et al.

As a result, the pseudo-operational matrices of fractional integration is gained as follows:

I xr1 G(x) x r1 Υ1r1 G(x), (24)

where
⎡ ⎤

1
r1 r1 
1
r1 r1 1 r1 r1
⎢ ρ1,k ωk,1 ρ1,k ωk,2 ··· ρ ω
k=0 1,k k,M ⎥
⎢ k=0 k=0 ⎥
⎢ ⎥
⎢ ⎥
⎢   2 ⎥
⎢ 2 ρ r 1 ωr 1 2
r1 r1
ρ2,k ωk,2 ··· r1 r1 ⎥
⎢ k=0 ρ2,k ωk,M ⎥
Υ1r1 ⎢
= ⎢ k=0
2,k k,1
⎥.
k=0 ⎥
⎢ ⎥
⎢ .. .. .. ⎥
⎢ .. ⎥
⎢ . . . . ⎥
⎢ ⎥
⎣M
r1 r1 
M
r1 r1 
M ⎦
ρ M,k ωk,1 ρ M,k ωk,2 ··· ρ rM,k
1 r1
ωk,M
k=0 k=0 k=0

Here, Υ1r1 is called the Genocchi polynomial pseudo-operational matrix of the fractional
integration of order r1 . According to the above relations, the pseudo-operational matrix of
the fractional integration of order r2 > 0 for Genocchi polynomials with respect to t is
defined as follows:

Itr2 G(t) t r2 Υ2r2 G(t). (25)

Next, the pseudo-operational matrix of fractional integration for 2D-Genocchi polynomi-

als of order r = (r1 , r2 ) > 0 is obtained as follows:
(r1 ,r2 )
Iθr G mn (x, t) = Iθ (G m (x)G n (t)) = (I xr1 G m (x))(Itr2 G n (t))

M
r1

N
r2
x r1 t r2 ξm,k,i G i (x) ξn,l, j G j (t)
i=1 j=1


M 
N
r1 r2
= x r1 t r2 ξm,k,i ξn,l, j G i (x)G j (t)
i=1 j=1


M 
N
= x r1 t r2 ηm,n,k,l,i,
r
j G i j (x, t), (26)
i=1 j=1

where
r1 r2
ηm,n,k,l,i,
r
j = ξm,k,i ξn,k, j , m = 1, 2, . . . , M, n = 1, 2, . . . , N .

Therefore, from the above equation, we conclude

Iθr G(x, t) x r1 t r2 Ψ r G(x, t), (27)

where

Ψ r = Υ1r1 ⊗ Υ2r2 , (28)

and Ψ r is called the pseudo-operational matrix of the fractional integration for 2D-Genocchi
polynomials of order r . It is worth to mention that the symbol ⊗ denotes the Kronecker
product (Schafer 1996). 

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 13 of 32 259

5.2 Dual pseudo-operational matrix of fractional order

In this section, we introduce the fractional dual pseudo-operational matrix of 2D-Genocchi

polynomials. First, we derive the fractional dual pseudo-operational matrix of order r1 for
Genocchi polynomials as follows:

I xr1 (G(x)G T (x)) x r1 Q r11 G(x). (29)

Using Eq. (3) and the properties of the operator I xr1 , we obtain each element of the fractional
dual pseudo-operational matrix as follows:

m   i  

m i
I xr1 (G m (x)G i (x)) = I xr1 ( gm−k x k gi− j x j )
k j
k=0 j=0

 i 
m   
m i
= gm−k gi− j I xr1 (x k+ j )
k j
k=0 j=0

 i   
m 
m i Γ (k + j + 1)
= gm−k gi− j x k+ j+r1
k j Γ (k + j + 1 + r1 )
k=0 j=0

 i   
m 
m i Γ (k + j + 1)
= x r1 gm−k gi− j x k+ j
k j Γ (k + j + 1 + r1 )
k=0 j=0


m 
i
= x r1 dmik j x k+ j , (30)
k=0 j=0

where
  
m i Γ (k + j + 1)
dmik j = gm−k gi− j .
k j Γ (k + j + 1 + r1 )

Let x k+ j be expanded in terms of Genocchi polynomials as


M
x k+ j ak js G s (x).
s=1

Now, by replacing the above equation in Eq. (30), we have


m 
i 
M 
M  m 
i
I xr1 (G m (x)G i (x)) x r1 dmik j ak js G s (x) = x r1 ( dmik j ak js )G s (x)
k=0 j=0 s=1 s=1 k=0 j=0


M
= x r1 qmis G s (x), (31)
s=1

where


m 
i
qmis = dmik j ak js .
k=0 j=0

123
259 Page 14 of 32 H. Dehestani et al.

Consequently, we get
⎡ ⎤
q111 q112 · · · q11M
⎢ q212 q212 · · · q21M ⎥
⎢ ⎥
Q r11 =⎢ . .. .. .. ⎥ .
⎣ .. . . . ⎦
q M M1 q M M2 · · · qM M M
Similarly, according to the above process, we have
Itr2 (G(t)G T (t)) t r2 Q r22 G(t). (32)
As a result, by using Eqs. (29) and (32), we deduce
Iθr (G(x, t)G T (x, t)) x r1 t r2 Q̂ r G(x, t) = x r1 t r2 (Q r11 ⊗ Q r21 )G(x, t), (33)

where Q̂ r is the fractional dual pseudo-operational matrix of 2D-Genocchi polynomials so

that each element of Q̂ r is obtained using the following formula:
(r1 ,r2 )
Iθr (Ĝ mn (x, t)Ĝ i j (x, t)) = Iθ (G m (x)G n (t)G i (x)G j (t))
= (I xr1 G m (x)G i (x))(Itr2 G n (t)G j (t))

M 
N
x r1 t r2 qmps qn jl G s (x)G l (t)
s=1 l=1

M 
N
x r1 t r2 qmps qn jl G sl (x, t). (34)
s=1 l=1

6 Error vector of the pseudo-operational matrix of fractional

integration

Lemma 1 Suppose that H is a Hilbert space and Y is a closed subspace of H such that
dim(Y ) < ∞ and {y1 , y2 , . . . , ym } is a basis for Y . Let z be an arbitrary element in H and
y ∗ be the unique best approximation to z out of Y . Then (Kreyszig 1978)
Gram(z, y1 , y2 , . . . , ym )
z − y ∗ 22 = , (35)
Gram(y1 , y2 , . . . , ym )
where


z, z
z, y1  · · ·
z, ym 


y1 , z
y1 , y1  · · ·
y1 , ym 

Gram(z, y1 , y2 , . . . , ym ) = . .. .. .. . (36)
.. . . .


ym , z
ym , y1  · · ·
ym , ym 

Lemma 2 Suppose φ ∈ L 2 [0, 1] is approximated by φ M as (Rivlin 1981)


M
φ(x) φ M (x) = τm G m (x).
m=1

Consider
L M (φ) = φ(x) − φ M (x)2 .

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 15 of 32 259

Then, we have

M → ∞, ⇒ L M (φ) → 0.

The error vector of the pseudo-operational matrix of fractional integration for Genocchi
polynomials with respect to variable x is defined as follows:

E xr1 = I xr1 G − Υ1r1 G, r1 > 0, (37)

where
 
E xr1 = erx,m
1 , m = 1, 2, . . . , M.

According to Eq. (21) and Lemma 1, we get


M   21
r1 Gram(x k , G 1 (x), . . . , G M (x))
x −
k
ωk, j G j (x)2 = . (38)
Gram(G 1 (x), G 2 (x), . . . , G M (x))
j=1

Due to the above equations and Eq. (22), we have


m
r1

M
r (x,t)
1  = I r1 G (x) − x r1
erx,m 2 x m ρm,k ωk,1 j G j (x)2
k=0 j=1
m  
  r M
m Γ (k + 1)
≤ |gm−k |x r1 2 x k − ωk,1 j G j (x)2
k Γ (k + r1 + 1)
k=0 j=1
m  
   21
m Γ (k + 1) Gram(x k , G 1 (x), . . . , G M (x))
≤ |gm−k |
k Γ (k + r1 + 1) Gram(G 1 (x), G 2 (x), . . . , G M (x))
k=0
= δm , m = 1, 2, . . . , M. (39)

By considering Eq. (39) and Lemma 2, we can deduce that by increasing the number of
Genocchi polynomials bases M, the error vector E xr1 tends to zero.
According to the above process, we introduce the error vector of the pseudo-operational
matrix of fractional integration for Genocchi polynomials with respect to variable t as follows:

E tr2 = Itr2 G − Υ2r2 G, r2 > 0, (40)

where
 r2 
E tr2 = et,n , n = 1, 2, . . . , N .

So, we have


N   21
r2 Gram(t k , G 1 (t), . . . , G N (t))
t −
k
ωk, j G j (x)2 = . (41)
Gram(G 1 (t), G 2 (t), . . . , G N (t))
j=1

123
259 Page 16 of 32 H. Dehestani et al.

As a result,

r2

n 
N
et,n 2 = Itr2 G n (t) − t r2 r2
ρn,k r2
ωk, j G j (x)2
k=0 j=1
n 
   r (x,t) N
n Γ (k + 1)
≤ |gn−k |t r2 2 t k − ωk,2 j G j (x)2
k Γ (k + r2 + 1)
k=0 j=1
n  
   21
n Γ (k + 1) Gram(t k , G 1 (t), . . . , G N (t))
≤ |gn−k |
k Γ (k + r2 + 1) Gram(G 1 (t), G 2 (t), . . . , G N (t))
k=0
= σ n, n = 1, 2, . . . , N . (42)
By considering Eq. (42) and Lemma 2, we can deduce that by increasing the number of
Genocchi polynomials bases N , the error vector E tr2 tends to zero.
Consequently, the error vector of the pseudo-operational matrix of fractional integration
for 2D-Genocchi polynomials is defined as follows:
Rθr = Iθr G − Ψ r G, r = (r1 , r2 ) > 0, (43)
where
 r 
Rθr = r̃ x,t,m,n , m = 1, 2, . . . , M, n = 1, 2, . . . , N .
So, by using Eqs. (26), (39), (42) and Lemma 2, we have


M 
N
r̃ x,t,m,n
r
2 = Iθr G mn (x, t) − ηm,n,k,l,i,
r
j G i j (x, t)
i=1 j=1 2
M ⎛ N ⎞
 
= Iθr G m (x)G n (t) − r1
ξm,k,i G i (x) ⎝ r2
ξn,l, j G j (t)
⎠
i=1 2 j=1
 ⎛ N ⎞

M 
= I xr1 G m (x)Itr2 G n (t) − r1
ξm,k,i G i (x) ⎝ r2
ξn,l, j G j (t)
⎠
i=1 j=1 2
⎛ ⎞

N
≤ I xr1 G m (x) ⎝ Itr2 G n (t) − r2
ξn,l, j G j (t)
⎠
j=1 2
 

M
r1

N
r2
+ I xr1 G m (t) − ξm,k,i G i (x) ξn,l, j G j (t)
i=1 j=1 2


N
≤ I xr1 G m (x)2 Itr2 G n (t) − r2
ξn,l, j G j (t)
j=1 2


M
r1

N
r2
+ I xr1 G m (t) − ξm,k,i G i (x)2  ξn,l, j G j (t)
i=1 j=1 2


N
r2
≤ σ n I xr1 G m (x)2 + δ m ξn,l, j G j (t) . (44)
j=1 2

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 17 of 32 259

Now, by using Eq. (3), we find that

⎛ 2 ⎞ 21
 m  
1 
m Γ (k + 1) k+r1
I xr1 G m (x)2 =⎝ gm−k x dx ⎠
0 k Γ (k + r1 + 1)
k=0
! m  2
"! m 
" 1
 m Γ (k + 1)2  1 2
≤ 2
gm−k x 2k+2r1
dx
k Γ (k + r1 + 1)2
k=0 k=0 0
 m   1
 m 2 Γ (k + 1)2 2
= 2
gm−k , (45)
k Γ (k + r1 + 1)2 (2k + 2r1 + 1)
k=0

and
⎛ 2 ⎞ 21
 N

N
⎜
1  ⎟
 r2
ξn,l, r2 dt ⎠
j G j (t)2 = ⎝ ξ G
n,l, j j (t)
j=1 0 j=1
⎛⎡ ⎤⎡ ⎤⎞ 1
N N  1
2

≤ ⎝⎣ r
|ξn,l,
2
j | 2 ⎦⎣ |G j (t)| dt ⎦⎠ .
2
(46)
j=1 j=1 0

On the other hand,

 1 j  2 2
 j g j−k
|G j (t)|2 dt ≤ . (47)
0 k 2k + 1
k=1

Then,
⎛⎡ ⎤⎡ ⎤⎞ 1

N 
N  j  2 2
N  2
j g j−k
r2
ξn,l, j G j (x) ≤ ⎝⎣ r2
|ξn,l, j|
2⎦ ⎣ ⎦⎠ . (48)
k 2k + 1
j=1 2 j=1 j=1 k=1

Thus, we achieve
 m  2
1
 m Γ (k + 1)2 2
r̃ x,t,m,n
r
2 ≤ 2
gm−k σn
k Γ (k + r1 + 1)2 (2k + 2r1 + 1)
k=0
⎛⎡ ⎤⎡ ⎤⎞ 1
N  j  2 2
N  2
(x,t) j g
+δ m ⎝ ⎣ r
|ξn,l, j |
2 2⎦ ⎣ j−k ⎦ ⎠ . (49)
k 2k + 1
j=1 j=1 k=1

By considering the above descriptions, Lemmas 1 and 2, we can conclude that by increasing
the terms of the 2D-Genocchi polynomials the error vector Rθr tends to zero.

7 Genocchi collocation method

In this section, we describe a numerical scheme for obtaining the approximate solution of
VO-FPIDEs by using 2D-Genocchi polynomials. For solving the problem in Eq. (1), we

123
259 Page 18 of 32 H. Dehestani et al.

expand
u(x, t) U T G(x, t), (50)
where
U = [u 11 , u 12 , . . . , u 1N , u 21 , u 22 , . . . , u 2N , . . . , u M1 , u M2 , . . . , u M N ]T .
By using the operational and pseudo-operational matrixes of derivative defined in Sect. 4,
we have
γ (x,t) γ (x,t) γ (x,t)
Dx u(x, t) U T Dx G(x, t) x q−γ (x,t) U T W1 G(x, t),
ν(x,t) ν(x,t) ν(x,t)
Dt u(x, t) U T
Dt G(x, t) t p−ν(x,t)
U T
W2 G(x, t), (51)
and
∂ i u(x, t) ∂ j u(x, t) j
U T L i1 G(x, t), U T L 2 G(x, t). (52)
∂x i ∂t j
Also, the kernel function k(x, t, ξ, η) can be approximated by the 2D-Genocchi polynomials
as follows:
k(x, t, ξ, η) G T (x, t)K G G(ξ, η). (53)
The block matrix K G is defined as follows:
 M
K G = K ∗ i,m=1 ,
where
 N
K ∗ = K i jmn j,n=1 , m, i = 1, 2, . . . , M.
Each component of this matrix is obtained by the following formula:

k(x, t, ξ, η), G mn (x, t), G i j (ξ, η)
K i jmn = ,
G mn (x, t)2 G i j (x, t)2
m, i = 1, 2, . . . , M, n, j = 1, 2, . . . , N .
Thus, we approximate an integral part of Eq. (1) by using dual pseudo-operational matrix
and Eqs. (50) and (53) as follows:
Iθr (k(x, t, ξ, η)u(ξ, η)) Iθr (G T (x, t)K G G(ξ, η)G T (ξ, η)U )
= G T (x, t)K G Iθr (G(ξ, η)G T (ξ, η))U
x r1 t r2 G T (x, t)K G Q̂ r G(x, t)U . (54)
Consequently, by substituting Eqs. (50)–(52) and (54) in Eq. (1), we obtain
γ (x,t) ν(x,t)
x q−γ (x,t) U T W1 G(x, t) + t p−ν(x,t) U T W2 G(x, t)

2 
2
j
+ μi (x, t)U T L i1 G(x, t) + ϑ j (x, t)U T L 2 G(x, t)
i=1 j=1

+U G(x, t) − g(x, t) − λx t G T (x, t)K G Q̂ r G(x, t)U = 0.

T r1 r2
(55)
Substituting nodal points of Newton–Cotes (Dehestani et al. 2018),
2s − 1 2h − 1
xs = , s = 1, 2, . . . , M, th = , h = 1, 2, . . . , N , (56)
2M 2N

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 19 of 32 259

in Eq. (55), we get a system of algebraic equations

q−γ (xs ,th ) γ (xs ,th ) p−ν(xs ,th ) ν(xs ,th )
xs U T W1 G(xs , th ) + th U T W2 G(xs , th )

2 
2
j
+ μi (xs , th )U T L i1 G(xs , th ) + ϑ j (xs , th )U T L 2 G(xs , th )
i=1 j=1

+U G(xs , th ) −
T
g(xs , th ) − λxsr1 thr2 G T (xs , th )K G Q̂ r G(xs , th )U = 0. (57)
Using Eqs. (2), (52) and nodal points defined in Eq. (56), we obtain the approximate of
supplementary conditions

1 
1
j
U T L i1 G(xs , 0) + U T L 2 G(xs , 0) = αk (xs ), k = 0, 1, s = 1, 2, . . . , M,
i=0 j=0


1 
1
j
U T L i1 G(0, th ) + U T L 2 G(0, th ) = βl (th ), l = 0, 1, h = 1, 2, . . . , N . (58)
i=0 j=0

Then, by substituting the systems defined in Eq. (58) in last rows of the system in Eq. (57),
we achieve the ultimate system. Finally, we solve the ultimate system to achieve the unknown
vector U .

8 Illustrative examples

In this section, we consider eight numerical examples to illustrate the applicability and
accuracy of the proposed method.

Example 1 For the first example, consider the time VO fractional mobile–immobile
advection–dispersion model is as follows (Zhang et al. 2013; Jiang and Liu 2017):
∂u(x, t) ν(x,t) ∂u(x, t) ∂ 2 u(x, t)
+ Dt u(x, t) = − + + g(x, t),
∂t ∂x ∂x2
with initial condition u(x, 0) = 10x 2 (1 − x)2 , 0 ≤ x ≤ 1 and boundary conditions
u(0, t) = u(1, t) = 0, 0 ≤ t ≤ 1,
with
t 1−ν(x,t)
g(x, t) = 10(1 + )x 2 (1 − x)2
Γ (2 − ν(x, t))
+10(4x 3 − 6x 2 + 2x − 12x 2 + 12x − 2)(t + 1),
and
ν(x, t) = 1 − 0.5 exp(−xt).
The exact solution of this example is u(x, t) = 10(t + 1)x 2 (1 − x)2 . The numerical results of
this example are shown in Table 1. Table 1 displays the absolute errors obtained by present
method for M = 5, N = 2 and other methods in Zhang et al. (2013) and Jiang and Liu
(2017) at t = 1. From the comparison with the methods in Zhang et al. (2013) and Jiang and
Liu (2017), it is illustrated that our method is more accurate than methods in Zhang et al.
(2013) and Jiang and Liu (2017).

123
 259 Page 20 of 32 H. Dehestani et al.

Table 1 Comparison of absolute errors with Refs. (Zhang et al. 2013; Jiang and Liu 2017) at t = 1 of
Example 1
xi Present method Method in Zhang et al. (2013) Method in Jiang and Liu (2017)
M = 5, N = 2 n = 100 n = 20

0.1 3.02 × 10−12 1.56 × 10−4 6.37 × 10−5

0.2 2.39 × 10−12 1.40 × 10−3 4.90 × 10−6
0.3 1.91 × 10−12 2.97 × 10−3 5.98 × 10−5
0.4 1.55 × 10−12 4.29 × 10−3 7.78 × 10−6
0.5 1.31 × 10−12 4.97 × 10−3 5.80 × 10−5
0.6 1.16 × 10−12 4.80 × 10−3 8.29 × 10−6
0.7 1.09 × 10−12 3.81 × 10−3 5.89 × 10−5
0.8 1.09 × 10−12 2.27 × 10−3 5.94 × 10−6
0.9 1.15 × 10−12 7.20 × 10−4 6.29 × 10−5

Example 2 We consider the VO fractional partial integro-differential equation in the following

form:
γ (x,t) ν(x,t)
Dx u(x, t) + Dt u(x, t)
% %
x t 1 1 (1,1)
= 4t + 4x − x 4 t 3 − x 3 t 2 + Iθ ((x 2 t + ξ )u(ξ, η)),
π π 2 3
(x, t) ∈ [0, 1] × [0, 1], 0 < γ (x, t), ν(x, t) ≤ 1,
with the supplementary conditions
u(x, 0) = 0, 0 ≤ x ≤ 1,
and
u(0, t) = 0, 0 ≤ t ≤ 1.
The exact solution of this example, when γ (x, t) = ν(x, t) = 0.5 is u(x, t) = 2xt. For
M = N = 3 and γ (x, t) = ν(x, t) = 0.5, we obtain
u 1 = 0.4999999999999998, u 2 = 0.4999999999999998,
u 3 = −4.769467547464007 × 10−16 , u 4 = 0.5,
u 5 = 0.5, u 6 = −9.507055209956305 × 10−17 ,
u 7 = −5.244209351627797 × 10−16 , u 8 = −5.244209351627797 × 10−16
u 9 = −8.94445347252264 × 10−16 .
Then, the approximate solution is gained as follows:
u(x, t) = −8.0500081252703761 × 10−15 x 2 t 2 + 7.47958481267299746062 × 10−15 xt 2
−1.14562860794051303921 × 10−15 t 2 + 4.90348251429369821661 × 10−15 t x 2
+2xt + 6.66133814775093924254 × 10−16 t.

Also, the absolute error gained for various choices of γ (x, t) = ν(x, t) with M = N = 3 is
listed in Table 2. From the results of Table 2, it can be understood that as γ (x, t) and ν(x, t)
approach 0.5, the absolute error tends to zero.

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 21 of 32 259

Table 2 The absolute errors gained by the present method for various choices of γ (x, t) = ν(x, t) with
M = N = 3 of Example 2
xi = ti γ (x, t) = ν(x, t)
0.5 − 0.05xt 0.5 − 0.03xt 0.5 − 0.01xt 0.5

0 3.1598 × 10−3 1.9142 × 10−3 6.4406 × 10−4 1.1102 × 10−15

0.1 2.3864 × 10−4 1.5511 × 10−4 5.5607 × 10−5 3.1907 × 10−16
0.2 2.8139 × 10−3 1.6756 × 10−3 5.5431 × 10−4 3.9419 × 10−15
0.3 5.9979 × 10−3 3.5778 × 10−3 1.1857 × 10−3 1.0296 × 10−15
0.4 9.3134 × 10−3 5.5517 × 10−3 1.8385 × 10−3 1.5870 × 10−15
0.5 1.2760 × 10−2 7.5971 × 10−3 2.5129 × 10−3 2.0666 × 10−15
0.6 1.6339 × 10−2 9.7141 × 10−3 3.2086 × 10−3 2.4682 × 10−15
0.7 2.0049 × 10−2 1.1903 × 10−2 3.9259 × 10−3 2.7920 × 10−15
0.8 2.3890 × 10−2 1.4163 × 10−2 4.6646 × 10−3 3.0379 × 10−15
0.9 2.7863 × 10−2 1.6495 × 10−2 5.4248 × 10−3 3.2058 × 10−15
1 3.1967 × 10−2 1.8898 × 10−2 6.2064 × 10−3 3.2959 × 10−15

Example 3 We consider the VO fractional partial integro-differential equation in the following

form:
ν(x,t)
Dt u(x, t) + u(x, t) = g(x, t) + Iθr ((x + ξ η)u(ξ, η)), (x, t) ∈ [0, 1] × [0, 1],
r = (r1 , r2 ) > 0, 0 < ν(x, t) ≤ 1,
with initial condition u(x, 0) = 0, 0 ≤ x ≤ 1, and
2x 2x 2+r1 t 2+r2 12x 2+r1 t 3+r2
g(x, t) = xt 2 + t 2−ν(x,t) − − .
Γ (2 − ν(x, t)) Γ (2 + r1 )Γ (3 + r2 ) Γ (3 + r1 )Γ (4 + r2 )
The exact solution of this example is u(x, t) = xt 2 . The numerical results of this example
are displayed in Table 3. Table 3 shows the absolute errors for various values of r1 , r2 and
various functions of ν(x, t) by applying the present method with M = 2, N = 3. From these
results, it is seen that the numerical method can be considered as an efficient method.

Example 4 We consider the partial integro-differential equation in the following form:

∂u ∂u
(x, t) + (x, t) = g(x, t) + Iθr (u(ξ, η)), r = (r1 , r2 ), (x, t) ∈ [0, 1] × [0, 1],
∂x ∂t
with initial condition u(x, 0) = exp(x), 0 ≤ x ≤ 1, and
g(x, t) = −1 + exp(x) + exp(t) + exp(x + t).
The exact solution of this example, when r1 , r2 = 1 is u(x, t) = exp(x + t). The numerical
results of this example are displayed in Table 4 and Figs. 1 and 2. Table 4 shows the absolute
errors for various values of M, N with r1 , r2 = 1. In Table 4, it is clear that by increasing the
number of Genocchi polynomials the absolute error tends to zero. The graphs of approximate
solutions for various values of r1 , r2 by applying the present method with M = N = 5 are
plotted in Fig. 1. From Fig. 1, it is seen that as r1 , r2 approaches 1, the numerical solutions
converge to the exact solution. Also, the absolute errors and approximate solution for M =
7, N = 5 with r1 = r2 = 1 are plotted in Fig. 2.

123
 259

123
Page 22 of 32

Table 3 The absolute error for different values of ν(x, t), r1 , r2 with M = 2, N = 3 of Example 3
(xi , ti ) r1 = 1 r1 = 0.25 r1 = 0.1 r1 = 100 r1 = 1.5 r1 = 20
r2 = 1 r2 = 0.75 r2 = 0.9 r2 = 1000 r2 = 2.5 r2 = 200
ν(x, t) = 1 ν(x, t) = 0.5 ν(x, t) = 0.2 ν(x, t) = 0.7 ν(x, t) = sin(x + t) ν(x, t) = 1 − exp(− xt
2)

(0.1, 0.1) 1.55 × 10−17 3.33 × 10−18 0 1.11 × 10−17 1.66 × 10−17 7.77 × 10−18
(0.2, 0.2) 1.77 × 10−17 1.33 × 10−17 2.86 × 10−42 2.22 × 10−17 1.11 × 10−17 8.88 × 10−18
(0.3, 0.3) 6.66 × 10−18 2.99 × 10−17 0 3.33 × 10−17 1.66 × 10−17 3.33 × 10−18
(0.4, 0.4) 1.77 × 10−17 5.32 × 10−17 0 4.44 × 10−17 6.66 × 10−17 8.88 × 10−18
(0.5, 0.5) 5.55 × 10−17 8.32 × 10−17 0 5.55 × 10−17 1.38 × 10−16 2.77 × 10−17
(0.6, 0.6) 1.06 × 10−16 1.19 × 10−16 4.59 × 10−41 6.66 × 10−17 2.33 × 10−16 5.32 × 10−17
(0.7, 0.7) 1.70 × 10−16 1.63 × 10−16 0 7.77 × 10−17 3.49 × 10−16 8.54 × 10−17
(0.8, 0.8) 2.48 × 10−16 2.13 × 10−16 1.83 × 10−40 8.88 × 10−17 4.88 × 10−16 1.24 × 10−16
(0.9, 0.9) 3.39 × 10−16 2.69 × 10−16 1.84 × 10−40 9.99 × 10−17 6.49 × 10−16 1.69 × 10−16
(1, 1) 4.44 × 10−16 3.33 × 10−16 0 1.11 × 10−16 8.32 × 10−16 2.22 × 10−16
H. Dehestani et al.
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 23 of 32 259

Table 4 The absolute error for (xi , ti ) M = 3, N = 3 M = 3, N = 5 M = 5, N = 5

different values of M, N with
r1 , r2 = 1 of Example 4 (0, 0) 6.51 × 10−3 6.49 × 10−3 3.79 × 10−5
(0.1, 0.1) 4.45 × 10−3 7.17 × 10−3 1.30 × 10−5
(0.2, 0.2) 3.25 × 10−3 7.08 × 10−3 2.42 × 10−5
(0.3, 0.3) 5.77 × 10−3 6.80 × 10−3 2.58 × 10−5
(0.4, 0.4) 6.73 × 10−3 6.87 × 10−3 2.32 × 10−5
(0.5, 0.5) 2.14 × 10−3 7.58 × 10−3 3.69 × 10−5
(0.6, 0.6) 1.33 × 10−2 8.78 × 10−3 6.41 × 10−5
(0.7, 0.7) 4.71 × 10−2 9.52 × 10−3 2.54 × 10−5
(0.8, 0.8) 1.09 × 10−1 7.56 × 10−3 3.15 × 10−4
(0.9, 0.9) 2.16 × 10−1 1.32 × 10−3 1.47 × 10−3
(1, 1) 3.87 × 10−1 2.42 × 10−2 4.46 × 10−3

7 Exact solution
r =r =1
1 2
6.5 r =r =1.2
1 2
r =r =1.4
1 2
6
r =r =1.6
1 2
5.5 r =r =1.8
1 2
r1=r2=2
5

4.5

3.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

Fig. 1 The approximate solution for various values of r1 and r2 with t = 1 and M = N = 5 of Example 4

Fig. 2 a The absolute error and b the approximate solution for r1 = r2 = 1 with M = 7, N = 5 of Example 4

123
259 Page 24 of 32 H. Dehestani et al.

Fig. 3 a The absolute error and b the approximate solution for r1 = r2 = 20 with M = N = 3 of Example 5

Example 5 In this example, consider the partial integral equation in the following form (Tari
et al. 2009):

u(x, t) = g(x, t) + Iθr ((xξ 2 + cos(η))u(ξ, η)), (x, t) ∈ [0, 1] × [0, 1], r = (r1 , r2 ) > 0,

with
 x  t
1
g(x, t) = x sin(t) − (x − ξ )r1 −1 (t − η)r2 −1 ξ sin(η)dξ dη.
Γ (r1 )Γ (r2 ) 0 0

The exact solution of this example is u(x, t) = x sin(t). The numerical results are displayed
in Table 5 and Fig. 3. The absolute errors for the various choices of M, N , r1 and r2 are
shown in Table 5. As seen in this table, the absolute error tends to zero when the number of
Genocchi polynomials increases. Figure 3 demonstrates the absolute error and approximate
solution obtained by the proposed method for r1 = r2 = 20 with M = N = 3.

Example 6 We consider the VO fractional partial integro-differential equation as follows:

ν(x,t) ∂u ∂ 2u ∂ 2u
Dt u(x, t) + (x, t) + 2 (x, t) + 2 (x, t) − u(x, t)
∂x ∂x ∂t
(1,1)
= g(x, t) + Iθ (cos(ξ η)u(ξ, η)), 0 < ν(x, t) ≤ 1,

with initial condition

∂u
u(x, 0) = x, (x, 0) = x 2 , 0 ≤ x ≤ 1,
∂t
and boundary conditions

u(0, t) = 0, u(1, t) = t 4 sin(1) + t + 1, 0 ≤ t ≤ 1,

where
Γ (5)t 2−ν(x,t) Γ (2)t 1−ν(x,t) 2
g(x, t) = sin(x) + x − t x 2 + (2t − 1)x + 2t + t 4 cos(x)
Γ (6 − ν(x, t)) Γ (2 − ν(x, t))
 x t
+12t 2 sin(x) − 2t 4 sin(x) + 1 − cos(ξ η)u(ξ, η)dξ dη.
0 0

123
 Table 5 The absolute error for different values of M, N , r1 and r2 of Example 5
(xi , ti ) r1 = r2 = 1 r1 = 0.5, r2 = 1.5 r1 = 20, r2 = 200
M = 3, N = 5 M = 3, N = 7 M = 3, N = 5 M = 3, N = 7 M = 3, N = 5 M = 3, N = 7

(0, 0) 1.00 × 10−7 6.00 × 10−9 6.00 × 10−7 4.80 × 10−8 1.11 × 10−16 0
(0.1, 0.1) 1.26 × 10−7 6.33 × 10−9 3.65 × 10−7 3.05 × 10−8 1.92 × 10−7 3.44 × 10−8
(0.2, 0.2) 3.92 × 10−7 5.10 × 10−9 1.11 × 10−7 8.53 × 10−8 4.10 × 10−7 5.93 × 10−8
(0.3, 0.3) 3.33 × 10−7 4.45 × 10−8 8.86 × 10−7 1.16 × 10−7 3.70 × 10−7 7.29 × 10−8
(0.4, 0.4) 4.93 × 10−8 1.11 × 10−7 6.82 × 10−7 1.22 × 10−7 1.54 × 10−7 7.57 × 10−8
(0.5, 0.5) 8.33 × 10−7 2.04 × 10−7 9.97 × 10−8 9.86 × 10−8 6.82 × 10−7 6.89 × 10−8
(0.6, 0.6) 5.39 × 10−7 3.22 × 10−7 1.25 × 10−6 4.73 × 10−8 3.72 × 10−7 4.74 × 10−8
The novel operational matrices based on 2D-Genocchi polynomials. . .

(0.7, 0.7) 4.34 × 10−6 4.50 × 10−7 4.36 × 10−6 4.13 × 10−8 2.38 × 10−6 1.37 × 10−8
(0.8, 0.8) 1.83 × 10−6 5.49 × 10−7 4.55 × 10−6 1.84 × 10−7 9.63 × 10−6 3.66 × 10−8
(0.9, 0.9) 5.74 × 10−5 6.67 × 10−7 6.76 × 10−5 2.65 × 10−7 7.98 × 10−5 2.75 × 10−7
(1, 1) 2.46 × 10−4 1.32 × 10−6 2.74 × 10−4 4.03 × 10−7 3.01 × 10−4 1.25 × 10−6

123
Page 25 of 32
259
259 Page 26 of 32 H. Dehestani et al.

Table 6 Errors for different values of M, N and ν(x, t) of Example 6

M N ν(x, t) e2 e∞

5 5 1 9.7883 × 10−5 8.2231 × 10−5

0.95 1.0090 × 10−4 8.5625 × 10−5
0.75 1.1193 × 10−4 9.6721 × 10−5
0.25 1.2024 × 10−4 1.0601 × 10−4
0.05 1.2347 × 10−4 1.0882 × 10−4
0.8 + 0.2 exp(−x) sin(t) 1.0778 × 10−4 9.2715 × 10−5
7 5 1 2.9500 × 10−6 2.0858 × 10−6
0.95 2.6241 × 10−6 2.1118 × 10−6
0.75 2.1118 × 10−6 2.1818 × 10−6
0.25 2.3594 × 10−6 1.8563 × 10−6
0.05 2.5614 × 10−6 1.7918 × 10−6
0.8 + 0.2 exp(−x) sin(t) 1.1583 × 10−6 9.4881 × 10−7

The exact solution of this example is u(x, t) = t 4 sin(x) + t x 2 + x. The numerical results are
displayed in Table 6 and Fig. 4. In Table 6, for various values of M, N the following errors
are calculated:

e∞ = max |u(x, t) − u M N (x, t)|,

(x,t)∈Ω
  21

M 
N
e2 = (u(xm , tn ) − u M N (xm , tn ))2 .
m=1 n=1

The calculated results in Table 6 show that the approximate solution has a good agreement
with the exact solution. Also, the absolute error for various choices of ν(x, t) and behavior
of the approximate solution with exact solution for M = 3, N = 5 at time t = 1 are given
in Fig. 4a, b, respectively.

Example 7 We consider the partial integro-differential equation as follows: (Tari et al. 2009):

∂u ∂ 2u ∂ 2u
(x, t) + 2 (x, t) + 2 (x, t) + u(x, t)
∂t ∂x ∂t
−Iθr ((xξ + t sin(η))u(ξ, η)) = g(x, t), r = (r1 , r2 ) > 0,

with initial condition

∂u
u(x, 0) = x, u(0, t) + (0, t) = 1 + t, 0 ≤ x ≤ 1.
∂t
g(x, t) is defined such that the exact solution of this example is u(x, t) = x cos(t) + t. The
numerical results are displayed in Tables 7 and 8. These tables show the absolute errors by
applying the proposed method. From these tables, we can see that the current method has
good accuracy.

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 27 of 32 259

10 -3 (a) (b)
20
3
(x,t)=0.95 Exact solution
(x,t)=0.75 approximate solution
15 (x,t)=0.5 2.5
(x,t)=0.25
(x,t)=0.8+0.2exp(-x)sin(t) 2
10
-3
10 1.5
4

5
3.5 1

3
0 0.2 0.25 0.5

0
-5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X

Fig. 4 a The absolute error for the different choices of ν(x, t), b the approximate solution and exact solution
for ν(x, t) = 1, with M = 3, N = 5 and t = 1 for Example 6

Table 7 Comparison of the absolute errors for different values of M, N and r1 = r2 = 1 with method in Ref.
(Tari et al. 2009) of Example 7
(xi , ti ) Present method Method in Tari et al. (2009)
M = 2, N = 5 M = 2, N = 7 M = 2, N = 9 N = 10

(0.2, 0.3) 4.5050 × 10−6 2.4254 × 10−7 1.8094 × 10−7 4.5567 × 10−7
(0.4, 0.6) 1.7046 × 10−5 2.2253 × 10−7 2.0896 × 10−7 9.6668 × 10−7
(0.5, 0.8) 7.3611 × 10−5 4.2005 × 10−8 3.0924 × 10−7 6.4149 × 10−7
(0.7, 0.5) 2.3236 × 10−5 5.4154 × 10−7 5.4592 × 10−7 2.9899 × 10−6
(1, 1) 6.6639 × 10−4 4.1389 × 10−6 5.4806 × 10−7 1.9872 × 10−6

Table 8 The absolute error for different values of r1 , r2 with M = 2, N = 7 of Example 7

(xi , ti ) r1 = 0.9, r2 = 0.9 r1 = 2.8, r2 = 0.9 r1 = 2.8, r2 = 2.8 r1 = 2.8, r2 = 6.9

(0.1, 0.1) 4.1786 × 10−10 6.5107 × 10−8 1.5217 × 10−7 5.0058 × 10−6
(0.3, 0.3) 1.0012 × 10−7 3.4469 × 10−7 3.3358 × 10−7 4.7096 × 10−6
(0.5, 0.5) 3.4238 × 10−7 7.8348 × 10−7 3.6157 × 10−7 1.1960 × 10−5
(0.7, 0.7) 7.2070 × 10−7 1.3264 × 10−6 2.3532 × 10−7 2.0872 × 10−5
(0.9, 0.9) 3.3940 × 10−7 4.1389 × 10−6 9.5909 × 10−7 2.9332 × 10−5

Example 8 We consider the partial integral system as follows (Kajani and Shehni 2011):
⎧
⎨ u 1 (x, t) + Iθr (−u 1 (ξ, η) − u 2 (ξ, η)) = x sin(t) − 21 x 2 + 21 x 2 cos(t) − 21 x 2 sin(t),
⎩ ∗
u 2 (x, t) + Iθr (−u 1 (ξ, η) + u 2 (ξ, η)) = x cos(t) − 21 x 2 + 21 x 2 sin(t) − 21 x 2 cos(t),
r = (r1 , r2 ) > 0, r ∗ = (r3 , r4 ) > 0.
The exact solution of this example, when r = r ∗ = 1, is u 1 (x, t) = x sin(t) and u 2 (x, t) =
x cos(t). The numerical results are displayed in Table 9 and Figs. 5 and 6. In Table 9, we

123
259 Page 28 of 32 H. Dehestani et al.

1.5
Exact solution
r1=r2=r3=r4=1
r1=r2=r3=r4=0.9
r1=r2=r3=r4=0.8
1
r1=0.8,r2=0.9,r3=0.8,r4=0.9
r =r =r =r =1.2
1 2 3 4
r =0.8,r =1.2,r =1.2,r =0.8
1 2 3 4

0.5

0
0 0.2 0.4 0.6 0.8 1
x

Fig. 5 The approximate solution of u 1 (x, t) for different values of r1 , r2 with M = 3, N = 7 and t = 1 for
Example 8

0.5

0.4

0.3 Exact solution

r =r =r =r =1
1 2 3 4

0.2
r =r =r =r =0.9
1 2 3 4
r =r =r =r =0.8
1 2 3 4

0.1 r =0.8,r =0.9,r =0.8,r =0.9

1 2 3 4
r1=r2=r3=r4=1.2
0 r1=0.8,r2=1.2,r3=1.2,r4=0.8

0 0.2 0.4 0.6 0.8 1

X

Fig. 6 The approximate solution of u 2 (x, t) for different values of r1 , r2 with M = 3, N = 7 and t = 1 for
Example 8

compare the absolute errors between the exact and approximate solutions using the present
method for various values of M, N and r1 = r2 = r3 = r4 = 1, with the method in (Kajani
and Shehni 2011). From Table 9, it is seen that the results obtained by the presented method
are superior to the method in (Kajani and Shehni 2011). Moreover, if higher values of M and
N are chosen, the absolute errors decay rapidly. Figures 5 and 6 demonstrate that if r and r ∗
are chosen near to 1, the computational results have better accuracy.

123
 Table 9 Comparison of the absolute error for different values of M, N and r = r ∗ = 1 with method in Ref. (Kajani and Shehni 2011) of Example 8
(xi , ti ) Present method Method in Kajani and Shehni (2011)
M = 3, N = 7 M = 3, N = 9 N = 10
e1 (x, t) e2 (x, t) e1 (x, t) e2 (x, t) e1 (x, t) e2 (x, t)

(0.1, 0.6) 2.01 × 10−9 1.02 × 10−7 8.43 × 10−14 2.09 × 10−14 0.0 1.0 × 10−9
(0.1, 0.9) 4.53 × 10−8 1.63 × 10−7 1.62 × 10−11 1.62 × 10−11 1 × 10−9 9.0 × 10−9
(0.2, 0.7) 2.13 × 10−8 2.39 × 10−7 3.57 × 10−13 2.67 × 10−13 1 × 10−9 2.0 × 10−9
(0.3, 0.8) 4.09 × 10−7 4.09 × 10−7 1.38 × 10−12 1.38 × 10−12 1 × 10−9 9.0 × 10−9
(0.3, 1) 4.96 × 10−7 7.29 × 10−7 9.20 × 10−10 5.13 × 10−10 8 × 10−9 8.2 × 10−8
The novel operational matrices based on 2D-Genocchi polynomials. . .

(0.4, 0.8) 5.46 × 10−7 2.09 × 10−7 3.99 × 10−12 1.79 × 10−12 1 × 10−9 1.2 × 10−8
(0.5, 0.8) 6.82 × 10−7 6.83 × 10−7 5.05 × 10−12 2.17 × 10−12 2 × 10−9 1.4 × 10−8
(1, 0.6) 1.02 × 10−7 1.02 × 10−7 1.85 × 10−12 1.18 × 10−12 0.0 2.0 × 10−9
(1, 1) 2.44 × 10−7 2.45 × 10−7 3.05 × 10−9 1.74 × 10−9 2.5 × 10−8 2.7 × 10−7

123
Page 29 of 32
259
 259 Page 30 of 32 H. Dehestani et al.

9 Conclusion

In this article, based on 2D-Genocchi polynomials, collocation method and the operational
and pseudo-operational matrices of integration and derivative, we approximate solution of
VO-FPIDEs. In addition, the dual pseudo-operational matrix of fractional order integration is
used to approximate an integral part of VO-FPIDEs. The technique of obtaining the pseudo-
operational matrices causes to achieve the approximate solution with high accuracy. Also,
the suggested method demonstrates that by applying a small number of the 2D-Genocchi
polynomials, we get a satisfactory result. Meanwhile, an upper bound for the error of pseudo-
operational matrix of fractional integration was derived by using Genocchi polynomials and
properties of the Riemann–Liouville fractional integral operator. Moreover, we presented
the numerical results and the graphical representations to verify the validity of the proposed
method.

Acknowledgements This work is supported by the National Elite Foundation and Alzahra University. We
express our sincere thanks to the anonymous referees for valuable suggestions that improved the final
manuscript.

Compliance with ethical standards

Conflict of interest The authors declare that they have no conflict of interest.

References
Aghazadeh N, Khajehnasiri AA (2013) Solving nonlinear two-dimensional Volterra integro-differential equa-
tions by block-pulse functions. J Math Sci 7(1):3
Araci S, Acikgoz M, Sen E (2013) On the extended Kims p-adic q-deformed fermionic integrals in the p-adic
integer ring. J Numb Theory 133(10):3348–3361
Avazzadeh Z, Rizi ZB, Ghaini FMM, Loghmani GB (2012) A numerical solution of nonlinear parabolic-
type Volterra partial integro-differential equations using radial basis functions. Eng Anal Bound Elem
36(5):881–893
Babaei A, Jafari H, Banihashemi S (2020) Numerical solution of variable order fractional nonlinear quadratic
integro-differential equations based on the sixth-kind Chebyshev collocation method. J Comput Appl
Math 377:112908
Babolian E, Dastani N (2011) Numerical solutions of two-dimensional linear and nonlinear Volterra integral
equations: homotopy perturbation method and differential transform method. J Indu Math 3(3):157–167
Beyer H, Kempfle S (1995) Definition of physically consistent damping laws with fractional derivatives. Z
Angew Math Mech 75(8):623–635
Bloom F (1981) Ill-posed problems for integrodifferential equations in mechanics and electromagnetic theory.
Siam 3
Chen Y, Liu L, Li B, Sun Y (2014) Numerical solution for the variable order linear cable equation with
Bernstein polynomials. Appl Math Comput 238:329–341
Coimbra CFM (2003) Mechanics with variable-order differential operators. Ann Phys 12(11–12):692–703
Dehestani H, Ordokhani Y, Razzaghi M (2018) Fractional-order Legendre-Laguerre functions and their appli-
cations in fractional partial differential equations. Appl Math Comput 336:433–453
Dehestani H, Ordokhani Y, Razzaghi M (2019) On the applicability of Genocchi wavelet method for different
kinds of fractional order differential equations with delay. Numer Linear Algebra Appl 2259(26):1–29
Dehestani H, Ordokhani Y, Razzaghi M (2019) Application of the modified operational matrices in multiterm
variable-order time-fractional partial differential equations. Math Meth Appl Sci 42:7296–7313
Dehestani H, Ordokhani Y, Razzaghi M (2019) numerical technique for solving various kinds of fractional
partial differential equations via Genocchi hybrid functions. RACSAM 113:3297–3321
Dehestani H, Ordokhani Y, Razzaghi M (2020) Pseudo-operational matrix method for the solution of variable-
order fractional partial integro-differential equations. Eng Comput. https://doi.org/10.1007/s00366-019-
00912-z

123
The novel operational matrices based on 2D-Genocchi polynomials. . . Page 31 of 32 259

Dehestani H, Ordokhani Y, Razzaghi M (2020) The computational method for generalized fractional Benjamin-
Bona-Mahony-Burgers equations arising from the propagation of water waves. Sadhana 45(1):1–20
Diaz G, Coimbra CFM (2009) Nonlinear dynamics and control of a variable order oscillator with application
to the van der Pol equation. Nonlinear Dyn 56(1–2):145–157
Engler H (1983) On some parabolic integro-differential equations: Existence and asymptotics of solutions.
Springer, Berlin, pp 161–167
Fakhar-Izadi F, Dehghan M (2011) The spectral methods for parabolic Volterra integro-differential equations.
J Comput Appl Math 235(14):4032–4046
Ganji RM, Jafari H, Baleanu D (2020) A new approach for solving multi variable orders differential equations
with Mittag-Leffler kernel. Chaos Solitons Fract 130:109405
Ganji RM, Jafari H (2020) A numerical approach for multi-variable orders differential equations using Jacobi
polynomials. Int J Appl Comput Math 5(2):34
Gorenflo R, Francesco M (1997) Fractional calculus. Springer, Berlin, pp 223–276
Grasselli M, Kabanikhin SI, Lorenzi A (1990) An inverse hyperbolic integrodifferential problem arising in
geophysics II. Nonlinear Anal Theory Methods Appl 15(3):283–298
Greenwell-Yanik CE, Fairweather G (1986) Analyses of spline collocation methods for parabolic and hyper-
bolic problems in two space variables. SIAM J Numer Anal 23(2):282–296
Habetler GJ, Schiffman RL (1970) A finite difference method for analyzing the compression of poro-
viscoelastic media. Computing 6(3–4):342–348
Heydari MH (2016) A new approach of the Chebyshev wavelets for the variable-order time fractional mobile-
immobile advection dispersion model. arXiv preprint arXiv:1605.06332
Ingman D, Suzdalnitsky J (2004) Control of damping oscillations by fractional differential operator with
time-dependent order. Comput Methods Appl Mech Eng 193(52):5585–5595
Isah A, Phang C (2018) Operational matrix based on Genocchi polynomials for solution of delay differential
equations. Ain Shams Eng 9(4):2123–2128
Isah A, Phang C (2019) New operational matrix of derivative for solving non-linear fractional differential
equations via Genocchi polynomials. J King Saud Univ Sci 31:1–7
Isah A, Phang C (2016) Genocchi Wavelet-like operational matrix and its application for solving non-linear
fractional differential equations. Open Phys 14(1):463–472
Jang B (2009) Comments on solving a class of two-dimensional linear and nonlinear Volterra integral equations
by the differential transform method. J Comput Appl Math 233(2):224–230
Jiang W, Liu N (2017) A numerical method for solving the time variable fractional order mobile-immobile
advection-dispersion model. Appl Numer Math 119:18–32
Kajani MT, Shehni NA (2011) Solutions of two-dimensional integral equation systems by using differential
transform method. Am Math 11:74–77
Khajehnasiri AA (2016) Numerical solution of nonlinear 2D Volterra-Fredholm integro-differential equations
by two-dimensional triangular function. J Appl Comput Math 2(4):575–591
Kreyszig E (1978) Introductory functional analysis with applications. Wiley, New York
Kulish VV, Lage JL (2002) Application of fractional calculus to fluid mechanics. J Fluids Eng 12(3):803–806
Larsson S, Thomée V, Wahlbin L (1998) Numerical solution of parabolic integro-differential equations by the
discontinuous Galerkin method. Math Comput 67(221):45–71
Loh JR, Phang C (2018) A new numerical scheme for solving system of Volterra integro-differential equation.
Alex Eng J 57:1117–1124
Li X, Wu B (2015) A numerical technique for variable fractional functional boundary value problems. Appl
Math Lett 43:108–113
Ma J (2007) Finite element methods for partial Volterra integro-differential equations on two-dimensional
unbounded spatial domains. Appl Math Comput 186(1):598–609
Magin RL, Ingo C, Colon-Perez L, Triplett W, Mareci TH (2013) Characterization of anomalous diffusion in
porous biological tissues using fractional order derivatives and entropy. Microporous Mesoporous Mater
178:39–43
Meleshko SV, Grigoriev YN, Ibragimov NK, Kovalev VF (2010) Symmetries of integro-differential equations:
with applications in mechanics and plasma physics. Springer, Berlin
Moghaddam BP, Machado JT (2017) A stable three-level explicit spline finite difference scheme for a class of
nonlinear time variable order fractional partial differential equations. Comput Math Appl 73(6):1262–
1269
Najafalizadeh S, Ezzati R (2016) Numerical methods for solving two-dimensional nonlinear integral equations
of fractional order by using two-dimensional block pulse operational matrix. Appl Math Comput 280:46–
56
Nemati S, Ordokhani Y (2013) Legendre expansion methods for the numerical solution of nonlinear 2D
Fredholm integral equations of the second kind. J Appl Math Inform 31(5–6):609–621

123
 259 Page 32 of 32 H. Dehestani et al.

Ostalczyk P, Rybicki T (2008) Variable-fractional-order dead-beat control of an electromagnetic servo. J Vib

Control 14(9–10):1457–1471
Pachpatte BG (1983) On a nonlinear diffusion system arising in reactor dynamics. J Math Anal Appl 94(2):501–
508
Pao CV, Payne L, Amann H (1979) Bifurcation analysis on a nonlinear diffusion system in reactor dynamics.
Appl Anal. 9(2):107–119
Pao CV (1974) Solution of a nonlinear integrodifferential system arising in nuclear reactor dynamics. J Math
Anal Appl 48(2):470–492
Ramirez LES, Coimbra CFM (2007) A variable order constitutive relation for viscoelasticity. Ann Phys 16(7–
8):543–552
Rivlin TJ (1981) An introduction to the approximation of functions. Courier Corporation, Chelmsford
Samko SG (1995) Fractional integration and differentiation of variable order. Anal Math 21(3):213–236
Sadeghi S, Jafari H, Nemati S (2020) Operational matrix for Atangana–Baleanu derivative based on Genocchi
polynomials for solving FDEs. Chaos Solitons Fract 135:109736
Schafer RD (1996) An introduction to nonassociative algebras. Dover, New York, p 12
Shu C, Ding H, Yeo KS (2004) Solution of partial differential equations by a global radial basis function-based
differential quadrature method. Eng Anal Bound Elem 28(10):1217–1226
Soon CM, Coimbra CFM, Kobayashi MH (2005) The variable viscoelasticity oscillator. Ann Phys 14(6):378–
389
Sun HG, Chen W, Wei H, Chen YQ (2011) A comparative study of constant-order and variable-order fractional
models in characterizing memory property of systems. Eur Phys J Spec Top 193(1):185
Tari A, Rahimi MY, Shahmorad S, Talati F (2009) Development of the tau method for the numerical solution
of two-dimensional linear Volterra integro-differential equations. Comput Methods Appl Math 9(4):421–
435
Torvik PJ, Bagley RL (1984) On the appearance of the fractional derivative in the behavior of real materials.
J Appl Mech 51(2):294–298
Vityuk AN, Golushkov AV (2004) Existence of solutions of systems of partial differential equations of frac-
tional order. Nonlinear Oscil 7(3):318–325
Yaghoobi S, Moghaddam BP, Ivaz K (2017) An efficient cubic spline approximation for variable-order frac-
tional differential equations with time delay. Nonlinear Dyn 87(2):815–826
Yan YI, Fairweather G (1992) Orthogonal spline collocation methods for some partial integro-differential
equations. SIAM J Numer Anal 29(3):755–768
Yanik EG, Fairweather G (1988) Finite element methods for parabolic and hyperbolic partial integro-
differential equations. Nonlinear Anal 12(8):785–809
Zhang H, Liu F, Phanikumar MS, Meerschaert MM (2013) A novel numerical method for the time variable
fractional order mobile-immobile advection-dispersion model. Comput Math Appl 66(5):693–701

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and
institutional affiliations.

123