Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article Special Issues

A spectral collocation method for the coupled system of nonlinear fractional differential equations

  • This paper analyzes the coupled system of nonlinear fractional differential equations involving the caputo fractional derivatives of order α(1,2) on the interval (0, T). Our method of analysis is based on the reduction of the given system to an equivalent system of integral equations, then the resulting equation is discretized by using a spectral method based on the Legendre polynomials. We have constructed a Legendre spectral collocation method for the coupled system of nonlinear fractional differential equations. The error bounds under the L2 and Lnorms is also provided, then the theoretical result is validated by a number of numerical tests.

    Citation: Xiaojun Zhou, Yue Dai. A spectral collocation method for the coupled system of nonlinear fractional differential equations[J]. AIMS Mathematics, 2022, 7(4): 5670-5689. doi: 10.3934/math.2022314

    Related Papers:

    [1] A. H. Tedjani, A. Z. Amin, Abdel-Haleem Abdel-Aty, M. A. Abdelkawy, Mona Mahmoud . Legendre spectral collocation method for solving nonlinear fractional Fredholm integro-differential equations with convergence analysis. AIMS Mathematics, 2024, 9(4): 7973-8000. doi: 10.3934/math.2024388
    [2] M. A. Zaky, M. Babatin, M. Hammad, A. Akgül, A. S. Hendy . Efficient spectral collocation method for nonlinear systems of fractional pantograph delay differential equations. AIMS Mathematics, 2024, 9(6): 15246-15262. doi: 10.3934/math.2024740
    [3] Younes Talaei, Sanda Micula, Hasan Hosseinzadeh, Samad Noeiaghdam . A novel algorithm to solve nonlinear fractional quadratic integral equations. AIMS Mathematics, 2022, 7(7): 13237-13257. doi: 10.3934/math.2022730
    [4] Imran Talib, Md. Nur Alam, Dumitru Baleanu, Danish Zaidi, Ammarah Marriyam . A new integral operational matrix with applications to multi-order fractional differential equations. AIMS Mathematics, 2021, 6(8): 8742-8771. doi: 10.3934/math.2021508
    [5] Khalid K. Ali, Mohamed A. Abd El Salam, Mohamed S. Mohamed . Chebyshev fifth-kind series approximation for generalized space fractional partial differential equations. AIMS Mathematics, 2022, 7(5): 7759-7780. doi: 10.3934/math.2022436
    [6] K. Ali Khalid, Aiman Mukheimer, A. Younis Jihad, Mohamed A. Abd El Salam, Hassen Aydi . Spectral collocation approach with shifted Chebyshev sixth-kind series approximation for generalized space fractional partial differential equations. AIMS Mathematics, 2022, 7(5): 8622-8644. doi: 10.3934/math.2022482
    [7] Yumei Chen, Jiajie Zhang, Chao Pan . Numerical approximation of a variable-order time fractional advection-reaction-diffusion model via shifted Gegenbauer polynomials. AIMS Mathematics, 2022, 7(8): 15612-15632. doi: 10.3934/math.2022855
    [8] Ahmed F. S. Aboubakr, Gamal M. Ismail, Mohamed M. Khader, Mahmoud A. E. Abdelrahman, Ahmed M. T. AbdEl-Bar, Mohamed Adel . Derivation of an approximate formula of the Rabotnov fractional-exponential kernel fractional derivative and applied for numerically solving the blood ethanol concentration system. AIMS Mathematics, 2023, 8(12): 30704-30716. doi: 10.3934/math.20231569
    [9] Dumitru Baleanu, Babak Shiri . Generalized fractional differential equations for past dynamic. AIMS Mathematics, 2022, 7(8): 14394-14418. doi: 10.3934/math.2022793
    [10] Chuanli Wang, Biyun Chen . An hp-version spectral collocation method for fractional Volterra integro-differential equations with weakly singular kernels. AIMS Mathematics, 2023, 8(8): 19816-19841. doi: 10.3934/math.20231010
  • This paper analyzes the coupled system of nonlinear fractional differential equations involving the caputo fractional derivatives of order α(1,2) on the interval (0, T). Our method of analysis is based on the reduction of the given system to an equivalent system of integral equations, then the resulting equation is discretized by using a spectral method based on the Legendre polynomials. We have constructed a Legendre spectral collocation method for the coupled system of nonlinear fractional differential equations. The error bounds under the L2 and Lnorms is also provided, then the theoretical result is validated by a number of numerical tests.



    Fractional-order derivatives arise in many applied fields, because the nonlocality for fractional calculus operator is very suitable for describing materials with memory and genetic properties. Such as control theory in power system [1,2], viscoelastic materials [3,4,5], information theory [6], electrical properties of materialsand [7], abnormal diffusion of ions in nerve cells [8,9,10], the modeling and analysis of various problems in bio-mathematical sciences [11,12], and etc.

    Inspired by the great popularity of the subject, many researchers turned to the further development of this branch of mathematical analysis. Coupled boundary conditions arise in the study of reaction-diffusion equations, Sturm-Liouville problems, mathematical biology [13,14,15] and so on; Many scholars have analyzed the existence of solutions of boundary value problems for coupled systems of nonlinear fractional differential equations [16,17,18,19,20,21,22], many of these are discussions of fractional α(1,2), but there are few numerical solutions.

    In general, there exists no method that yields an exact solution for the coupled system of nonlinear fractional differential equations. Only approximate solutions can be derived using linearization or perturbation methods. Z. Odibat et al. [23] presented He's homotopy perturbation method; H. Jafari et al. [24] employed Adomians decomposition method to give approximate solutions; S. Momani et al. [25] implemented variation iteration method to obtain approximate solutions; Zhou et al. [26] constructed a high order schemes for the numerical solution.

    These methods Most of the existing methods are discussed on α(0,1) and rarely on α(1,2). Moreover, Either these methods are based on local operations, and the effect of these methods is not good for nonlocality and weak singularity problems, or the convergence region of the corresponding results is rather small. Therefore, we construct a numerical method to solve them For the coupled fractional differential equations of α(1,2), this method converges rapidly when the solution is smooth, and still has good convergence when the solution is weakly singular.

    A general technique to construct legendre spectral collocation method for the numerical solution of the nonlinear fractional boundary value problems has been presented in [27]. In this paper we will extend this legendre spectral collocation method that is mentioned to the coupled system of nonlinear fractional differential equations. Due to the influence of coupling system and nonlinear term, the convergence analysis of spectral collocation method becomes very difficult. For this purpose, we use two kinds of polynomial interpolation, namely Legendre-Gauss interpolation and Jacobian-Gauss interpolation.

    The outline of the paper is as follows: In Sect. 2, We introduce some definitions and lemmas that will be used later. In Sect. 3, we transform the Eq (3.1) into an equivalent Volterra Fredholm integral equations, and replace the equations with a variable to get an equations defined on the interval (1,1). Then, we obtain a numerical scheme for problem (3.6) by using Legendre spectral collocation method. In Sect. 4, we derive the error analysis of the numerical scheme (3.9). In Sect. 5, some numerical experiments are provided to support the theoretical statement. Finally, some concluding remarks are given in the final section.

    Definition 2.1.(cf. p. 70 [28]) Let t(0,1), the left-sided Caputo derivative of order α, n1<α<n, nN+, are defined as:

    C0Dαtu(t)=1Γ(nα)t0(tτ)nα1u(n)(τ)dτ,

    where Γ() denotes Gamma function.

    Definition 2.2.(cf. [29]) Let Jα,βn(x),xΛ be the standard Jacobi polynomial of degree n. The set of Jacobi polynomials is a complete L2ωα,β(Λ)-orthogonal system, i.e.

    11Jα,βm(x)Jα,βn(x)ωα,β(x)dx=γα,βnδm,n, (2.1)

    where δm,n is the Kronecker function, and

    ωα,β(x)=(1x)α(1+x)β,forα,β>1.
    γα,βn={2α+β+1Γ(α+1)Γ(β+1)Γ(α+β+2),n=0,2α+β+1(2n+α+β+1)Γ(n+α+1)Γ(n+β+1)n!Γ(n+α+β+1),n1.

    Specially,

    Jα,β0(x)=1,Jα,β1(x)=12(α+β+2)x+12(αβ).

    Definition 2.3.(cf. [27]) For any of the Gauss-type quadratures defined above with the points and weights {xα,βn,ωα,βn}Nn=0(N0), we can define a discrete inner product in interval Λ :

    Λϕ(x)ωα,β(x)dxNn=0ϕ(xα,βn)ωα,βn. (2.2)

    Lemma 2.1.(cf. [29]) Let PN be the space of all polynomials of degree at most N, which is exact for any ϕ(x)P2N+1. Particularly,

    Nn=0Jα,βp(xα,βn)Jα,βq(xα,βn)ωα,βn=γα,βpδp,q,0p+q2N+1. (2.3)

    Lemma 2.2.(Lemma 2.1. [27]) For boundary value problem C0Dαty(t)=f(t,y(t)),t(0,T) with y(0)=y(T)=0, Let 1<α<2. Assume that y(t) is a function with an absolutely continuous first derivative, and f:[0,T]×RR is continuous. Then we have that yC1[0,T] is a solution of the problem if and only if it is a solution of the Fredholm integral equation:

    y(t)=1Γ(α)t0(tτ)α1f(τ,y(τ))dτtTΓ(α)T0(Tτ)α1f(τ,y(τ))dτ.

    We consider the coupled system of nonlinear fractional differential equations:

    {C0Dα1ty1(t)=f1(t,y1(t),y2(t)),C0Dα2ty2(t)=f2(t,y1(t),y2(t)),t(0,T),y1(0)=y1(T)=0,y2(0)=y2(T)=0, (3.1)

    where f1,f2:[0,T]×RR is continuous, and C0Dα1t, C0Dα2t is the left-sided Caputo derivative of order α1,α2 (1, 2).

    By Lemma 2.2, it has been proved easily that the problem (3.1) is equivalent to the following Fredholm integral equations when αi,yi(t),fi satisfy the condition of Lemma 2.2:

    yi(t)=1Γ(αi)t0(tτ)αi1fi(τ,y1(τ),y2(τ))dτtTΓ(αi)T0(Tτ)αi1fi(τ,y1(τ),y2(τ))dτ, (3.2)

    for t[0,T],i=1,2.

    Lemma 3.1. Let αi,yi(t),fi satisfy the condition of Lemma 2.2. Futhermore, let fi satisfy a Lipschitz condition with the Lipschitz constant L<Γ(α+1)4Tα. Γ(α+1)Tα=max(Γ(α1+1)Tα1,Γ(α2+1)Tα2). Then the system (3.1) has a unique solution.

    Proof. Let y(t)=(y1(t),y2(t)), define the operator Ay(t)=(A1y(t),A2y(t)). where

    Aiy(t)=1Γ(αi)t0(tτ)αi1fi(τ,y(τ))dτtTΓ(αi)T0(Tτ)αi1fi(τ,y(τ))dτ.

    Similar to Lemma 2.3 of [27], we have the following formula

    |Aiy(t)Aiˆy(t)|2LTαiΓ(αi+1)

    Let \|y-\widehat{y}\|_{L^{\infty}(0, T)} = \max(\|y_1-\widehat{y_1}\|, \|y_2-\widehat{y_2}\|) , we get

    \begin{eqnarray*} |\mathcal{A}_iy_i(t)-\mathcal{A}_i\widehat{y}_i(t)| \leq \cfrac{4LT^{\alpha}}{\Gamma(\alpha+1)}\|y-\widehat{y}\|_{L^{\infty}(0,T)}. \end{eqnarray*}

    We know C^1[0, T] \subset L^{\infty}(0, T) , this means that \mathcal{A}_i and \mathcal{A} is contraction. Then, by Banach's fixed point theorem, we know \mathcal{A} is a unique fixed point.

    Let \Lambda = [-1, 1] , then use the change of variable t = \cfrac{1}{2}T(x+1), x \in \Lambda . We transfer the problem (3.2) to an equivalent problem defined in \Lambda , then arrive at the following schema:

    \begin{eqnarray} y_i(\frac{1}{2}T(x+1))& = &\frac{1}{\Gamma(\alpha_i)}\int_0^{\frac{1}{2}T(x+1)}{ (\frac{1}{2}T(x+1)-\tau)^{\alpha_i-1} f_i(\tau,y_1(\tau),y_2(\tau))}d\tau\\ &&-\cfrac{x+1}{2\Gamma(\alpha_i)}\int_0^{T}(T-\hat{\tau})^{\alpha_i-1} f_i(\hat{\tau},y_1(\hat{\tau}),y_2(\hat{\tau}))d\hat{\tau}, \end{eqnarray} (3.3)

    where \tau = \cfrac{1}{2}T(\xi+1) and \hat{\tau} = \cfrac{1}{2}T(\lambda+1) . Furthermore, transfer the interval (0, \cfrac{1}{2}T(x+1)) to (-1, x) and (0, T) to (-1, 1) .

    For the convenience, let

    \begin{eqnarray*} &&Y_i(x) = y_i(\frac{1}{2}T(x+1)),\\ &&F_i(\xi, Y_1(\xi), Y_2(\xi)) = f_i(\frac{1}{2}T(\xi+1), y_1(\frac{1}{2}T(\xi+1)), y_2(\frac{1}{2}T(\xi+1))). \end{eqnarray*}

    Using the abbreviation, (3.3) can be read to

    \begin{eqnarray} Y_i(x)& = &\cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i)}\int_{-1}^{x}(x-\xi)^{\alpha_i-1} F_i(\xi,Y_1(\xi),Y_2(\xi))d\xi \\ &&-\cfrac{T^{\alpha_i}(x+1)}{2^{\alpha_{i}+1}\Gamma(\alpha_i)}\int_{-1}^1(1-\lambda)^{\alpha_i-1} F_i(\lambda,Y_1(\lambda),Y_2(\lambda))d\lambda. \end{eqnarray} (3.4)

    Finally, under the linear transformation

    \begin{eqnarray} \xi = \xi(x,\theta) : = \cfrac{x+1}{2}\theta + \cfrac{x-1}{2}, \; \theta \in \Lambda. \end{eqnarray} (3.5)

    In summary, we obtain

    \begin{eqnarray} Y_i(x)& = &\cfrac{T^{\alpha_i}(x+1)^{\alpha_i}}{4^{\alpha_i}\Gamma(\alpha_i)}\int_{-1}^1(1-\theta)^{\alpha_i-1} F_i(\xi(x,\theta),Y_1(\xi(x,\theta)),Y_2(\xi(x,\theta)))d\theta\\ &&-\cfrac{T^{\alpha_i}(x+1)}{2^{\alpha_{i}+1}\Gamma(\alpha_i)}\int_{-1}^1(1-\lambda)^{\alpha_i-1} F_i(\lambda,Y_1(\lambda),Y_2(\lambda))d\lambda. \end{eqnarray} (3.6)

    In the following, we will give a Legendre spectral collocation method for solving the system (3.6).

    For \forall v \in C(\Lambda) , we denote the Jacobi-Gauss interpolation operator in the x -direction: \mathcal{I}_{x, N}^{\alpha, \beta}:C(\Lambda) \rightarrow \mathcal{P}_N , such that

    \begin{eqnarray} \mathcal{I}_{x,N}^{\alpha,\beta}v(x_n^{\alpha,\beta}) = v(x_n^{\alpha,\beta}),\; 0 \leq n \leq N. \end{eqnarray} (3.7)

    Obviously

    \begin{eqnarray} \mathcal{I}_{x,N}^{\alpha,\beta}v(x) = \sum\limits_{p = 0}^Nv_p^{\alpha,\beta}\mathcal{J}_p^{\alpha,\beta}(x), \; \; where\; \; v_p^{\alpha,\beta} = \cfrac{1}{\gamma_p^{\alpha,\beta}}\sum\limits_{n = 0}^Nv(x_n)\mathcal{J}_p^{\alpha,\beta}(x_n)\omega_n^{\alpha,\beta}. \end{eqnarray} (3.8)

    When \alpha = \beta = 0 , the Jacobi polynomial is equivalent to the Legendre polynomial L_k(x) . Moreover, we read x_n = x_n^{0, 0}, \omega_n = \omega_n^{0, 0} and \mathcal{I}_{x, N} = \mathcal{I}_{x, N}^{0, 0} .

    Then we construct the schema for the next steps. We want to derive U_i(x) \in \mathcal{P}_N(\Lambda) with N \geq 1 , such that

    \begin{eqnarray} U_i(x)& = &\cfrac{T^{\alpha_i}}{4^{\alpha_i}\Gamma(\alpha_i)}\mathcal{I}_{x,N}[(x+1)^{\alpha_i} \int_{-1}^1(1-\theta_i)^{\alpha_i-1}\mathcal{I}_{\theta_i,N}^{\alpha_i-1,0}F_i(\xi(x,\theta_i), U_1(\xi(x,\theta_i)),U_2(\xi(x,\theta_i)))d\theta_i]\\ &&-\cfrac{T^{\alpha_i}(x+1)}{2^{\alpha_{i}+1}\Gamma(\alpha_i)}\int_{-1}^1(1-\lambda_i)^{\alpha_i-1} \mathcal{I}_{\lambda_i,N}^{\alpha_i-1,0}F_i(\lambda_i,U_1(\lambda_i),U_2(\lambda_i))d\lambda_i, \end{eqnarray} (3.9)

    where x = x^{0, 0}, \theta_i = \theta^{\alpha_i-1, 0}, \lambda_i = \lambda^{\alpha_i-1, 0} .

    The above formula is an implicit format. (3.9) has unique solution if F_i satisfies the Lipschitz condition with the Lipschitz constant L < \cfrac{\Gamma(\alpha+1)}{4T^{\alpha}} .

    Remark 1. The proof are similar to Appendix in [27] (reference therein), it is easy to the proof, so we omit the details.

    Next, we want to derive an approximation of scheme (3.9). We set

    \begin{eqnarray} &&U_i(x) = \sum\limits_{p = 0}^Nu_{i,p}L_p(x),\\ &&\mathcal{I}_{x,N}\mathcal{I}_{\theta_i,N}^{\alpha_i-1,0}((x+1)^{\alpha_i}F_i(\xi(x,\theta_i), U_1(\xi(x,\theta_i),U_2(\xi(x,\theta_i)))\\ && = \sum\limits_{p = 0}^{N}\sum\limits_{p' = 0}^{N}d_{i,p,p'}L_p(x)\mathcal{J}_{p'}^{\alpha_i-1,0}(\theta_i). \end{eqnarray} (3.10)

    Using (3.10) and (2.1), we arrive at the following schema:

    \begin{eqnarray} &&\cfrac{T^{\alpha_i}}{4^{\alpha_i}\Gamma(\alpha_i)}\int_{-1}^1(1-\theta_i)^{\alpha_i-1} \mathcal{I}_{x,N}\mathcal{I}_{\theta_i,N}^{\alpha_i-1,0}((x+1)^{\alpha_i}F_i(\xi(x,\theta_i), U_1(\xi(x,\theta_i)),U_2(\xi(x,\theta_i)))d\theta_i\\ && = \cfrac{T^{\alpha_i}}{4^{\alpha_i}\Gamma(\alpha_i)}\sum\limits_{p = 0}^{N}\sum\limits_{p' = 0}^{N}d_{i,p,p'}L_p(x) \int_{-1}^1(1-\theta_i)^{\alpha_i-1}\mathcal{J}_{p'}^{\alpha_i-1,0}(\theta_i)d\theta_i \\ && = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}\sum\limits_{p = 0}^{N}d_{i,p,0}L_p(x), \end{eqnarray} (3.11)

    where

    \begin{eqnarray} d_{i,p,0}& = &\cfrac{\alpha_i(2p+1)}{2^{1+\alpha_i}}\sum\limits_{m = 0}^{N}\sum\limits_{n = 0}^{N}(x_{m}+1)^{\alpha_i} F_i(\xi(x_{m},\theta_{n}^{\alpha_i-1,0}),\\ &&U_1(\xi(x_{m},\theta_{n}^{\alpha_i-1,0})), U_2(\xi(x_{m},\theta_{n}^{\alpha_i-1,0})))L_p(x_{m})\omega_{m}\omega_{n}^{\alpha_i-1,0}. \end{eqnarray} (3.12)

    Futhermore, by (2.2) we obtain

    \begin{eqnarray} &&\int_{-1}^1(1-\lambda_i)^{\alpha_i-1}\mathcal{I}_{\lambda_i,N}^{\alpha_i-1,0} F_i(\lambda_i,U_1(\lambda_i),U_2(\lambda_i))d\lambda_i \\ && = \sum\limits_{n = 0}^{N}F_i(\lambda_{n}^{\alpha_i-1,0},U_1(\lambda_{n}^{\alpha_i-1,0}),U_2(\lambda_{n}^{\alpha_i-1,0}))\omega_{n}^{\alpha_i-1,0}. \end{eqnarray} (3.13)

    To summarize, that is by combining (3.9)-(3.13), we arrive at the following overall schema

    \begin{eqnarray} \sum\limits_{p = 0}^Nu_{i,p}L_p(x) & = & \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}\sum\limits_{p = 0}^{N}d_{i,p,0}L_p(x) - \cfrac{T^{\alpha_i}(x+1)}{2^{\alpha_{i}+1}\Gamma(\alpha_i)} \sum\limits_{n = 0}^{N}F_i(\lambda_{n}^{\alpha_i-1,0},\\ &&U_1(\lambda_{n}^{\alpha_i-1,0}),U_2(\lambda_{n}^{\alpha_i-1,0}))\omega_{n}^{\alpha_i-1,0}. \end{eqnarray} (3.14)

    The coefficients of (3.14) yields are expanded and compared, we have

    \begin{eqnarray} \begin{cases} u_{i,p} = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,p,0}- \cfrac{T^{\alpha_i}}{2^{\alpha_{i}+1}\Gamma(\alpha_i)}\sum\limits_{n = 0}^{N}F_i(\lambda_{n}^{\alpha_i-1,0},\\ \; \; \; \; \; \; \; \; \; U_1(\lambda_{n}^{\alpha_i-1,0}),U_2(\lambda_{n}^{\alpha_i-1,0}))\omega_{n}^{\alpha_i-1,0},\; \; \; for \; p = 0,1,\\ u_{i,p} = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,p,0} \; \; \; for \; 2 \leq p \leq N. \end{cases} \end{eqnarray} (3.15)

    Remark 2. For the sake of economy of exposition, we focus on the coupled system ^C_0D^{\alpha_{i}}_{t}y_i(t) = f_i(t, y_1(t), y_2(t)), t \in (0, T) with y_i(0) = y_i(T) = 0 in which i = 1, 2 . In fact, Similar to the provied method in this paper, we can easily get the numerical scheme and error analysis of the system ^C_0D^{\alpha_{i}}_{t}y_i(t) = f_i(t, y_1(t), y_2(t), \cdots, y_z(t)), t \in (0, T) with y_i(0) = y_i(T) = 0 in which i = 1, 2, \cdots, z . In Example 5.3, we also give the numerical experiment when z = 3 .

    Remark 3. If y_i(0)\neq 0, y_i(T)\neq 0 , That is, the initial value problem, with the initial value is not zero.

    \begin{eqnarray} \left\{\begin{array}{ll} ^C_0D^{\alpha_{1}}_{t}y_1(t) = f_1(t,y_1(t),y_2(t)),\\ ^C_0D^{\alpha_{2}}_{t}y_2(t) = f_2(t,y_1(t),y_2(t)),&t \in (0,T),\\ y_1^{(k)}(0) = c_1^k,y_2^{(k)}(0) = c_2^k. \end{array}\right. \end{eqnarray} (3.16)

    By Volterra integral equations the problem (3.16) is equivalent to the following Fredholm integral equations

    \begin{align} \left\{\begin{array}{ll} y_1(t) = g_1(t)+\cfrac{1}{\Gamma(\alpha_1)}\int_0^{t}(t-\tau)^{\alpha_1-1}f_1(\tau, y_1(\tau), y_2(\tau))d\tau, \\ y_2(t) = g_2(t)+\cfrac{1}{\Gamma(\alpha_1)}\int_0^{t}(t-\tau)^{\alpha_1-1}f_2(\tau, y_1(\tau), y_2(\tau))d\tau, \end{array}\right. \end{align} (3.17)

    where g_i(t) = \sum\limits_{k = 0}^{n_i-1}c_i^k\cfrac{t^k}{k!}, i = 1, 2, \ldots, n.

    Similarly, we can get

    \begin{eqnarray} \begin{cases} u_{i,p} = c_i^{0}+c_i^{1}\cfrac{T}{2}+\cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,p,0},\; \; \; for \; p = 0,\\ u_{i,p} = c_i^{1}\cfrac{T}{2}+\cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,p,0},\; \; \; for \; p = 1,\\ u_{i,p} = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,p,0} \; \; \; for \; 2 \leq p \leq N. \end{cases} \end{eqnarray} (3.18)

    In the next sections, we will give a error analysis for the above schema under the space L^2(\Lambda) and L^{\infty}(\Lambda) , respectively.

    We introduce the Jacobi-weighted Sobolev space and the norm and semi-norm.

    \begin{eqnarray*} H^l_{\omega^{\alpha,\beta}}(\Lambda) = \{v:\|v\|_{H^l_{\omega^{\alpha,\beta}}} < \infty\},\; \; l \geq 0, \end{eqnarray*}
    \begin{eqnarray*} \|v\|_{H^l_{\omega^{\alpha,\beta}}} = (\sum\limits_{k = 0}^l|v|_{H^l_{\omega^{\alpha,\beta}}})^{\frac{1}{2}},\; \; \; |v|_{H^k_{\omega^{\alpha,\beta}}} = \|\partial_x^kv\|_{\omega^{\alpha+k,\beta+k}}, \end{eqnarray*}

    where \|\cdot\|_{\omega^{\alpha, \beta}} is the weighted L^2_{\omega^{\alpha, \beta}}(\Lambda) -norm. Especially, L^2(\Lambda) = H^0_{\omega^{0, 0}}(\Lambda) , \|\cdot\| = \|\cdot\|_{L^2(\Lambda)} and \|\cdot\|_{\infty} = \|\cdot\|_{L^{\infty}(\Lambda)} .

    Now, we analyze the errors of numerical scheme (3.9). Let e_i(x) = Y_i(x)-U_i(x), (i = 1, 2) and using \mathcal{I} to represent identity operator. Observely,

    \begin{eqnarray} &&\|e_i\| \leq \|Y_i -\mathcal{I}_{x,N}Y_i\| + \|\mathcal{I}_{x,N}Y_i - U_i\|. \end{eqnarray} (4.1)

    Moreover, let e(x) = {\bf{Y}}(x)-{\bf{U}}(x) . where {\bf{Y}}(x) = (Y_1(x), Y_2(x)) , {\bf{U}}(x) = (U_1(x), U_2(x)) .

    Lemma 4.1. For N \geq 1 , We can get the following inequality

    \begin{eqnarray*} \|e_i\| \leq \sum\limits_{n = 1}^{5}\|B_{i,n}\|, \; \; for\; i = 1, 2, \end{eqnarray*}

    where

    \begin{eqnarray*} B_{i,1}(x) & = & Y_i(x)-\mathcal{I}_{x,N}Y_i(x), \nonumber\\ B_{i,2}(x) & = & \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i)}\mathcal{I}_{x,N} \int_{-1}^{x}(x-\xi_i)^{\alpha_i-1}(\mathcal{I}-_{x}\widetilde{\mathcal{I}}_{\xi_i,N}^{\alpha_i-1,0}) F_i(\xi_i,{\bf{Y}}(\xi_i))d\xi_i, \nonumber\\ B_{i,3}(x) & = & \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i)}\mathcal{I}_{x,N} \int_{-1}^{x}(x-\xi_i)^{\alpha_i-1}{_{x}\widetilde{\mathcal{I}}_{\xi_i,N}^{\alpha_i-1,0}} (F_i(\xi_i,{\bf{Y}}(\xi_i))-F_i(\xi_i,{\bf{U}}(\xi_i)))d\xi_i, \nonumber\\ B_{i,4}(x) & = & \cfrac{T^{\alpha_i}(x+1)}{2^{\alpha_i+1}\Gamma(\alpha_i)}\int_{-1}^{1}(1-\lambda_i)^{\alpha_i-1} \mathcal{I}_{\lambda_i,N}^{\alpha_i-1,0} (F_i(\lambda_i,{\bf{U}}(\lambda_i))-F_i(\lambda_i,{\bf{Y}}(\lambda_i)))d\lambda_i, \nonumber\\ B_{i,5}(x) & = & \cfrac{T^{\alpha_i}(x+1)}{2^{\alpha_i+1}\Gamma(\alpha_i)}\int_{-1}^{1}(1-\lambda_i)^{\alpha_i-1} (\mathcal{I}_{\lambda_i,N}^{\alpha_i-1,0}-\mathcal{I}) F_i(\lambda_i,{\bf{Y}}(\lambda_i))d\lambda_i. \end{eqnarray*}

    Proof. Similar to Lemma 4.3 of [27], we can easily prove it.

    Lemma 4.2. Assume the hypothesis of Lemma 4.1. We have the following inequality

    \begin{eqnarray*} \|e\| \leq \sum\limits_{n = 1}^{5}\|B_{1,n}\| + \sum\limits_{n = 1}^{5}\|B_{2,n}\|. \end{eqnarray*}

    Proof. For e(x) = {\bf{Y}}(x)-{\bf{U}}(x) . we have

    \begin{eqnarray*} e(x)& = &{\bf{Y}}(x)-{\bf{U}}(x)\\ & = & (Y_1(x), Y_2(x))-(U_1(x), U_2(x)) \\ & = & (Y_1(x)-U_1(x), Y_2(x)-U_2(x))\\ & = & (e_1(x), e_2(x)). \end{eqnarray*}

    Then, combination with Lemma 4.1

    \begin{eqnarray*} \|e\| & = & \|(e_1, e_2)\| = \sqrt{\|e_1\|^2 + \|e_2\|^2} \leq \|e_1\| + \|e_2\| \\ &\leq& \sum\limits_{n = 1}^{5}\|B_{1,n}\| + \sum\limits_{n = 1}^{5}\|B_{2,n}\|. \end{eqnarray*}

    We set the Nemytskii operator

    \begin{eqnarray*} \mathbb{F}:H_{\omega^{l,l}}^{l}(\Lambda) \rightarrow H_{\omega^{\alpha_i+l-1,l}}^{l}(\Lambda) (l \in {\bf{N}}, 1 \leq l \leq N+1). \end{eqnarray*}

    Let {\bf{Y}}(x) as the solution of system (3.6) and {\bf{U}}(x) as the solution of system (3.9). We define that

    \begin{eqnarray*} \mathbb{F}({\bf{Y}})(x): = F(x,{\bf{Y}}(x)), \end{eqnarray*}
    \begin{eqnarray*} \|\partial_{x}^{l}Y\|_{\omega^{l,l}} = \max(\|\partial_{x}^{l}Y_1\|_{\omega^{l,l}}, \|\partial_{x}^{l}Y_2\|_{\omega^{l,l}}), \end{eqnarray*}
    \begin{eqnarray*} \|\partial_{x}^{l}F(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}} = \max(\|\partial_{x}^{l}F_1(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_1+l-1,l}}, \|\partial_{x}^{l}F_2(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_2+l-1,l}}). \end{eqnarray*}

    Now we prove the convergence of the spectral collocation method in the space L^2(\Lambda) .

    Lemma 4.3. Suppose that \alpha_i \in (1, 2) , Y_i \in H_{\omega^{l, l}}^{l}(\Lambda), i = 1, 2 , F_i fulfills the Lipschitz condition with the Lipschitz constant L < \cfrac{\Gamma(\alpha+1)}{4T^{\alpha}} . Then we can obtain

    \begin{eqnarray*} \|Y_i-U_i\| \leq cN^{-l}(\|\partial_{x}^{l}Y\|_{\omega^{l,l}}+\|\partial_{x}^{l} F(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_i+l-1,l}}). \end{eqnarray*}

    Proof. Similar to Lemma 4.4 of [27], we can easily get

    \begin{eqnarray} &&\|B_{i,1}\| = \|Y_i - \mathcal{I}_{x,N}Y_i\| \leq cN^{-l}\|\partial_{x}^{l}Y_i\|_{\omega^{l,l}},\\ &&\|B_{i,2}\| \leq cN^{-l}\|\partial_{\xi}^lF_i(\cdot,{\bf{Y}}(\cdot)\|_{\omega^{\alpha_i+l-1,l}},\\ &&\|B_{i,3}\| \leq \frac{1}{2}(cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_i+l-1,l}}+\sqrt{\frac{\alpha_i}{3\times2^{\alpha_i}}}\|e\|),\\ &&\|B_{i,4}\| \leq \frac{1}{2}(\sqrt{\frac{\alpha_i}{12}}\|e\|+cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_i+l-1,l}}),\\ &&\|B_{i,5}\| \leq cN^{-l}\|\partial_{x}^lF_i(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_i+l-1,l}}. \end{eqnarray} (4.2)

    Then by Lemma 4.1, we get

    \begin{eqnarray} \|e_i\| &\leq& cN^{-l}\|\partial_{x}^{l}Y_i\|_{\omega^{l,l}}+ cN^{-l}\|\partial_{\xi}^lF_i(\cdot,{\bf{Y}}(\cdot)\|_{\omega^{\alpha_i+l-1,l}}\\ &&+\frac{1}{2}(cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_i+l-1,l}}+\sqrt{\frac{\alpha_i}{3\times2^{\alpha_i}}}\|e\|)\\ &&+\frac{1}{2}(\sqrt{\frac{\alpha_i}{12}}\|e\|+cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_i+l-1,l}})\\ &&+ cN^{-l}\|\partial_{x}^lF_i(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_i+l-1,l}}\\ &\leq& cN^{-l}(\|\partial_{x}^{l}Y\|_{\omega^{l,l}} +\|\partial_{x}^lF(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}})\\ &&+\frac{1}{2}(\sqrt{\frac{\alpha_i}{3\times 2^{\alpha_i}}}+\sqrt{\frac{\alpha_i}{12}})\|e\|. \end{eqnarray} (4.3)

    Then, by (4.3) and Lemma 4.2, we have

    \begin{eqnarray*} \|{\bf{Y}}-{\bf{U}}\|&\leq& cN^{-l}(\|\partial_{x}^{l}Y\|_{\omega^{l,l}} +\|\partial_{x}^lF(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}}) \\ &&+\frac{1}{2}(\sqrt{\frac{\alpha_1}{3\times 2^{\alpha_1}}}+\sqrt{\frac{\alpha_1}{12}})\|e\| +\frac{1}{2}(\sqrt{\frac{\alpha_2}{3\times 2^{\alpha_2}}}+\sqrt{\frac{\alpha_2}{12}})\|e\|. \end{eqnarray*}

    Clearly,

    \begin{eqnarray*} \sqrt{\frac{\alpha_i}{3\times 2^{\alpha_i}}}+\sqrt{\frac{\alpha_i}{12}} < 1,\; \; \; \; \; \forall \alpha_i \in (1,2). \end{eqnarray*}

    Then, we easily get

    \begin{eqnarray} \|{\bf{Y}}-{\bf{U}}\|&\leq& cN^{-l}(\|\partial_{x}^{l}Y\|_{\omega^{l,l}} +\|\partial_{x}^lF(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}}). \end{eqnarray} (4.4)

    Thus

    \begin{eqnarray*} \|Y_i-U_i\| &\leq& cN^{-l}(\|\partial_{x}^{l}Y\|_{\omega^{l,l}} +\|\partial_{x}^lF(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}}). \end{eqnarray*}

    Let u_i(t): = U_i(\frac{2t}{T}-1)(i = 1, 2) be the approximate solution obtained by using the Legendre spectral collocation method with t \in (0, T) , \chi^{\alpha, \beta}(t): = (T-t)^{\alpha}t^{\beta} is defined as a weighting function. Furthermore, define the Nemytskii operator

    \begin{eqnarray*} \mathbb{K}:H_{\chi^{l,l}}^{l}(0,T)\rightarrow H_{\chi^{\alpha_i+l-1,l}}^{l}(0,T)(l \in {\bf{N}}, 1 \leq l \leq N+1). \end{eqnarray*}

    Let {\bf{y}}(t) = (y_1(t), y_2(t)) as the exact solution to the system (3.1), and {\bf{u}}(t) = (u_1(t), u_2(t)) be the approximate solution. We define that

    \begin{eqnarray*} \mathbb{K}(y_i)(t): = f_i(t,y_1(t),y_2(t)),i = 1,2, \end{eqnarray*}
    \begin{eqnarray*} \|\partial_{t}^{l}y\|_{L_{\chi^{l,l}}^{2}(0,T)} = \max(\|\partial_{t}^{l}y_1\|_{L_{\chi^{l,l}}^{2}(0,T)}, \|\partial_{t}^{l}y_2\|_{L_{\chi^{l,l}}^{2}(0,T)}), \end{eqnarray*}
    \begin{eqnarray*} \|\partial_{t}^{l}f(\cdot,{\bf{y}}(\cdot))\|_{L_{\chi^{\alpha+l-1,l}}^{2}(0,T)} = \max(\|\partial_{t}^{l}f_1(\cdot,{\bf{y}}(\cdot))\|_{L_{\chi^{\alpha_1+l-1,l}}^{2}(0,T)}, \|\partial_{t}^{l}f_2(\cdot,{\bf{y}}(\cdot))\|_{L_{\chi^{\alpha_2+l-1,l}}^{2}(0,T)}). \end{eqnarray*}

    Then, according to the above Lemma, we can get the following theorem.

    Theorem 4.1. Suppose that \alpha_i \in (1, 2) , y_i \in H_{\chi^{l, l}}^{l}(0, T) . Moreover, f_i fulfills the Lipschitz condition with the Lipschitz constant L < \cfrac{\Gamma(\alpha_i+1)}{4T^{\alpha_i}} . Then we get

    \begin{eqnarray*} \|y_i-u_i\| \leq cN^{-l}(\|\partial_{t}^{l}y\|_{L_{\chi^{l,l}}^{2}(0,T)}+\|\partial_{t}^{l} f(\cdot,{\bf{y}}(\cdot))\|_{L_{\chi^{\alpha+l-1,l}}^{2}(0,T)}). \end{eqnarray*}

    Next, we estimate the error in function space L^{\infty}(\Lambda) .

    Lemma 4.4. Suppose that \alpha_i \in (1, 2) , Y_i \in L^{\infty}(\Lambda)\bigcap H^{l}(\Lambda) , F_i fulfills the Lipschitz condition with the Lipschitz constant L < \cfrac{\Gamma(\alpha+1)}{4T^{\alpha}} . Then we get

    \begin{eqnarray*} \|Y_i-U_i\| \leq cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y\|_{\omega^{l,l}} +cN^{\frac{1}{2}-l}\|\partial_{x}^{l}F(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_i+l-1,l}},i = 1,2. \end{eqnarray*}

    Proof. Similarly, we get that

    \begin{eqnarray} \|e_i\|_{\infty} \leq \|Y_i-\mathcal{I}_{x,N}Y_i\|_{\infty} + \|\mathcal{I}_{x,N}Y_i-U_i\|_{\infty} \leq \sum\limits_{n = 1}^5\|B_{i,n}\|_{\infty}. \end{eqnarray} (4.5)

    Then, similar to Lemma 4.5 of [27], we can easily get

    \begin{eqnarray} &&\|B_{i,1}\|_{\infty} \leq cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y_i\|,\\ &&\|B_{i,2}\|_{\infty} \leq cN^{\frac{1}{2}-l}\|\partial_{\xi}^{l}F_i(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_i+l-1,l}},\\ &&\|B_{i,3}\|_{\infty} \leq \frac{1}{2}(cN^{\frac{1}{2}-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_i+l-1,l}}+cN^{\frac{1}{2}}\|e\|),\\ &&|B_{i,4}|_{\infty} \leq \frac{1}{2}(\frac{1}{2}\|e\|_{\infty} + cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_i+l-1,l}}),\\ &&|B_{i,5}|_{\infty} \leq cN^{-l}\|\partial_{x}^{l}F_i(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_i+l-1,l}}. \end{eqnarray} (4.6)

    By (4.6) and Lemma 4.2, we have

    \begin{eqnarray*} \|{\bf{Y}}-{\bf{U}}\|_{\infty} &\leq& (cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y_1\|)+(cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y_2\|) \nonumber\\ &&+(cN^{\frac{1}{2}-l}\|\partial_{\xi}^{l}F_1(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_1+l-1,l}}) + (cN^{\frac{1}{2}-l}\|\partial_{\xi}^{l}F_2(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_2+l-1,l}})\nonumber\\ &&+\frac{1}{2}cN^{\frac{1}{2}-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_1+l-1,l}} +\frac{1}{2}cN^{\frac{1}{2}-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_2+l-1,l}}+cN^{\frac{1}{2}}\|e\| \nonumber\\ &&+\frac{1}{2}\|e\|_{\infty} + \frac{1}{2}cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_1+l-1,l}} + \frac{1}{2}cN^{-l}\|\partial_{x}^{l}{\bf{Y}}\|_{\omega^{\alpha_2+l-1,l}}\nonumber\\ &&+cN^{-l}\|\partial_{x}^{l}F_1(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_1+l-1,l}} +cN^{-l}\|\partial_{x}^{l}F_2(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha_2+l-1,l}}\nonumber\\ &\leq& cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y\| +cN^{\frac{1}{2}-l}\|\partial_{x}^lF(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}} +cN^{\frac{1}{2}}\|e\|+\frac{1}{2}\|e\|_{\infty}. \end{eqnarray*}

    Then, combination with (4.4) we can obtain

    \begin{eqnarray*} \|{\bf{Y}}-{\bf{U}}\|_{\infty} &\leq& cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y\| +cN^{\frac{1}{2}-l}\|\partial_{x}^lF(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}}. \end{eqnarray*}

    Thus

    \begin{eqnarray*} \|Y_i-U_i\| &\leq& cN^{\frac{3}{4}-l}\|\partial_{x}^{l}Y\| +cN^{\frac{1}{2}-l}\|\partial_{x}^lF_i(\cdot,{\bf{Y}}(\cdot))\|_{\omega^{\alpha+l-1,l}}. \end{eqnarray*}

    Then, according to the Lemma 4.5, we can get the following theorem.

    Theorem 4.2. Suppose that \alpha_i \in (1, 2) , y_i \in L^{\infty}(0, T)\bigcap H^{l}(0, T) . Futhermore, f_i fulfills the Lipschitz condition with the Lipschitz constant L < \cfrac{\Gamma(\alpha+1)}{4T^{\alpha}} . Then we have

    \begin{eqnarray*} \|y_i-u_i\|_{L^{\infty}(0,T)} \leq cN^{\frac{3}{4}-l}\|\partial_{t}^{l}y\|_{L^2(0,T)} +cN^{\frac{1}{2}-l}\|\partial_{t}^{l}f(\cdot,{\bf{y}}(\cdot))\|_{L_{\chi^{\alpha+l-1,l}}^2(0,T)}. \end{eqnarray*}

    In this section, we carry out some numerical experiments to verify the theoretical results of the previous method.

    Example 5.1. We consider the coupled system of nonlinear problems with weakly singular solutions as follows:

    \begin{eqnarray*} \begin{cases} ^C_0D^{\alpha_1}_{t}y_1(t) = y_1^2(t)+y_2^2(t)+g_1(t),\\ ^C_0D^{\alpha_2}_{t}y_2(t) = y_1^2(t)-y_2^2(t)+g_2(t),&t \in (0,1),\\ y_1(0) = y_1(1) = 0,y_2(0) = y_2(1) = 0, \end{cases} \end{eqnarray*}

    where g_1(t) = -\Gamma(2+\alpha_1)t-(t-t^{1+\alpha_1})^2-(t-t^{1+\frac{\alpha_2}{2}})^2 , g_2(t) = -\cfrac{\Gamma(2+\frac{\alpha_2}{2})}{\Gamma(2-\frac{\alpha_2}{2})}t^{1-\frac{\alpha_2}{2}} -(t-t^{1+\alpha_1})^2+(t-t^{1+\frac{\alpha_2}{2}})^2 . It can be verified that the exact solution is y_1(t) = t-t^{1+\alpha_1} , y_2(t) = t-t^{1+\frac{\alpha_2}{2}} , which is weakly singular at the endpoint t = 0.

    When N = 2 , by (3.10) we have

    \begin{eqnarray} U_i(x) = u_{i,0}L_0(x) + u_{i,1}L_1(x) + u_{i,2}L_2(x), \end{eqnarray} (5.1)

    where L_0(x), L_1(x), L_2(x) known, u_{i, 0}, u_{i, 1}, u_{i, 2} by (3.15) we can obtain

    \begin{eqnarray*} \begin{cases} u_{i,0} = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,0,0}-\cfrac{T^{\alpha_i}}{2^{\alpha_{i}+1}\Gamma(\alpha_i)} \sum\limits_{n = 0}^{N}F_i(\lambda_{n}^{\alpha_i-1,0},U_1(\lambda_{n}^{\alpha_i-1,0}),U_2(\lambda_{n}^{\alpha_i-1,0}))\omega_{n}^{\alpha_i-1,0},\\ u_{i,1} = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,1,0}-\cfrac{T^{\alpha_i}}{2^{\alpha_{i}+1}\Gamma(\alpha_i)} \sum\limits_{n = 0}^{N}F_i(\lambda_{n}^{\alpha_i-1,0},U_1(\lambda_{n}^{\alpha_i-1,0}),U_2(\lambda_{n}^{\alpha_i-1,0}))\omega_{n}^{\alpha_i-1,0},\\ u_{i,2} = \cfrac{T^{\alpha_i}}{2^{\alpha_i}\Gamma(\alpha_i+1)}d_{i,2,0}. \end{cases} \end{eqnarray*}

    d_{i, p, 0} is given in (3.12). Substitute the above results into (5.1) to obtain U_i(x) , by (3.2), (3.3), (3.6), (3.7), (3.9), we know

    \begin{eqnarray*} U_i(x) = Y_i(x) = y_i(\frac{1}{2}T(x+1)) = y_i(t). \end{eqnarray*}

    In Table 1, we show The value of u_{i, p} when \alpha_1 = \alpha_2 = 5/4 , N = 2 . In Table 2, we show substitute the results of tab 1 into (5.1) to obtain numerical solution and the exact solution. In Table 3, we show the error of y_1, y_2 varies with N when \alpha_1 = \alpha_2 = 5/4 . Moreover, In Figures 1 and 2, we show the error between the approximate solution and the exact solution in L^{\infty}(\Lambda) and L^2(\Lambda) space, respectively. It is observed that for all \alpha \in (1, 2) , the method has good convergence even though the solution is weakly singular.

    Table 1.  The value of u_{i, p} when \alpha_1 = \alpha_2 = 5/4 , N = 2 .
    i u_{i, 0} u_{i, 1} u_{i, 2}
    1 0.1923 0.0115 -0.1951
    2 0.1201 -0.0129 -0.1145

     | Show Table
    DownLoad: CSV
    Table 2.  Numerical solution and exact solution when \alpha_1 = \alpha_2 = 5/4 , N = 2 .
    x U_1(x)^{1} U_2(x)^{2} t ^{3} y_1(t) y_2(t)
    -1.0000 -0.0144 0.0186 0 0 0
    -0.9195 0.0318 0.0441 0.0402 0.0395 0.0348
    -0.7388 0.1216 0.0932 0.1306 0.1204 0.0940
    -0.4779 0.2175 0.1443 0.2610 0.2123 0.1483
    -0.1653 0.2800 0.1748 0.4174 0.2774 0.1756
    0.1653 0.2838 0.1706 0.5826 0.2861 0.1669
    0.4779 0.2285 0.1320 0.7390 0.2327 0.1273
    0.7388 0.1387 0.0741 0.8694 0.1395 0.0728
    0.9195 0.0530 0.0203 0.9598 0.0480 0.0243
    1.0000 0.0087 -0.0073 1.0000 0 0
    Note: ^{1} U_1(x)=\sum\limits_{p=0}^2u_{1, p}L_p(x). ^{2} U_2(x)=\sum\limits_{p=0}^2u_{2, p}L_p(x). ^{3} t=\cfrac{1}{2}T(x+1).

     | Show Table
    DownLoad: CSV
    Table 3.  Numerical solution and exact solution when \alpha_1 = \alpha_2 = 5/4 , N = 2 .
    N L^{\infty}-error of y_1 L^{2}-error of y_1 L^{\infty}-error of y_2 L^{2}-error of y_2
    2 1.4362e-02 4.3702e-03 1.8581e-02 4.3702e-03
    4 7.9055e-04 4.5786e-04 3.0695e-03 4.5786e-04
    6 1.6147e-04 1.3548e-04 1.0241e-03 1.3548e-04
    8 5.1592e-05 6.0834e-05 4.5738e-04 6.0834e-05
    10 2.1007e-05 3.3689e-05 2.4107e-04 3.3689e-05
    12 9.9849e-06 2.1058e-05 1.4155e-04 2.1058e-05
    14 5.2881e-06 1.4235e-05 8.9695e-05 1.4235e-05
    16 3.0344e-06 1.0160e-05 6.0167e-05 1.0160e-05
    18 1.8525e-06 7.5343e-06 4.2182e-05 7.5343e-06
    20 1.1881e-06 5.7290e-06 3.0635e-05 5.7290e-06

     | Show Table
    DownLoad: CSV
    Figure 1.  L^{\infty}-errors and L^{2}-errors of y_1 for Expample 5.1.
    Figure 2.  L^{\infty}-errors and L^{2}-errors of y_2 for Expample 5.1.

    Example 5.2. We consider the coupled system of nonlinear problems with smooth solutions as follows:

    \begin{eqnarray*} \begin{cases} ^C_0D^{\alpha_1}_{t}y_1(t) = y_1^2(t)+y_2^2(t)+g_1(t),\\ ^C_0D^{\alpha_2}_{t}y_2(t) = y_1^2(t)-y_2^2(t)+g_2(t),&t \in (0,1),\\ y_1(0) = y_1(1) = 0,y_2(0) = y_2(1) = 0, \end{cases} \end{eqnarray*}

    where g_1(t) = -\cfrac{\Gamma(17/4)}{\Gamma(17/4-\alpha_1)}t^{13/4-\alpha_1}-(t-t^{13/4})^2-(t-t^{15/4})^2 , g_2(t) = -\cfrac{\Gamma(19/4)}{\Gamma(19/4-\alpha_2)}t^{15/4-\alpha_2}-(t-t^{13/4})^2+(t-t^{15/4})^2 . It can be verified that the exact solution is y_1(t) = t-t^{\frac{13}{4}} , y_2(t) = t-t^{\frac{15}{4}} , which is smooth on the interval [0, 1].

    In Figures 3 and 4, we show the error between the approximate solution and the exact solution in L^{\infty}(\Lambda) and L^2(\Lambda) space, respectively. It is observed that for all \alpha \in (1, 2) , the method converges rapidly when the exact solution is very smooth.

    Figure 3.  L^{\infty}-errors and L^{2}-errors of y_1 for Expample 5.2.
    Figure 4.  L^{\infty}-errors and L^{2}-errors of y_2 for Expample 5.2.

    Example 5.3. We consider the system of nonlinear problems as follows:

    \begin{eqnarray*} \begin{cases} ^C_0D^{\alpha_1}_{t}y_1(t) = y_1^2(t)+y_2^2(t)+y_3^2(t)+g_1(t),\\ ^C_0D^{\alpha_2}_{t}y_2(t) = y_1^2(t)-y_2^2(t)+y_3^2(t)+g_2(t),\\ ^C_0D^{\alpha_3}_{t}y_3(t) = y_1^2(t)+y_2^2(t)-y_3^2(t)+g_3(t),&t \in (0,1),\\ y_1(0) = y_1(1) = 0,y_2(0) = y_2(1) = 0,y_3(0) = y_3(1) = 0, \end{cases} \end{eqnarray*}

    where g_1(t) = -\Gamma(2+\alpha_1)t-(t-t^{1+\alpha_1})^2-(t-t^{1+\frac{\alpha_2}{2}})^2-(t-t^{5/4})^2 , g_2(t) = -\cfrac{\Gamma(2+\frac{\alpha_2}{2})}{\Gamma(2-\frac{\alpha_2}{2})}t^{1-\frac{\alpha_2}{2}} -(t-t^{1+\alpha_1})^2+(t-t^{1+\frac{\alpha_2}{2}})^2-(t-t^{5/4})^2 , g_3(t) = -\cfrac{\Gamma(17/4)}{\Gamma(17/4-\alpha_3)}t^{13/4-\alpha_3}-(t-t^{1+\alpha_1})^2-(t-t^{1+\frac{\alpha_2}{2}})^2+(t-t^{13/4})^2 . It can be verified that the exact solution is y_1(t) = t-t^{1+\alpha_1} , y_2(t) = t-t^{1+\frac{\alpha_2}{2}} , y_3(t) = t-t^{13/4} .

    In Figures 57, we show the error between the approximate solution and the exact solution in L^{\infty}(\Lambda) and L^2(\Lambda) space, respectively. The results show that this method still has good convergence when z = 3 .

    Figure 5.  L^{\infty}-errors and L^{2}-errors of y_1 for Expample 5.3.
    Figure 6.  L^{\infty}-errors and L^{2}-errors of y_2 for Expample 5.3.
    Figure 7.  L^{\infty}-errors and L^{2}-errors of y_3 for Expample 5.3.

    Example 5.4. We consider the coupled system of nonlinear problems with initial value is not zero as follows:

    \begin{eqnarray*} \begin{cases} ^C_0D^{\alpha_1}_{t}y_1(t) = y_1^2(t)+y_2^2(t)+g_1(t),\\ ^C_0D^{\alpha_2}_{t}y_2(t) = y_1^2(t)-y_2^2(t)+g_2(t),&t \in (0,1),\\ y_1^{(k)}(0) = c_1^k,y_2^{(k)}(0) = c_2^k, \end{cases} \end{eqnarray*}

    where g_1(t) = \Gamma(2+\alpha_1)t-(t^{1+\alpha_1}-3t+2)^2-(t^{13/4}-5t+2)^2 , g_2(t) = \cfrac{\Gamma(17/4)}{\Gamma(17/4-\alpha_2)}t^{13/4-\alpha_2}-(t^{1+\alpha_1}-3t+2)^2+(t^{13/4}-5t+2)^2 . It can be verified that the exact solution is y_1(t) = t^{1+\alpha_1}-3t+2 , y_2(t) = t^{13/4}-5t+2 . y_1(0) = 2, y_1'(0) = -3, y_1(1) = 0, y_2(0) = 2, y_2'(0) = -5, y_2(1) = -2 .

    In Figures 8 and 9, we show the error between the approximate solution and the exact solution in L^{\infty}(\Lambda) and L^2(\Lambda) space, respectively, where y_i(0)\neq 0 . The results show that this method still has good convergence when the initial value is not zero.

    Figure 8.  L^{\infty}-errors and L^{2}-errors of y_1 for Expample 5.4.
    Figure 9.  L^{\infty}-errors and L^{2}-errors of y_2 for Expample 5.2.

    We presented a Legendre spectral collocation method for the system of nonlinear fractional differential equations. We established an error estimate for the numerical solution, and showing that the proposed schema is converges. The carried out numerical tests confirmed the theoretical prediction.

    The authors would like to thank the referees and editor for their valuable comments and suggestions which helped to improve the results of this paper. This work was supported by the Ph.D. Research-Starting Foundation of Guizhou Normal University, China ([2016] Numerical Method of Anomalous Diffusion and Its Related Problems)

    The authors declare that there are no conflicts of interest.



    [1] H. H. Sun, A. A. Abdelwahab, B. Onaral, Linear approximation of transfer function with a pole of fractional power, IEEE T. Automat. Contr., 29 (1984), 441–444. https://doi.org/10.1109/TAC.1984.1103551 doi: 10.1109/TAC.1984.1103551
    [2] A. Oustaloup, La dérivation non entière: Thorie, synthèse et applications, Hermes Science Publications, 1995.
    [3] R. C. Koeller, Applications of fractional calculus to the theory of viscoelasticity, J. Appl. Mech., 51 (1984), 299–307. https://doi.org/10.1115/1.3167616 doi: 10.1115/1.3167616
    [4] R. L. Bagley, R. A. Calico, Fractional order state equations for the control of viscoelasticallydamped structures, J. Guid. Control Dynam., 14 (1991), 304–311. https://doi.org/10.2514/3.20641 doi: 10.2514/3.20641
    [5] F. Amblard, A. C. Maggs, B. Yurke, A. N. Pargellis, S. Leibler, Subdiffusion and anomalous local viscoelasticity in actin networks, Phys. Rev. Lett., 77 (1996), 4470–4473. https://doi.org/10.1103/PhysRevLett.77.4470 doi: 10.1103/PhysRevLett.77.4470
    [6] B. Mandelbrot, Some noises with 1/f spectrum, a bridge between direct current and white noise, IEEE T. Inf. Theory., 13 (1967), 289–298. https://doi.org/10.1109/TIT.1967.1053992 doi: 10.1109/TIT.1967.1053992
    [7] G. E. Carlson, C, Halijak, Approximation of fractional capacitors by a regular newton process, IEEE T. Circuit Theory, 11 (1964), 210–213. https://doi.org/10.1109/TCT.1964.1082270 doi: 10.1109/TCT.1964.1082270
    [8] J. P. Bouchaud, A. Georges, Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep., 195 (1990), 127–293. https://doi.org/10.1016/0370-1573(90)90099-N doi: 10.1016/0370-1573(90)90099-N
    [9] I. Goychuk, E. Heinsalu, M. Patriarca, G. Schmid, P. Hanggi, Current and universal scaling in anomalous transport, Phys. Rev. E, 73 (2006), 020101. https://doi.org/10.1103/PhysRevE.73.020101 doi: 10.1103/PhysRevE.73.020101
    [10] R. Klages, G. Radons, I. M. Sokolov, Anomalous transport: Foundations and applications, Wiley-VCH Verlag, 2008. https://doi.org/10.1002/9783527622979
    [11] K. M. Saad, M. Alqhtani, Numerical simulation of the fractal-fractional reaction diffusion equations with general nonlinear, AIMS Math., 6 (2021), 3788–3804. https://doi.org/10.3934/math.2021225 doi: 10.3934/math.2021225
    [12] H. M. Srivastava, K. M. Saad, M. Khader, An efficient spectral collocation method for the dynamic simulation of the fractional epidemiological model of the Ebola virus, Chaos Solitons Fractals, 140 (2020), 110174. https://doi.org/10.1016/j.chaos.2020.110174 doi: 10.1016/j.chaos.2020.110174
    [13] H. Amann, Parabolic evolution equations and nonlinear boundary conditions, J. Differ. Equations, 72 (1988), 201–269. https://doi.org/10.1016/0022-0396(88)90156-8 doi: 10.1016/0022-0396(88)90156-8
    [14] C. V. Pao, Finite difference reaction diffusion systems with coupled boundary conditions and time delays, J. Math. Anal. Appl., 272 (2002), 407–434. https://doi.org/10.1016/S0022-247X(02)00145-2 doi: 10.1016/S0022-247X(02)00145-2
    [15] S. Wang, Doubly nonlinear degenerate parabolic systems with coupled nonlinear boundary conditions, J. Differ. Equations, 182 (2002), 431–469. https://doi.org/10.1006/jdeq.2001.4101 doi: 10.1006/jdeq.2001.4101
    [16] X. W. Su, Boundary value problem for a coupled system of nonlinear fractional differential equations, Appl. Math. Lett., 22 (2009), 64–69. https://doi.org/10.1016/j.aml.2008.03.001 doi: 10.1016/j.aml.2008.03.001
    [17] M. J. Li, Y. L. Liu, Existence and uniqueness of positive solutions for a coupled system of nonlinear fractional differential equations, Open J. Appl. Sci., 3 (2013), 53–61. https://doi.org/10.4236/ojapps.2013.31B1011 doi: 10.4236/ojapps.2013.31B1011
    [18] Y. Wang, S. Liang, Q. Wang, Existence results for fractional differential equations with integral and multi-point boundary conditions, Bound. Value Probl., 2018 (2018), 4. https://doi.org/10.1186/s13661-017-0924-4 doi: 10.1186/s13661-017-0924-4
    [19] B. Ahmad, S. Hamdan, A. Alsaedi, S. K. Ntouyas, A study of a nonlinear coupled system of three fractional differential equations with nonlocal coupled boundary conditions, Adv. Differ. Equ., 2021 (2021), 1–21. https://doi.org/10.1186/s13662-021-03440-7 doi: 10.1186/s13662-021-03440-7
    [20] J. H. Wang, H. J. Xiang, Z. G. Liu, Positive solution to nonzero boundary values problem for a coupled system of nonlinear fractional differential equations, Int. J. Differ. Equ., 2010 (2010), 186928. https://doi.org/10.1155/2010/186928 doi: 10.1155/2010/186928
    [21] W. G. Yang, Positive solutions for a coupled system of nonlinear fractional differential equations with integral boundary conditions, Comput. Math. Appl., 63 (2012), 288–297. https://doi.org/10.1016/j.camwa.2011.11.021 doi: 10.1016/j.camwa.2011.11.021
    [22] Y. Wang, L. S. Liu, Y. H. Wu, Positive solutions for a class of higher-order singular semipositone fractional differential systems with coupled integral boundary conditions and parameters, Adv. Differ. Equ., 2014 (2014), 268. https://doi.org/10.1186/1687-1847-2014-268 doi: 10.1186/1687-1847-2014-268
    [23] Z. Odibat, S. Momani, Modified homotopy perturbation method: Application to quadratic riccati differential equation of fractional order, Chaos Solitons Fractals, 36 (2008), 167–174. https://doi.org/10.1016/j.chaos.2006.06.041 doi: 10.1016/j.chaos.2006.06.041
    [24] H. Jafari, V. Daftardar-Gejji, Solving a system of nonlinear fractional differential equations using adomian decomposition, J. Comput. Appl. Math., 196 (2006), 644–651. https://doi.org/10.1016/j.cam.2005.10.017 doi: 10.1016/j.cam.2005.10.017
    [25] S. Momani, Z. Odibat, Numerical approach to differential equations of fractional order, J. Comput. Appl. Math., 207 (2007), 96–110. https://doi.org/10.1016/j.cam.2006.07.015 doi: 10.1016/j.cam.2006.07.015
    [26] X. J. Zhou, C. J. Xu, Numerical solution of the coupled system of nonlinear fractional ordinary differential equations, Adv. Appl. Math. Mech., 9 (2017), 574–595. https://doi.org/10.4208/aamm.2015.m1054 doi: 10.4208/aamm.2015.m1054
    [27] C. L. Wang, Z. Q. Wang, L. L. Wang, A spectral collocation method for nonlinear fractional boundary value problems with a caputo derivative, J. Sci. Comput., 76 (2018), 166–188. https://doi.org/10.1007/s10915-017-0616-3 doi: 10.1007/s10915-017-0616-3
    [28] A. Kilbas, H. M. Srivastava, J. J. Trujillo, Theory and applications of fractional differential equations, Elsevier, 2006.
    [29] J. Shen, T. Tang, Z. H. Teng, Spectral and high-order methods with applications, 2006.
  • This article has been cited by:

    1. Hamid Reza Marzban, Atiyeh Nezami, Analysis of nonlinear fractional optimal control systems described by delay Volterra–Fredholm integral equations via a new spectral collocation method, 2022, 162, 09600779, 112499, 10.1016/j.chaos.2022.112499
    2. S Yaghoubi, H Aminikhah, K Sadri, A spectral shifted gegenbauer collocation method for fractional pantograph partial differential equations and its error analysis, 2023, 48, 0973-7677, 10.1007/s12046-023-02270-5
    3. Waleed Mohamed Abd-Elhameed, Youssri Hassan Youssri, Amr Kamel Amin, Ahmed Gamal Atta, Eighth-Kind Chebyshev Polynomials Collocation Algorithm for the Nonlinear Time-Fractional Generalized Kawahara Equation, 2023, 7, 2504-3110, 652, 10.3390/fractalfract7090652
  • 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(2599) PDF downloads(128) Cited by(3)

Figures and Tables

Figures(9)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog