# A New Efficient Method for Solving System of Weakly Singular Fractional Integro-Differential Equations by Shifted Sixth-Kind Chebyshev Polynomials

Authors: S. Yaghoubi; H. Aminikhah; K. Sadri
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this paper, a new approach for solving the system of fractional integro-differential equation with weakly singular kernels is introduced. The method is based on a class of symmetric orthogonal polynomials called shifted sixth-kind Chebyshev polynomials. First, the operational matrices are constructed, and after that, the method is described. This method reduces a system of weakly singular fractional integro-differential equations (WSFIDEs) by the collocation method into a system of algebraic equations. Thereupon, an upper error bound for the proposed method is determined. Finally, some numerical examples are prepared to test the accuracy and efficiency of the presented method.


---

## Body

## 1. Introduction

The study of fractional calculus has applications and popularity in various and wide fields of biology, physics, and fluid mechanics. Fractional calculus is actually integration and differentiation of arbitrary orders [1–4]. In various problems of physics and engineering, the fractional differential equations have been proved to be valuable tools in modeling of many phenomena [5, 6]. As we know, many mathematical models of real phenomena (arising in engineering and physics) are described as linear or nonlinear systems. It is worth mentioning that with the development of fractional calculus, the behavior of many systems can be described using the fractional differential and fractional integro-differential system [7, 8]. In recent years, systems of the fractional differential and integral equations are the subjects of extensive study due to their frequent appearance in many engineering and scientific disciplines [9–11]. However, most of the fractional-order equations and integral equations do not have analytic solutions or are hard to find. So, it is essential to find numerical methods to get approximate or exact solutions of a system of integro-differential equations. So far, researchers have utilized diverse numerical methods for solving a system of fractional integro-differential equations. In [12], the homotopy perturbation method was proposed for solving linear and nonlinear systems of fractional integro-differential equations. Heydari et al. have used the Chebyshev wavelet method for solving a class of systems of nonlinear fractional singular Volterra integro-differential equation in [13]. In the next year, for the first time, the hybrid functions composed of the Block-pulse functions and Bernoulli polynomials were applied for problems with fractional-order differential equations in [14]. Also, a novel technique based on iterative refinement was presented to analytically approximate a system of linear fractional integro-differential equations [15]. In 2018, Hesameddini and Shahbazi developed the concept of [14] and used it to solve the FDIE system in [16]. Also, Xie and Yi presented the simple and fast method based on the Block-pulse function to solve a nonlinear system of fractional Volterra-Fredholm integro-differential equations in the same year [17]. In [18], the authors implemented the new Jacobi operational matrices to reduce the complexity of calculations to solve WSFDIEs. Next, the Haar wavelet method was employed to solve a coupled system of FIDEs, and the Muntz-Legendre wavelets were introduced to solve FIDVFEs in [19, 20]. Also, the other authors applied the Chebyshev Pseudo spectral method for solving fractional-order nonlinear system of Volterra integro-differential equations and a least square collocation Chebyshev technique for solving a system of linear fractional integro-differential equations [21, 22]. In this paper, we consider the following system of weakly singular integro-differential equations:
(1)Dνiuix=Fix,u1x,u2x,⋯,umx+∑j=1mθij∫0xKijx,zGijz,ujzx−zαijdz+fix,i=1,⋯,m,uik0+=u0ik,k=0,1,⋯,r−1,r−1<νi≤r,whereuix∈Cr0,1,i=1,⋯,m, are the unknown functions, Fi:0,1×Cr0,1m⟶ℝ, Gij:0,1×Cr0,1⟶ℝ, and Kijx,z:0,12⟶ℝ are continuous operators and functions that satisfy Lipschitz conditions, and Dνi is the Caputo fractional derivative operator where νi∈0,1. The parameters θij,αij∈ℝ such that θij≤1, 0<αij<1, and i,j=1,2,⋯,m. Moreover, fix,i=1,⋯,m, are known and sufficiently smooth functions.As usual, a way for solving functional equations is to express the solution as a linear combination of the so-called basis functions. In most researches, various polynomials such as the Legendre, Chebyshev, Taylor, Hermit, and Bernstein are used as basis functions. Among all of them, the Chebyshev polynomials are the most important in the analysis and numerical analysis. Chebyshev polynomials are orthogonal on the interval−1,1 and have good properties that are used widely in the approximation of the functions. For this reason, many studies are done based on the different kinds of Chebyshev polynomials. In [23], Masjed-Jamei introduced two half-trigonometric orthogonal Chebyshev polynomials, and he named them as the Chebyshev polynomials of the fifth and sixth kinds. The basic formulas and properties of this class of polynomials are displayed in [24, 25]. Up to now, many researchers have used various kinds of Chebyshev polynomials for the fractional-order differential and integro-differential equations (see [26, 27]). However, there are only a few works that have used the sixth-kind Chebyshev polynomials. The main aim of this work is to introduce these polynomials as a new basis function for solving WSFDIEs. In the current paper, we apply the orthogonal shifted sixth-kind Chebyshev polynomials together with the collocation method for solving a system of weakly singular integro-differential equations with fractional derivatives that, to the best of our knowledge, is proposed here for the first time. For solving these equations, we derive the fractional operational matrices of fractional and integer orders and the product operational matrix, as well. Also, we introduce an operational matrix to approximate the integral term that has the weakly singular kernel in Equation (1). As far as we can tell, this operational matrix is presented for the first time. By substituting appropriate approximations in Equation (1), the original equations convert into algebraic equations that each of the equations of algebraic systems is collocated at N+1 roots of the N+1th shifted sixth-kind Chebyshev polynomials (SSKCPs). By solving these algebraic systems, the approximate solutions of the original system are obtained. Although the calculation of the operational matrices may be complicated, we show that the obtained results are equal to other methods or are even more accurate. Implementing these matrices leads to a decrease in the number of required computations, and therefore, the computation time will be reduced.The rest of the paper is organized as follows. In section2, some essential preliminaries are mentioned briefly. Section 3 is devoted to constructing the operational matrices of SSKCPs. The proposed numerical procedure is described in Section 4. The error analysis of the proposed method is discussed in Section 5. Some numerical applications are indicated in Section 6, and conclusions are presented in Section 7.


## 2. Preliminaries and Notations

In this section, we recall some definitions and properties of fractional integral and derivative operators which will be used later. After that, some necessary definitions and fundamental properties of the shifted sixth-kind Chebyshev polynomials are reviewed briefly.


### 2.1. Some Essentials of the Fractional Calculus

Definition 1.
The Riemann-Liouville fractional integral operatorJα of order α is given by [2]
(2)Jαfx=1Γα∫0xx−zα−1fzdz,α>0,x>0.Definition 2.
Letα∈ℝ, n−1<α≤n, n∈ℕ, and fx be a real-valued continuous function defined on 0,∞. Then, the Caputo fractional derivative of order α>0 is defined by [2]
(3)0Dxαfx=1Γn−α∫0xfnzx−zα+1−ndz,fnz,α=n,where Γx is the Gamma function as
(4)Γx=∫0∞e−zzx−1dz,Rez>0,Γx+1=xΓx,Bu,v=∫01zu−11−zv−1dz=ΓuΓvΓu+v,Reu>0,Rev>0.The last integral is often called the Beta integral. Forα1,α2>0, the Riemann-Liouville integral and Caputo fractional derivative operators satisfy the following properties:
(1)
Jα1Jα2fx=Jα2Jα1fx=Jα1+α2fx(2)
Jαλ1fx+λ2gx=λ1Jαfx+λ2Jαgx(3)
JαDαfx=fx−∑i=0n−1fi0xi/i!,n−1<α≤n,x>0(4)
Dαxγ=0,α>γ,Γγ+1/Γγ−α+1xγ−α,otherwise(5)
Jαxv=Γv+1/Γv+α+1xv+α,v>−1


### 2.2. Shifted Sixth-Kind Chebyshev Polynomials and Their Properties (SSKCPs)

Definition 3.
The sixth-kind Chebyshev polynomials are orthogonal functions on the interval−1,1 and can be determined with the following recursive formula [23, 24]:
(5)Sjx=xSj−1x−jj+1+−1j2j+1+14jj+1Sj−2x,j≥2,S0x=1,S1x=x.Definition 4.
The shifted sixth-kind Chebyshev polynomials on0,1 is defined by [23, 24]
(6)Sj∗x=Sj2x−1,j=0,1,2,⋯.
These polynomials have the following explicit analytic form:(7)Sj∗x=∑k=0jρkjxk,where
(8)ρkj=22k−j2k+1!∑j/2i=k+1/2−1j/2+i+k2i+k+1!2i−k!,jeven,22k−j+12k+1!j+1∑j−1/2i=k/2−1j+1/2+i+ki+12i+k+2!2i−k+1!,jodd.Moreover, the shifted polynomialsSj∗x are orthogonal on 0,1 with respect to the weight function Vx=2x−12x−x2 in the sense that
(9)∫01Si∗xSj∗xVxdx=λiδij,(10)λi=π22i+5,ieven,πi+322i+5i+1,iodd.Now, lethx∈L20,1; then, hx can be approximated in terms of Sj∗x as
(11)hx≈∑j=0NϱjSj∗x=FTSx=STxF,where
(12)Sx=S0∗x,S1∗x,⋯,SN∗xT,F=ϱ0,ϱ1,⋯,ϱNT,where the coefficients ϱj are given by
(13)ϱj=1λj∫01hxSj∗xVxdx,and λj is defined in Equation (10). Similarly, any continuous two-variable function, Fx,z, defined on 0,1×0,1 can be approximated by means of the double-shifted sixth-kind Chebyshev polynomials as
(14)Fx,z≈∑j=0N∑i=0NFijSi∗xSj∗z=STxFSz,where F is a N+1×N+1 matrix, and its entries are given by
(15)Fij=1λiλj∫01∫01Fx,zSi∗xSj∗zVxVzdxdz.


## 2.1. Some Essentials of the Fractional Calculus

Definition 1.
The Riemann-Liouville fractional integral operatorJα of order α is given by [2]
(2)Jαfx=1Γα∫0xx−zα−1fzdz,α>0,x>0.Definition 2.
Letα∈ℝ, n−1<α≤n, n∈ℕ, and fx be a real-valued continuous function defined on 0,∞. Then, the Caputo fractional derivative of order α>0 is defined by [2]
(3)0Dxαfx=1Γn−α∫0xfnzx−zα+1−ndz,fnz,α=n,where Γx is the Gamma function as
(4)Γx=∫0∞e−zzx−1dz,Rez>0,Γx+1=xΓx,Bu,v=∫01zu−11−zv−1dz=ΓuΓvΓu+v,Reu>0,Rev>0.The last integral is often called the Beta integral. Forα1,α2>0, the Riemann-Liouville integral and Caputo fractional derivative operators satisfy the following properties:
(1)
Jα1Jα2fx=Jα2Jα1fx=Jα1+α2fx(2)
Jαλ1fx+λ2gx=λ1Jαfx+λ2Jαgx(3)
JαDαfx=fx−∑i=0n−1fi0xi/i!,n−1<α≤n,x>0(4)
Dαxγ=0,α>γ,Γγ+1/Γγ−α+1xγ−α,otherwise(5)
Jαxv=Γv+1/Γv+α+1xv+α,v>−1


## 2.2. Shifted Sixth-Kind Chebyshev Polynomials and Their Properties (SSKCPs)

Definition 3.
The sixth-kind Chebyshev polynomials are orthogonal functions on the interval−1,1 and can be determined with the following recursive formula [23, 24]:
(5)Sjx=xSj−1x−jj+1+−1j2j+1+14jj+1Sj−2x,j≥2,S0x=1,S1x=x.Definition 4.
The shifted sixth-kind Chebyshev polynomials on0,1 is defined by [23, 24]
(6)Sj∗x=Sj2x−1,j=0,1,2,⋯.
These polynomials have the following explicit analytic form:(7)Sj∗x=∑k=0jρkjxk,where
(8)ρkj=22k−j2k+1!∑j/2i=k+1/2−1j/2+i+k2i+k+1!2i−k!,jeven,22k−j+12k+1!j+1∑j−1/2i=k/2−1j+1/2+i+ki+12i+k+2!2i−k+1!,jodd.Moreover, the shifted polynomialsSj∗x are orthogonal on 0,1 with respect to the weight function Vx=2x−12x−x2 in the sense that
(9)∫01Si∗xSj∗xVxdx=λiδij,(10)λi=π22i+5,ieven,πi+322i+5i+1,iodd.Now, lethx∈L20,1; then, hx can be approximated in terms of Sj∗x as
(11)hx≈∑j=0NϱjSj∗x=FTSx=STxF,where
(12)Sx=S0∗x,S1∗x,⋯,SN∗xT,F=ϱ0,ϱ1,⋯,ϱNT,where the coefficients ϱj are given by
(13)ϱj=1λj∫01hxSj∗xVxdx,and λj is defined in Equation (10). Similarly, any continuous two-variable function, Fx,z, defined on 0,1×0,1 can be approximated by means of the double-shifted sixth-kind Chebyshev polynomials as
(14)Fx,z≈∑j=0N∑i=0NFijSi∗xSj∗z=STxFSz,where F is a N+1×N+1 matrix, and its entries are given by
(15)Fij=1λiλj∫01∫01Fx,zSi∗xSj∗zVxVzdxdz.


## 3. Operational Matrices of SSKCPs

In this section, the formulas of operational matrices with the fractional order will be derived for the sixth-kind Chebyshev polynomials. The following are the required lemmas.Lemma 5.
Ifr≥l, l∈ℕ, then we have
(16)∫01xrSl∗xVxdx=∑m=0lρmlπΓr+m+3/22Γr+m+5m2+m+r2+2rm+3+r.Proof.
From the properties of the orthogonal polynomials, ifr≤l, we have
(17)∫01xrSl∗xVxdx=0.
Hence, we supposer≥l. The lemma can be easily proved by the integration of the analytic form of SSKCPs in Equation (7).Theorem 6.
LetSx be the SSKCP vector given by Equation (12) and μ∈ℝ; then,
(18)JμSx≈PμSx,where Pμ is the N+1×N+1 operational matrix of the fractional integration of the order μ in the Riemann-Liouville sense which is defined by
(19)Pμ=a~00a~01⋯a~0Na~10a~11⋯a~1N⋮⋮⋱⋮a~N0a~N1⋯a~NN,(20)a~ij=∑l=0iϖijl,i=0⋯N,j=0⋯N.ϖijl are given by
(21)ϖijl=ρliΓl+1λjΓl+μ+1∑k=0jρkjπΓk+μ+l+3/22Γk+μ+l+5×k+μ2+2l+1k+μ+ll+1+3,i=0⋯N,j=0⋯N.Proof.
By applying the Riemann-Liouville integral operator to the SSKCPs’ analytic form, we have(22)JμSi∗x=∑l=0iρliΓl+1Γl+μ+1xμ+l.Now, we can express xμ+l in terms of the shifted sixth-kind Chebyshev polynomials as follows:
(23)xμ+l≈∑j=0NC~ljSj∗x,where the coefficients C~lj are given by
(24)C~lj=1λj∫01xμ+lSj∗xVxdx.
According to Lemma5, we can rewrite Equation (22) as
(25)JμSj∗x≈∑j=0N∑l=0iρliΓl+1λjΓl+μ+1×∑k=0jρkjπΓk+μ+l+3/22Γk+μ+l+5k+μ2+2l+1k+μ+ll+1+3Sj∗x=∑j=0Na~ijSj∗x,where a~ij is given in Equation (20). Rewriting the last relation in the vector form gives
(26)JμSj∗x=a~i0,a~i1,⋯,a~iNSj∗x,i=0,1,⋯,N.
This leads to the desired result.In the following, some useful and applicable lemmas are presented to get the Chebyshev operational matrix of product.Lemma 7.
IfSj∗x and Si∗x are jth and ith shifted sixth-kind Chebyshev polynomials, respectively, then we can write the product of Sj∗x and Si∗x as
(27)Qi+jx=Si∗xSj∗x=∑k=0i+jχki,jxk.Proof.
See [18].Lemma 8.
IfSi∗x, Sj∗x, and Sk∗x are ith, jth, and kth shifted sixth-kind Chebyshev polynomials, then
(28)dijk=∫01Si∗xSj∗xSk∗xVxdx=∑r=0j+k∑l=0iπρliχki,jΓr+l+3/22Γr+l+5l+rl+r+1+3,where χki,j is obtained by Lemma 7.Proof.
According to Lemma7, we can write
(29)Qj+kx=Sj∗xSk∗x=∑r=0j+kχrj,kxr.Then,
(30)dijk=∫01Si∗x∑r=0j+kχrj,kxrVxdx=∑r=0j+kχrj,k∫01Si∗xxrVxdx.
The value of the last integral is obtained by Lemma5.Assuming thatE is a N+1×1 vector, we have
(31)SxSTxE≃E~Sx,where E~ is a N+1×N+1 matrix called the product operational matrix. The next theorem presents a general form for entries of the matrix E~.Theorem 9.
The entries of the matrixE~ in Equation (31) are as follows:
(32)E~jk=1λk∑i=0NEidijk,j,k=0,1,⋯,N,where dijk is obtained by Lemma 8 and Ei is the element of the vector E.Proof.
See [18].In the following, we get an approximation for the integral part with the singular kernel in Equation (1). Before that, we present a theorem.Theorem 10.
The following relation is determined for0<κ<1:
(33)∫0xzrx−zκdz=Γr+1Γ1−κΓr−κ+2xr−κ+1,r=0,1,2,⋯.Proof.
By performing Equation (33) and the substitution of z=ξx into Equation (33) and then using the definition of the Beta function, we obtain
(34)xr−κ+1∫011−ξ−κξrdξ=Γ1−κΓr+1Γr−κ+2xr−κ+1,r=0,1,2,⋯.Now, we present an approximation for the integral part with a weakly singular kernel. For this purpose, see the following theorem.Theorem 11.
Suppose thatux is a continuous function on the interval 0,1 and κ∈0,1 and ux≈STxF=FTSx where Sx and F are defined by Equation (12), then we have
(35)∫0xuzx−zκdz≈FTIκSx,where Iκ is a N+1×N+1 matrix as follows:
(36)Iκ=σ00σ01⋯σ0Nσ10σ11⋯σ1N⋮⋮⋱⋮σN0a~N1⋯σNN,and its entries are determined as follows:
(37)Iijκ=σij=∑l=0iρliΓl+1Γ1−κΓl−κ+2C¯jl−κ+1,i,j=0,1,⋯,N,where the quantities ρli and C¯jl−κ+1 are introduced by relation (8) and Lemma 5.Proof.
By the definition of vectorSx and Lemma 5, we can write
(38)STx=S0∗x,S1∗x,⋯,SN∗x=∑l=00ρl0xl,⋯,∑l=0NρlNxl.By applying Theorem 10, we have
(39)∫0xSTzx−zκdz=∑l=00∫0xρl0zlx−zκdz,⋯,∑l=0NρlN∫0xzlx−zκdz=∑l=00ρl0Γl+1Γ1−κΓl−κ+2xl−κ+1,⋯,∑l=0NρlNΓl+1Γ1−κΓl−κ+2xl−κ+1.
Now, we approximatexl−κ+1 in terms of SSKCPs as follows:
(40)xl−κ+1≈∑j=0NC¯jl−κ+1Sj∗x,C¯jl−κ+1=1λj∫01xl−κ+1Sj∗xVxdx,j=0,1,⋯,N,where C¯jl−κ+1 is obtained by applying Lemma 5. Thus,
(41)∑l=0iρliΓl+1Γ1−κΓl−κ+2xl−κ+1≈∑j=0N∑l=0iρliΓl+1Γ1−κC¯jl−κ+1Γl−κ+2Sj∗x=∑j=0NσijSj∗x,i=0,1,⋯,N.
Equation (39) is then represented as follows:
(42)∫0xSTzz−xκdz≈σ00σ01⋯σ0Nσ10σ11⋯σ1N⋮⋮⋱⋮σN0a~N1⋯σNNS0∗xS1∗x⋮SN∗x=IκSx.


## 4. Numerical Procedure

We consider the system of fractional integro-differential equation with weakly singular kernels described in Equation (1). To solve the system, we approximate the function Dνiu in a matrix form
(43)Dνiuix≈STxFi,Fi=ϱ0i,ϱ1i,⋯,ϱNiT,i=1,2,⋯,m.Then, according to the initial conditions of the problem, we can approximate the known functionHix=∑k=0r−1u0ikxk/Γk+1,k=νi as follows:
(44)Hix≈STxCi,i=1,2,⋯,m.Using Equations (43) and (44), we compute an approximate for uix:
(45)uix≈STxPνiFi+STxCi=STxUi,i=1,2,⋯,m,where Pνi is the integral operational matrix presented in Theorem 6. We have the following approximations for the rest of the system:
(46)Fix,u1x,u2x,⋯,umx≈XiTSx,Gijx,ujx≈STxVij,i,j=1,2,⋯,m,Kijx,z≈STxKijSz.Using Theorems9 and 11, we have
(47)∫0xKijx,zGijz,yjzx−zαijdz≈∫0xSTxKijSzSTzVijx−zαijdz≈STxKijV~ij∫0xSzx−zαijdz≈STxKijV~ijIαijSx.Now, by substituting Equations (43)-(47) into Equation (1), we obtain
(48)STxFi−XiTSx−∑j=0mθijSTxKijV~ijIαijSx−fix≈0,i=1,2,⋯,m.Each equation of algebraic system (48) is collocated at N+1 roots of the N+1th shifted sixth-kind Chebyshev polynomials. Thus, an algebraic system, including mN+1 equations, is acquired. By solving the resultant algebraic system, we can obtain an approximation for unknown vectors Fi,i=1,2,⋯,m, and by substituting the vector Fi into Equation (45), we obtain an approximation for uix,i=1,2,⋯,m.


## 5. Error Analysis

In this section, we prove some theorems. Then, we obtain an upper error bound for the approximation error. For this aim, we need the following norms:(49)fL2I=∫Ifx2Vxdx1/2,X1=∑i=0nxi,where f∈L2I is a square integrable function on the interval I=0,1 and X=x0,x1,⋯,xnT is a vector.Theorem 12.
Suppose thatYNx=∑i=0NEiSi∗x is an approximation in SSKCPs to the continuous function Yx on the interval 0,1. Then, the coefficients Ei, for i=0,1,⋯,N, are bounded as
(50)Ei≤MYλi∑m=0iρmiπΓm+3/22Γm+5m2+m+3,where MY denotes the maximum value of Yx on the interval 0,1.Proof.
Using Equations (7) and (9) for i=0,1,⋯,N, we have
(51)Ei=1λi∫01YxSi∗xVxdx=1λi∫01Yx∑m=0iρmixmVxdx=1λi∑m=0iρmi∫01YxxmVxdx.
SinceYx is a continuous function on the interval 0,1, so it is bounded and there is a constant MY such that
(52)∀x∈0,1,Yx≤MY.Using Equations (51) and (52), inequality (50) is deduced.Theorem 13.
Suppose thatYx is a continuous function and YNx is an approximation to Yx in terms of SSKCPs. Then, the error bound can be achieved as follows:
(53)Yx−YNxL2≤∑i=N+1∞Ωi1/2=ΩN,where
(54)Ωi=MY2λi∑m=0iρmiπΓm+3/22Γm+5m2+m+32.Proof.
AssumeYx is an arbitrary function. So, Yx and YNx have the following forms using SSKCP series:
(55)Yx=∑i=0∞EiSi∗x,YNx=∑i=0NEiSi∗x,so,
(56)Yx−YNx=∑i=N+1∞EiSi∗x.
Using Equations (9) and (56) and Theorem 12, we have
(57)Yx−YNxL22=∫01Yx−YNx2Vxdx=∫01∑i=N+1∞EiSi∗x2Vxdx=∫01∑j=N+1∞∑i=N+1∞EiEjSi∗xSj∗xVxdx=∑i=N+1∞Ei2λi≤∑i=N+1∞Ωi.Theorem 14.
Suppose that the continuous two-variable functionHx,y is approximated on the interval 0,1×0,1 in terms of SSKCPs as HNx,y=∑i=0N∑j=0NHijSi∗xSj∗y; then, the coefficients Hij, for i,j=0,1,⋯,N, can be bounded as follows:
(58)Hij≤MHπ4λiλj∑m=0iρmiΓm+3/2Γm+5m2+m+3∑r=0jρrjΓr+3/2Γr+5r2+r+3,where MH denotes the maximum value of Hx,y on the interval 0,1×0,1.Proof.
Using Equations (7) and (15), we have
(59)Hij=1λiλj∫01∫01Hx,ySi∗xSj∗yVxVydxdy=1λiλj∫01∑m=0iρmixmVx∫01Hx,y∑r=0jρrjyrVydydx=1λiλj∑m=0iρmi∑r=0jρrj∫01∫01xmHx,yyrVxVydxdy.
SinceHx,y is a bounded and continuous function on the interval 0,1×0,1, so there is a constant MH such that
(60)∀x,y∈0,1×0,1,Hx,y≤MH.
Using Equations (59) and (60), Theorem 14 is proved.Theorem 15.
Suppose thatHx,y is a continuous function with two variables and HNx,y is the approximation to Hx,y using SSKCPs. Then, the error bound can be obtained as
(61)Hx,y−HNx,yL2≤∑i=0N∑j=N+1∞ςij2λiλj1/2+∑i=N+1∞∑j=0∞ςij2λiλj1/2=ΛH,where
(62)ςij=MHπ4λiλj∑m=0iρmiΓm+3/2Γm+5m2+m+3∑r=0jρrjΓr+3/2Γr+5r2+r+3.Proof.
Suppose thatHx,y is an arbitrary function. SSKCP series of Hx,y and its approximation in terms of SSKCPs have the following form:
(63)Hx,y=∑i=0∞∑j=0∞HijSi∗xSj∗y,HNx,y=∑i=0N∑j=0NHijSi∗xSj∗y.Thus,
(64)Hx,y−HNx,y=∑i=0N∑j=N+1∞HijSi∗xSj∗y+∑i=N+1∞∑j=0∞HijSi∗xSj∗y.Using Equations (9) and (64) and Theorem 14, we conclude
(65)Hx,yx−HNx,yL2≤∑i=0N∑j=N+1∞HijSi∗xSj∗yL2+∑i=N+1∞∑j=0∞HijSi∗xSj∗yL2=∫01∫01∑i=0N∑j=N+1∞HijSi∗xSj∗y2VxVydydx1/2+∫01∫01∑i=N+1∞∑j=0∞HijSi∗xSj∗y2VxVydydx1/2=∑i=0N∑j=N+1∞Hij2λiλj1/2+∑i=N+1∞∑j=0∞Hij2λiλj1/2≤∑i=0N∑j=N+1∞ςij2λiλj1/2+∑i=N+1∞∑j=0∞ςij2λiλj1/2.In the following theorem, we obtain an upper error bound of the proposed method. First, suppose that for theuix,yix∈Cr0,1,i=1,2,⋯,m, there exist positive constants ξri,ηij>0 such that the following Lipschitz conditions hold.
(1)
Fix,u1x,u2x,⋯,umx−Fix,y1x,y2x,⋯,ymxL2≤ξ1iu1x−y1xL2+ξ2iu2x−y2xL2+⋯+ξmiumx−ymxL2,i=1,2,⋯,m(2)
Gijz,ujz−Gijz,yjzL2≤ηijujz−yjzL2,i,j=1,2,⋯,mTheorem 16.
Suppose thatU~x=u~1x,u~2x,⋯,u~mx is a set of approximate solutions obtained from the SSKCP collocation method, Ux=u1x,u2x,⋯,umx is the set of exact solutions of system (1), E=u1x−u~1xL2,⋯,umx−u~mxL2 is the error vector of approximate solutions, and also rix,i=1,2,⋯,m denote the residual functions associated to the approximate solutions that are named perturbation terms. Assume that Hypotheses (1) and (2) are satisfied; then, a bound for the method error can be achieved as
(66)E1≤ΓN∗+∑i=1m∑j=1mγji1−ΔN∗,0<ΔN∗<1.Proof.
First, we apply the Riemann-Liouville integral operator on Equation (1) and obtain the following equation:
(67)uix=gix+1Γνi∫0xx−zνi−1Fiz,u1z,u2z,⋯,umzdz+∑j=1mθijΓ1−αijΓνi−αij+1∫0xx−zνi−αijKijx,zGijz,ujzdz,where
(68)gix=∑k=0r−1u0ikΓk+1xk+1Γνi∫0xx−zνi−1fizdz.
We can write the approximate equation of Equation (67) as follows:
(69)u~ix=gix+1Γνi∫0xx−zνi−1Fiz,u~1z,u~2z,⋯,u~mzdz+∑j=1mθijΓ1−αijΓνi−αij+1∫0xx−zνi−αijK~ijx,zGijz,u~jzdz+rix,where rix is the perturbation term. We subtract Equation (69) from Equation (67) and obtain the following result:
(70)rix=uix−u~ix+1Γνi∫0xx−zνi−1Fiz,u1z,u2z,⋯,umz−Fiz,u~1z,u~2z,⋯,u~mzdz+∑j=1mθijΓ1−αijΓνi−αij+1∫0xx−zνi−αijKijx,zGijz,ujz−K~ijx,zGijz,u~jzdz.
First, we obtain a bound for the perturbation term, so by taking theL2-norm on Equation (70), we get the following inequality:
(71)rixL2≤uix−u~ixL2+1Γνi∫0xx−zνi−1Fiz,u1z,⋯,umz−Fiz,u~1z,⋯,u~mzdzL2+∑j=1mθijΓ1−αijΓνi−αij+1×∫0xx−zνi−αijKijx,zGijz,ujz−K~ijx,zGijz,u~jzdzL2.
Using Hypothesis 1 and Theorem13, we have
(72)1Γνi∫0xx−zνi−1Fiz,u1z,u2z,⋯,umz−Fiz,u~1z,u~2z,⋯,u~mzdzL2≤Ψνi∑j=1mξjiujz−u~jzL2≤Ψνi∑j=1mξjiΩNj,where
(73)Ψνi=1Γνi4Γ7/2Γ2νi−1/2Γ2νi+3−4Γ5/2Γ2νi−1/2Γ2νi+2+Γ3/2Γ2νi−1/2Γ2νi+11/2.
Noting thatKijx,z and Gijz,ujz are continuous and known functions, thus, there are constants MKij and NGij such that
(74)Kijx,zL2≤MKij,Gijz,ujzL2≤NGij.
From Hypothesis 2, Theorem13, Theorem 15, and Equation (74), the following inequality is obtained:
(75)∑j=1mθijΓ1−αijΓνi−αij+1∫0xx−zνi−αijKijx,zGijz,ujz−K~ijx,zGijz,u~jzdzL2≤∑j=1mθijΨνi−αijMKijηjiujz−u~jzL2+NGijΛKij≤∑j=1mθijΨνi−αijMKijηjiΩNj+NGijΛKij,where
(76)Ψνi−αij=Γ1−αijΓνi−αij+1×4Γ7/2Γ2νi−αij+3/2Γ2νi−αij+5−4Γ5/2Γ2νi−αij+3/2Γ2νi−αij+4+Γ3/2Γ2νi−αij+3/2Γ2νi−αij+31/2.
Here,Δji and γji can be introduced as below:
(77)Δji=Ψνiξji+Ψνi−αijθijMKijηji,(78)γji=Ψνi−αijθijNGijΛKij.
From Equations (72)-(78), we can get the following upper bound for rix:
(79)rixL2≤ΩNi+∑j=1mΔjiΩNj+∑j=1mγji,i=1,2,⋯,m.
By adding the abovem inequalities, we have
(80)∑i=1mrixL2≤∑i=1mΩNi+∑i=1m∑j=1mΔjiΩNj+∑i=1m∑j=1mγji=ΓN∗.
We define the vectorR as
(81)R=r1L2,r2L2,⋯,rmL2T.Thus, we have
(82)R1=∑i=1mriL2≤ΓN∗.
Again, we consider Equation (70). So we have
(83)uix−u~ixL2≤rixL2+∑j=1mΔjiuix−u~ixL2+∑j=1mγji,i=1,2,⋯,m.
Adding the abovem inequalities leads to the following inequality:
(84)∑i=1muix−u~ixL2≤∑i=1mrixL2+∑i=1m∑j=1mΔjiuix−u~ixL2+∑i=1m∑j=1mγji=R1+Δ11u1x−u~1xL2+⋯+Δm1umx−u~mxL2+Δ12u1x−u~1xL2+⋯+Δm2umx−u~mxL2+⋯+Δ1mu1x−u~1xL2+⋯+Δmmumx−u~mxL2+∑i=1m∑j=1mγji=R1+Δ11+Δ12+⋯+Δ1mu1x−u~1xL2+⋯+Δm1+Δm2+⋯+Δmmumx−u~mxL2+∑i=1m∑j=1mγji.
By defining the following quantities(85)ΔN∗=Max∑i=1mΔji,j=1,⋯,m,E=u1x−u~1xL2,u2x−u~2xL2,⋯,umx−u~mxL2T,we obtain the following upper bound for the method error:
(86)∑i=1muix−u~ixL2≤R1+ΔN∗∑i=1muix−u~ixL2+∑i=1m∑j=1mγji,E1≤ΓN∗+∑i=1m∑j=1mγji1−ΔN∗,0<ΔN∗<1.


## 6. Numerical Applications

In this section, three linear and nonlinear examples are presented to illustrate the practical implementation of our numerical method. Also, the comparison of results obtained from the proposed method with those of the other methods is shown. All calculations are done with mathematical software Maple 18.Example 1.
Let us consider the following linear system of WSFIDE:(87)Dν1u1x+u1x+∫0xu1zx−z1/2dz−12∫0xxu2zx−z1/2dz=f1x,Dν2u2x+u2x+13∫0xx2zu1zx−z1/2dz+13∫0xu2zx−z1/2dz=f2x,where
(88)f1x=2x+x2+25x5/2,f2x=1+x+32105x11/2+49x3/2,and 0<νi≤1,i=1,2, and the initial conditions are u10=u20=0. The exact solutions are u1x=x2 and u2x=x if ν1=ν2=1. According to the procedure presented in Section 4, we reach the following approximations:
(89)Dν1u1x≈STxF1,u1x≈STxPν1TF1=STxU1,Dν2u2x≈STxF2,u2x≈STxPν2TF2=STxU2,x≈STxK12Sz,x2z≈STxK21Sz,∫0xu1zx−z1/2dz≈∫0xU1TSzx−z1/2dz=U1T∫0xSzx−z1/2dz≈U1TI1/2Sx,∫0xxu2zx−z1/2dz≈∫0xSTxK12SzSTzU2x−z1/2dz=STxK12∫0xSzSTzU2x−z1/2dz≈STxK12U~2∫0xSzx−z1/2dz≈STxK12U~2I1/2Sx,∫0xx2zu1zx−z1/2dz≈∫0xSTxK21SzSTzU1x−z1/2dz=STxK21∫0xSzSTzU1x−z1/2dz=STxK21U~1∫0xSzx−z1/2dz=STxK21U~1I1/2Sx,∫0xu2zx−z1/2dz≈∫0xU2TSzx−z1/2dz≈U2T∫0xSzx−z1/2dz≈U2TI1/2Sx,where U~1 and U~2 are the operational matrices of the product, corresponding to the vectors U1 and U2, respectively. Setting above approximations in system (87) leads to the following linear algebraic system.
(90)STxF1+STxU1+U1TI1/2Sx−12STxK12U~2I1/2Sx−f1x≈0,STxF2+STxU2+13STxK21U~1I1/2Sx+13U2TI1/2Sx−f2x≈0.Figures1 and 2 show a comparison between the exact and numerical solutions and absolute error functions for ν1=ν2=1 and N=20. The maximum absolute errors are listed in Table 1 in versus of N. It can be seen that our proposed method shows good consistency between the numerical results and analytic solutions and also this method can achieve a higher convergence result when N increases. Figure 3 shows the behavior of the numerical solutions for N=20 and νi=0.8,0.9,1,i=1,2. As seen from Figure 3, as νi⟶1, the approximate solutions are u~ix⟶uix for i=1,2.Figure 1
The graphs of the exact and approximate solutions and the absolute error function ofu~1x for N=20 and ν1=ν2=1 of Example 1.Figure 2
The graphs of the exact and approximate solutions and the absolute error function ofu~2x for N=20 and ν1=ν2=1 of Example 1.Table 1
Maximum absolute errors for different values ofN for Example 1.
Nu1xu2xCPU time32.7867×10−53.0371×10−53.01083.7450×10−64.7250×10−68.642112.3054×10−77.7804×10−718.330162.1369×10−73.0206×10−768.048208.0541×10−81.1680×10−7166.843Figure 3
Some illustrations for the exact and approximate solutions using different values ofνi: νi=0.8,0.9,1,i=1,2 in Example 1.Here, we calculate a numerical error bound for the first example. This result could confirm the correctness of the analytical error bound. This validation could be accomplished similarly for the other cases, but the calculations are long.(91)u1x−u~1xL2=8.93650018×10−9,u2x−u~2xL2=1.81941246×10−8,E=8.93650018×10−9,1.81941246×10−8T,E1=2.71305687×10−8,ξ11=0.4,ξ21=0.3,ξ12=0.8,ξ22=0.62,η11=η12=η21=η22=1,Ψν1=Ψν2=0.31332853,Ψν1−α11=Ψν1−α12=0.31332834,Ψν2−α21=Ψν2−α22=0.31332834,Δ11=0.43865995,Δ12=0.25392268,Δ21=0.12405832,Δ22=0.29870654,ΔN∗=Max∑i=12Δji,j=1,2=0.69258263,ΛK11=ΛK12=ΛK21=ΛK22=14.31641021,γ11=0.72969033,γ21=−0.43034788,γ12=0.24323011,γ22=0.28689859,ΓN∗=2.71305877×10−8,ΓN∗+∑i=1m∑j=1mγji1−ΔN∗=2.69819226,E1=2.71305687×10−8≤2.69819226.Example 2.
In this example, consider the following system of linear WSFIDEs [13, 18]:
(92)Dν1u1x−u3x−∫0xxzu1zx−z1/2dz−∫0xu2zx−z1/2dz=f1x,Dν2u2x−u1x−∫0xu2zx−z1/3dz−∫0xu3zx−z1/3dz=f2x,Dν3u3x−u3x−∫0xu1zx−z1/4dz−∫0xx2zu2zx−z1/4dz=f3x,where 0≤x≤1 and
(93)f1x=2x−1−x3+1615x7/2−1615x5/2−3235x9/2,f2x=2x−xx−1−2740x8/3−243440x11/3,f3x=3x2−x3+1621x7/4−128231x11/4−5121155x23/4,and 0<νi≤1,i=1,2,3, the initial conditions are u10=u20=u30=0. The exact solutions are u1x=xx−1,u2x=x2,u3x=x3 if ν1=ν2=ν3=1. By assuming that Dν1u1x≈STxF1, Dν2u2x≈STxF2, and Dν3u3x≈STxF3 and using proper operational matrices, the following approximations are obtained:
(94)u1x≈STxPν1TF1=STxU1,u2x≈STxPν2TF2=STxU2,u3x≈STxPν3TF3=STxU3,xz≈STxK11Sz,x2z≈STxK32Sz,∫0xxzu1zx−z1/2dz+∫0xu2zx−z1/2dz≈STxK11U~1I1/2Sx+U2TI1/2Sx,∫0xu2zx−z1/3dz+∫0xu3zx−z1/3dz≈U2TI1/3Sx+U3TI1/3Sx,∫0xu1zx−z1/4dz+∫0xx2zu2zx−z1/4dz≈U1TI1/4Sx+STxK32U~2I1/4Sx,where U~1 and U~2 are the operational matrices of product, corresponding to the vectors U1 and U2, respectively. By substituting the above approximations into Equation (92), we achieve the following algebraic system
(95)STxF1−STxU3−STxK11U~1I1/2Sx−U2TI1/2Sx−f1x≈0,STxF2−STxU1−U2TI1/3Sx−U3TI1/3Sx−f2x≈0,STxF3−STxU3−U1TI1/4Sx−STxK32U~2I1/4Sx−f3x≈0.
Solving the above system by the collocation method forN=9, we can determine the unknown vectors Fi,i=1,2,3. Table 2 displays the maximum absolute errors for various N. The data in this table show that the numerical solutions get close to the analytical solutions with the increase of values of N. Table 3 shows a comparison between the SSKCP collocation and Jacobi collocation methods. This table shows that the results are approximately the same as reported by [18]. Figures 4–6 display a graphical comparison between approximate solutions and exact solutions and also absolute error functions for N=9. Figure 7 shows the behavior of the numerical solutions for N=9 and νi=0.85,0.90,0.95,1,i=1,2,3. From Figure 7, it can be seen that as νi⟶1, the approximate solutions are u~ix⟶uix,i=1,2,3.Example 3.
Consider the following nonlinear system of WSFIDEs [13, 18]:
(96)D3/4u1x−u2x−∫0xu1zx−z1/2dz−∫0xu22zx−z1/2dz=f1x,D1/2u2x−u13x−∫0xu1zx−z2/3dz−∫0xu2zx−z2/3dz=f2x,
where0≤x≤1 and
(97)f1x=2Γ9/4x5/4−x3−Γ3Γ1/2Γ7/2x5/2−Γ7Γ1/2Γ5/2x13/2,f2x=Γ4Γ7/2x5/2−x6−Γ3Γ1/3Γ10/3x7/3−Γ4Γ1/3Γ13/3x10/3.
The initial conditions areu10=u20=0 and the exact solutions for this example are u1x=x2 and u2x=x3. We use the following approximations and convert the system (96) into an algebraic system, which will be described below:
(98)D3/4u1x≈STxF1,u1x≈STxP3/4TF1=STxU1,D1/2u2x≈STxF2,u2x≈STxP1/2TF2=STxU2,u13x≈STxU3,U3≈U~12TU1,u22x≈STxU4,U4≈U~2TU2,∫0xu1zx−z1/2dz+∫0xu22zx−z1/2dz≈U1TI1/2Sx+U4TI1/2Sx,∫0xu1zx−z2/3dz+∫0xu2zx−z2/3dz≈U1TI2/3Sx+U2TI2/3Sx,where U~1 and U~2 are operational matrices of the product, corresponding to the vectors U1 and U2, respectively. So we have
(99)STxF1−STxU2−U1TI1/2Sx−U4TI1/2Sx−f1x≈0,STxF2−STxU3−U1TI2/3Sx−U2TI2/3Sx−f2x≈0.
Table4 shows a comparison between the absolute errors of the SSKCP collocation method and the Jacobi collocation method at equally spaced points xi=0.1i,i=0,1,⋯,10 for N=12, which shows that the absolute errors of the Jacobi collocation method are less than the presented method, but with increasing N, the errors of the proposed method decrease. Maximum absolute errors for different values of N and numerical results for N=25 reported in Tables 5 and 6 confirm that the results are close to those reported by [18]. Figures 8 and 9 show the comparison between the numerical results and the exact solutions and also the absolute error functions of u1x,u2x for N=25, respectively.Table 2
Maximum absolute errors for different values ofN for Example 2.
Nu1xu2xu3xCPU time48.1514×10−56.7763×10−51.0743×10−45.38292.3334×10−61.8107×10−63.3112×10−617.457126.8859×10−75.4537×10−71.1140×10−636.941171.2296×10−79.5306×10−82.0253×10−7115.316205.9278×10−84.6454×10−81.0049×10−7216.638Table 3
Maximum absolute errors obtained by SSKCP collocation and Jacobi collocation methods forN=9, α=−1/2, and β=1/2 in Example 2.
Methodu1xu2xu3xSSKCP collocation2.3334×10−61.8107×10−63.3112×10−6Jacobi collocation [18]2.0034×10−61.8376×10−62.0763×10−6Figure 4
The graphs of the exact and approximate solutions and the absolute error function ofu~1x for N=9 and ν1=ν2=ν3=1 of Example 2.Figure 5
The graphs of the exact and approximate solutions and the absolute error function ofu~2x for N=9 and ν1=ν2=ν3=1 of Example 2.Figure 6
The graphs of the exact and approximate solutions and the absolute error function ofu~3x for N=9 and ν1=ν2=ν3=1 of Example 2.Figure 7
Some illustrations for the exact and approximate solutions using different values ofνi: νi=0.85,0.90,0.95,1,i=1,2,3 for N=9 in Example 2.Table 4
Comparison between absolute errors of SSKCP collocation and Jacobi collocation methods at equally spaced points forN=12 in Example 3.
xiSSKCP collocation methodJacobi collocation method [18]Erroru1Erroru2Erroru1Erroru20.01.7182×10−59.8329×10−79.2840×10−66.2830×10−70.15.5657×10−63.3268×10−61.5202×10−61.6472×10−70.25.5935×10−66.9555×10−67.0504×10−79.8923×10−80.38.1095×10−61.2777×10−58.7678×10−72.9284×10−70.41.1355×10−52.2448×10−54.8193×10−74.5120×10−70.52.0566×10−53.9797×10−58.8367×10−77.7271×10−70.63.2387×10−57.0527×10−53.7500×10−71.3807×10−60.75.7063×10−51.2877×10−41.0664×10−62.4661×10−60.81.1121×10−42.5191×10−42.2537×10−64.8902×10−60.92.3722×10−45.4688×10−44.7116×10−61.0475×10−51.05.7984×10−41.3758×10−31.2704×10−52.6197×10−5Table 5
Maximum absolute errors obtained by SSKCP collocation method for different values ofN in Example 3.
NErroru1Erroru2CPU time82.8334×10−36.7948×10−35.772125.7819×10−41.3597×10−323.587161.8841×10−44.4719×10−497.048207.8006×10−51.8536×10−4164.409252.9931×10−56.1196×10−5483.681Table 6
Values of absolute errors at equally spaced points forN=25 in Example 3.
xiErroru1Erroru20.01.3659×10−63.2865×10−80.12.6535×10−71.7914×10−70.23.1410×10−73.8266×10−70.34.3165×10−77.0524×10−70.46.3587×10−71.2387×10−60.59.6524×10−72.1315×10−60.61.5052×10−63.6310×10−60.72.8039×10−66.5226×10−60.85.5119×10−61.2576×10−60.91.1848×10−52.7085×10−51.02.8779×10−56.9756×10−5Figure 8
The graphs of the exact and approximate solutions and the absolute error function ofu~1x for N=25 of Example 3.Figure 9
The graphs of the exact and approximate solutions and the absolute error function ofu~2x for N=25 of Example 3.


## 7. Conclusion

In this paper, the shifted sixth-kind Chebyshev polynomials together with the collocation method were used to solve a class of the system of fractional integro-differential equations with weakly singular kernels. For this purpose, the integral and product operational matrices were calculated, and using the obtained approximations, the original system of equations was transformed into a corresponding linear and nonlinear system of algebraic equations that are easier for solving. Choosing an appropriate value ofN, each of algebraic equations was collocated in the roots of SN+1∗x, and finally, for obtaining the unknown vectors Fi,i=1,2,⋯,m, an algebraic system involving mN+1 algebraic equations was solved. To eliminate the singularity of the kernels of the equations under study, an operational matrix was derived. Also, an error bound was determined for the proposed method. To show the ability and efficiency of the proposed method, three examples were presented and the maximum absolute errors were calculated for different N, and graphs of the absolute error functions and numerical solutions were plotted and numerical results showed a good agreement between the approximate and exact solutions. When the order of the fractional derivative ν was uncertain, the numerical solutions for the various values of ν, 0<ν≤1, were approached to the exact solutions as ν⟶1. The comparison of the proposed method with the Jacobi collocation method [18] showed good implementation of SSKCP collocation method for solving a system of linear fractional integro-differential equations with weakly singular kernels and for a system of nonlinear WSFIDEs; the error decreased with increasing N. CPU times were computed for all examples. According to the numerical results which were obtained in the relevant instances and compared with exact solutions and those obtained from the Jacobi collocation method, it can be concluded that the SSKCP collocation method is very helpful to look for approximate solutions of a system of WSFIDEs. This method can be applied to the linear and nonlinear systems of fractional-order Volterra-Fredholm integro-differential equations with weakly singular kernels, but additional operational matrices are required.



---

*Source: 9087359-2022-10-17.xml*