Time-Stepping Error Estimates of Linearized Grünwald–Letnikov Difference Schemes for Strongly Nonlinear Time-Fractional Parabolic Problems

Abstract

A fully discrete scheme is proposed for numerically solving the strongly nonlinear time-fractional parabolic problems. Time discretization is achieved by using the Grünwald–Letnikov (G–L) method and some linearized techniques, and spatial discretization is achieved by using the standard second-order central difference scheme. Through a Grönwall-type inequality and some complementary discrete kernels, the optimal time-stepping error estimates of the proposed scheme are obtained. Finally, several numerical examples are given to confirm the theoretical results.

1. Introduction

In this paper, we present time-stepping error estimates of the Grünwald–Letnikov (G–L) difference scheme for the following strongly nonlinear time-fractional parabolic problem

$${\partial }_t^{\alpha }u=q\left(u\right)\sum _{i=1}^d\frac{{\partial }^{2}u}{\partial x_i^2}+f\left(u\right),(\mathit{x},t)\in \Omega \times (0,T]$$

(1)

with homogeneous Dirichlet boundary conditions and the initial condition $u(\mathit{x},0)=u_0\left(\mathit{x}\right)$, where $\mathit{x}\in \Omega \subset {\mathbb{R}}^d$$(d=1,2,3)$, $q\left(u\right)\ge 0$, $f\left(u\right)\in C\left(\mathbb{R}\right)$ is the nonlinear function, and ${\partial }_t^{\alpha }u$ is the Riemann–Liouville fractional derivative with $0<\alpha <1$. The fractional equations have been widely used as mathematical models in the fields of physics [1,2,3,4], economics [5,6], biology [7,8], and so on.

In recent years, there have been plenty of numerical results on the error estimates of different schemes, e.g., L1 schemes [9,10], L2-type schemes [11,12], fast time-stepping schemes [13,14], and so on. These convergence results usually give maximum errors for time-fractional problems. Recently, many numerical methods have given a lower convergence order as t tends to zero and a higher convergence order as t is far away from the beginning. This results in investigations into the time-stepping error estimates of different numerical algorithms. For example, Gracia et al. [15] studied the time-stepping error estimates of the L1 scheme on uniform time mesh to solve the linear fractional convection–diffusion equation and proved that the time error estimate is proportional to $\tau t_n^{\alpha -1}$, where $t_n=n\tau$, with $\tau$ being the time step. Kopteva [16] further considered the time-stepping error estimate of the L1 scheme on quasi-uniform temporal mesh. Li et al. [17] developed some fractional discrete Grönwall-type inequalities and obtained a time-stepping error estimate of the L1 scheme under general regularity assumptions. Yan et al. [18] obtained time-stepping error estimates of a modified L1 scheme. Jin et al. [19] obtained time-stepping error estimates of the k-th backward differentiation formula (BDF) convolution quadrature method for solving time-fractional equations. Zhang et al. [20] studied the time-stepping error estimates of the fractional BDF method and the fractional Crank–Nicolson-type method. Li et al. [21] further improved the results by using a novel discrete fractional Grönwall-type inequality. Santra and Mohapatra [22] solved the multi-term time-fractional differential equation by using the L1 scheme and obtained related time-stepping error estimates. Recently, Chen et al. [23] studied the time-stepping error estimates of the G–L scheme for the equations ${\partial }_t^{\alpha }u+c\left(t\right)u\left(t\right)=f\left(t\right)$, where $c\left(t\right)>0$. In the case of $c\left(t\right)>0$, the errors decrease in time and the principle of induction can be applied. The analytical methods do not work for the general nonlinear fractional differential Equation (1). To the best of the authors’ knowledge, the time-stepping error estimates of the G–L scheme with maximum norm are still missing for nonlinear time-fractional problems.

In this paper, we propose a fully discrete linearized numerical scheme to solve the nonlinear fractional differential Equation (1). In the temporal direction, the G–L method and the linearized method are applied. In the spatial direction, the second-order central difference scheme is used to approximate the spatial derivatives. A sharp time-stepping error estimate of the proposed scheme is obtained. The results imply that the temporal convergence order is of ${\tau }^{\alpha }$ as t tends to 0 and of 1 as t is far away from $t=0$. Numerical results are given to confirm the theoretical results.

The key to the proof of the convergence results is the complementary discrete kernels, which were proposed to investigate the maximum errors of L1-type schemes for time-fractional problems [10,24]. In this study, we continue to study the complementary discrete kernels for the G–L scheme and obtain time-stepping truncation errors for the time-fractional problems. Thanks to the time-stepping truncation errors and a Grönwall-type inequality, we obtain the optimal time-stepping error estimates of the G–L scheme for the nonlinear problems.

The paper is organized as follows: In Section 2, we construct the fully discrete numerical scheme for problems (1). In Section 3, the convergence results of the scheme are rigorously proved. In Section 4, several numerical examples are provided to confirm the theoretical results in the paper. Throughout the paper, we declare that C and $C_i$$(i=1,2)$ have different values in different places.

2. Numerical Scheme

In this section, we derive a full discrete scheme for solving Equation (1) numerically. Without loss of generality, we set $d=2$ and $\Omega ={(a,b)}^2$.

Let $\tau =T/N$ be the time step and $h_x=(b-a)/M_x$, $h_y=(b-a)/M_y$ the meshsize, where N, $M_x$, and $M_y$ are positive integers. Denote $t_n=n\tau ,n=0,1,2,\dots ,N$, $x_i=a+i{h}_x,i=0,1,2,\dots ,M_x$, $y_j=a+j{h}_y,j=0,1,2,\dots ,M_y$, ${\Omega }_{\tau }=\{t_n|0\le n\le N\}$, ${\overline{\Omega }}_h=\{(x_i,y_j)|0\le i\le M_x,0\le j\le M_y\}$, ${\Omega }_h={\overline{\Omega }}_h\cap \Omega$, $\partial {\Omega }_h={\overline{\Omega }}_h\cap \partial \Omega$, $\omega =\{(i,j)|(x_i,y_j)\in {\Omega }_h\}$, $\partial \omega =\{(i,j)|(x_i,y_j)\in \partial {\Omega }_h\}$, $\overline{\omega }=\omega \cup \partial \omega$. Let ${\mathcal{U}}_h=\{v=\left(v_{i,j}\right)|0\le i\le M_x,0\le j\le M_y\}$ be the grid function space which is defined on ${\overline{\Omega }}_h$. For $v\in {\mathcal{U}}_h$, we introduce some notations:

$${\delta }_{▵}^2{v}_{i,j}=\frac{1}{h_x^2}(v_{i+1,j}-2v_{i,j}+v_{i-1,j})+\frac{1}{h_y^2}(v_{i,j+1}-2v_{i,j}+v_{i,j-1}){,\parallel v\parallel }_{\infty }=\underset{\begin{array}{c}1\le i\le M_x\\ 1\le j\le M_y\end{array}}{max}|v_{i,j}|.$$

As in ([25], p. 6, (16)), we assume that the solution $u(x,y,t)$ of (1) has the following form:

$$\begin{array}{c}\hfill u(x,y,t)=\sum _{{(j,k)}_{\alpha }}{\gamma }_{j,k}{t}^{j+k\alpha }\phi (x,y)+Y(x,y,t),0\le t\le T,\end{array}$$

(2)

where ${(j,k)}_{\alpha }:=\{(j,k):j,k\in {\mathbb{N}}_0,j+k\alpha <2\}$, ${\mathbb{N}}_0=\{0,1,2,\cdots \}$, ${\gamma }_{j,k}$ are some constants, $Y(x,y,t)\in C^2[0,T]$. Here, the functions $\phi (x,y)$ and $Y(x,y,t)$ are smooth in the spatial direction.

The G–L approximation ([23], p. 54, (2.5)) of ${\partial }_t^{\alpha }u(·,t_n)$ is as follows:

$$\begin{array}{ccc}\hfill {\partial }_t^{\alpha }u(·,t_n)& \approx & A_{0}u(·,t_n)-\sum _{k=1}^{n-1}(A_{n-k-1}-A_{n-k})u(·,t_k)-A_{n-1}u(·,t_0)\hfill \\ & =& \sum _{k=1}^n{A}_{n-k}\left(u(·,t_k)-u(·,t_{k-1})\right)\hfill \\ & :=& D_{\tau }^{\alpha }u(·,t_n),\hfill \end{array}$$

(3)

where $A_{k-1}=\frac{1}{{\tau }^{\alpha }}{d}_k^{\left(\alpha \right)}$ with $d_k^{\left(\alpha \right)}=\frac{\Gamma (k-\alpha )}{\Gamma (1-\alpha )\Gamma \left(k\right)}$ for $k\ge 1$. The truncation errors of the G–L approximation is stated below.

Lemma 1 

([23], p. 56, Lemma 3.2). Under the assumption (2), it holds that

$$|{\partial }_t^{\alpha }u(·,t_n)-D_{\tau }^{\alpha }u(·,t_n)|\le C(n^{-2}+{\tau }^{\gamma -\alpha }{n}^{-(1+\alpha -\gamma )}),\gamma =min\{1,2\alpha \},1\le n\le N.$$

Now, we begin to present the fully discrete scheme of (1). First, we consider Equation (1) at the grid point $(x_i,y_j,t_n)$. For $(i,j)\in \omega$ and $1\le n\le N$, it holds that

$$D_{\tau }^{\alpha }u(x_i,y_j,t_n)=q\left(u(x_i,y_j,t_{n-1})\right){\delta }_{▵}^{2}u(x_i,y_j,t_n)+f\left(u(x_i,y_j,t_{n-1})\right)+R_{i,j}^n,$$

(4)

where $R_{i,j}^n$ is the truncation error. Using Lemma 1, one can obtain

$$R_{i,j}^n\le C\left((h_x^2+h_y^2)+\tau +n^{-2}+{\tau }^{\gamma -\alpha }{n}^{-(1+\alpha -\gamma )}+\tau (h_x^2+h_y^2)\right).$$

(5)

With the notation $U_{i,j}^n=u(x_i,y_j,t_n)$, Equation (4) can be rewritten as

$$D_{\tau }^{\alpha }{U}_{i,j}^n=q\left(U_{i,j}^{n-1}\right){\delta }_{\Delta }^2{U}_{i,j}^n+f\left(U_{i,j}^{n-1}\right)+R_{i,j}^n.$$

(6)

Omitting small term $R_{i,j}^n$ and replacing $U_{i,j}^n$ with $u_{i,j}^n$, we obtain the following numerical scheme: for $(i,j)\in \partial \omega ,0\le n\le N$,

$$\begin{array}{cc}\hfill D_{\tau }^{\alpha }{u}_{i,j}^n=q\left(u_{i,j}^{n-1}\right){\delta }_{▵}^2{u}_{i,j}^n+f\left(u_{i,j}^{n-1}\right),& (i,j)\in \omega ,1\le n\le N,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill u_{i,j}^0=u_0(x_i,y_j),& (i,j)\in \overline{\omega },\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill u_{i,j}^n=0,& (i,j)\in \partial \omega ,0\le n\le N.\hfill \end{array}$$

(9)

3. Numerical Analysis

In this section, the time-stepping error estimates of scheme (7)–(9) are rigorously proved. For the G–L approximation, one can check that

$$A_0\ge A_1\ge A_2\ge \cdots \ge A_{m-1}>0,1\le m\le M.$$

Now, we construct a family of complementary discrete kernels $P_{n-j}$ satisfying

$$\sum _{k=m}^j{P}_{j-k}{A}_{k-m}\equiv 1,1\le m\le j\le M,$$

(10)

which implies, for given coefficients $A_{k-m}$,

$$\begin{array}{c}\hfill P_0:=\frac{1}{A_0},P_j:=\frac{1}{A_0}\sum _{k=0}^{j-1}(A_{j-k-1}-A_{j-k})P_k.\end{array}$$

(11)

Below, we present several lemmas for the complementary discrete kernels.

Lemma 2 

([23], p. 57, Lemma 4.1, Corollary 4.1). It holds that

$$\begin{array}{c}\hfill \frac{{\tau }^{\alpha }{(n+1)}^{\alpha -1}}{\Gamma \left(\alpha \right)}<P_n=\frac{{\tau }^{\alpha }\Gamma (n+\alpha )}{\Gamma \left(\alpha \right)\Gamma (n+1)}<\frac{{\tau }^{\alpha }{n}^{\alpha -1}}{\Gamma \left(\alpha \right)},n\ge 1.\end{array}$$

(12)

Lemma 3. 

It holds that

$$\begin{array}{c}\hfill \sum _{k=1}^n|R_{i,j}^k|P_{n-k}\le C\left(t_n^{\alpha }(h_x^2+h_y^2)+\tau t_n^{\alpha -1}\right).\end{array}$$

(13)

Proof. 

The desired result is an immediate consequence of Lemma 4.3 in [23] and inequality (5). □

Lemma 4 

([26], p. 625, Theorem 2.1). Let ${\left\{w_i\right\}}_{i=0}^N$ be a sequence of non-negative real numbers satisfying

$$w_i\le \eta +\frac{\phi }{{\left(i\tau \right)}^{1-\alpha }}+M{\tau }^{\alpha }\sum _{j=0}^{i-1}\frac{w_j}{{(i-j)}^{1-\alpha }},$$

where $0<\alpha <1$, φ and η are non-negative constants and M is a positive constant independent of τ. Then, it holds

$$w_i\le \eta E_{\alpha }(M\Gamma \left(\alpha \right){\left(i\tau \right)}^{\alpha })+\frac{\phi \Gamma \left(\alpha \right)}{{\left(i\tau \right)}^{1-\alpha }}\sum _{n=0}^{\infty }\frac{{(M\Gamma \left(\alpha \right){\left(i\tau \right)}^{\alpha })}^n}{\Gamma \left(\alpha \right(n+1\left)\right)},0\le i\le N,$$

where $E_{\alpha }\left(z\right)={\sum }_{k=0}^{\infty }\frac{z^k}{\Gamma (1+k\alpha )}$ is the Mittag-Leffler function.

Theorem 1. 

Let $U^m=\left(U_{i,j}^m\right)=u(x_i,y_j,t_n)$ be the exact solutions of problelms (1). Let $u^m=\left(u_{i,j}^m\right)\in {\mathcal{U}}_h$ be numerical solution of scheme (7)–(9). Then, there exists a positive constant ${\tau }^{*}$ such that when $\tau \le {\tau }^{*}$, the following error estimate holds

$${\parallel U^m-u^m\parallel }_{\infty }\le C\left(\tau t_m^{\alpha -1}+t_n^{\alpha }(h_x^2+h_y^2)\right),0\le m\le N.$$

(14)

Proof. 

We prove the main results by using mathematical induction. First, the time-stepping error holds for $m=0$. Suppose that the error estimates (14) hold for $m=0,1,\cdots ,n-1$. Then, substituting (6) from (7) and denoting $e_{i,j}^m=u_{i,j}^m-U_{i,j}^m$, we can obtain

$$\begin{array}{cc}& D_{\tau }^{\alpha }{e}_{i,j}^m=q\left(u_{i,j}^{m-1}\right){\delta }_{▵}^2{u}_{i,j}^m-q\left(U_{i,j}^{m-1}\right){\delta }_{\Delta }^2{U}_{i,j}^m+f\left(u_{i,j}^{m-1}\right)-f\left(U_{i,j}^{m-1}\right)+R_{i,j}^m,(i,j)\in \omega .\end{array}$$

(15)

For the first two terms in the right-hand side of (15), we have

$$\begin{array}{cc}\hfill & q\left(u_{i,j}^{m-1}\right){\delta }_{▵}^2{u}_{i,j}^m-q\left(U_{i,j}^{m-1}\right){\delta }_{\Delta }^2{U}_{i,j}^m\hfill \\ \hfill =& q\left(u_{i,j}^{m-1}\right){\delta }_{▵}^2{u}_{i,j}^m-q\left(u_{i,j}^{m-1}\right){\delta }_{\Delta }^2{U}_{i,j}^m+q\left(u_{i,j}^{m-1}\right){\delta }_{\Delta }^2{U}_{i,j}^m-q\left(U_{i,j}^{m-1}\right){\delta }_{\Delta }^2{U}_{i,j}^m\hfill \\ \hfill =& q\left(u_{i,j}^{m-1}\right){\delta }_{▵}^2{e}_{i,j}^m+\left(q\left(u_{i,j}^{m-1}\right)-q\left(U_{i,j}^{m-1}\right)\right){\delta }_{\Delta }^2{U}_{i,j}^m.\hfill \end{array}$$

Thus, Equation (15) can be further rewritten as

$$D_{\tau }^{\alpha }{e}_{i,j}^m=q\left(u_{i,j}^{m-1}\right){\delta }_{▵}^2{e}_{i,j}^m+\left(q\left(u_{i,j}^{m-1}\right)-q\left(U_{i,j}^{m-1}\right)\right){\delta }_{\Delta }^2{U}_{i,j}^m+f\left(u_{i,j}^{m-1}\right)-f\left(U_{i,j}^{m-1}\right)+R_{i,j}^m.$$

(16)

Let $e^m=\left(e_{i,j}^m\right)\in {\mathcal{U}}_h$ and suppose that $e_{i_m,j_m}^m={\parallel e^m\parallel }_{\infty }$. At the grid point $(x_{i_m},y_{j_m},t_m)$, we have

$$D_{\tau }^{\alpha }{e}_{i_m,j_m}^m=q\left(u_{i_m,j_m}^{m-1}\right){\delta }_{▵}^2{e}_{i_m,j_m}^m+\left(q\left(u_{i_m,j_m}^{m-1}\right)-q\left(U_{i_m,j_m}^{m-1}\right)\right){\delta }_{\Delta }^2{U}_{i_m,j_m}^m+f\left(u_{i_m,j_m}^{m-1}\right)-f\left(U_{i_m,j_m}^{m-1}\right)+R_{i_m,j_m}^m,$$

which implies

$$\begin{array}{cc}\hfill & \left(A_0+\frac{2q\left(u_{i_m,j_m}^{m-1}\right)}{h_x^2}+\frac{2q\left(u_{i_m,j_m}^{m-1}\right)}{h_y^2}\right)e_{i_m,j_m}^m\hfill \\ \hfill =& q\left(u_{i_m,j_m}^{m-1}\right)\frac{e_{i_m+1,j_m}^m+e_{i_m-1,j_m}^m}{h_x^2}+q\left(u_{i_m,j_m}^{m-1}\right)\frac{e_{i_m,j_m+1}^m+e_{i_m,j_m-1}^m}{h_y^2}\hfill \\ \hfill & +\sum _{k=1}^{m-1}(A_{m-k-1}-A_{m-k})e_{i_m,j_m}^k+\left(q\left(u_{i_m,j_m}^{m-1}\right)-q\left(U_{i_m,j_m}^{m-1}\right)\right){\delta }_{\Delta }^2{U}_{i_m,j_m}^m\hfill \\ \hfill & +f\left(u_{i_m,j_m}^{m-1}\right)-f\left(U_{i_m,j_m}^{m-1}\right)+R_{i_m,j_m}^m.\hfill \end{array}$$

From the fact that $A_k>A_{k+1}$, it follows that

$$\begin{array}{cc}\hfill & \left(A_0+\frac{2q(u_{i_m,j_m}^{m-1})}{h_x^2}+\frac{2q(u_{i_m,j_m}^{m-1})}{h_y^2}\right)|e_{i_m,j_m}^m|\hfill \\ \hfill \le & q(u_{i_m,j_m}^{m-1})\left|\frac{e_{i_m+1,j_m}^m+e_{i_m-1,j_m}^m}{h_x^2}\right|+q(u_{i_m,j_m}^{m-1})\left|\frac{e_{i_m,j_m+1}^m+e_{i_m,j_m-1}^m}{h_y^2}\right|\hfill \\ \hfill & +\sum _{k=1}^{m-1}(A_{m-k-1}-A_{m-k})|e_{i_m,j_m}^k|+\left|\left(q(u_{i_m,j_m}^{m-1})-q(U_{i_m,j_m}^{m-1})\right){\delta }_{\Delta }^2{U}_{i_m,j_m}^m\right|\hfill \\ \hfill & +|f(u_{i_m,j_m}^{m-1})-f(U_{i_m,j_m}^{m-1})|+|R_{i_m,j_m}^{m-1}|.\hfill \end{array}$$

Noting that $|e_{i_m+1,j_m}^m|+|e_{i_m-1,j_m}^m|\le 2|e_{i_m,j_m}^m|$, $|e_{i_m,j_m+1}^m|+|e_{i_m,j_m-1}^m|\le 2|e_{i_m,j_m}^m|$, we have

$$\begin{array}{c}\hfill A_0|e_{i_m,j_m}^m|\le \sum _{k=1}^{m-1}(A_{m-k-1}-A_{m-k})|e_{i_m,j_m}^k|+|f(u_{i_m,j_m}^{m-1})-f(U_{i_m,j_m}^{m-1})|+C|e_{i_m,j_m}^{m-1}|+|R_{i_m,j_m}^m|.\end{array}$$

Now, by the assumptions that (14) hold when $n=m+1$, we conclude that

$$|U_{i_m,j_m}^{n-1}|\le |u_{i_m,j_m}^{n-1}|+1.$$

whenever $\tau$, $h_x$ and $h_y$ are sufficiently small. As a result, there exists a constant C such that

$$|f(u_{i_m,j_m}^{m-1})-f(U_{i_m,j_m}^{m-1})|<C|u_{i_m,j_m}^{m-1}-u_{i_m,j_m}^{m-1}|=C{e}_{i_m,j_m}^{m-1}.$$

Therefore,

$$\begin{array}{c}\hfill A_0\parallel e^m{\parallel }_{\infty }\le \sum _{k=1}^{m-1}(A_{m-k-1}-A_{m-k})\parallel e^k{\parallel }_{\infty }+C\parallel e^{m-1}{\parallel }_{\infty }+|R_{i_m,j_m}^m|.\end{array}$$

Multiplying both sides of the above inequality by $P_{n-m}$, summing over for m from 1 to n and using equality (3), we arrive at

$$\begin{array}{c}\hfill \sum _{m=1}^n{P}_{n-m}\sum _{k=1}^m{A}_{m-k}(\parallel e^k{\parallel }_{\infty }-\parallel e^{k-1}{\parallel }_{\infty })\le \sum _{m=1}^n{P}_{n-m}\left(C\parallel e^{m-1}{\parallel }_{\infty }+|R_{i_m,j_m}^m|\right).\end{array}$$

(17)

Applying identity (10), we rewrite the left-hand side of inequality (17) as

$$\begin{array}{cc}\hfill \sum _{m=1}^n{P}_{n-m}\sum _{k=1}^m{A}_{m-k}(\parallel e^k{\parallel }_{\infty }-\parallel e^{k-1}{\parallel }_{\infty })& =\sum _{k=1}^n(\parallel e^k{\parallel }_{\infty }-\parallel e^{k-1}{\parallel }_{\infty })\sum _{m=k}^n{P}_{n-m}{A}_{m-k}\hfill \\ & = \parallel e^n{\parallel }_{\infty }-{\parallel e^0\parallel }_{\infty }.\hfill \end{array}$$

(18)

Substituting (18) into inequality (17) leads to

$$\begin{array}{c}\hfill {\parallel e^n\parallel }_{\infty }-{\parallel e^0\parallel }_{\infty }\le \sum _{m=1}^n{P}_{n-m}\left(C{\parallel e^{m-1}\parallel }_{\infty }+|R_{i_m,j_m}^m|\right).\end{array}$$

(19)

Using Lemmas 2 and 3, and noticing the fact $\frac{1}{{(n-m)}^{1-\alpha }}\le \frac{2^{1-\alpha }}{{(n-m+1)}^{1-\alpha }}$ for $n>m$, one can obtain the following from (19)

$$\begin{array}{c}\hfill \parallel e^n{\parallel }_{\infty }\le C_1{\tau }^{\alpha }\left(\sum _{m=1}^n\frac{\parallel e^m{\parallel }_{\infty }}{{(n-m+1)}^{1-\alpha }}+\sum _{m=1}^n\frac{\parallel e^{m-1}{\parallel }_{\infty }}{{(n-m+1)}^{1-\alpha }}\right)+C_2\left(t_n^{\alpha }(h_x^2+h_y^2)+\tau t_n^{\alpha -1}\right).\end{array}$$

Thanks to Lemma 4, we can derive

$$\begin{array}{c}\hfill {\parallel e^n\parallel }_{\infty }\le C\left(\tau t_n^{\alpha -1}+t_n^{\alpha }(h_x^2+h_y^2)\right).\end{array}$$

Therefore, the error estimates hold for $m=n$. This closes the mathematical induction and completes the proof. □

4. Numerical Experiments

In this section, we conduct several numerical examples to verify the theoretical results in the paper. All the computations are performed by MATLAB 2018b.

Example 1. 

Consider the following one-dimensional nonlinear fractional differential equation:

$$\begin{array}{cc}\hfill {\partial }_t^{\alpha }u=u_{xx}+u-u^2+g(x,t),& 0<x<1,0<t\le T,\hfill \end{array}$$

where $g(x,t)$, initial boundary conditions are given by the exact solution

$$u(x,t)=t^{\alpha }sin\left(\pi x\right).$$

In order to test the convergence order in time direction, we set $M_x=1000$ and $T=1$. The errors and convergence orders at the last time with $N=20,40,80,160,320$ are displayed in

Table 1. The errors and convergence orders at $t=1$ in the time direction for Example 1.

N	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
20	1.9119E-02	–		2.8577E-02	–		3.8115E-02	–
40	9.4595E-03	1.0152		1.4165E-02	1.0152		1.8921E-02	1.0104
80	4.7039E-03	1.0079		7.0520E-03	1.0063		9.4227E-03	1.0051
160	2.3448E-03	1.0044		3.5180E-03	1.0044		4.7049E-03	1.0026
320	1.1702E-03	1.0027		1.7567E-03	1.0019		2.3500E-03	1.0015

. The maximum errors and convergence orders at the last time with $N=320,M_x=1000$, $T={10}^{-3},{10}^{-4},{10}^{-5},{10}^{-6},{10}^{-7},{10}^{-8},{10}^{-9}$ are listed in

Table 2. The errors and convergence orders at the last time in the time direction for Example 1.

T	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
${10}^{-3}$	6.2044E-04	–		1.1365E-05	–		1.4471E-06	–
${10}^{-4}$	1.7885E-05	0.5402		1.8707E-06	0.7836		1.6928E-07	0.9319
${10}^{-5}$	5.3146E-06	0.5270		3.9938E-07	0.6706		2.5282E-08	0.8258
${10}^{-6}$	1.7576E-06	0.4806		9.5724E-08	0.6204		3.9679E-09	0.8042
${10}^{-7}$	6.3740E-07	0.4405		2.3752E-08	0.6053		6.2788E-10	0.8007
${10}^{-8}$	2.4350E-07	0.4179		5.9478E-09	0.6013		9.9488E-11	0.8001
${10}^{-9}$	9.5287E-08	0.4075		1.4929E-09	0.6003		1.5767E-11	0.8000

. These results show that the convergence order of the scheme is of 1 when t is far away from the initial time and of $\alpha$ when t is close to the initial time. We test the convergence order in the spatial direction with $N=\mathrm{10,000}$ and $M_x=4,8,16,32,64$. The maximum errors and convergence orders at $t=1$ are listed in

Table 3. The errors and convergence orders at $t=1$ in the spatial direction for Example 1.

N	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
4	4.4921E-02	–		4.4905E-02	–		4.4776E-02	–
8	1.1055E-02	2.0227		1.1036E-02	2.0246		1.0990E-02	2.0266
16	2.7260E-03	2.0198		2.7072E-03	2.0274		2.6813E-03	2.0351
32	6.5285E-04	2.0620		6.3405E-04	2.0941		6.1335E-04	2.1282
64	1.3513E-04	2.2724		1.1635E-04	2.4462		9.6931E-05	2.6617

. The results indicate that the convergence order is of 2 in the spatial direction. These numerical results verify the theoretical results in the paper.

Example 2. 

Consider the following two-dimensional nonlinear fractional differential equation:

$${\partial }_t^{\alpha }u=(1+u^2)(u_{xx}+u_{yy})+u-u^3+g(x,y,t),(x,y)\in {(0,1)}^2,t\in (0,T],$$

subject to homogeneous boundary conditions. The initial condition and $g(x,y,t)$ are determined by using the exact solution

$$u(x,y,t)=(x-x^2)(y-y^2)(t^{\alpha }+t^2).$$

In order to test the convergence order in time direction, we set $M_x=M_y=40$ and $T=1$. The errors and convergence orders at the last time with $N=20,40,80,160,320$ are displayed in

Table 4. The errors and convergence orders at $t=1$ in the time direction for Example 2.

N	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
20	7.5426E-03	–		8.1545E-03	–		8.7658E-03	–
40	3.8026E-03	0.9881		4.1081E-03	0.9891		4.4155E-03	0.9893
80	1.9090E-03	0.9942		2.0617E-03	0.9947		2.2158E-03	0.9948
160	9.5638E-04	0.9971		1.0327E-03	0.9974		1.1099E-03	0.9974
320	4.7865E-04	0.9986		5.1682E-04	0.9987		5.5543E-04	0.9987

. The maximum errors and convergence orders at the last time with $N=10,M_x=M_y=40$, $T={10}^{-4},{10}^{-5},{10}^{-6},{10}^{-7},{10}^{-8},{10}^{-9}$ are listed in

Table 5. The errors and convergence orders as $t\to 0$ in the time direction for Example 2.

T	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
${10}^{-4}$	4.5550E-05	–		4.2083E-06	–		3.4867E-07	–
${10}^{-5}$	1.2815E-05	0.5508		8.2348E-07	0.7085		5.0240E-08	0.8414
${10}^{-6}$	3.8886E-06	0.5179		1.9169E-07	0.6330		7.8361E-09	0.8069
${10}^{-7}$	1.3330E-06	0.4650		4.7188E-08	0.6088		1.2388E-09	0.8011
${10}^{-8}$	4.9542E-07	0.4299		1.1792E-08	0.6022		1.9625E-10	0.8002
${10}^{-9}$	1.9157E-07	0.4126		2.9583E-09	0.6006		3.1102E-11	0.8000

. Again, we can see that the scheme is of order 1 when t is far away from the initial time and of order $\alpha$ when t is near the initial time.

Example 3. 

Consider the following one-dimensional nonlinear fractional differential equation with homogeneous boundary conditions:

$$\begin{array}{c}\hfill {\partial }_t^{\alpha }u=(1+u^4)u_{xx}+u-u^3,0<x<1,0<t\le T,\end{array}$$

(20)

whose initial value is taken as

$$\begin{array}{c}\hfill u(x,0)=2xsin\left(2\pi x\right),0\le x\le 1.\end{array}$$

In this example, we set $h_x=0.1$. The numerical solutions computed by the proposed scheme with $N=2048$ and $\tau =T/N$ are used as the reference solutions in the accuracy test. To verify the convergence order far away from initial time, the errors and convergence orders at time $t=1$ are contained in

Table 6. The errors and convergence orders at $t=1$ in the time direction for Example 3.

N	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
8	3.8667E-05	–		5.2839E-05	–		5.2703E-05	–
16	1.8713E-05	1.0471		2.5169E-05	1.0699		5.4107E-05	1.1284
32	9.1542E-06	1.0315		1.2221E-05	1.0423		1.1499E-05	1.0679
64	4.4738E-06	1.0329		5.9512E-06	1.0381		5.5528E-06	1.0502
128	2.1575E-06	1.0521		2.8649E-06	1.0547		2.6622E-06	1.0606

. To test the convergence order near initial time, the errors and convergence orders at the time $t={10}^{-4},{10}^{-5},{10}^{-6},{10}^{-7},{10}^{-8},{10}^{-9}$ with $N=10$ are displayed in

Table 7. The errors and convergence orders as $t\to 0$ in the time direction for Example 3.

T	$\mathit{\alpha }$ = 0.4			$\mathit{\alpha }$ = 0.6			$\mathit{\alpha }$ = 0.8
	$\mathit{errors}$	Orders		$\mathit{errors}$	Orders		$\mathit{errors}$	Orders
${10}^{-4}$	6.9689E-05	–		3.3337E-05	–		1.6673E-05	–
${10}^{-5}$	3.8249E-05	0.2605		1.9890E-05	0.2243		2.1532E-06	0.8889
${10}^{-6}$	2.1335E-05	0.2535		6.7966E-06	0.4663		2.8157E-07	0.8835
${10}^{-7}$	1.7114E-05	0.0958		1.7214E-06	0.5964		4.2738E-08	0.8188
${10}^{-8}$	1.0875E-05	0.1969		4.2407E-07	0.6085		6.7243E-09	0.8032
${10}^{-9}$	5.6239E-06	0.2864		1.0576E-07	0.6031		1.0645E-09	0.8005

. Clearly, the proposed scheme is of order 1 when t far away from initial time and of order $\alpha$ when t is near initial time.

5. Conclusions

In this paper, we present a fully discrete linearized numerical scheme for solving the nonlinear fractional differential equation. The Caputo fractional derivative is discretized by the G–L method. The spatial discretization is achieved by using the standard second-order finite difference scheme. With the help of the discrete fractional Grönwall inequality, the optimal error estimates of the scheme are rigorously analyzed. Finally, several numerical examples are shown to confirm our theoretical results.