Processing math: 56%
Research article

Markov random fields model and applications to image processing

  • Received: 18 October 2021 Revised: 15 December 2021 Accepted: 15 December 2021 Published: 22 December 2021
  • MSC : 60G20, 60H10, 60K35, 81T15, 81T18

  • Markov random fields (MRFs) are well studied during the past 50 years. Their success are mainly due to their flexibility and to the fact that they gives raise to stochastic image models. In this work, we will consider a stochastic differential equation (SDE) driven by Lévy noise. We will show that the solution Xv of the SDE is a MRF satisfying the Markov property. We will prove that the Gibbs distribution of the process Xv can be represented graphically through Feynman graphs, which are defined as a set of cliques, then we will provide applications of MRFs in image processing where the image intensity at a particular location depends only on a neighborhood of pixels.

    Citation: Boubaker Smii. Markov random fields model and applications to image processing[J]. AIMS Mathematics, 2022, 7(3): 4459-4471. doi: 10.3934/math.2022248

    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
  • Markov random fields (MRFs) are well studied during the past 50 years. Their success are mainly due to their flexibility and to the fact that they gives raise to stochastic image models. In this work, we will consider a stochastic differential equation (SDE) driven by Lévy noise. We will show that the solution Xv of the SDE is a MRF satisfying the Markov property. We will prove that the Gibbs distribution of the process Xv can be represented graphically through Feynman graphs, which are defined as a set of cliques, then we will provide applications of MRFs in image processing where the image intensity at a particular location depends only on a neighborhood of pixels.



    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, a11=0, or a22=0, respectively. In order to achieve these, we need the following lemmas involving bivariate Mittag-Leffler function:

    Lemma 1. [1] For α,β>0 and γ1>α, we have

    dαdtα[tγ1Eα,β,γ(λ1tα,λ2tβ)]=tγα1Eα,β,γα(λ1tα,λ2tβ), (4.1)

    for any t,α,β,γ,λ1,λ2R.

    Proof: See Lemma 2.2 in [1].

    Lemma 2. [1] For α,β>0, we have

    1+a11tαEα,β,α+1(a11tα,a22tβ)+a22tβEα,β,β+1(a11tα,a22tβ)=Eα,β,1(a11tα,a22tβ),tα1Γ(α)+a11t2α1Eα,β,2α(a11tα,a22tβ)+a22tα+β1Eα,β,α+β(a11tα,a22tβ)=tα1Eα,β,α(a11tα,a22tβ),tβ1Γ(β)+a11tα+β1Eα,β,α+β(a11tα,a22tβ)+a22t2β1Eα,β,2β(a11tα,a22tβ)=tβ1Eα,β,β(a11tα,a22tβ), (4.2)

    for any t,α,βR.

    Proof: See [1].

    In this case, we have the hypergeometric function with A=1, i.e. a11a22=a12a21. The following identities are important for finding the explicit analytical solution of system (1.1).

    2F1(p,q;1;1)=Γ(1)Γ(p+q+1)Γ(p+1)Γ(q+1)=Γ(p+q+1)Γ(p+1)Γ(q+1)=(p+qq)(p+1)2F1(p,1q;2;1)=(p+1)Γ(2)Γ(p+q+1)Γ(p+2)Γ(q+1)=Γ(p+q+1)Γ(p+1)Γ(q+1)=(p+qq). (4.3)

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

    x1(t)=x01Eα(a11tα)+x11tEα,2(a11tα)+t0(tτ)α1Eα,α(a11(tτ)α)g1(τ)dτ+x01a11p=0q=1ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+qq)+x11a11p=0q=1ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+qq)+x02a12p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+qq)+x12a12p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+qq)+a11p=0q=1ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)(p+qq)+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)(p+qq). (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

    x1(t)=x01+x01p=0ap+111tpα+αΓ(pα+α+1)+x01a11p=0q=1ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+qq)+x02a12p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+qq)+x11t+x11tp=0ap+111tpα+αΓ(pα+α+2)+x11a11p=0q=1ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+qq)+x12a12p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+qq)+1Γ(α)t0(tτ)α1g1(τ)dτ+p=0ap+111t0(tτ)pα+2α1g1(τ)dτΓ(pα+2α)+a11p=0q=1ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)(p+qq)+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)(p+qq)=x01+(x01a11+x02a12)p=0q=0ap11aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+qq)+x11t+(x11a11+x12a12)p=0q=0ap11aq22tpα+qβ+α+1Γ(pα+qβ+α+2)(p+qq)+t0(tτ)α1g1(τ)dτΓ(α)+a11p=0q=0ap11aq22t0(tτ)pα+qβ+2α1g1(τ)dτΓ(pα+qβ+2α)(p+qq)+a12p=0q=0ap11aq22t0(tτ)pα+qβ+α+β1g2(τ)dτΓ(pα+qβ+α+β)(p+qq). (4.5)

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

    x1(t)=x01+(x01a11+x02a12)tαEα,β,α+1(a11tα,a22tβ)+x11t+(x11a11+x12a12)tα+1Eα,β,α+2(a11tα,a22tβ)+1Γ(α)t0(tτ)α1g1(τ)dτ+a11t0(tτ)2α1Eα,β,2α(a11(tτ)α,a22(tτ)β)g1(τ)dτ+a12t0(tτ)α+β1Eα,β,α+β(a11(tτ)α,a22(tτ)β)g2(τ)dτ. (4.6)

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

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

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

    with initial conditions x1(0)=x01,x2(0)=x02,x1(0)=x11,x2(0)=x12 and constant A=a12a21a11a22=1 has the following solutions given by:

    x1(t)=x01+(x01a11+x02a12)tαEα,β,α+1(a11tα,a22tβ)+x11t+(x11a11+x12a12)tα+1Eα,β,α+2(a11tα,a22tβ)+1Γ(α)t0(tτ)α1g1(τ)dτ+a11t0(tτ)2α1Eα,β,2α(a11(tτ)α,a22(tτ)β)g1(τ)dτ+a12t0(tτ)α+β1Eα,β,α+β(a11(tτ)α,a22(tτ)β)g2(τ)dτ, (4.7)
    x2(t)=x02+(x01a21+x02a22)tβEα,β,β+1(a11tα,a22tβ)+x12t+(x11a21+x12a22)tβ+1Eα,β,β+2(a11tα,a22tβ)+1Γ(β)t0(tτ)β1g2(τ)dτ+a21t0(tτ)α+β1Eα,β,α+β(a11(tτ)α,a22(tτ)β)g1(τ)dτ+a22t0(tτ)2β1Eα,β,2β(a11(tτ)α,a22(tτ)β)g2(τ)dτ. (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:

    CDαx1(t)=(x01a11+x02a12)Eα,β,1(a11tα,a22tβ)+(x11a11+x12a12)tEα,β,2(a11tα,a22tβ)+g1(t)+a11t0(tτ)α1Eα,β,α(a11(tτ)α,a22(tτ)β)g1(τ)dτ+a12t0(tτ)β1Eα,β,β(a11(tτ)α,a22(tτ)β)g2(τ)dτ, (4.9)
    a11x1+a12x2+g1(t)=a11x01+a11(x01a11+x02a12)tαEα,β,α+1(a11tα,a22tβ)+x11a11t+a11(x11a11+x12a12)tα+1Eα,β,α+2(a11tα,a22tβ)+a11Γ(α)t0(tτ)α1g1(τ)dτ+a11a11t0(tτ)2α1Eα,β,2α(a11(tτ)α,a22(tτ)β)g1(τ)dτ+a11a12t0(tτ)α+β1Eα,β,α+β(a11(tτ)α,a22(tτ)β)g2(τ)dτ+a12x02+a12(x01a21+x02a22)tβEα,β,β+1(a11tα,a22tβ)+x12a12t+a12(x11a21+x12a22)tβ+1Eα,β,β+2(a11tα,a22tβ)+a12a21t0(tτ)α+β1Eα,β,α+β(a11(tτ)α,a22(tτ)β)g1(τ)dτ+a12Γ(β)t0(tτ)β1g2(τ)dτ+a12a22t0(tτ)2β1Eα,β,2β(a11(tτ)α,a22(tτ)β)g2(τ)dτ+g1(t). (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 x01, x02, x11, x12, g1(τ) and g2(τ). First, we take part of (4.10) involving x01 to prove its equivalence with the corresponding x01 term in (4.9). Since A=1, then a11a22=a12a21 which yields

    x01(a11+a211tαEα,β,α+1(a11tα,a22tβ)+a12a21tβEα,β,β+1(a11tα,a22tβ))=x01a11(1+a11tαEα,β,α+1(a11tα,a22tβ)+a22tβEα,β,β+1(a11tα,a22tβ))=x01a11(1+p=0q=0ap+111aq22tpα+qβ+αΓ(pα+qβ+α+1)(p+qq)+p=0q=0ap11aq+122tpα+qβ+βΓ(pα+qβ+β+1)(p+qq))=x01a11(1+p=1q=0ap11aq22tpα+qβΓ(pα+qβ+1)(p+q1q)+p=0q=1ap11aq22tpα+qβΓ(pα+qβ+1)(p+q1q1))=x01a11(1+p=1ap11tpαΓ(pα+1)+q=1aq22tqβΓ(qβ+1)+p=1q=1ap11aq22tpα+qβΓ(pα+qβ+1)[(p+q1q)+(p+q1q1)])=x01a11(p=0q=0ap11aq22tpα+qβΓ(pα+qβ+1)(p+qq))=x01a11Eα,β,1(a11tα,a22tβ). (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] S. Albeverio, L. Dipersio, E. Mastrogiacomo, B. Smii, Invariant measures for SDEs driven by Lévy noise: A case study for dissipative nonlinear drift in infinite dimension, Commun. Math. Sci., 15 (2017), 957–983. https://doi.org/10.4310/CMS.2017.v15.n4.a3
    [2] S. Albeverio, L. Dipersio, E. Mastrogiacomo, B. Smii, A class of Lévy driven SDEs and their explicit invariant measures, Potential Anal., 45 (2016), 229–259. https://doi.org/10.1007/s11118-016-9544-3 doi: 10.1007/s11118-016-9544-3
    [3] S. Albeverio, E. Mastrogiacomo, B. Smii, Small noise asymptotic expansions for stochastic PDE's driven by dissipative nonlinearity and Lévy noise, Stoch. Proc. Appl., 123 (2013), 2084–2109. https://doi.org/10.1016/j.spa.2013.01.013 doi: 10.1016/j.spa.2013.01.013
    [4] S. Albeverio, B. Smii, Borel summation of the small time expansion of some SDE's driven by Gaussian white noise, Asymptotic Anal., 114 (2019), 211–223. https://doi.org/10.3233/ASY-191525 doi: 10.3233/ASY-191525
    [5] S. Albeverio, B. Smii, Asymptotic expansions for SDE's with small multiplicative noise, Stoch. Proc. Appl., 125 (2015), 1009–1031. https://doi.org/10.1016/j.spa.2014.09.009 doi: 10.1016/j.spa.2014.09.009
    [6] J. Besag, Spatial interaction and the statistical analysis of lattice systems, J. Royal Stat. Soc. Ser. B, 36 (1974), 192–236.
    [7] A. Blake, P. Kohli, C. Rother, Markov random fields for vision and image processing, The MIT Press, 2011.
    [8] P. Brémaud, Markov chains, Gibbs fields, Monte Carlo simulation and queues, Springer-Verlag, 1999.
    [9] Z. Brze\acute{z}niak, E. Hausenblas, Uniqueness in law of the It\hat{o} integral with respect to Lévy noise, In: R. Dalang, M. Dozzi, F. Russo, Seminar on stochastic analysis, random fields and applications VI, Vol. 63, Basel: Springer, 2011. https://doi.org/10.1007/978-3-0348-0021-1_3
    [10] Y. Fu, Y. Kang, G. Chen, Stochastic resonance based visual perception using spiking neural networks, Front. Comput. Neurosci., 14 (2020), 24. https://doi.org/10.3389/fncom.2020.00024 doi: 10.3389/fncom.2020.00024
    [11] S. Geman, D. Geman, Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images, IEEE T. Pattern Anal., PAMI-6 (1984), 721–741. https://doi.org/10.1109/TPAMI.1984.4767596 doi: 10.1109/TPAMI.1984.4767596
    [12] D. Griffeath, Introduction to random fields, In: Denumerable Markov chains, Graduate Texts in Mathematics, New York: Springer, 1976. https://doi.org/10.1007/978-1-4684-9455-6_12
    [13] H. Gottschalk, B. Smii, How to determine the law of the solution to a SPDE driven by a Lévy space-time noise, J. Math. Phys., 43 (2007), 1–22.
    [14] P. Hartman, A. Wintner, On the infinitisimal generators of integral convolutions, Am. J. Math., 64 (1942), 273–298. https://doi.org/10.2307/2371683 doi: 10.2307/2371683
    [15] K. He, Y. Li, S. Soundarajanc, J. E. Hopcroft, Hidden community detection in social networks, Inform. Sciences, 425 (2018), 92–106. https://doi.org/10.1016/j.ins.2017.10.019 doi: 10.1016/j.ins.2017.10.019
    [16] R. Kinderman, J. L. Snell, Markov random fields and their applications, Contemporary Mathematics, 1980. http://dx.doi.org/10.1090/conm/001
    [17] P. Kr\ddot{a}henb\ddot{a}hl, V. Koltun, Efficient inference in fully connected CRFs with gaussian edge potentials, Adv. Neural Inf. Process. Syst., 24 (2011), 109–117.
    [18] D. Mugnolo, Gaussian estimates for a heat equation on a network, Netw. Heterog. Media, 2 (2007), 55–79. https://doi.org/10.3934/nhm.2007.2.55 doi: 10.3934/nhm.2007.2.55
    [19] D. Mugnolo, S. Romanelli, Dynamic and generalized Wentzell node conditions for network equations, Math. Methods Appl. Sci., 30 (2007), 681–706. https://doi.org/10.1002/mma.805 doi: 10.1002/mma.805
    [20] K. I. Sato, Lévy processes and infinitely divisible distributions, Cambridge University Press, 1999.
    [21] B. Smii, Asymptotic expansion of the transition density of the semigroup associated to a SDE driven by Lévy noise, Asymptotic Anal., 124 (2021), 51–68. https://doi.org/10.3233/ASY-201640 doi: 10.3233/ASY-201640
    [22] C. Turchetti, Stochastic models of neural networks, IOS Press, 2004.
    [23] V. K. Pandey, H. Agarwal, A. K. Aggarwal Image solution of stochastic differential equation of diffusion type driven by Brownian motion, Singapore: Springer, 2021,542–553. https://doi.org/10.1007/978-981-16-1092-9_46
  • 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
  • © 2022 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(2867) PDF downloads(134) Cited by(3)

Figures and Tables

Figures(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog