Processing math: 39%
Research article

An automatic density peaks clustering based on a density-distance clustering index

  • Received: 06 August 2023 Revised: 07 October 2023 Accepted: 23 October 2023 Published: 25 October 2023
  • MSC : 91C20

  • The density peaks clustering (DPC) algorithm plays an important role in data mining by quickly identifying cluster centers using decision graphs to identify arbitrary clusters. However, the decision graph introduces uncertainty in determining the cluster centers, which can result in an incorrect number of clusters. In addition, the cut-off distance parameter relies on prior knowledge, which poses a limitation. To address these issues, we propose an improved automatic density peaks clustering (ADPC) algorithm. First, a novel clustering validity index called density-distance clustering (DDC) is introduced. The DDC index draws inspiration from the density and distance characteristics of cluster centers, which is applicable to DPC and aligns with the general definition of clustering. Based on the DDC index, the ADPC algorithm automatically selects the suitable cut-off distance and acquires the optimal number of clusters without additional parameters. Numerical experimental results validate that the introduced ADPC algorithm successfully automatically determines the optimal number of clusters and cut-off distance, significantly outperforming DPC, AP and DBSCAN algorithms.

    Citation: Xiao Xu, Hong Liao, Xu Yang. An automatic density peaks clustering based on a density-distance clustering index[J]. AIMS Mathematics, 2023, 8(12): 28926-28950. doi: 10.3934/math.20231482

    Related Papers:

    [1] Yanping Yang, Muhammad Shoaib Saleem, Waqas Nazeer, Ahsan Fareed Shah . New Hermite-Hadamard inequalities in fuzzy-interval fractional calculus via exponentially convex fuzzy interval-valued function. AIMS Mathematics, 2021, 6(11): 12260-12278. doi: 10.3934/math.2021710
    [2] Jamshed Nasir, Saber Mansour, Shahid Qaisar, Hassen Aydi . Some variants on Mercer's Hermite-Hadamard like inclusions of interval-valued functions for strong Kernel. AIMS Mathematics, 2023, 8(5): 10001-10020. doi: 10.3934/math.2023506
    [3] Hari Mohan Srivastava, Soubhagya Kumar Sahoo, Pshtiwan Othman Mohammed, Bibhakar Kodamasingh, Kamsing Nonlaopon, Khadijah M. Abualnaja . Interval valued Hadamard-Fejér and Pachpatte Type inequalities pertaining to a new fractional integral operator with exponential kernel. AIMS Mathematics, 2022, 7(8): 15041-15063. doi: 10.3934/math.2022824
    [4] Miguel Vivas-Cortez, Muhammad Aamir Ali, Artion Kashuri, Hüseyin Budak . Generalizations of fractional Hermite-Hadamard-Mercer like inequalities for convex functions. AIMS Mathematics, 2021, 6(9): 9397-9421. doi: 10.3934/math.2021546
    [5] Manar A. Alqudah, Artion Kashuri, Pshtiwan Othman Mohammed, Muhammad Raees, Thabet Abdeljawad, Matloob Anwar, Y. S. Hamed . On modified convex interval valued functions and related inclusions via the interval valued generalized fractional integrals in extended interval space. AIMS Mathematics, 2021, 6(5): 4638-4663. doi: 10.3934/math.2021273
    [6] Iqra Nayab, Shahid Mubeen, Rana Safdar Ali, Faisal Zahoor, Muath Awadalla, Abd Elmotaleb A. M. A. Elamin . Novel fractional inequalities measured by Prabhakar fuzzy fractional operators pertaining to fuzzy convexities and preinvexities. AIMS Mathematics, 2024, 9(7): 17696-17715. doi: 10.3934/math.2024860
    [7] Zehao Sha, Guoju Ye, Dafang Zhao, Wei Liu . On interval-valued K-Riemann integral and Hermite-Hadamard type inequalities. AIMS Mathematics, 2021, 6(2): 1276-1295. doi: 10.3934/math.2021079
    [8] Thabet Abdeljawad, Muhammad Aamir Ali, Pshtiwan Othman Mohammed, Artion Kashuri . On inequalities of Hermite-Hadamard-Mercer type involving Riemann-Liouville fractional integrals. AIMS Mathematics, 2021, 6(1): 712-725. doi: 10.3934/math.2021043
    [9] Muhammad Bilal Khan, Muhammad Aslam Noor, Thabet Abdeljawad, Bahaaeldin Abdalla, Ali Althobaiti . Some fuzzy-interval integral inequalities for harmonically convex fuzzy-interval-valued functions. AIMS Mathematics, 2022, 7(1): 349-370. doi: 10.3934/math.2022024
    [10] Muhammad Bilal Khan, Pshtiwan Othman Mohammed, Muhammad Aslam Noor, Abdullah M. Alsharif, Khalida Inayat Noor . New fuzzy-interval inequalities in fuzzy-interval fractional calculus by means of fuzzy order relation. AIMS Mathematics, 2021, 6(10): 10964-10988. doi: 10.3934/math.2021637
  • The density peaks clustering (DPC) algorithm plays an important role in data mining by quickly identifying cluster centers using decision graphs to identify arbitrary clusters. However, the decision graph introduces uncertainty in determining the cluster centers, which can result in an incorrect number of clusters. In addition, the cut-off distance parameter relies on prior knowledge, which poses a limitation. To address these issues, we propose an improved automatic density peaks clustering (ADPC) algorithm. First, a novel clustering validity index called density-distance clustering (DDC) is introduced. The DDC index draws inspiration from the density and distance characteristics of cluster centers, which is applicable to DPC and aligns with the general definition of clustering. Based on the DDC index, the ADPC algorithm automatically selects the suitable cut-off distance and acquires the optimal number of clusters without additional parameters. Numerical experimental results validate that the introduced ADPC algorithm successfully automatically determines the optimal number of clusters and cut-off distance, significantly outperforming DPC, AP and DBSCAN algorithms.



    System of fractional differential equations with incommensurate order derivatives have received increasing attention recently as this incommensurate order derivative is better in describing the real phenomena, such as financial system [2,3], circuit simulation [4], eco-epidemiological model[5], HIV model [6] and modeling glucose-insulin regulatory system [7]. In this research direction, many works had been done to study stability analysis [8,9,10,11], synchronization [12] and other rich dynamical behaviour [13,14].

    Due to the emerging of cross-discipline research in this incommensurate fractional order system, finding the solution of the incommensurate fractional order system is becoming more and more important. In this case, numerical methods, such as the predictor-corrector scheme [15,16], are always used to obtain the solution for the incommensurate fractional order system. Apart from this, some algorithms are developed to obtain the approximation solution for incommensurate fractional order systems, such as the Adomian decomposition algorithm [17], reduced-order model approximation via genetic algorithm [18]. However, not much research was done to find the analytical solution or exact solution for this incommensurate fractional order system. Until recently, Huseynov et al. in [1] successfully derive the analytical solution for the incommensurate fractional order 0<α,β<1 by converting the system into a corresponding Volterra integral equation. Besides that, Ahmadova et al. [19] found the analytical solution for this incommensurate fractional order system via trivariate Mittag-Leffer functions. However, their proposed methods are only limited to incommensurate fractional order 0<α,β<1. Hence, this motivates us to derive the analytical solution for a higher order of incommensurate fractional order system.

    In this paper, we extend the work by Huseynov et al. in [1], which is limited for incommensurate fractional order system for 0<α,β<1. We intend to derive the analytical solution of higher order incommensurate fractional order system. Specifically, for α,β(1,2), we consider the incommensurate fractional order system as follows:

    CDαx1(t)=a11x1(t)+a12x2(t)+g1(t),CDβx2(t)=a21x1(t)+a22x2(t)+g2(t), (1.1)

    with initial condition x1(0)=x01, x2(0)=x02, x1(0)=x11 and x2(0)=x12. The physical meaning of such an incommensurate fractional order system as well as the advantages of using incommensurate models over the classical one (compare to commensurate models) are shown in [2–14]. The fractional derivatives are defined with Caputo sense and the initial value problems to be solved for x1,x2C1[0,). Similar to the works in [1], we convert the system in (1.1) to Volterra integral equations and Picard's successive approximations were used to derive the analytical expression of the solution for an incommensurate fractional order system for 1<α,β<2. Similar to [1], we use Picard's successive approximations to solve the Volterra integral equations arise because this method is based on the Banach fixed point theorem. In order to obtain the fixed point of a functional operator, start with an arbitrary function (i.e. the zeroth approximation) and apply the operator repeatedly to obtain a sequence of successive approximations which should converge towards the fixed point. This method has been applied to derive the explicit analytical solution of incommensurate fractional differential equation systems with fractional order 0<α,β<1 [1]. The solution will be simplified via some combinatorial concepts and bivariate Mittag-Leffler function. In short, this paper aims to contribute to analytical method that gives new explicit solutions to a certain class of fractional differential systems. These kind of explicit solutions for solving fractional differential equations or systems have been increasingly investigated by researchers in these research areas, such as in [20,21,22,23,24]. In short, we hope to contribute in obtaining an explicit analytical solution for fractional calculus problems, which is relatively less investigated compared to numerical solution, such as in [25,26,27,28,29,30,31,32].

    The rest of this paper is structured as follows. Section 2 is devoted to some preliminaries regarding some important definitions, concepts and notations in fractional calculus and special functions. Section 3 is devoted to presenting the derivation of analytical solutions for the incommensurate fractional order system in higher order. Moreover, some special cases will be discussed in Section 4. Sections 5 and 6 are devoted to presenting some examples and conclusion of this paper, respectively.

    In this section, we briefly explain some important definitions, concepts and notations in fractional calculus and special functions, which is important for obtaining the analytical solution for this incommensurate fractional order system.

    Definition 1. Let α>0, n=[α]+1 if αN, n=α if αN and x>0. The left Caputo fractional derivative of a function of order α, denoted by CDαxf(x) is

    CDαx f(x)=1Γ(nα)x0f(n)(τ)(xτ)αn+1dτ, (2.1)

    with n1α<n.

    For Caputo fractional derivative, we have this important expression:

    CDαx xβ=Γ(β+1)Γ(β+1α)xβα,forβ>α. (2.2)

    Definition 2. For Re(α),Re(β)>0, the classical Mittag-Leffler function (i.e. one parameter) and two-parameter Mittag-Leffler function are defined as

    Eα(t)=k=0tkΓ(αk+1),Eα,β(t)=k=0tkΓ(αk+β). (2.3)

    Definition 3. [33] For Re(α),Re(β),Re(γ)>0, the three-parameter version of bivariate Mittag-Leffler function can be defined as:

    Eα,β,γ(x,y)=k=0l=0(k+l)!xkylΓ(αk+βl+γ)k!l!. (2.4)

    The convergence of this bivariate Mittag-Leffler function was shown in Section 2, the new bivariate Mittag-Leffler function in [33]. The Mittag-Leffler function is used as the solution of system of fractional differential equations as this Mittag-Leffler function is the generalization of the exponential function, which exponential function is widely used to express the solution of integer order system of differential equations. The Mittag-Leffler function is a series which the terms are up to infinity. Hence, to calculate these Mittag-Leffler functions, ones can refer the numerical algorithm such as in [34,35,36]. With Caputo fractional derivative, we have this important expression for the fractional derivative involving Mittag-Leffler function:

    CDαx(Eα(λxα))=λ(Eα(λxα)). (2.5)

    In this paper, we will use some important integration with was introduced in [1] as follows:

    t0(tτ)a1τb1dτ=Γ(a)Γ(b)Γ(a+b)ta+b1,fora>0,b>0. (2.6)
    tu(tτ)a1(τu)b1dτ=Γ(a)Γ(b)Γ(a+b)(tu)a+b1,fora>0,b>0. (2.7)
    t0τ0(tτ)a1(τu)b1f(u)dudτ=Γ(a)Γ(b)Γ(a+b)t0(tu)a+b1f(u)du,fora>0,b>0. (2.8)

    Remarks: We can also write Γ(a)Γ(b)Γ(a+b)=B(a,b), where B(a,b) is the Beta function.

    For the f(τ)=τv, where v>0, using Eq (2.6), we have the following integration involving Mittag-Leffler function.

    t0(tτ)aEα,β(λ(tτ)b)τvdτ=k=0t0λk(tτ)bk+aτvdτΓ(αk+β)=k=0λkΓ(bk+a+1)Γ(v+1)tbk+a+v+1Γ(αk+β)Γ(bk+a+v+2). (2.9)

    If a=β1 and b=α, from Eq (2.9), we obtain

    t0(tτ)β1τvEα,β(λ(tτ)α)dτ=Γ(v+1)tβ+vEα,β+v+1(λtα).

    The lower incomplete gamma function is defined for Re(α)>0,Re(z)>0 as follows:

    γ(α,z)=z0ettα1dt. (2.10)

    Definition 4. Hypergeometric functions 2F1(a1,a2;b;z) and 1F2(a;b1,b2;z) are defined by the series

    2F1(a1,a2;b;z)=k=0(a1)k(a2)k(b)kzkk!,|z|<1,1F2(a;b1,b2;z)=k=0(a)k(b1)k(b2)kzkk!, (2.11)

    where the pochhammer symbol, (a)k=Γ(a+k)Γ(a).

    For the sake of simplicity, throughout the writing, we use n1,n2,,nk=0 to represent multiple series n1=0n2=0nk=0.

    In order to derive the analytical solutions for the incommensurate fractional differential equation systems with order 1<α,β<2 as in Eq (1.1), basically one can follow the following steps:

    Step 1: Write the system in Volterra integral equations of second kind.

    Step 2: Perform the Picard's successive approximations.

    Step 3: Simplify the solution by using some combinatorial formulae.

    Step 4: Verify the solution by using substitution.

    Here, we will derive the inhomogeneous case. By setting g1(t)=0 and g2(t), the Eq (1.1) will reduce to the homogeneous case. Similarly, if we take the value of α=β, the incommensurate fractional differential equation systems with fractional order 1<α,β<2 will be reduced to commensurate fractional differential equation systems with fractional order 1 to 2.

    Step 1: Write the system in Volterra integral equations of second kind.

    Using the result from Theorem 5.15 in [37], we obtain the single fractional differential equation for σ(1,2) in Caputo sense as follows:

    CDσy(t)=λy(t)+h(t),y(0)=y0,y(0)=y1. (3.1)

    We have the following solution:

    y(t)=y0Eσ(λtσ)+y1tEσ,2(λtσ)+t0(tτ)σ1h(τ)Eσ,σ(λ(tτ)σ)dτ. (3.2)

    Using Eq (3.2), the Volterra integral equation of second kind for the equation in (1.1) can be written as

    x1(t)=x01Eα(a11tα)+x11tEα,2(a11tα)+t0(tτ)α1[a12x2(τ)+g1(τ)]Eα,α(a11(tτ)α)dτ,x2(t)=x02Eβ(a22tβ)+x12tEβ,2(a22tβ)+t0(tτ)β1[a21x1(τ)+g2(τ)]Eβ,β(a22(tτ)β)dτ. (3.3)

    Substituting x2(t) into the first equation in (3.3) and x1(t) into the second equation in (3.3), we obtain the following:

    (3.4)

    Using identity as in Eq (2.6) to Eq (2.8), we obtain

    x1(t)=x01Eα(a11tα)+x11tEα,2(a11tα)+a12x02n1,n2=0an111an222tn1α+n2β+αΓ(n1α+n2β+α+1)+a12x12n1,n2=0an111an222tn1α+n2β+α+1Γ(n1α+n2β+α+2)+n1=0an111Γ(n1α+α)t0(tτ)n1α+α1g1(τ)dτ+a12n1,n2=0an111an222Γ(n1α+n2β+α+β)t0(tτ)n1α+n2β+α+β1(a21x1(τ)+g2(τ))dτ. (3.5)

    By using a similar approach, we obtain the expression for x2(t) as follows:

    x2(t)=x02Eβ(a22tβ)+x12tEβ,2(a22tβ)+a21x01n1,n2=0an122an211tn1β+n2α+βΓ(n1β+n2α+β+1)+a21x11n1,n2=0an122an211tn1β+n2α+β+1Γ(n1β+n2α+β+2)+a21n1,n2=0an122an211Γ(n1β+n2α+α+β)t0(tτ)n1β+n2α+α+β1(a12x2(τ)+g1(τ))dτ+n1=0an122Γ(n1β+β)t0(tτ)n1β+β1g2(τ)dτ. (3.6)

    Step 2: Perform the Picard's successive approximation.

    Using Picard's successive approximation, the solution of the Volterra integral equations as in Eqs (3.5) and (3.6) can be obtained via setting

    (3.7)

    and

    (3.8)

    For m=1, using the identities (2.6)–(2.8), we have

    x1,1(t)=x1,0(t)+a12a21x01n1,n2,n3=0an1+n311an222t(n1+n3)α+n2β+α+βΓ((n1+n3)α+n2β+α+β+1)+a12a21x11n1,n2,n3=0an1+n311an222t(n1+n3)α+n2β+α+β+1Γ((n1+n3)α+n2β+α+β+2)+a212a21x02n1,n2,n3,n4=0an1+n311an2+n422t(n1+n3)α+(n2+n4)β+2α+βΓ((n1+n3)α+(n2+n4)β+2α+β+1)+a212a21x12n1,n2,n3,n4=0an1+n311an2+n422t(n1+n3)α+(n2+n4)β+2α+β+1Γ((n1+n3)α+(n2+n4)β+2α+β+2)+a12a21n1,n2,n3=0an1+n311an222t0(tτ)(n1+n3)α+n2β+2α+β1g1(τ)dτΓ((n1+n3)α+n2β+2α+β)+a212a21n1,n2,n3,n4=0an1+n311an2+n422t0(tτ)(n1+n3)α+(n2+n4)β+2α+2β1g2(τ)dτΓ((n1+n3)α+(n2+n4)β+2α+2β). (3.9)

    Meanwhile, for m=2, we have

    (3.10)

    In general, after some algebraic manipulation, we obtain

    (3.11)

    where n1++n2m+1 and n2++n2m denote n1+n3++n2m+1 and n2+n4++n2m, respectively.

    When m, we can rewrite the solution of x1(t) as follows:

    (3.12)

    Similarly, by symmetry, we have successive approximations for x2(t)=limmx2,m(t) as follows:

    (3.13)

    Step 3: Simplify the solution by using some combinatorial formulae.

    We write j as all the odd-indexed terms together and m as all the even-indexed appear together, i.e. j=n1+n3++n2k+1, m=n2+n4++n2k or m=n2+n4++n2k+2, we obtain,

    x1(t)=x01n1=0an111tn1αΓ(n1α+1)+x01k=1j=0m=0n1,n2,,n2k+1:n1+n3++n2k+1=j,n2+n4++n2k=mak12ak21aj11am22t(j+k)α+(m+k)βΓ((j+k)α+(m+k)β+1)+x11n1=0an111tn1α+1Γ(n1α+2)+x11k=1j=0m=0n1,n2,,n2k+1:n1+n3++n2k+1=j,n2+n4++n2k=mak12ak21aj11am22t(j+k)α+(m+k)β+1Γ((j+k)α+(m+k)β+2)+x02k=0j=0m=0n1,n2,,n2k+2:n1+n3++n2k+1=j,n2+n4++n2k+2=mak+112ak21aj11am22t(j+k+1)α+(m+k)βΓ((j+k+1)α+(m+k)β+1)+x12k=0j=0m=0n1,n2,,n2k+2:n1+n3++n2k+1=j,n2+n4++n2k+2=mak+112ak21aj11am22t(j+k+1)α+(m+k)β+1Γ((j+k+1)α+(m+k)β+2)+n1=0an111t0(tτ)n1α+α1g1(τ)dτΓ(n1α+α)+k=1j=0m=0n1,n2,,n2k+1:n1+n3++n2k+1=j,n2+n4++n2k=mak12ak21aj11am22t0(tτ)(j+k+1)α+(m+k)β1g1(τ)dτΓ((j+k+1)α+(m+k)β)+k=0j=0m=0n1,n2,,n2k+2:n1+n3++n2k+1=j,n2+n4++n2k+2=mak+112ak21aj11am22t0(tτ)(j+k+1)α+(m+k+1)β1g2(τ)dτΓ((j+k+1)α+(m+k+1)β). (3.14)

    Then, we have the simple combinatorial identity as follows for any k, j and m:

    n1,n2,,n2k+1:n1+n3++n2k+1=j,n2+n4++n2k=m(1)=|{(n1,n3,,n2k+1):=j}||{n2,n4,,n2k):=m}|=(k+j)!k!j!(k+m1)!(k1)!m!=(k+jk)(k+m1k1) (3.15)
    n1,n2,,n2k+2:n1+n3++n2k+1=j,n2+n4++n2k+2=m(1)=|{(n1,n3,,n2k+1):=j}||{n2,n4,,n2k+2):=m}|=(k+j)!k!j!(k+m)!k!m!=(k+jk)(k+mk) (3.16)

    Applying Eqs (3.15) and (3.16) to Eq (3.14) yields

    x1(t)=x01n1=0an111tn1αΓ(n1α+1)+x01k=1j=0m=0ak12ak21aj11am22t(j+k)α+(m+k)β(k+jk)(k+m1k1)Γ((j+k)α+(m+k)β+1)+x11n1=0an111tn1α+1Γ(n1α+2)+x11k=1j=0m=0ak12ak21aj11am22t(j+k)α+(m+k)β+1(k+jk)(k+m1k1)Γ((j+k)α+(m+k)β+2)+x02k=0j=0m=0ak+112ak21aj11am22t(j+k+1)α+(m+k)β(k+jk)(k+mk)Γ((j+k+1)α+(m+k)β+1)+x12k=0j=0m=0ak+112ak21aj11am22t(j+k+1)α+(m+k)β+1(k+jk)(k+mk)Γ((j+k+1)α+(m+k)β+2)+n1=0an111t0(tτ)n1α+α1g1(τ)dτΓ(n1α+α)+k=1j=0m=0ak12ak21aj11am22t0(tτ)(j+k+1)α+(m+k)β1g1(τ)dτ(k+jk)(k+m1k1)Γ((j+k+1)α+(m+k)β)+k=0j=0m=0ak+112ak21aj11am22t0(tτ)(j+k+1)α+(m+k+1)β1g2(τ)dτ(k+jk)(k+mk)Γ((j+k+1)α+(m+k+1)β). (3.17)

    By writing some of the terms in Mittag-Leffler function and let all the summations start from 0, we have

    x1(t)=x01Eα(a11tα)+x01k=0j=0m=0ak+112ak+121aj11am22t(j+k+1)α+(m+k+1)β(k+j+1k+1)(k+mk)Γ((j+k+1)α+(m+k+1)β+1)+x11tEα,2(a11tα)+x11k=0j=0m=0ak+112ak+121aj11am22t(j+k+1)α+(m+k+1)β+1(k+j+1k+1)(k+mk)Γ((j+k+1)α+(m+k+1)β+2)+x02k=0j=0m=0ak+112ak21aj11am22t(j+k+1)α+(m+k)β(k+jk)(k+mk)Γ((j+k+1)α+(m+k)β+1)+x12k=0j=0m=0ak+112ak21aj11am22t(j+k+1)α+(m+k)β+1(k+jk)(k+mk)Γ((j+k+1)α+(m+k)β+2)+t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ+k=0j=0m=0ak+112ak+121aj11am22t0(tτ)(j+k+2)α+(m+k+1)β1g1(τ)dτ(k+j+1k+1)(k+mk)Γ((j+k+2)α+(m+k+1)β)+k=0j=0m=0ak+112ak21aj11am22t0(tτ)(j+k+1)α+(m+k+1)β1g2(τ)dτ(k+jk)(k+mk)Γ((j+k+1)α+(m+k+1)β). (3.18)

    In the same manner, for x2(t), we obtain the following expression

    (3.19)

    Using (k+j+1k+1)(k+mk)=(k+j+1)!(k+1)!j!(k+m)!k!m! and (k+jk)(k+mk)=(k+j)!k!j!(k+m)!k!m! and assuming p=j+k and q=m+k and a11,a220, we obtain

    (3.20)

    The above equation can also be rewritten as follows:

    x1(t)=x01Eα(a11tα)+x11tEα,2(a11tα)+t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ+x01a12a21a22p=0q=1ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)min(p,q1)k=0(a12a21a11a22)k(p+1)!(k+1)!(pk)!(q1)!k!(qk1)!+x11a12a21a22p=0q=1ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)min(p,q1)k=0(a12a21a11a22)k(p+1)!(k+1)!(pk)!(q1)!k!(qk1)!+x02a12p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!+x12a12p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!+a12a21a22p=0q=1ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)×min(p,q1)k=0(a12a21a11a22)k(p+1)!(k+1)!(pk)!(q1)!k!(qk1)!+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)×min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!. (3.21)

    Similarly, by symmetric, we obtain

    x2(t)=x02Eβ(a22tβ)+x12tEβ,2(a22tβ)+t0(tτ)β1Eβ,β(a22(tτ)β)g2(τ)dτ+x02a12a21a11p=1q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)min(p1,q)k=0(a12a21a11a22)k(p1)!k!(pk1)!(q+1)!(k+1)!(qk)!+x12a12a21a11p=1q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)min(p1,q)k=0(a12a21a11a22)k(p1)!k!(pk1)!(q+1)!(k+1)!(qk)!+x01a21p=0q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!+x11a21p=0q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!+a21p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g1(τ)dτΓ(pα+qβ+α+β)min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!+a12a21a11p=1q=0ap11aq22t0(tτ)pα+qβ+2β1g2(τ)dτΓ(pα+qβ+2β)×min(p1,q)k=0(a12a21a11a22)k(p1)!k!(pk1)!(q+1)!(k+1)!(qk)!. (3.22)

    Letting A=a12a21a11a22 with a11,a220, we can simplify the inner series in (3.21) into hypergeometric function expression using the following identities.

    min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!=2F1(p,q;1;A),min(p,q1)k=0(a12a21a11a22)k(p+1)!(k+1)!(pk)!(q1)!k!(qk1)!=(p+1)2F1(p,1q;2;A). (3.23)

    Hence, we have x1(t) as follows:

    x1(t)=x01Eα(a11tα)+x11tEα,2(a11tα)+t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ+x01a11Ap=0q=1ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+1)2F1(p,1q;2;A)+x11a11Ap=0q=1ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+1)2F1(p,1q;2;A)+x02a12p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)2F1(p,q;1;A)+x12a12p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)2F1(p,q;1;A)+a11Ap=0q=1ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)(p+1)2F1(p,1q;2;A)+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)2F1(p,q;1;A). (3.24)

    The same is applied for x2(t), where we can simplify the inner series in (3.22) into a hypergeometric function expression using the following identities.

    min(p,q)k=0(a12a21a11a22)kp!k!(pk)!q!k!(qk)!=2F1(p,q;1;A),min(p1,q)k=0(a12a21a11a22)k(p1)!k!(pk1)!(q+1)!(k+1)!(qk)!=(q+1)2F1(1p,q;2;A). (3.25)

    Hence, we obtain x2(t) as follows:

    x2(t)=x02Eβ(a22tβ)+x12tEβ,2(a22tβ)+t0(tτ)β1Eβ,β(a22(tτ)β)g2(τ)dτ+x02a22Ap=1q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)(q+1)2F1(1p,q;2;A)+x12a22Ap=1q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)(q+1)2F1(1p,q;2;A)+x01a21p=0q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)2F1(p,q;1;A)+x11a21p=0q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)2F1(p,q;1;A)+a22Ap=1q=0ap11aq22t0(tτ)pα+qβ+2β1g2(τ)dτΓ(pα+qβ+2β)(q+1)2F1(1p,q;2;A)+a21p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g1(τ)dτΓ(pα+qβ+α+β)2F1(p,q;1;A). (3.26)

    Substituting q=0 in the double series with q=1 in Eq (3.24) (i.e. the 4th, 5th and 8th terms in the RHS of Eq (3.24)), and using 2F1(p,1;2;A)=1(1A)p+1A(p+1), we obtain the following new expression for the 4th, 5th and 8th terms in the RHS of Eq (3.24).

    x01a11Ap=0ap11t(p+1)αΓ((p+1)α+1)(p+1)2F1(p,1;2;A)=x01Eα(a11tα)x01Eα(a11(1A)tα),x11a11Ap=0ap11t(p+1)α+1Γ((p+1)α+2)(p+1)2F1(p,1;2;A)=x11tEα,2(a11tα)x11tEα,2(a11(1A)tα),a11Ap=0ap11t0(tτ)(p+1)α+α1g1(τ)dτΓ((p+1)α+α)(p+1)2F1(p,1;2;A)=t0(tτ)α1[Eα,α(a11(tτ)α)Eα,α(a11(1A)(tτ)α)]g1(τ)dτ. (3.27)

    Hence, we have the final solutions for the x1(t) to the system (1.1) as follows:

    x1(t)=x01Eα(a11(1A)tα)+x11tEα,2(a11(1A)tα)+t0(tτ)α1Eα,α(a11(1A)(tτ)α)g1(τ)dτ+x01a11Ap=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+1)2F1(p,1q;2;A)+x11a11Ap=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+1)2F1(p,1q;2;A)+x02a12p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)2F1(p,q;1;A)+x12a12p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)2F1(p,q;1;A)+a11Ap=0q=0ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)(p+1)2F1(p,1q;2;A)+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)2F1(p,q;1;A). (3.28)

    Using similar approach, the final solutions for the x2(t) to the system (1.1) is given as follows:

    x2(t)=x02Eβ(a22(1A)tβ)+x12tEβ,2(a22(1A)tβ)+t0(tτ)β1Eβ,β(a22(1A)(tτ)β)g2(τ)dτ+x02a22Ap=0q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)(q+1)2F1(1p,q;2;A)+x12a22Ap=0q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)(q+1)2F1(1p,q;2;A)+x01a21p=0q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)2F1(p,q;1;A)+x11a21p=0q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)2F1(p,q;1;A)+a21p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g1(τ)dτΓ(pα+qβ+α+β)2F1(p,q;1;A)+a22Ap=0q=0ap11aq22t0(tτ)pα+qβ+2β1g2(τ)dτΓ(pα+qβ+2β)(q+1)2F1(1p,q;2;A). (3.29)

    Step 4: Verify the solution by using substitution.

    Finally, we can verify the solutions by substituting (3.28) (i.e. x1(t)) and (3.29) (i.e. x2(t)) into CDαx1(t)=a11x1(t)+a12x2(t)+g1(t), which is the first equation of incommensurate fractional order system (1.1). Hence, the right-hand-side of the first equation of (1.1) is given by

    a11x1(t)+a12x2(t)+g1(t)=x01a11Eα(a11tα)+x11a11tEα,2(a11tα)+[a11t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ+g1(t)]+x02a12Eβ(a22tβ)+x12a12tEβ,2(a22tβ)+a12t0(tτ)β1Eβ,β(a22(tτ)β)g2(τ)dτ+x01a11Ap=0q=1ap11aq22tpα+qβΓ(pα+qβ+1)(p2F1(1p,1q;2;A)+2F1(p,1q;1;A))+x11a11Ap=0q=1ap11aq22tpα+qβ+1Γ(pα+qβ+2)(p2F1(1p,1q;2;A)+2F1(p,1q;1;A))+x02a12p=1q=0ap11aq22tpα+qβΓ(pα+qβ+1)(2F1(1p,q;1;A)+Aq2F1(1p,1q;2;A))+x12a12p=1q=0ap11aq22tpα+qβ+1Γ(pα+qβ+2)(2F1(1p,q;1;A)+Aq2F1(1p,1q;2;A))+a11Ap=0q=1ap11aq22t0(tτ)pα+qβ+α1g1(τ)dτΓ(pα+qβ+α)(p2F1(1p,1q;2;A)+2F1(p,1q;1;A))+a12p=1q=0ap11aq22t0(tτ)pα+qβ+β1g2(τ)dτΓ(pα+qβ+β)(2F1(1p,q;1;A)+Aq2F1(1p,1q;2;A)), (3.30)

    while the left-hand-side of the equation is given by

    CDαx1(t)=a11x01Eα(a11tα)+a11x11tEα,2(a11tα)+CDα(t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ)+x01a11Ap=0q=1ap11aq22tpα+qβΓ(pα+qβ+1)(p+1)2F1(p,1q;2;A)+x11a11Ap=0q=1ap11aq22tpα+qβ+1Γ(pα+qβ+2)(p+1)2F1(p,1q;2;A)+x02a12Eβ(a22tβ)+x02a12p=1q=0ap11aq22tpα+qβΓ(pα+qβ+1)2F1(p,q;1;A)+x12a12tEβ,2(a22tβ)+x12a12p=1q=0ap11aq22tpα+qβ+1Γ(pα+qβ+2)2F1(p,q;1;A)+a11Ap=0q=1ap11aq22t0(tτ)pα+qβ+α1g1(τ)dτΓ(pα+qβ+α)(p+1)2F1(p,1q;2;A)+a12t0(tτ)β1Eβ,β(a22(tτ)β)g2(τ)dτ+a12p=1q=0ap11aq22t0(tτ)pα+qβ+β1g2(τ)dτΓ(pα+qβ+β)2F1(p,q;1;A). (3.31)

    Using the Lemma 2.1 in [1], all the terms in Eqs (3.31) and (3.30) are equivalent, except for the third term of Eqs (3.31) and (3.30). Hence, we show here using some algebraic manipulation that the third term of Eq (3.31) is indeed equivalent to the third term in Eq (3.30). We have

    CDα(t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ)=CDα(t0(tτ)α1Γ(α)g1(τ)dτ+k=1t0ak11(tτ)kα+α1Γ(kα+α)g1(τ)dτ)=CDα(Iαg1(t))+CDα(k=1t0ak11(tτ)kα+α1Γ(kα+α)g1(τ)dτ)=g1(t)+k=1t0ak11(tτ)(k1)α+α1Γ((k1)α+α)g1(τ)dτ=g1(t)+a11k=0t0ak11(tτ)kα+α1Γ(kα+α)g1(τ)dτ=g1(t)+a11t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ. (3.32)

    By applying Lemma 2.1 in [1] and Eq (3.32), the first equation in (1.1) holds true, while the second equation in (1.1) can be verified using a similar approach. Hence, we have the theorem as follows:

    Theorem 1. The incommensurate fractional differential equation systems with fractional order 1<α,β<2 are given by:

    CDαx1(t)=a11x1(t)+a12x2(t)+g1(t),CDβx2(t)=a21x1(t)+a22x2(t)+g2(t), (3.33)

    with initial conditions x1(0)=x01,x2(0)=x02,x1(0)=x11,x2(0)=x12 and constant A=a12a21a11a22(1), a11,a220 have the solutions as follows:

    x1(t)=x01Eα(a11(1A)tα)+x11tEα,2(a11(1A)tα)+t0(tτ)α1Eα,α(a11(1A)(tτ)α)g1(τ)dτ+x01a11Ap=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+1)2F1(p,1q;2;A)+x11a11Ap=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+1)2F1(p,1q;2;A)+x02a12p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)2F1(p,q;1;A)+x12a12p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)2F1(p,q;1;A)+a11Ap=0q=0ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)(p+1)2F1(p,1q;2;A)+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)2F1(p,q;1;A), (3.34)
    x2(t)=x02Eβ(a22(1A)tβ)+x12tEβ,2(a22(1A)tβ)+t0(tτ)β1Eβ,β(a22(1A)(tτ)β)g2(τ)dτ+x02a22Ap=0q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)(q+1)2F1(1p,q;2;A)+x12a22Ap=0q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)(q+1)2F1(1p,q;2;A)+x01a21p=0q=0ap11aq22tpα+qβ+βΓ(pα+qβ+β+1)2F1(p,q;1;A)+x11a21p=0q=0ap11aq22tpα+qβ+β+1Γ(pα+qβ+β+2)2F1(p,q;1;A)+a21p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g1(τ)dτΓ(pα+qβ+α+β)2F1(p,q;1;A)+a22Ap=0q=0ap11aq22t0(tτ)pα+qβ+2β1g2(τ)dτΓ(pα+qβ+2β)(q+1)2F1(1p,q;2;A). (3.35)

    In this section, we will present some special cases of Theorem 1, which including the case when A = 1 , a_{11} = 0 , or a_{22} = 0 , respectively. In order to achieve these, we need the following lemmas involving bivariate Mittag-Leffler function:

    Lemma 1. [1] For \alpha, \beta > 0 and \gamma-1 > \lfloor \alpha \rfloor , we have

    \begin{equation} \frac{{\rm d}^{\alpha}}{{\rm d}t^{\alpha}} \left[ t^{\gamma-1}E_{\alpha,\beta,\gamma}(\lambda_{1}t^{\alpha},\lambda_{2}t^{\beta}) \right] = t^{\gamma-\alpha-1}E_{\alpha,\beta,\gamma-\alpha}(\lambda_{1}t^{\alpha},\lambda_{2}t^{\beta}), \end{equation} (4.1)

    for any t, \alpha, \beta, \gamma, \lambda_{1}, \lambda_{2}\in \mathbb{R} .

    Proof: See Lemma 2.2 in [1].

    Lemma 2. [1] For \alpha, \beta > 0 , we have

    \begin{equation} \begin{split} 1+a_{11}t^{\alpha}E_{\alpha,\beta, \alpha+1}(a_{11}t^{\alpha},a_{22}t^{\beta})+a_{22}t^{\beta}E_{\alpha,\beta, \beta+1}(a_{11}t^{\alpha},a_{22}t^{\beta})& = E_{\alpha,\beta, 1}(a_{11}t^{\alpha},a_{22}t^{\beta}),\\ \frac{t^{\alpha-1}}{\Gamma(\alpha)}+a_{11}t^{2\alpha-1}E_{\alpha,\beta, 2\alpha}(a_{11}t^{\alpha},a_{22}t^{\beta})+a_{22}t^{\alpha+\beta-1}E_{\alpha,\beta,\alpha+\beta}(a_{11}t^{\alpha},a_{22}t^{\beta})& = t^{\alpha-1}E_{\alpha,\beta, \alpha}(a_{11}t^{\alpha},a_{22}t^{\beta}),\\ \frac{t^{\beta-1}}{\Gamma(\beta)}+a_{11}t^{\alpha+\beta-1}E_{\alpha,\beta, \alpha+\beta}(a_{11}t^{\alpha},a_{22}t^{\beta})+a_{22}t^{2\beta-1}E_{\alpha,\beta,2\beta}(a_{11}t^{\alpha},a_{22}t^{\beta})& = t^{\beta-1}E_{\alpha,\beta,\beta}(a_{11}t^{\alpha},a_{22}t^{\beta}), \end{split} \end{equation} (4.2)

    for any t, \alpha, \beta \in \mathbb{R} .

    Proof: See [1].

    In this case, we have the hypergeometric function with A = 1 , i.e. a_{11}a_{22} = a_{12}a_{21} . The following identities are important for finding the explicit analytical solution of system (1.1).

    \begin{equation} \begin{split} _2F_1(-p,-q;1;1) = &\; \frac{\Gamma(1)\Gamma(p+q+1)}{\Gamma(p+1)\Gamma(q+1)} = \frac{\Gamma(p+q+1)}{\Gamma(p+1)\Gamma(q+1)} = \begin{pmatrix} p+q\\q \end{pmatrix}\\ (p+1)_2F_1(-p,1-q;2;1) = &\; (p+1)\frac{\Gamma(2)\Gamma(p+q+1)}{\Gamma(p+2)\Gamma(q+1)} = \frac{\Gamma(p+q+1)}{\Gamma(p+1)\Gamma(q+1)} = \begin{pmatrix} p+q\\q \end{pmatrix}. \end{split} \end{equation} (4.3)

    Using Eqs (3.24), (3.26), A = 1 and the identities in Eq (4.3), we can express x_1(t) as follows:

    \begin{array}{l} x_1(t) = x^0_1 E_\alpha(a_{11}t^\alpha)+x^1_1 tE_{\alpha,2}(a_{11}t^{\alpha})+\int_0^t (t-\tau)^{\alpha-1}E_{\alpha,\alpha}(a_{11}(t-\tau)^{\alpha})g_1(\tau)\; {\rm d}\tau \\ + x^0_1a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty} \frac{a_{11}^p a_{22}^q t^{p\alpha+q\beta+\alpha}}{\Gamma(p\alpha+q\beta+\alpha+1)} \begin{pmatrix} p+q \nonumber\\ q \end{pmatrix}+x^1_1a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha+1}}{\Gamma(p\alpha+q\beta+\alpha+2)}\begin{pmatrix} p+q \nonumber\\ q \end{pmatrix} \\ +x^0_2a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha}}{\Gamma(p\alpha+q\beta+\alpha+1)}\begin{pmatrix} p+q \nonumber\\ q \end{pmatrix}+x^1_2a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha+1}}{\Gamma(p\alpha+q\beta+\alpha+2)}\begin{pmatrix} p+q \nonumber\\ q \end{pmatrix} \\ +a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty} \frac{a_{11}^{p}a_{22}^{q} \int_0^t (t-\tau)^{p\alpha+q\beta+2\alpha-1}g_1(\tau)\; {\rm d}\tau}{\Gamma(p\alpha+q\beta+2\alpha)}\begin{pmatrix} p+q\\q \end{pmatrix}\\ +a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{a_{11}^{p}a_{22}^{q}\int_0^t (t-\tau)^{p\alpha+q\beta+\alpha+\beta-1} g_2(\tau) \; {\rm d}\tau }{\Gamma(p\alpha+q\beta+\alpha+\beta)}\begin{pmatrix} p+q \nonumber\\ q \end{pmatrix}. \end{array} (4.4)

    Expanding the Mittag-Leffler function and bivariate Mittag-Leffler function in the first three terms of Eq (4.4) and rearranging the terms in RHS of (4.4) yields

    \begin{array}{l} x_1(t) = \; x^0_1+x^0_1 \sum^\infty_{p = 0} \frac{a_{11}^{p+1}t^{p\alpha+\alpha}}{\Gamma(p\alpha+\alpha+1)}+ x^0_1a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty} \frac{a_{11}^p a_{22}^q t^{p\alpha+q\beta+\alpha}}{\Gamma(p\alpha+q\beta+\alpha+1)}\begin{pmatrix} p+q\nonumber \\ q \end{pmatrix} \\+x^0_2a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha}}{\Gamma(p\alpha+q\beta+\alpha+1)}\begin{pmatrix} p+q \nonumber\\ q \end{pmatrix} \\ +x^1_1 t+x^1_1 t \sum^\infty_{p = 0} \frac{a_{11}^{p+1}t^{p\alpha+\alpha}}{\Gamma(p\alpha+\alpha+2)}+x^1_1a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha+1}}{\Gamma(p\alpha+q\beta+\alpha+2)}\begin{pmatrix} p+q \nonumber\\ q \end{pmatrix} \\ +x^1_2a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha+1}}{\Gamma(p\alpha+q\beta+\alpha+2)}\begin{pmatrix} p+q\nonumber\\ q \end{pmatrix} \\ +\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau+\sum^\infty_{p = 0} \frac{a_{11}^{p+1} \int_0^t (t-\tau)^{p\alpha+2\alpha-1}g_1(\tau)\; {\rm d}\tau}{\Gamma(p\alpha+2\alpha)}\\ +a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty}\frac{a_{11}^{p}a_{22}^{q} \int_0^t (t-\tau)^{p\alpha+q\beta+2\alpha-1}g_1(\tau)\; {\rm d}\tau}{\Gamma(p\alpha+q\beta+2\alpha)}\begin{pmatrix} p+q\\q \end{pmatrix} \\ +a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{a_{11}^{p}a_{22}^{q}\int_0^t (t-\tau)^{p\alpha+q\beta+\alpha+\beta-1} g_2(\tau) \; {\rm d}\tau }{\Gamma(p\alpha+q\beta+\alpha+\beta)}\begin{pmatrix} p+q \nonumber\\q \end{pmatrix} \\ = \; x^0_1+\left( x^0_1a_{11}+x^0_2a_{12} \right)\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{a_{11}^p a_{22}^q t^{p\alpha+q\beta+\alpha}}{\Gamma(p\alpha+q\beta+\alpha+1)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix} \\ +x^1_1 t+\left(x^1_1a_{11}+x^1_2a_{12}\right)\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta+\alpha+1}}{\Gamma(p\alpha+q\beta+\alpha+2)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix} \\ +\frac{\int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau}{\Gamma(\alpha)}+a_{11}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^{q} \int_0^t (t-\tau)^{p\alpha+q\beta+2\alpha-1}g_1(\tau)\; {\rm d}\tau}{\Gamma(p\alpha+q\beta+2\alpha)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix} \\ +a_{12}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{a_{11}^{p}a_{22}^{q}\int_0^t (t-\tau)^{p\alpha+q\beta+\alpha+\beta-1} g_2(\tau) \; {\rm d}\tau }{\Gamma(p\alpha+q\beta+\alpha+\beta)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix}. \end{array} (4.5)

    Rewriting some terms in the above equation using bivariate Mittag-Leffler function yields

    \begin{align} x_1(t) = &\; x^0_1+\left(x^0_1 a_{11}+x^0_2 a_{12}\right)t^{\alpha}E_{\alpha,\beta,\alpha+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+x^1_1 t+\left(x^1_1 a_{11}+x^1_2 a_{12}\right) t^{\alpha+1}E_{\alpha,\beta,\alpha+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau \\ &+ a_{11} \int^t_0 (t-\tau)^{2\alpha-1} E_{\alpha,\beta,2\alpha}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau \\ &+ a_{12}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{align} (4.6)

    Using the similar approach for x_2(t) , we then have the following theorem:

    Theorem 2. The incommensurate fractional differential equation systems with fractional order 1 < \alpha, \beta < 2

    \begin{equation*} \label{maintheorem2} \begin{split} ^{C}{\rm D}^{\alpha}x_1(t)& = a_{11}x_1(t)+a_{12}x_2(t)+g_1(t),\\ ^{C}{\rm D}^{\beta}x_2(t)& = a_{21}x_1(t)+a_{22}x_2(t)+g_2(t), \end{split} \end{equation*}

    with initial conditions x_1(0) = x^0_1, \; x_2(0) = x^0_2, \; x'_1(0) = x^1_1, \; x'_2(0) = x^1_2 and constant A = \frac{a_{12}a_{21}}{a_{11}a_{22}} = 1 has the following solutions given by:

    \begin{equation} \begin{split} x_1(t) = &\; x^0_1+\left(x^0_1 a_{11}+x^0_2 a_{12}\right)t^{\alpha}E_{\alpha,\beta,\alpha+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+x^1_1 t+\left(x^1_1 a_{11}+x^1_2 a_{12}\right) t^{\alpha+1}E_{\alpha,\beta,\alpha+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau\\ &+ a_{11} \int^t_0 (t-\tau)^{2\alpha-1} E_{\alpha,\beta,2\alpha}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau, \end{split} \end{equation} (4.7)
    \begin{equation} \begin{split} x_2(t) = &\; x^0_2+\left(x^0_1 a_{21}+x^0_2 a_{22}\right)t^{\beta}E_{\alpha,\beta,\beta+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+x^1_2 t+\left(x^1_1 a_{21}+x^1_2 a_{22}\right) t^{\beta+1}E_{\alpha,\beta,\beta+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+\frac{1}{\Gamma(\beta)} \int^t_0 (t-\tau)^{\beta-1} g_2(\tau)\; {\rm d}\tau\\ &+ a_{21}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{22} \int^t_0 (t-\tau)^{2\beta-1} E_{\alpha,\beta,2\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{split} \end{equation} (4.8)

    Proof: The solution is proved when the first equation of (1.1) is satisfied. Hence, using Eqs (4.7) and (4.8), the LHS of (1.1) (after taking the fractional derivative for Eq (4.7) with Lemma 1) and RHS of (1.1) are shown as follows:

    \begin{align} ^C{\rm D}^\alpha x_1(t) = &\; \left(x^0_1 a_{11}+x^0_2 a_{12}\right)E_{\alpha,\beta,1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+\left(x^1_1 a_{11}+x^1_2 a_{12}\right) tE_{\alpha,\beta,2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+g_1(t)+ a_{11} \int^t_0 (t-\tau)^{\alpha-1} E_{\alpha,\beta,\alpha}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}\int^t_0 (t-\tau)^{\beta-1} E_{\alpha,\beta,\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau, \end{align} (4.9)
    \begin{align} a_{11}x_1+a_{12}x_2+g_1(t) = &\; a_{11}x^0_1+a_{11}\left(x^0_1 a_{11}+x^0_2 a_{12}\right)t^{\alpha}E_{\alpha,\beta,\alpha+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+x^1_1 a_{11} t+a_{11}\left(x^1_1 a_{11}+x^1_2 a_{12}\right) t^{\alpha+1}E_{\alpha,\beta,\alpha+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+\frac{a_{11}}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau\\ &+ a_{11}a_{11} \int^t_0 (t-\tau)^{2\alpha-1} E_{\alpha,\beta,2\alpha}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{11}a_{12}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau \\ &+a_{12}x^0_2+a_{12}\left(x^0_1 a_{21}+x^0_2 a_{22}\right)t^{\beta}E_{\alpha,\beta,\beta+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+x^1_2 a_{12}t+a_{12}\left(x^1_1 a_{21}+x^1_2 a_{22}\right) t^{\beta+1}E_{\alpha,\beta,\beta+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ &+ a_{12}a_{21}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+\frac{a_{12}}{\Gamma(\beta)} \int^t_0 (t-\tau)^{\beta-1} g_2(\tau)\; {\rm d}\tau\\ &+a_{12} a_{22} \int^t_0 (t-\tau)^{2\beta-1} E_{\alpha,\beta,2\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau+g_1(t). \end{align} (4.10)

    In order to verify the LHS (i.e. Eq (4.9)) is equal to RHS (i.e. Eq (4.10)) for the first equation in (1.1), we will compare the terms containing x^0_1 , x^0_2 , x^1_1 , x^1_2 , g_1(\tau) and g_2(\tau) . First, we take part of (4.10) involving x^0_1 to prove its equivalence with the corresponding x^0_1 term in (4.9). Since A = 1 , then a_{11}a_{22} = a_{12}a_{21} which yields

    \begin{array}{l} x^0_1\Big(a_{11}+a^2_{11}t^{\alpha}E_{\alpha,\beta,\alpha+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+a_{12}a_{21}t^{\beta}E_{\alpha,\beta,\beta+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\Big)\\ = x^0_1 a_{11}\Big(1+a_{11}t^{\alpha}E_{\alpha,\beta,\alpha+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+a_{22}t^{\beta}E_{\alpha,\beta,\beta+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\Big)\\ = x^0_1 a_{11}\Bigg(1+\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p+1}a_{22}^qt^{p\alpha+q\beta+\alpha}}{\Gamma(p\alpha+q\beta+\alpha+1)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix}+\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^{q+1}t^{p\alpha+q\beta+\beta}}{\Gamma(p\alpha+q\beta+\beta+1)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix} \Bigg) \\ = x^0_1 a_{11}\Bigg(1+\sum\limits_{p = 1}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^qt^{p\alpha+q\beta}}{\Gamma(p\alpha+q\beta+1)}\begin{pmatrix} p+q-1\nonumber\\q \end{pmatrix}+\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 1}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta}}{\Gamma(p\alpha+q\beta+1)}\begin{pmatrix} p+q-1\nonumber\\q-1 \end{pmatrix} \Bigg) \\ = x^0_1 a_{11}\Bigg(1+\sum\limits_{p = 1}^{\infty}\frac{a_{11}^{p}t^{p\alpha}}{\Gamma(p\alpha+1)}+\sum\limits_{q = 1}^{\infty}\frac{a_{22}^qt^{q\beta}}{\Gamma(q\beta+1)} +\sum\limits_{p = 1}^{\infty}\sum\limits_{q = 1}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta}}{\Gamma(p\alpha+q\beta+1)}\Bigg[\begin{pmatrix} p+q-1\nonumber\\q \end{pmatrix}+\begin{pmatrix} p+q-1\nonumber\\q-1 \end{pmatrix}\Bigg] \Bigg) \\ = x^0_1 a_{11}\Bigg(\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{a_{11}^{p}a_{22}^{q}t^{p\alpha+q\beta}}{\Gamma(p\alpha+q\beta+1)}\begin{pmatrix} p+q\nonumber\\q \end{pmatrix} \Bigg) \\ = x^0_1a_{11}E_{\alpha,\beta,1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right). \end{array} (4.11)

    Indeed, the above expression can be obtained via Lemma 2. Hence, using a similar approach, we have the proof for the terms with x^0_2 as follows:

    \begin{equation} \begin{split} &x^0_2a_{12} \Big(1+a_{11}t^{\alpha}E_{\alpha,\beta,\alpha+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)+a_{22}t^{\beta}E_{\alpha,\beta,\beta+1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\Big)\\ & = x^0_2a_{12} E_{\alpha,\beta,1}\left(a_{11}t^{\alpha},a_{22}t^\beta\right). \end{split} \end{equation} (4.12)

    Similarity, for the terms with x^1_1 and x^1_2 , the LHS is equal to the RHS since we have

    \begin{equation} \begin{split} &x^1_1 a_{11} t+a_{11}\left(x^1_1 a_{11}+x^1_2 a_{12}\right) t^{\alpha+1}E_{\alpha,\beta,\alpha+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right) +x^1_2 a_{12}t\\ &+a_{12}\left(x^1_1 a_{21}+x^1_2 a_{22}\right) t^{\beta+1}E_{\alpha,\beta,\beta+2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right)\\ = &\left(x^1_1 a_{11}+x^1_2 a_{12}\right) tE_{\alpha,\beta,2}\left(a_{11}t^{\alpha},a_{22}t^\beta\right). \end{split} \end{equation} (4.13)

    For the terms containing g_1(\tau) and g_2(\tau) , using Lemma 2 and A = 1 , we show that it is equivalent via the following equations.

    \begin{equation} \begin{split} &\frac{a_{11}}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau\\ &+ a^2_{11} \int^t_0 (t-\tau)^{2\alpha-1} E_{\alpha,\beta,2\alpha}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}a_{21}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ & = a_{11}\int^t_0 (t-\tau)^{\alpha-1} E_{\alpha,\beta,\alpha}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau, \end{split} \end{equation} (4.14)

    and

    \begin{equation} \begin{split} &a_{12}\Bigg( \frac{1}{\Gamma(\beta)} \int^t_0 (t-\tau)^{\beta-1} g_2(\tau)\; {\rm d}\tau\\ &+a_{11}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha,\beta,\alpha+\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau \\ &+a_{22} \int^t_0 (t-\tau)^{2\beta-1} E_{\alpha,\beta,2\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau\Bigg)\\ & = a_{12}\int^t_0 (t-\tau)^{\beta-1} E_{\alpha,\beta,\beta}\left(a_{11}(t-\tau)^{\alpha},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{split} \end{equation} (4.15)

    Theorem 2 is verified since the Eq (4.9) is equivalent to Eq (4.10) for each of the terms containing x^0_1 , x^0_2 , x^1_1 , x^1_2 , g_1(\tau) and g_2(\tau) .

    We have emphasized that a_{11} and a_{22} are not equal to zero in Theorem 1. However, we can still make assumption for these special cases. For the case a_{11} = 0 , we consider \alpha, \beta \in (1, 2) . The incommensurate fractional differential equation system (1.1) is now given by

    \begin{align} ^{C}{\rm D}^{\alpha}x_1(t)& = a_{12}x_2(t)+g_1(t),\\ ^{C}{\rm D}^{\beta}x_2(t)& = a_{21}x_1(t)+a_{22}x_2(t)+g_2(t), \end{align} (4.16)

    with the same initial conditions as (1.1).

    Since a_{11} = 0 , we use Eq (3.18) which makes it double series when j = 0 , yielding

    \begin{align} x_1(t) = &\; x^0_1+ x^0_1\sum\limits_{k = 0}^{\infty}\sum\limits_{m = 0}^{\infty} \frac{a^{k+1}_{12}a^{k+1}_{21}a_{22}^m t^{(k+1)\alpha+(m+k+1)\beta}}{\Gamma((k+1)\alpha+(m+k+1)\beta+1)}\frac{(k+m)!}{k!m!}\\ &+x^1_1 t+x^1_1\sum\limits_{k = 0}^{\infty}\sum\limits_{m = 0}^{\infty}\frac{a^{k+1}_{12}a^{k+1}_{21}a_{22}^{m}t^{(k+1)\alpha+(m+k+1)\beta+1}}{\Gamma((k+1)\alpha+(m+k+1)\beta+2)}\frac{(k+m)!}{k!m!}\\ &+x^0_2\sum\limits_{k = 0}^{\infty}\sum\limits_{m = 0}^{\infty} \frac{a^{k+1}_{12}a^k_{21}a_{22}^{m}t^{(k+1)\alpha+(m+k)\beta}}{\Gamma((k+1)\alpha+(m+k)\beta+1)}\frac{(k+m)!}{k!m!}\\ &+x^1_2\sum\limits_{k = 0}^{\infty}\sum\limits_{m = 0}^{\infty}\frac{a^{k+1}_{12}a^k_{21}a_{22}^{m}t^{(k+1)\alpha+(m+k)\beta+1}}{\Gamma((k+1)\alpha+(m+k)\beta+2)}\frac{(k+m)!}{k!m!}\\ &+\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}g_1(\tau) \; {\rm d}\tau\\ &+\sum\limits_{k = 0}^{\infty}\sum\limits_{m = 0}^{\infty} \frac{a^{k+1}_{12}a^{k+1}_{21}a_{22}^{m} \int_0^t (t-\tau)^{(k+2)\alpha+(m+k+1)\beta-1}g_1(\tau)\; {\rm d}\tau}{\Gamma((k+2)\alpha+(m+k+1)\beta)}\frac{(k+m)!}{k!m!}\\ &+\sum\limits_{k = 0}^{\infty}\sum\limits_{m = 0}^{\infty} \frac{a^{k+1}_{12}a^k_{21}a_{22}^{m}\int_0^t (t-\tau)^{(k+1)\alpha+(m+k+1)\beta-1} g_2(\tau) \; {\rm d}\tau }{\Gamma((k+1)\alpha+(m+k+1)\beta)}\frac{(k+m)!}{k!m!}. \end{align} (4.17)

    Rewriting the above equation using bivariate Mittag-Leffler function yields

    \begin{align} x_1(t) = &\; x^0_1+x^0_1 a_{12}a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\beta,\alpha+\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^1_1 t+x^1_1 a_{12}a_{21}t^{\alpha+\beta+1}E_{\alpha+\beta,\beta,\alpha+\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2 a_{12}t^{\alpha}E_{\alpha+\beta,\beta,\alpha+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right) +x^1_2 a_{12}t^{\alpha+1}E_{\alpha+\beta,\beta,\alpha+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau \\ &+ a_{12}a_{21} \int^t_0 (t-\tau)^{2\alpha+\beta-1} E_{\alpha+\beta,\beta,2\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau \\ &+ a_{12}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{align} (4.18)

    Since a_{11} = 0 , the x_2(t) solution can be obtained directly from the first equation of (4.16). By rearranging the first equation of (4.16), we obtain

    \begin{align} x_2(t) = &\; \frac{^{C}{\rm D}^{\alpha}x_1(t)}{a_{12}}-\frac{g_1(t)}{a_{12}}\\ = &\; \frac{1}{a_{12}}\Bigg(x^0_1 a_{12}a_{21}t^{\beta}E_{\alpha+\beta,\beta,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_1 a_{12}a_{21}t^{\beta+1}E_{\alpha+\beta,\beta,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2 a_{12}E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_2 a_{12}tE_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+g_1(t)\\ &+ a_{12}a_{21} \int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau \\ &+ a_{12}\int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau\Bigg)-\frac{g_1(t)}{a_{12}}. \end{align} (4.19)

    Using Lemma 1, and since our \alpha, \beta > 1 , we obtain the x_2(t) as follows:

    \begin{align} x_2(t) = &\; x^0_1 a_{21}t^{\beta}E_{\alpha+\beta,\beta,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_1 a_{21}t^{\beta+1}E_{\alpha+\beta,\beta,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2 E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_2 tE_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+ a_{21} \int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau \\ &+ \int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{align} (4.20)

    Hence, we can obtain the following theorem.

    Theorem 3. For special case a_{11} = 0 , the system (1.1) can be written as

    \begin{equation*} \label{main_a11 = 0new} \begin{split} ^{C}{\rm D}^{\alpha}x_1(t)& = a_{12}x_2(t)+g_1(t),\\ ^{C}{\rm D}^{\beta}x_2(t)& = a_{21}x_1(t)+a_{22}x_2(t)+g_2(t), \end{split} \end{equation*}

    and the explicit analytical solution of the above system with initial conditions x_1(0) = x^0_1 , x_2(0) = x^0_2 , x'_1(0) = x^1_1 , x'_2(0) = x^1_2 is given by:

    \begin{equation} \begin{split} x_1(t) = &\; x^0_1+x^0_1 a_{12}a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\beta,\alpha+\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^1_1 t+x^1_1 a_{12}a_{21}t^{\alpha+\beta+1}E_{\alpha+\beta,\beta,\alpha+\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2 a_{12}t^{\alpha}E_{\alpha+\beta,\beta,\alpha+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^1_2 a_{12}t^{\alpha+1}E_{\alpha+\beta,\beta,\alpha+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}a_{21} \int^t_0 (t-\tau)^{2\alpha+\beta-1} E_{\alpha+\beta,\beta,2\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau, \end{split} \end{equation} (4.21)
    \begin{equation} \begin{split} x_2(t) = &\; x^0_1 a_{21}t^{\beta}E_{\alpha+\beta,\beta,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_1 a_{21}t^{\beta+1}E_{\alpha+\beta,\beta,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2 E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_2 t E_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+a_{21} \int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+\int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{split} \end{equation} (4.22)

    Proof: The theorem can be checked by substituting the solutions into the second equation of (4.16). For the LHS, take the fractional derivative for Eq (4.22) and use Lemma 1, which yields

    \begin{equation} \begin{split} ^C{\rm D}^\beta x_2(t) = &\; x^0_1 a_{21}E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_1 a_{21}tE_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2 \; ^C{\rm D}^\beta \left[E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right]+x^1_2 \; ^C{\rm D}^\beta \left[t E_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right]\\ &+a_{21} \int^t_0 (t-\tau)^{\alpha-1} E_{\alpha+\beta,\beta,\alpha}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+{^C{\rm D}^\beta} \left[\int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau\right].\\ \end{split} \end{equation} (4.23)

    For the RHS, substitute Eqs (4.21) and (4.22) into the second equation of (4.16) yields

    \begin{equation} \begin{split} &a_{21}x_1(t)+a_{22}x_2(t)+g_2(t)\\ = &\; x^0_1a_{21}\left(1+a_{12}a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\beta,\alpha+\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right)\\ &+x^1_1 a_{21}t\left(1+a_{12} a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\beta,\alpha+\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right)\\ &+x^0_2 a_{12}a_{21}t^{\alpha}E_{\alpha+\beta,\beta,\alpha+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^1_2 a_{12}a_{21}t^{\alpha+1}E_{\alpha+\beta,\beta,\alpha+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+\frac{a_{21}}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}a^2_{21} \int^t_0 (t-\tau)^{2\alpha+\beta-1} E_{\alpha+\beta,\beta,2\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+ a_{12}a_{21}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau\\ &+x^0_1 a_{21}a_{22}t^{\beta}E_{\alpha+\beta,\beta,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_1 a_{21}a_{22}t^{\beta+1}E_{\alpha+\beta,\beta,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+x^0_2a_{22} E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+x^1_2 a_{22}t E_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\\ &+a_{21}a_{22} \int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau\\ &+a_{22}\int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{split} \end{equation} (4.24)

    Similar approach will be employed in proving Theorem 2, where we will be comparing one by one of the terms containing x^0_1 , x^0_2 , x^1_1 , x^1_2 , g_1(\tau) and g_2(\tau) . First, for the parts involving x^0_1 , x^1_1 and g_1(t) , by using Lemma 2, we prove Eq (4.24) to be equivalent as those corresponding parts in Eq (4.23) as follows:

    \begin{align} &x^0_1a_{21}\left(1+a_{12}a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\beta,\alpha+\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+a_{22}t^{\beta}E_{\alpha+\beta,\beta,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right)\\ = &\; x^0_1a_{21}\left[1+\sum^{\infty}_{m = 1}\sum^{\infty}_{n = 0} \frac{a_{12}^m a_{21}^m a_{22}^n t^{m\alpha+m\beta+n\beta}{\binom{m+n-1}{n}}}{\Gamma(m\alpha+m\beta+n\beta+1)}+\sum^\infty_{m = 0}\sum^\infty_{n = 1} \frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta}{\binom{m+n-1}{n-1}}}{\Gamma(m\alpha+m\beta+n\beta+1)}\right]\\ = &\; x^0_1 a_{21}E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right), \end{align} (4.25)
    \begin{align} &x^1_1 a_{21}t\left(1+a_{12} a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\beta,\alpha+\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+a_{22}t^{\beta}E_{\alpha+\beta,\beta,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right)\\ = &\; x^1_1a_{21}t\left[1+\sum^\infty_{m = 1}\sum^\infty_{n = 0} \frac{a_{12}^m a_{21}^m a_{22}^n t^{m\alpha+m\beta+n\beta}\binom{m+n-1}{n}}{\Gamma(m\alpha+m\beta+n\beta+2)} +\sum^\infty_{m = 0}\sum^\infty_{n = 1} \frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta}\binom{m+n-1}{n-1}}{\Gamma(m\alpha+m\beta+n\beta+2)}\right]\\ = &\; x^1_1 a_{21}tE_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right), \end{align} (4.26)
    \begin{align} &a_{21}\Bigg[\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau \\ &+a_{12}a_{21} \int^t_0 (t-\tau)^{2\alpha+\beta-1} E_{\alpha+\beta,\beta,2\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau \\ &+a_{22} \int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau \Bigg] \\ = &\; a_{21}\Bigg[\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau \\ &+\sum^\infty_{m = 0}\sum^\infty_{n = 0} \frac{a^{m+1}_{12}a^{m+1}_{21}a^n_{22}\int_0^t (t-\tau)^{m\alpha+m\beta+n\beta+2\alpha+\beta-1}g_1(\tau)\; {\rm d}\tau \; \binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta+2\alpha+\beta)}\\ &+\sum^\infty_{m = 0}\sum^\infty_{n = 0} \frac{a^m_{12}a^m_{21}a^{n+1}_{22}\int_0^t (t-\tau)^{m\alpha+m\beta+n\beta+\alpha+\beta-1}g_1(\tau)\; {\rm d}\tau \; \binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta+\alpha+\beta)} \Bigg]\\ = &\; a_{21}\Bigg[\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau \\ &+\sum^\infty_{m = 1}\sum^\infty_{n = 0} \frac{a^m_{12}a^m_{21}a^n_{22}\int_0^t (t-\tau)^{m\alpha+m\beta+n\beta+\alpha-1}g_1(\tau)\; {\rm d}\tau \; \binom{m+n-1}{n}}{\Gamma(m\alpha+m\beta+n\beta+\alpha)}\\ &+\sum^\infty_{m = 0}\sum^\infty_{n = 1} \frac{a^m_{12}a^m_{21}a^n_{22}\int_0^t (t-\tau)^{m\alpha+m\beta+n\beta+\alpha-1}g_1(\tau)\; {\rm d}\tau \; \binom{m+n-1}{n-1} }{\Gamma(m\alpha+m\beta+n\beta+\alpha)} \Bigg]\\ = &\; a_{21}\Bigg[\frac{1}{\Gamma(\alpha)} \int^t_0 (t-\tau)^{\alpha-1} g_1(\tau)\; {\rm d}\tau \\ &+\sum^\infty_{m = 1} \frac{a^m_{12}a^m_{21}\int_0^t (t-\tau)^{m\alpha+m\beta+\alpha-1}g_1(\tau)\; {\rm d}\tau}{\Gamma(m\alpha+m\beta+\alpha)}+\sum^\infty_{n = 1} \frac{a^n_{22}\int_0^t (t-\tau)^{n\beta+\alpha-1}g_1(\tau)\; {\rm d}\tau}{\Gamma(n\beta+\alpha)}\\ &+\sum^\infty_{m = 1}\sum^\infty_{n = 1} \frac{a^m_{12}a^m_{21}a^n_{22}\int_0^t (t-\tau)^{m\alpha+m\beta+n\beta+\alpha-1}g_1(\tau)\; {\rm d}\tau \; \left[ \binom{m+n-1}{n}+ \binom{m+n-1}{n-1}\right] }{\Gamma(m\alpha+m\beta+n\beta+\alpha)}\Bigg]\\ = &\; a_{21}\int^t_0 (t-\tau)^{\alpha-1}E_{\alpha+\beta,\beta,\alpha}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_1(\tau)\; {\rm d}\tau. \end{align} (4.27)

    Meanwhile, for x^0_2 , x^1_2 and g_2(t) parts, we proceed the proving from (4.23) as follows:

    \begin{align} &x^0_2 \; ^C{\rm D}^\beta \left[E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right]\\ = &\; x^0_2 \underset{(m,n)\neq(0,0)}{\sum^\infty_{m,n = 0}}\frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta-\beta}\binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta-\beta+1)} \\ = &\; x^0_2 \Bigg[\sum^\infty_{m = 1}\sum^\infty_{n = 0}\frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta-\beta}\binom{m+n-1}{n}}{\Gamma(m\alpha+m\beta+n\beta-\beta+1)} +\sum^\infty_{m = 0}\sum^\infty_{n = 1}\frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta-\beta}\binom{m+n-1}{n-1}}{\Gamma(m\alpha+m\beta+n\beta-\beta+1)} \Bigg]\\ = &\; x^0_2 \Bigg[\sum^\infty_{m = 0}\sum^\infty_{n = 0}\frac{a^{m+1}_{12}a^{m+1}_{21}a^n_{22}t^{(m+1)\alpha+m\beta+n\beta}\binom{m+n}{n}}{\Gamma((m+1)\alpha+m\beta+n\beta+1)}+\sum^\infty_{m = 0}\sum^\infty_{n = 0}\frac{a^m_{12}a^m_{21}a^{n+1}_{22}t^{m\alpha+m\beta+n\beta}\binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta+1)} \Bigg]\\ = &\; x^0_2\left(a_{12}a_{21}t^\alpha E_{\alpha+\beta,\beta,\alpha+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+a_{22}E_{\alpha+\beta,\beta,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right), \end{align} (4.28)
    \begin{align} &x^1_2 \; ^C{\rm D}^\beta \left[t E_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right]\\ = &\; x^1_2 \underset{(m,n)\neq(0,0)}{\sum^\infty_{m,n = 0}}\frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta-\beta+1}\binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta-\beta+2)}\\ = &\; x^1_2 \Bigg[\sum^\infty_{m = 1}\sum^\infty_{n = 0}\frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta-\beta+1}\binom{m+n-1}{n}}{\Gamma(m\alpha+m\beta+n\beta-\beta+2)} +\sum^\infty_{m = 0}\sum^\infty_{n = 1}\frac{a^m_{12}a^m_{21}a^n_{22}t^{m\alpha+m\beta+n\beta-\beta+1}\binom{m+n-1}{n-1}}{\Gamma(m\alpha+m\beta+n\beta-\beta+2)} \Bigg]\\ = &\; x^1_2 \Bigg[\sum^\infty_{m = 0}\sum^\infty_{n = 0}\frac{a^{m+1}_{12}a^{m+1}_{21}a^n_{22}t^{(m+1)\alpha+m\beta+n\beta+1}\binom{m+n}{n}}{\Gamma((m+1)\alpha+m\beta+n\beta+2)} +\sum^\infty_{m = 0}\sum^\infty_{n = 0}\frac{a^m_{12}a^m_{21}a^{n+1}_{22}t^{m\alpha+m\beta+n\beta+1}\binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta+2)} \Bigg]\\ = &\; x^1_2 \left(a_{12}a_{21}t^{\alpha+1} E_{\alpha+\beta,\beta,\alpha+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)+a_{22}tE_{\alpha+\beta,\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{22}t^\beta\right)\right), \end{align} (4.29)
    \begin{align} &^C{\rm D}^\beta \left[\int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau\right]\\ = &\underset{(m,n)\neq(0,0)}{\sum^\infty_{m,n = 0}}\frac{a^m_{12}a^m_{21}a^n_{22}\int^t_0 (t-\tau)^{m\alpha+m\beta+n\beta-1} g_2(\tau)\; {\rm d}\tau \; \binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta)} \\ = &\sum^\infty_{m = 1}\sum^\infty_{n = 0}\frac{a^m_{12}a^m_{21}a^n_{22}\int^t_0 (t-\tau)^{m\alpha+m\beta+n\beta-1} g_2(\tau)\; {\rm d}\tau \; \binom{m+n-1}{n}}{\Gamma(m\alpha+m\beta+n\beta)}\\ &+\sum^\infty_{m = 0}\sum^\infty_{n = 1}\frac{a^m_{12}a^m_{21}a^n_{22}\int^t_0 (t-\tau)^{m\alpha+m\beta+n\beta-1} g_2(\tau)\; {\rm d}\tau \; \binom{m+n-1}{n-1}}{\Gamma(m\alpha+m\beta+n\beta)} \\ = &\sum^\infty_{m = 0}\sum^\infty_{n = 0}\frac{a^{m+1}_{12}a^{m+1}_{21}a^n_{22}\int^t_0 (t-\tau)^{m\alpha+m\beta+n\beta+\alpha+\beta-1} g_2(\tau)\; {\rm d}\tau\; \binom{m+n}{n}}{\Gamma((m\alpha+m\beta+n\beta+\alpha+\beta)}\\ & +\sum^\infty_{m = 0}\sum^\infty_{n = 0}\frac{a^m_{12}a^m_{21}a^{n+1}_{22}\int^t_0 (t-\tau)^{m\alpha+m\beta+n\beta+\beta-1} g_2(\tau)\; {\rm d}\tau\; \binom{m+n}{n}}{\Gamma(m\alpha+m\beta+n\beta+\beta)}\\ = &\; a_{12}a_{21}\int^t_0 (t-\tau)^{\alpha+\beta-1} E_{\alpha+\beta,\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau\\ &+a_{22}\int^t_0 (t-\tau)^{\beta-1} E_{\alpha+\beta,\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{22}(t-\tau)^\beta\right) g_2(\tau)\; {\rm d}\tau. \end{align} (4.30)

    Since all the terms in Eqs (4.23) and (4.24) are equivalent, the solution of the system (4.16) is verified.

    For the case a_{22} = 0 , we consider \alpha, \beta \in (1, 2) . The incommensurate fractional differential equation system (1.1) is now given by

    \begin{equation} \begin{split} ^{C}{\rm D}^{\alpha}x_1(t)& = a_{11}x_1(t)+a_{12}x_2(t)+g_1(t),\\ ^{C}{\rm D}^{\beta}x_2(t)& = a_{21}x_1(t)+g_2(t), \end{split} \end{equation} (4.31)

    with the same initial conditions as in Eq (1.1).

    For the case a_{22} = 0 , using similar approach from the previous subsection, we obtain the following solution

    \begin{align} x_1(t) = &\; x^0_1 E_{\alpha+\beta,\alpha,1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)+x^1_1 t E_{\alpha+\beta,\alpha,2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)\\ &+x^0_2 a_{12}t^{\alpha}E_{\alpha+\beta,\alpha,\alpha+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)+x^1_2 a_{12}t^{\alpha+1}E_{\alpha+\beta,\alpha,\alpha+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)\\ &+\int^t_0 (t-\tau)^{\alpha-1}g_1(\tau) E_{\alpha+\beta,\alpha,\alpha}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{11}(t-\tau)^\alpha\right) \; {\rm d}\tau\\ &+a_{12} \int^t_0 (t-\tau)^{\alpha+\beta-1} g_2(\tau) E_{\alpha+\beta,\alpha,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{11}(t-\tau)^\alpha\right) \; {\rm d}\tau, \end{align} (4.32)
    \begin{align} x_2(t) = &\; x^0_1 a_{21}t^{\beta}E_{\alpha+\beta,\alpha,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)+x^1_1 a_{21}t^{\beta+1}E_{\alpha+\beta,\alpha,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)\\ &+x^0_2+x^0_2 a_{12}a_{21}t^{\alpha+\beta}E_{\alpha+\beta,\alpha,\alpha+\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)\\ &+x^1_2 t+x^1_2 a_{12}a_{21}t^{\alpha+\beta+1}E_{\alpha+\beta,\alpha,\alpha+\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta},a_{11}t^\alpha\right)\\ &+ a_{21}\int^t_0 (t-\tau)^{\alpha+\beta-1} g_1(\tau) E_{\alpha+\beta,\alpha,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{11}(t-\tau)^\alpha\right) \; {\rm d}\tau\\ &+\frac{1}{\Gamma(\beta)} \int^t_0 (t-\tau)^{\beta-1} g_2(\tau)\; {\rm d}\tau\\ &+a_{12}a_{21} \int^t_0 (t-\tau)^{\alpha+2\beta-1}g_2(\tau) E_{\alpha+\beta,\alpha,\alpha+2\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta},a_{11}(t-\tau)^\alpha\right) \; {\rm d}\tau. \end{align} (4.33)

    For the case a_{11} = 0 and a_{22} = 0 , we consider \alpha, \beta \in (1, 2) . The incommensurate fractional differential equation system (1.1) is now given by

    \begin{equation} \begin{split} ^{C}{\rm D}^{\alpha}x_1(t)& = a_{12}x_2(t)+g_1(t),\\ ^{C}{\rm D}^{\beta}x_2(t)& = a_{21}x_1(t)+g_2(t), \end{split} \end{equation} (4.34)

    with the same initial conditions as in Eq (1.1).

    For the case a_{11} = 0 and a_{22} = 0 , using similar approach from the previous subsection, we obtain the following solution,

    \begin{equation} \begin{split} x_1(t) = &\; x^0_1 E_{\alpha+\beta}\left(a_{12}a_{21}t^{\alpha+\beta}\right)+x^1_1 t E_{\alpha+\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta}\right)+x^0_2 a_{12}t^{\alpha}E_{\alpha+\beta,\alpha+1}\left(a_{12}a_{21}t^{\alpha+\beta}\right)\\ &+x^1_2 a_{12}t^{\alpha+1}E_{\alpha+\beta,\alpha+2}\left(a_{12}a_{21}t^{\alpha+\beta}\right)+\int^t_0 (t-\tau)^{\alpha-1}g_1(\tau) E_{\alpha+\beta,\alpha}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta}\right) \; {\rm d}\tau\\ &+a_{12}\int^t_0 (t-\tau)^{\alpha+\beta-1}g_2(\tau) E_{\alpha+\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta}\right) \; {\rm d}\tau, \end{split} \end{equation} (4.35)
    \begin{equation} \begin{split} x_2(t) = &\; x^0_1 a_{21}t^{\beta}E_{\alpha+\beta,\beta+1}\left(a_{12}a_{21}t^{\alpha+\beta}\right)+x^1_1 a_{21}t^{\beta+1}E_{\alpha+\beta,\beta+2}\left(a_{12}a_{21}t^{\alpha+\beta}\right)+x^0_2 E_{\alpha+\beta}\left(a_{12}a_{21}t^{\alpha+\beta}\right)\\&+x^1_2 t E_{\alpha+\beta,2}\left(a_{12}a_{21}t^{\alpha+\beta}\right)+a_{21} \int^t_0 (t-\tau)^{\alpha+\beta-1} g_1(\tau) E_{\alpha+\beta,\alpha+\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta}\right) \; {\rm d}\tau\\ &+\int^t_0 (t-\tau)^{\beta-1} g_2(\tau) E_{\alpha+\beta,\beta}\left(a_{12}a_{21}(t-\tau)^{\alpha+\beta}\right) \; {\rm d}\tau. \end{split} \end{equation} (4.36)

    This section illustrates an explicit analytical solution for the incommensurate fractional differential system for order \alpha, \beta \in (1, 2) using the theorems that we had derived in previous sections. Four examples will be presented using Theorems 1–3, respectively, for the case of a_{11} = 0 and a_{22} = 0 .

    Example 1. Consider the following incommensurate linear fractional differential equation system

    \begin{equation} \begin{pmatrix} \frac{{\rm d}^{3/2}x_1}{{\rm d}t^{3/2}} \\ \frac{{\rm d}^{4/3}x_2}{{\rm d}t^{4/3}} \end{pmatrix} = \begin{pmatrix} 3&2 \\ 1&2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} + \begin{pmatrix} 2t \\ -4t \end{pmatrix}, \end{equation} (5.1)

    with respect to the initial conditions x_1(0) = 3/2 , x_1'(0) = 2 , x_2(0) = 2 and x_2'(0) = 3 .

    Solution: Using Theorem 1, the explicit analytical solution is obtained as follows:

    \begin{align} x_1(t) = &\; \frac{3}{2} E_{\frac{3}{2}}\left(2t^{\frac{3}{2}}\right)+2 tE_{\frac{3}{2},2}\left(2t^{\frac{3}{2}}\right)+2 t^{\frac{5}{2}}E_{\frac{3}{2},\frac{7}{2}}\left(2t^{\frac{3}{2}}\right)\\ &+\frac{3}{2}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^p 2^q t^{\frac{3p}{2}+\frac{4q}{3}+\frac{3}{2}}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{5}{2}\right)} (p+1) _2F_1\left(-p,1-q;2;\frac{1}{3}\right)\\ &+2\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{5}{2}}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{7}{2}\right)}(p+1) _2F_1\left(-p,1-q;2;\frac{1}{3}\right)\\ &+4\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{3}{2}}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{5}{2}\right)}\; _2F_1\left(-p,-q;1;\frac{1}{3}\right)\\ &+6\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{5}{2}}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{7}{2}\right)}\; _2F_1\left(-p,-q;1;\frac{1}{3}\right)\\ &+2\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+4}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+5\right)}(p+1) _2F_1\left(-p,1-q;2;\frac{1}{3}\right) \\ &-8\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{23}{6}}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{29}{6}\right)}\; _2F_1\left(-p,-q;1;\frac{1}{3}\right), \end{align} (5.2)
    \begin{align} x_2(t) = &\; 2 E_\frac{4}{3}\left(\frac{4}{3}t^\frac{4}{3}\right)+3 tE_{\frac{4}{3},2}\left(\frac{4}{3}t^{\frac{4}{3}}\right)-4t^{\frac{7}{3}}E_{\frac{4}{3},\frac{10}{3}}\left(\frac{4}{3}t^{\frac{4}{3}}\right)\\ &+\frac{4}{3}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^p 2^q t^{\frac{3p}{2}+\frac{4q}{3}+\frac{4}{3}}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{7}{3}\right)} (q+1) _2F_1\left(1-p,-q;2;\frac{1}{3}\right)\\ &+2\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{7}{3}}}{\Gamma(\frac{3p}{2}+\frac{4q}{3}+\frac{10}{3})}(q+1) _2F_1\left(1-p,-q;2;\frac{1}{3}\right)\\ &+\frac{3}{2}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{4}{3}}}{\Gamma(\frac{3p}{2}+\frac{4q}{3}+\frac{7}{3})}\; _2F_1\left(-p,-q;1;\frac{1}{3}\right)\\ &+2\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty}\frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{7}{3}}}{\Gamma(\frac{3p}{2}+\frac{4q}{3}+\frac{10}{3})}\; _2F_1\left(-p,-q;1;\frac{1}{3}\right)\\ &+2\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^{p}2^{q}t^{\frac{3p}{2}+\frac{4q}{3}+\frac{23}{6}}}{\Gamma(\frac{3p}{2}+\frac{4q}{3}+\frac{29}{6})}\; _2F_1\left(-p,-q;1;\frac{1}{3}\right)\\ &-\frac{8}{3}\sum\limits_{p = 0}^{\infty}\sum\limits_{q = 0}^{\infty} \frac{3^{p}2^{q} t^{\frac{3p}{2}+\frac{4q}{3}+\frac{11}{3}}}{\Gamma(\frac{3p}{2}+\frac{4q}{3}+\frac{14}{3})}(q+1) _2F_1\left(1-p,-q;2;\frac{1}{3}\right). \end{align} (5.3)

    Example 2. Consider the following incommensurate linear fractional differential equation system

    \begin{equation} \begin{pmatrix} \frac{{\rm d}^{3/2}x_1}{{\rm d}t^{3/2}} \\ \frac{{\rm d}^{4/3}x_2}{{\rm d}t^{4/3}} \end{pmatrix} = \begin{pmatrix} 3/2&3 \\ 1&2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} + \begin{pmatrix} \cos(t) \\ e^{t} \end{pmatrix}, \end{equation} (5.4)

    with respect to the initial conditions x_1(0) = 3 , x_1'(0) = -1 , x_2(0) = 1 and x_2'(0) = 2 .

    Solution: This example is the case when A = 1 . Hence, using Theorem 2, we have the following explicit analytical solution

    \begin{align} x_1(t) = &\; 3-t+\frac{15}{2}t^{\frac{3}{2}}E_{\frac{3}{2},\frac{4}{3},\frac{5}{2}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+\frac{9}{2} t^{\frac{5}{2}}E_{\frac{3}{2},\frac{4}{3},\frac{7}{2}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+\frac{1}{\Gamma(\frac{3}{2})} \int^t_0 (t-\tau)^{\frac{1}{2}} \cos(\tau)\; {\rm d}\tau \\ &+ \frac{3}{2} \int^t_0 (t-\tau)^{2} E_{\frac{3}{2},\frac{4}{3},3}\left(\frac{3}{2}(t-\tau)^{\frac{3}{2}},2(t-\tau)^{\frac{4}{3}}\right) \cos(\tau)\; {\rm d}\tau \\ &+3\int^t_0 (t-\tau)^{\frac{11}{6}} E_{\frac{3}{2},\frac{4}{3},\frac{17}{6}}\left(\frac{3}{2}(t-\tau)^{\frac{3}{2}},2(t-\tau)^{\frac{4}{3}}\right) e^\tau\; {\rm d}\tau \\ = &\; 3-t+\frac{15}{2}t^{\frac{3}{2}}E_{\frac{3}{2},\frac{4}{3},\frac{5}{2}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+\frac{9}{2} t^{\frac{5}{2}}E_{\frac{3}{2},\frac{4}{3},\frac{7}{2}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right) \\ &+\frac{\sqrt{t}}{3\Gamma(\frac{3}{2})} \left[3\sin(t) \; _1F_2 \left(\frac{1}{4};\frac{1}{2},\frac{5}{4};\frac{-t^2}{4}\right)-t \cos(t) \; _1F_2 \left(\frac{3}{4};\frac{3}{2},\frac{7}{4};\frac{-t^2}{4}\right)\right] \\ &+ \frac{3}{2} \sum^{\infty}_{p,q = 0}\frac{\left(\frac{3}{2}\right)^{p}\left(2\right)^{q}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+4\right)}\frac{(p+q)!}{p!q!}t^{\frac{3p}{2}+\frac{4q}{3}+3}\; _1F_2 \left(1;\frac{3p}{4}+\frac{2q}{3}+2,\frac{3p}{4}+\frac{2q}{3}+\frac{5}{2};\frac{-t^2}{4}\right) \\ &+3e^t\sum\limits_{p,q = 0}^{\infty}\frac{\left(\frac{3}{2}\right)^{p}\left(2\right)^{q}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{17}{6}\right)}\; \gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{17}{6},t\right), \end{align} (5.5)
    \begin{align} x_2(t) = &\; 1+2t+5t^{\frac{4}{3}}E_{\frac{3}{2},\frac{4}{3},\frac{7}{3}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+3t^{\frac{7}{3}}E_{\frac{3}{2},\frac{4}{3},\frac{10}{3}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+\frac{1}{\Gamma(\frac{4}{3})} \int^t_0 (t-\tau)^{\frac{1}{3}}e^{\tau}\; {\rm d}\tau \\ &+ \int^t_0 (t-\tau)^{\frac{11}{6}} E_{\frac{3}{2},\frac{4}{3},\frac{17}{6}}\left(\frac{3}{2}(t-\tau)^{\frac{3}{2}},2(t-\tau)^{\frac{4}{3}}\right) \cos(\tau)\; {\rm d}\tau \\ &+ 2\int^t_0 (t-\tau)^{\frac{5}{3}} E_{\frac{3}{2},\frac{4}{3},\frac{8}{3}}\left(\frac{3}{2}(t-\tau)^{\frac{3}{2}},2(t-\tau)^{\frac{4}{3}}\right) e^{\tau}\; {\rm d}\tau \\ = &\; 1+2t+5t^{\frac{4}{3}}E_{\frac{3}{2},\frac{4}{3},\frac{7}{3}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+3t^{\frac{7}{3}}E_{\frac{3}{2},\frac{4}{3},\frac{10}{3}}\left(\frac{3}{2}t^{\frac{3}{2}},2t^{\frac{4}{3}}\right)+\frac{e^t}{\Gamma(\frac{4}{3})} \; \gamma\left(\frac{4}{3},t\right)\\ &+\sum^{\infty}_{p,q = 0}\left(1.5\right)^{p}\left(2\right)^{q}t^{\frac{9p+8q+17}{6}}\frac{(p+q)!}{p!q!}\Bigg[\frac{1}{\Gamma\left(\frac{9p+8q+23}{6}\right)}-\frac{t^2\; _1F_2 \left(1;\frac{9p+8q+35}{12},\frac{9p+8q+41}{12};\frac{-t^2}{4}\right)}{\Gamma\left(\frac{9p+8q+35}{6}\right)}\Bigg]\\ &+2e^t\sum\limits_{p,q = 0}^{\infty}\frac{\left(\frac{3}{2}\right)^{p}\left(2\right)^{q}}{\Gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{8}{3}\right)}\; \gamma\left(\frac{3p}{2}+\frac{4q}{3}+\frac{8}{3},t\right). \end{align} (5.6)

    Example 3. Consider the following incommensurate linear fractional differential equation system

    \begin{equation} \begin{pmatrix} \frac{{\rm d}^{1.2}x_1}{{\rm d}t^{1.2}} \\ \frac{{\rm d}^{1.5}x_2}{{\rm d}t^{1.5}} \end{pmatrix} = \begin{pmatrix} 0&1 \\ 5&1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} + \begin{pmatrix} t+e^t \\ e^{2t} \end{pmatrix}, \end{equation} (5.7)

    with respect to the initial conditions x_1(0) = 1 , x_1'(0) = \Gamma(1.8) , x_2(0) = 1 and x_2'(0) = 3 .

    Solution: Since a_{11} = 0 , Theorem 3 will be applied. For sake of simplicity, we present this solution using decimal numbers. This example have the following explicit analytical solution

    \begin{equation} \begin{split} x_1(t) = &\; 1+t\Gamma(1.8)+5t^{2.7}E_{2.7,1.5,3.7}\left(5t^{2.7},t^{1.5}\right)+5\Gamma(1.8)t^{3.7}E_{2.7,1.5,4.7}\left(5t^{2.7},t^{1.5}\right)\\ &+t^{1.2}E_{2.7,1.5,2.2}\left(5t^{2.7},t^{1.5}\right)+3t^{2.2}E_{2.7,1.5,3.2}\left(5t^{2.7},t^{1.5}\right)+\frac{1}{\Gamma(1.2)} \int^t_0 (t-\tau)^{0.2} \left(\tau+e^\tau\right)\; {\rm d}\tau\\ &+ 5 \int^t_0 (t-\tau)^{2.9} E_{2.7,1.5,3.9}\left(5(t-\tau)^{2.7},(t-\tau)^{1.5}\right) \left(\tau+e^\tau\right)\; {\rm d}\tau\\ &+ \int^t_0 (t-\tau)^{1.7} E_{2.7,1.5,2.7}\left(5(t-\tau)^{2.7},(t-\tau)^{1.5}\right) e^{2\tau}\; {\rm d}\tau\\ = &\; 1+t\Gamma(1.8)+5t^{2.7}E_{2.7,1.5,3.7}\left(5t^{2.7},t^{1.5}\right)+5\Gamma(1.8)t^{3.7}E_{2.7,1.5,4.7}\left(5t^{2.7},t^{1.5}\right)\\&+t^{1.2}E_{2.7,1.5,2.2}\left(5t^{2.7},t^{1.5}\right)+3t^{2.2}E_{2.7,1.5,3.2}\left(5t^{2.7},t^{1.5}\right)\\ &+\frac{t^{2.2}}{\Gamma(3.2)}+\frac{e^t \gamma(1.2,t)}{\Gamma(1.2)}+ 5t^{4.9} E_{2.7,1.5,5.9}\left(5t^{2.7},t^{1.5}\right)\\ &+ 5e^t\sum^\infty_{p,q = 0}\frac{5^{p}\; \gamma(2.7p+1.5q+3.9,t)}{\Gamma(2.7p+1.5q+3.9)} + e^{2t}\sum^\infty_{p,q = 0}\frac{5^{p}\gamma(2.7p+1.5q+2.7,2t)}{2^{2.7p+1.5q+2.7}\Gamma(2.7p+1.5q+2.7)}, \end{split} \end{equation} (5.8)
    \begin{equation} \begin{split} x_2(t) = & \; 5t^{1.5}E_{2.7,1.5,2.5}\left(5t^{2.7},t^{1.5}\right)+5\Gamma(1.8)t^{2.5}E_{2.7,1.5,3.5}\left(5t^{2.7},t^{1.5}\right)\\ &+E_{2.7,1.5,1}\left(5t^{2.7},t^{1.5}\right)+3t E_{2.7,1.5,2}\left(5t^{2.7},t^{1.5}\right)\\ &+5 \int^t_0 (t-\tau)^{1.7} E_{2.7,1.5,2.7}\left(5(t-\tau)^{2.7},(t-\tau)^{1.5}\right) \left(\tau+e^\tau\right)\; {\rm d}\tau\\ &+\int^t_0 (t-\tau)^{0.5} E_{2.7,1.5,1.5}\left(5(t-\tau)^{2.7},(t-\tau)^{1.5}\right) e^{2\tau}\; {\rm d}\tau\\ = &\; 5t^{1.5}E_{2.7,1.5,2.5}\left(5t^{2.7},t^{1.5}\right)+5\Gamma(1.8)t^{2.5}E_{2.7,1.5,3.5}\left(5t^{2.7},t^{1.5}\right)\\ &+E_{2.7,1.5,1}\left(5t^{2.7},t^{1.5}\right)+3t E_{2.7,1.5,2}\left(5t^{2.7},t^{1.5}\right)+5 t^{3.7} E_{2.7,1.5,4.7}\left(5t^{2.7},t^{1.5}\right) \\&+ 5e^t\sum^\infty_{p,q = 0}\frac{5^{p}\; \gamma(2.7p+1.5q+2.7,t)}{\Gamma(2.7p+1.5q+2.7)}+e^{2t}\sum^\infty_{p,q = 0}\frac{ 5^{p}\gamma(2.7p+1.5q+1.5,2t)}{2^{2.7p+1.5q+1.5}\Gamma(2.7p+1.5q+1.5)}. \end{split} \end{equation} (5.9)

    Example 4. Consider the following incommensurate linear fractional differential equation system

    \begin{equation} \begin{pmatrix} \frac{{\rm d}^{1.2}x_1}{{\rm d}t^{1.2}} \\ \frac{{\rm d}^{1.6}x_2}{{\rm d}t^{1.6}} \end{pmatrix} = \begin{pmatrix} 0&3 \\ -1&0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} + \begin{pmatrix} 2.5t^{0.8}-9t \\ \Gamma\left(0.8\right)t^2 \end{pmatrix}, \end{equation} (5.10)

    with respect to the initial conditions x_1(0) = 2 , x_1'(0) = 0 , x_2(0) = -6 and x_2'(0) = 3 .

    Solution: Since a_{11} = a_{22} = 0 , applying the result presented in subsection 4.4 and with the help from Eq (2.9), we have the following explicit analytical solution

    \begin{equation} \begin{split} x_1(t) = &\; 2 E_{2.8}\left(-3t^{2.8}\right)-18t^{1.2}E_{2.8,2.2}\left(-3t^{2.8}\right)+9t^{2.2}E_{2.8,3.2}\left(-3t^{2.8}\right)\\ &+2.5\Gamma(1.8) t^{2} E_{2.8,3}\left(-3t^{2.8}\right)-9 t^{2.2} E_{2.8,3.2}\left(-3t^{2.8}\right)+3 \Gamma(0.8) t^{4.8}E_{2.8,5.8}\left(-3t^{2.8}\right), \end{split} \end{equation} (5.11)
    \begin{equation} \begin{split} x_2(t) = &-2 t^{1.6}E_{2.8,2.6}\left(-3t^{2.8}\right)-6 E_{2.8}\left(-3t^{2.8}\right)+3 t E_{2.8,2}\left(-3t^{2.8}\right)\\ &-2.5\Gamma(1.8) t^{3.6} E_{2.8,4.6}\left(-3t^{2.8}\right)+9 t^{3.8} E_{2.8,4.8}\left(-3t^{2.8}\right)+\Gamma(0.8) t^{3.6}E_{2.8,4.6}\left(-3t^{2.8}\right). \end{split} \end{equation} (5.12)

    The solution of this example is shown in Figure 1. For the purpose of validate the solution (i.e. LHS equal to RHS of the problem), we can find the LHS via fractional derivative of these x_1(t) and x_2(t) for the desired order (i.e. in this example, are 1.2 and 1.6 respectively). Meanwhile for the RHS, substitute the solution x_1(t) and x_2(t) in the RHS of problem. If the analytical expression is too lengthy, we suggest to plot the both sides up to desired power. We use Maple to perform all the computation.

    Figure 1.  Solution x_1(t) and x_2(t) for Example 4.

    This paper has successfully derived the explicit analytical solution of linear incommensurate fractional differential equation systems with fractional order 1 < \alpha, \beta < 2 . Using the new theorems, analytical solutions are obtained, and we presented them via some examples. This paper serves as an extension of the similar result recently achieved in [1,19], which limited to fractional order 0 < \alpha, \beta < 1 . Moreover, the analytical solution obtained in this paper may enable us to investigate more rigorously the stability analysis and asymptotic stability for incommensurate fractional differential equation systems with fractional order 1 < \alpha, \beta < 2 , especially when this kind of incommensurate system may be more suitable to represent the real-world applications such as COVID-19 [38], cancer modelling, fluid flows problems. It may also be extended to higher order in the future. Explicit analytical solution for higher order (i.e. \alpha, \beta > 2 ) incommensurate fractional differential equation systems may be obtained using a similar approach.

    Communication of this research is made possible through monetary assistance by Universiti Tun Hussein Onn Malaysia and the UTHM Publisher's Office via Publication Fund E15216. The first writer wish to thank UTHM for the financial support during her study through GPPS H049.

    Authors declare that there is no conflict of interests regarding the publication of the paper.



    [1] H. Kim, Geospatial data-driven assessment of earthquake-induced liquefaction impact mapping using classifier and cluster ensemble, Appl. Soft Comput., 140 (2023), 110266. https://doi.org/10.1016/j.asoc.2023.110266 doi: 10.1016/j.asoc.2023.110266
    [2] E. Ivannikova, H. Park, T. Hämäläinen, K. Lee, Revealing community structures by ensemble clustering using group diffusion, Inform. Fusion, 42 (2018), 24–36. https://doi.org/10.1016/j.inffus.2017.09.013 doi: 10.1016/j.inffus.2017.09.013
    [3] X. Zeng, A. Chen, M. Zhou, Color perception algorithm of medical images using density peak based hierarchical clustering, Biomed. Signal Proces., 48 (2019), 69–79. https://doi.org/10.1016/j.bspc.2018.09.013 doi: 10.1016/j.bspc.2018.09.013
    [4] Y. Slimen, S. Allio, J. Jacques, Model-based co-clustering for functional data, Neurocomputing, 291 (2018), 97–108. https://doi.org/10.1016/j.neucom.2018.02.055 doi: 10.1016/j.neucom.2018.02.055
    [5] Q. Zhang, C. Zhu, L. Yang, Z. Chen, L. Zhao, P. Li, An incremental cfs algorithm for clustering large data in industrial internet of things, IEEE T. Ind. Inform., 13 (2017), 1193–1201. https://doi.org/10.1109/TII.2017.2684807 doi: 10.1109/TII.2017.2684807
    [6] A. Fahad, N. Alshatri, Z. Tari, A. Alamri, I. Khalil, A. Zomaya, et al., A survey of clustering algorithms for big data: taxonomy and empirical analysis, IEEE T. Emerg. Top. Com., 2 (2014), 267–279. https://doi.org/10.1109/TETC.2014.2330519 doi: 10.1109/TETC.2014.2330519
    [7] D. Wang, T. Li, P. Deng, F. Zhang, W. Huang, P. Zhang, et al., A generalized deep learning clustering algorithm based on non-negative matrix factorization, ACM T. Knowl. Discov. D., 17 (2023), 99. https://doi.org/10.1145/3584862 doi: 10.1145/3584862
    [8] M. Shahzad, S. Riazul Islam, M. Hossain, M. Abdullah-Al-Wadud, A. Alamri, M. Hussain, Gafor: genetic algorithm based fuzzy optimized re-clustering in wireless sensor networks, Mathematics, 9 (2021), 43. https://doi.org/10.3390/math9010043 doi: 10.3390/math9010043
    [9] W. Zhao, C. Deng, C. Ngo, k-means: a revisit, Neurocomputing, 291 (2018), 195–206. https://doi.org/10.1016/j.neucom.2018.02.072 doi: 10.1016/j.neucom.2018.02.072
    [10] Y. Zhu, K. Ting, M. Carman, Density-ratio based clustering for discovering clusters with varying densities, Pattern Recogn., 60 (2016), 983–997. https://doi.org/10.1016/j.patcog.2016.07.007 doi: 10.1016/j.patcog.2016.07.007
    [11] Chaomurilige, How klfcm works—convergence and parameter analysis for klfcm clustering algorithm, Mathematics, 11 (2023), 2285. https://doi.org/10.3390/math11102285 doi: 10.3390/math11102285
    [12] H. Ling, J. Wu, Y. Zhou, W. Zheng, How many clusters? A robust pso-based local density mode, Neurocomputing, 207 (2016), 264–275. https://doi.org/10.1016/j.neucom.2016.03.071 doi: 10.1016/j.neucom.2016.03.071
    [13] A. Rodriguez, A. Laio, Clustering by fast search and find of density peaks, Science, 344 (2014), 1492–1496. https://doi.org/10.1126/science.1242072 doi: 10.1126/science.1242072
    [14] R. Liu, H. Wang, X. Yu, Shared-nearest-neighbor-based clustering by fast search and find of density peaks, Inform. Sciences, 450 (2018), 200–226. https://doi.org/10.1016/j.ins.2018.03.031 doi: 10.1016/j.ins.2018.03.031
    [15] X. Xu, S. Ding, Y. Wang, L. Wang, W. Jia, A fast density peaks clustering algorithm with sparse search, Inform. Sciences, 554 (2021), 61–83. https://doi.org/10.1016/j.ins.2020.11.050 doi: 10.1016/j.ins.2020.11.050
    [16] J. Xu, G. Wang, T. Li, W. Deng, G. Gou, Fat node leading tree for data stream clustering with density peaks, Knowl.-Based Syst., 120 (2017), 99–117. https://doi.org/10.1016/j.knosys.2016.12.025 doi: 10.1016/j.knosys.2016.12.025
    [17] S. Ding, M. Du, T. Sun, X. Xu, Y. Xue, An entropy-based density peaks clustering algorithm for mixed type data employing fuzzy neighborhood, Knowl.-Based Syst., 133 (2017), 294–313. https://doi.org/10.1016/j.knosys.2017.07.027 doi: 10.1016/j.knosys.2017.07.027
    [18] M. Karaayvaz, S. Cristea, S. Gillespie, A. Patel, R. Mylvaganam, C. Luo, et al., Unravelling subclonal heterogeneity and aggressive disease states in tnbc through single-cell rna-seq, Nat. Commun., 9 (2018), 3588. https://doi.org/10.1038/s41467-018-06052-0 doi: 10.1038/s41467-018-06052-0
    [19] X. Li, K. Wong, Evolutionary multiobjective clustering and its applications to patient stratification, IEEE T. Cybernetics, 49 (2019), 1680–1693. https://doi.org/10.1109/TCYB.2018.2817480 doi: 10.1109/TCYB.2018.2817480
    [20] T. Xu, J. Jiang, A graph adaptive density peaks clustering algorithm for automatic centroid selection and effective aggregation, Expert Syst. Appl., 195 (2022), 116539. https://doi.org/10.1016/j.eswa.2022.116539 doi: 10.1016/j.eswa.2022.116539
    [21] L. Bai, X. Cheng, J. Liang, H. Shen, Y. Guo, Fast density clustering strategies based on the k-means algorithm, Pattern Recogn., 71 (2017), 375–386. https://doi.org/10.1016/j.patcog.2017.06.023 doi: 10.1016/j.patcog.2017.06.023
    [22] J. Xu, G. Wang, W. Deng, Denpehc: density peak based efficient hierarchical clustering, Inform. Sciences, 373 (2016), 200–218. https://doi.org/10.1016/j.ins.2016.08.086 doi: 10.1016/j.ins.2016.08.086
    [23] J. Chen, H. He, A fast density-based data stream clustering algorithm with cluster centers self-determined for mixed data, Inform. Sciences, 345 (2016), 271–293. https://doi.org/10.1016/j.ins.2016.01.071 doi: 10.1016/j.ins.2016.01.071
    [24] Y. Liu, Z. Ma, F. Yu, Adaptive density peak clustering based on k-nearest neighbors with aggregating strategy, Knowl.-Based Syst., 133 (2017), 208–220. https://doi.org/10.1016/j.knosys.2017.07.010 doi: 10.1016/j.knosys.2017.07.010
    [25] M. Masud, J. Huang, C. Wei, J. Wang, I. Khan, M. Zhong, I-nice: a new approach for identifying the number of clusters and initial cluster centres, Inform. Sciences, 466 (2018), 129–151. https://doi.org/10.1016/j.ins.2018.07.034 doi: 10.1016/j.ins.2018.07.034
    [26] M. D'Errico, E. Facco, A. Laio, A Rodriguez, Automatic topography of high-dimensional data sets by non-parametric density peak clustering, Inform. Sciences, 560 (2021), 476–492. https://doi.org/10.1016/j.ins.2021.01.010 doi: 10.1016/j.ins.2021.01.010
    [27] P. Rousseeuw, Silhouettes: a graphical aid to the interpretation and validation of cluster analysis, J. Comput. Appl. Math., 20 (1987), 53–65. https://doi.org/10.1016/0377-0427(87)90125-7 doi: 10.1016/0377-0427(87)90125-7
    [28] L. Lovmar, A. Ahlford, M. Jonsson, A. Syvanen, Silhouette scores for assessment of SNP genotype clusters, BMC Genomics, 6 (2005), 35. https://doi.org/10.1186/1471-2164-6-35 doi: 10.1186/1471-2164-6-35
    [29] X. Xu, S. Ding, Z. Shi, An improved density peaks clustering algorithm with fast finding cluster centers, Knowl.-Based Syst., 158 (2018), 65–74. https://doi.org/10.1016/j.knosys.2018.05.034 doi: 10.1016/j.knosys.2018.05.034
    [30] J. Xie, H. Gao, W. Xie, X. Liu, P. Grant, Robust clustering by detecting density peaks and assigning points based on fuzzy weighted k-nearest neighbors, Inform. Sciences, 354 (2016), 19–40. https://doi.org/10.1016/j.ins.2016.03.011 doi: 10.1016/j.ins.2016.03.011
    [31] M. Du, S. Ding, H. Jia, Study on density peaks clustering based on k-nearest neighbors and principal component analysis, Knowl.-Based Syst., 99 (2016), 135–145. https://doi.org/10.1016/j.knosys.2016.02.001 doi: 10.1016/j.knosys.2016.02.001
    [32] S. Ding, C. Li, X. Xu, L. Ding, J. Zhang, L. Guo, et al., A sampling-based density peaks clustering algorithm for large-scale data, Pattern Recogn., 136 (2023), 109238. https://doi.org/10.1016/j.patcog.2022.109238 doi: 10.1016/j.patcog.2022.109238
    [33] Z. Liang, P. Chen, Delta-density based clustering with a divide-and-conquer strategy: 3dc clustering, Pattern Recogn. Lett., 73 (2016), 52–59. https://doi.org/10.1016/j.patrec.2016.01.009 doi: 10.1016/j.patrec.2016.01.009
    [34] M. Chen, L. Li, B. Wang, J. Cheng, L. Pan, X. Chen, Effectively clustering by finding density backbone based-on knn, Pattern Recogn., 60 (2016), 486–498. https://doi.org/10.1016/j.patcog.2016.04.018 doi: 10.1016/j.patcog.2016.04.018
    [35] M. Wang, F. Min, Z. Zhang, Y. Wu, Active learning through density clustering, Expert Syst. Appl., 85 (2017), 305–317. https://doi.org/10.1016/j.eswa.2017.05.046 doi: 10.1016/j.eswa.2017.05.046
    [36] B. Wu, B. Wilamowski, A fast density and grid based clustering method for data with arbitrary shapes and noise, IEEE T. Ind. Inform., 13 (2017), 1620–1628. https://doi.org/10.1109/TII.2016.2628747 doi: 10.1109/TII.2016.2628747
    [37] Z. Li, Y. Tang, Comparative density peaks clustering, Expert Syst. Appl., 95 (2018), 236–247. https://doi.org/10.1016/j.eswa.2017.11.020 doi: 10.1016/j.eswa.2017.11.020
    [38] K. Ting, Y. Zhu, M. Carman, Y. Zhu, Z. Zhou, Overcoming key weaknesses of distance-based neighbourhood methods using a data dependent dissimilarity measure, Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, 2016, 1205–1214. https://doi.org/10.1145/2939672.2939779 doi: 10.1145/2939672.2939779
    [39] S. Ding, W. Du, X. Xu, T. Shi, Y. Wang, C. Li, An improved density peaks clustering algorithm based on natural neighbor with a merging strategy, Inform. Sciences, 624 (2023), 252–276. https://doi.org/10.1016/j.ins.2022.12.078 doi: 10.1016/j.ins.2022.12.078
    [40] F. Samaria, A. Harter, Parameterisation of a stochastic model for human face identification, Proceedings of 1994 IEEE Workshop on Applications of Computer Vision, 1994,138–142. https://doi.org/10.1109/ACV.1994.341300 doi: 10.1109/ACV.1994.341300
    [41] B. Frey, D. Dueck, Clustering by passing messages between data points, Science, 315 (2007), 972–976. https://doi.org/10.1126/science.1136800 doi: 10.1126/science.1136800
    [42] D. Ienco, G. Bordogna, Fuzzy extensions of the DBScan clustering algorithm, Soft Comput., 22 (2018), 1719–1730. https://doi.org/10.1007/s00500-016-2435-0 doi: 10.1007/s00500-016-2435-0
    [43] J. Jiang, X. Yan, Z. Yu, J. Guo, W. Tian, A Chinese expert disambiguation method based on semi-supervised graph clustering, Int. J. Mach. Learn. Cyber., 6 (2015), 197–204. https://doi.org/10.1007/s13042-014-0255-z doi: 10.1007/s13042-014-0255-z
    [44] H. Jia, S. Ding, M. Du, Y. Xue, Approximate normalized cuts without eigen-decomposition, Inform. Sciences, 374 (2016), 135–150. https://doi.org/10.1016/j.ins.2016.09.032 doi: 10.1016/j.ins.2016.09.032
    [45] N. Vinh, J. Epps, J. Bailey, Information theoretic measures for clusterings comparison: is a correction for chance necessary? Proceedings of the 26th Annual International Conference on Machine Learning, 2009, 1073–1080. https://doi.org//10.1145/1553374.1553511 doi: 10.1145/1553374.1553511
    [46] M. Sampat, Z. Wang, S. Gupta, A. Bovik, M. Markey, Complex wavelet structural similarity: a new image similarity index, IEEE T. Image Process., 18 (2009), 2385–2401. https://doi.org/10.1109/TIP.2009.2025923 doi: 10.1109/TIP.2009.2025923
  • This article has been cited by:

    1. Mohammed A. Almalahi, Satish K. Panchal, Fahd Jarad, Phang Chang, Results on Implicit Fractional Pantograph Equations with Mittag-Leffler Kernel and Nonlocal Condition, 2022, 2022, 2314-4785, 1, 10.1155/2022/9693005
    2. Ahmad Qazza, Aliaa Burqan, Rania Saadeh, Fahd Jarad, Application of ARA-Residual Power Series Method in Solving Systems of Fractional Differential Equations, 2022, 2022, 1563-5147, 1, 10.1155/2022/6939045
    3. Thongchai Botmart, Ravi P. Agarwal, Muhammed Naeem, Adnan Khan, Rasool Shah, On the solution of fractional modified Boussinesq and approximate long wave equations with non-singular kernel operators, 2022, 7, 2473-6988, 12483, 10.3934/math.2022693
    4. A. I. Ahmed, T. A. Al-Ahmary, Sagheer Abbas, Fractional-Order Chelyshkov Collocation Method for Solving Systems of Fractional Differential Equations, 2022, 2022, 1563-5147, 1, 10.1155/2022/4862650
    5. Soon Hock Gan, Chang Phang, 2023, Chapter 12, 978-981-99-2849-1, 131, 10.1007/978-981-99-2850-7_12
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1541) PDF downloads(88) Cited by(4)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog