A computational method for a class of systems of nonlinear variable-order fractional quadratic integral equations

Applied Numerical Mathematics 153 (2020) 164–178

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A computational method for a class of systems of nonlinear variable-order fractional quadratic integral equations M.H. Heydari Department of Mathematics, Shiraz University of Technology, Shiraz, Iran

a r t i c l e

i n f o

Article history: Received 20 May 2019 Received in revised form 2 February 2020 Accepted 13 February 2020 Available online 19 February 2020 Keywords: Quadratic integral equation (QIE) Chebyshev cardinal functions (CCFs) Operational matrix (OM) Convergence and error analysis

a b s t r a c t This paper develops an accurate computational method based on the shifted Chebyshev cardinal functions (CCFs) for a new class of systems of nonlinear variable-order fractional quadratic integral equations (QIEs). In this way, a new operational matrix (OM) of variable-order fractional integration is obtained for these cardinal functions. In the proposed method, the unknown functions of a system of nonlinear variable-order fractional QIEs are approximated by the shifted CCFs with undetermined coefficients. Then, these approximations are substituted into the system. Next, the OM of variable-order fractional integration and the cardinal property of the shifted CCFs are utilized to reduce the system into an equivalent system of nonlinear algebraic equations. Finally, by solving this algebraic system an approximate solution for the problem is obtained. The main idea behind this approach is to reduce such problems to solving systems of nonlinear algebraic equations, which greatly simplifies the problem. Convergence of the presented method is investigated theoretically and numerically. Furthermore, the proposed approach is numerically evaluated by solving some test problems. © 2020 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Over the last decades, fractional integral equations have achieved great popularity among the mathematicians and physicists in modeling applied problems, for instance see [9,16], and since it is difficult to solve integral equations analytically, numerical methods have been extensively used on these equations, e.g. [18,24]. Quadratic integral equations (QIEs) have been appeared in a variety of areas such as kinetic theory of gases, the theory of radiative transfer and the traffic theory [3]. In recent years, several numerical methods have been successfully used to numerically solve this important class of equations, for instance see [5,6]. Recently, many researchers have studied fractional form of the QIEs, e.g. [7,8]. Variable-order fractional calculus (i.e. the theory of differential and integral operators with variable fractional orders) can be a useful mathematical framework for a deeper and more precise analysis of dynamical phenomena [20,22]. Note that the mathematical systems including variable-order fractional operators are too complex to be analytically solved and so numerical approaches are applied for obtaining their approximate solutions. For instance, some effective numerical methods are investigated in [2,4,19,23,25] for solving the systems involved with the variable-order fractional derivatives. It should be noted that there are not various methods for solution of variable-order fractional integral equations. In [17], Khane Keshi et al. have used a discretization scheme for numerical solution of a type of variable-order fractional integral equations.

E-mail address: [email protected]. https://doi.org/10.1016/j.apnum.2020.02.011 0168-9274/© 2020 IMACS. Published by Elsevier B.V. All rights reserved.

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

165

Recently, Heydari [11] has proposed a wavelet method for the numerical solution of a new category of nonlinear variableorder fractional QIEs. We remind that interpolation is one of the basic concepts in numerical analysis for finding approximate solutions of mathematical problems. In the interpolation methods, by replacing the interpolants instead of unknown functions and substituting the collocation points an integral/differential equation can be reduced into a system of algebraic equations. In particular, interpolants defined on cardinal basis functions can greatly simplify the equations including differential/integral operators because the final system of algebraic equations are notably sparse. Some successful approximation techniques which depend on particular properties of cardinal functions have been recently published [11–13]. To the best of our knowledge, there is not any previous study focused on the cardinal methods for systems of variableorder fractional QIEs. The main purpose of this paper is to propose a high accurate method by utilizing the shifted Chebyshev cardinal functions (CCFs) for solution of a novel class of systems of nonlinear variable-order fractional QIEs in the form

⎛ uk (t ) = gk (t )+

1

(αk (t )) 1

(βk (t ))

⎞

t

⎝ (t − τ )αk (t )−1 Lk1 (t , τ , u 1 (τ ), u 2 (τ ), . . . , um (τ )) dτ ⎠ 0

⎞ ⎛ t  ⎝ (t − τ )βk (t )−1 Lk2 (t , τ , u 1 (τ ), u 2 (τ ), . . . , um (τ )) dτ ⎠ , k = 1, 2, . . . , m,

(1.1)

0

where t ∈ [0, t max ], and for k = 1, 2, . . . m, gk ’s are known functions, αk ’s and βk ’s are positive continuous functions, and Lk1 and Lk2 are Lipschitz continuous functions. To this end, first a new operational matrix (OM) of variable-order fractional integration is extracted for the shifted CCFs. In the suggested approach, the unknown functions in Eq. (1.1) are expanded in the shifted CCFs. Then, the cardinal property of the shifted CCFs, together with the derived OM are utilized to reduce the system into a system of algebraic equations. Eventually, the numerical solutions are computed by solving the algebraic system. The convergence of the proposed approach is theoretically and numerically examined to show its reliability. It is worth noting that Eq. (1.1) can be used for expressing the variable-order fractional model of the problems arisen in the radiative transfer. This study is arranged as follows: The main definitions and some preliminaries of the variable-order fractional calculus are given in Section 2. Section 3 contains some definitions and propositions about the shifted CCFs, and the new relevant results are expressed in Section 4. The proposed approach is formulated in Section 5, and its convergence and error analysis are investigated in Section 6. Some test problems are studied in Section 7. Eventually, a brief conclusion is given in Section 8. 2. Variable-order fractional calculus Herein, some definitions about the variable-order fractional calculus theory are briefly reviewed. Definition 2.1. ([2]) Let α : [0, τmax ] −→ [0, ∞) is a continuous function. The fractional integration of order function u (τ ) in the Riemann-Liouville sense is defined by

I α (τ ) u (τ ) =

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩

τ

1

(τ − ζ )α (τ )−1 u (ζ )dζ, α (τ ) > 0,

(α (τ )) 0

α (τ ) = 0.

u (τ ),

Corollary 2.2. The above definition results in

I α (τ ) τ k =

k!

(k + α (τ ) + 1)

τ k +α ( τ ) ,

k ∈ N ∪ {0}.

Definition 2.3. ([21]) The Mittag-Leffler function of two positive parameters λ1 and λ2 is expressed by

Eλ1 ,λ2 (t ) =

∞  j =0

tj

( j λ1 + λ2 )

,

t ∈ R.

Corollary 2.4. Let α : [0, τmax ] −→ R+ is a continuous function. Then, we have

I α (τ ) e τ = τ α (τ ) E1,α (τ )+1 (τ ),

I α (τ ) e −τ = τ α (τ ) E1,α (τ )+1 (−τ ),

I α (τ ) sin(τ ) = τ α (τ )+1 E2,α (τ )+2 (−τ 2 ),

I α (τ ) cos(τ ) = τ α (τ ) E2,α (τ )+1 (−τ 2 ).

α (τ ) of the

166

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

Proof. By using the Taylor’s series of the above mentioned functions, Corollary 2.2 and Definition 2.3, the proof is straightforward.  3. Shifted Chebyshev cardinal functions (CCFs) The shifted CCFs of order n are expressed on the interval [0, τmax ] by [10] n +1

ψi (τ ) =

k =1 k=i



τ − τk , τi − τk

i = 1, 2, . . . , n + 1,

(3.1)

  (2i −1)π 1 − ξ and ξ = cos for i = 1, 2, . . . , n + 1 is the i-th root of the Chebyshev polynomial of (n + 1)( ) i i 2 2(n+1) th degree on [−1, 1] [1]. A function u ∈ C [0, τmax ] can be expanded by the above basis functions as follows where

τi =

τmax

u (τ ) 

n +1 

c i ψi (τ )  C T n (τ ),

(3.2)

i =1

where

C = [c 1 c 2 . . . cm ] T ,

(3.3)

n (τ ) = [ψ1 (τ ) ψ2 (τ ) . . . ψn+1 (τ )] T ,

and c i = u (τi ). Similarly, any continuous function u (t , τ ) defined over [0, τmax ] × [0, τmax ] can be approximated by the shifted CCFs as follows [14]

u (t , τ ) 

n +1 n +1  

u i j ψi (t )ψ j (τ )  n (t ) T Un (τ ),

(3.4)

i =1 j =1

where U = [u i j ] and u i j = u (τi , τ j ) for i = 1, 2, . . . , n + 1 and j = 1, 2, . . . , n + 1. Lemma 3.1. Suppose that L : [0, τmax ]2 × Rm −→ R be a continuous operator and CkT n (τ ) for k = 1, 2, . . . , m be the approxima-



tions of the functions uk (τ ) by the shifted CCFs, where Ck = ck1 ck2 . . . ck(n+1)

T

. Then, we have

L (t , τ , u 1 (τ ), u 2 (τ ), . . . , um (τ ))  n (t )T Vn (τ ),

(3.5)

where V = [ v i j ] is an (n + 1)-order square matrix, and its entires are given by

vij = L





τi , τ j , c1 j , c2 j , . . . , cmj .

Proof. According to the cardinal property of the shifted CCFs, the proof is clear.



4. The OM of variable-order fractional integration The aim of this section is to show that the fractional integration of order be approximated as follows

α (τ ) > 0 of n (τ ) introduced in Eq. (3.3) can

I α (τ ) n (τ )  P(α ) n (τ ), where the (n + 1)-order square matrix P(α ) is the OM of fractional integration of variable-order

α (τ ).

Remark 1. ([15]) The basis functions introduced in Eq. (3.1) can be represented as follows n 1 

ψi (τ ) =

γi

b ik τ n−k ,

i = 1, 2, . . . , n + 1,

where for i = 1, 2, . . . , n + 1

γi =

n +1

(τi − τl ),

l =1 l=i

(4.1)

k =0

b ik =

⎧ 1 ⎪ ⎪ ⎨ 1 ⎪ ⎪ ⎩ −k

k = 0, k  l =1

ail b ik−l ,

k = 1, 2, . . . , n ,

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

167

in which

ail =

n +1 

τrl ,

l = 1, 2, . . . , k .

r =1 r =i

Theorem 4.1. Assume that n (τ ) be the vector introduced in Eq. (3.3) and α : [0, τmax ] −→ R+ be a continuous function. The fractional integration of order α (τ ) of this vector can be approximated as follows

I α (τ ) n (τ )  P(α ) n (τ ),

(4.2)

(α )

where P(α ) = [ p i j ] is an (n + 1)-order square matrix, and its elements are computed by (α )

pi j =

n 1 

γi

n−k+α (τ j )

(n − k)!τ j

b ik

(n − k + α (τ j ) + 1)

k =0

,

in which the coefficients b ik were introduced in Remark 1. Proof. From the analytical form expressed in Eq. (4.1) and the linearity of the variable-order fractional integration, we have

 I α (τ ) ψi (τ ) = I α (τ )

n 1 

γi

 b ik τ

n−k

=

k =0

n 1 

γi

b ik I α (τ ) τ n−k ,

i = 1, 2, . . . , n + 1.

k =0

Moreover from Corollary 2.2, the above relation results in n 1 

I α (τ ) ψi (τ ) =

γi

b ik

k =0

(n − k)! τ n−k+α (τ ) , (n − k + α (τ ) + 1)

i = 1, 2, . . . , n + 1.

Now, by expanding the right-hand side of the above relation in terms of the shifted CCFs, we obtain

I α (τ ) ψi (τ ) 

n +1 

(α )

i = 1, 2, . . . , n + 1,

p i j ψ j (τ ),

j =1

where (α )

pi j =

n 1 

γi

n−k+α (τ j )

b ik

k =0

(n − k)! τ j

(n − k + α (τ j ) + 1)

which completes the proof.

,

i = 1, 2, . . . , n + 1, j = 1, 2, . . . , n + 1,



5. Explanation of the proposed method To solve the system introduced in Eq. (1.1), the known and unknown functions therein are approximated by the shifted CCFs with determined and undetermined coefficients, respectively. Then, some useful properties of the shifted CCFs are utilized to transform the system into an equivalent system of algebraic equations that is easier to solve. Eventually, an approximate solution is obtained for the problem by solving the yielded algebraic system. Accordingly, suppose that

uk (t )  CkT n (t ),



k = 1, 2, . . . , m ,

(5.1)

T

where Ck = ck1 ck2 . . . ck(n+1) and n (t ) is the vector introduced in Eq. (3.3). Also, the shifted CCFs are applied to approximate the functions gk (t ) for k = 1, 2, . . . , m as follows

gk (t )  GkT n (t ),

(5.2)

where Gk for k = 1, 2, . . . , m is an (n + 1) dimensional column vector with known elements. Replacing from Eqs. (5.1) and (5.2) into Eq. (1.1), yields

168

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

CkT n (t )



GkT n (t )+

1

(αk (t )) 1

(βk (t ))

⎞ ⎛ t    α ( t )− 1 T T T ⎝ (t − τ ) k Lk1 t , τ , C1 n (τ ), C2 n (τ ), . . . , Cm n (τ ) dτ ⎠ 0

⎞ ⎛ t    T ⎝ (t − τ )βk (t )−1 Lk2 t , τ , C1T n (τ ), C2T n (τ ), . . . , Cm n (τ ) dτ ⎠ ,

(5.3)

0

where k = 1, 2, . . . , m. Moreover, applying Lemma 3.1, results in

  T Lk1 t , τ , C1T n (τ ), C2T n (τ ), . . . , Cm n (τ )  n (t )T Vk1 n (τ ),   T Lk2 t , τ , C1T n (τ ), C2T n (τ ), . . . , Cm n (τ )  n (t )T Vk2 n (τ ),

(5.4)

where Vk1 and Vk2 for k = 1, 2, . . . , m are (n + 1)-order square matrices with elements

Vk1 = [ v k1 ]i j = Lk1 Vk2 = [ v k2 ]i j = Lk2

 



τi , τ j , c1 j , c2 j , . . . , cmj , 

(5.5)

τi , τ j , c1 j , c2 j , . . . , cmj .

Substituting from Eq. (5.4) into Eq. (5.3) yields

⎛ CkT n (t )



GkT n (t )+ ⎝n (t ) T Vk1

1

(αk (t ))

⎛ ⎝n (t )T Vk2

1

(βk (t ))

t

⎞ (t − τ )αk (t )−1 n (τ )dτ ⎠

0

t

⎞ (t − τ )βk (t )−1 n (τ )dτ ⎠ ,

(5.6)

0

where k = 1, 2, . . . , m. By applying Theorem 4.1, the above relation can be represented in the form





CkT n (t )  GkT n (t ) + n (t ) T Vk1 P(αk ) n (t )

 n (t )T Vk2 P(βk ) n (t ) ,

k = 1, 2, . . . , m .

(5.7)

Moreover, by applying the cardinal property of the basis functions for k = 1, 2, . . . , m, we get

  T n (t )T Vk1 P(αk ) n (t )  diag Vk1 P(αk ) n (t )  k1 n (t ),   T n (t )T Vk2 P(βk ) n (t )  diag Vk2 P(βk ) n (t )  k2 n (t ),

(5.8)

and consequently



  n (t )T Vk1 P(αk ) n (t ) n (t )T Vk2 P(βk ) n (t )  kT (t ),

(5.9)

where k = k1  k2 , and  indicates the point-wise product of k1 and k2 . Eqs. (5.7)-(5.9) result in





CkT − kT n (t )  GkT n (t ),

k = 1, 2, . . . , m ,

(5.10)

and the cardinal property of the shifted CCFs implies that

CkT − kT = GkT ,

k = 1, 2, . . . , m .

(5.11)

Note that the system produced in Eq. (5.11) can be easily solved by suitable algorithms to compute Ck for k = 1, 2, . . . , m. Finally, the approximate solution of the main problem is obtained by replacing Ck into Eq. (5.1). To this end the “fsolve”

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

169

command in MAPLE is employed for solution of the system achieved in Eq. (5.11). The method is expressed in the following algorithm: Algorithm of the presented method. Input: The numbers τmax ∈ R+ , m, n ∈ N ; the functions gk : [0, τmax ] −→ R, αk , βk : [0, τmax ] −→ R+ and Lk1 , Lk2 : [0, τmax ]2 × Rm −→ R for k = 1, 2, . . . , m.  T  T T Output: The approximatesolution:  u 1 (t ), u 2 (t ), . . . , um (t )  C1 n (t ), C2 n (t ), . . . , Cm n (t ) . (2i −1)π 2(n+1)

Step 1: Generate ξi = cos

and

τi =

τmax 2

(1 − ξi ) for i = 1, 2, . . . , n + 1.

Step 2: Define the functions ψi (τ ) for i = 1, 2, . . . , n + 1 by Eq. (3.1). Step 3: Construct the shifted CCFs vector n (τ ) by Eq. (3.3). Step 4: Compute the operational matrices P(αk ) and P(βk ) for k = 1, 2, . . . , m using Theorem 4.1. Step 5: Define the vectors Ck for k = 1, 2, . . . , m with undetermined elements. Step 6: Compute the vectors Gk in Eq. (5.2) for k = 1, 2, . . . , m similar to Eq. (3.2). Step 7: Compute the matrices Lk1 and Lk2 in Eq. (5.4) for k = 1, 2, . . . , m by Eq. (5.5). Step 8: Compute the vectors k1 and k2 in Eq. (5.8), and k in Eq. (5.9) for k = 1, 2, . . . , m. Step 9: Put CkT − kT = GkT for k = 1, 2, . . . , m. Step 10: Solve the system achieved in Step 9 and evaluate the vectors Ck for k = 1, 2, . . . , m.

6. Convergence and error analysis Here, the convergence and error analysis of the presented method are theoretically examined. In the beginning, some necessary preliminaries and definitions are briefly reviewed. Remark 2.([10]) The shifted CCFs introduced in Eq. (3.1) are orthogonal over [0, τmax ] with respect to the weight function 2 ω(τ ) = 1/ 1 − ( τmax τ − 1)2 , and their orthogonality property is given by



τmax

 ψi (τ ), ψ j (τ ) ω =

ψi (τ )ψ j (τ )ω(τ )dτ = 0

πτmax 2(n + 1)

δi j ,

where δ indicates the Kronecker delta. ˆ ω(τ ) be as defined in Remark 2 and mˆ ∈ N ∪ {0}. The Sobolev space H m ω (0, τmax ) is defined by  ˆ 2 ˆ , u ( j ) ∈ L 2ω (0, τmax ) . Hm ω (0, τmax ) = u ∈ L ω (0, τmax ) : for j = 0, 1, . . . , m

Definition 6.1. ([1]) Let



Note that the introduced space is equipped with the weighted inner product

u , v mˆ ,ω =

ˆ τmax m 

u ( j ) (τ ) v ( j ) (τ )ω(τ )dτ ,

(6.1)

j =0 0

and is a Hilbert space with the introduced norm

u  H mˆ (0,τ w

max )

⎛ ˆ m  =⎝ u ( j ) 22

⎞1/2

L ω (0,τmax )

j =0

⎠

.

It is worth noting that the semi-norm associated with the above norm is given by

⎛ |u | H mˆ ;m˜ (0,τ w

max )

=⎝

ˆ m  ˆ ,m ˜ +1 ) j =min(m

⎞1/2 u ( j ) 2L 2 (0,τ ω

max )

ˆ Theorem 6.2. ([10]) Let u ∈ H m ω (0, τmax ) and P n u (τ ) =

⎠

n +1 

.

c i ψi (τ ), where c i = u (τi ). Then, we have

i =1

u − P n u L ∞ (0,τmax ) ≤ mˆ n1/2−mˆ

2

τmax

1/2 |u | H mˆ ;n (0,τ w

max )

,

(6.2)

170

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

where u  L ∞ (0,τmax ) =

sup

0≤τ ≤τmax

|u (τ )| and

⎛ |u | H mˆ ;n (0,τ w

max )

τ

ˆ m 

=⎝

max

⎞1/2

2 j

u ( j ) 2L 2 (0,τ

2

ˆ ,n+1) j =min(m

ω

max )

⎠

,

ˆ and independent on n. in which mˆ is a positive constant, dependent on m Remark 3. The following assumptions are imposed on Eq. (1.1):

• The kernel functions Lk1 and Lk2 for k = 1, 2, . . . , m are continuous on [0, τmax ]2 × Rm and satisfy the following Lipschitz conditions: m       Lk1 (t , τ , u 1 , u 2 , . . . , um ) − Lk1 t , τ , u˜ 1 , u˜ 2 , . . . , u˜ m  ≤ ρkl ul − u˜ l  , l =1

   Lk2 (t , τ , u 1 , u 2 , . . . , um ) − Lk2 t , τ , u˜ 1 , u˜ 2 , . . . , u˜ m  ≤

m 





σkl ul − u˜ l  ,

(6.3)

l =1

max |Lk1 (t , τ , 0, 0, . . . , 0)| = μk1 ,

I = [0, τmax ] × [0, τmax ] ,

max |Lk2 (t , τ , 0, 0, . . . , 0)| = μk2 ,

I = [0, τmax ] × [0, τmax ] ,

(t ,τ )∈ I (t ,τ )∈ I

where ρkl and σkl for l = 1, 2, . . . , m, and μk1 and μk2 are positive constants. • There are positive constants α¯ k and β¯k such that



max

t ∈[0,τmax ]





¯ k, t αk (t ) /(αk (t ) + 1) = α

max

t ∈[0,τmax ]



t βk (t ) /(βk (t ) + 1) = β¯k ,

k = 1, 2, . . . , m .

(6.4)

ˆ m

ˆ m

Theorem 6.3. Suppose that in Eq. (1.1), gk ∈ H ω 1 (0, τmax ) and uk ∈ H ω 2 (0, τmax ) for k = 1, 2, . . . , m. Let g˜ k (t ) be the shifted CCFs expansion of gk (t ) and u˜ k (t ) be the approximate solution of Eq. (1.1) obtained by the presented method. Then, we have

  uk − u˜ k 

L ∞ (0,τmax )

where

 Bk = k1 + α¯ k β¯k

≤ Bk ,

m 

(6.5)



σkl l2

μk1 +

l =1

m 



 + α¯ k β¯k

ρkl ul L ∞ (0,τmax )

l =1



μk2 +

m 



ρkl l2

l =1 m 



σkl l2 + ul L ∞ (0,τmax )



 ,

(6.6)

l =1

in which

k1 = mˆ 1 n1/2−mˆ 1 k2 = mˆ 2 n

ˆ2 1/2−m



1/2

2

τmax 2

| gk | 1/2 |uk |

τmax

ˆ ;n m

H w1 (0,τmax )

ˆ ;n m

H w2 (0,τmax )

, (6.7)

.

Proof. By applying the presented method, we have for k = 1, 2, . . . , m

⎛

u˜ k (t )  g˜ k (t )+

1

(αk (t )) 1

(βk (t ))

t

⎝ (t − τ )αk (t )−1 Lk1

⎞  t , τ , u˜ 1 (τ ), u˜ 2 (τ ), . . . , u˜ m (τ ) dτ ⎠



0

⎞ ⎛ t    ⎝ (t − τ )βk (t )−1 Lk2 t , τ , u˜ 1 (τ ), u˜ 2 (τ ), . . . , u˜ m (τ ) dτ ⎠ . 0

Subtracting the above equation from Eq. (1.1) results in

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

uk (t ) − u˜ k (t )  gk (t ) − g˜ k (t ) +

1

(βk (t )) +

0

0

(βk (t ))

(αk (t ))

(αk (t ))

⎞ ⎛ t  α ( t )− 1 ⎝ (t − τ ) k Lk1 (t , τ , u 1 (τ ), . . . , um (τ )) dτ ⎠

⎞ ⎛ t     ⎝ (t − τ )βk (t )−1 Lk2 (t , τ , u 1 (τ ), . . . , um (τ )) − Lk2 t , τ , u˜ 1 (τ ), . . . , u˜ m (τ ) dτ ⎠

1

1

1

⎞ ⎛ t    ⎝ (t − τ )βk (t )−1 Lk2 t , τ , u˜ 1 (τ ), . . . , u˜ m (τ ) dτ ⎠ 0

⎞ ⎛ t     ⎝ (t − τ )αk (t )−1 Lk1 (t , τ , u 1 (τ ), . . . , um (τ )) − Lk1 t , τ , u˜ 1 (τ ), . . . , u˜ m (τ ) dτ ⎠ . 0

Thus, we have

    uk (t ) − u˜ k (t )   gk (t ) − g˜ k (t ) + 1

(βk (t )) +

171

0

0

(βk (t ))

(αk (t ))

(αk (t ))

⎞ ⎛ t  α ( t )− 1 ⎝ (t − τ ) k |Lk1 (t , τ , u 1 (τ ), . . . , um (τ ))| dτ ⎠

⎞ ⎛ t     β ( t )− 1 Lk2 (t , τ , u 1 (τ ), . . . , um (τ )) − Lk2 t , τ , u˜ 1 (τ ), . . . , u˜ m (τ )  dτ ⎠ ⎝ (t − τ ) k

1

1

1

⎞ ⎛ t     ⎝ (t − τ )βk (t )−1 Lk2 t , τ , u˜ 1 (τ ), . . . , u˜ m (τ )  dτ ⎠ 0

⎞ ⎛ t     ⎝ (t − τ )αk (t )−1 Lk1 (t , τ , u 1 (τ ), . . . , um (τ )) − Lk1 t , τ , u˜ 1 (τ ), . . . , u˜ m (τ )  dτ ⎠ . 0

Now according to Remark 3, the above relation yields

    uk (t ) − u˜ k (t )   gk (t ) − g˜ k (t ) + 1

(βk (t )) +

l =1

0

l =1

0

1

1

(αk (t ))

⎛ t  ⎞   m  α ( t )− 1 ⎝ (t − τ ) k μk1 + ρkl |ul (t )| dτ ⎠

⎛ t m  ⎞     β ( t )− 1 ⎝ (t − τ ) k σkl ul (t ) − u˜ l (t ) dτ ⎠

(βk (t ))

(αk (t ))

1

⎛ t   ⎞  m    β ( t )− 1 ⎝ (t − τ ) k μk2 + σkl u˜ l (t ) dτ ⎠ l =1

0

⎛ t m  ⎞     α ( t )− 1 ⎝ (t − τ ) k ρkl ul (t ) − u˜ l (t ) dτ ⎠ . l =1

0

The above inequality together with Eq. (6.4) result in

  uk − u˜ k 

L ∞ (0,τmax )

    gk − g˜ k  

μk1 + 

μk2 +

L ∞ (0,τmax

m  l =1 m  l =1



+ α¯ k β¯k )

m 

  σkl ul − u˜ l 

l =1





ρkl ul L ∞ (0,τmax ) + α¯ k β¯k   σkl ul − u˜ l 

m 

 L ∞ (0,τmax )

  ρkl ul − u˜ l 

l =1 L ∞ (0,τmax )

+ ul L ∞ (0,τmax )



 L ∞ (0,τmax )

 (6.8)

.

Now according to inequalities expressed in Eqs. (6.8) and (6.2), the proof is complete.



172

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

Corollary 6.4. Let U (t ) = [u 1 (t ) u 2 (t ) . . . um (t )] T be the analytical solution of the system introduced in Eq. (1.1), and U (t ) =



T

u˜ 1 (t ) u˜ 2 (t ) . . . u˜ m (t ) be the approximate solution obtained by the proposed method, and assume the assumption of Theorem 6.3 is satisfied. Then, we have

  U − U  ∞ L (0,τ

max )

 max {B1 , B2 , . . . , Bm } ,

(6.9)

where Bk for k = 1, 2, . . . , m is defined in Eq. (6.6). Proof. According to the properties of infinity-norm, the proof is clear.



7. Numerical examples Herein, some test problems are investigated to show the precision of the presented method. All the computations are done via Maple 17 with 65 digits accuracy on a 3.5 GHz Core i5 personal computer with 8 GB of RAM. Meanwhile, the convergence order (C-O) of the presented method is computed by

 

CO = log

 !

ε2  ε1 

log

n2 + 1 n1 + 1

 ,

where ε1 and ε2 are the first and the second values of the maximum absolute error (MA-E) appeared by the proposed method, respectively. Moreover, (ni + 1) for i = 1, 2 is the number of the shifted CCFs applied in the first and the second approximations, respectively. Example 1. First, consider the nonlinear variable-order fractional QIE [11]

⎛ u (t ) = g (t ) + ⎝

t

1

 (α (t ))

⎞⎛ (t − τ )α (t )−1 et −τ u 2 (τ )dτ ⎠ ⎝

1

where

β(t )−1

(t − τ )

 (β(t ))

0

⎞

t sin(t )e

u (τ )dτ ⎠ ,

−4 τ 3

0





g (t ) = et 1 − sin(t )t α (t )+β(t ) E1,α (t )+1 (t )E1,β(t )+1 (−t ) . The analytical solution of this integral equation is u (t ) = et . This example is solved on the domain [0, 1] by using the presented method with some values for n, where

α (t ) = 1 + 2.5e−t ,

β(t ) = 1 + t + 3 sin(t ).

A comparison between the numerical results obtained by the presented method and those reported in [11] is given in Table 1. The results of this table indicate that the proposed method is more accurate than the one used in our previous work [11]. Moreover, whenever the number of the shifted CCFs increases, the precision of the numerical solution improves. Diagrams of the absolute error (AE) function with two values n = 6 and n = 12 are shown in Fig. 1, which confirm the high accuracy of the presented approach all over the domain. Note that we have computed the number of the basis functions utilized in the method proposed in [11] by n = 2k M where M indicates the order of the CCFs used in definition of the Chebyshev cardinal wavelets. The values of the MA-E and C-O of the proposed method are reported in Table 5. From this Table, it can be seen that by increasing the number of the basis functions, the approximate solutions tend to the exact solution of the problem with high order of convergence. Example 2. Consider the nonlinear variable-order fractional QIE [11]

⎛ u (t ) = g (t )+ ⎝

1

 (α (t )) 0

⎛ ⎝

t

1

t

 (β(t ))

⎞  (t − τ )α (t )−1 sin(t − τ ) sec2 (τ ) 1 − u 2 (τ ) dτ ⎠ 

⎞ (t − τ )β(t )−1 e −t τ 3 arcsin(u (τ ))dτ ⎠ ,

0

where



g (t ) = sin(t ) − e −t sin(t )t α (t ) E2,α (t )+1 (−t )2 − cos(t )t α (t )+1 E2,α (t )+2 (−t )2

 24 t β(t )+4 . (β(t ) + 5)

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

173

Table 1 The AE values of the presented method and the method applied in [11] for Example 1. t

0.1 0.3 0.5 0.7 0.9

The method applied in [11]

The proposed method

n=6

n = 12

n = 14

n = 28

n=6

n=8

n = 10

n = 12

7.5564E-04 4.8192E-04 4.6060E-04 7.7092E-04 1.3333E-03

5.2076E-05 1.0884E-04 2.2185E-05 1.4488E-04 1.1195E-04

1.5305E-09 8.3581E-09 4.2994E-09 2.8275E-07 5.6064E-07

2.7443E-11 5.1160E-10 9.6441E-10 1.3677E-08 1.2088E-07

7.8953E-09 1.0113E-08 2.8438E-13 1.0629E-08 8.7329E-09

2.9529E-11 1.8203E-11 5.0964E-16 1.8946E-11 3.1985E-11

1.3401E-14 1.9123E-14 7.2198E-20 1.9771E-14 1.4323E-14

3.7697E-18 6.2765E-18 5.5619E-23 6.4583E-18 3.9915E-18

Fig. 1. Diagrams of the AE function with n = 6 (left side) and n = 12 (right side) in Example 1.

Table 2 The AE values of the presented method and the method applied in [11] for Example 2. t

0.1 0.3 0.5 0.7 0.9

The method applied in [11]

The proposed method

n=6

n = 12

n = 14

n = 28

n=6

n=8

n = 10

n = 12

5.9390E-04 3.5694E-04 1.7880E-04 2.7031E-04 4.3479E-04

4.5869E-05 7.1328E-05 9.7385E-06 6.1693E-05 3.0831E-05

2.9467E-10 4.8973E-09 1.2538E-09 1.8669E-08 6.9271E-09

1.2282E-12 1.3175E-10 1.8894E-10 2.9257E-10 2.0024E-10

4.4695E-09 5.5225E-08 1.0296E-14 5.3736E-09 4.2315E-09

1.6512E-11 9.8861E-12 1.9282E-18 9.6723E-12 1.5805E-11

7.4321E-15 1.0348E-14 2.0384E-21 1.0161E-14 7.1661E-15

2.0784E-18 3.3881E-18 1.9171E-25 3.3356E-18 2.0145E-18

The analytical solution of this example is u (t ) = sin(t ). The presented method with some values of n is applied for numerical solution of this example on the domain [0, 1], where

α (t ) = 0.75 + 2 sinh(t ),

β(t ) = 0.5 + 2 cosh(t ).

In Table 2, a numerical comparison has been performed between the numerical results obtained by the presented method and the ones obtained [11]. The information of this table confirms that the proposed method is more accurate than the one presented in [11]. Moreover, it can be seen that the high accurate results are achieved for even a few number of the shifted CCFs. Furthermore, as the number of the basis functions increases, accuracy of the method improves rapidly. Diagrams of the AE function with two values n = 6 and n = 12 are shown in Fig. 2 which express the high and global accuracy of the presented method. The values of the MA-E and C-O of the presented scheme are listed in Table 5. The reported results confirm that the approximate solutions tend to the exact solution of the problem with high order of convergence. Example 3. Consider the following system of nonlinear variable-order fractional QIEs:

⎛ u 1 (t ) = g 1 (t )+ ⎝

1

 (α1 (t ))

⎛ ⎝

1

 (β1 (t ))

t 0

t 0

⎞   (t − τ )α1 (t )−1 et −τ u 21 (τ ) + cos2 (τ ) u 22 (τ ) + sin2 (τ ) dτ ⎠ 

⎞ (t − τ )β1 (t )−1 e −t u 1 (τ )u 2 (τ )dτ ⎠ ,

174

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

Fig. 2. Diagrams of the AE function with n = 6 (left side) and n = 12 (right side) in Example 2. Table 3 The AE values of the presented method with different values of n in Example 3. t

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 1.7 1.9

u 1 (t )

u 2 (t )

n=5

n=7

n=9

n = 11

n=5

n=7

n=9

n = 11

1.2061E-04 9.0686E-05 6.9763E-05 9.8526E-05 6.2553E-07 1.1755E-04 1.5455E-04 9.5510E-05 1.8580E-05 4.7852E-05

1.6368E-06 1.2959E-06 4.5058E-07 1.8326E-06 1.4524E-06 4.2531E-07 4.9428E-07 1.5339E-06 1.8870E-06 3.0241E-06

6.5733E-09 6.2941E-09 1.6882E-08 3.7725E-09 8.3734E-09 1.4795E-08 1.1023E-08 6.9793E-09 1.3984E-08 4.4058E-10

3.0079E-11 8.5431E-11 4.3905E-11 4.9492E-12 1.0036E-10 5.7386E-11 3.9897E-11 7.8675E-11 3.3516E-11 5.6748E-11

3.2156E-05 1.0344E-05 2.0942E-05 6.3792E-06 2.5030E-05 4.1639E-05 4.2374E-05 5.0573E-05 7.9078E-05 7.8377E-05

1.6186E-07 1.6282E-07 3.0696E-08 2.1007E-07 3.8112E-07 5.8160E-07 6.9854E-07 5.7795E-07 5.4473E-07 8.1365E-07

2.4669E-10 4.2703E-10 1.4328E-09 1.6696E-09 9.8037E-10 1.7841E-09 5.0212E-09 6.8575E-09 4.0301E-09 5.0250E-09

1.0289E-12 3.5676E-12 2.8621E-12 1.7385E-12 1.4086E-11 2.2721E-11 1.4431E-11 2.7736E-11 3.2765E-11 2.0185E-11

⎛ u 2 (t ) = g 2 (t )+ ⎝

1

 (α2 (t ))

⎛ ⎝

1

 (β2 (t ))

t

⎞ (t − τ )α2 (t )−1 e cos(τ )−sin(τ ) sin(t ) cos(τ )e u 1 (τ )−u 2 (τ ) dτ ⎠

0

t

⎞ (t − τ )β2 (t )−1 (sin(t )u 1 (τ ) + cos(t )u 2 (τ )) dτ ⎠ ,

0

where

g 1 (t ) = sin(t ) − t α1 (t )+β1 (t )+1 E1,α1 (t )+1 (−t ) E2,β1 (t )+2 (−4t 2 ),





g 2 (t ) = cos(t ) − t α2 (t ) sin(t ) E2,α2 (t )+1 (−t 2 ) t β2 (t )+1 sin(t ) E2,β2 (t )+2 (−t 2 ) + t β2 (t ) cos(t ) E2,β2 (t )+1 (−t 2 ) . The exact solution of this system is u 1 (t ), u 2 (t ) = sin(t ), cos(t ) . The method proposed in Section 5 with some values of n is applied for numerical solution of this example on the domain [0, 2], where

α1 (t ) = 0.55 + 0.25 sin(t ), α2 (t ) = 0.45 + 0.35 sin(t ), β1 (t ) = 1 − 0.25e−t , β2 (t ) = 1 − 0.75e−2t . The obtained results are illustrated in Table 3 and Fig. 3, which show that the presented method is highly accurate in solving such system. Meanwhile, using more number of the basis functions increases the accuracy rapidly. It is worth noting that although the investigated system is a weakly singular system of nonlinear variable-order fractional QIEs, the method is successfully applied. The values of the MA-E and C-O of the presented scheme are shown in Table 5. The reported results show that the approximate solutions tend to the exact solution of the problem with high order of convergence. Example 4. Consider the following system of nonlinear variable-order fractional QIEs:

⎛

u 1 (t ) = g 1 (t )+ ⎝

1

 (α1 (t ))

t 0

⎞

(t − τ )α1 (t )−1 et −τ e τ u 1 (τ )u 2 (τ ) u 3 (τ )dτ ⎠

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

175

Fig. 3. Diagrams of the AE function with n = 11 for u 1 (t ) (left side) and u 2 (t ) (right side) in Example 3.

⎛ ⎝

t

1

 (β1 (t ))

⎛ u 2 (t ) = g 2 (t )+ ⎝

 (α2 (t ))

⎛ ⎝

 (β2 (t ))

⎛ u 3 (t ) = g 3 (t )+ ⎝

 (α3 (t ))

⎛ ⎝

(t − τ )α2 (t )−1 sin(t − τ )e −tu 1 (τ )u 2 (τ ) dτ ⎠ ⎞ (t − τ )β2 (t )−1t 3 τ 2 u 1 (τ )u 2 (τ ) arcsin(u 3 (τ ))dτ ⎠ ,

0

t

1

⎞

0

t

1

(t − τ )β1 (t )−1 cos(t − τ )u 1 (τ )u 2 (τ )dτ ⎠ ,

0

t

1

⎞

⎞ (t − τ )α3 (t )−1 et −sin(τ ) u 1 (τ )u 2 (τ )e u 3 (τ ) dτ ⎠

0

t

1

 (β3 (t ))

where

(t − τ )β3 (t )−1 0

e −(t +τ ) 2 + sin(τ )

⎞ ln (1 + u 1 (τ )u 2 (τ )) (u 3 (τ ) + 2) dτ ⎠ ,





g 1 (t ) = e −t − t α1 (t )β1 (t )+1 et E2,α1 (t )+2 (−t 2 ) cos(t )E2,β1 (t )+1 (−t 2 ) + t sin(t )E2,β1 (t )+2 (−t 2 ) ,



6

g 2 (t ) = et −



t α2 (t )+β2 (t )+6 e −t sin(t )E2,α2 (t )+1 (−t 2 ) − t cos(t )E2,α2 (t )+2 (−t 2 ) ,

(β2 (t ) + 4) ln(2) g 3 (t ) = sin(t ) − t α3 (t )+β3 (t ) E1,β3 (t )+1 (−t ). (α3 (t ) + 1) 



The exact solution of this system is u 1 (t ), u 2 (t ), u 3 (t ) = e −t , et , sin(t ) . The presented method with some values of n is applied for numerical solution of this example on the domain [0, 1.5], where

α1 (t ) = 3.25 + 0.25 sin(t ), α2 (t ) = 2.55 + 0.55 sin(t ), α3 (t ) = 1.75 + 1.35 sin(t ), β1 (t ) = 0.65 − 0.25e−t ,

β2 (t ) = 1.55 − 0.65e−2t ,

β3 (t ) = 2.75 − 0.25e−3t .

The numerical results are shown in Table 4 and Figs. 4–6 which expresses the high accuracy of the presented approach. The values of the MA-E and C-O of the established method are given in Table 6. The reported results declare that the approximate solutions tend to the exact solution of the problem with high order of convergence. Example 5. Finally, consider the nonlinear variable-order fractional QIE

⎛ u (t ) = g (t ) + ⎝

1

t

 (α (t )) 0

⎞⎛ (t − τ )α (t )−1 e −t τ 2 u 2 (τ )dτ ⎠ ⎝

1

t

 (β(t )) 0

⎞ √ (t − τ )β(t )−1 sin(t ) τ u 3 (τ )dτ ⎠ ,

176

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

Table 4 The AE values of the presented method with different values of n in Example 4. t

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5

u 1 (t )

u 2 (t )

u 3 (t )

n=7

n=9

n = 11

n=7

n=9

n = 11

n=7

n=9

n = 11

5.1344E-09 4.1645E-09 8.8643E-09 8.1721E-09 3.5691E-10 7.3253E-09 8.1551E-09 5.3939E-09

7.4533E-12 1.4752E-11 1.4131E-11 1.1171E-11 6.1977E-12 2.1758E-12 6.5034E-12 1.7212E-11

1.6317E-14 2.1152E-15 9.4241E-15 1.0891E-14 1.1636E-14 1.3475E-14 1.4203E-14 1.1451E-14

1.9728E-08 1.6709E-08 3.6894E-08 3.5602E-08 1.6933E-09 3.1880E-08 4.2175E-08 4.4947E-08

2.9928E-11 6.1583E-11 6.1332E-11 5.0716E-11 2.8171E-11 9.5286E-12 2.5256E-11 6.9959E-11

6.5694E-14 8.7745E-15 4.0117E-14 4.7855E-14 5.2229E-14 6.3613E-14 6.4768E-14 7.2872E-14

6.0990E-09 5.2010E-09 1.1543E-08 1.1170E-08 5.4194E-10 1.0025E-08 1.3206E-08 1.4065E-08

9.3275E-12 1.9287E-11 1.9280E-11 1.5976E-11 8.8907E-12 2.9910E-12 7.9576E-12 2.1884E-11

2.0584E-14 2.7594E-15 1.2653E-14 1.5118E-14 1.6524E-14 2.0111E-14 2.0482E-14 2.3036E-14

Fig. 4. Diagrams of the AE function for u 1 (t ) with n = 7 (left side) and n = 11 (right side) in Example 4.

Fig. 5. Diagrams of the AE function for u 2 (t ) with n = 7 (left side) and n = 11 (right side) in Example 4.

Table 5 The values of the MA-E and C-O of the proposed method in Examples 1-3. n

6 7 8 9 10 11

Example 1

Example 2

u (t )

u (t )

Example 3 u 1 (t )

u 2 (t )

MA-E

C-O

MA-E

C-O

MA-E

C-O

MA-E

C-O

4.2819E-08 1.3271E-09 3.5146E-11 9.1155E-13 1.9899E-14 4.1200E-16

– 26.0162 30.8299 34.6630 40.1266 44.5620

2.1547E-08 3.9650E-10 1.8442E-11 2.7117E-13 1.0470E-14 1.1882E-16

– 29.9204 26.0482 40.0495 34.1435 51.4721

7.2172E-06 2.0489E-06 9.5564E-08 1.8361E-08 7.1028E-10 1.0922E-10

– 09.4297 26.0246 15.6564 34.1235 21.5178

2.9534E-06 6.9856E-07 2.3534E-08 6.9065E-08 1.2241E-10 3.7165E-11

– 10.7966 28.7866 10.2182 66.4716 13.6994

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

177

Fig. 6. Diagrams of the AE function for u 3 (t ) with n = 7 (left side) and n = 11 (right side) in Example 4.

Fig. 7. Diagrams of the AE function with n = 6 (left side) and n = 11 (right side) in Example 5. Table 6 The values of the MA-E and C-O of the proposed method in Examples 4 and 5. n

Example 4

Example 5

u 1 (t )

6 7 8 9 10 11

u 2 (t )

u 3 (t )

u (t )

MA-E

C-O

MA-E

C-O

MA-E

C-O

MA-E

C-O

2.1088E-07 9.9610E-08 3.8943E-10 1.7212E-11 5.0207E-13 1.5662E-14

– 05.6168 47.0724 29.6038 37.0854 39.8511

9.6662E-07 4.4947E-08 1.7536E-09 6.9959E-11 2.2290E-12 6.8995E-14

– 22.9782 27.5405 30.5761 36.1593 39.9404

3.1401E-07 1.4065E-08 5.9175E-10 2.1884E-11 7.5479E-13 2.1803E-14

– 23.2584 26.8999 31.2956 35.3275 40.7348

9.9314E-03 7.7857E-03 4.0455E-03 2.7104E-03 2.4461E-03 4.6684E-04

– 01.8228 05.5583 03.8013 01.0764 19.0350

where

g (t ) =

√

t −e

 −t

sin(t )

6 t 3+α (t )



(α (t ) + 4)

The exact solution of this example is u (t ) = [0, 1], where

α (t ) = 0.75 + 2 sinh(t ),

2 t 2+β(t )

(β(t ) + 3)

 .

√

t. The method expressed in Section 5 is used for this example on the domain

β(t ) = 0.5 + 2 cosh(t ).

Diagrams of the AE function with two values n = 6 and n = 11 are shown in Fig. 7. Also, the values of the MA-E and C-O of the proposed method are reported in Table 6. From Fig. 7 and Table 6, one can see that applying more terms of the basis functions can not provide numerical results with high accuracy. Meanwhile, from Table 6, it can be observed that by increasing the number of the shifted CCFs the approximate solutions tend to the exact solution of the problem with

178

M.H. Heydari / Applied Numerical Mathematics 153 (2020) 164–178

low order of convergence. These limitations are occurred due to the fact that solution of the problem is not sufficiently differentiable. 8. Conclusion In this study, a new operational matrix (OM) method based on the shifted Chebyshev cardinal functions (CCFs) has been developed and applied for numerical solution of a novel class of systems of nonlinear variable-order fractional quadratic integral equations (QIEs). In this way, a new OM of variable-order fractional integration has been extracted for the shifted CCFs and used to establish the presented method. The method transforms such systems into equivalent systems of algebraic equations. The convergence and error analysis of the established method have been investigated theoretically and numerically. Several test problems have been provided to show the accuracy of the presented approach. The achieved results confirm that only a few number of the shifted CCFs are necessary to obtain a high accurate solution for such systems. References [1] C. Canuto, M. Hussaini, A. Quarteroni, T. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1988. [2] Y. Chen, L. Liu, B. Li, Y. Sun, Numerical solution for the variable order linear cable equation with Bernstein polynomials, Appl. Math. Comput. 238 (2014) 329–341. [3] K. Deimling, Nonlinear Functional Analysis, Springer-Verlag, Berlin, 1985. [4] A.A. El-Sayed, P. Agarwal, Numerical solution of multiterm variable-order fractional differential equations via shifted Legendre polynomials, Math. Methods Appl. Sci. 4 (11) (2019) 3978–3991. [5] A.M.A. El-Sayed, M.M. Saleh, E.A.A. Ziada, Numerical and analytic solution for nonlinear quadratic integral equations, Math. Sci. Res. J. 12 (8) (2008) 183–191. [6] A.M.A. El-Sayed, H.H.G. Hashem, E.A.A. Ziada, Picard and Adomian methods for quadratic integral equation, J. Comput. Appl. Math. 29 (2010) 447–463. [7] A.M.A. El-Sayed, M.Sh. Mohamed, F.F.S. Mohamed, Existence of positive continuous solution of a quadratic integral equation of fractional orders, Math. Sci. Lett. 1 (9) (2011) 1–7. [8] A.M.A. El-Sayed, H.H.G. Hashem, Y.M.Y. Omar, Positive continuous solution of a quadratic integral equation of fractional orders, Math. Sci. Lett. 2 (1) (2013) 19–27. [9] R.M. Evans, U.N. Katugampola, D.A. Edwards, Applications of fractional calculus in solving Abel-type integral equations: surface-volume reaction problem, J. Comput. Appl. Math. 73 (2017) 1–17. [10] M.H. Heydari, A new direct method based on the Chebyshev cardinal functions for variable-order fractional optimal control problems, J. Franklin Inst. 355 (2018) 4970–4995. [11] M.H. Heydari, Chebyshev cardinal wavelets for nonlinear variable-order fractional quadratic integral equations, Appl. Numer. Math. 144 (2019) 190–203. [12] M.H. Heydari, M.R. Mahmoudi, A. Shakiba, Z. Avazzadeh, Chebyshev cardinal wavelets and their application in solving nonlinear stochastic differential equations with fractional Brownian motion, Commun. Nonlinear Sci. Numer. Simul. 64 (2018) 98–121. [13] M.H. Heydari, Z. Avazzadeh, M.R. Mahmoudi, Chebyshev cardinal wavelets for nonlinear stochastic differential equations driven with variable-order fractional Brownian motion, Chaos Solitons Fractals 124 (2019) 105–124. [14] M.H. Heydari, Z. Avazzadeh, Y. Yang, A computational method for solving variable-order fractional nonlinear diffusion-wave equation, Appl. Math. Comput. 352 (2019) 235–248. [15] M.H. Heydari, A. Atangana, Z. Avazzadeh, Y. Yang, Numerical treatment of the strongly coupled nonlinear fractal-fractional Schrödinger equations through the shifted Chebyshev cardinal functions, Alex. Eng. J. (2020), https://doi.org/10.1016/j.aej.2019.12.039, in press. [16] R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific, 2000. [17] F. Khane Keshi, B. Parsa Moghaddam, A. Aghili, A numerical approach for solving a class of variable-order fractional functional integral equations, J. Comput. Appl. Math. 37 (2018) 4821–4834. [18] X. Ma, C. Huang, Numerical solution of fractional integro-differential equations by a hybrid collocation method, Appl. Math. Comput. 219 (12) (2013) 6750–6760. [19] A.M. Nagy, N.H. Sweilam, A.A. El-Sayed, New operational matrix for solving multi-term variable order fractional differential equations, J. Comput. Nonlinear Dyn. 13 (2018) 11001–11007. [20] H.T.C. Pedro, M.H. Kobayashi, J.M.C. Pereira, C.F.M. Coimbra, Variable order modeling of diffusive-convective effects on the oscillatory flow past a sphere, J. Vib. Control 14 (2008) 1569–1672. [21] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [22] L.E.S. Ramirez, C.F.M. Coimbra, On the variable order dynamics of the nonlinear wake caused by a sedimenting particle, Physica D 240 (2011) 1111–1118. [23] R. Roohi, M.H. Heydari, O. Bavi, H. Emdad, Chebyshev polynomials for generalized Couette flow of fractional Jeffrey nanofluid subjected to several thermochemical effects, Eng. Comput. (2019), https://doi.org/10.1007/s00366-019-00843-9, in press. [24] N.H. Sweilam, M.M. Khader, A Chebyshev pseudo-spectral method for solving fractional order integro-differential equations, ANZIAM J. 51 (2010) 464–475. [25] M. Zayernouri, G.E. Karniadakis, Fractional spectral collocation methods for linear and nonlinear variable order FPdEs, J. Comput. Phys. 80 (1) (2015) 312–338.