High Order Difference Schemes for Space Riesz Variable-Order Nonlinear Fractional Diffusion Equations

1. Introduction

Fractional differential equations have been successfully applied in many fields [2,3,4,36] and attracted many scientists. And many researchers investigated different numerical methods. Commonly used discretization methods include finite difference method, finite element method, and finite volume method [14,20,34]. In genearal, the discretization algebra systems are solved by certain preconditioned iteration methods; see for instance [18,35].

Since constant-order (CO) fractional differential equations fail to describe some complex transport diffusion processes [6,13], the variable-order fractional derivatives have be proposed and the variable-order fractional differential equations have received more and more attention [7,15,19,28,37,43]. Samko and Ross firstly gave the definition of variable-order (VO) fractional derivatives [32]. Lorenzo and Hartley derived different type of VO fractional calculus [21]. For applications of VO fractional differential equations, we refer readers to [15,37,43].

In recent years, many scholars have been studying approximate schemes for variable-order fractional differential equations. In [19], the relationship between the Riemann-Liouville VO fractional derivative and Grünwald-Letnikov expansion has been developed, and in [45], the first order explicit and implicit Euler methods are used to discretize variable-order fractional differential equations. Du, Alikhanov and Sun proposed a second order approximation scheme for time Caputo VO derivative by using special points on each time interval [10]. By utilizng the relationship between VO and CO fractional derivative, VO-WSGD operators for Riemann-Liouville VO fractional derivatives is derived, theory results suggest that the proposed schemes is greater than or equal to second order accuracy [17]. Here WSGD is the abbreviation of weighted-shifted-Grünwald-difference. Kheirkhah, Hajipour and Baleanu proposed a third-order WSGD formula to Caputo derivative for one-dimensional (1D) and two-dimensional (2D) time-fractional advection-reaction-subdiffusion equations with variable-order $\alpha (x,t)\in (0,1)$ [16]. Wang, She, Lao and Lin proposed second-order FCD and fourth-order WSFCD operators for space Riesz VO nonlinear fractional diffusion equations (SRVONFDEs) [42]. However, the fourth-order WSFCD may cause stability problem if

$$\alpha (x,t)\notin \left(1,-3+\sqrt[3]{48+\sqrt{{\displaystyle \frac{62207}{27}}}}+\sqrt[3]{48-\sqrt{{\displaystyle \frac{62207}{27}}}}\right)\approx (1,1.6516).$$

The aim of this work is to develop new difference schemes for the Riesz VO fractional derivative by 4-point weighted-shifted-Grünwald-difference (4WSGD) operators. Similar to [17,42], by applying the relationship between the Riesz CO and VO fractional derivatives, we obtain VO-4WSGD operators for Riesz VO fractional derivative. Then we employ the VO-4WSGD operators to SRVONFDEs.

In this paper, we first consider the following initial-boundary value (IBV) problem for a 1D SRVONFDE:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle \frac{\partial u(x,t)}{\partial t}=K(x,t)\frac{{\partial }^{\alpha (x,t)}u(x,t)}{\partial {\left|x\right|}^{\alpha (x,t)}}+f(u,x,t),}\hfill & (x,t)\in (a,b)\times (0,T],\hfill \\ u(x,0)=u_0\left(x\right),\hfill & x\in (a,b),\hfill \\ u(a,t)=u_a\left(t\right)=0,u(b,t)=u_b\left(t\right)=0,\hfill & t\in [0,T],\hfill \end{array}\right.\end{array}$$

(1)

where $0<\underline{K_C}\le K_C(x,t)=K(x,t)C_{\alpha (x,t)}\le \overline{K_C}$ for certain positive constants $\underline{K_C}$ and $\overline{K_C}$, $f(u,x,t)$ is a nonlinear source term satisfying the Lipschitz condition

$$|f(u_1,x,t)-f(u_2,x,t)|\le L|u_1-u_2|$$

for a certain positive constant L, ${\displaystyle \frac{{\partial }^{\alpha (x,t)}u(x,t)}{\partial {\left|x\right|}^{\alpha (x,t)}}}$ denotes the $\alpha (x,t)$th order Riesz VO fractional derivative (with regard to x) with $1<\underline{\alpha }\le \alpha (x,t)\le \overline{\alpha }<2$. Here, the Riesz VO fractional derivative is defined by [30,45]

$${\displaystyle \frac{{\partial }^{\alpha (x,t)}u(x,t)}{\partial {\left|x\right|}^{\alpha (x,t)}}}=C_{\alpha (x,t)}\left[{}_a{}^{R}D_x^{\alpha (x,t)}{+}_x^R{D}_b^{\alpha (x,t)}\right]u(x,t),C_{\alpha (x,t)}=-{\displaystyle \frac{1}{2cos\frac{\pi \alpha (x,t)}{2}}},$$

where ${}_a{}^{R}D_x^{\alpha (x,t)}$ and ${}_x{}^{R}D_b^{\alpha (x,t)}$ are the left-sided and right-sided Riemann-Liouville VO fractional derivatives defined by

$${}_a{}^{R}D_x^{\alpha (x,t)}u(x,t)={\displaystyle \frac{1}{\mathsf{\Gamma }(2-\alpha (x,t\left)\right)}}{\left[{\displaystyle \frac{d^2}{d{\xi }^2}}{\int }_a^{\xi }{\displaystyle \frac{u(\eta ,t)}{{(\xi -\eta )}^{\alpha (x,t)-1}}}d\eta \right]}_{\xi =x},x\in (a,b),$$

$${}_x{}^{R}D_b^{\alpha (x,t)}u(x,t)={\displaystyle \frac{{(-1)}^2}{\mathsf{\Gamma }(2-\alpha (x,t\left)\right)}}{\left[{\displaystyle \frac{d^2}{d{\xi }^2}}{\int }_{\xi }^b{\displaystyle \frac{u(\eta ,t)}{{(\eta -\xi )}^{\alpha (x,t)-1}}}d\eta \right]}_{\xi =x},x\in (a,b).$$

The 2D problem will be discussed in Section 5.

The rest of paper is organised as below. In Section 2, one gains VO-4WSGD operators for Riemann-Liouville VO fractional derivatives and analyze their accuracy. In Section 3, we obtain IE-4WSGD schemes for 1D SRVONFDE and derive stability and convergence. In Section 4, two banded preconditioners are proposed for the discrete systems of linear equations. The spectral radius and the condition number of the preconditioned systems are derived. Extension of the IE-4WSGD scheme to 2D SRVONFDE is discussed in Section 5. Numerical experiments are carried out in Section 6 to verify our theoretical results and concluding remarks are given in Section 7.

2. 4WSGD Operators for Riemann-Liouville VO Fractional Derivatives

Denote $\alpha >0$ and $p\in {\mathbb{N}}^{*}$, the shifted Grünword operator for the left-sided fractional derivative ${}_{-\infty }{}^{R}D_x^{\alpha }u\left(x\right)$ is

$${}_{L}A_{h,p}^{\alpha }u\left(x\right)={\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=0}^{\infty }{g}_k^{\left(\alpha \right)}u(x-(k-p)h),$$

where $g_k^{\left(\alpha \right)}={(-1)}^{k}inom\alpha k$, i.e.,

$$g_0^{\left(\alpha \right)}=1,g_k^{\left(\alpha \right)}={\displaystyle \frac{{(-1)}^k}{k!}}\prod _{j=0}^{k-1}(\alpha -j),k=1,2,\dots .$$

(2)

As a matter of fact, $g_k^{\left(\alpha \right)}$ satisfies

$${(1-z)}^{\alpha }=\sum _{k=0}^{\infty }{g}_k^{\left(\alpha \right)}{z}^k,\left|z\right|<1.$$

To approximate the left-sided Riemann-Liouville fractional derivative, the ${}_{L}A_{h,p}^{\alpha }$ has first-order accuracy:

$${}_{L}A_{h,p}^{\alpha }u\left(x\right){=}_{-\infty }^R{D}_x^{\alpha }u\left(x\right)+\mathcal{O}\left(h\right);$$

(3)

see [22]. Similarly, for the right-sided Riemann-Liouville fractional derivative, the shifted Grünword operator is

$${}_{R}A_{h,p}^{\alpha }u\left(x\right)={\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=0}^{\infty }{g}_k^{\left(\alpha \right)}u(x+(k-p)h)$$

and

$${}_{R}A_{h,p}^{\alpha }u\left(x\right)={}_x{}^{R}D_{\infty }^{\alpha }u\left(x\right)+O\left(h\right).$$

In [34], a class of third-order 4WSGD operators are derived. ${\gamma }_p(p=-1,0,1,2)$ represents the weights, and satisfies

$$\left\{\begin{array}{c}{\gamma }_2=\gamma ,\hfill \\ {\gamma }_1={\displaystyle \frac{3{\alpha }^2+5\alpha }{24}}-3\gamma ,\hfill \\ {\gamma }_0=1+3\gamma +{\displaystyle \frac{\alpha -3{\alpha }^2}{12}},\hfill \\ {\gamma }_{-1}={\displaystyle \frac{3{\alpha }^2-7\alpha }{24}}-\gamma ,\hfill \end{array}\right.$$

(4)

where $\gamma$ is a free parameter. The left 4WSGD operator is defined by

$$\begin{array}{ccc}\hfill {}_{L}D_{h,\gamma }^{\alpha }u\left(x\right)& =& {\gamma }_2·{}_{L}A_{h,2}^{\alpha }u\left(x\right)+{\gamma }_1·{}_{L}A_{h,1}^{\alpha }u\left(x\right)+{\gamma }_0·{}_{L}A_{h,0}^{\alpha }u\left(x\right)+{\gamma }_{-1}·{}_{L}A_{h,-1}^{\alpha }u\left(x\right)\hfill \\ & & ={\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left(\alpha \right)}u(x-(k-2)h),\hfill \end{array}$$

(5)

where

$${\omega }_{k,\gamma }^{\left(\alpha \right)}={\gamma }_2{g}_k^{\left(\alpha \right)}+{\gamma }_1{g}_{k-1}^{\left(\alpha \right)}+{\gamma }_0{g}_{k-2}^{\left(\alpha \right)}+{\gamma }_{-1}{g}_{k-3}^{\left(\alpha \right)}$$

with $g_{-1}^{\left(\alpha \right)}=g_{-2}^{\left(\alpha \right)}=g_{-3}^{\left(\alpha \right)}=0$ and $g_k^{\left(\alpha \right)}$ for $k\ge 0$ being defined by (2). Similarly, the right 4WSGD operator is defined by

$$\begin{array}{c}\hfill {}_{R}D_{h,\gamma }^{\alpha }u\left(x\right)={\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left(\alpha \right)}u(x+(k-2)h).\end{array}$$

(6)

The following theorem indicates that if the weights are given by (4) and $u\left(x\right)$ is sufficiently smooth, then the the corresponding 4WSGD scheme has at least third-order accuracy.

Theorem 1.

[34] Assume $u\left(x\right)\in L^1\left(\mathbb{R}\right)\cap C^2\left(\mathbb{R}\right)$, and ${}_{-\infty }{}^{R}D_x^{\alpha +4}u\left(x\right)$, ${}_x{}^{R}D_{\infty }^{\alpha +4}u\left(x\right),$$\mathcal{F}\left[{}_{-\infty }D_x^{\alpha +4}u\left(x\right)\right]\left(\mathbb{R}\right),$ and $\mathcal{F}\left[{}_x{}^{R}D_{\infty }^{\alpha +4}u\left(x\right)\right]\left(\mathbb{R}\right)\in L^1\left(\mathbb{R}\right)$, and ${\gamma }_{-1}$, ${\gamma }_0$, ${\gamma }_1$, ${\gamma }_2$ are given by (4), then

$$\begin{array}{c}\hfill {}_{-\infty }{}^{R}D_x^{\alpha }u\left(x\right)={}_{L}D_{h,\gamma }^{\alpha }u\left(x\right)+{\eta }_L,{}_x{}^{R}D_{\infty }^{\alpha }u\left(x\right)={}_{R}D_{h,\gamma }^{\alpha }u\left(x\right)+{\eta }_R,\end{array}$$

where

$$|{\eta }_L|\le C_L|\gamma -{\displaystyle \frac{{\alpha }^3-{\alpha }^2-4\alpha }{48}}|h^3+O\left(h^4\right),\left|{\eta }_R\right|\le C_R|\gamma -{\displaystyle \frac{{\alpha }^3-{\alpha }^2-4\alpha }{48}}|h^3+O\left(h^4\right)$$

with $C_L$ and $C_R$ being positive constants.

Remark 1.

From the inequalities for $|{\eta }_L|$ and $|{\eta }_R|$, we know that the 4WSGD operators corresponding to $\gamma ={\displaystyle \frac{{\alpha }^3-{\alpha }^2-4\alpha }{48}}$ has fourth order accuracy. But the 4WSGD operator with 4th order may not suitable for variable-order fractional differential equations since it often causes stability problem.

Now we present the 4WSGD operator for the Riemann-Liouville VO fractional derivatives. According to the concept of Riemann-Liouville VO fractional derivatives, for each $\alpha (x_n,t_m)$, we can rewrite the definition as

$$\begin{array}{ccc}\hfill {}_a{}^{R}D_x^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}& =& {\left[{}_a{}^{R}D_{\xi }^{\alpha (x_n,t_m)}u(\xi ,t_m)\right]}_{\xi =x_n},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {}_x{}^{R}D_b^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}& =& {\left[{}_{\xi }{}^{R}D_b^{\alpha (x_n,t_m)}u(\xi ,t_m)\right]}_{\xi =x_n}.\hfill \end{array}$$

(8)

The formulas (7)–(8) inspire an idea to develop 4WSGD operators for the Riemann-Liouville VO fractional derivatives. From (5) and (7), we get 4WSGD operators for ${}_a{}^{R}D_x^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}$:

$$\begin{array}{ccc}\hfill {}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}u(x_n,t_m)& =& {\left[{}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}u(\xi ,t_m)\right]}_{\xi =x_n}\hfill \\ & =& {\left[{\displaystyle \frac{1}{h^{\alpha (x_n,t_m)}}}\sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}u(\xi -(k-2)h,t_m)\right]}_{\xi =x_n}\hfill \\ & =& {\displaystyle \frac{1}{h^{\alpha (x_n,t_m)}}}\sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}u(x_n-(k-2)h,t_m),\hfill \end{array}$$

(9)

where ${\omega }_k^{\left(\alpha (x_n,t_m)\right)}$ are given by

$${\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}={\gamma }_2{g}_k^{\left(\alpha (x_n,t_m)\right)}+{\gamma }_1{g}_{k-1}^{\left(\alpha (x_n,t_m)\right)}+{\gamma }_0{g}_{k-2}^{\left(\alpha (x_n,t_m)\right)}+{\gamma }_{-1}{g}_{k-3}^{\left(\alpha (x_n,t_m)\right)},$$

(10)

with

$$\left\{\begin{array}{c}{\gamma }_2=\gamma ,\hfill \\ {\gamma }_1=-3\gamma +{\displaystyle \frac{3\alpha {(x_n,t_m)}^2+5\alpha (x_n,t_m)}{24}},\hfill \\ {\gamma }_0=3\gamma +{\displaystyle \frac{\alpha (x_n,t_m)-3\alpha {(x_n,t_m)}^2}{12}}+1,\hfill \\ {\gamma }_{-1}=-\gamma +{\displaystyle \frac{3\alpha {(x_n,t_m)}^2-7\alpha (x_n,t_m)}{24}},\hfill \end{array}\right.$$

(11)

for $k=0,1,2,\dots$ . Here

$$\begin{array}{ccc}\hfill g_{-1}^{\left(\alpha (x_n,t_m)\right)}& =& g_{-2}^{\left(\alpha (x_n,t_m)\right)}=g_{-3}^{\left(\alpha (x_n,t_m)\right)}=0,g_0^{\left(\alpha (x_n,t_m)\right)}=1,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill g_k^{\left(\alpha (x_n,t_m)\right)}& =& {\displaystyle \frac{{(-1)}^k}{k!}}\prod _{j=0}^{k-1}(\alpha (x_n,t_m)-j),k=1,2,\dots .\hfill \end{array}$$

(13)

We refer readers to [17] for the idea of deriving formula (9). Similarly, 4WSGD operators for ${}_x{}^{R}D_{\infty }^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}$ are

$$\begin{array}{c}\hfill {}_{R}D_{h,\gamma }^{\alpha (x_n,t_m)}u(x_n,t_m)={\displaystyle \frac{1}{h^{\alpha (x_n,t_m)}}}\sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}u(x_n+(k-2)h,t_m).\end{array}$$

(14)

We call the formulas (9) and (14) left and right VO-4WSGD operator, respectively.

The zero extension of $u(x,t)$ to $\mathbb{R}\times [0,T]$:

$$\tilde{u}(x,t)=\left\{\begin{array}{cc}u(x,t),\hfill & (x,t)\in (a,b)\times [0,T],\hfill \\ 0,\hfill & \mathrm{else}.\hfill \end{array}\right.$$

Theorem 2.

Let $1<\underline{\alpha }\le \alpha (x,t)\le \overline{\alpha }<2$ and $u(x,t)\in L^1\left(\mathbb{R}\right)\cap C^2\left(\mathbb{R}\right)$, and ${}_{-\infty }{}^{R}D_x^{\overline{\alpha }+4}\tilde{u}(x,t)$, ${}_x{}^{R}D_{\infty }^{\overline{\alpha }+4}\tilde{u}(x,t)$ and their Fourier transforms belong to $L^1\left(\mathbb{R}\right)$. Let ${}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}\tilde{u}(x_n,t_m)$ and ${}_{R}D_{h,\gamma }^{\alpha (x_n,t_m)}\tilde{u}(x_n,t_m)$ be the VO-4WSGD operators given by (9) and (14), respectively. Then

$$\begin{array}{ccc}\hfill {}_{-\infty }{}^{R}D_x^{\alpha (x_n,t_m)}\tilde{u}(x,t_m){|}_{x=x_n}& =& {}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}\tilde{u}(x_n,t_m)+{\eta }_L(x_n,t_m),\hfill \\ \hfill {}_x{}^{R}D_{\infty }^{\alpha (x_n,t_m)}\tilde{u}(x,t_m){|}_{x=x_n}& =& {}_{R}D_{h,\gamma }^{\alpha (x_n,t_m))}\tilde{u}(x_n,t_m))+{\eta }_R(x_n,t_m)),\hfill \end{array}$$

where

$$|{\eta }_L|\le C_L|\gamma -{\displaystyle \frac{\alpha {(x_n,t_m)}^3-\alpha {(x_n,t_m)}^2-4\alpha (x_n,t_m)}{48}}|h^3+O\left(h^4\right),$$

$$|{\eta }_R|\le C_R|\gamma -{\displaystyle \frac{\alpha {(x_n,t_m)}^3-\alpha {(x_n,t_m)}^2-4\alpha (x_n,t_m)}{48}}|h^3+O\left(h^4\right)$$

with $C_L$ and $C_R$ being positive constants.

Proof. 

Since $\alpha (x,t)\in [\underline{\alpha },\overline{\alpha }]\subset (1,2)$, for given x and t, we have that ${}_{-\infty }{}^{R}D_{\xi }^{\alpha (x_n,t_m)+4}\tilde{u}(\xi ,t_m)$, ${}_x{}^{R}D_{\infty }^{\alpha (x_n,t_m)+4}\tilde{u}(\xi ,t)$ and their Fourier transforms belong to $L^1\left(\mathbb{R}\right)$. Thus, for $\alpha (x_n,t_m)$, one gets by Theorem 1,

$$\begin{array}{c}\hfill {}_{-\infty }{}^{R}D_{\xi }^{\alpha (x_n,t_m)}\tilde{u}(\xi ,t_m)={}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}\tilde{u}(\xi ,t_m)+{\eta }_L(\xi ,t_m).\end{array}$$

In particular, by (7) we have

$$\begin{array}{ccc}\hfill {}_{-\infty }{}^{R}D_x^{\alpha (x_n,t_m)}\tilde{u}(x,t_m){|}_{x=x_n}& =& {\left[{}_{-\infty }{}^{R}D_{\xi }^{\alpha (x_n,t_m)}\tilde{u}(\xi ,t_m)\right]}_{\xi =x_n}\hfill \\ & =& {\left[{}_{L}D_{h,\gamma }^{\alpha (x,t)}\tilde{u}(\xi ,t_m)+{\eta }_L(\xi ,t)\right]}_{\xi =x_n}\hfill \\ & =& {}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}\tilde{u}(x_n,t_m)+{\eta }_L(x_n,t_m),\hfill \end{array}$$

where $|{\eta }_L(x_n,t_m)|\le C_L|\gamma -{\displaystyle \frac{\alpha {(x_n,t_m)}^3-\alpha {(x_n,t_m)}^2-4\alpha (x_n,t_m)}{48}}|h^3+O\left(h^4\right).$

Similarly, the relation between the right-sided Riemann-Liouville VO fractional derivative and the right-sided VO-4WSGD operator is obtained by using Theorem 1 and (8).□

Remark 2.

(i) Let $u(x,t)$ be a function defined in $[a,b]\times [0,T]$ with $u(a,t)=u(b,t)=0$ for $t\in [0,T]$, then ${}_a{}^{R}D_x^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}$ can be approximated by

$$\begin{array}{ccc}\hfill {}_{L}D_{h,\gamma }^{\alpha (x_n,t_m)}u(x_n,t_m)& =& {\displaystyle \frac{1}{h^{\alpha (x_n,t_m)}}}\sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}u(x_n-(k-2)h,t_m)\hfill \\ & =& {\displaystyle \frac{1}{h^{\alpha (x_n,t_m)}}}\sum _{k=0}^{\lfloor (x_n-a)/h\rfloor +2}{\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}u(x_n-(k-2)h,t_m),\hfill \end{array}$$

where $\lfloor x\rfloor$ denotes the largest integer that is not greater than x. Similarly, ${}_x{}^{R}D_b^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}$ can be approximated by

$${}_{R}D_{h,\gamma }^{\alpha (x_n,t_m)}u(x_n,t_m)={\displaystyle \frac{1}{h^{\alpha (x_n,t_m)}}}\sum _{k=0}^{\lfloor (b-x_n)/h\rfloor +2}{\omega }_{k,\gamma }^{\left(\alpha (x_n,t_m)\right)}u(x_n+(k-2)h,t_m).$$

(ii) Typically, for a discretization scheme that has good properties, the free parameter $\gamma$, which dependents on the value of $\alpha (x,t)$, is of great importance.

3. IE-4WSGD Scheme for 1D SRVONFDE

3.1. IE-4WSGD Scheme

In this subsection, we derive IE-4WSGD schemes for the IBV problem of 1D SRVONFDE Eq (1). Let $N,M\in {\mathbb{N}}^{*}$, and let $h=(b-a)/N$ and $\tau =T/M$, which represent the space grid size and the time step length, respectively. And denote

$$x_n=a+nh,n=0,1,2,\cdots N,t_m=m\tau ,m=0,1,2,\cdots M.$$

The following symbols are utilized

$$\begin{array}{ccc}& & u_n^m=u(x_n,t_m),{\alpha }_n^m=\alpha (x_n,t_m),f_{e,n}^m=f(u(x_n,t_m),x_n,t_m),d_n^m={\displaystyle \frac{\tau K_{C,n}^m}{h^{{\alpha }_n^m}}},\hfill \\ & & {{(}_a^R{D}_x^{{\alpha }_n^m}u)}_n^m={}_a{}^{R}D_x^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n},{{(}_x^R{D}_b^{{\alpha }_n^m}u)}_n^m={}_x{}^{R}D_b^{\alpha (x_n,t_m)}u(x,t_m){|}_{x=x_n}.\hfill \end{array}$$

Applying the implicit Euler (IE) formula

$${\displaystyle \frac{\partial u(x_n,t_m)}{\partial t}}={\displaystyle \frac{u(x_n,t_m)-u(x_n,t_{m-1})}{\tau }}+O\left(\tau \right)$$

to discretize ${\displaystyle \frac{\partial u(x,t)}{\partial t}}$, one gains

$$\begin{array}{c}\hfill {\displaystyle \frac{u_n^m-u_n^{m-1}}{\tau }}=K_{C,n}^m\left[{{(}_a^R{D}_x^{{\alpha }_n^m}u)}_n^m+{{(}_x^R{D}_b^{{\alpha }_n^m}u)}_n^m\right]+f_{e,n}^m+O({}_{\tau }).\end{array}$$

Approximating ${{(}_a^R{D}_x^{{\alpha }_n^m}u)}_n^m$ and ${{(}_x^R{D}_b^{{\alpha }_n^m}u)}_n^m$ by the left- and right-4WSGD operators respectively, and approximating $f(u_n^m,x_n,t_m)$ by the Taylor expansion, we get the following systems:

$$\begin{array}{ccc}\hfill u_n^m& =& u_n^{m-1}+{\displaystyle \frac{\tau K_{C,n}^m}{h^{{\alpha }_n^m}}}\left[\sum _{k=0}^{n+2}{\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}{u}_{n-k+2}^m+\sum _{k=0}^{N-n+3}{\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}{u}_{n+k-2}^m\right]+\tau f_{e,n}^{m-1}+\tau R_n^m,\hfill \end{array}$$

where $n=1,2,\cdots ,N-1,m=1,2,\cdots ,M$, $|R_n^m|\le C(\tau +h^l)$ for a certain constant $C>0$ with $l\in \{3,4\}$.

Let the numerical approximation of $u_n^m$ by the IE-4WSGD scheme be denoted by $U_n^m$ and let $f_n^m=f(U_n^m,x_n,t_m)$. By using $u_0^m=u_N^m=0$, $u_{-1}^m=u_{N+1}^m=0$, we obtain the linear equations about $\{U_n^m:n=1,2,\cdots ,N-1,m=1,2,\cdots ,M\}$:

$$\begin{array}{c}\hfill U_n^m-d_n^m\left[\sum _{k=0}^{n+2}{\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}{U}_{n-k+2}^m+\sum _{k=0}^{N-n+3}{\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}{U}_{n+k-2}^m\right]=U_n^{m-1}+\tau f_n^{m-1}.\end{array}$$

(15)

Let

$$\begin{array}{ccc}\hfill W_L^m& =& \left(\begin{array}{ccccccc}{\omega }_{2,\gamma }^{\left({\alpha }_1^m\right)}& {\omega }_{1,\gamma }^{\left({\alpha }_1^m\right)}& {\omega }_{0,\gamma }^{\left({\alpha }_1^m\right)}& 0& \cdots & \cdots & 0\\ {\omega }_{3,\gamma }^{\left({\alpha }_2^m\right)}& {\omega }_{2,\gamma }^{\left({\alpha }_2^m\right)}& {\omega }_{1,\gamma }^{\left({\alpha }_2^m\right)}& {\omega }_{0,\gamma }^{\left({\alpha }_2^m\right)}& 0& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& ⋮& \ddots & \ddots & \ddots & \ddots & 0\\ ⋮& ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ {\omega }_{N-1,\gamma }^{\left({\alpha }_{N-2}^m\right)}& {\omega }_{N-2,\gamma }^{\left({\alpha }_{N-2}^m\right)}& {\omega }_{N-3,\gamma }^{\left({\alpha }_{N-2}^m\right)}& {\omega }_{N-4,\gamma }^{\left({\alpha }_{N-2}^m\right)}& \ddots & {\omega }_{2,\gamma }^{\left({\alpha }_{N-2}^m\right)}& {\omega }_{1,\gamma }^{\left({\alpha }_{N-2}^m\right)}\\ {\omega }_{N,\gamma }^{\left({\alpha }_{N-1}^m\right)}& {\omega }_{N-1,\gamma }^{\left({\alpha }_{N-1}^m\right)}& {\omega }_{N-2,\gamma }^{\left({\alpha }_{N-1}^m\right)}& {\omega }_{N-3,\gamma }^{\left({\alpha }_{N-1}^m\right)}& \cdots & {\omega }_{3,\gamma }^{\left({\alpha }_{N-1}^m\right)}& {\omega }_{2,\gamma }^{\left({\alpha }_{N-1}^m\right)}\end{array}\right),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill W_R^m& =& \left(\begin{array}{ccccccc}{\omega }_{2,\gamma }^{\left({\alpha }_1^m\right)}& {\omega }_{3,\gamma }^{\left({\alpha }_1^m\right)}& {\omega }_{4,\gamma }^{\left({\alpha }_1^m\right)}& \cdots & \cdots & {\omega }_{N-1,\gamma }^{\left({\alpha }_1^m\right)}& {\omega }_{N,\gamma }^{\left({\alpha }_1^m\right)}\\ {\omega }_{1,\gamma }^{\left({\alpha }_2^m\right)}& {\omega }_{2,\gamma }^{\left({\alpha }_2^m\right)}& {\omega }_{3,\gamma }^{\left({\alpha }_2^m\right)}& \cdots & \cdots & {\omega }_{N-2,\gamma }^{\left({\alpha }_2^m\right)}& {\omega }_{N-1,\gamma }^{\left({\alpha }_2^m\right)}\\ {\omega }_{0,\gamma }^{\left({\alpha }_3^m\right)}& {\omega }_{1,\gamma }^{\left({\alpha }_3^m\right)}& {\omega }_{2,\gamma }^{\left({\alpha }_3^m\right)}& \cdots & \cdots & {\omega }_{N-3,\gamma }^{\left({\alpha }_3^m\right)}& {\omega }_{N-2,\gamma }^{\left({\alpha }_3^m\right)}\\ 0& {\omega }_{0,\gamma }^{\left({\alpha }_4^m\right)}& {\omega }_{1,\gamma }^{\left({\alpha }_4^m\right)}& \ddots & \cdots & {\omega }_{N-4,\gamma }^{\left({\alpha }_4^m\right)}& {\omega }_{N-3,\gamma }^{\left({\alpha }_4^m\right)}\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮& ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & {\omega }_{2,{\gamma }_{N-2}^m}^{\left({\alpha }_{N-2}^m\right)}& {\omega }_{3,{\gamma }_{N-2}^m}^{\left({\alpha }_{N-2}^m\right)}\\ 0& \cdots & \cdots & 0& {\omega }_{0,\gamma }^{\left({\alpha }_{N-1}^m\right)}& {\omega }_{1,\gamma }^{\left({\alpha }_{N-1}^m\right)}& {\omega }_{2,\gamma }^{\left({\alpha }_{N-1}^m\right)}\end{array}\right),\hfill \end{array}$$

(17)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\omega }_{0,\gamma }^{\left({\alpha }_j^m\right)}={\gamma }_j^m,{\omega }_{1,\gamma }^{\left({\alpha }_j^m\right)}=-(3+{\alpha }_j^m){\gamma }_j^m+{\displaystyle \frac{3{\left({\alpha }_j^m\right)}^2+5{\alpha }_j^m}{24}},\hfill \\ {\omega }_{2,\gamma }^{\left({\alpha }_j^m\right)}={\displaystyle \frac{{\left({\alpha }_j^m\right)}^2+5{\alpha }_j^m+6}{2}}{\gamma }_j^m-{\displaystyle \frac{{\alpha }_j^m(3{\left({\alpha }_j^m\right)}^2+11{\alpha }_j^m-2)}{24}},\hfill \\ {\omega }_{k,\gamma }^{\left({\alpha }_j^m\right)}={\gamma }_j^m{g}_k^{\left({\alpha }_j^m\right)}+(-3{\gamma }_j^m+{\displaystyle \frac{3{\left({\alpha }_j^m\right)}^2+5{\alpha }_j^m}{24}})g_{k-1}^{\left({\alpha }_j^m\right)}+(3{\gamma }_j^m+{\displaystyle \frac{{\alpha }_j^m-3{\left({\alpha }_j^m\right)}^3}{12}}+1)g_{k-2}^{\left({\alpha }_j^m\right)}\hfill \\ +(-{\gamma }_j^m+{\displaystyle \frac{3{\left({\alpha }_j^m\right)}^2-7{\alpha }_j^m}{24}})g_{k-3}^{\left({\alpha }_j^m\right)},k\ge 3.\hfill \end{array}\end{array}$$

(18)

Denote

$$D^m=\mathrm{diag}\left(d_1^m,d_2^m,\cdots ,d_{N-1}^m\right),U^m={(U_1^m,U_2^m,\dots ,U_{N-1}^m)}^T,F^m={(f_1^m,f_2^m,\dots ,f_{N-1}^{m-1})}^T.$$

We get the matrix form of Eq. (15):

$$A^m{U}^m=U^{m-1}+\tau F^{m-1},m=1,2,\dots ,M,$$

(19)

where

$$\begin{array}{c}\hfill A^m=I-D^m{W}^m,\end{array}$$

(20)

where $W^m=W_L^m+W_R^m.$ Obviously, $W^m={\left[w_{ij}^m\right]}_{i,j=1}^{N-1}$ satisfy

$$w_{ij}^m=\left\{\begin{array}{cc}2{\omega }_{2,\gamma }^{\left({\alpha }_i^m\right)},\hfill & j=i,\hfill \\ {\omega }_{1,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_i^m\right)},\hfill & j=i\pm 1,\hfill \\ {\omega }_{0,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_i^m\right)},\hfill & j=i\pm 2,\hfill \\ {\omega }_{|i-j|+2,\gamma }^{\left({\alpha }_i^m\right)},\hfill & |i-j|>2.\hfill \end{array}\right.$$

(21)

3.2. Stability and Convergence Analysis

In this subsection, we analyze the stability and convergence of the IE-4WSGD scheme (15). We note that $W_L^m$, $W_R^m$ and $W^m$ are not Toeplitz matrices since the variable order $\alpha (x,t)$ depends on the space variable x (and the time variable t). Recalling that $D^m$ are diagonal matrices with nonnegative entries, we put forward conditions for $W^m$ to be diagonally dominant with negative diagonal entries, which is sufficient for the IE-4WSGD scheme to be stable.

Lemma 1.

[27,41] For $n=0,1,2,...,N,m=1,2,...,M$ and ${\alpha }_n^m\in (1,2),$ the coefficient ${\left\{g_k^{\left({\alpha }_n^m\right)}\right\}}_{k\ge 0}$ satisfy

(i)

$g_0^{\left({\alpha }_n^m\right)}=1,g_1^{\left({\alpha }_n^m\right)}=-{\alpha }_n^m,g_k^{\left({\alpha }_n^m\right)}>0k\ge 2;$

(ii)

$\sum _{k=0}^{\infty }{g}_k^{\left({\alpha }_n^m\right)}=0,\sum _{k=0}^q{g}_k^{\left({\alpha }_n^m\right)}<0,q\ge 1;$

(iii)

${\displaystyle \frac{{\alpha }_n^m({\alpha }_n^m-1)(2-{\alpha }_n^m)(3-{\alpha }_n^m)}{180}}{\left({\displaystyle \frac{4}{k}}\right)}^{{\alpha }_n^m+1}<g_k^{\left({\alpha }_n^m\right)}\le {\displaystyle \frac{{\alpha }_n^m({\alpha }_n^m-1)}{2}}{\left({\displaystyle \frac{3}{k+1}}\right)}^{{\alpha }_n^m+1},k\ge 2;$

(iv)

${\displaystyle \frac{({\alpha }_n^m-1)(2-{\alpha }_n^m)(3-{\alpha }_n^m)}{45}}{\left({\displaystyle \frac{4}{k}}\right)}^{{\alpha }_n^m}<\sum _{q=k}^{\infty }{g}_q^{\left({\alpha }_n^m\right)}\le {\displaystyle \frac{({\alpha }_n^m-1)3^{{\alpha }_n^m+1}}{2k^{{\alpha }_n^m}}},k\ge 2;$

(v)

$g_k^{\left({\alpha }_n^m\right)}=O\left(k^{-({\alpha }_n^m+1)}\right),k\ge 2$.

Lemma 2.

Assume $1<{\alpha }_n^m<2$ and ${\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}$, $n=0,1,2,...,N,m=1,2,...,M$, are given by (18). Then we have

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\omega }_{0,\gamma }^{\left({\alpha }_n^m\right)}={\gamma }_n^m,{\omega }_{1,\gamma }^{\left({\alpha }_n^m\right)}=-(3+{\alpha }_n^m){\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2+5{\alpha }_n^m}{24}},\hfill \\ {\omega }_{2,\gamma }^{\left({\alpha }_n^m\right)}={\displaystyle \frac{{\left({\alpha }_n^m\right)}^2+5{\alpha }_n^m+6}{2}}{\gamma }_n^m-{\displaystyle \frac{{\alpha }_n^m(3{\left({\alpha }_n^m\right)}^2+11{\alpha }_n^m-2)}{24}},\hfill \\ {\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}=({\displaystyle \frac{-{\left({\alpha }_n^m\right)}^3-6{\left({\alpha }_n^m\right)}^2-11{\alpha }_n^m-6}{k(k-1)(k-2)}}{\gamma }_n^m\hfill \\ +{\displaystyle \frac{24k^2+(-12{\left({\alpha }_n^m\right)}^2-36{\alpha }_n^m-96)k+(3{\left({\alpha }_n^m\right)}^4+14{\left({\alpha }_n^m\right)}^3+33{\left({\alpha }_n^m\right)}^2+46{\alpha }_n^m+72)}{24(k-1)(k-2)}})g_{k-3}^{\left({\alpha }_n^m\right)},k\ge 3,\hfill \\ \sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}=0.\hfill \end{array}\right.\end{array}$$

(22)

Proof. 

The formulas for ${\omega }_{0,\gamma }^{\left({\alpha }_i^m\right)},{\omega }_{1,\gamma }^{\left({\alpha }_i^m\right)}$, and ${\omega }_{2,\gamma }^{\left({\alpha }_i^m\right)}$ can be easily obtained from (11) and (18).

For $k\ge 3,$ utilizing (13), we get

$$\begin{array}{ccc}\hfill g_k^{\left({\alpha }_n^m\right)}& =& {\displaystyle \frac{k-1-{\alpha }_n^m}{k}}{g}_{k-1}^{\left({\alpha }_n^m\right)}={\displaystyle \frac{(k-1-{\alpha }_n^m)(k-2-{\alpha }_n^m)}{k(k-1)}}{g}_{k-2}^{\left({\alpha }_n^m\right)}\hfill \\ & =& {\displaystyle \frac{(k-1-{\alpha }_n^m)(k-2-{\alpha }_n^m)(k-3-{\alpha }_n^m)}{k(k-1)(k-2)}}{g}_{k-3}^{\left({\alpha }_n^m\right)}.\hfill \end{array}$$

Furthermore, by (10), we get

$$\begin{array}{ccc}& & {\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}\hfill \\ & =& {\gamma }_n^m{g}_k^{\left({\alpha }_n^m\right)}+\left(-3{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2+5{\alpha }_n^m}{24}}\right)g_{k-1}^{\left({\alpha }_n^m\right)}+\left(3{\gamma }_n^m+{\displaystyle \frac{{\alpha }_n^m-3{\left({\alpha }_n^m\right)}^3}{12}}+1\right)g_{k-2}^{\left({\alpha }_n^m\right)}\hfill \\ & & +\left(-{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2-7{\alpha }_n^m}{24}}\right)g_{k-3}^{\left({\alpha }_n^m\right)}\hfill \\ & =& ({\gamma }_n^m{\displaystyle \frac{(k-1-{\alpha }_n^m)(k-2-{\alpha }_n^m)(k-3-{\alpha }_n^m)}{k(k-1)(k-2)}}\hfill \\ & & +\left(-3{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2+5{\alpha }_n^m}{24}}\right){\displaystyle \frac{(k-2-{\alpha }_n^m)(k-3-{\alpha }_n^m)}{(k-1)(k-2)}})g_{k-3}^{\left({\alpha }_n^m\right)}\hfill \\ & & +\left(\left(3{\gamma }_n^m+{\displaystyle \frac{{\alpha }_n^m-3{\left({\alpha }_n^m\right)}^3}{12}}+1\right){\displaystyle \frac{k-3-{\alpha }_n^m}{k-2}}+\left(-{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2-7{\alpha }_n^m}{24}}\right)\right)g_{k-3}^{\left({\alpha }_n^m\right)}\hfill \\ & =& ({\displaystyle \frac{-{\left({\alpha }_n^m\right)}^3-6{\left({\alpha }_n^m\right)}^2-11{\alpha }_n^m-6}{k(k-1)(k-2)}}{\gamma }_n^m\hfill \\ & & +{\displaystyle \frac{24k^2+(-12{\left({\alpha }_n^m\right)}^2-36{\alpha }_n^m-96)k+(3{\left({\alpha }_n^m\right)}^4+14{\left({\alpha }_n^m\right)}^3+33{\left({\alpha }_n^m\right)}^2+46{\alpha }_n^m+72)}{24(k-1)(k-2)}})g_{k-3}^{\left({\alpha }_n^m\right)}.\hfill \end{array}$$

Finally, by Lemma 1 and Eq. (12), we get

$$\begin{array}{ccc}& & \sum _{k=0}^{\infty }{\omega }_{k,\gamma }^{\left({\alpha }_n^m\right)}\hfill \\ & =& {\gamma }_n^m\sum _{k=0}^{\infty }{g}_k^{\left({\alpha }_n^m\right)}+\left(-3{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2+5{\alpha }_n^m}{24}}\right)\sum _{k=0}^{\infty }{g}_{k-1}^{\left({\alpha }_n^m\right)}+\left(3{\gamma }_n^m+{\displaystyle \frac{{\alpha }_n^m-3{\left({\alpha }_n^m\right)}^3}{12}}+1\right)\sum _{k=0}^{\infty }{g}_{k-2}^{\left({\alpha }_n^m\right)}\hfill \\ & & +\left(-{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2-7{\alpha }_n^m}{24}}\right)\sum _{k=0}^{\infty }{g}_{k-3}^{\left({\alpha }_n^m\right)}\hfill \\ & =& {\gamma }_n^m\sum _{k=0}^{\infty }{g}_k^{\left({\alpha }_n^m\right)}+\left(-3{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2+5{\alpha }_n^m}{24}}\right)\sum _{k=0}^{\infty }{g}_k^{\left({\alpha }_n^m\right)}+\left(3{\gamma }_n^m+{\displaystyle \frac{{\alpha }_n^m-3{\left({\alpha }_n^m\right)}^3}{12}}+1\right)\sum _{k=0}^{\infty }{g}_k^{\left({\alpha }_n^m\right)}\hfill \\ & & +\left(-{\gamma }_n^m+{\displaystyle \frac{3{\left({\alpha }_n^m\right)}^2-7{\alpha }_n^m}{24}}\right)\sum _{k=0}^{\infty }{g}_k^{\left({\alpha }_n^m\right)}=0.\hfill \end{array}$$

□

Lemma 3.

Let ${\gamma }_i^m$ be chosen such that ${\omega }_{k,\gamma }^{\left({\alpha }_i^m\right)}$, $k=0,1,2,\dots$, satisfy

$${\omega }_{2,\gamma }^{\left({\alpha }_i^m\right)}\le 0,{\omega }_{1,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_i^m\right)}\ge 0,{\omega }_{0,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_i^m\right)}\ge 0,{\omega }_{k,\gamma }^{\left({\alpha }_i^m\right)}\ge 0,k\ge 5.$$

(23)

Then $W^m$ are diagonally dominant with negative diagonal entries.

Proof. 

Consider the ith row of $W^m$. From (21) and (23), it can be seen that all diagonal entries $w_{ii}^m\le 0$ and all off-diagonal entries $w_{ij}^m\ge 0,i\ne j$. Thus, by Lemma 2, it is easily checked that

$$\begin{array}{ccc}& & \sum _{j=1}^{N-1}{\left(W^m\right)}_{i,j}=\sum _{j=1}^{N-1}{w}_{ij}^m\hfill \\ & \le & 2{\omega }_{2,\gamma }^{\left({\alpha }_i^m\right)}+2({\omega }_{1,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_i^m\right)})+2({\omega }_{0,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_i^m\right)})+2\sum _{j=5}^{\infty }{\omega }_{j,\gamma }^{\left({\alpha }_i^m\right)}\hfill \\ & =& 2\sum _{j=0}^{\infty }{\omega }_{j,\gamma }^{\left({\alpha }_i^m\right)}=0.\hfill \end{array}$$

Equivalently,

$$\sum _{j=1,j\ne i}^{\infty }{w}_{ij}^m\le -w_{ii}^m=\left|w_{ii}^m\right|.$$

Therefore, $W^m$ is diagonally dominant with negative diagonal entries. □

Theorem 3.

Let ${\alpha }_i^m\in (1,2)$, $i=1,2,...,N-1,m=1,2,...,M$. If ${\gamma }_i^m$ satisfies

$$\begin{array}{c}\hfill {\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_i^m\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\}\end{array}$$

(24)

where

$$\left\{\begin{array}{c}{\tilde{\gamma }}_2\left({\alpha }_i^m\right)={\displaystyle \frac{3{\left({\alpha }_i^m\right)}^4+14{\left({\alpha }_i^m\right)}^3+3{\left({\alpha }_i^m\right)}^2-52{\alpha }_i^m}{8({\left({\alpha }_i^m\right)}^3+6{\left({\alpha }_i^m\right)}^2+17{\alpha }_i^m+24)}},\hfill \\ {\tilde{\gamma }}_3\left({\alpha }_i^m\right)={\displaystyle \frac{{\alpha }_i^m[3{\left({\alpha }_i^m\right)}^4+14{\left({\alpha }_i^m\right)}^3-15{\left({\alpha }_i^m\right)}^2-98{\alpha }_i^m+72]}{6({\left({\alpha }_i^m\right)}^4+6{\left({\alpha }_i^m\right)}^3+11{\left({\alpha }_i^m\right)}^2+6{\alpha }_i^m+24)}},\hfill \\ {\tilde{\gamma }}_4\left({\alpha }_i^m\right)={\displaystyle \frac{5[3{\left({\alpha }_i^m\right)}^4+14{\left({\alpha }_i^m\right)}^3-27{\left({\alpha }_i^m\right)}^2-134{\alpha }_i^m+192]}{24({\left({\alpha }_i^m\right)}^3+6{\left({\alpha }_i^m\right)}^2+11{\alpha }_i^m+6)}}.\hfill \end{array}\right.$$

(25)

Then (23) holds, i.e., ${\omega }_{2,\gamma }^{\left({\alpha }_i^m\right)}\le 0,{\omega }_{1,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_i^m\right)}\ge 0,{\omega }_{0,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_i^m\right)}\ge 0,{\omega }_{k,\gamma }^{\left({\alpha }_i^m\right)}\ge 0,k\ge 5.$

Moreover, Eq. (24) has a solution for ${\alpha }_i^m\in (1,{\alpha }^{*})$, where

$${\alpha }^{*}\approx 1.83384978$$

(26)

is the unique solution of

$$3{\mu }^8+32{\mu }^7+42{\mu }^6-556{\mu }^5-1377{\mu }^4+2204{\mu }^3+2484{\mu }^2-12048\mu +23040=0.$$

in the interval $(1,2)$.

Proof. 

For the sake of simplicity, we denote ${\alpha }_i^m$ by $\mu$. According to Lemma 2, it can be seen that

$${\omega }_{2,\gamma }^{\left(\mu \right)}={\displaystyle \frac{{\mu }^2+5\mu +6}{2}}{\gamma }_i^m-{\displaystyle \frac{\mu (3{\mu }^2+11\mu -2)}{24}}\le 0,$$

if and only if

$${\gamma }_i^m\le {\displaystyle \frac{\mu (3{\mu }^2+11\mu -2)}{12(\mu +2)(\mu +3)}}\triangleq {\tilde{\gamma }}_1\left(\mu \right).$$

Notice that

$$\begin{array}{ccc}& & {\omega }_{1,\gamma }^{\left(\mu \right)}+{\omega }_{3,\gamma }^{\left(\mu \right)}\hfill \\ & =& -(3+\mu ){\gamma }_i^m+{\displaystyle \frac{3{\mu }^2+5\mu }{24}}+{\displaystyle \frac{-{\mu }^3-6{\mu }^2-11\mu -6}{6}}{\gamma }_i^m+{\displaystyle \frac{3{\mu }^4+14{\mu }^3-3{\mu }^2-62\mu }{48}}\hfill \\ & =& -{\displaystyle \frac{{\mu }^3+6{\mu }^2+17\mu +24}{6}}{\gamma }_i^m+{\displaystyle \frac{3{\mu }^4+14{\mu }^3+3{\mu }^2-52\mu }{48}},\hfill \end{array}$$

it follows that ${\omega }_{1,\gamma }^{\left(\mu \right)}+{\omega }_{3,\gamma }^{\left(\mu \right)}\ge 0$ if and only if

$${\gamma }_i^m\le {\displaystyle \frac{3{\mu }^4+14{\mu }^3+3{\mu }^2-52\mu }{8({\mu }^3+6{\mu }^2+17\mu +24)}}\triangleq {\tilde{\gamma }}_2\left(\mu \right).$$

By

$$\begin{array}{ccc}& & {\omega }_{0,\gamma }^{\left(\mu \right)}+{\omega }_{4,\gamma }^{\left(\mu \right)}\hfill \\ & =& {\gamma }_i^m+{\displaystyle \frac{\mu ({\mu }^3+6{\mu }^2+11\mu +6)}{24}}{\gamma }_i^m-{\displaystyle \frac{\mu (3{\mu }^4+14{\mu }^3-15{\mu }^2-98\mu +72)}{144}}\hfill \\ & =& {\displaystyle \frac{{\mu }^4+6{\mu }^3+11{\mu }^2+6\mu +24}{24}}{\gamma }_i^m-{\displaystyle \frac{\mu (3{\mu }^4+14{\mu }^3-15{\mu }^2-98\mu +72)}{144}},\hfill \end{array}$$

we see that ${\omega }_{0,\gamma }^{\left(\mu \right)}+{\omega }_{4,\gamma }^{\left(\mu \right)}\ge 0$ if and only if

$${\gamma }_i^m\ge {\displaystyle \frac{\mu (3{\mu }^4+14{\mu }^3-15{\mu }^2-98\mu +72)}{6({\mu }^4+6{\mu }^3+11{\mu }^2+6\mu +24)}}\triangleq {\tilde{\gamma }}_3\left(\mu \right).$$

For $k\ge 5,$ we have

$$\begin{array}{ccc}\hfill {\omega }_{k,\gamma }^{\left(\mu \right)}& =& ({\displaystyle \frac{-{\mu }^3-6{\mu }^2-11\mu -6}{k(k-1)(k-2)}}{\gamma }_i^m\hfill \\ & & +{\displaystyle \frac{24k^2+(-12{\mu }^2-36\mu -96)k+3{\mu }^4+14{\mu }^3+33{\mu }^2+46\mu +72}{24(k-1)(k-2)}})g_{k-3}^{\left(\mu \right)}.\hfill \end{array}$$

Notice that $g_j^{\left(\mu \right)}>0$ for $j\ge 2$, we have that $w_{k,\gamma }^{\left(\mu \right)}\ge 0$ if and only if

$${\gamma }_i^m\le {\displaystyle \frac{k[24k^2+(-12{\mu }^2-36\mu -96)k+3{\mu }^4+14{\mu }^3+33{\mu }^2+46\mu +72]}{24({\mu }^3+6{\mu }^2+11\mu +6)}}.$$

Let

$$\begin{array}{ccc}\hfill g\left(x\right)& =& 24x^3+(-12{\mu }^2-36\mu -96)x^2+(3{\mu }^4+14{\mu }^3+33{\mu }^2+46\mu +72)x.\hfill \end{array}$$

It’s not difficult to verify that $g\left(x\right)$ monotonically increases. Thus $g\left(k\right)\ge g\left(5\right)$ for $k\ge 5$. In the other hand, for $k\ge 5,$ we have

$$\begin{array}{ccc}\hfill {\gamma }_i^m& \le & {\displaystyle \frac{5[600-5(12{\mu }^2+36\mu +96)+3{\mu }^4+14{\mu }^3+33{\mu }^2+46\mu +72]}{24({\mu }^3+6{\mu }^2+11\mu +6)}}\hfill \\ & =& {\displaystyle \frac{5[3{\mu }^4+14{\mu }^3-27{\mu }^2-134\mu +192]}{24({\mu }^3+6{\mu }^2+11\mu +6)}}\triangleq {\tilde{\gamma }}_4\left(\mu \right).\hfill \end{array}$$

In the following, we will show that

$${\tilde{\gamma }}_1\left(\mu \right)\ge {\tilde{\gamma }}_2\left(\mu \right)\ge {\tilde{\gamma }}_3\left(\mu \right)$$

and derive the condition under which ${\tilde{\gamma }}_3\left(\mu \right)\le {\tilde{\gamma }}_4\left(\mu \right)$ (and therefore (24) holds).

We have

$$\begin{array}{ccc}& & {\tilde{\gamma }}_1\left(\mu \right)-{\tilde{\gamma }}_2\left(\mu \right)\hfill \\ & =& {\displaystyle \frac{\mu (3{\mu }^2+11\mu -2)}{12(\mu +2)(\mu +3)}}-{\displaystyle \frac{3{\mu }^4+14{\mu }^3+3{\mu }^2-52\mu }{8({\mu }^3+6{\mu }^2+17\mu +24)}}\hfill \\ & =& {\displaystyle \frac{2\mu (3{\mu }^2+11\mu -2)({\mu }^3+6{\mu }^2+17\mu +24)}{24(\mu +2)(\mu +3)({\mu }^3+6{\mu }^2+17\mu +24)}}-{\displaystyle \frac{3(3{\mu }^4+14{\mu }^3+3{\mu }^2-52\mu )({\mu }^2+5\mu +6)}{24(\mu +2)(\mu +3)({\mu }^3+6{\mu }^2+17\mu +24)}}\hfill \\ & =& {\displaystyle \frac{-3{\mu }^6-29{\mu }^5-43{\mu }^4+353{\mu }^3+1186{\mu }^2+840\mu }{24({\mu }^4+8{\mu }^3+29{\mu }^2+58\mu +48)(\mu +3)}}\hfill \\ & =& {\displaystyle \frac{(-3{\mu }^5-20{\mu }^4+17{\mu }^3+302{\mu }^2+280\mu )(\mu +3)}{24({\mu }^4+8{\mu }^3+29{\mu }^2+58\mu +48)(\mu +3)}}\hfill \\ & =& {\displaystyle \frac{-3{\mu }^5-20{\mu }^4+17{\mu }^3+302{\mu }^2+280\mu }{24({\mu }^4+8{\mu }^3+29{\mu }^2+58\mu +48)}}.\hfill \end{array}$$

Let

$${\tilde{\gamma }}_{12}\left(\mu \right)=-3{\mu }^5-20{\mu }^4+17{\mu }^3+302{\mu }^2+280\mu ,\mu \in (1,2).$$

Then

$${\tilde{\gamma }}_{12}^{\prime }\left(\mu \right)=-15{\mu }^4-80{\mu }^3+51{\mu }^2+604\mu +280.$$

It’s easy to check that ${\tilde{\gamma }}_{12}^{\prime }\left(\mu \right)\ge 0$ for $\mu \in (1,2).$ Thus ${\tilde{\gamma }}_{12}\left(\mu \right)$ monotonically increases. It follows that ${\tilde{\gamma }}_{12}\left(\mu \right)\ge {\tilde{\gamma }}_{12}\left(1\right)=576>0$ for $\mu \in (1,2),$ and therefore ${\tilde{\gamma }}_1\left(\mu \right)\ge {\tilde{\gamma }}_2\left(\mu \right).$

Next, we prove that ${\tilde{\gamma }}_2\left(\mu \right)\ge {\tilde{\gamma }}_3\left(\mu \right)$. We have

$$\begin{array}{ccc}& & {\tilde{\gamma }}_2\left(\mu \right)-{\tilde{\gamma }}_3\left(\mu \right)\hfill \\ & =& {\displaystyle \frac{3{\mu }^4+14{\mu }^3+3{\mu }^2-52\mu }{8({\mu }^3+6{\mu }^2+17\mu +24)}}-{\displaystyle \frac{\mu (3{\mu }^4+14{\mu }^3-15{\mu }^2-98\mu +72)}{6({\mu }^4+6{\mu }^3+11{\mu }^2+6\mu +24)}}\hfill \\ & =& {\displaystyle \frac{{\tilde{\gamma }}_{23}\left(\mu \right)}{24({\mu }^7+12{\mu }^6+64{\mu }^5+198{\mu }^4+391{\mu }^3+510{\mu }^2+552\mu +576)}},\hfill \end{array}$$

where

$$\begin{array}{c}\hfill {\tilde{\gamma }}_{23}\left(\mu \right)=-3{\mu }^8-32{\mu }^7-120{\mu }^6-74{\mu }^5+1371{\mu }^4+5722{\mu }^3+3792{\mu }^2-10656\mu .\end{array}$$

The derivative of ${\tilde{\gamma }}_{23}\left(\mu \right)$ is

$$\begin{array}{c}\hfill {\tilde{\gamma }}_{23}^{\prime }\left(\mu \right)=-24{\mu }^7-224{\mu }^6-720{\mu }^5-370{\mu }^4+5484{\mu }^3+17166{\mu }^2+7584\mu -10656.\end{array}$$

It’s easy to check that ${\tilde{\gamma }}_{23}^{\prime }\left(\mu \right)\ge 0,$ $\mu \in (1,2).$ Therefore,

$${\tilde{\gamma }}_{23}\left(\mu \right)\ge {\tilde{\gamma }}_{23}\left(1\right)=18240>0,\mu \in (1,2).$$

Hence, ${\tilde{\gamma }}_2\left(\mu \right)\ge {\tilde{\gamma }}_3\left(\mu \right).$

In order that (24) has a solution, it is required that ${\tilde{\gamma }}_4\left(\mu \right)\ge {\tilde{\gamma }}_3\left(\mu \right)$. By some computation, we get

$$\begin{array}{ccc}& & {\tilde{\gamma }}_4\left(\mu \right)-{\tilde{\gamma }}_3\left(\mu \right)\hfill \\ & =& {\displaystyle \frac{5(3{\mu }^4+14{\mu }^3-27{\mu }^2-134\mu +192)}{24({\mu }^3+6{\mu }^2+11\mu +6)}}-{\displaystyle \frac{\mu (3{\mu }^4+14{\mu }^3-15{\mu }^2-98\mu +72)}{6({\mu }^4+6{\mu }^3+11{\mu }^2+6\mu +24)}}\hfill \\ & =& {\displaystyle \frac{{\tilde{\gamma }}_{34}\left(\mu \right)}{24({\mu }^7+12{\mu }^6+58{\mu }^5+144{\mu }^4+217{\mu }^3+276{\mu }^2+300\mu +144)}},\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill {\tilde{\gamma }}_{34}\left(\mu \right)& =& 3{\mu }^8+32{\mu }^7+42{\mu }^6-556{\mu }^5-1377{\mu }^4+2204{\mu }^3+2484{\mu }^2-12048\mu +23040.\hfill \end{array}$$

It is easy to get

$$\begin{array}{c}\hfill {\tilde{\gamma }}_{34}^{\prime }\left(\mu \right)=24{\mu }^7+224{\mu }^6+252{\mu }^5-2780{\mu }^4-5508{\mu }^3+6612{\mu }^2+4968\mu -12048.\end{array}$$

We can verify that ${\tilde{\gamma }}_{34}^{\prime }\left(\mu \right)$ decreases for $\mu \in (1,2)$ and it follows that ${\tilde{\gamma }}_{34}^{\prime }\left(\mu \right)\le {\tilde{\gamma }}_{34}^{\prime }\left(1\right)=-8256<0.$ Therefore, ${\tilde{\gamma }}_{34}\left(\mu \right)$ also decreases. Notice that

$${\tilde{\gamma }}_{34}\left(1\right)=13824>0,{\tilde{\gamma }}_{34}\left(2\right)=-5760<0,$$

we have that there exists a unique ${\alpha }^{*}\in (1,2)$ such that ${\tilde{\gamma }}_{34}\left({\alpha }^{*}\right)=0.$ By using Matlab, we get ${\alpha }^{*}\approx 1.83384978$. Thus, ${\tilde{\gamma }}_4\left(\mu \right)\ge {\tilde{\gamma }}_3\left(\mu \right)$ for $\mu \in (1,{\alpha }^{*}),$ and (24) has solution(s). □

Lemma 4.

[1] (The barrier lemma) Assume $M_A$ is a a monotone matrix with $M_A\in {\mathbb{R}}^{n\times n}.$ If there is a unit vector $\overrightarrow{\varrho }$, i.e., $\parallel \overrightarrow{\varrho }{\parallel }_{\infty }=1$, and $\underset{1\le i\le n}{min}{\left(M_A\overrightarrow{\varrho }\right)}_i\ge a_M$ with $a_M>0$ a positive constant. Then one gains $\parallel M_A^{-1}{\parallel }_{\infty }\le 1/a_M.$

Define $\mathbf{e}={(1,1,\dots ,1)}^T$. For any vectors $\mathbf{x},\mathbf{y}\in {\mathbb{R}}^n$, denote $\mathbf{x}\ge (\mathrm{or}>)\mathbf{y}$ if ${[\mathbf{x}-\mathbf{y}]}_i\ge 0$ (or ${[\mathbf{x}-\mathbf{y}]}_i>0$) $i=1,2,\dots ,n$.

Lemma 5.

Suppose $1<{\alpha }_i^m\le {\alpha }^{*}$ with ${\alpha }^{*}$ being given by (26) and ${\gamma }_i^m$ is chosen such that ${\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_i^m\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\}$, $i=1,2,\dots ,N-1,$ where ${\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_3\left({\alpha }_i^m\right)$ and ${\tilde{\gamma }}_4\left({\alpha }_i^m\right)$ are given by (25). Then $\parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }\le 1$.

Proof. 

Notice that $A^m=I-D^m{W}^m$ with $W^m$ given by (21), we are aware of that $A^m$ is an L-matrix. By Lemmas 2 and 3, we have

$$\begin{array}{ccc}& & \sum _{j=1}^{N-1}{\left[A^m\right]}_{ij}=1-d_i^m\sum _{j=1}^{N-1}{w}_{ij}^{\left({\alpha }_i^m\right)}\hfill \\ & \ge & 1-\left[2{\omega }_{2,\gamma }^{\left({\alpha }_i^m\right)}+2({\omega }_{1,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_i^m\right)})+2({\omega }_{0,\gamma }^{\left({\alpha }_i^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_i^m\right)})+2\sum _{j=5}^{\infty }{\omega }_{j,\gamma }^{\left({\alpha }_i^m\right)}\right]\hfill \\ & \ge & 1-2\sum _{j=0}^{\infty }{\omega }_{j,\gamma }^{\left({\alpha }_i^m\right)}\ge 1,i=1,2,\dots ,N-1.\hfill \end{array}$$

Hence, $A^m\mathbf{e}\ge \mathbf{e}>\mathbf{0}.$ It follows that $A^m$ is an M-matrix [12]. Furthermore, $\parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }\le \underset{1\le i\le N}{min}{\left[\mathbf{e}\right]}_i=1$ from Lemma 4. □

Lemma 6.

[45] (Discrete Gronwall Inequality) Let $\left\{{\zeta }^n\right|n\ge 0\}$ be a sequence with ${\zeta }^n>0$, and ${\varrho }^n\ge 0.$ If ${\zeta }^{n+1}\le (1+c\tau ){\zeta }^n+\tau {\varrho }^n,n=0,1,2,\dots ,$ then ${\zeta }^n$ satisfies

$${\zeta }^{n+1}\le e^{cn\tau }\sum _{k=0}^n\tau {\varrho }^k,$$

where ${\zeta }^0=0,$ c is a nonnegative constant.

Theorem 4.

If $1<{\alpha }_i^m\le {\alpha }^{*}$ with ${\alpha }^{*}$ being given by (26) and ${\gamma }_i^m$ is chosen such that ${\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_i^m\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\},$ $i=1,2,\dots ,N-1,m=1,2\dots ,M,$ where ${\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_3\left({\alpha }_i^m\right)$ and ${\tilde{\gamma }}_4\left({\alpha }_i^m\right)$ are given by (25). Then the IE-4WSGD scheme (19) is stable.

Proof. 

Suppose ${\tilde{U}}^m$ is the approximation solution of $U^m$ and ${\tilde{F}}^m={({\tilde{f}}_1^m,{\tilde{f}}_2^m,\dots ,{\tilde{f}}_{N-1}^{m-1})}^T$. Let $E^m=U^m-{\tilde{U}}^m.$ Then we have

$$A^m{E}^m=E^{m-1}+\tau (F^{m-1}-{\tilde{F}}^{m-1}),m=1,2,\dots ,M,$$

that is,

$$E^m={\left(A^m\right)}^{-1}{E}^{m-1}+\tau {\left(A^m\right)}^{-1}(F^{m-1}-{\tilde{F}}^{m-1}).$$

From Lemma 5, we get

$$\begin{array}{ccc}\hfill \parallel E^m{\parallel }_{\infty }& \le & \parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }\parallel E^{m-1}{\parallel }_{\infty }+\tau L\parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }{\parallel E^{m-1}\parallel }_{\infty }\hfill \\ & \le & (1+\tau L)\parallel E^{m-1}{\parallel }_{\infty }\le {(1+\tau L)}^m\parallel E^0{\parallel }_{\infty }\le e^{m\tau L}\parallel E^0{\parallel }_{\infty }\le e^{TL}{\parallel E^0\parallel }_{\infty }.\hfill \end{array}$$

This can lead to the conclusion. □

Theorem 5.

Denote $u_n^m$ the exact solution of IBV problem (1) and $U_n^m$ the solution of (19). If $1<{\alpha }_i^m\le {\alpha }^{*}$ and ${\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_i^m\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\},$ $i=1,2,\dots ,N-1,m=1,2,\dots ,M.$ Then for $1\le m\le M$, it holds

$$\parallel U^m-u^m{\parallel }_{\infty }\le C{e}^{TL}T(\tau +h^l),$$

where $C>0$ and $l\in \{3,4\}$.

Proof. 

Denote $e^m=U^m-u^m$. From Eqs. (1) and (19), we get

$$A^m{e}^m=e^{m-1}+\tau (F_e^{m-1}-F^{m-1})+\tau R^m,m=1,2,\dots ,M,$$

where $\parallel R^m{\parallel }_{\infty }\le C(\tau +h^l),l\in \{3,4\}.$ That is,

$$e^m={\left(A^m\right)}^{-1}{e}^{m-1}+\tau {\left(A^m\right)}^{-1}(F_e^{m-1}-F^{m-1})+\tau {\left(A^m\right)}^{-1}{R}^m.$$

Notice that $e^0=0,$ from Lemmas 5 and 6, we get

$$\begin{array}{ccc}\hfill \parallel e^m{\parallel }_{\infty }& \le & \parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }\parallel e^{m-1}{\parallel }_{\infty }+\tau L\parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }\parallel e^{m-1}{\parallel }_{\infty }+\tau \parallel {\left(A^m\right)}^{-1}{\parallel }_{\infty }{\parallel R^m\parallel }_{\infty }\hfill \\ & \le & (1+\tau L)\parallel e^{m-1}{\parallel }_{\infty }+\tau {\parallel R^m\parallel }_{\infty }\le e^{m\tau L}Cm\tau (\tau +h^l)\le C{e}^{TL}T(\tau +h^l).\hfill \end{array}$$

This can lead to the conclusion. □

In the following, we discuss how to choose the parameters ${\gamma }_i^m$ such that the IE-4WSGD scheme is stable and has the optimal error bound.

Theorem 6.

Let ${\gamma }_{obj}\left({\alpha }_i^m\right)={\displaystyle \frac{{\left({\alpha }_i^m\right)}^3-{\left({\alpha }_i^m\right)}^2-4{\alpha }_i^m}{48}},i=1,2,\dots ,N-1,m=1,2,\dots ,M$ and let

$${\alpha }^{**}\approx 1.79968326$$

(27)

be the unique solution of $-{\mu }^6-5{\mu }^5+29{\mu }^4+169{\mu }^3-220{\mu }^2-1316\mu +1920=0$ in the interval $(1,2)$. Then for ${\alpha }_i^m\in \left(1,{\alpha }^{**}\right)$, it holds

$${\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_{obj}\left({\alpha }_i^m\right)\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\},$$

(28)

where ${\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_3\left({\alpha }_i^m\right)$ and ${\tilde{\gamma }}_4\left({\alpha }_i^m\right)$ are given by (25). That is, if $\alpha (x,t)\in (0,{\alpha }^{**})$, then we can select the free parameter of the 4WSGD scheme so that the corresponding scheme has fourth order of accuracy.

Proof. 

We first compare ${\gamma }_{obj}\left(\mu \right)$ and ${\tilde{\gamma }}_3\left(\mu \right)$. We have

$$\begin{array}{ccc}& & {\gamma }_{obj}\left(\mu \right)-{\tilde{\gamma }}_3\left(\mu \right)\hfill \\ & =& {\displaystyle \frac{{\mu }^3-{\mu }^2-4\mu }{48}}-{\displaystyle \frac{\mu (3{\mu }^4+14{\mu }^3-15{\mu }^2-98\mu +72)}{6({\mu }^4+6{\mu }^3+11{\mu }^2+6\mu +24)}}\hfill \\ & =& {\displaystyle \frac{{\gamma }_{obj3}\left(\mu \right)}{48({\mu }^4+6{\mu }^3+11{\mu }^2+6\mu +24)}},\hfill \end{array}$$

where

$${\gamma }_{obj3}\left(\mu \right)={\mu }^7+5{\mu }^6-23{\mu }^5-141{\mu }^4+94{\mu }^3+736{\mu }^2-672\mu .$$

It follows that

$${\gamma }_{obj3}^{\prime }\left(\mu \right)=7{\mu }^6+30{\mu }^5-115{\mu }^4-564{\mu }^3+282{\mu }^2+1472\mu -672.$$

It is easily checked that ${\gamma }_{obj3}^{\prime }\left(\mu \right)$ increases first and then decreases with ${\gamma }_{obj3}^{\prime }\left(1\right)>0$ and ${\gamma }_{obj3}^{\prime }\left({\alpha }^{*}\right)<0$ for $(1,{\alpha }^{*})$. Thus ${\gamma }_{obj3}\left(\mu \right)$ also increases first and then decreases and there exists a minimum at $\mu =1$ or $\mu ={\alpha }^{*}.$ Since ${\gamma }_{obj3}\left(1\right)=0,$ ${\gamma }_{obj3}\left({\alpha }^{*}\right)>0,$ therefore ${\gamma }_{obj3}\left(\mu \right)\ge 0,$ that is ${\gamma }_{obj}\left(\mu \right)\ge {\tilde{\gamma }}_3\left(\mu \right)$ for $(1,{\alpha }^{*})$.

Next we compare ${\tilde{\gamma }}_2\left(\mu \right)$ and ${\gamma }_{obj}\left(\mu \right)$. We have

$$\begin{array}{ccc}& & {\tilde{\gamma }}_2\left(\mu \right)-{\gamma }_{obj}\left(\mu \right)\hfill \\ & =& {\displaystyle \frac{3{\mu }^4+14{\mu }^3+3{\mu }^2-52\mu }{8({\mu }^3+6{\mu }^2+17\mu +24)}}-{\displaystyle \frac{{\mu }^3-{\mu }^2-4\mu }{48}}\hfill \\ & =& {\displaystyle \frac{{\gamma }_{obj2}\left(\mu \right)}{48({\mu }^3+6{\mu }^2+17\mu +24)}},\hfill \end{array}$$

where

$${\gamma }_{obj2}\left(\mu \right)=-{\mu }^6-5{\mu }^5+11{\mu }^4+101{\mu }^3+110{\mu }^2-216\mu .$$

Then it is easy to get

$${\gamma }_{obj2}^{\prime }\left(\mu \right)=-6{\mu }^5-25{\mu }^4+44{\mu }^3+303{\mu }^2+220\mu -216,$$

and one can verify that ${\gamma }_{obj2}^{\prime }\left(\mu \right)$ increases monotonically for $\mu \in (1,{\alpha }^{*})$. Notice that ${\gamma }_{obj2}^{\prime }\left(1\right)=320>0$, we have that ${\gamma }_{obj2}\left(\mu \right)\ge {\gamma }_{obj2}\left(1\right)\ge 0.$ Equivalently, ${\tilde{\gamma }}_2\left(\mu \right)\ge {\gamma }_{obj}\left(\mu \right)$ for $(1,{\alpha }^{*})$.

Finally, we compare ${\tilde{\gamma }}_4\left(\mu \right)$ and ${\gamma }_{obj}\left(\mu \right)$. We have

$$\begin{array}{ccc}& & {\tilde{\gamma }}_4\left(\mu \right)-{\gamma }_{obj}\left(\mu \right)\hfill \\ & =& {\displaystyle \frac{5(3{\mu }^4+14{\mu }^3-27{\mu }^2-134\mu +192)}{24({\mu }^3+6{\mu }^2+11\mu +6)}}-{\displaystyle \frac{{\mu }^3-{\mu }^2-4\mu }{48}}\hfill \\ & =& {\displaystyle \frac{{\gamma }_{obj4}\left(\mu \right)}{48({\mu }^3+6{\mu }^2+11\mu +6)}},\hfill \end{array}$$

where

$${\gamma }_{obj4}\left(\mu \right)=-{\mu }^6-5{\mu }^5+29{\mu }^4+169{\mu }^3-220{\mu }^2-1316\mu +1920.$$

It can be checked that

$${\gamma }_{obj4}^{\prime }\left(\mu \right)=-6{\mu }^5-25{\mu }^4+116{\mu }^3+507{\mu }^2-440\mu -1316$$

increases monotonically and ${\gamma }_{obj4}^{\prime }\left({\alpha }^{*}\right)<0.$ Therefore, ${\gamma }_{obj4}\left(\mu \right)$ decreases monotonically. Notice that

$${\gamma }_{obj4}\left(1\right)>0\mathrm{and}{\gamma }_{obj4}\left({\alpha }^{*}\right)\le 0,$$

we see that there exists a unique ${\alpha }^{**}\in (1,{\alpha }^{*})$ such that ${\gamma }_{obj4}\left({\alpha }^{**}\right)=0$. By using Matlab, we get an approximate value of ${\alpha }^{**}$, i.e., ${\alpha }^{**}\approx 1.79968326$. That is ${\gamma }_{obj4}\left(\mu \right)>0$ for ${\alpha }_i^m\in (1,{\alpha }^{**})$, and ${\gamma }_{obj4}\left({\alpha }^{*}\right)<0$ for ${\alpha }_i^m\in ({\alpha }^{**},{\alpha }^{*})$.

In summary,

$${\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_{obj}\left({\alpha }_i^m\right)\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\}{\alpha }_i^m\in \left(1,{\alpha }^{**}\right).$$

□

To end this section, we discuss how to choose the parameter ${\gamma }_i^m$ for ${\alpha }^{**}<{\alpha }_i^m\le {\alpha }^{*}$, where ${\alpha }^{*}$ and ${\alpha }^{**}$ are given by (26) and (27) respectively. The minimization problem will take into account

$$\begin{array}{ccc}& & min|{\gamma }_i^m-{\displaystyle \frac{{\left({\alpha }_i^m\right)}^3-{\left({\alpha }_i^m\right)}^2-4{\alpha }_i^m}{48}}|\hfill \\ \hfill s.t.& & {\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_i^m\le min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4\left({\alpha }_i^m\right)\},{\alpha }^{**}\le {\alpha }_i^m\le {\alpha }^{*}.\hfill \end{array}$$

To seek above-mentioned problem, we draw the curves below (see Figure 1):

• Objective: ${\gamma }_{obj}\left({\alpha }_i^m\right),$
• Lower bound: ${\tilde{\gamma }}_3\left({\alpha }_i^m\right),$
• Upper bound: $min\{{\tilde{\gamma }}_2\left({\alpha }_i^m\right),{\tilde{\gamma }}_4{\left(\alpha \right)}_i^m\}.$

Combining Theorem 6 and the curves in Figure 1, it can be easily seen that ${\tilde{\gamma }}_4\left({\alpha }_i^m\right)\le {\gamma }_{obj}\left({\alpha }_i^m\right),{\tilde{\gamma }}_3\left({\alpha }_i^m\right)\le {\gamma }_{obj}\left({\alpha }_i^m\right)$ for $[{\alpha }^{**},{\alpha }^{*})$, and ${\tilde{\gamma }}_4\left({\alpha }_i^m\right)$ is closer to the objective ${\tilde{\gamma }}_{obj}\left({\alpha }_i^m\right)$. Therefore, we can choose ${\gamma }_i^m={\tilde{\gamma }}_4\left({\alpha }_i^m\right)$ for $\left[{\alpha }^{**},{\alpha }^{*}\right).$ In summary, select ${\gamma }_i^m$ as follows:

$${\gamma }_i^m=\left\{\begin{array}{c}{\gamma }_{obj}\left({\alpha }_i^m\right),{\alpha }_i^m\in \left(1,{\alpha }^{**}\right)\hfill \\ {\tilde{\gamma }}_4\left({\alpha }_i^m\right),{\alpha }_i^m\in \left[{\alpha }^{**},{\alpha }^{*}\right).\hfill \end{array}\right.$$

(29)

4. Preconditioning Techniques

In this section, based on Section 3, we consider applying preconditioned GMRES methods to solve the linear systems in Eq. (19). Although the coefficient matrices are dense and do not have a Toeplitz-like structure, each $W^m$ has off-diagonal decay property, hence we employ banded matrices as preconditioners for the linear systems to speed up the relevant GMRES methods. Similar discussion can be found in the literature; see for instance [18,42].

Firstly, a banded preconditioner with diagonal compensation will be take into account. Letting $Q={\left(q_{ij}\right)}_{i,j=1}^N$, we denote $\mathrm{diag}(\sum _{j=1}^N{q}_{1j},\sum _{j=1}^N{q}_{2j},\cdots ,\sum _{j=1}^N{q}_{Nj})$ by $\mathrm{diag}\left(\mathrm{sum}\right(Q,2\left)\right)$. Let $k_b$ be an integer that is much less than N. Split $W^m$ as

$$W^m=W_{k_b}^m+(W^m-W_{k_b}^m),$$

where

$$W_{k_b}^m={\tilde{W}}_{k_b}^m+\mathrm{diag}\left(\mathrm{sum}(W^m-{\tilde{W}}_k^m,2)\right)$$

with

$${\tilde{W}}_{k_b}^m=\left(\begin{array}{cccccc}{w}_{1,1}^{\left(m\right)}& \cdots & w_{1,k_b}^{\left(m\right)}& 0& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & \ddots & ⋮\\ w_{k_b,1}^{\left(m\right)}& \cdots & w_{k_b,k_b}^{\left(m\right)}& \ddots & \ddots & 0\\ 0& \ddots & \ddots & \ddots & \ddots & w_{N-k_b,N-1}^{\left(m\right)}\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& w_{N-1,N-k_b}^{\left(m\right)}& \cdots & w_{N-1,N-1}^{\left(m\right)}\end{array}\right).$$

(30)

Then we use

$$P_{k_b}^m=I-D^m{W}_{k_b}^m$$

as preconditioner.

Lemma 7.

For $k_b\ge 3,$ we have

$$\begin{array}{c}\hfill \parallel W^m-{\tilde{W}}_{k_b}^m{\parallel }_{\infty }\le r_{\alpha ,k_b}^m,{\parallel W^m-W_{k_b}^m\parallel }_{\infty }\le 2r_{\alpha ,k_b}^m,\end{array}$$

where $r_{\alpha ,k_b}^m=4\underset{1\le n\le N-1}{max}\left(\right|\gamma |,|{\gamma }_1|,|{\gamma }_0|,|{\gamma }_{-1}\left|\right)\underset{0\le n\le N}{max}\left({\displaystyle \frac{({\alpha }_n^m-1)3^{{\alpha }_n^m+1}}{{(k_b-1)}^{{\alpha }_n^m}}}\right).$

Proof. 

From (21) and Lemma 1, we get

$$\begin{array}{ccc}& & \parallel W^m-{\tilde{W}}_{k_b}^m{\parallel }_{\infty }\hfill \\ & \le & 2\underset{1\le n\le N-1}{max}\left(\sum _{q=k_b+2}^{\infty }{\omega }_{q,\gamma }^{\left({\alpha }_n^m\right)}\right)\hfill \\ & \le & 8\underset{1\le n\le N-1}{max}\left(\right|\gamma |,|{\gamma }_1|,|{\gamma }_0|,|{\gamma }_{-1}\left|\right)\underset{1\le i\le N-1}{max}\left(\sum _{q=k_b-1}^{\infty }{g}_q^{\left({\alpha }_n^m\right)}\right)\hfill \\ & \le & 8\underset{1\le n\le N-1}{max}\left(\right|\gamma |,|{\gamma }_1|,|{\gamma }_0|,|{\gamma }_{-1}\left|\right)\underset{1\le n\le N-1}{max}\left({\displaystyle \frac{({\alpha }_n^m-1)3^{{\alpha }_n^m+1}}{2{(k_b-1)}^{{\alpha }_n^m}}}\right)\hfill \\ & \triangleq & r_{\alpha ,k_b}^m.\hfill \end{array}$$

Similarly, we can get

$$\begin{array}{c}\hfill \parallel W^m-W_{k_b}^m{\parallel }_{\infty }\le 2r_{\alpha ,k_b}^m.\end{array}$$

□

Theorem 7.

Let $d_{max}^m=\underset{1\le n\le N-1}{max}\left(d_n^m\right)$ and $r_{\alpha ,k_b}^m$ be the same as in Lemma 7. For $k_b\ge 3,$ we have following conclusions:

(i)

The spectrum radius of ${\left(P_{k_b}^m\right)}^{-1}{A}^m-I$ is bounded as follows:

$$\rho \left({\left(P_{k_b}^m\right)}^{-1}{A}^m-I\right)\le 2d_{max}^m{r}_{\alpha ,k_b}^m;$$

(ii)

The condition number of the preconditioned matrix is bounded as follows:

$${\mathrm{cond}}_{\infty }\left({\left(P_{k_b}^m\right)}^{-1}\left(A^m\right)\right)\le {(1+2d_{max}^m{r}_{\alpha ,k_b}^m)}^2,m=1,2,\dots ,M.$$

Proof. 

Similar to Lemma 5, we have ${∥{\left(I-D^m{W}_{k_b}^m\right)}^{-1}∥}_{\infty }\le 1$ by using Lemma 4.

(i) Utilizing Lemma 7, we get

$$\begin{array}{ccc}& & {∥{\left(P_{k_b}^m\right)}^{-1}{A}^m-I∥}_{\infty }\hfill \\ & =& {∥{(I-D^m{W}_{k_b}^m)}^{-1}(I-D^m{W}^m)-I∥}_{\infty }\hfill \\ & =& {∥{(I-D^m{W}_{k_b}^m)}^{-1}\left[I-D^m\left(W_{k_b}^m+(W^m-W_{k_b}^m)\right)\right]-I∥}_{\infty }\hfill \\ & =& {∥{(I-D^m{W}_{k_b}^m)}^{-1}{D}^m(W^m-W_{k_b}^m)∥}_{\infty }\hfill \\ & \le & {∥{(I-D^m{W}_{k_b}^m)}^{-1}∥}_{\infty }{∥D^m(W^m-W_{k_b}^m)∥}_{\infty }\hfill \\ & \le & \underset{1\le n\le N-1}{max}\left(d_n^m\right){∥W^m-W_{k_b}^m∥}_{\infty }\le 2d_{max}^m{r}_{\alpha ,k_b}^m.\hfill \end{array}$$

Therefore, we get $\rho \left({\left(P_{k_b}^m\right)}^{-1}\left(A^m\right)-I\right)\le {∥{\left(P_{k_b}^m\right)}^{-1}\left(A^m\right)-I∥}_{\infty }\le 2d_{max}^m{r}_{\alpha ,k_b}^m.$

(ii) We have

$$\begin{array}{ccc}\hfill {∥{\left(P_{k_b}^m\right)}^{-1}{A}^m∥}_{\infty }& =& {∥{\left(P_{k_b}^m\right)}^{-1}{A}^m-I+I∥}_{\infty }\hfill \\ & \le & 1+{∥{\left(P_{k_b}^m\right)}^{-1}{A}^m-I∥}_{\infty }\hfill \\ & \le & 1+2d_{max}^m{r}_{\alpha ,k_b}^m.\hfill \end{array}$$

Similarly, we have

$$\begin{array}{ccc}& & {∥{\left[{\left(P_{k_b}^m\right)}^{-1}{A}^m\right]}^{-1}∥}_{\infty }\hfill \\ & =& {∥{(I-D^m{W}^m)}^{-1}(I-D^m{W}_{k_b}^m)∥}_{\infty }\hfill \\ & =& {∥I+{(I-D^m{W}^m)}^{-1}{D}^m(W^m-W_{k_b}^m)∥}_{\infty }\hfill \\ & \le & 1+{∥{(I-D^m{W}^m)}^{-1}∥}_{\infty }{∥D^m(W^m-W_{k_b}^m)∥}_{\infty }\hfill \\ & \le & 1+2d_{max}^m{r}_{\alpha ,k_b}^m.\hfill \end{array}$$

The conclusion of the theorem follows. □

Now, a banded preconditioner without diagonal compensation will be take into account.

$${\tilde{P}}_{k_b}^m=I-D^m{\tilde{W}}_{\alpha ,k_b}^m.$$

Theorem 8.

Let $d_{max}^m=\underset{1\le n\le N-1}{max}\left(d_n^m\right)$ and $r_{\alpha ,k_b}^m$ be the same as in Lemma 7. For $k_b\ge 3,$ we have following conclusions:

(i)

The spectrum radius of ${\left({\tilde{P}}_{k_b}^m\right)}^{-1}{A}^m-I$ is bounded as follows:

$$\rho \left({\left({\tilde{P}}_{k_b}^m\right)}^{-1}{A}^m-I\right)\le d_{max}^m{r}_{\alpha ,k_b}^m.$$

(ii)

The condition number of the preconditioned matrix is bounded as follows:

$${\mathrm{cond}}_{\infty }\left({\left({\tilde{P}}_{k_b}^m\right)}^{-1}\left(A^m\right)\right)\le {(1+d_{max}^m{r}_{\alpha ,k_b}^m)}^2,m=1,2,\dots ,M.$$

Proof. 

By using Lemma 7, the theorem can be proved similar to the proof of Theorem 7. □

5. IE-4WSGD Scheme for 2D SRVONFDE

Nowadays, we take into account the IBV problem for a 2D SRVONFDE:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{\partial u(x,y,t)}{\partial t}=K_{\alpha }(x,y,t)\frac{{\partial }^{\alpha (x,y,t)}u(x,y,t)}{\partial {\left|x\right|}^{\alpha (x,y,t)}}+K_{\beta }(x,y,t)\frac{{\partial }^{\beta (x,y,t)}u(x,y,t)}{\partial {\left|y\right|}^{\beta (x,y,t)}}+f(u,x,y,t),}\hfill \\ (x,y)\in \Omega =(x_L,x_R)\times (y_L,y_R),t\in (0,T],\hfill & \\ u(x,y,0)=u_0(x,y),(x,y)\in \overline{\Omega },\hfill \\ u(x,y,t)=0,(x,y)\in \partial \Omega ,t\in [0,T],\hfill \end{array}\right.\end{array}$$

(31)

where ${\displaystyle \frac{{\partial }^{\alpha (x,y,t)}u(x,y,t)}{\partial {\left|x\right|}^{\alpha (x,y,t)}}}$ and ${\displaystyle \frac{{\partial }^{\beta (x,y,t)}u(x,y,t)}{\partial {\left|y\right|}^{\beta (x,y,t)}}}$ are Riesz VO fractional derivatives, $0<\underline{K_{C\alpha }}\le K_{C\alpha }(x,y,t)=K_{\alpha }(x,y,t)C_{\alpha }(x,y,t)\le \overline{K_{\alpha }}$ for certain positive constants $\underline{K_{C\alpha }}$ and $\overline{K_{C\alpha }}$, $0<\underline{K_{C\beta }}\le K_{C\beta }(x,y,t)=K_{\beta }(x,y,t)C_{\beta }(x,y,t)\le \overline{K_{C\beta }}$ for certain positive constants $\underline{K_{C\beta }}$ and $\overline{K_{C\beta }}$, $f(u,x,y,t)$ is the source term that satisfies the Lipschitz condition.

We discuss discretization of the IBV problem firstly. $N_x,N_y,M\in {\mathbb{N}}^{*}$ Let $h_x={\displaystyle \frac{x_R-x_L}{N_x}}$, $h_y={\displaystyle \frac{y_R-y_L}{N_y}}$ and $\tau ={\displaystyle \frac{T}{M}}$. Let $x_i=x_L+i{h}_x,0\le i\le N_x$, $y_j=y_L+j{h}_y,0\le j\le N_y,$ and $t_m=m\tau ,m=0,1,2,\cdots ,M.$ $g(x,y,t)$ is a given function, let $g_{ij}^m$ or $g_{i,j}^m$ represent $g(x_i,y_j,t_m)$ . The IE-4WSGD scheme presented in Section 3 can also be applied to 2D RVONFDE, i.e., Eq. (31). The IE-4WSGD scheme for Eq. (31) is

$$\begin{array}{ccc}\hfill u_{ij}^m& =& u_{ij}^{m-1}+{\displaystyle \frac{\tau K_{C\alpha ,i,j}^m}{h_x^{{\alpha }_{i,j}^{m-1}}}}\left[\sum _{k=0}^{i+2}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{u}_{i-k+2,j}^m+\sum _{k=0}^{N_x-i+3}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{u}_{i+k-2,j}^m\right]\hfill \\ & & +{\displaystyle \frac{\tau K_{C\beta ,i,j}^m}{h_y^{{\beta }_{i,j}^{m-1}}}}\left[\sum _{k=0}^{j+2}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{u}_{i,j-k+2}^m+\sum _{k=0}^{N_y-j+3}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{u}_{i,j+k-2}^m\right]+\tau f_{e,ij}^{m-1}+\tau R_{ij}^m,\hfill \end{array}$$

(32)

where $f_{e,ij}^m=f(u_{ij}^m,x_i,y_j,t_m)$ and $|R_{ij}^m|\le C_{xy}(\tau +h_x^l+h_y^l),l\in \{3,4\}.$ Let $U_{ij}^m$ be an approximation of $u_{ij}^m$, we have the following systems of linear equations:

$$\begin{array}{ccc}\hfill U_{ij}^m& -& {\displaystyle \frac{\tau K_{C\alpha ,i,j}^m}{h_x^{{\alpha }_{i,j}^{m-1}}}}\left[\sum _{k=0}^{i+2}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{U}_{i-k+2,j}^m+\sum _{k=0}^{N_x-i+3}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{U}_{i+k-2,j}^m\right]\hfill \\ & -& {\displaystyle \frac{\tau K_{C\beta ,i,j}^m}{h_y^{{\beta }_{i,j}^{m-1}}}}\left[\sum _{k=0}^{j+2}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{U}_{i,j-k+2}^m+\sum _{k=0}^{N_y-j+3}{\omega }_{k,\gamma }^{\left({\alpha }_{i,j}^m\right)}{U}_{i,j+k-2}^m\right]=U_{ij}^{m-1}+\tau f_{ij}^{m-1}.\hfill \end{array}$$

(33)

Now, we take into account the matrix form of the discretization linear system (33). Denote

$$U^m={(U_{1,1}^m,U_{2,1}^m,\cdots ,U_{N_x-1,1}^m,\cdots ,U_{1,N_y-1}^m,\cdots ,U_{N_x-1,N_y-1}^m)}^T.$$

For $i=1,2,\dots ,N_x-1$, $j=1,2\cdots ,N_y-1$, let $f_{ij}^m=f(U_{ij}^m,x_i,y_j,t_m)$ and then define

$$F^m={(f_{1,1}^m,f_{2,1}^m,\cdots ,f_{N_x-1,1}^m,\cdots ,f_{1,N_y-1}^m,\cdots ,f_{N_x-1,N_y-1}^m)}^T.$$

For $j=1,2,\dots ,N_y-1$, denote

$$D_{x,j}^m=\mathrm{diag}\left({\displaystyle \frac{\tau K_{C\alpha ,1,j}^m}{h_x^{{\alpha }_{1,j}^m}}},{\displaystyle \frac{\tau K_{C\alpha ,2,j}^m}{h_x^{{\alpha }_{2,j}^m}}},\cdots ,{\displaystyle \frac{\tau K_{C\alpha ,N_x-1,j}^m}{h_x^{{\alpha }_{N_x-1,j}^m}}}\right),$$

$$D_{y,j}^m=\mathrm{diag}\left({\displaystyle \frac{\tau K_{C\beta ,1,j}^m}{h_y^{{\beta }_{1,j}^m}}},{\displaystyle \frac{\tau K_{C\beta ,2,j}^m}{h_y^{{\beta }_{2,j}^m}}},\cdots ,{\displaystyle \frac{\tau K_{C\beta ,N_x-1,j}^m}{h_y^{{\beta }_{N_x-1,j}^m}}}\right).$$

Then define

$$D_x^m=\mathrm{diag}\left(D_{x,1}^m,D_{x,2}^m,\dots ,D_{x,N_y-1}^m\right),D_y^m=\mathrm{diag}\left(D_{y,1}^m,D_{y,2}^m,\dots ,D_{y,N_y-1}^m\right).$$

Similar to Eq. (21), let $W_{x,j}^m={\left[w_{x,i_1,i_2}^{m,j}\right]}_{i_1,i_2=1}^{N_x-1}$, where

$$w_{x,i_1,i_2}^{m,j}=\left\{\begin{array}{cc}2{\omega }_{2,\gamma }^{\left({\alpha }_{i_1,j}^m\right)},\hfill & i_2=i_1,\hfill \\ {\omega }_{1,\gamma }^{\left({\alpha }_{i_1,j}^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_{i_1,j}^m\right)},\hfill & i_2=i_1\pm 1,\hfill \\ {\omega }_{0,\gamma }^{\left({\alpha }_{i_1,j}^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_{i_1,j}^m\right)},\hfill & i_2=i_1\pm 2,\hfill \\ {\omega }_{|i_1-i_2|+2,\gamma }^{\left({\alpha }_{i_1,j}^m\right)},\hfill & |i_1-i_2|>2.\hfill \end{array}\right.$$

(34)

Then let

$$W_x^m=\mathrm{diag}\left(W_{x,1}^m,W_{x,2}^m,\dots ,W_{x,N_y-1}^m\right).$$

Let $W_{y,i}^m={\left[w_{y,j_1,j_2}^{m,i}\right]}_{j_1,j_2=1}^{N_y-1}$, where

$$w_{y,j_1,j_2}^{m,i}=\left\{\begin{array}{cc}2{\omega }_{2,\gamma }^{\left({\beta }_{i,j_1}^m\right)},\hfill & j_2=j_1,\hfill \\ {\omega }_{1,\gamma }^{\left({\beta }_{i,j_1}^m\right)}+{\omega }_{3,\gamma }^{\left({\beta }_{i,j_1}^m\right)},\hfill & j_2=j_1\pm 1,\hfill \\ {\omega }_{0,\gamma }^{\left({\beta }_{i,j_1}^m\right)}+{\omega }_{4,\gamma }^{\left({\beta }_{i,j_1}^m\right)},\hfill & j_2=j_1\pm 2,\hfill \\ {\omega }_{|j_1-j_2|+2,\gamma }^{\left({\beta }_{i,j_1}^m\right)},\hfill & |j_1-j_2|>2.\hfill \end{array}\right.$$

(35)

Let

$$P=\left({\overrightarrow{\ell }}_1,{\overrightarrow{\ell }}_{1+(N_y-1)},\dots ,{\overrightarrow{\ell }}_{1+(N_x-2)(N_y-1)},\dots ,{\overrightarrow{\ell }}_{N_y-1},{\overrightarrow{\ell }}_{2(N_y-1)},\dots ,{\overrightarrow{\ell }}_{(N_x-1)(N_y-1)}\right),$$

where ${\overrightarrow{\ell }}_j$ mean the jth column of the identity square matrix and the matrix belongs to $(N_x-1)(N_y-1)$. Then let

$$W_y^m=P^T\mathrm{diag}\left(W_{y,1}^m,W_{y,2}^m,\dots ,W_{y,N_x-1}^m\right)P.$$

Using the above notations, we obtain the linear system of (33):

$$\begin{array}{c}\hfill \left(I-D_x^m{W}_x^m-D_y^m{W}_y^m\right)U^m=U^{m-1}+\tau F^{m-1},m=1,2,\dots ,M.\end{array}$$

(36)

The following theoretical results require the conditions:

$$\left\{\begin{array}{c}1<{\alpha }_{ij}^m\le {\alpha }^{*},{\tilde{\gamma }}_3\left({\alpha }_{ij}^m\right)\le {\gamma }_{x,ij}^m\le min\{{\tilde{\gamma }}_2\left({\alpha }_{ij}^m\right),{\tilde{\gamma }}_4\left({\alpha }_{ij}^m\right)\};\\ 1<{\beta }_{ij}^m\le {\alpha }^{*},{\tilde{\gamma }}_3\left({\beta }_{ij}^m\right)\le {\gamma }_{y,ij}^m\le min\{{\tilde{\gamma }}_2\left({\beta }_{ij}^m\right),{\tilde{\gamma }}_4\left({\beta }_{ij}^m\right)\};\end{array}\right.$$

(37)

for $i=1,2,\dots ,N_x-1,j=1,2,\dots ,N_y-1,$ where ${\alpha }^{*}\approx 1.83384978$ is given by (26), and ${\tilde{\gamma }}_2(·),{\tilde{\gamma }}_3(·)$ and ${\tilde{\gamma }}_4(·)$ are given by (25).

Lemma 8.

Suppose (37) holds, then

$${∥{\left(I-D_x^m{W}_x^m-D_y^m{W}_y^m\right)}^{-1}∥}_{\infty }\le 1.$$

Proof. 

Based on the given conditions, $I-D_x^m{W}_x^m-D_y^m{W}_y^m$ is an L-matrix. We get the sum of $((j-1)(N_x-1)+i)$th row of $I-D_x^m{W}_x^m-D_y^m{W}_y^m$ is

$$\begin{array}{ccc}& & 1-{\displaystyle \frac{\tau K_{C\alpha ,i,j}^m}{h_x^{{\alpha }_{i,j}^m}}}\sum _{j=1}^{(N_x-1)(N_y-1)}{\left[W_x^m\right]}_{ij}-{\displaystyle \frac{\tau K_{C\beta ,i,j}^m}{h_y^{{\beta }_{i,j}^m}}}\sum _{j=1}^{(N_x-1)(N_y-1)}{\left[W_y^m\right]}_{ij}\hfill \\ & \ge & 1-{\displaystyle \frac{\tau K_{C\alpha ,i,j}^m}{h_x^{{\alpha }_{i,j}^m}}}\left[2{\omega }_{2,\gamma }^{\left({\alpha }_{i,j}^m\right)}+2({\omega }_{1,\gamma }^{\left({\alpha }_{i,j}^m\right)}+{\omega }_{3,\gamma }^{\left({\alpha }_{i,j}^m\right)})+2({\omega }_{0,\gamma }^{\left({\alpha }_{i,j}^m\right)}+{\omega }_{4,\gamma }^{\left({\alpha }_{i,j}^m\right)})+2\sum _{j=5}^{\infty }{\omega }_{j,\gamma }^{\left({\alpha }_{i,j}^m\right)}\right]\hfill \\ & & -{\displaystyle \frac{\tau K_{C\beta ,i,j}^m}{h_x^{{\beta }_{i,j}^m}}}\left[2{\omega }_{2,\gamma }^{\left({\beta }_{i,j}^m\right)}+2({\omega }_{1,\gamma }^{\left({\beta }_{i,j}^m\right)}+{\omega }_{3,\gamma }^{\left({\beta }_{i,j}^m\right)})+2({\omega }_{0,\gamma }^{\left({\beta }_{i,j}^m\right)}+{\omega }_{4,\gamma }^{\left({\beta }_{i,j}^m\right)})+2\sum _{j=5}^{\infty }{\omega }_{j,\gamma }^{\left({\beta }_{i,j}^m\right)}\right]\hfill \\ & \ge & 1-{\displaystyle \frac{2\tau K_{C\alpha ,i,j}^m}{h_x^{{\alpha }_{i,j}^m}}}\sum _{j=0}^{\infty }{\omega }_{j,\gamma }^{\left({\alpha }_{i,j}^m\right)}-{\displaystyle \frac{2\tau K_{C\beta ,i,j}^m}{h_x^{{\beta }_{i,j}^m}}}\sum _{j=0}^{\infty }{\omega }_{j,\gamma }^{\left({\beta }_{i,j}^m\right)}\hfill \\ & \ge & 1.\hfill \end{array}$$

Therefore, $(I-D_x^m{W}_x^m-D_y^m{W}_y^m)\mathbf{e}\ge \mathbf{e}>\mathbf{0}.$ It follows that $I-D_x^m{W}_x^m-D_y^m{W}_y^m$ is an M-matrix. Furthermore, ${∥{\left(I-D_x^m{W}_x^m-D_y^m{W}_y^m\right)}^{-1}∥}_{\infty }\le 1$ by Lemma 4. □

Theorem 9.

Suppose (37) holds, then the difference scheme (36) is stable.

Proof. 

Define $E^m=U^m-{\tilde{U}}^m$. Then we have

$$\begin{array}{c}\hfill (I-D_x^m{W}_x^m-D_y^m{W}_y^m)E^m=E^{m-1}+\tau (F^{m-1}-{\tilde{F}}^{m-1}),m=1,2,\dots ,M.\end{array}$$

That is

$$\begin{array}{c}\hfill E^m={(I-D_x^m{W}_x^m-D_y^m{W}_y^m)}^{-1}{E}^{m-1}+\tau {(I-D_x^m{W}_x^m-D_y^m{W}_y^m)}^{-1}(F^{m-1}-{\tilde{F}}^{m-1}).\end{array}$$

By Lemma 8, we have

$$\begin{array}{c}\hfill \parallel E^m{\parallel }_{\infty }\le (1+\tau L)\parallel E^{m-1}{\parallel }_{\infty }\le {(1+\tau L)}^m\parallel E^0{\parallel }_{\infty }\le e^{m\tau L}\parallel E^0{\parallel }_{\infty }\le e^{TL}{\parallel E^0\parallel }_{\infty }.\end{array}$$

The conclusion of the theorem follows. □

Theorem 10.

Let the exact solution to IBV problem (31) at the grid points and the solution of (36) be denoted by $u_{ij}^m$ and $U_{ij}^m$, respectively. Suppose (37) holds, then for $m=1,2,\cdots ,M$, it holds

$$\parallel u^m-U^m{\parallel }_{\infty }\le CT(\tau +h_x^l+h_y^l),$$

where $C>0,$ $l\in \{3,4\}$.

Proof. 

Denote $e^m=u^m-U^m$. Utilizing Eqs. (31) and (36), we get

$$(I-D_x^m{W}_x^m-D_y^m{W}_y^m)e^m=e^{m-1}+\tau (F_e^{m-1}-F^{m-1})+\tau R^m,m=1,2,\dots ,M,$$

where $\parallel R^m{\parallel }_{\infty }\le C(\tau +h_x^l+h_y^l),l\in \{3,4\}.$

That is,

$$e^m={(I-D_x^m{W}_x^m-D_y^m{W}_y^m)}^{-1}{e}^{m-1}+\tau {(I-D_x^m{W}_x^m-D_y^m{W}_y^m)}^{-1}((F_e^{m-1}-F^{m-1})+R^m).$$

Notice that $e^0=0,$ and from Lemmas 6 and 8, we get

$$\begin{array}{c}\hfill \parallel e^m{\parallel }_{\infty }\le (1+\tau L)\parallel e^{m-1}{\parallel }_{\infty }+\tau {\parallel R^m\parallel }_{\infty }\le e^{m\tau L}Cm\tau (\tau +h_x^l+h_y^l)\le C{e}^{TL}T(\tau +h_x^l+h_y^l).\end{array}$$

The conclusion of the theorem follows. □

Remark 3.

For preconditioners, we extend 1D case to 2D case. Let

$$W_{x,j}^m=W_{x,j,k_b}^m+(W_{x,j}^m-W_{x,j,k_b}^m),j=1,2,\cdots ,Ny-1,$$

$$W_{y,i}^m=W_{y,i,k_b}^m+(W_{y,i}^m-W_{y,i,k_b}^m),i=1,2,\cdots ,Nx-1,$$

where

$$W_{x,j,k_b}^m={\tilde{W}}_{x,j,k_b}^m+\mathrm{diag}\left(\mathrm{sum}(W_{x,j,}^m-{\tilde{W}}_{x,j,k}^m,2)\right),W_{y,i,k_b}^m={\tilde{W}}_{y,i,k_b}^m+\mathrm{diag}\left(\mathrm{sum}(W_{y,i,}^m-{\tilde{W}}_{y,i,k}^m,2)\right),$$

wich ${\tilde{W}}_{x,j,k_b}^m$ and ${\tilde{W}}_{y,i,k_b}^m$ being similar to (30). Let

$$W_{x,k_b}^m=\mathrm{diag}\left(W_{x,1,k_b}^m,W_{x,2,k_b}^m,\dots ,W_{x,N_y-1,k_b}^m\right).$$

$$W_{y,k_b}^m=P^T\mathrm{diag}\left(W_{y,1,k_b}^m,W_{y,2,k_b}^m,\dots ,W_{y,N_x-1,k_b}^m\right)P.$$

Then we use

$$P_{k_b}^m=I-D_x^m{W}_{x,k_b}^m-D_y^m{W}_{y,k_b}^m$$

or

$${\tilde{P}}_{k_b}^m=I-D_x^m{\tilde{W}}_{x,k_b}^m-D_y^m{\tilde{W}}_{y,k_b}^m$$

as a preconditioner.

Spectral analysis of preconditioners is similar to Section 4.

6. Numerical Results

We verify the validity of the developed IE-4WSGD schemes with two examples though numerical experiments in this section. We utilize Gauss elimination (GE) and the preconditioned GMRES (PGMRES) method to solve the relevant linear systems. Moreover, their CPU time is compared. The CPU time cover the whole solving process, involving the original matrix, corresponding preconditioner, and solving linear systems. For the restarted PGMRES method, we take the initial guess:

$$U^{m+1,0}=\left\{\begin{array}{cc}{U}^m,\hfill & m=0,\hfill \\ 2{\widehat{U}}^m-{\widehat{U}}^{m-1},\hfill & m>0.\hfill \end{array}\right.$$

The stopping criterion is

$${\displaystyle \frac{\parallel {{\mathbf{F}}_{\mathbf{r}}}_j{\parallel }_2}{\parallel {{\mathbf{F}}_{\mathbf{b}}}^{m+1}{\parallel }_2}}<{10}^{-5},$$

where ${{\mathbf{F}}_{\mathbf{r}}}_j$ represents the residual vector after j iterations and ${{\mathbf{F}}_{\mathbf{b}}}^{m+1}$ represents the right-sided vector of the corresponding linear equations.

For 2D problem, we take $N_x=N_y=N$. We use a parenthesis (Iter, CPU time) to represent the average number of iterations for solving discrete linear systems (the first component) and the CPU time in seconds (the second component). For the Gauss elimination (GE) method, the first component is simply replaced by “–”. The order of spatial accuracy is denoted by ${\mathrm{“rate}}_h$”, which is given by

$${rate}_h={\displaystyle \frac{{log}_2(\mathrm{Error}(N_1,M_1)/\mathrm{Error}(N_2,M_2))}{{log}_2(N_2/N_1)}}$$

with $M_i=O\left(N_i^4\right)$ and

$$\mathrm{Error}(N,M)={∥{\tilde{\mathbf{u}}}^M-{\mathbf{u}}_T∥}_{\infty },$$

where ${\tilde{\mathbf{u}}}^M$ denote the approximate solution vector at the final time step and ${\mathbf{u}}_T$ denote the exact solution vector at $t=T$. All computations are performed via Matlab version 2021a on a DESKTOP-1VC6MHV (Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz 2.30GHz, 12.0GB RAM).

Example 1.

Consider the IBV problem for a 1D SRVONFDE with $a=0$, $b=2$, $T=1$, $\alpha (x,t)={\displaystyle \frac{7}{5}}+{\displaystyle \frac{1}{2}}{e}^{-x^2{t}^2-1}$ or ${\displaystyle \frac{545-{\left(xt\right)}^4}{300}}$, and $K(x,t)=-2{(x+1)}^2{e}^{x}cos{\displaystyle \frac{\pi \alpha (x,t)}{2}}$. Let $u(x,t)=e^{-t}{x}^4{(2-x)}^4$, we can get

$${\displaystyle \frac{\partial u(x,t)}{\partial t}}=-e^{-t}{x}^4{(2-x)}^4,{\displaystyle \frac{{\partial }^{\alpha (x,t)}u(x,t)}{\partial {\left|x\right|}^{\alpha (x,t)}}}=h(x,t),$$

where

$$h(x,t)=-e^{-t}{\displaystyle \frac{sec\left(\frac{\pi \alpha (x,t)}{2}\right)}{2}}\left(\sum _{i=5}^9{\displaystyle \frac{{(-1)}^{i-1}{2}^{9-i}{C}_4^{i-5}\mathsf{\Gamma }\left(i\right)}{\mathsf{\Gamma }(i-\alpha (x,t\left)\right)}}\left(x^{i-1-\alpha (x,t)}+{(2-x)}^{i-1-\alpha (x,t)}\right)\right).$$

Let

$$f(u,x,t)=-e^{-t}{x}^4{(2-x)}^4-K(x,t)h(x,t)+u^3-e^{-3t}{x}^{12}{(2-x)}^{12}.$$

It can be checked that $u(x,t)=e^{-t}{x}^4{(2-x)}^4$ is the exact solution of the relevant IBV problem.

Table 1 and Table 2 demonstrate the spatial convergence rate for Example 1. We take different variable order functions ${\alpha }_1(x,t)=1.4+0.5e^{-{\left(xt\right)}^2-1}\in (1,{\alpha }^{**})$ and ${\alpha }_2(x,t)={\displaystyle \frac{545-{\left(xt\right)}^4}{300}}\in (1,{\alpha }^{*})$, where ${\alpha }^{*}\approx 1.83384978$ and ${\alpha }^{**}\approx 1.79968326$ are given by (26) and (27) respectively. Table 1 shows that the IE-4WSGD scheme has spatial fourth accuracy. For ${\alpha }_2(x,t)={\displaystyle \frac{545-{\left(xt\right)}^4}{300}}$, we choose ${\gamma }_i^m$ by using (29). From Table 2, it can also be seen that the IE-4WSGD scheme achieve spatial fourth accuracy for ${\alpha }_2(x,t)$. In fact, although some values of ${\alpha }_2(x,t)$ are larger than ${\alpha }^{**}$, the difference ${\alpha }_2(x,t)-{\alpha }^{**}$ is small, which leads to high order spatial accuracy. For the PGMRES method, we choose $k_b=5$. We can see that the PGMRES methods require less CPU times than the GE method for sufficiently large N. Although the average number of iterations between preconditioners $P_{k_b}^m$ and ${\tilde{P}}_{k_b}^m$ is barely no difference for large N, the former requires more CPU times than the latter because of higher costs of constructing $P_{k_b}^m$ at each time step. Note that for the $N=2^6,M=2^2$ and $N=2^7,M=2^6$, the PGMRES method with ${\tilde{P}}_{k_b}^m$ need much more iterations than the one with $P_{k_b}^m$, which indicates that $P_{k_b}^m$ is more efficient in the sense of the number of iterations.

Example 2.

Consider the IBV problem for 2D SRVONFDE (31) with following conditions:

$$\begin{array}{ccc}& & (x,y,t)\in (1,3)\times (2,4)\times (0,1),u(x,y,0)=0,\hfill \\ & & K_{\alpha }(x,y,t)=-10e^{xy+1}cos\left({\displaystyle \frac{\pi \alpha (x,y,t)}{2}}\right),K_{\beta }(x,y,t)=-16e^{xy+2}cos\left({\displaystyle \frac{\pi \beta (x,y,t)}{2}}\right),\hfill \end{array}$$

and

$$\begin{array}{c}\hfill f(u,x,y,t)=cos\left(t\right){(x-1)}^4{(3-x)}^4{(y-2)}^4{(4-y)}^4+(e^{sin\left(t\right){(x-1)}^4{(3-x)}^4{(y-2)}^4{(4-y)}^4}-e^u)+h(x,y,t)\end{array}$$

with

$$\begin{array}{ccc}& & h(x,y,t)\hfill \\ & =& {\displaystyle \frac{sin\left(t\right)K_{\alpha }(x,y,t)}{2cos\left(\pi \alpha \right(x,y,t)/2)}}\hfill \\ & & \times \left(\sum _{i=5}^9{\displaystyle \frac{{(-1)}^{i-1}{C}_4^{i-5}\mathsf{\Gamma }\left(i\right)}{\mathsf{\Gamma }(i-\alpha (x,y,t\left)\right)}}\left({(x-1)}^{i-1-\alpha (x,y,t)}+{(3-x)}^{i-1-\alpha (x,y,t)}\right)\right){(y-2)}^4{(4-y)}^4\hfill \\ & & +{\displaystyle \frac{sin\left(t\right)K_{\beta }(x,y,t)}{2cos\left(\pi \beta \right(x,y,t)/2)}}\hfill \\ & & \times \left(\sum _{i=5}^9{\displaystyle \frac{{(-1)}^{i-1}{C}_4^{i-5}\mathsf{\Gamma }\left(i\right)}{\mathsf{\Gamma }(i-\beta (x,y,t\left)\right)}}\left({(y-2)}^{i-1-\beta (x,y,t)}+{(4-y)}^{i-1-\beta (x,y,t)}\right)\right){(x-1)}^4{(3-x)}^4.\hfill \end{array}$$

Then the exact solution is

$$u(x,y,t)=sin\left(t\right){(x-1)}^4{(3-x)}^4{(y-2)}^4{(4-y)}^4.$$

The order of the fractional orders are specified in Table 3 and Table 4.

We observe the spatial accuracy of the IE-4WSGD scheme for the 2D problem (31). In Table 3 and Table 5, the orders $\alpha (x,y,t)$, $\beta (x,y,t)$ are in $(1,{\alpha }^{**})$, so we can choose the corresponding parameters to achieve fourth order spatial accuracy. Numerical results in Table 3 and Table 4 verify that the derived scheme has 4th order accuracy. In Table 4 and Table 6, the orders $\alpha (x,y,t)$, $\beta (x,y,t)$ are in $(1,{\alpha }^{*})$. Similar to one-dimensional example, numerical results illustrate that the proposed scheme can achieve theoretical accuracy presented in Section 5 for $\alpha (x,y,t)$ and $\beta (x,y,t)$ that satisfy the conditions. For PGMRES method, we choose $k_b=3$ in this example. Similar to Example 1, PGMRES methods are more efficient than the GE method for large N. Moreover, the PGMRES method with ${\tilde{P}}_{k_b}^m$ need sightly more the average number of iterations than the one with $P_{k_b}^m$ while the former requires more CPU time.

7. Concluding Remarks

In this paper, we proposed 4WSGD schemes for the Riesz VO fractional derivative, and derived the convergence and stability of IE-4WSGD schemes for 1D and 2D SRVONFDEs. Numerical experiments show that the presented schemes and the PGMRES methods are very efficient. In the future, we will study the application of the schemes and algorithms for solving other SRVONFDEs, e.g., the movement of groundwater pollution and control.