# An Operational Matrix Technique Based on Chebyshev Polynomials for Solving Mixed Volterra-Fredholm Delay Integro-Differential Equations of Variable-Order

Authors: Kamal R. Raslan; Khalid K. Ali; Emad M. H. Mohamed; Jihad A. Younis; Mohamed A. Abd El salam
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

In this work, an algorithm for finding numerical solutions of linear fractional delay-integro-differential equations (LFDIDEs) of variable-order (VO) is introduced. The operational matrices are used as discretization technique based on shifted Chebyshev polynomials (SCPs) of the first kind with the spectral collocation method. The proposed VO-LFDIDEs have multiterms of integer, fractional-order derivatives for delayed or nondelayed and mixed Volterra-Fredholm integral terms. The introduced model is a more general form of linear fractional VO pantograph, neutral, and mixed Fredholm–Volterra integro-differential equations with delay arguments. Caputo’s VO fractional derivative operator is used to generate the matrices of the derivative. Operational matrices are presented for all terms. The reliability and efficiency of the proposed scheme are demonstrated by some numerical experiments. Also, some examples are included to improve the validity and applicability of the techniques. Finally, comparisons between the proposed method and other methods were used to solve this kind of equation.


---

## Body

## 1. Introduction

In recent few decades, fractional calculus has shown to be a useful implement to formulate many problems in science and engineering where the fractional derivatives and integrals can be used for the description of the properties of various real materials in different branches of science [1–7]. In 1993, Samko and Ross [8] have introduced the variable order derivative (VOD) operator just as a generalization of the fractional-order derivative, and they provided some of its main properties. In this introduced operator, the order of the derivative is a function of the independent variables such as time and space variables. With this generalization, many applications have been established in mechanics, physics, signal processing, and control [9–12]. One of the most important tools in studying the prime numbers is the Riemann zeta-function which embodies both additive and multiplicative structures in a single function, the fractional generalization of it is also found in [13–17].Consequently, a generalized kind of differential equation appears named VO fractional differential equations (VOFDEs) [18–24]. In VOFDEs, the differential operator is VO, and the derivative changes in general concerning the independent variable or concerning an external functional behavior [8, 25]. Analytical solutions for the VOFDEs are difficult to obtain because the kernel of the VO’s operator has a variable exponent, hence, there will be increasingly rapid developments in numerical approaches to VOFDEs [26–29]. The delay, neutral delay FDEs, and fractional integro-differential equations (FIDEs) are considered as a generalization and development of FDEs, and dealing with them analytically in most of the cases is difficult [30–34]. The VO fractional delay differential equations (VOFDDEs) are a kind of generalization for the fractional delay differential equations (FDDEs) [35–38]. VOFDDEs did not receive the attention that the FDDEs accomplished; however, the potential to describe complicated behavior by the VO of differentiation or integration is clear, all this was a motivation to start studying this type of equation. Soon after, a variety of definitions have been offered for variable order integral and derivative operators such as Riemann–Liouville (RL) [39], Grünwald–Letnikov [40], Caputo, and Fabrizio [41] derivatives.In this work, we use operational matrices discretization technique based on shifted Chebyshev polynomials (SCPs) of the first kind with the spectral collocation method to solve the following type of linear fractional delay-integro-differential equations of variable-order (LFDIDEs):(1)∑i=0n1QixDνixyx+∑j=0n2PjxDαjxyx−τj+∑s=0n3Rsxysx−εs=gx+∫abKx,tytdt+∫axVx,tytdt,x∈a,b,under the conditions
(2)yiρi=μi,i=0,1,2⋯,m−1,where 0<αjx≤1, 0<νix≤1, and τj,εs are positive real delay arguments, and the known functions Qix,Pjx, Rsx, Kx,t, Vx,t, and gx are well defined, additionally, yx is an unknown function to be determined. The symbole D denotes the variable-order derivative (VOD) operator in the Caputo sense.Corollary 1.
The independent variablex of (1) belongs to a,b, which is the intersection of the intervals of the different delayed arguments and 0,h, i.e., x∈a,b=τj,h+τj⋂εs,h+εs⋂0,h, and τj,εs<h.The introduced model (1) is a more general form of linear fractional VO pantograph, neutral, and mixed Fredholm–Volterra integro-differential equations with delay arguments. Several methods have been presented for solving VO differential equations [42–46]. For example, Ganji et al. applied the first-kind CPs to obtain the solution of variable order differential equations [47]. Doha et al. solved variable-order fractional Volterra integro-differential equations by shifted Legendre–Gauss–Lobatto collocation method [39]. A numerical method based on the Jacobi polynomials for solving variable-order differential equations has been proposed in [48]. Furthermore, Ganji et al. have introduced a numerical scheme based on the Bernstein polynomials to solve variable order diffusion-wave equations [49]. As a notation, the proposed model (1) is a generalization of our previous reports [50–54] and other work [55–57] as well.


## 2. Basic Definitions

In this preliminary section, the definitions and properties of Caputo’s VO fractional derivative and the required mathematical tools will be introduced. Also, briefly, the SCPs of the first kind, we will present what we need from their properties.


### 2.1. Caputo Variable-Order’s Operator

Definition 2.
The Caputo’sVO fractional derivative operator of order n−1<αx≤n operates a function gx∈Cn0,h in [42, 43] as
(3)Dαxgx=1Γn−αx∫0xgsnx−sαx−n+1ds,forn−1<αx<n,gnx,forαx=n,and n∈ℵ.Remark 3.
In this work, we write the Variable-Order Caputo’s fractional symbolDαx instead of 0CDxαx for short.Remark 4.
The following useful property according to Definition2:
(4)Dαxxm=Γm+1Γm+1−αxxm−αx,form≥n,0,otherwise.where
(5)n−1<αx≤n,m∈ℵ0andℵ0=ℵ∪0.Definition 5.
The Caputo’s VO fractional derivative operator of orderαx,0<αx≤1; operates gx, where gx∈Cn0,h can be defined as [42, 58]
(6)Dαxgx=1Γ1−αx∫0xg′xx−tαxdt+g0+−g0−Γ1−αxx−αx,and 0<αx≤1. Definition 5 can be reformulated if we assume that; the starting time in a perfect situation, then:
(7)Dαxgx=1Γ1−αx∫0xg′xx−tαxdt,0<αx<1.Remark 6.
By (7), we can get the following formulas:
(i)
The Caputo’s VO operator is linear(8)Dαxμfx+λgx=μDαxfx+λDαxgx,where μ and λ are any two constants.
(ii)
As the ordinary differentiation operator the Caputo’s VO operator as well, such thatDαxK=0, and K is constant.


### 2.2. Shifted First Kind of Chebyshev Polynomials

Next, let us introduce some properties of the SCPs [59]. It is well known that the classical CPs are defined on −1,1 by the three-term recurrence relation:
(9)Tn+1t=2tTnt−Tn−1t,T0t=1,T1t=t,n=1,2,⋯.Let the SCPsTn2t/h−1 be denoted by Tnhx, which are defined on 0,h. Then, Tnhx can be generated by using the following recurrence relation:
(10)Tn+1hx=22xh−1Tnhx−Tn−1hx,n=1,2,⋯,where
(11)T0hx=1,T1hx=2xh−1.The analytic form of the SCPsTih of degree i is given by:
(12)Tihx=∑k=0i−1i−jii+j−1!22ji−j!2j!hjxk.The orthogonality condition of the SCPs is(13)∫0hTihxTjhxwhdx=δj,where
(14)wh=1/hx−x2,δj=0,fori≠j,ϕj/2π,fori=jϕ0=2,ϕj=1,j≥1.A functionyx, square integrable in 0,h, may be expressed in terms of shifted Chebyshev polynomials as
(15)yx=∑i=0∞aiTihx,where the coefficients ai are given by
(16)ai=1hi∫0hyxTihxwhdx,i=0,1,2,⋯.In practice, only the first1+N-terms shifted Chebyshev polynomials are considered. The special values
(17)Ti,jh=−1i−jii+j−1!22ji−j!2j!hj,j≤i,will be of important use later.In the approximation theory, the series in (15) can be approximated by taking the first N+1 terms as follows:
(18)yx≃yNx=∑i=0NaiTihx=ATThx,where A=a0,a1,⋯,aNT is a vector, Thx=T0hxT1hx⋯TNhxT and if we assume that Xx=x0x1⋯xN.From (12) and (17), Thx can be written as the following form:
(19)Thx=XxMT,where M is square lower triangle matrix with size N+1×N+1 given by:
(20)M=m0000⋯0m10m110⋯0m20m21m22⋯0⋮⋮⋮⋮⋮mN+1,0mN+1,1mN+1,2⋯mN+1,N+1.and its entries elements are given by:
(21)mij=1,ifi=j=0,−1i−jii+j−1!22ji−j!2j!hj,ifj≤i,0,ifj>i.For example, ifN=4, then the square matrix M is given by:
(22)M=10000−120001−8800−118−483201−32160−256128.Now, from (19), we can obtain the kth derivative of the matrix Tx as
(23)Tkx=XkxMT,k=0,1,2,⋯.


## 2.1. Caputo Variable-Order’s Operator

Definition 2.
The Caputo’sVO fractional derivative operator of order n−1<αx≤n operates a function gx∈Cn0,h in [42, 43] as
(3)Dαxgx=1Γn−αx∫0xgsnx−sαx−n+1ds,forn−1<αx<n,gnx,forαx=n,and n∈ℵ.Remark 3.
In this work, we write the Variable-Order Caputo’s fractional symbolDαx instead of 0CDxαx for short.Remark 4.
The following useful property according to Definition2:
(4)Dαxxm=Γm+1Γm+1−αxxm−αx,form≥n,0,otherwise.where
(5)n−1<αx≤n,m∈ℵ0andℵ0=ℵ∪0.Definition 5.
The Caputo’s VO fractional derivative operator of orderαx,0<αx≤1; operates gx, where gx∈Cn0,h can be defined as [42, 58]
(6)Dαxgx=1Γ1−αx∫0xg′xx−tαxdt+g0+−g0−Γ1−αxx−αx,and 0<αx≤1. Definition 5 can be reformulated if we assume that; the starting time in a perfect situation, then:
(7)Dαxgx=1Γ1−αx∫0xg′xx−tαxdt,0<αx<1.Remark 6.
By (7), we can get the following formulas:
(i)
The Caputo’s VO operator is linear(8)Dαxμfx+λgx=μDαxfx+λDαxgx,where μ and λ are any two constants.
(ii)
As the ordinary differentiation operator the Caputo’s VO operator as well, such thatDαxK=0, and K is constant.


## 2.2. Shifted First Kind of Chebyshev Polynomials

Next, let us introduce some properties of the SCPs [59]. It is well known that the classical CPs are defined on −1,1 by the three-term recurrence relation:
(9)Tn+1t=2tTnt−Tn−1t,T0t=1,T1t=t,n=1,2,⋯.Let the SCPsTn2t/h−1 be denoted by Tnhx, which are defined on 0,h. Then, Tnhx can be generated by using the following recurrence relation:
(10)Tn+1hx=22xh−1Tnhx−Tn−1hx,n=1,2,⋯,where
(11)T0hx=1,T1hx=2xh−1.The analytic form of the SCPsTih of degree i is given by:
(12)Tihx=∑k=0i−1i−jii+j−1!22ji−j!2j!hjxk.The orthogonality condition of the SCPs is(13)∫0hTihxTjhxwhdx=δj,where
(14)wh=1/hx−x2,δj=0,fori≠j,ϕj/2π,fori=jϕ0=2,ϕj=1,j≥1.A functionyx, square integrable in 0,h, may be expressed in terms of shifted Chebyshev polynomials as
(15)yx=∑i=0∞aiTihx,where the coefficients ai are given by
(16)ai=1hi∫0hyxTihxwhdx,i=0,1,2,⋯.In practice, only the first1+N-terms shifted Chebyshev polynomials are considered. The special values
(17)Ti,jh=−1i−jii+j−1!22ji−j!2j!hj,j≤i,will be of important use later.In the approximation theory, the series in (15) can be approximated by taking the first N+1 terms as follows:
(18)yx≃yNx=∑i=0NaiTihx=ATThx,where A=a0,a1,⋯,aNT is a vector, Thx=T0hxT1hx⋯TNhxT and if we assume that Xx=x0x1⋯xN.From (12) and (17), Thx can be written as the following form:
(19)Thx=XxMT,where M is square lower triangle matrix with size N+1×N+1 given by:
(20)M=m0000⋯0m10m110⋯0m20m21m22⋯0⋮⋮⋮⋮⋮mN+1,0mN+1,1mN+1,2⋯mN+1,N+1.and its entries elements are given by:
(21)mij=1,ifi=j=0,−1i−jii+j−1!22ji−j!2j!hj,ifj≤i,0,ifj>i.For example, ifN=4, then the square matrix M is given by:
(22)M=10000−120001−8800−118−483201−32160−256128.Now, from (19), we can obtain the kth derivative of the matrix Tx as
(23)Tkx=XkxMT,k=0,1,2,⋯.


## 3. Operational Matrices

In this section, we introduce the generalized operational matrices forThkx, Thsx−τs, DνixThx, and DαixThx−τi according to fractional calculus.Lemma 7.
Thekth order derivative of the row vector Thx can be in the following relation form [50, 60]:
(24)Thkx=XxBkMT,where B is the N+1×N+1 operational matrix of derivative for Xx and can be obtained from
(25)B=εij=0,forotherwise,j+1,fori=j+1,,j=0,1,⋯,N.Proof (see [50, 60]).
The row vectorTx−τ represents in terms of the vector Xx in the following form [50, 60]:
(26)Tx−τ=XxB−τMT,where B−τ is square upper triangle matrix with size N+1×N+1 given by:
(27)B−τ=βij=ji−τj−i,forj≥i,0,forotherwise,,i,j=0,1,⋯,N.Corollary 8.
Thesth order derivative of the delay row vector Thx−τs can be represented as [50, 60]
(28)Thsx−τ=Xsx−τsMT=XxBsB−τsMT.According to the previous lemmas with the fractional calculus by using the Caputo’s variable-order fractional derivative, we introduce the following theorem.Theorem 9.
Theνixth variable-order fractional derivative of the vector Tx can be written as
(29)DνixThx=XνixxBνixxMT,where
(30)Xνixx=0,x1−νix,x2−νix,⋯,xN−νix,0<νix<1,where Bνixx diagonal matrix with size N+1×N+1, its elements λij given by:
(31)λij=Γi+1Γi+1−νix,fori=j≠0.,0,forotherwise.,i=0,1,⋯,N.Proof.
(32)DνixThx=DνixXxMT=Dνix1,x,x2,⋯,xNMT,0<νix<1,=0,Γ2Γ2−νixx1−νix,Γ3Γ3−νixx2−νix,⋯,ΓN+1ΓN+1−νixxN−νixMT=XνixxBνixMT.Corollary 10.
Relation (29) satisfies for n−1<νix<n according to Definition 2 by induction with
(33)Xνixx=0,0,0,⋯,xn−νix,⋯,xN−νix,n−1<νix<n,and
(34)Bνixx=λij=Γi+1Γi+1−νix,fori=j≥n,0,forotherwise.,i=0,1,⋯,N,n<N.Corollary 11.
Theαixth variable-order fractional derivative of the delay vector Tx−τi can be written in the following form:
(35)DαixThx−τi=XαixxBαixxB−τiMT.
Form (18), we get
(36)yNx=ThxA,(37)yNkx=ThkxA,(38)DνixyNx=ThνixxA,and
(39)DαixyNx−τi=Tαixx−τiA,i=0,1,⋯,N.
By puttingx⟶x−εs in the relation (37), we obtain the matrix form
(40)yNsx−εs=Thsx−εsA.
Consequently, by substituting the matrix form (24) into (37), we have the matrix relation
(41)yNkx=XkxBkMTA.
By substituting the matrix forms (29) into (38), we have the matrix relation
(42)DνixyNx=XνixxBνixxMTA.
By using (35) and (39), we get
(43)DαjxyNx−τj=XαjxxBαjxxB−τjMTA.
Similar form (28) and (40), we obtain
(44)yNsx−εs=XxBsB−εsMTA.


### 3.1. Matrix Representation for Fredholm Integral Terms

Now, we try to find the matrix form corresponding to the integral term. Assume thatKx,t can be expanded to univariate Chebyshev series concerning t, as follows
(45)Kx,t≅∑r=0NurxTrt.Then the matrix representation of the kernel functionKx,t is given by
(46)Kx,t≅UxTTt,where
(47)Ux=u0xu1x⋯uNx.Substituting the relations (29) and (46) in the present integral part, we obtained
(48)∫abKx,tytdt=∫abUxTTtTtAdt=∫abUxMXTtXtMTAdt=∫abUxMt0t1⋯tNTt0t1⋯tNMTAdt=UxM∫abtptqdtMTA=UxM∫abtp+qdtMTA=UxMZMTA,p,q=0,1,⋯,N,where
(49)Z=∫abtp+qdt,p,q=0,1,⋯,N,or
(50)Z=zpq=bp+q+1−ap+q+1p+q+1,p,q=0,1,⋯,N.So, the present integral term can be written as(51)∫abKx,tytdt=UxMZMTA=UxMZMTA.


### 3.2. Matrix Representation for Volterra Integral Terms

Now, we try to find the matrix form corresponding to the integral term. By the same way,Vcx,t can be expanded as (45)
(52)Vx,t≅∑r=0NfrxTrt.Then, the matrix representation of the kernel functionVcx,t is given by
(53)Vx,t≅FxTTt,where
(54)Fx=f0xf1x⋯fNx.Substituting the relations (29) and (53) in the present integral part, we obtained
(55)∫axVx,tytdt=∫axFxTTtTtAdt=∫axFxMXTtXtMTAdt=∫axFxMt0t1⋯tNTt0t1⋯tNMTAdt=FxM∫axtptqdtMTA=FxM∫axtp+qdtMTA=FxMExMTA,p,q=0,1,⋯,N,where
(56)Ex=∫axtp+qdt,p,q=0,1,⋯,N,or
(57)Ex=epqx=xp+q+1−ap+q+1p+q+1,p,q=0,1,⋯,N.So, the present integral term can be written as:(58)∫axVx,tytdt=FxMExMTA.


### 3.3. Conditions’ Matrix Relation

Finally, we can obtain the matrix form for condition (2) by using (36) on the form:
(59)XbiBiMTA=μi,i=0,1,2⋯,m−1,Via, Equations (41), (42), (43), (44), (51), and (58), then, Equation (1) converted to
(60)∑i=0n1QixXνixxBνixxMT+∑j=0n2PjxXαjxxBαjxxB−τjMT+∑s=0n3RsxXxBsB−εsMT−UxMZMT−FxMExMTA=gx.We can write (60) in the form:
(61)OA=G,orO;G,where
(62)O=∑i=0n1QixXνixxBνixxMT+∑j=0n2PjxXαjxxBαjxxB−τjMT+∑s=0n3RsxXxBsB−εsMT−UxMZMT−FxMExMT,The collocation points are defined in this form:(63)xl=h2l+12N+2,l=0,1,2,⋯,N.By substituting (63) in (60), then (60) can be turned into the following system:
(64)∑i=0n1QixlX¯νixlxlB¯νixlxlMT+∑j=0n2PjxlX¯αjxlxlB¯αjxlxlB−τjMT+∑s=0n3RsxlXxlBsB−εsMT−UxlMZMT−F¯ME¯MTA=G,whereX¯νixlxl,B¯νixlxl,X¯αjxlxl, and B¯αjxlxl are block matrices, also F¯ and E¯ are block matrices defined in this form:
(65)F¯=fx000⋯00fx10⋯000fx2⋯0⋮⋮⋮⋮000⋯fxN,E¯=ex000⋯00ex10⋯000ex2⋯0⋮⋮⋮⋮000⋯exN.G=gx1gx2⋮gxN.Here, the system in (64) in addition to the supplementary conditions can be solved numerically (the inverse matrix method is preferred when O be invertible) to determine the unknown vector A. Hence, the approximate analytical solution defined in series (18) can be calculated.


### 3.4. Convergence Analysis

Theorem 12.
Assume that a functionyx∈0,h be nth times continuously differentiable and ynx be the best square approximation of yx defined in (15). Then, we have
(66)yx−yNx≤MHN+1√lN+1!,where
(67)M=maxx∈0,hyn+1x,H=maxh−x0,x0andl=∫0h1xh−xdt.Proof.
Using Taylor expansion foryx as follows:
(68)yx=yx0+x−x0y′x0+⋯+x−x0NN!ynx0+x−x0N+1N+1!yN+1β,where x0∈0,h and β∈x0,x. Assume
(69)y¯Nx=yx0+x−x0y′x0+⋯+x−x0NN!yNx0+x−x0N+1N+1!yN+1β,then
(70)yx−yNx=x−x0N+1N+1!yN+1β.
According to,yNx that given in (18), we obtain
(71)yx−yNx2≤yx−y¯Nx2=∫0hωxyx−y¯Nx2dx=∫0hωxx−x0N+1N+1!yN+1β2dx≤M2N+1!2∫0hωxx−x0N+12dx.
Now, letH=maxh−x0,x0, thus (71) rewritten as
(72)yx−yNx2≤M2HN+12N+1!2∫0hωxdx.
Since,ω=1/xh−x,
Then(73)yx−yNx2≤M2HN+12N+1!2∫0h1xh−xdx.
Hence, the proof is completed.


## 3.1. Matrix Representation for Fredholm Integral Terms

Now, we try to find the matrix form corresponding to the integral term. Assume thatKx,t can be expanded to univariate Chebyshev series concerning t, as follows
(45)Kx,t≅∑r=0NurxTrt.Then the matrix representation of the kernel functionKx,t is given by
(46)Kx,t≅UxTTt,where
(47)Ux=u0xu1x⋯uNx.Substituting the relations (29) and (46) in the present integral part, we obtained
(48)∫abKx,tytdt=∫abUxTTtTtAdt=∫abUxMXTtXtMTAdt=∫abUxMt0t1⋯tNTt0t1⋯tNMTAdt=UxM∫abtptqdtMTA=UxM∫abtp+qdtMTA=UxMZMTA,p,q=0,1,⋯,N,where
(49)Z=∫abtp+qdt,p,q=0,1,⋯,N,or
(50)Z=zpq=bp+q+1−ap+q+1p+q+1,p,q=0,1,⋯,N.So, the present integral term can be written as(51)∫abKx,tytdt=UxMZMTA=UxMZMTA.


## 3.2. Matrix Representation for Volterra Integral Terms

Now, we try to find the matrix form corresponding to the integral term. By the same way,Vcx,t can be expanded as (45)
(52)Vx,t≅∑r=0NfrxTrt.Then, the matrix representation of the kernel functionVcx,t is given by
(53)Vx,t≅FxTTt,where
(54)Fx=f0xf1x⋯fNx.Substituting the relations (29) and (53) in the present integral part, we obtained
(55)∫axVx,tytdt=∫axFxTTtTtAdt=∫axFxMXTtXtMTAdt=∫axFxMt0t1⋯tNTt0t1⋯tNMTAdt=FxM∫axtptqdtMTA=FxM∫axtp+qdtMTA=FxMExMTA,p,q=0,1,⋯,N,where
(56)Ex=∫axtp+qdt,p,q=0,1,⋯,N,or
(57)Ex=epqx=xp+q+1−ap+q+1p+q+1,p,q=0,1,⋯,N.So, the present integral term can be written as:(58)∫axVx,tytdt=FxMExMTA.


## 3.3. Conditions’ Matrix Relation

Finally, we can obtain the matrix form for condition (2) by using (36) on the form:
(59)XbiBiMTA=μi,i=0,1,2⋯,m−1,Via, Equations (41), (42), (43), (44), (51), and (58), then, Equation (1) converted to
(60)∑i=0n1QixXνixxBνixxMT+∑j=0n2PjxXαjxxBαjxxB−τjMT+∑s=0n3RsxXxBsB−εsMT−UxMZMT−FxMExMTA=gx.We can write (60) in the form:
(61)OA=G,orO;G,where
(62)O=∑i=0n1QixXνixxBνixxMT+∑j=0n2PjxXαjxxBαjxxB−τjMT+∑s=0n3RsxXxBsB−εsMT−UxMZMT−FxMExMT,The collocation points are defined in this form:(63)xl=h2l+12N+2,l=0,1,2,⋯,N.By substituting (63) in (60), then (60) can be turned into the following system:
(64)∑i=0n1QixlX¯νixlxlB¯νixlxlMT+∑j=0n2PjxlX¯αjxlxlB¯αjxlxlB−τjMT+∑s=0n3RsxlXxlBsB−εsMT−UxlMZMT−F¯ME¯MTA=G,whereX¯νixlxl,B¯νixlxl,X¯αjxlxl, and B¯αjxlxl are block matrices, also F¯ and E¯ are block matrices defined in this form:
(65)F¯=fx000⋯00fx10⋯000fx2⋯0⋮⋮⋮⋮000⋯fxN,E¯=ex000⋯00ex10⋯000ex2⋯0⋮⋮⋮⋮000⋯exN.G=gx1gx2⋮gxN.Here, the system in (64) in addition to the supplementary conditions can be solved numerically (the inverse matrix method is preferred when O be invertible) to determine the unknown vector A. Hence, the approximate analytical solution defined in series (18) can be calculated.


## 3.4. Convergence Analysis

Theorem 12.
Assume that a functionyx∈0,h be nth times continuously differentiable and ynx be the best square approximation of yx defined in (15). Then, we have
(66)yx−yNx≤MHN+1√lN+1!,where
(67)M=maxx∈0,hyn+1x,H=maxh−x0,x0andl=∫0h1xh−xdt.Proof.
Using Taylor expansion foryx as follows:
(68)yx=yx0+x−x0y′x0+⋯+x−x0NN!ynx0+x−x0N+1N+1!yN+1β,where x0∈0,h and β∈x0,x. Assume
(69)y¯Nx=yx0+x−x0y′x0+⋯+x−x0NN!yNx0+x−x0N+1N+1!yN+1β,then
(70)yx−yNx=x−x0N+1N+1!yN+1β.
According to,yNx that given in (18), we obtain
(71)yx−yNx2≤yx−y¯Nx2=∫0hωxyx−y¯Nx2dx=∫0hωxx−x0N+1N+1!yN+1β2dx≤M2N+1!2∫0hωxx−x0N+12dx.
Now, letH=maxh−x0,x0, thus (71) rewritten as
(72)yx−yNx2≤M2HN+12N+1!2∫0hωxdx.
Since,ω=1/xh−x,
Then(73)yx−yNx2≤M2HN+12N+1!2∫0h1xh−xdx.
Hence, the proof is completed.


## 4. Numerical Examples and Results Analysis

In this section, we introduce some numerical examples for VO fractional delay-differential equation to illustrate the above results. All results are coded and obtained by using Mathematica package programming. Moreover, the absolute maximum errorEmax will be used in our computational results for the comparison between the exact and approximate solutions, where Emax is defined as the below formula:
(74)Emax=maxl=0,1,2,⋯,Nyxl−yNxl,xl=h2l+12N+2,l=0,1,2,⋯,N,where x, yx, and yNx are the space vectors, exact, and numerical solutions, respectively.Example 1.
Consider the following linear variable-order fractional delay integro-differential equation:(75)Dνxyx+Dαxyx−0.001+y′′x−0.002−yx=gx+∫0hx+tytdt+∫0x3t−2xytdt,x∈0,h.
The initial conditions arey0=0, y′0=−1, and the exact solution is yx=x2−x where h=1,νx=1−Sinx/9, αx=2−Cosx/7, Q0x=1, R0x=−1,R2x=1, P0x=1, and gx=25/12+7x/6−x2−x4/12+x1−αx1.002/−1.+αx+2.x/2−3αx+αx2/Gamma1−αx+x1−υx−2+2x+υx/Gamma1−υx−2+υx−1+υx. The fundamental matrix equation of the problem (75) at N=3 is defined by:
(76)Q0X¯νxlB¯νxlMT+P0X¯αxlB¯αxlB−τMT+R2XB2BεMT−R0XMT−UMZMT−F¯ME¯MTA=G,where
(77)Q0=R0=P0=idintity,B=0100002000030000,M=1000−12001−880−118−4832,X=10.1250.0156250.0019531310.3750.1406250.052734410.6250.3906250.24414110.8750.7656250.669922,G=1.92111.851111.944912.15722,B¯νx=000001.0387400001.0918300001.12841000001.0285500001.0660800001.0917000001.0189800001.0430200001.0593000001.0107600001.0239900001.03288,B¯αx=000001.0554600001.1373300001.19466000001.0584600001.1460100001.2075000001.0641400001.162900001.232701.0718700001.187100001.26924,X¯νx=00000.1530180000.01912730000.00239091000000000.4018140000.150680000.0565052000000000.8780230000.3991810000.29992000000000.9318990000.768270000.672236T,X¯αx=00000.2881080000.03601350000.00450169000000000.5627810000.2110430000.079141000000000.7677950000.4798720000.29992000000000.9318990000.8154110000.713485T,Z=1121314121314151314151614151617,B0.001=1−0.0011.1×10−6−1.0×10−901−0.0023.0×10−6001−0.0030001.In Table1, we list the absolute errors between the exact solutions and the approximate solutions using the proposed method with νx=1−Sinx/9, αx=2−Cosx/7 and three choices of N. In Figure 1 is the behavior of the absolute errors for the proposed method of different N.Example 2.
Consider the following linear variable-order fractional delay integro-differential equation [44]:
(78)Dαxyx−τ0=gx+∫0hx−tytdt+∫0xt+xytdt,x∈0,h.with the initial condition y0=1, y′0=1. Problem (78) found in [44] with the exact solution given as yx=ex at τ0=0. Now, by applying the suggested technique introduced in the previous section with finite terms N+1 to (78), then, we have
(79)P0XαxBαxB−τ0MT−UxMZMT−FxMExMTA=gx,
Hence, using the collocation method with the collection pointsxl that is given in the previous section. Then, the previous equation can be rewritten as follows:
(80)P0X¯αxlB¯αxlB−τ0MT−UMZMT−F¯ME¯MTA=G,Table 1
The comparisons between the absolute errors of differentN for Example 1.
xAbsolute errorsOur methodN=3Our methodN=4Our methodN=50.21.66×10−162.77×10−161.94×10−160.41.66×10−163.33×10−162.77×10−160.61.94×10−162.77×10−162.22×10−160.81.11×10−163.88×10−162.77×10−1611.38×10−163.40×10−161.75×10−16Figure 1
The behavior of the absolute errors for the proposed method withN=3,4,5 for Example 1.The numerical solutions obtained by different values ofαx together with the exact solution are displayed in Figure 2. From this figure, we see that the numerical solution converges to the exact one by increasing the number of basis functions. In Table 2, the comparisons between the absolute errors of the method are given in [44] and the presented method for different values of α and N=7 and choices of τ0=0 and h=1.Figure 2
Exact solution and numerical solutions (a), absolute error at two different cases ofαx (b) for Example 2 of N=7, τ=0, and h=1.
(a)(b)Table 2
The comparisons between the absolute errors of the method given in [44] the presented method with different values of αx, τ0=0 and N=7 for Example 2.
xAbsolute errorsOur methodMethod [44]Our methodMethod [44]αx=xαx=xαx=1−cos2xαx=1−cos2x0.13.97×10−95.04×10−65.08×10−97.34×10−60.35.43×10−91.37×10−66.27×10−92.86×10−60.55.53×10−91.63×10−65.64×10−93.55×10−60.76.32×10−99.62×10−75.94×10−93.17×10−70.92.01×10−88.21×10−72.09×10−86.10×10−7In Table3 shows the absolute errors between the exact and approximate solutions for ν=1 with N=6 and various choices of τ0. In Figure 3 is the exact solution and numerical solutions (a), absolute error (b) at N=6, τ0=0.001,τ0=0.0001, and h=1.Example 3.
Consider the following linear variable-order fractional delay integro-differential equation [45]:
(81)Dαxyx−τ=gx+13∫0hx+tytdt+110∫0xtxytdt,x∈0,h.where αx=x/2, gx=−17/216−5x/56−1/720x99+8x+720x6−αx7+7x−αx/Gamma8−αx and the initial condition is equal to y0=0. The exact solution as given in [45] is given by yx=x7+x6. Now, by applying the suggested technique that introduced in the previous section with finite terms N+1 to (53) then we have
(82)P0XαxBαxB−τMT−UxMZMT−FxMExMTA=gx,hence, using the collocation method with the collection points xl that given in the previous section. Then, the previous equation can be rewritten as follows:
(83)P0X¯αxlB¯αxlB−τMT−UMZMT−F¯ME¯MTA=G.Table 3
The comparisons between the absolute errors of the presented method for Example2 with different values of τ and N=7 in case αx=1−cos2x.
xAbsolute errorsOur methodτ=0.02Our methodτ=0.002Our methodτ=0.00020.15.73×10−95.14×10−95.09×10−90.36.86×10−96.33×10−96.28×10−90.56.47×10−96.54×10−96.70×10−90.76.76×10−95.98×10−95.94×10−90.92.22×10−82.13×10−82.09×10−8Figure 3
Exact solution and numerical solutions (a), absolute error (b) for Example2 with three different values of τ and N=7 in case αx=1−cos2x and h=1.
(a)(b)In Table4, the comparison is the maximum absolute error of the method given in [45] and the presented method for different values of N and h=1 in αx=x/2. The absolute error when h=1 for different values of N is shown in Figure 4.Example 4.
Consider the following linear variable-order fractional delay integro-differential equation [45]:
(84)Dαxyx−τ0=gx+∫01sinxtytdt+∫0xt−xytdt,x∈0,h,where P0=1, gx=−16x27/4/621−25x41/5/1476+x19/4−αxGamma23/4/Gamma23/4−αx+x29/20Gamma36/5/Gamma36/5−αx−299Sinx/1107, αx=x and the initial condition is equal to y0=0. The exact solution of this equation is given by yx=x19/4+x31/5 at τ0=0 (see [45]). Now, by applying the suggested technique that introduced in the previous section with finite terms N+1 to (8), then, we have
(85)P0XαxBαxB−τ0MT−UxMZMT−FxMExMTA=gx,
Hence, using the collocation method with the collection pointsxl that given in the previous section. Then, the previous equation can be rewritten as follows:
(86)P0X¯αxlB¯αxlB−τ0MT−UMZMT−F¯ME¯MTA=G.Table 4
Comparison of the maximum absolute error of the method given in [45] and the presented method for Example 3 with different values of N and h=1 in αx=x/2 for Example 3.
xAbsolute errorsOur method atN=7[45] at N=7.Our method atN=5[45] at N=50.19.05×10−162.10×10−61.56×10−42.88×10−30.31.80×10−155.45×10−65.58×10−44.72×10−40.53.49×10−153.50×10−66.27×10−41.12×10−40.73.49×10−155.25×10−46.94×10−42.01×10−30.99.54×10−154.03×10−24.18×10−34.39×10−2Figure 4
The absolute error whenh=1 for three different values of N for Example 3.By solving the algebraic system (86), we can obtain the vector A=a0,a1,⋯,anT . Subsequently, numerical solution according to (18) is obtained. Table 5 gives the comparison of the absolute error at some points obtained by the present method and method of [45] with different values of N, and it is seen that our method gives more accurate results than the method of [45] at τ0=0. In Figure 5, the exact solution and numerical solutions given (a), where the absolute error shown (b), are all with different values of N and h=1.Example 5.
Consider the following linear variable-order fractional integro-differential equation [39, 46]:
(87)Dνxyx+sinxy′x+x−1yx=gx−∫0xexytdt,x∈0,1,where Q0=1, R0x=x−1, R1x=sinx, and gx=−1+x1+−1+xx+1/6exx6+x−3+2x+−1+2xSinx+x1−υx−2+2x+υx/Gamma3−υx and the initial conditions is equal to y0=0,y′0=−1. The exact solution of (87) according to [39, 46] is given by yx=x2−x+1 at νx=x/3. Now, by applying the suggested technique that introduced in the previous section with finite terms N+1 to (8), then, we have
(88)Q0XνxBνxMT+R0XxMT+R1XxBMT+FxMExMTA=gx,
Hence, using the present collocation method with the collection pointsxl that given in the previous section. Then, the previous equation can be rewritten as follows:
(89)Q0X¯νxlB¯νxlMT+R0XxlMT+R1XxlBMT+F¯ME¯MTA=G.Table 5
Comparison of the absolute error of the method given in [45] and the presented method for Example 4 with different values of N, h=1, and τ0=0.
xAbsolute errorsOur method atN=7[45] at N=7Our method atN=5[45] at N=50.12.95×10−68.24×10−48.95×10−61.31×10−30.32.31×10−69.11×10−51.89×10−51.28×10−40.53.76×10−61.25×10−32.39×10−53.72×10−30.75.88×10−65.25×10−49.57×10−32.64×10−20.92.64×10−59.57×10−38.65×10−21.81×10−1Figure 5
Exact solution and numerical solutions (a), absolute error (b) for Example4 with different values of N and h=1 in case of h=1.
(a)(b)By solving the algebraic system (89), we can obtain the vector A=a0,a1,⋯,anT . Subsequently, numerical solution (18) is obtained. In Table 6, a comparison between the AEs is obtained in [39, 46], and the results are obtained in this work for N=4 of for example (12). In Figure 6, we see the matching of the AEs values in this figure and that in Table 6. While Figure 7 displays the curves of exact and approximate solutions for example (12) with values of parameters listed in their caption.Table 6
The comparisons between the absolute errors of the presented method for Example5 with value of N=4.
xAbsolute errorsOur methodMethod [46]Method [39]0.11.11×10−163.52×10−134.99×10−160.31.11×10−162.33×10−131.22×10−150.501.35×10−131.05×10−150.704.86×10−146.66×10−160.902.99×10−145.55×10−16Figure 6
Exact solution and numerical solutions of Example4-case 1 for different values of N,h.Figure 7
Comparison between the approximate solutions and the exact solutions of Example5 for N=4.


## 5. Conclusion

In this work, the general form of mixed Fredholm–Volterra integro-differential equations with delay arguments of variable-order based on the operational matrix method with the Shifted Chebyshev polynomials (CPs) of the first kind is presented. The spectral collocation technique with the aid of CPs is used as an operational matrix method for solving the proposed model, which is reduced by the operational matrices to the matrix form. Caputo’s VO fractional derivative operator is used to generate the matrices of the derivative. The accuracy of the proposed technique is obtained by many numerical examples. Finally, we used codes written with the Mathematica package to calculate our numerical results and graphs.



---

*Source: 6203440-2022-04-09.xml*