Crank-Nicolson Quasi-Compact Scheme for the Nonlinear Two-Sided Spatial Fractional Advection-Diffusion Equations

Dechao Gao, Zeshan Qiu^*, Lizan Wang, Jianxin Li

Department of Basic Education, Xinjiang University of Political Science and Law, Tumushuke, China.

DOI: 10.4236/jamp.2024.124068                

Abstract

The higher-order numerical scheme of nonlinear advection-diffusion equations is studied in this article, where the space fractional derivatives are evaluated by using weighted and shifted Grünwald difference operators and combining the compact technique, in the time direction is discretized by the Crank-Nicolson method. Through the energy method, the stability and convergence of the numerical scheme in the sense of L₂-norm are proved, and the convergence order is . Some examples are given to show that our numerical scheme is effective.

Keywords

Crank-Nicolson Quasi-Compact Scheme, Fractional Advection-Diffusion Equations, Nonlinear, Stability and Convergence

1. Introduction

In the last two decades, many fractional differential models have been studied, and these models [1] - [8] have been widely applied in many fields of science and technology, and a lot of research results have been obtained. Among them, Meerschaert and Tadjeran [4] proposed the shifted Grünwald difference operator, which combined with the Crank-Nicolson method to derive the numerical scheme of the advection-dispersion equation. The scheme is unconditionally stable and convergent with order $O\left({\tau }^2+h\right)$ . Later, Meerschaert and Tadjeran [5] utilized the shifted Grünwald difference operators to approximate the left and right Riemann-Liouville fractional derivatives, in time direction is discretized by the explicit Euler method, numerical scheme is unconditionally stable and convergent with order $O\left(\tau +h\right)$ . Tadjeran et al. [9] adopted Crank-Nicolson method to discrete time partial derivative, and one-order shifted Grünwald difference operator to approximate space fractional derivative, derived numerical scheme with convergence order $O\left({\tau }^2+h\right)$ for solving the one-sided space fractional diffusion equation with variable coefficients. By weighting the shifted Grünwald difference operator, Tian et al. [10] combined Crank-Nicolson time discretization to construct a class of second-order numerical scheme for solving the two-sided space fractional diffusion equations. The schemes are proven to achieve convergence accuracy $O\left({\tau }^2+h^2\right)$ . Zhou et al. [11] used the compact technique on the basis of weighted and shifted Grünwald difference operator and combined Crank-Nicolson time discretization to derive quasi-compact scheme for solving the two-sided space fractional diffusion equations, which is proved to be unconditionally stable and convergence with order $O\left({\tau }^2+h^3\right)$ . Hao et al. [12] further derived numerical scheme with convergence accuracy $O\left({\tau }^2+h^4\right)$ by using the compact technique. Haghi et al. [7] proposed a high-order compact numerical scheme for solving the two-dimensional nonlinear time-fractional fourth-order reaction-diffusion equation, the unique solvability of the numerical method is proved in detail. In addition, Li and Deng [13] proposed tempered weighted and shifted Grünwald difference operators for the Riemann-Liouville tempered fractional derivatives, and then a class of second-order numerical schemes is proposed for solving two-sided space tempered fractional diffusion equation. Numerical schemes show convergence order is $O\left({\tau }^2+h^2\right)$ . More work on tempered fractional models reference [14] - [23] .

In this paper, the following nonlinear fractional advection-diffusion equations are considered:

$(\begin{array}{ll}\frac{\partial u\left(x,t\right)}{\partial t}=-\kappa \frac{\partial u\left(x,t\right)}{\partial x}+l_a{D}_x^{\alpha }u\left(x,t\right)+r_x{D}_b^{\alpha }u\left(x,t\right)+f\left(u,x,t\right),\hfill & \left(x,t\right)\in \left(a,b\right)\times \left(0,T\right],\hfill \\ u\left(x,0\right)=\phi \left(x\right),\hfill & x\in \left[a,b\right],\hfill \\ u\left(a,t\right)={\psi }_l\left(t\right),u\left(b,t\right)={\psi }_r\left(t\right),\hfill & t\in \left[0,T\right],\hfill \end{array}$ (1)

where $1<\alpha <2$ , $\kappa$ is the mean advective velocity, non-negative constants l and r denote the diffusion coefficients, which satisfy that $l+r\ne 0$ , ${\psi }_l\left(t\right)\equiv 0$ if $l\ne 0$ , and ${\psi }_r\left(t\right)\equiv 0$ if $r\ne 0$ , the nonlinear source term $f\left(u\mathrm{,}x\mathrm{,}t\right)$ satisfies Lipschitz condition:

$\left|f\left(u\mathrm{,}x\mathrm{,}t\right)-f\left(v\mathrm{,}x\mathrm{,}t\right)\right|\le L\left|u-v\right|\mathrm{,}\forall u\mathrm{,}v\in ℝ\mathrm{,}$ (2)

${}_{a}D{}_x^{\alpha }u\left(x\mathrm{,}t\right)$ and ${}_{x}D{}_b^{\alpha }u\left(x\mathrm{,}t\right)$ represent the left and right Riemann-Liouville fractional derivatives respectively, which are defined as [24] :

${}_{a}D{}_x^{\alpha }u\left(x\mathrm{,}t\right)=\frac{1}{\Gamma \left(2-\alpha \right)}\frac{{\partial }^2}{\partial x^2}\left({\displaystyle {\int }_a^x}\frac{u\left(\tau \mathrm{,}t\right)}{{\left(\tau -x\right)}^{\alpha -1}}\text{d}\tau \right)\mathrm{,}$

${}_{x}D{}_b^{\alpha }u\left(x\mathrm{,}t\right)\mathrm{<mo>\; =\; </mo>}\frac{1}{\Gamma \left(2-\alpha \right)}\frac{{\partial }^2}{\partial x^2}\left({\displaystyle {\int }_x^b}\frac{u\left(\tau \mathrm{,}t\right)}{{\left(\tau -x\right)}^{\alpha -1}}\text{d}\tau \right)\mathrm{.}$

The sections of this article are set as follows. The quasi-compact difference approximations of fractional derivatives are introduced in Section 2. The derivation process of numerical scheme for solving problem (1) is given in Section 3. Section 4 gives a detailed proof of the stability and convergence of the numerical scheme. In Section 5, numerical experiments are given to show that the numerical scheme constructed is effective. Section 6 gives a brief summary of the work.

2. Quasi-Compact Difference Approximations for the Fractional Derivatives

$S^{n+\alpha }\left(ℝ\right)=\left\{\nu |\nu \in L_1\left(ℝ\right)\text{and}{\displaystyle {\int }_{ℝ}}\left(1+\left|w\right|\right){|}^{n+\alpha }\left|\hat{\nu }\left(w\right)\right|\text{d}w<\infty \right\},$

is a fractional Sobolev space $S^{n+\alpha }\left(ℝ\right)$ , where $\hat{\nu }\left(w\right)={\displaystyle {\int }_{ℝ}}\nu \left(x\right){\text{e}}^{-iwx}\text{d}x$ is the Fourier transform of $\nu \left(x\right)$ .

Lemma 2.1. [9] Let $\nu \left(x\right)\in S^{n+\alpha }\left(ℝ\right)$ , $1<\alpha <2$ , the shift number p is an integer. The shifted Grünwald difference operators are defined as:

$A_{h,p}^{\alpha }\nu \left(x\right)=\frac{1}{h^{\alpha }}{\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}}{g}_k^{\alpha }\nu \left(x-\left(k-p\right)h\right),$ (3)

${\hat{A}}_{h,p}^{\alpha }\nu \left(x\right)=\frac{1}{h^{\alpha }}{\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}}{g}_k^{\alpha }\nu \left(x+\left(k-p\right)h\right),$ (4)

then

$A_{h,p}^{\alpha }\nu \left(x\right)={}_{-\infty }D{}_x^{\alpha }\nu \left(x\right)+{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}{a}_p^{\alpha ,k}{}_{-\infty }D{}_x^{\alpha +k}\nu \left(x\right)h^k+O\left(h^n\right),$ (5)

${\hat{A}}_{h,p}^{\alpha }\nu \left(x\right)={}_{x}D{}_{+\infty }^{\alpha }\nu \left(x\right)+{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}{a}_p^{\alpha ,k}{}_{x}D{}_{+\infty }^{\alpha +k}\nu \left(x\right)h^k+O\left(h^n\right),$ (6)

where $g_k^{\alpha }={\left(-1\right)}^k\left(\begin{array}{c}\alpha \\ k\end{array}\right)$ ( $k\ge 0$ ) denotes the normalized Grünwald weights, $a_p^{\alpha \mathrm{,}k}$ are the power series expansion coefficients of the function $W_p\left(s\right)={\text{e}}^{ps}{\left(\frac{1-{\text{e}}^{-s}}{s}\right)}^{\alpha }=1+\left(p-\frac{\alpha }{2}\right)s+\left(\frac{p^2}{2}-\frac{p\alpha }{2}+\frac{\alpha }{6}+\frac{\alpha \left(\alpha -1\right)}{8}\right)s^2+O\left({\left|s\right|}^3\right)$ .

Lemma 2.2. [25] Let $\nu \left(x\right)\in S^{3+\alpha }\left(ℝ\right)$ , $1<\alpha <2$ , if two difference operators are defined as:

$\begin{array}{c}{B}_h^{\alpha }\nu \left(x\right)={\gamma }_1{A}_{h,-1}^{\alpha }\nu \left(x\right)+{\gamma }_2{A}_{h,0}^{\alpha }\nu \left(x\right)+{\gamma }_3{A}_{h,1}^{\alpha }\nu \left(x\right)\\ =\frac{1}{h^{\alpha }}{\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}}{w}_k^{\alpha }\nu \left(x-\left(k-1\right)h\right),\end{array}$ (7)

$\begin{array}{c}{\hat{B}}_h^{\alpha }\nu \left(x\right)={\gamma }_1{\hat{A}}_{h,-1}^{\alpha }\nu \left(x\right)+{\gamma }_2{\hat{A}}_{h,0}^{\alpha }\nu \left(x\right)+{\gamma }_3{\hat{A}}_{h,1}^{\alpha }\nu \left(x\right)\\ =\frac{1}{h^{\alpha }}{\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}}{w}_k^{\alpha }\nu \left(x+\left(k-1\right)h\right).\end{array}$ (8)

then

${\left(I+\frac{1}{6}{h}^2{}_{-\infty }{D}_x^2\right)}_{-\infty }{D}_x^{\alpha }\nu \left(x\right)=B_h^{\alpha }\nu \left(x\right)+O\left(h^3\right)\mathrm{,}$ (9)

${\left(I+\frac{1}{6}{h}^2{}_x{D}_{+\infty }^2\right)}_x{D}_{+\infty }^{\alpha }\nu \left(x\right)={\hat{B}}_h^{\alpha }\nu \left(x\right)+O\left(h^3\right)\mathrm{.}$ (10)

where $w_k^{\alpha }={\gamma }_1{g}_{k-2}^{\alpha }+{\gamma }_2{g}_{k-1}^{\alpha }+{\gamma }_3{g}_k^{\alpha }\left(g_{-2}^{\alpha }=g_{-1}^{\alpha }=0\right)$ , ${\gamma }_1=\frac{1}{24}\left(3{\alpha }^2-7\alpha +4\right)$ , ${\gamma }_2=\frac{1}{12}\left(-3{\alpha }^2+\alpha +8\right)$ , ${\gamma }_3=\frac{1}{24}\left(3{\alpha }^2+5\alpha +4\right)$ .

Noticed, ${}_{-\infty }D{}_x^{2+\alpha }\nu \left(x\right)=\frac{{\text{d}}^2}{\text{d}{x}^2}\left({}_{-\infty }D{}_x^{\alpha }\nu \left(x\right)\right)$ , ${}_{x}D{}_{+\infty }^{2+\alpha }\nu \left(x\right)=\frac{{\text{d}}^2}{\text{d}{x}^2}\left({}_{x}D{}_{+\infty }^{\alpha }\nu \left(x\right)\right)$ [24] , specially,

$\begin{array}{c}\left(I+\frac{1}{6}{h}^2\frac{{\text{d}}^2}{\text{d}{x}^2}\right)\frac{\text{d}\nu \left(x\right)}{\text{d}x}=\frac{1}{2h}\left(\nu \left(x+h\right)-\nu \left(x-h\right)\right)+O\left(h^3\right)\\ =B_h^1\nu \left(x\right)+O\left(h^3\right)\mathrm{,}\end{array}$ (11)

$\begin{array}{c}\left(I+\frac{1}{6}{h}^2\frac{{\text{d}}^2}{\text{d}{x}^2}\right)\nu \left(x\right)=\frac{1}{6}\left(\nu \left(x+h\right)+4\nu \left(x\right)+\nu \left(x-h\right)\right)\mathrm{<mo\; stretchy="false">\; )\; </mo>}+O\left(h^4\right)\\ =C_h\nu \left(x\right)+O\left(h^4\right)\mathrm{.}\end{array}$ (12)

3. Numerical Scheme

The time interval $\left[\mathrm{0,}T\right]$ and the space interval $\left[a\mathrm{,}b\right]$ are divided into equidistant grids, and the time stepsize is denoted as $\tau =T/N$ and the space stepsize is denoted as $h=\left(b-a\right)/M$ respectively, $t_n=n\tau$ , $0\le n\le N$ , $x_i=a+ih$ , $0\le i\le M$ .

Denoting $t_{n+1/2}=\frac{t_n+t_{n+1}}{2}$ , ${\delta }_t{u}_i^{n+1/2}=\frac{u_i^{n+1}-u_i^n}{\tau }$ , $u_i^{n+1/2}=\frac{u_i^n+u_i^{n+1}}{2}$ , $f_i^{n+1/2}=f\left(x_i\mathrm{,}{t}_{n+1/2}\mathrm{,}\frac{u_i^n+u_i^{n+1}}{2}\right)$ , ${\hat{f}}_i^{n+1/2}=f\left(x_i\mathrm{,}{t}_{n+1/2}\mathrm{,}\frac{U_i^n+U_i^{n+1}}{2}\right)$ .

The function $u\left(x\mathrm{,}\cdot \right)$ in problem (1) belongs to $S^{3+\alpha }\left(ℝ\right)$ after zero extension.

Applying $\left(I+\frac{1}{6}{h}^2\frac{{\text{d}}^2}{\text{d}{x}^2}\right)$ to both sides of the equation in problem (1), andthen discretizing time partial derivative at point $\left(x_i\mathrm{,}{t}_{n+1/2}\right)$ by Crank-Nicolson method, and from Lemma 2.2, we can see:

$C_h{\delta }_t{u}_i^{n+1/2}=-\kappa B_h^1{u}_i^{n+1/2}+l{B}_h^{\alpha }{u}_i^{n+1/2}+r{\hat{B}}_h^{\alpha }{u}_i^{n+1/2}+C_h{f}_i^{n+1/2}+R_i^n\mathrm{,0}\le n\le N-\mathrm{1,}$ (13)

$C_h\frac{u_i^{n+1}-u_i^n}{\tau }=-\kappa B_h^1{u}_i^{n+1/2}+l{B}_h^{\alpha }{u}_i^{n+1/2}+r{\hat{B}}_h^{\alpha }{u}_i^{n+1/2}+C_h{f}_i^{n+1/2}+R_i^n\mathrm{,}$ (14)

where $R_i^n=O\left({\tau }^2+h^3\right)$ is the local truncation error.

Denoting $U_i^{n+1/2}=\frac{U_i^n+U_i^{n+1}}{2}$ , the corresponding numerical scheme is obtained by eliminating the local truncation error in (14):

$C_h\frac{U_i^{n+1}-U_i^n}{\tau }=-\kappa B_h^1{U}_i^{n+1/2}+l{B}_h^{\alpha }{U}_i^{n+1/2}+r{\hat{B}}_h^{\alpha }{U}_i^{n+1/2}+C_h{\hat{f}}_i^{n+1/2}\mathrm{,}1\le i\le M-1.$ (15)

The corresponding matrix form is:

$C\left(U^{n+1}-U^n\right)=P{U}^{n+1/2}+\tau C{\hat{f}}^{n+1/2}+\tau F^n\mathrm{,}$ (16)

where $U^n={\left(U_1^n\mathrm{,}{U}_2^n\mathrm{,}\cdots \mathrm{,}{U}_{M-1}^n\right)}^{\text{T}}$ , ${\hat{f}}^{n+1/2}={\left({\hat{f}}_1^{n+1/2}\mathrm{,}{\hat{f}}_2^{n+1/2}\mathrm{,}\cdots \mathrm{,}{\hat{f}}_{M-1}^{n+1/2}\right)}^{\text{T}}$ , $P=\tau \left(-\kappa B^1+l{B}^{\alpha }+r{B}^{\alpha }{}^T\right)$ , $C=tridiag\left\{\frac{1}{6}\mathrm{,}\frac{2}{3}\mathrm{,}\frac{1}{6}\right\}$ , $B^1=\frac{1}{2h}tridiag\left\{-\mathrm{1,0,1}\right\}$ ,

$B^{\alpha }=\frac{1}{h^{\alpha }}\left(\begin{array}{cccccc}{w}_1^{\alpha }& w_0^{\alpha }& \mathrm{<mi>\; </mi>}& \mathrm{<mi>\; </mi>}& \mathrm{<mi>\; </mi>}& \mathrm{<mi>\; </mi>}\\ w_2^{\alpha }& w_1^{\alpha }& w_0^{\alpha }& \mathrm{<mi>\; </mi>}& \mathrm{<mi>\; </mi>}& \mathrm{<mi>\; </mi>}\\ ⋮& ⋮& \ddots & \ddots & \mathrm{<mi>\; </mi>}& \mathrm{<mi>\; </mi>}\\ ⋮& ⋮& ⋮& \ddots & \ddots & \mathrm{<mi>\; </mi>}\\ w_{M-2}^{\alpha }& w_{M-3}^{\alpha }& w_{M-4}^{\alpha }& \dots & w_1^{\alpha }& w_0^{\alpha }\\ w_{M-1}^{\alpha }& w_{M-2}^{\alpha }& w_{M-3}^{\alpha }& \dots & w_2^{\alpha }& w_1^{\alpha }\end{array}\right)\mathrm{,}$ (17)

$\begin{array}{c}{F}^n=\frac{1}{6\tau }{\left(U_0^n-U_0^{n+1},0,\cdots ,0,U_M^n-U_M^{n+1}\right)}^{\text{T}}+\frac{1}{6}{\left({\hat{f}}_0^{n+1/2},0,\cdots ,0,{\hat{f}}_M^{n+1/2}\right)}^{\text{T}}\\ +\frac{\kappa }{2h}{\left(U_0^{n+1/2},0,\cdots ,0,-U_M^{n+1/2}\right)}^{\text{T}}+\frac{U_0^{n+1/2}}{h^{\alpha }}{\left(l{w}_2^{\alpha }+r{w}_0^{\alpha },l{w}_3^{\alpha },\cdots ,l{w}_{M-1}^{\alpha },l{w}_M^{\alpha }\right)}^{\text{T}}\\ +\frac{U_M^{n+1/2}}{h^{\alpha }}{\left(r{w}_M^{\alpha },r{w}_{M-1}^{\alpha },\cdots ,r{w}_3^{\alpha },r{w}_2^{\alpha }+l{w}_0^{\alpha }\right)}^{\text{T}}.\end{array}$

4. Theoretical Analysis of the Numerical Scheme

Before giving a detailed theoretical analysis, let’s first give some lemmas that will be used.

Lemma 4.1. [26] Real matrix A of order n is negative definite if and only if $D=\frac{A+A^{\text{T}}}{2}$ is negative definite.

Lemma 4.2. $C=tridiag\left\{\frac{1}{6}\mathrm{,}\frac{2}{3}\mathrm{,}\frac{1}{6}\right\}$ , for all given nonzero real column vectors $\epsilon$ , C satisfies that:

$\frac{1}{3}{\epsilon }^{\text{T}}\epsilon \le {\epsilon }^{\text{T}}C\epsilon \le {\epsilon }^{\text{T}}\epsilon \mathrm{.}$

Proof. By Gerschgorin disk theorem [27] , it is easy to check $\frac{1}{3}\le \lambda \left(C\right)\le 1$ , then $\frac{1}{3}{\epsilon }^{\text{T}}\epsilon \le {\epsilon }^{\text{T}}C\epsilon \le {\epsilon }^{\text{T}}\epsilon$ .¨

Lemma 4.3. Matrix P is given in (16) is negative definite.

Proof. Let $D=\frac{P+P^{\text{T}}}{2}=\frac{\left(l+r\right)\tau }{2}\left(B^{\alpha }+B^{\alpha }{}^{\text{T}}\right)$ , from [17] , $B^{\alpha }$ is negative definite, by Lemma 4.1, it is known that the matrix P is negative definite.¨

Define $U_h=\{u\mathrm{<mo>\; |\; </mo>}u=\left\{u_i\right\}$ is a grid function defined on ${\left\{x_i=a+ih\right\}}_{i=1}^{M-1}\}$ . For $u\in U_h$ , the corresponding discrete $L_2$ -norm is defined as ${\Vert u\Vert }_{L_2}={\left(h{\displaystyle \underset{i=1}{\overset{M-1}{\sum }}}{u}_i^2\right)}^{1/2}$ .

Theorem 4.1. The numerical scheme (15) is stable.

Proof. Let $U_i^n$ and $V_i^n$ represent the numerical solutions obtained by scheme (15) under different initial conditions respectively. Denoting ${\epsilon }_i^n=U_i^n-V_i^n$ , ${\epsilon }^n={\left({\epsilon }_1^n\mathrm{,}{\epsilon }_2^n\mathrm{,}\cdots \mathrm{,}{\epsilon }_{M-1}^n\right)}^{\text{T}}$ , then

$C\left({\epsilon }^{n+1}-{\epsilon }^n\right)=P{\epsilon }^{n+1/2}+\tau C\left({\hat{f}}^{n+1/2}-{\bar{f}}^{n+1/2}\right),$ (18)

where ${\bar{f}}_i^{n+1/2}=f\left(x_i\mathrm{,}{t}_{n+1/2}\mathrm{,}\frac{V_i^n+V_i^{n+1}}{2}\right)$ .

Let’s multiply both sides of (18) by $h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}$ , we get:

$h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}C\left({\epsilon }^{n+1}-{\epsilon }^n\right)=h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}P{\epsilon }^{n+1/2}+\tau h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}C\left({\hat{f}}^{n+1/2}-{\bar{f}}^{n+1/2}\right),$ (19)

because

$\begin{array}{c}h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}C\left({\epsilon }^{n+1}-{\epsilon }^n\right)=\frac{h}{2}{\left({\epsilon }^{n+1}+{\epsilon }^n\right)}^{\text{T}}C\left({\epsilon }^{n+1}-{\epsilon }^n\right)\\ =\frac{h}{2}\left({\left({\epsilon }^{n+1}\right)}^{\text{T}}C{\epsilon }^{n+1}-{\left({\epsilon }^n\right)}^{\text{T}}C{\epsilon }^n\right),\end{array}$ (20)

arrange the formula to obtain:

$\begin{array}{l}\frac{h}{2}\left({\left({\epsilon }^{n+1}\right)}^{\text{T}}C{\epsilon }^{n+1}-{\left({\epsilon }^n\right)}^{\text{T}}C{\epsilon }^n\right)\\ =h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}P{\epsilon }^{n+1/2}+\tau h{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}C\left({\hat{f}}^{n+1/2}-{\bar{f}}^{n+1/2}\right).\end{array}$ (21)

The following inequalities are established:

$\left|{\epsilon }_i^{n+1/2}\right|=\left|\frac{1}{2}\left({\epsilon }_i^n+{\epsilon }_i^{n+1}\right)\right|\le \frac{1}{2}\left(\left|{\epsilon }_i^n\right|+\left|{\epsilon }_i^{n+1}\right|\right)\mathrm{,}$

$\left|{\hat{f}}_i^{n+1/2}-{\bar{f}}_i^{n+1/2}\right|\le L\left|{\epsilon }_i^{n+1/2}\right|\le \frac{L}{2}\left(\left|{\epsilon }_i^n\right|+\left|{\epsilon }_i^{n+1}\right|\right)\mathrm{,}$

$\begin{array}{l}{\left({\epsilon }^{n+1/2}\right)}^{\text{T}}C\left({\hat{f}}^{n+1/2}-{\bar{f}}^{n+1/2}\right)\\ \le \frac{L}{4}\left(\left|{\epsilon }_1^n\right|+\left|{\epsilon }_1^{n+1}\right|,\cdots ,\left|{\epsilon }_{M-1}^n\right|+\left|{\epsilon }_{M-1}^{n+1}\right|\right)C{\left(\left|{\epsilon }_1^n\right|+\left|{\epsilon }_1^{n+1}\right|,\cdots ,\left|{\epsilon }_{M-1}^n\right|+\left|{\epsilon }_{M-1}^{n+1}\right|\right)}^{\text{T}}\\ \le \frac{L}{2}(\left(\left|{\epsilon }_1^n\right|,\cdots ,\left|{\epsilon }_{M-1}^n\right|\right)C{\left(\left|{\epsilon }_1^n\right|,\cdots ,\left|{\epsilon }_{M-1}^n\right|\right)}^{\text{T}}\\ +\left(\left|{\epsilon }_1^{n+1}\right|,\cdots ,\left|{\epsilon }_{M-1}^{n+1}\right|\right)C{\left(\left|{\epsilon }_1^{n+1}\right|,\cdots ,\left|{\epsilon }_{M-1}^{n+1}\right|\right)}^{\text{T}})\\ \le \frac{L}{2h}\left({\Vert {\epsilon }^n\Vert }_{L_2}^2+{\Vert {\epsilon }^{n+1}\Vert }_{L_2}^2\right).\end{array}$

Denoting $E^{n+1}=h{\left({\epsilon }^{n+1}\right)}^{\text{T}}C{\epsilon }^{n+1}$ , from Lemma 4.3 and (21), then

$\begin{array}{c}{E}^{n+1}\le E^n+\tau L\left({\Vert {\epsilon }^n\Vert }_{L_2}^2+{\Vert {\epsilon }^{n+1}\Vert }_{L_2}^2\right)\\ \le E^0+\tau L\left({\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert {\epsilon }^k\Vert }_{L_2}^2+{\displaystyle \underset{k=1}{\overset{n+1}{\sum }}}{\Vert {\epsilon }^k\Vert }_{L_2}^2\right)\\ =E^0+\tau L{\Vert {\epsilon }^0\Vert }_{L_2}^2+\tau L{\Vert {\epsilon }^{n+1}\Vert }_{L_2}^2+2\tau L{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\Vert {\epsilon }^k\Vert }_{L_2}^2.\end{array}$ (22)

Applying $\frac{1}{3}{\Vert {\epsilon }^k\Vert }_{L_2}^2\le E^k\le {\Vert {\epsilon }^k\Vert }_{L_2}^2$ , then for all given $\mu \in \left(\mathrm{0,1}\right)$ , when $0<\tau \le \frac{1-\mu }{3L}$ , we obtain:

$\begin{array}{c}{\Vert {\epsilon }^{n+1}\Vert }_{L_2}^2\le \frac{3+3\tau L}{1-3\tau L}{\Vert {\epsilon }^0\Vert }_{L_2}^2+\frac{6\tau L}{1-3\tau L}{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\Vert {\epsilon }^k\Vert }_{L_2}^2\\ \le \frac{4-\mu }{\mu }{\Vert {\epsilon }^0\Vert }_{L_2}^2+\frac{6\tau L}{\mu }{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\Vert {\epsilon }^k\Vert }_{L_2}^2.\end{array}$ (23)

By the discrete Gronwall inequalities,

${\Vert {\epsilon }^{n+1}\Vert }_{L_2}^2\le \frac{4-\mu }{\mu }{e}^{\frac{6n\tau L}{\mu }}{\Vert {\epsilon }^0\Vert }_{L_2}^2\le \frac{4-\mu }{\mu }{e}^{\frac{6TL}{\mu }}{\Vert {\epsilon }^0\Vert }_{L_2}^2\mathrm{.}$ (24)

Therefore,

${\Vert {\epsilon }^{n+1}\Vert }_{L_2}\le C_1{\Vert {\epsilon }^0\Vert }_{L_2}\mathrm{,}0\le n\le N-\mathrm{1,}$ (25)

where $C_1=\sqrt{\frac{4-\mu }{\mu }{\text{e}}^{\frac{6TL}{\mu }}}$ .¨

Theorem 4.2. The numerical scheme (15) is convergent.

${\Vert e^n\Vert }_{L_2}\le C_2\left({\tau }^2+h^3\right)\mathrm{,1}\le n\le N\mathrm{,}$

where $e^n={\left(e_1^n\mathrm{,}{e}_2^n\mathrm{,}\cdots \mathrm{,}{e}_{M-1}^n\right)}^{\text{T}}$ , $e_i^n=u_i^n-U_i^n$ , $C_2$ is existed constant.

Proof. Subtracting (15) from (14), we know that:

$C_h\frac{e_i^{n+1}-e_i^n}{\tau }=-\kappa B_h^1{e}_i^{n+1/2}+l{B}_h^{\alpha }{e}_i^{n+1/2}+r{\hat{B}}_h^{\alpha }{u}_i^{n+1/2}+C_h\left(f_i^{n+1/2}-{\hat{f}}_i^{n+1/2}\right)+R_i^n\mathrm{.}$ (26)

The corresponding matrix form is:

$C\left(e^{n+1}-e^n\right)=P{e}^{n+1/2}+\tau C\left(f^{n+1/2}-{\hat{f}}^{n+1/2}\right)+\tau R^n\mathrm{,}$ (27)

where $R^n={\left(R_1^n\mathrm{,}{R}_2^n\mathrm{,}\cdots \mathrm{,}{R}_{M-1}^n\right)}^{\text{T}}$ .

Let’s multiply both sides of (27) by $h{\left(e^{n+1/2}\right)}^{\text{T}}$ ,

$\begin{array}{l}h{\left(e^{n+1/2}\right)}^{\text{T}}C\left(e^{n+1}-e^n\right)\\ =h{\left(e^{n+1/2}\right)}^{\text{T}}P{e}^{n+1/2}+\tau h{\left(e^{n+1/2}\right)}^{\text{T}}C\left(f^{n+1/2}-{\hat{f}}^{n+1/2}\right)+\tau h{\left(e^{n+1/2}\right)}^{\text{T}}{R}^n\mathrm{.}\end{array}$ (28)

Denoting $E^{n+1}=h{\left(e^{n+1}\right)}^{\text{T}}C{e}^{n+1}$ , similar to the proof of Theorem 4.1, then

$\begin{array}{c}{E}^{n+1}\le E^n+\tau L\left({\Vert e^n\Vert }_{L_2}^2+{\Vert e^{n+1}\Vert }_{L_2}^2\right)+2\tau h{\left(e^{n+1/2}\right)}^{\text{T}}{R}^n\\ =E^n+\tau L\left({\Vert e^n\Vert }_{L_2}^2+{\Vert e^{n+1}\Vert }_{L_2}^2\right)+\tau h{\displaystyle \underset{i=1}{\overset{M-1}{\sum }}}\left(e_i^n+e_i^{n+1}\right)R_i^n\\ \le E^n+\tau L\left({\Vert e^n\Vert }_{L_2}^2+{\Vert e^{n+1}\Vert }_{L_2}^2\right)+\frac{1}{2}\tau h{\displaystyle \underset{i=1}{\overset{M-1}{\sum }}}\left({\left|e_i^n\right|}^2+{\left|e_i^{n+1}\right|}^2+2{\left|R_i^n\right|}^2\right)\\ =E^n+\tau \frac{1+2L}{2}\left({\Vert e^n\Vert }_{L_2}^2+{\Vert e^{n+1}\Vert }_{L_2}^2\right)+\tau {\Vert R^n\Vert }_{L_2}^2\\ \le E^0+\tau \frac{1+2L}{2}\left({\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert e^k\Vert }_{L_2}^2+{\displaystyle \underset{k=1}{\overset{n+1}{\sum }}}{\Vert e^k\Vert }_{L_2}^2\right)+\tau {\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert R^n\Vert }_{L_2}^2\\ =E^0+\tau \frac{1+2L}{2}{\Vert e^0\Vert }_{L_2}^2+\tau \frac{1+2L}{2}{\Vert e^{n+1}\Vert }_{L_2}^2+\left(1+2L\right)\tau {\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\Vert e^k\Vert }_{L_2}^2+\tau {\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert R^n\Vert }_{L_2}^2.\end{array}$ (29)

Applying $\frac{1}{3}{\Vert e^k\Vert }_{L_2}^2\le E^k\le {\Vert e^k\Vert }_{L_2}^2$ , and noticing that ${\Vert e^0\Vert }_{L_2}^2=0$ , then for all given $\mu \in \left(\mathrm{0,1}\right)$ , when $0<\tau \le \frac{2-6\mu }{3+6L}$ , we obtain:

$\begin{array}{c}{\Vert e^{n+1}\Vert }_{L_2}^2\le \frac{6\tau }{2-\left(3+6L\right)\tau }{\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert R^k\Vert }_{L_2}^2+\frac{6\left(1+2L\right)\tau }{2-\left(3+6L\right)\tau }{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\Vert e^k\Vert }_{L_2}^2\\ \le \frac{\tau }{\mu }{\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert R^k\Vert }_{L_2}^2+\frac{\left(1+2L\right)\tau }{\mu }{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\Vert e^k\Vert }_{L_2}^2.\end{array}$ (30)

By the discrete Gronwall inequalities,

${\Vert e^{n+1}\Vert }_{L_2}^2\le \frac{\tau }{\mu }{e}^{\frac{n\left(1+2L\right)\tau }{\mu }}{\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert R^k\Vert }_{L_2}^2\le \frac{\tau }{\mu }{e}^{\frac{\left(1+2L\right)T}{\mu }}{\displaystyle \underset{k=0}{\overset{n}{\sum }}}{\Vert R^k\Vert }_{L_2}^2.$ (31)

${\Vert e^{n+1}\Vert }_{L_2}\le C_2\left({\tau }^2+h^3\right)\mathrm{.}$ (32)

¨

5. Numerical Experiments

In this section, we give some numerical experiments to prove the accuracy and validity of the numerical scheme.

$\text{Order}={\log }_m\left(\frac{{\Vert e\Vert }_{L_2,h}}{{\Vert e\Vert }_{L_2,h/m}}\right),$

is the order of measurement.

Example 5.1. Consider the following nonlinear fractional advection-diffusion equations.

$(\begin{array}{ll}\frac{\partial u\left(x,t\right)}{\partial t}=-\kappa \frac{\partial u\left(x,t\right)}{\partial x}+l_a{D}_x^{\alpha }u\left(x,t\right)+r_x{D}_b^{\alpha }u\left(x,t\right)+f\left(u,x,t\right),\hfill & \left(x,t\right)\in \left(0,1\right)\times \left(0,1\right],\hfill \\ u\left(0,t\right)=0,u\left(1,t\right)=0,\hfill & t\in \left[0,1\right],\hfill \\ u\left(x,0\right)=x^4{\left(1-x\right)}^4,\hfill & x\in \left[0,1\right],\hfill \end{array}$

where $1<\alpha <2$ , and the nonlinear source term is:

$\begin{array}{c}f\left(u,x,t\right)={\left(u\left(x,t\right)\right)}^{\frac{3}{2}}+4\kappa e^{\alpha t}\left(1-2x\right){\left(1-x\right)}^3{x}^3+\alpha e^{\alpha t}{x}^4{\left(1-x\right)}^4\\ -e^{\alpha t}[l{\displaystyle \underset{k=0}{\overset{4}{\sum }}}{\left(-1\right)}^k\left(\begin{array}{c}4\\ k\end{array}\right)\frac{\Gamma \left(5+k\right)}{\Gamma \left(5+k-\alpha \right)}{x}^{4+k-\alpha }\\ +r{\displaystyle \underset{k=0}{\overset{4}{\sum }}}{\left(-1\right)}^k\left(\begin{array}{c}4\\ k\end{array}\right)\frac{\Gamma \left(5+k\right)}{\Gamma \left(5+k-\alpha \right)}{\left(1-x\right)}^{4+k-\alpha }+e^{\frac{\alpha t}{2}}{x}^6{\left(1-x\right)}^6].\end{array}$

The analytical solution is $u\left(x\mathrm{,}t\right)=e^{\alpha t}{x}^4{\left(1-x\right)}^4$ .

Choose different $\alpha$ , the mean advective velocity $\kappa$ , diffusion coefficients l and r, the proposed method is used to solve Example 5.1. Let $\tau =h^{\frac{3}{2}}$ , the error results and measurement order results obtained by the numerical method are displayed in Table 1 and Table 2. From Table 1 and Table 2, we can see that the numerical scheme reaches the third-order precision in the spatial direction, and the results are in complete agreement with the conclusion of theoretical analysis. In Table 3 and Table 4, we choose different time stepsizes and fixed $h=\frac{1}{1000}$ , and obtain the errors and the time measurement order results. From Table 3 and Table 4, we verify the numerical scheme is second-order in time.

6. Conclusion

The novelty of this paper is that the higher-order numerical scheme of the

Table 1. Errors and the corresponding space measurement order results at $t=1$ , $\kappa =1$ , $l=1$ , $r=2$ .

Table 2. Errors and the corresponding space measurement order results at $t=1$ , $\kappa =\frac{1}{10}$ , $l=2$ , $r=5$ .

Table 3. Errors and the corresponding time measurement order results at $t=1$ , $h=\frac{1}{1000}$ , $\kappa =1$ , $l=1$ , $r=2$ .

Table 4. Errors and the corresponding time measurement order results at $t=1$ , $h=\frac{1}{1000}$ , $\kappa =\frac{1}{10}$ , $l=2$ , $r=5$ .

nonlinear fractional advection-diffusion equations is studied, and the spatial convergence accuracy reaches the third order. Firstly, Crank-Nicolson quasi-compact scheme is constructed. Secondly, the numerical scheme is analyzed theoretically by the energy method in the sense of L₂-norm. Finally, the effectiveness of the numerical scheme is verified in numerical experiments.

Acknowledgments

The authors would like to express thanks to the editors and anonymous reviewers for their valuable comments and suggestions. This research is supported by the Xinjiang University of Political Science and Law President Fund (XZSK2022039).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Baeumer, B. and Meerschaert, M. (2010) Tempered Stable Lévy Motion and Transient Super-Diffusion. Journal of Computational and Applied Mathematics, 233, 2438-2448. https://doi.org/10.1016/j.cam.2009.10.027
[2]	Liu, F., Anh, V. and Turner, I. (2004) Numerical Solution of the Space Fractional Fokker-Plank Equation. Journal of Computational and Applied Mathematics, 166, 209-219. https://doi.org/10.1016/j.cam.2003.09.028
[3]	Liu, F., Zhuang, P. and Burrage, K. (2012) Numerical Methods and Analysis for a Class of Fractional Advection-Dispersion Models. Computers & Mathematics with Applications, 64, 2990-3007. https://doi.org/10.1016/j.camwa.2012.01.020
[4]	Meerschaert, M. and Tadjeran, C. (2004) Finite Difference Approximations for Fractional Advection-Dispersion Equations. Journal of Computational and Applied Mathematics, 172, 65-77. https://doi.org/10.1016/j.cam.2004.01.033
[5]	Meerschaert, M. and Tadjeran, C. (2006) Finite Difference Approximations for Two-Sided Space-Fractional Partial Differential Equations. Applied Numerical Mathematics, 56, 80-90. https://doi.org/10.1016/j.apnum.2005.02.008
[6]	Ortigueira, M. (2006) Riesz Potential Operators and Inverses via Fractional Centred Derivatives. International Journal of Mathematics and Mathematical Sciences, 2006, Article ID: 048391. https://doi.org/10.1155/IJMMS/2006/48391
[7]	Haghi, M., Ilati, M. and Dehghan, M. (2023) A Fourth-Order Compact Difference Method for the Nonlinear Time-Fractional Fourth-Order Reaction-Diffusion Equation. Engineering with Computers, 39, 1329-1340. https://doi.org/10.1007/s00366-021-01524-2
[8]	Zamani-Gharaghoshi, H., Dehghan, M. and Abbaszadeh, M. (2024) A Meshless Collocation Method Based on Pascal Polynomial Approximation and Implicit Closest Point Method for Solving Reaction-Diffusion Systems on Surfaces. Engineering with Computers, 40, 313-322. https://doi.org/10.1007/s00366-023-01794-y
[9]	Tadjeran, C., Meerschaert, M. and Scheffler (2006) A Second-Order Accurate Numerical Approximation for the Fractional Diffusion Equation. Journal of Computational Physics, 213, 205-213. https://doi.org/10.1016/j.jcp.2005.08.008
[10]	Tian, W., Zhou, H. and Deng, W. (2015) A Class of Second Order Difference Approximations for Solving Space Fractional Diffusion Equations. Mathematics of Computation, 84, 1703-1727. https://doi.org/10.1090/S0025-5718-2015-02917-2
[11]	Zhou, H., Tian, W. and Deng, W. (2013) Quasi-Compact Finite Difference Schemes for Space Fractional Diffusion Equations. Journal of Scientific Computing, 56, 45-66. https://doi.org/10.1007/s10915-012-9661-0
[12]	Hao, Z., Sun, Z. and Cao, W. (2015) A Fourth-Order Approximation of Fractional Derivatives with Its Applications. Journal of Computational Physics, 281, 787-825. https://doi.org/10.1016/j.jcp.2014.10.053
[13]	Li, C. and Deng, W. (2016) High Order Schemes for the Tempered Fractional Diffusion Equations. Advances in Computational Mathematics, 42, 543-572. https://doi.org/10.1007/s10444-015-9434-z
[14]	Celik, C. and Duman, M. (2017) Finite Element Method for a Symmetric Tempered Fractional Diffusion Equation. Applied Numerical Mathematics, 120, 270-286. https://doi.org/10.1016/j.apnum.2017.05.012
[15]	Dehghan, M., Abbaszadeh, M. and Deng, W. (2017) Fourth-Order Numerical Method for the Space-Time Tempered Fractional Diffusion-Wave Equation. Applied Mathematics Letters, 73, 120-127. https://doi.org/10.1016/j.aml.2017.04.011
[16]	Hu, D. and Cao, X. (2019) The Implicit Midpoint Method for Riesz Tempered Fractional Diffusion Equation with a Nonlinear Source Term. Advances in Difference Equations, 2019, Article No. 66. https://doi.org/10.1186/s13662-019-1990-y
[17]	Yu, Y., Deng, W., Wu, Y. and Wu, J. (2017) Third Order Difference Schemes (without Using Points Outside of the Domain) for One Sided Space Tempered Fractional Partial Differential Equations. Applied Numerical Mathematics, 112, 126-145. https://doi.org/10.1016/j.apnum.2016.10.011
[18]	Zhang, Y., Li, Q. and Ding, H. (2018) High-Order Numerical Approximation Formulas for Riemann-Liouville (Riesz) Tempered Fractional Derivatives: Construction and Application (I). Applied Mathematics and Computation, 329, 432-443. https://doi.org/10.1016/j.amc.2018.02.023
[19]	Abbaszadeh, M. and Dehghan, M. (2021) Fourth-Order Alternating Direction Implicit Difference Scheme to Simulate the Space-Time Riesz Tempered Fractional Diffusion Equation. International Journal of Computer Mathematics, 98, 2137-2155. https://doi.org/10.1080/00207160.2020.1841175
[20]	Deng, W. and Zhang, Z. (2017) Numerical Schemes of the Time Tempered Fractional Feynman-Kac Equation. Computers & Mathematics with Applications, 73, 1063-1076. https://doi.org/10.1016/j.camwa.2016.12.017
[21]	Hanert, E. and Piret, C. (2014) A Chebyshev Pseudospectral Method to Solve the Space-Time Tempered Fractional Diffusion Equation. SIAM Journal of Scientific Computing, 36, A1797-A1812. https://doi.org/10.1137/130927292
[22]	Hu, D. and Cao, X. (2019) A Fourth-Order Compact ADI Scheme for Two-Dimensional Riesz Space Fractional Nonlinear Reaction-Diffusion Equation. International Journal of Computer Mathematics, 97, 1928-1948. https://doi.org/10.1080/00207160.2019.1671587
[23]	Yang, Q., Liu, F. and Turner, I. (2010) Numerical Methods for Fractional Partial Differential Equations with Riesz Space Fractional Derivatives. Applied Mathematical Modelling, 34, 200-218. https://doi.org/10.1016/j.apm.2009.04.006
[24]	Podlubny, I. (1999) Fractional Differential Equations. Academic Press, New York. https://lib.ugent.be/catalog/rug01:002178612
[25]	Yu, Y., Deng, W. and Wu, Y. (2017) High-Order Quasi-Compact Difference Schemes for Fractional Diffusion Equations. Communications in Mathematical Sciences, 15, 1183-1209. https://doi.org/10.4310/CMS.2017.v15.n5.a1
[26]	Bhatia, R. (2009) Positive Definite Matrices. Princeton University Press, Princeton. https://doi.org/10.1515/9781400827787
[27]	Varga, R. (2009) Matrix Iterative Analysis. Springer Science and Business Media, Berlin.