# Fractional-Order Chelyshkov Collocation Method for Solving Systems of Fractional Differential Equations

Authors: A. I. Ahmed; T. A. Al-Ahmary
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper presents an efficient method for solving systems of fractional differential equations by using the fractional-order Chelyshkov functions (FCHFs). The fractional derivative and the fractional integral are considered in the Caputo sense and the Riemann–Liouville sense, respectively. The proposed method is based on using the operational matrix of fractional integration for FCHFs, together with the spectral collocation method to transform the system of fractional differential equations into a system of algebraic equations. The convergence of the presented method is demonstrated. The performance of the presented method is tested through various examples of systems of fractional differential equations. A comparison with some existing methods shows that the presented method can be successfully used to solve systems of fractional differential equations with accuracy and efficiency.


---

## Body

## 1. Introduction

Fractional differential equations arise in many applications, such as fluid mechanics, chemistry, physics, and engineering [1–3]. Therefore, many researchers have focused their attention on finding solutions to fractional differential equations. The existence and uniqueness of solutions to fractional differential equations have been studied in several papers [3,4]. Most of these equations do not have exact analytical solutions, so many studies have been devoted to develop new methods to find numerical or approximate solutions, such as the collocation method [5], Tau method [6], B-Spline operational matrix method [7], Galerkin method [8], space–time Chebyshev collocation method [9], Euler functions method [10], differential Transformation method [11], and homotopy analysis method [12].In general, the discretization methods that give numerical solutions to differential equations can be classified into spectral methods, finite differences methods, and finite element methods. Spectral methods are more preferred than others since they are accurate for a given number of unknowns [13] and exhibit exponential rates of convergence/spectral accuracy for smooth problems [14]. The idea of spectral methods is to approximate the solution of a differential equation as an expansion of certain basic functions and then to determine the coefficients in that expansion in order to satisfy the differential equation. The most applied spectral versions are the Galerkin, Tau, and collocation methods [15, 16].Orthogonal functions have been considered as a powerful tool to find and interpret the solutions to several types of differential equations. These polynomials simplify the treatment of differential equations and can be generated by applying the Gram-Schmidt orthonormalization procedure. Many researchers have used orthogonal functions to solve differential equations, such as Chebyshev wavelets [17], fractional-order Legendre wavelets [6], improved block–pulse basis functions [18], shifted second kind Chebyshev polynomials [19], shifted Jacobi polynomials [20], etc. Orthogonal Chelyshkov polynomials were introduced in [21] and used for solving several types of differential equations [22–24].The common and efficient methods that find solutions to fractional differential equations are based on employing series expansion of the form∑k=0nqkxkα,α>0 such as Adomian’s decomposition method [25] and He’s variational iteration method [26]. Thus, it may be useful to use orthogonal functions of fractional order to solve fractional differential equations to increase efficiency.Systems of fractional differential equations represent an appropriate tool in modeling many processes, including the COVID-19 pandemic [27, 28]. The fractional model is preferred because it dominates the effects by adding some parameters like memory and correlations which lead to better analysis and prediction [29]. Various numerical and analytical solutions have been presented for some nonlinear systems of fractional differential equations [29–39]. However, the solution of systems of fractional differential equations is still in progress and far from completeness.In this paper, we consider the following general form of the system of nonlinear fractional differential equations:(1)Dαixit=fit,x1,x2,…,xm,xij0=λij,i=1,2,…,m,j=0,1,…,αi−1,where αi is the smallest integer greater than or equal to αi and Dαi denotes the Caputo fractional derivative of order αi. The existence, uniqueness, and stability of solutions to the system of nonlinear fractional differential equations (1) have been discussed in [40, 41]. Many methods have been introduced for solving this system [31–33, 38]. To the best knowledge of the author, there is no article about solving system (1) using FCHFs. Therefore, the motivation of this paper is to provide a method for solving system of nonlinear fractional differential equations (1) by using FCHFs. To achieve this aim, we construct the FCHFs based on Chelyshkov polynomials and then derive the operational matrix of fractional integration for FCHFs and use it to present a new technique for solving system (1). Our method possesses some beneficial properties as follows:(i)
The present method can be applied without using linearization or perturbation unlike some other methods such as the homotopy perturbation method.(ii)
It uses orthogonal functions of fractional order which overcome the drawbacks of the poor rate of convergence and the required high number of basic functions that occur when solutions to the system of fractional differential equations contain terms with fractional powers.(iii)
It reduces the system of fractional differential equations to a system of algebraic equations.(iv)
The obtained numerical results show that the presented method is accurate, helpful, and easy to implement.The paper is arranged as follows: In the next section, we introduce some definitions and preliminaries of fractional calculus. In Section3, we present the FCHFs and their properties and derive the operational matrix of fractional integration. Section 4 demonstrates the convergence of the fractional-order Chelyshkov method. In Section 5, the proposed technique for solving systems of nonlinear fractional differential equations is described. Some examples are shown in Section 6 with satisfactory numerical results. Finally, some conclusions are drawn in Section 7.


## 2. Preliminaries

This section presents some definitions and properties of fractional integration and differentiation which will be used in the subsequent sections [42,43]. The literature containany different definitions of fractional integration of order α≥0. The Riemann–Liouville definition is the more widely used, which can be given as follows:Definition 1.
The Riemann–Liouville fractional integral operator of orderα≥0 is defined by(2)Iαxt=1Γα∫0tt−τα−1xτdτ,α>0,t>0,xt,α=0.WhereΓα is the gamma function. For the operator Iα we can write(3)Iαtμ=Γμ+1Γμ+α+1tμ+α,μ>−1,(4)Iα1Iα2xt=Iα1+α2xt=Iα2Iα1xt,α1,α2≥0,(5)Iαη1x1t+η2x2t=η1Iαx1t+η2Iαx2t,where η1 and η2 are constants.Definition 2.
The fractional derivative of orderα>0 in the Caputo sense is given by the following equation:(6)Dαxt=Iα−αDαxt=1Γα−α∫0tt−τα−α−1xατdτ,α−1<α≤α,t>0.For Caputo fractional derivative operatorDα, we have(7)Dαtμ=0,μ∈N0,μ<α,Γμ+1Γμ−α+1tμ−α,μ∈N0,μ≥αorμ∉N,μ>α,(8)Dαη=0,ηisaconstant,(9)DαIαxt=xt,(10)IαDαxt=xt−∑j=0α−1xj0+j!tj,t>0,(11)Dαρ1x1t+ρ2x2t=ρ1Dαx1t+ρ2Dαx2t,where α denotes the largest integer less than or equal to α, N=1,2,…, N0=0,1,2,… and ρ1, ρ2 are constants.Definition 3.
(Generalized Taylor’s formula). Assume thatDkμxt∈C0,1 for k=0,1,…,n+1 and 0<μ≤1, then(12)xt=∑l=0ntlμΓlμ+1Dlμx0++tn+1μΓn+1μ+1Dn+1μxτ,0<τ≤t,∀t∈0,1,where Dlμ=DμDμ⋯Dμ︸ltimes.In case ofμ=1, the generalized Taylor’s formula reduces to the classical Taylor’s formula.


## 3. Fractional Chelyshkov Functions and Their Properties

### 3.1. Fractional Chelyshkov Functions

The family of Chelyshkov polynomialsFn=Fnii=0n is a set of orthogonal polynomials over the interval 0,1 with respect to the weight function ωτ=1 and is given as follows:(13)Fniτ=∑l=in−1l−in−il−in+l+1n−iτl,i=0,1,…n.These polynomials satisfy the orthogonality relationship(14)∫01FnsτFnrτdτ=0,s≠r,1s+r+1,s=r,s,r=0,1,…,n.From equations (13), we can say that every polynomial in the family Fn has degree n, this is the main difference between Chelyshkov polynomials and the other families of orthogonal polynomials. The relationship between Chelyshkov polynomials Fniτ and Jacobi polynomials Jmα,βτ, where α,β>−1 and m≥0 can be expressed by(15)Fniτ=−1n−iτiJn−i0,2i+12τ−1,i=0,1,…,n.To construct the fractional-order Chelyshkov functions (FCHFs), we substitutetμ, μ>0 for τ in equations (13). Let the FCHFs Fnitμ be denoted by Fniμt. Then, Fniμt can take the following form:(16)Fniμt=∑l=in−1l−in−il−in+l+1n−itlμ,i=0,1,…n.The FCHFs are orthogonal over the interval0,1 with respect to the weight function ωμt=tμ−1 and satisfy the orthogonality condition(17)∫01FnsμtFnrμtωμtdt=0,s≠r,1s+r+1μ,s=r,s,r=0,1,…,n.Let us introduce some notation.(18)Fμt=Fn0μt,Fn1μt,…,FnnμtT,(19)Tμt=1,tμ,t2μ,…,tnμT,(20)E=eili,l=0n,eileil=0,0≤l<i,−1l−in−il−in+l+1n−i,i≤l≤n,where Fμt denotes the FCHFs vector and can be written in the form(21)Fμt=ETμt.


### 3.2. Approximation of Functions

Let the weighted spaceLωμ2υ, υ=0,1 be defined as follows:(22)Lωμ2υ=ξ:υ⟶ℝ;ξis measurable onυ&∫01ξt2ωμtdt<∞.and its inner product and norm are given as follows:(23)ξ,ζωμ=∫01ξtζtωμtdt,(24)ξωμ=ξ,ξωμ1/2.LetFn=spanFn0μt,Fn1μt,…,Fnnμt denotes a finite dimensional and closed subspace of Lωμ2υ, then for each xt∈Lωμ2υ there is a unique best approximation xnt∈Fn such that(25)x−xnωμ≤x−σωμfor allσ∈Fn.Since the vector spaceFn contains the function xnt, then it can be approximated in terms of FCHFs as(26)xnt=∑i=0nciFniμt=CTFμt,where the coefficients ci are given by(27)ci=x,FniμωμFniμ,Fniμωμ,i=0,1,…,n,(28)C=c0,c1,…,cnT.


### 3.3. Operational Matrix of Fractional Integration

In the following section, we aim to build the operational matrix of fractional integration of FCHFs, which will be used to reduce the system of fractional differential equations into a set of algebraic equations.Theorem 1.
LetFμt be the FCHFs vector defined by equation (21), then(29)IαFμt≃PαFμt,where Pα=piki,k=0n is the n+1×n+1 operational matrix of fractional integration of order α>0 in the Riemann–Liouville sense and its elements take the following form(30)pik=∑l=in∑s=kn−1l+s−i−kμ2k+1Γlμ+1μl+s+1+αΓlμ+α+1n−il−in+l+1n−in−ks−kn+s+1n−k.Proof.
The integration of equation (21) gives (31)IαFμt=IαETμt=EIαTμt.
From equations (3) and (19), we can write(32)IαTμt=1Γα+1tα,Γμ+1Γμ+α+1tμ+α,Γ2μ+1Γ2μ+α+1t2μ+α,…,Γnμ+1Γnμ+α+1tnμ+αT=E¯T¯μt,where(33)E¯=diag1Γα+1,Γμ+1Γμ+α+1,Γ2μ+1Γ2μ+α+1,…,Γnμ+1Γnμ+α+1,(34)T¯μt=tα,tμ+α,t2μ+α,…,tnμ+αT.
By approximating the components of the vectorT¯μt in terms of FCHFs, we get(35)tlμ+α≃∑k=0ne^lkFnkμt,l=0,1,…,n,where(36)e^lk=μ2k+1∫01tlμ+αFnkμtωμtdt=μ2k+1∑s=kn−1s−kn−ks−kn+s+1n−k∫01tl+s+1μ+α−1dt=μ2k+1∑s=kn−1s−kμl+s+1+αn−ks−kn+s+1n−k,
thus(37)T¯μt≃E^Fμt,where the elements of the matrix E^=e^lkl,k=0n are given by equation (36). Employing equations (31), (32), and (37), leads to(38)IαFμt≃EE¯E^Fμt,and by multiplying the matrices E, E¯, and E^ described by equations (20), (33), and (36), we get(39)IαFμt≃PαFμt.Remark 1.
Whenμ=1, we get the same operational matrix of fractional integration in [23], which means that the operational matrix of fractional integration in [23] is a special case of the operational matrix of fractional integration in the present paper.


## 3.1. Fractional Chelyshkov Functions

The family of Chelyshkov polynomialsFn=Fnii=0n is a set of orthogonal polynomials over the interval 0,1 with respect to the weight function ωτ=1 and is given as follows:(13)Fniτ=∑l=in−1l−in−il−in+l+1n−iτl,i=0,1,…n.These polynomials satisfy the orthogonality relationship(14)∫01FnsτFnrτdτ=0,s≠r,1s+r+1,s=r,s,r=0,1,…,n.From equations (13), we can say that every polynomial in the family Fn has degree n, this is the main difference between Chelyshkov polynomials and the other families of orthogonal polynomials. The relationship between Chelyshkov polynomials Fniτ and Jacobi polynomials Jmα,βτ, where α,β>−1 and m≥0 can be expressed by(15)Fniτ=−1n−iτiJn−i0,2i+12τ−1,i=0,1,…,n.To construct the fractional-order Chelyshkov functions (FCHFs), we substitutetμ, μ>0 for τ in equations (13). Let the FCHFs Fnitμ be denoted by Fniμt. Then, Fniμt can take the following form:(16)Fniμt=∑l=in−1l−in−il−in+l+1n−itlμ,i=0,1,…n.The FCHFs are orthogonal over the interval0,1 with respect to the weight function ωμt=tμ−1 and satisfy the orthogonality condition(17)∫01FnsμtFnrμtωμtdt=0,s≠r,1s+r+1μ,s=r,s,r=0,1,…,n.Let us introduce some notation.(18)Fμt=Fn0μt,Fn1μt,…,FnnμtT,(19)Tμt=1,tμ,t2μ,…,tnμT,(20)E=eili,l=0n,eileil=0,0≤l<i,−1l−in−il−in+l+1n−i,i≤l≤n,where Fμt denotes the FCHFs vector and can be written in the form(21)Fμt=ETμt.


## 3.2. Approximation of Functions

Let the weighted spaceLωμ2υ, υ=0,1 be defined as follows:(22)Lωμ2υ=ξ:υ⟶ℝ;ξis measurable onυ&∫01ξt2ωμtdt<∞.and its inner product and norm are given as follows:(23)ξ,ζωμ=∫01ξtζtωμtdt,(24)ξωμ=ξ,ξωμ1/2.LetFn=spanFn0μt,Fn1μt,…,Fnnμt denotes a finite dimensional and closed subspace of Lωμ2υ, then for each xt∈Lωμ2υ there is a unique best approximation xnt∈Fn such that(25)x−xnωμ≤x−σωμfor allσ∈Fn.Since the vector spaceFn contains the function xnt, then it can be approximated in terms of FCHFs as(26)xnt=∑i=0nciFniμt=CTFμt,where the coefficients ci are given by(27)ci=x,FniμωμFniμ,Fniμωμ,i=0,1,…,n,(28)C=c0,c1,…,cnT.


## 3.3. Operational Matrix of Fractional Integration

In the following section, we aim to build the operational matrix of fractional integration of FCHFs, which will be used to reduce the system of fractional differential equations into a set of algebraic equations.Theorem 1.
LetFμt be the FCHFs vector defined by equation (21), then(29)IαFμt≃PαFμt,where Pα=piki,k=0n is the n+1×n+1 operational matrix of fractional integration of order α>0 in the Riemann–Liouville sense and its elements take the following form(30)pik=∑l=in∑s=kn−1l+s−i−kμ2k+1Γlμ+1μl+s+1+αΓlμ+α+1n−il−in+l+1n−in−ks−kn+s+1n−k.Proof.
The integration of equation (21) gives (31)IαFμt=IαETμt=EIαTμt.
From equations (3) and (19), we can write(32)IαTμt=1Γα+1tα,Γμ+1Γμ+α+1tμ+α,Γ2μ+1Γ2μ+α+1t2μ+α,…,Γnμ+1Γnμ+α+1tnμ+αT=E¯T¯μt,where(33)E¯=diag1Γα+1,Γμ+1Γμ+α+1,Γ2μ+1Γ2μ+α+1,…,Γnμ+1Γnμ+α+1,(34)T¯μt=tα,tμ+α,t2μ+α,…,tnμ+αT.
By approximating the components of the vectorT¯μt in terms of FCHFs, we get(35)tlμ+α≃∑k=0ne^lkFnkμt,l=0,1,…,n,where(36)e^lk=μ2k+1∫01tlμ+αFnkμtωμtdt=μ2k+1∑s=kn−1s−kn−ks−kn+s+1n−k∫01tl+s+1μ+α−1dt=μ2k+1∑s=kn−1s−kμl+s+1+αn−ks−kn+s+1n−k,
thus(37)T¯μt≃E^Fμt,where the elements of the matrix E^=e^lkl,k=0n are given by equation (36). Employing equations (31), (32), and (37), leads to(38)IαFμt≃EE¯E^Fμt,and by multiplying the matrices E, E¯, and E^ described by equations (20), (33), and (36), we get(39)IαFμt≃PαFμt.Remark 1.
Whenμ=1, we get the same operational matrix of fractional integration in [23], which means that the operational matrix of fractional integration in [23] is a special case of the operational matrix of fractional integration in the present paper.


## 4. Convergence Analysis

In this section, we prove the convergence of the FCHFs with respect to theLωμ2-norm and derive an upper bound of the error for the operational matrix of fractional integration.Theorem 2.
Suppose thatDiμxt∈C0,1 for i=0,1,…,n+1. If xnt is the best approximation to x out of Fn, then the error bound can be given by(40)xt−xntωμ≤MμΓn+1μ+12n+3μ,where Mμ=supDn+1μxt,t∈0,1.Proof.
In view of Definition3 and by supposing that ρt is the generalized Taylor’s formula of xt we have(41)ρt=∑l=0ntlμΓlμ+1Dlμx0+,t∈0,1,with an error bound given as follows:(42)xt−ρt=tn+1μΓn+1μ+1Dn+1μxτ,0<τ≤t≤Mμtn+1μΓn+1μ+1.
SinceFn contains both ρt and xnt, and xn is the best approximation to x from Fn then(43)xt−xntωμ2≤xt−ρtωμ2=∫01xt−ρt2ωμtdt≤∫01Mμtn+1μΓn+1μ+12ωμtdt=Mμ2Γn+1μ+12∫01t2n+3μ−1dt=Mμ2Γn+1μ+122n+3μ,which leads to the desired result.Theorem2 shows that for a fixed μ, the best approximation xn converges to x as n⟶∞. Now, we find an upper bound for the error of Pα. The error vector Eα of the operational matrix Pα is given by the following equation:(44)Eα=IαFμt−PαFμt,Eα=ε0,ε1,…,εnT.The functiontlμ+α was approximated in equation (35) as ∑k=0ne^lkFnkμt, so from Theorem 2, we can write(45)tlμ+α−∑k=0ne^lkFnkμtωμ≤M¯μΓn+1μ+12n+3μ,where(46)M¯μ=supDn+1μtlμ+α,t∈0,1.Then, by using equations (16), (30) and (45), we get(47)εiωμ=IαFniμt−∑k=0npikFnkμtωμ=∑l=in−1l−iΓlμ+1Γlμ+α+1n−il−in+l+1n−itlμ+α−∑k=0n∑l=in−1l−ie^lkΓlμ+1Γlμ+α+1n−il−in+l+1n−iFnkμtωμ≤∑l=inΓlμ+1Γlμ+α+1n−il−in+l+1n−iM¯μΓn+1μ+12n+3μ,i=0,1,…,n.From Inequality (47), it can be seen that by increasing the number of FCHFs, the error vector Eα tends to zero.


## 5. Solution Method of Systems of Fractional Differential Equations

In this section, the fractional Chelyshkov matrix method which is based on computing Chelyshkov coefficients is applied to solve the system of fractional differential equations (1).First, we approximateDαixit,i=1,2,…,m as follows:(48)Dαixit≃∑k=0ncikFnkμt=CiTFμt,where Ci,i=1,2,…,m are unknown vectors. By integrating equations (48) and applying equation (27), we obtain(49)IαiDαixit≃CiTIαiFμt≃CiTPαiFμt,i=1,2,…,m.From equations (10) and (49), one can write(50)xit−∑j=0αi−1xij0j!tj≃CiTPαiFμt,i=1,2,…,m.Employing the initial conditions gives(51)xit≃CiTPαiFμt+∑j=0αi−1λijj!tj,i=1,2,…,m.Now, we approximate∑j=0αi−1λij/j!tj for i=1,2,…,m in terms of FCHFs as follows:(52)∑j=0αi−1λijj!tj≃∑r=0nairFnrμt=AiTFμt,where the vectors Ai=ai0,ai1,…,ainT,i=1,2,…,m are given by(53)air=μ2r+1∑j=0αi−1λijj!∫01tjFnrμtωμtdt=μ2r+1∑j=0αi−1∑l=rn−1l−rλijj!n−rl−rn+l+1n−r∫01tμl+1+j−1dt=μ2r+1∑j=0αi−1∑l=rn−1l−rλijμl+1+jj!n−rl−rn+l+1n−r,i=1,2,…,m,r=0,1,…,n.Therefore, equation (51) becomes(54)xit≃CiTPαi+AiTFμt,i=1,2,…,m.By substituting the approximations forDαixit and xit given by equations (51) and (54) in equations (1), the nonlinear system of fractional differential equations (1) is reduced into the following system of algebraic equations:(55)CiTFμt−fit,C1TPα1+A1TFμt,…,CmTPαm+AmTFμt≃0,i=1,2,…,m.Collocating equations (55) at the equidistant nodes tk=k/n,k=0,1,…,n yields(56)CiTFμtk−fitk,C1TPα1+A1TFμtk,…,CmTPαm+AmTFμtk≃0,i=1,2,…,m.Equation (56) is a system of mn+1 nonlinear algebraic equations which can be solved to find the unknown vectors Ci,i=1,2,…,m. Thus, approximate solutions of xit,i=1,2,…,m can be determined.


## 6. Numerical Examples

To show the efficiency and accuracy of the presented method, we use it to solve some linear and nonlinear systems of fractional differential equations.Example 1.
Consider the following system of fractional differential equations [44]:(57)D3/2x1t=x2t,D1/2x2t=−x2t−x1t+t+1,
subject to(58)x10=x1′0=1,x20=0.
The exact solutions arex1t=t+1, x2t=0. By applying the presented method, we get(59)D3/2x1t=C1TFμt,D1/2x2t=C2TFμt,where Ci=ci0,ci1,…,cinT, i=1,2. Then(60)x1t=C1TP3/2+A1TFμt,x2t=C2TP1/2+A2TFμt.
By using equations (59) and (60), we rewrite equation (57) as follows:(61)C1T−C2TP1/2−A2TFμt=0,C2T+C2TP1/2+C1TP3/2+A2T+A1TFμt−t−1=0.
Collocating equations (61) at the nodes tk=k/n,k=0,1,…,n yields(62)C1T−C2TP1/2−A2TFμtk=0,C2T+C2TP1/2+C1TP3/2+A2T+A1TFμtk−tk−1=0.
Now, takingn=2,μ=1 and using equations (30) and (53), we obtain(63)A1=13,54,3512T.(64)A2=0,0,0T.(65)P1/2=433π3177931058521−2083154878098105−8158.(66)P3/2=8273π−1372355109333−3255321528811833−72558.
Using equations (63)–(66), we rewrite equations (62) as(67)10395πc10−5580c20+832c21−96c22=0,3465π2c10−3c11−c12+8937c20+1616c21+328c22=0,c10−c11+c12−43465π97c20+144c21+904c22=0,−8165c10+288c11+960c12+585231π+124c20−10816c21+1248c22=0,96656c10+44544c11+4736c12+97448−90090πc20+168064+135135πc21+34112+45045πc22=0,25928c10+30976c11+12608c12+5044+45045πc20+7488−45045πc21+47008+45045πc22=0.
By solving the linear system (67), we get(68)cij=0,i=1,2,j=0,1,2.
Hence, the solutions are(69)x1t=A1TF1t=13,54,35123−12t+10t24t−5t2t2=t+1,x2t=A2TF1t=0,0,03−12t+10t24t−5t2t2=0,which are the exact solutions. Also, for the same n=2 with μ=1/2 we can get the exact solutions as follows:(70)A1=13,1,83T.(71)A2=0,0,0T.(72)P1/2=13π2312335−108π6−132−50+18π3115−25103.(73)P3/2=135π9−3π−44+15π7003−75π−103+π15−5π−1753+25π23−3353.
Similarly, by substituting from equations (70)–(73) in equation (62) and solving the resulting system, we get(74)cij=0,i=1,2,j=0,1,2.Therefore,(75)x1t=A1TF1/2t=13,1,833−12t+10t4t−5tt=t+1,x2t=A2TF1/2t=0,0,03−12t+10t4t−5tt=0.Example 2.
Consider the following nonlinear system of fractional differential equations [44]:(76)D3/2x1t=x22t,D1/2x2t=−x12t−t7/2x2t+1+4πt4+2,
subject to(77)x10=x1′0=0,x20=0.
The exact solutions to this system arex1t=t2, x2t=4t/π. Applying the proposed method yields(78)D3/2x1t=C1TFμt,D1/2x2t=C2TFμt,x1t=C1TP3/2+A1TFμt,x2t=C2TP1/2+A2TFμt.
Equations (76) can be written after using equations (78) in the following form(79)C1T−C2TP1/2−A2TFμt=0,C2TFμt+C1TP3/2+A1TFμt2+t7/2C2TP1/2+A2TFμt−1+4πt4−2=0.
By collocating equations (79) at the nodes tk=k/n,k=0,1,…,n, we get(80)C1T−C2TP1/2−A2TFμtk=0,C2TFμtk+C1TP3/2+A1TFμtk2+t7/2C2TP1/2+A2TFμtk−1+4πtk4−2=0.
Withn=4 and μ=1/2 and solving the resulting system from equations (80) after calculating A1,A2,P1/2 and P3/2 from equations (30) and (53), we get(81)c10=−9.66486918194092×10−16,c11=0.1128379167095563,c12=0.5641895835477753,c13=1.57973083393373,c14=3.385137501286536,c20=0.4,c21=1.200000000000042,c22=2.000000000000043,c23=2.8,c24=3.59999999999999.
Therefore, the solutions are as follows:(82)x1t=0.999999999999985t2+1.903338597841753×10−14t3/2+1.821459649775647×10−15t−1.266348137463069×10−15t+5.637851296924623×10−17,≃t2,x2t=8.17124146124115×10−14t2+5.329070518200751×10−14t3/2−2.700062395888381×10−13t+2.256758334191171t−4.830120678422922×10−15,≃4πt,which can be seen as exact solutions. Table 1 describes the absolute errors of x1t and x2t. The symbols e1 and e2 that are used in Table 1 denote the absolute errors for x1t and x2t, respectively. Also, Table 1 contains a comparison of the computational results obtained by the introduced method and those in [44] in terms of the absolute errors. The results for the compared method have been executed by the authors in [44] and we use these results to make a direct comparison with the presented method. The symbol “–” means that the result is unavailable for that method. From Table 1, we can see that the absolute errors achieved by the presented method are close to zero and less than those in [44] for almost all values of t. This means that the presented method is more coincidence with exact solutions than [44] for this system and can find solutions with high accuracy within acceptable computational costs.Table 1
Comparison of the results for Example2.
Our methodAtabakzadeh’s [44]n=4, μ=1/2n=4te1e2e1e20.05.64×10−174.83×10−15—0.00.29.71×10−161.42×10−14—8.8×10−40.42.44×10−156.00×10−15—7.1×10−30.63.72×10−150.00—3.5×10−30.84.66×10−150.00—2.8×10−21.05.11×10−155.77×10−15—5.4×10−2Example 3.
Consider the fractional system of the pollution problem of three lakes interconnected by channels flowing between them that has been modeled in [45] and has been solved numerically in [46,47]:(83)Dα1x1t=381180x3t−202900x1t−182900x1t+sint+1,Dα2x2t=182900x1t−18850x2t,Dα3x3t=202900x1t+18850x2t−381180x3t.
subject to(84)x10=0,x20=0,x30=0,where 0<αi≤1,i=1,2,3 and xit denote the amount of pollution in lake i at time t. By applying the method described in Section 5, we obtain(85)Dα1x1t=C1TFμt,Dα2x2t=C2TFμt,Dα3x3t=C3TFμt,x1t=C1TPα1+A1TFμt,x2t=C2TPα2+A2TFμt,x3t=C3TPα3+A3TFμt.
By using equation (85), equation (83) take the form(86)C1T−381180C3TPα3+A3T+382900C1TPα1+A1TFμt−sint−1=0,C2T−182900C1TPα1+A1T+18850C2TPα2+A2TFμt=0,C3T−202900C1TPα1+A1T−18850C2TPα2+A2T+381180C3TPα3+A3TFμt=0.
Considerαi=1,i=1,2,3, n=3, and μ=1. After collocating equations (86) at the nodes tk=k/n,k=0,1,…,n and solving the resulting system, we obtain(87)x1t=−7.985745491227481×10−2t3+0.5413124555867391t2+0.989767107680799t+5.113128579998585×10−4,x2t=8.40375201373044×10−4t3+3.237087889257668×10−3t2−3.354110606453802×10−5t+1.833796936640798×10−6,x3t=9.46943417407036×10−4t3+3.59433215804088×10−3t2−3.674536458210798×10−5t+2.012284555870566×10−6.Similarly, the solutions at αi=0.75,i=1,2,3, n=3 and μ=0.75 take the form(88)x1t=3.21736742031109×10−2t9/4+0.5829745706710028t3/2+1.023986404601416t3/4+3.20165796818566×10−3,x2t=1.973507604532053×10−3t9/4+4.30731995528251×10−3t3/2+4.081331088971815×10−5t3/4−9.58440146269299×10−7,x3t=2.227008738973193×10−3t9/4+4.777099386813366×10−3t3/2+4.734654981000435×10−5t3/4−1.165144594284984×10−6.
The approximate solutions to System (83) for various values of αi, i=1,2,3 are displayed in Figures 1–3. The accuracy of the approximate solutions (87) and (88) can be checked by substituting them in equations (83), thus the resulting equations must be satisfied approximately, that is,(89)E1t=Dα1x1t−381180x3t+202900x1t+182900x1t−sint−1≃0,E2t=Dα2x2t−182900x1t+18850x2t≃0,E3t=Dα3x3t−202900x1t−18850x2t+381180x3t≃0,where t∈0,1. Tables 2 and 3 show some numerical values of the residual functions Eit obtained by equations (89) for different values of n,μ, and αi where i=1,2,3. By comparing the residual functions in these tables, it can be seen that the best results are corresponding to n=6,μ=1. This means that when n increases, the errors are reduced, and the approximate solutions approach to the exact solutions. Nonetheless, larger n will need more computational processes and this may produce a computational error. Table 4 describes a comparison of the results obtained by the proposed method and those in [47] in terms of the residual functions. It can be seen that the residual functions obtained by the presented method are less than those in [47] on 15 out of all 18 values of t. This means that the presented method is more accurate than [47] for this system.Figure 1
Approximate solutionx1t to Example 3 for n=3 and various values of αi=μ, i=1,2,3.Figure 2
Approximate solutionx2t to Example 3 for n=3 and various values of αi=μ, i=1,2,3.Figure 3
Approximate solutionx3t to Example 3 for n=3 and various values of αi=μ, i=1,2,3.Table 2
Residual functionsEit of Example 3 for different values of αi,i=1,2,3.
n=3,μ=0.75n=3,μ=0.9αi=0.75,i=1,2,3αi=0.9,i=1,2,3tE1tE2tE3tE1tE2tE3t0.05.89×10−21.76×10−52.14×10−52.67×10−22.83×10−53.06×10−50.24.87×10−39.05×10−71.1×10−64.01×10−31.93×10−62.09×10−60.41.43×10−41.33×10−71.61×10−71.58×10−32.25×10−62.43×10−60.67.23×10−32.12×10−62.58×10−63.85×10−33.67×10−63.97×10−60.81.86×10−36.23×10−77.57×10−76.89×10−41.5×10−61.62×10−61.02.58×10−27.74×10−69.41×10−61.93×10−22.05×10−52.22×10−5Table 3
Residual functionsEit of Example 3 when αi=1,i=1,2,3.
n=3, μ=1n=6, μ=1tE1tE2tE3tE1tE2tE3t0.01.02×10−23.67×10−54.02×10−51.58×10−61.71×10−81.91×10−80.29.11×10−42.64×10−62.9×10−62.44×10−81.44×10−101.6×10−100.41.31×10−34.99×10−65.47×10−61.46×10−99.75×10−111.09×10−100.61.48×10−34.99×10−65.47×10−61.64×10−89.75×10−111.09×10−100.85.18×10−42.64×10−62.9×10−62.54×10−91.44×10−101.6×10−101.01.02×10−23.67×10−54.02×10−51.58×10−61.71×10−81.91×10−8Table 4
Comparison of the results for Example3 when αi=1,i=1,2,3.
Our methodYüzbaşı’s [47]n=6, μ=1n=6tE1tE2tE3tE1tE2tE3t0.01.58×10−61.71×10−81.91×10−80.00.00.00.22.44×10−81.44×10−101.6×10−104.23×10−86.08×10−106.78×10−100.41.46×10−99.75×10−111.09×10−104.12×10−85.54×10−106.19×10−100.61.64×10−89.75×10−111.09×10−106.57×10−88.32×10−109.28×10−100.82.54×10−91.44×10−101.6×10−102.03×10−72.43×10−92.71×10−91.01.58×10−61.71×10−81.91×10−81.05×10−51.19×10−71.33×10−7Example 4.
Consider the fractional system of the space distribution of vibrations [46,48,49]:(90)Dα1x1t=x2t−x13t,Dα2x2t=x1t−0.54x15t,subject to(91)x10=1.2,x1′0=0,x20=1.6412,x2′0=0,where 0<αi≤2,i=1,2. Equations (90) represent a system of two coupled fractional differential equations with a fifth-order strong non-linearity. We can find the approximate solutions by applying the method described in Section 5. We have eliminated the details and have displayed the approximate solutions to System (90) for various values of αi, i=1,2 in Figures 4 and 5. The residual functions Eit,i=1,2 of system (90) take the following form:(92)E1t=Dα1x1t−x2t+x13t≃0,E2t=Dα2x2t−x1t+0.54x15t≃0.
Some numerical values of the residual functionsEit obtained by equations (92) for different values of n,μ, and αi, i=1,2 are shown in Tables 5–7. It can be seen that the presented method can find solutions with good accuracy within acceptable computational costs and the accuracy of the approximate solutions can be improved with the increment of n. Table 8 contains a comparison of the results obtained by the introduced method and those in [46] in terms of the residual functions. We have computed the residual functions for [46] from the given solutions x1 and x2 in that paper. From Table 8, we can see that the errors achieved by the presented method are less than those in [46] for all values of t. Therefore, we can say that the proposed method is more efficient and superior to [46] for this system.Figure 4
Approximate solutionx1t to Example 4 for n=3 and various values of αi=μ, i=1,2.Figure 5
Approximate solutionx2t to Example 4 for n=3 and various values of αi=μ, i=1,2.Table 5
Residual functionsEit of Example 4 at αi=0.5,i=1,2.
n=3, μ=0.5n=6, μ=0.5tE1tE2tE1tE2t0.01.7×10−22.92×10−29.87×10−41.79×10−30.22.33×10−52.69×10−52.62×10−54.82×10−50.41.26×10−32.17×10−31.05×10−51.89×10−50.61.71×10−32.95×10−33.42×10−56.21×10−50.83.91×10−46.62×10−41.39×10−52.52×10−51.02.91×10−35.0×10−31.3×10−42.37×10−4Table 6
Residual functionsEit of Example 4 at αi=1,i=1,2.
n=3, μ=1n=6, μ=1tE1tE2tE1tE2t0.05.81×10−31.04×10−25.28×10−51.21×10−40.25.44×10−41.07×10−31.72×10−76.62×10−70.47.43×10−41.29×10−36.65×10−71.67×10−60.68.3×10−41.51×10−31.85×10−81.7×10−70.83.37×10−45.39×10−48.57×10−72.12×10−61.05.81×10−31.04×10−25.28×10−51.21×10−4Table 7
Residual functionsEit of Example 4 at αi=2,i=1,2.
n=3, μ=1n=6, μ=1tE1tE2tE1tE2t0.01.9×10−23.17×10−21.09×10−42.39×10−40.23.89×10−36.33×10−38.46×10−61.85×10−50.42.14×10−33.56×10−35.08×10−61.1×10−50.61.44×10−32.25×10−35.1×10−61.1×10−50.83.67×10−35.93×10−38.63×10−61.85×10−51.01.1×10−21.7×10−21.14×10−42.42×10−4Table 8
Comparison of the results for Example4 at αi=1,i=1,2.
Our methodKhader’s [46]n=3, μ=1n=3tE1tE2tE1tE2t0.05.81×10−31.04×10−21.730.140.25.44×10−41.07×10−30.890.650.47.43×10−41.29×10−31.221.030.68.3×10−41.51×10−33.742.250.83.37×10−45.39×10−412.4217.261.05.81×10−31.04×10−249.34261.96Example 5.
Consider the incommensurate system of fractional differential equations [37]:(93)D6/5x1t=x2t+t+et,D3/2x2t=5x1t+x2t+e2t,subject to(94)x10=1,x1′0=Γ1.8,x20=1,x2′0=3.
The approximate solutions can be obtained by applying the presented method in Section5. We have eliminated the details and displayed the approximate solutions to System (93) for various values of n in Figures 6 and 7. The residual functions Eit,i=1,2 of System (93) can be written as follows:(95)E1t=D6/5x1t−x2t−t−et≃0,E2t=D3/3x2t−5x1t−x2t−e2t≃0.
Numerical values of the residual functionsEit,i=1,2 obtained by equations (95) for different values of n are presented in Table 9. System (93) has been solved in [37]. They presented solutions in terms of three-parameter Mittag-Leffler functions and some infinite series. From Table 9, it can be seen that the presented method can find solutions that have acceptable accuracy with respect to the computational costs and the accuracy of the approximate solutions can be improved with the increment of n. By comparing our solutions with those in [37], we can say that the presented method gives solutions with simple structure and can be computed easily compared with those in [37].Figure 6
Approximate solutionx1t to Example 5 for various values of n.Figure 7
Approximate solutionx2t to Example 5 for various values of n.Table 9
Residual functionsEit,i=1,2 of Example 5.
n=6, μ=1n=9, μ=1tE1tE2tE1tE2t0.21.101.010.970.770.40.970.780.840.520.60.880.590.780.430.80.840.560.730.391.00.760.220.730.55


## 7. Conclusions

In this paper, an efficient method has been presented to handle systems of fractional differential equations. The proposed method can be applied without using linearization or perturbation unlike some other methods such as the homotopy perturbation method. It reduces the system of fractional differential equations into a system of algebraic equations and can be adapted to implement on computers. The reliability and accuracy of the presented method are tested on five different linear and nonlinear systems. From the results, we can see that the obtained solutions are either in good agreement with the exact solutions if the exact solutions are known or satisfy the system equations if the exact solutions are unknown. A comparison with other methods reveals that our method is more accurate and needs a small number of FCHFs to achieve good results. Finally, we can observe that using basis functions of fractional order gives more accurate results than integer-order basis functions, which demonstrates the significance of this method for solving fractional systems.



---

*Source: 4862650-2022-04-16.xml*