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

1.   Introduction

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 $\alpha \in (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 $\alpha \in (0, 1)$ and rarely on $\alpha \in (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 $\alpha \in (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.

2.   Preliminaries

Definition 2.1.(cf. p. 70 ^[28]) Let $t \in (0, 1)$, the left-sided Caputo derivative of order $\alpha$, $n-1 < \alpha < n$, $n \in \mathbb{N^{+}}$, are defined as:

where $\Gamma(\cdot)$ denotes Gamma function.

Definition 2.2.(cf. ^[29]) Let $\mathcal{J}_n^{\alpha, \beta}(x), x \in \Lambda$ be the standard Jacobi polynomial of degree $n$. The set of Jacobi polynomials is a complete $L^2_{\omega^{\alpha, \beta}}(\Lambda)$-orthogonal system, i.e.

where $\delta_{m, n}$ is the Kronecker function, and

Specially,

Definition 2.3.(cf. ^[27]) For any of the Gauss-type quadratures defined above with the points and weights $\{x_n^{\alpha, \beta}, \omega_n^{\alpha, \beta}\}_{n = 0}^N(N \geq 0)$, we can define a discrete inner product in interval $\Lambda$ :

Lemma 2.1.(cf. ^[29]) Let $\mathcal{P}_N$ be the space of all polynomials of degree at most $N$, which is exact for any $\phi(x) \in \mathcal{P}_{2N+1}$. Particularly,

Lemma 2.2.(Lemma 2.1. ^[27]) For boundary value problem $^C_0D^{\alpha}_{t}y(t) = f(t, y(t)), \; t \in (0, T)$ with $y(0) = y(T) = 0$, Let $1 < \alpha < 2$. Assume that $y(t)$ is a function with an absolutely continuous first derivative, and $f:[0, T] \times \mathbb{R} \rightarrow \mathbb{R}$ is continuous. Then we have that $y \in C^1[0, T]$ is a solution of the problem if and only if it is a solution of the Fredholm integral equation:

3.   A kinds of numerical schemes for the coupled system of nonlinear fractional differential equations

3.1.   Reformulation of the coupled system

We consider the coupled system of nonlinear fractional differential equations:

where $f_1, f_2:[0, T] \times \mathbb{R} \rightarrow \mathbb{R}$ is continuous, and $^C_0D^{\alpha_1}_{t}$, $^C_0D^{\alpha_2}_{t}$ is the left-sided Caputo derivative of order $\alpha_1, \alpha_2 \in$ (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 $\alpha_i, y_i(t), f_i$ satisfy the condition of Lemma 2.2:

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

Lemma 3.1. Let $\alpha_i, y_i(t), f_i$ satisfy the condition of Lemma 2.2. Futhermore, let $f_i$ satisfy a Lipschitz condition with the Lipschitz constant $L < \cfrac{\Gamma(\alpha+1)}{4T^{\alpha}}$. $\cfrac{\Gamma(\alpha+1)}{T^{\alpha}} = \max(\cfrac{\Gamma(\alpha_1+1)}{T^{\alpha_1}}, \cfrac{\Gamma(\alpha_2+1)}{T^{\alpha_2}})$. Then the system (3.1) has a unique solution.

Proof. Let $\mathrm{y}(t) = (y_1(t), y_2(t))$, define the operator $\mathcal{A}\mathrm{y}(t) = (\mathcal{A}_1\mathrm{y}(t), \mathcal{A}_2\mathrm{y}(t))$. where

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

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

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:

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

Using the abbreviation, (3.3) can be read to

Finally, under the linear transformation

In summary, we obtain

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

3.2.   The Legendre spectral collocation scheme

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

Obviously

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

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

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

where

Futhermore, by (2.2) we obtain

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

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

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.

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

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

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.

4.   Error analysis

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

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,

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

where

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

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

Then, combination with Lemma 4.1

We set the Nemytskii operator

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

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

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

Then by Lemma 4.1, we get

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

Clearly,

Then, we easily get

Thus

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

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

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

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

Proof. Similarly, we get that

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

By (4.6) and Lemma 4.2, we have

Then, combination with (4.4) we can obtain

Thus

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

5.   Numerical experiments

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:

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

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

$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

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.

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

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.

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

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 5–7, 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$.

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

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.

6.   Conclusions

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.

Acknowledgments

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)

Conflict of interest

The authors declare that there are no conflicts of interest.