The Finite Difference Method and Analysis for Simulating the Unsteady Generalized Maxwell Fluid with a Multi-Term Time Fractional Derivative

Abstract

The finite difference method is used to solve a new class of unsteady generalized Maxwell fluid models with multi-term time-fractional derivatives. The fractional order range of the Maxwell model index is from 0 to 2, which is hard to approximate with general methods. In this paper, we propose a new finite difference scheme to solve such problems. Based on the discrete $H^1$ norm, the stability and convergence of the considered discrete scheme are discussed. We also prove that the accuracy of the method proposed in this paper is $O(\tau +h^2)$. Finally, some numerical examples are provided to further demonstrate the superiority of this method through comparative analysis with other algorithms.

1. Introduction

In the past several years, non-Newtonian fluids have received widespread attention from mathematicians and physicists. Compared with Newtonian fluids, non-Newtonian fluids were widely used in industry due to their complexity and special properties [1,2]. Viscoelastic fluid is a special type of non-Newtonian fluid, as the name suggests, that possesses both the viscous properties of a liquid and the elastic properties of a solid. It can produce an instantaneous elastic response and continuous stress when a load is applied to it. Due to the widespread applications of viscoelastic fluids in medical [3], chemical [4], biological [5], physical [6], and other fields, various models have been proposed to study their properties [7,8,9]. The Maxwell liquid model is one of them, which is composed of multiple continuous Maxwell units in parallel. Due to its potential advantage of not relying on the complex effects of fluid viscosity, it has attracted the attention of a large number of researchers in recent years, and many works have been published. The heat transfer behavior of Maxwell fluids in a stretching oscillating wall channel during motion was studied under the Cattaneo–Christov heat flux model, and the influence of the parameters on fluid behavior was considered with the help of the homotopy deformation method and series analysis method [10]. The Cattaneo–Christov double diffusion theory model was applied to analyze the thermal effects and transport flow phenomena of Maxwell fluids based on it [11]. The dynamics of Maxwell fluids in rectangular containers under periodic pressure with simple harmonic motion characteristics were studied with the finite difference method, and thus the trend of resonance change was obtained [12]. The MHD flow of Maxwell fluid on its surface was discussed when the carrier medium had certain stratification and thermophoresis effects. The experimental results showed that stratification affects the temperature change rate of the fluid [13]. Compared with these research results, a more complex situation was discussed in [14]. Taking into account the characteristics of various flow parameters, six second-order partial differential equations were used to conduct theoretical and numerical analysis of Maxwell fluids with non-homogeneous boundary conditions. At the same time, on the basis of [14], the scholar in [15] considered the influence of the Joule heating phenomenon and Lorentz force on the experimental results, although the flow of Maxwell fluids still occurs in porous interstitial spaces.

However, it is regrettable that the above article only models the dynamic behavior of fluids in the sense of integer-order partial differential concepts. In practical industrial applications, fluid behavior is much more complex than we can imagine, and integer-order models cannot effectively describe the impact of the flow rate on stress. Fortunately, fractional calculus has become a powerful tool, and it has been several hundred years since its concept was first introduced [16]. It allows for noninteger order differential operations, which could be found in the communication between L’hospital and Leibniz back then. In recent years, fractional calculus has received extensive research in describing complex nonlinear systems due to its unique genetic and memory effects [17,18,19,20]. Fractional calculus is considered by many scholars as an effective means of describing the rheological effects of viscoelastic fluids. In this field, there are several high-quality articles worth exploring. The work in [21] may be one of the earliest representative works in this field, as the author systematically introduced the theoretical calculation process of fractional-order Maxwell models with the classical discrete Laplace transform. Three years later, the scope of application of this method was expanded, and analytical solutions were obtained using the Weber transform and Laplace transform under different flow conditions in [22]. In [23], the Laplace transform and Fourier transform were used to deal with the complex behavior of Maxwell fluids between vertical sidewalls. Theoretical studies have shown that integer-order models are the limited case of the generalized Maxwell model. The authors of [24] improved the method in [23] by mixing the Hankel transform and Laplace transform techniques, effectively reducing the difficulty of theoretical derivation. The results showed that parameters such as the order and Prandtl number have a significant impact on heat transfer. In [25], the process of using Maxwell fluids to clean residual oil from dead corners was studied, which meant that the stress genetic effect of viscoelastic fluid had to be taken into account. The semi-analytical solution was highly consistent with the image simulation results. Akyildiz and Siginer discussed the seepage flow phenomenon of Maxwell fluids in a triangular region and obtained an accurate solution to the momentum equation under the Caputo–Fabrizio definition [26]. At the same time, the (semi-) analytical solutions under the Caputo–Fabrizio and Atangana–Baleanu definition were obtained using the convolution theorem and Laplace transform method in [27], and the drawback of the relatively single-fractional order definition was improved in [26]. In [28], the unstable flow of Maxwell fluids under two different definitions in a heating environment was considered, and the research result included the integer order (fractional-order parameter was one) as a special case.

Nevertheless, almost all of the above articles are devoted to the study of analytical solutions, and the methods include the Laplace transform, Fourier transform, and other integral transformation methods. These methods are only suitable for some linear Maxwell fluids. But, in the actual simulation of heat transfer and the viscoelastic fluid flow model, because of the unpredictable complexity and the influence of physical factors, we will face nonlinear models in most cases. Seeking numerical solutions as an effective approach, there has been some effective numerical methods proposed, such as the finite difference method, finite element method, and spectral method. Interested readers can refer to [29,30,31,32] and the cited articles therein. Nevertheless, they all have some drawbacks. For example, in [29,30], the authors only considered fractional derivatives in the temporal dimension and neglected the spatial dimension. The authors of [31,32] both paid attention to the spatiotemporal dimension. Their algorithms produced significant truncation errors, which were $\mathcal{O}(\Delta t^{2-\alpha }+\Delta x)$. Such accuracy is insufficient in practical problems. Hence, an effective numerical method is important for the Maxwell fluid model with a fractional-order derivative.

In [33], an unsteady flow of a generalized fractional-order Maxwell fluid model was proposed:

$$(1+\lambda D_t^{\alpha })\frac{\partial u(x,t)}{\partial t}=v\frac{{\partial }^{2}u(x,t)}{\partial x^2},$$

(1)

where $v=\mu /\rho$ is the kinematic viscosity of the fluid. In [22], the authors proposed the unsteady rotating flows of a viscoelastic fluid with the fractional Maxwell model between coaxial cylinders as follows:

$$(1+\lambda D_t^{\alpha })\frac{\partial u(x,t)}{\partial t}=v(\frac{{\partial }^{2}u(x,t)}{\partial x^2}+\frac{1}{x}\frac{\partial u(y,t)}{\partial x}-\frac{u(x,t)}{x^2}).$$

(2)

In the following, we mainly study a new multi-term time fractional-order non-Newtonian diffusion model, shown below, which is more general:

$$a_1{D}_t^{\alpha }u(x,t)+a_2\frac{\partial u(x,t)}{\partial t}=a_3\frac{{\partial }^{2}u(x,t)}{\partial x^2}-a_{4}u(x,t)-a_5{D}_t^{\beta }u(x,t)+f(x,t),(x,t)\in \Omega ,$$

(3)

with the boundary conditions

$$u(0,t)=0,u(L,t)=0,0\le t\le T,$$

(4)

and the initial conditions

$$u(x,0)={\varphi }_1\left(x\right),u_t(x,0)={\varphi }_2\left(x\right),0\le x\le L,$$

(5)

where $a_1,a_2,a_3,a_4\ge 0$, $0<\beta <1$, $1<\alpha <2$, and $\Omega =(0,L)\times (0,T)$.

At the beginning of the simulation, the GL definition may generate certain errors, while the RL definition is mainly used for theoretical research of fractional integral equations. The Caputo definition is more relevant to modern physical applications. Meanwhile, the initial conditions for Caputo fractional differential equations are the same as those for integer-order equations. Consequently, we only study the Caputo fractional derivative in this paper. Then, $D_t^{\alpha }u(x,t)$ and $D_t^{\beta }u(x,t)$ are given as follows [34]:

$$D_t^{\alpha }u(x,t)u(x,t)=\frac{1}{\Gamma (2-\alpha )}{\int }_0^t{(t-\tau )}^{1-\alpha }\frac{{\partial }^{2}u(x,\tau )}{\partial {\tau }^2}d\tau ,1<\alpha <2,$$

(6)

and

$$D_t^{\beta }u(x,t)u(x,t)=\frac{1}{\Gamma (1-\beta )}{\int }_0^t{(t-\tau )}^{-\beta }\frac{\partial u(x,\tau )}{\partial \tau }d\tau ,0<\beta <1.$$

(7)

In the following, we will study the application of the finite difference method to these generalized Maxwell fluid models with multi-term fractional derivatives. The basic idea of the finite difference method is to discretize continuous partial differential equations into discrete algebraic equations, numerical solutions are obtained by solving this algebraic equation system. In the finite difference method, the solution region is divided into a series of grid points and regions. Then, the derivative of the original equation at each grid point is approximated. Thus, a discrete system of equations will be reached. A new finite difference scheme is proposed to approximate this model. The stability and convergence analysis of this scheme is proposed based on the discrete $H^1$ norm.

The general fractional order non-Newtonian fluid model not only has a multi-term fractional time derivative, but the range of the fractional order is from 0 to 2, which is challenging to approximate. A new finite difference scheme is proposed to approximate the multi-term fractional order Maxwell model. The main contributions of this paper are as follows. (1) Our method is proposed for the new Maxwell model, which has a multi-term fractional time derivative. We also propose the L2 scheme for the fractional order term $D_t^{\alpha }$ with the first-order accuracy. Meanwhile, we propose some useful and important lemma, and we find that these lemma can be extended to other multi-term fractional-order diffusion models. (2) The stability and convergence analysis of this scheme is proposed based on the discrete $H^1$ norm. We prove that the accuracy of the new finite difference method is $O(\tau +h^2)$. (3) A numerical example is used to illustrate the validity and rationality of the method with different fractional orders. Our method is effective enough to be generalized to solve the generalized fractional-order Oldroyd-B fluid model, the generalized fractional-order Burgers fluid model and other generalized non-Newtonian fluid models.

The structure of this article is as follows. Section 2 introduces some preliminary knowledge of the finite difference method, and we propose the numerical scheme for discretizing the time fractional derivative. In Section 3 and Section 4, the proof of solvability, stability, and convergence of our new numerical algorithm are given by the energy method. Then, the effectiveness of the algorithm is discussed when the time fractional order of the scheme is changed. In Section 5, the simulation results of some examples are given to verify the superiority of the proposed algorithm. Finally, the conclusions and future work are stated in Section 6.

2. Preliminary Knowledge of the Finite-Difference Method

In the following, some notations and properties are given for discretizing the time fractional order derivative.

Lemma 1.

For $1<\alpha <2$, we can define $b_k^{\left(\alpha \right)}={(k+1)}^{2-\alpha }-k^{2-\alpha }$, $k=0,1,2,\cdots$, and then

$$1.b_k^{\left(\alpha \right)}>0,b_0^{\left(\alpha \right)}=1,b_k^{\left(\alpha \right)}>b_{k+1}^{\left(\alpha \right)},$$

(8)

$$2.(2-\alpha ){(k+1)}^{1-\alpha }\le b_k^{\left(\alpha \right)}\le (2-\alpha )k^{1-\alpha },$$

(9)

$$3.\underset{k\to \infty }{lim}{b}_k^{\left(\alpha \right)}=0.$$

(10)

Proof. 

It is easy to prove that $b_k^{\left(\alpha \right)}>0,b_0^{\left(\alpha \right)}=1$. Let $g\left(x\right)=x^{2-\alpha }$ and $f\left(x\right)=g(x+1)-g\left(x\right)$, $1<\alpha <2$. Then, we have

$$f^{\prime }\left(x\right)=g^{\prime }(x+1)-g^{\prime }\left(x\right)=(2-\alpha )(\frac{1}{{(x+1)}^{\alpha -1}}-\frac{1}{{(x+1)}^{\alpha -1}})<0.$$

This indicates that $f\left(x\right)$ is a monotonically decreasing function. Then, we have

$$f\left(k\right)=b_k^{\left(\alpha \right)}>b_{k+1}^{\left(\alpha \right)}=f(k+1).$$

By taking the derivative of the function $g\left(x\right)$, we have

$$g^{\prime }\left(x\right)=(2-\alpha )x^{1-\alpha },$$

and

$$g^{″}\left(x\right)=(2-\alpha )(1-\alpha )x^{-\alpha }.$$

For all $x>0$, we can obtain

$$g^{″}\left(x\right)=(2-\alpha )(1-\alpha )x^{-\alpha }<0,$$

which indicates that $g^{\prime }\left(x\right)$ is a monotonically decreasing function. It follows from Lagrange’s mean value theorem that

$$g(k+1)-g\left(k\right)=g^{\prime }\left(\xi \right),k<\xi <k+1.$$

Since $g^{\prime }\left(x\right)$ is a monotonically decreasing function, we have

$$g^{\prime }(k+1)<g(k+1)-g\left(k\right)<g^{\prime }\left(k\right),$$

In other words, we have

$$(2-\alpha ){(k+1)}^{1-\alpha }\le b_k^{\left(\alpha \right)}\le (2-\alpha )k^{1-\alpha }.$$

This is easily obtained by the squeeze theorem. □

Lemma 2.

For $0<\beta <1$, we also define $b_k^{\left(\beta \right)}={(k+1)}^{1-\beta }-k^{1-\beta }$, $k=0,1,2,\cdots$, and similar properties can be obtained:

$$1.b_k^{\left(\beta \right)}>0,b_0^{\left(\beta \right)}=1,b_k^{\left(\beta \right)}>b_{k+1}^{\left(\beta \right)};$$

(11)

$$2.(1-\beta ){(k+1)}^{-\beta }\le b_k^{\left(\beta \right)}\le (1-\beta )k^{-\beta };$$

(12)

$$3.\underset{k\to \infty }{lim}{b}_k^{\left(\beta \right)}=0;$$

(13)

$$4.b_{k+1}^{\left(\beta \right)}-2b_k^{\left(\beta \right)}+b_{k-1}^{\left(\beta \right)}\ge 0.$$

(14)

Proof. 

The proof is similar to that of Lemma 1. □

To discretize the multi-term time fractional order (Equation (3)), we define the mesh points $x_i=ih$ and $t_n=h\tau$, $i=0,1,\cdots ,M$, $n=0,1,\cdots ,N$, where $\tau =T/N$ and $h=L/M$ are the uniform temporal step size and spatial step size, respectively. Let $u^n={[u_1^n,u_2^n,\cdots ,u_{M-1}^n]}^T$ be the numerical solution of ${\tilde{u}}^n={[u(x_1,t_n),u(x_2,t_n),\cdots ,u(x_{M-1},t_n)]}^T$. The following notations are introduced:

$$\begin{array}{ccc}\hfill & {\nabla }_t{u}_i^n=\frac{u_i^n-u_i^{n-1}}{\tau },& u_i^{n-\frac{1}{2}}=\frac{u_i^n+u_i^{n-1}}{2},\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& {\nabla }_x{u}_i^n=\frac{u_i^n-u_{i-1}^n}{h},& {\delta }_x^2{u}_i^n=\frac{u_{i-1}^n-2u_i^n+u_{i+1}^n}{h^2}.\hfill \end{array}$$

(16)

Then, let $V_h=\left\{v\right|v$ is a grid function on ${\Omega }_h$, and $v_0=v_M=0\}$. We can define the discrete inner products and induced norms for any $v,w\in V_h$ as follows:

$$\begin{array}{cc}\hfill & (v,w)=h{\sum }_{j=1}^{M-1}{v}_i{w}_i,\langle \nabla v,\nabla w\rangle =h{\sum }_{i=1}^M{\nabla }_x{v}_i·{\nabla }_x{w}_i,\end{array}$$

(17)

$$\begin{array}{cc}& {\left|\right|v\left|\right|}_0=\sqrt{(v,v)}{,\left|\right|v\left|\right|}_{\infty }={max}_{0\le i\le M}|v_i{|,|v|}_1=\sqrt{\langle \nabla v,\nabla v\rangle },\end{array}$$

(18)

We can easily check that

$$\begin{array}{cc}\hfill & ({\delta }_x^2{v}_k,v_n)=-\langle {\nabla }_x{v}^k,{\nabla }_x{v}^n\rangle ,\end{array}$$

(19)

$$\begin{array}{cc}& ({\delta }_x^2{v}_k,{\nabla }_t{v}_n)=-\frac{1}{\tau }\langle {\nabla }_x{v}^k,{\nabla }_x{v}^n-{\nabla }_x{v}^{n-1}\rangle .\end{array}$$

(20)

Then, we use the following result in [35] to discretize the time fractional-order derivative $D_t^{\alpha }u(x,t)$$(1<\alpha <2)$. For all $f\left(t\right)\in C^2[0,t_n]$, we have

$$\begin{array}{cc}\hfill & \left|{\int }_0^{t_n}\frac{f^{\prime }\left(t\right)dt}{{(t_n-t)}^{\alpha -1}}-\frac{{\tau }^{1-\alpha }}{2-\alpha }\left[b_0^{\left(\alpha \right)}f\left(t_n\right)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})f\left(t_k\right)-b_{n-1}^{\left(\alpha \right)}f\left(t_0\right)\right]\right|\hfill \\ \hfill & \le C\underset{0\le t\le t_n}{max}\left|f^{″}\left(t\right)\right|{\tau }^{3-\alpha },\hfill \end{array}$$

(21)

where $b_k^{\alpha }={(k+1)}^{2-\alpha }-k^{2-\alpha }$, $k=0,1,2,\cdots$ are defined in Lemma 1.

Due to Equation (21), we can obtain the following equation for $f\left(t\right)\in C^3[0,t_n]$:

$$\begin{array}{cc}\hfill & \left|{\int }_0^{t_n}\frac{f^{″}\left(t\right)dt}{{(t_n-t)}^{\alpha -1}}-\frac{{\tau }^{1-\alpha }}{2-\alpha }\left[b_0^{\left(\alpha \right)}{f}^{\prime }\left(t_n\right)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})f^{\prime }\left(t_k\right)-b_{n-1}^{\left(\alpha \right)}{f}^{\prime }\left(t_0\right)\right]\right|\hfill \\ \hfill & \le C\underset{0\le t\le t_n}{max}\left|f^{‴}\left(t\right)\right|{\tau }^{3-\alpha }.\hfill \end{array}$$

(22)

If we approximate $f^{\prime }\left(t\right)={▽}_{t}f\left(t_n\right)+C_1\tau$, hten we can obtain

$$\begin{array}{cc}\hfill & \frac{{\tau }^{1-\alpha }}{2-\alpha }\left[b_0^{\left(\alpha \right)}{C}_1\tau -\sum _{k=1}^{n-1}\left(b_{n-k-1}^{\left(\alpha \right)}\right)C_1\tau \right]=\frac{C_1{\tau }^{2-\alpha }}{2-\alpha }\left[b_0^{\left(\alpha \right)}-\sum _{k=1}^{n-1}\left(b_{n-k-1}^{\left(\alpha \right)}\right)\right]\hfill \\ & =\frac{C_1{\tau }^{2-\alpha }}{2-\alpha }{b}_{n-1}^{\alpha }\le C_1{\tau }^{2-\alpha }{(n-1)}^{1-\alpha }\le C_1{T}^{1-\alpha }\tau \le C_2\tau .\hfill \end{array}$$

(23)

It follows from Equations (22) and (23) that

$$\begin{array}{cc}\hfill & \left|{\int }_0^{t_n}\frac{f^{″}\left(t\right)dt}{{(t_n-t)}^{\alpha -1}}-\frac{{\tau }^{1-\alpha }}{2-\alpha }\left[b_0^{\left(\alpha \right)}{▽}_{t}f\left(t_n\right)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)}){▽}_{t}f\left(t_k\right)-b_{n-1}^{\left(\alpha \right)}{f}^{\prime }\left(t_0\right)\right]\right|\hfill \\ & \le C\underset{0\le t\le t_n}{max}\left|f^{‴}\left(t\right)\right|{\tau }^{3-\alpha }+C_2\tau \le C_3\tau .\hfill \end{array}$$

(24)

Then, the new discrete scheme for the fractional-order derivative $D_t^{\alpha }u(x,t)$$(1<\alpha <2)$ is obtained at mesh points $(x_i,t_n)$:

$$D_t^{\alpha }u(x_i,t_n)=\frac{{\tau }^{1-\alpha }}{\Gamma (3-\alpha )}[b_0^{\left(\alpha \right)}{\nabla }_{t}u(x_i,t_n)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})u(x_i,t_k)-b_{n-1}^{\left(\alpha \right)}{\varphi }_2\left(x_i\right)]+R_1\left(O\left({\tau }^{3-\alpha }\right)\right),$$

(25)

where $b_k^{\left(\alpha \right)}={(k+1)}^{2-\alpha }-k^{2-\alpha }$ and $|R_1|\le C\tau$. The new discrete scheme (called formula L1) for the fractional-order derivative $D_t^{\beta }(x,t)(0<\beta <1)$ is obtained at mesh points $(x_i,t_n)$:

$$D_t^{\beta }u(x_i,t_n)=\frac{{\tau }^{-\beta }}{\Gamma (2-\beta )}[b_0^{\left(\beta \right)}u(x_i,t_n)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u(x_i,t_k)-b_{n-1}^{\left(\beta \right)}u(x_i,t_0)]+R_2,$$

(26)

where $b_k^{\left(\beta \right)}={(k+1)}^{1-\beta }-k^{1-\beta }$ and $|R_2|\le C\underset{0<t\le T}{max}\left|\frac{{\partial }^{2}u(x,t)}{\partial t^2}\right|({\tau }^{2-\beta }+h^2)$.

Furthermore, the implicit finite difference scheme of Equation (3) is obtained:

$$a_1{D}_t^{\alpha }u(x_i,t_n)+a_2\frac{\partial u(x_i,t_n)}{\partial t}=a_3\frac{{\partial }^{2}u(x_i,t_n)}{\partial x^2}-a_{4}u(x_i,t_n)-a_5{D}_t^{\beta }u(x_i,t_n)+f(x_i,t_n).$$

(27)

It follows from Equations (25) and (26) that

$$\begin{array}{cc}\hfill & \frac{a_1{\tau }^{1-\alpha }}{\Gamma (3-\alpha )}[b_0^{\left(\alpha \right)}{\nabla }_{t}u(x_i,t_n)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})u(x_i,t_k)-b_{n-1}^{\left(\alpha \right)}{\varphi }_2\left(x_i\right)]+a_2{\nabla }_{t}u(x_i,t_n)\hfill \\ \hfill & =a_3{\delta }_x^{2}u(x_i,t_n)-\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}[b_0^{\left(\beta \right)}u(x_i,t_n)-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u(x_i,t_k)-b_{n-1}^{\left(\beta \right)}u(x_i,t_0)]\hfill \\ \hfill & -a_{4}u(x_i,t_n)+f(x_i,t_n)+R_i^n,\hfill \end{array}$$

(28)

where $|R_i^n|\le C(\tau +h^2)$.

If we omit the error term, then one can obtain the implicit finite difference scheme of Equation (3) as follows:

$$\begin{array}{cc}\hfill & \frac{a_1{\tau }^{1-\alpha }}{\Gamma (3-\alpha )}[b_0^{\left(\alpha \right)}{\nabla }_t{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})u_i^k-b_{n-1}^{\left(\alpha \right)}{\varphi }_2\left(x_i\right)]+a_2{\nabla }_t{u}_i^n\hfill \\ \hfill & =a_3{\delta }_x^2{u}_i^n-a_4{u}_i^n-\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}[b_0^{\left(\beta \right)}{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u_i^k-b_{n-1}^{\left(\beta \right)}{u}_i^0]+f_i^n,\hfill \end{array}$$

(29)

which can be recast into

$$\begin{array}{cc}\hfill & c_1{u}_i^n-c_2(u_{i-1}^n-2u_i^n+u_{i+1}^n)=d_1{u}_i^{n-1}+\frac{a_1{\tau }^{1-\alpha }}{\Gamma (3-\alpha )}[\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})u_i^k+b_{n-1}^{\left(\alpha \right)}{\varphi }_2\left(x_i\right)]\hfill \\ \hfill & +\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}[\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u_i^k+b_{n-1}^{\left(\beta \right)}{u}_i^0]+f_i^n,\hfill \end{array}$$

(30)

where $c_1=\frac{a_1{\tau }^{-\alpha }}{\Gamma (3-\alpha )}+\frac{a_2}{\tau }+a_4+\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}$, $c_2=\frac{a_3}{h^2}$ and $d_1=\frac{a_1{\tau }^{-\alpha }}{\Gamma (3-\alpha )}+\frac{a_2}{\tau }$.

Theorem 1.

The finite difference scheme (Equation (30)) for Equation (3) is uniquely solvable.

Proof. 

Due to the difference scheme in Equation (30), the coefficient matrix D is linear tridiagonal at each time level:

$$D=\left(\begin{array}{cccccc}{c}_1+2c_2& -c_2& 0& \cdots & 0& 0\\ -c_2& c_1+2c_2& -c_2& \cdots & 0& 0\\ 0& -b_2& c_1+2c_2& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & c_1+2c_2& -c_2\\ 0& 0& 0& \cdots & -c_2& c_1+2c_2\end{array}\right),$$

(31)

where $c_1=\frac{a_1{\tau }^{-\alpha }}{\Gamma (3-\alpha )}+\frac{a_2}{\tau }+a_4+\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}>0$ and $c_2=\frac{a_3}{h^2}>0$. It follows from Equation (31) that D is not only a strictly diagonally dominant matrix but also a nonsingular matrix, which implies that the finite difference scheme in Equation (30) is uniquely solvable. □

3. Stability Analysis of the Finite Difference Scheme

Lemma 3.

For any $G=\{g_0,g_1,g_2,\cdots ,g_N\}$, it holds that

$$\sum _{n=1}^N\left[b_0{g}_n-\sum _{k=1}^{n-1}(b_{n-k-1}-b_{n-k})g_k-b_{n-1}{g}_0\right]g_n\ge \frac{1}{2}{b}_{N-1}\sum _{n=1}^N{g}_n^2-\frac{1}{2}\sum _{n=1}^N{b}_{n-1}{g}_0^2,$$

(32)

where $b_k$ is defined as $b_k^{\left(\alpha \right)}$ (or $b_k^{\left(\beta \right)}$). If $b_k=b_k^{\left(\alpha \right)}$, then it follows from Lemma 1 that

$$\frac{1}{2}{b}_{N-1}^{\left(\alpha \right)}\sum _{n=1}^N{g}_n^2-\frac{1}{2}\sum _{n=1}^N{b}_{n-1}^{\left(\alpha \right)}{g}_0^2\ge \frac{(2-\alpha )N^{1-\alpha }}{2}\sum _{n=1}^N{g}_n^2-\frac{N^{2-\alpha }}{2}{g}_0^2.$$

(33)

If $b_k=b_k^{\left(\beta \right)}$, then it follows from Lemma 2 that

$$\frac{1}{2}{b}_{N-1}^{\beta }\sum _{n=1}^N{g}_n^2-\frac{1}{2}\sum _{n=1}^N{b}_{n-1}^{\left(\beta \right)}{g}_0^2\ge \frac{(1-\beta )N^{-\beta }}{2}\sum _{n=1}^N{g}_n^2-\frac{N^{1-\beta }}{2}{g}_0^2.$$

(34)

Lemma 4.

Let N be any positive integer and real vector $G={(g_0,g_1,\cdots ,g_N)}^T\in {\mathbb{R}}^{N+1}$. Then, we have

$$\sum _{n=1}^N\left[b_0^{\left(\beta \right)}{g}_n-\sum _{k=0}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k}^{(n-k)})g_k-b_{n-1}^{\left(\beta \right)}{g}_0\right](g_n-g_{n-1})\ge -g_0^2,$$

(35)

where $b_k^{\left(\beta \right)}$ is defined as shown in Lemma 2.

Proof. 

Let ${\omega }_0^n={\omega }_0=1$, ${\omega }_k^n={\omega }_k=b_k^{\left(\beta \right)}-b_{k-1}^{\left(\beta \right)}$, $k=1,2,\cdots ,n-1$, and ${\omega }_n^n=-b_{n-1}^{\left(\beta \right)}$$(n>0)$. Then, Equation (35) can be recast into the form

$$\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k(g_n-g_{n-1})\ge -g_0^2.$$

(36)

Firstly, we will prove

$$\sum _{n=0}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_n-\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_{n-1}\ge 0.$$

(37)

Since

$$\sum _{n=0}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_n=G^T\left(\begin{array}{cccccc}{\omega }_0& 0& 0& \cdots & 0& 0\\ {\omega }_1^1& {\omega }_0& 0& \cdots & 0& 0\\ {\omega }_2^2& {\omega }_1& {\omega }_0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\omega }_{N-1}^{N-1}& {\omega }_{N-2}& {\omega }_{N-3}& \cdots & {\omega }_0& 0\\ {\omega }_N^N& {\omega }_{N-1}& {\omega }_{N-2}& \cdots & {\omega }_1& {\omega }_0\end{array}\right)G,$$

(38)

and

$$\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_{n-1}=G^T\left(\begin{array}{cccccc}{\omega }_1^1& {\omega }_0& 0& \cdots & 0& 0\\ {\omega }_2^2& {\omega }_1& {\omega }_0& \cdots & 0& 0\\ {\omega }_3^3& {\omega }_2& {\omega }_1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\omega }_N^N& {\omega }_{N-1}& {\omega }_{N-2}& \cdots & {\omega }_1& {\omega }_0\\ 0& 0& 0& \cdots & 0& 0\end{array}\right)G,$$

(39)

then

$$\sum _{n=0}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_n-\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_{n-1}=G^{T}WG,$$

(40)

where

$$\begin{array}{c}\hfill W=\left(\begin{array}{cccccc}{\omega }_0-{\omega }_1^1& -{\omega }_0& 0& \cdots & 0& 0\\ {\omega }_1^1-{\omega }_2^2& {\omega }_0-{\omega }_1& -{\omega }_0& \cdots & 0& 0\\ {\omega }_2^2-{\omega }_3^3& {\omega }_1-{\omega }_2& {\omega }_0-{\omega }_1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\omega }_{N-1}^{N-1}-{\omega }_N^N& {\omega }_{N-2}-{\omega }_{N-1}& {\omega }_{N-3}-{\omega }_{N-2}& \cdots & {\omega }_0-{\omega }_1& -{\omega }_0\\ {\omega }_N^N& {\omega }_{N-1}& {\omega }_{N-2}& \cdots & {\omega }_1& {\omega }_0\end{array}\right).\end{array}$$

(41)

It follows from ${\omega }_1^1=-b_0^{\left(\beta \right)}=-{\omega }_0$ and ${\omega }_k^k-{\omega }_{k+1}^{k+1}={\omega }_k$ that

$$\begin{array}{c}\hfill W=\left(\begin{array}{cccccc}2{\omega }_0& -{\omega }_0& 0& \cdots & 0& 0\\ {\omega }_1& {\omega }_0-{\omega }_1& -{\omega }_0& \cdots & 0& 0\\ {\omega }_2& {\omega }_1-{\omega }_2& {\omega }_0-{\omega }_1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\omega }_{N-1}& {\omega }_{N-2}-{\omega }_{N-1}& {\omega }_{N-3}-{\omega }_{N-2}& \cdots & {\omega }_0-{\omega }_1& -{\omega }_0\\ {\omega }_N^N& {\omega }_{N-1}& {\omega }_{N-2}& \cdots & {\omega }_1& {\omega }_0\end{array}\right).\end{array}$$

(42)

Then, proving Equation (37) is equivalent to proving that the matrix W is positive semidefinite. As we know, the quadric form has the property $G^{T}WG=G^T\frac{W^{\prime }+W}{2}G$. Then, we just need to prove that $H=\frac{W^{\prime }+W}{2}$ is positive semidefinite, which means proving that the eigenvalues of H are nonnegative. It follows from Equation (37) that

$$H=\left(\begin{array}{ccccccc}2{\omega }_0& \frac{{\omega }_1-{\omega }_0}{2}& \frac{{\omega }_2}{2}& \frac{{\omega }_3}{2}& \cdots & \frac{{\omega }_{N-1}}{2}& \frac{{\omega }_N^N}{2}\\ \frac{{\omega }_1-{\omega }_0}{2}& {\omega }_0-{\omega }_1& \frac{{\omega }_1-{\omega }_2-{\omega }_0}{2}& \frac{{\omega }_2-{\omega }_3}{2}\cdots & \frac{{\omega }_{N-2}-{\omega }_{N-1}}{2}& \frac{{\omega }_{N-1}}{2}\\ \frac{{\omega }_2}{2}& \frac{{\omega }_1-{\omega }_2-{\omega }_0}{2}& {\omega }_0-{\omega }_1& \frac{{\omega }_1-{\omega }_2-{\omega }_0}{2}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ \frac{{\omega }_{N-1}}{2}& \frac{{\omega }_{N-2}-{\omega }_{N-1}}{2}& \frac{{\omega }_{N-3}-{\omega }_{N-2}}{2}& \frac{{\omega }_{N-4}-{\omega }_{N-3}}{2}& \cdots & {\omega }_0-{\omega }_1& \frac{{\omega }_1-{\omega }_0}{2}\\ \frac{{\omega }_N^N}{2}& \frac{{\omega }_{N-1}}{2}& \frac{{\omega }_{N-2}}{2}& \frac{{\omega }_{N-3}}{2}& \cdots & \frac{{\omega }_1-{\omega }_0}{2}& {\omega }_0\end{array}\right).$$

(43)

Due to Lemma 2, we have

$${\omega }_k^n<0,{\omega }_i^n\le {\omega }_{i+1}^n,k=1,2,\cdots ,n,i=1,2,\cdots ,n-2,$$

(44)

which implies that

$$h_{ii}>0,\forall i=1,2,\cdots ,N+1,$$

(45)

and

$$h_{ij}<0,\forall i\ne j,i,j=1,2,\cdots N+1.$$

(46)

It follows from Equation (42) and Lemma 2 that H is diagonally dominant; in other words, we have

$$h_{ii}\ge \sum _{j=1,j\ne i}^{N+1}\left|h_{ij}\right|,i=1,2,\cdots ,N+1.$$

(47)

Since H is a symmetric matrix, the eigenvalues of H are real numbers. Let ${\lambda }^{*}$ be any eigenvalue of H. It follows from Gershgorin’s theorem that

$$|{\lambda }^{*}-h_{ii}|\le r_i=\sum _{j=1,j\ne i}^{N+1}\left|h_{ij}\right|,$$

(48)

which implies

$$h_{i}i-\sum _{j=1,j\ne i}^{N+1}|h_{ij}|\le {\lambda }^{*}\le h_{ii}+\sum _{j=1,j\ne i}^{N+1}|h_{ij}|.$$

(49)

Then, we have

$${\lambda }^{*}\ge h_{ii}-\sum _{j=1,j\ne i}^{N+1}\left|h_{ij}\right|\ge 0,$$

(50)

which demonstrates that Equation (37) holds; in other words, this means that

$$\sum _{n=0}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_n-\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_{n-1}=\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_n-\sum _{n=1}^N\sum _{k=0}^n{\omega }_{n-k}^n{g}_k{g}_{n-1}+{\omega }_0{g}_0^2\ge 0.$$

(51)

Therefore, we have

$$\sum _{n=1}^N\left[b_0^{\left(\beta \right)}{g}_n-\sum _{k=0}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k}^{(n-k)})g_k-b_{n-1}^{\left(\beta \right)}{g}_0\right](g_n-g_{n-1})\ge -g_0^2.$$

(52)

□

Theorem 2.

The finite difference scheme (Equation (30)) for Equation (3) is unconditionally stable.

Proof. 

We will use the energy method to analyze the stability of the scheme in Equation (30). First, by multiplying Equation (30) by $h\tau {\nabla }_t{u}_i^n$, we sum i from 1 to $M-1$ and sum n from 1 to N, and then we have

$$\begin{array}{cc}\hfill & \frac{a_1{\tau }^{2-\alpha }}{\Gamma (3-\alpha )}\sum _{n=1}^N\sum _{i=1}^{M-1}h[b_0^{\left(\alpha \right)}{\nabla }_t{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})u_i^k-b_{n-1}^{\left(\alpha \right)}{\varphi }_2\left(x_i\right)]{\nabla }_t{u}_i^n\hfill \\ \hfill & +\frac{a_5{\tau }^{1-\beta }}{\Gamma (2-\beta )}\sum _{n=1}^N\sum _{i=1}^{M-1}h[b_0^{\left(\beta \right)}{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u_i^k-b_{n-1}^{\left(\beta \right)}{u}_i^0]{\nabla }_t{u}_i^n\hfill \\ \hfill & +a_2\tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{\left({\nabla }_t{u}_i^n\right)}^2+a_4\tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{u}_i^n{\nabla }_t{u}_i^n=a_3\tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{\delta }_x^2{u}_i^n{\nabla }_t{u}_i^n\hfill \\ \hfill & +\tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{f}_i^n{\nabla }_t{u}_i^n.\hfill \end{array}$$

(53)

It follows from Lemma 3 that

$$\begin{array}{cc}\hfill & \frac{a_1{\tau }^{2-\alpha }}{\Gamma (3-\alpha )}\sum _{n=1}^N\sum _{i=1}^{M-1}h\left[b_0^{\left(\alpha \right)}{\nabla }_t{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})u_i^k-b_{n-1}^{\left(\alpha \right)}{\varphi }_2\left(x_i\right)\right]{\nabla }_t{u}_i^n\hfill \\ \hfill & \ge \frac{a_1{\tau }^{2-\alpha }}{\Gamma (3-\alpha )}\sum _{i=1}^{M-1}h\left[\frac{(2-\alpha )N^{1-\alpha }}{2}\sum _{i=1}^N{\left({\nabla }_t{u}_i^n\right)}^2-\frac{N^{2-\alpha }}{2}{\varphi }_2^2\left(x_i\right)\right]\hfill \\ \hfill & \ge \frac{a_1\tau T^{1-\alpha }}{2\Gamma (2-\alpha )}\sum _{n=1}^N\left|\right|{\nabla }_t{u}^n{\left|\right|}_0^2-\frac{a_1{T}^{2-\alpha }}{2\Gamma (3-\alpha )}\left|\right|{\varphi }_2{\left|\right|}_0^2.\hfill \end{array}$$

(54)

Using Lemma 4, we have

$$\begin{array}{cc}\hfill & \frac{a_5{\tau }^{1-\beta }}{\Gamma (2-\beta )}\sum _{n=1}^N\sum _{i=1}^{M-1}h[b_0^{\left(\beta \right)}{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u_i^k-b_{n-1}^{\left(\beta \right)}{u}_i^0]{\nabla }_t{u}_i^n\hfill \\ \hfill & =\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}\sum _{n=1}^N\sum _{i=1}^{M-1}h[b_0^{\left(\beta \right)}{u}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})u_i^k-b_{n-1}^{\left(\beta \right)}{u}_i^0](u_i^n-u_i^{n-1})\hfill \\ \hfill & \ge -\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}\sum _{i=1}^{M-1}h{\left(u_i^0\right)}^2\hfill \\ \hfill & =-\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}\left|\right|u^0{\left|\right|}_0^2,\hfill \end{array}$$

(55)

and

$$\begin{array}{cc}\hfill a_4\tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{u}_i^n{\nabla }_t{u}_i^n& =a_4\sum _{n=1}^N\sum _{i=1}^{M-1}h{u}_i^n(u_i^n-u_i^{n-1})\hfill \\ \hfill & =a_4\sum _{i=1}^{M-1}h\sum _{n=1}^N{u}_i^n(u_i^n-u_i^{n-1})\hfill \\ \hfill & \ge \frac{a_4}{2}\sum _{i=1}^{M-1}h\sum _{n=1}^N({\left(u_i^n\right)}^2-{\left(u_i^{n-1}\right)}^2)\hfill \\ \hfill & =\frac{a_4}{2}\sum _{i=1}^{M-1}h({\left(u_i^N\right)}^2-{\left(u_i^0\right)}^2)\hfill \\ \hfill & =\frac{a_4}{2}\left|\right|u^N{\left|\right|}_0^2-\frac{a_4}{2}\left|\right|u^0{\left|\right|}_0^2.\hfill \end{array}$$

(56)

Therefore, the left-hand side of Equation (56) is bounded by

$$L\ge \tau (\frac{a_1{T}^{1-\alpha }}{2\Gamma (2-\alpha )}+a_2)\sum _{n=1}^N\left|\right|{\nabla }_t{u}^n{\left|\right|}_0^2+\frac{a_4}{2}\left|\right|u^N{\left|\right|}_0^2-\frac{a_1{T}^{2-\alpha }}{2\Gamma (3-\alpha )}\left|\right|{\varphi }_2{\left|\right|}_0^2-(\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}+\frac{a_4}{2})\left|\right|u^0{\left|\right|}_0^2.$$

(57)

On the other hand, under the relation of the inner and norms (Equation (17)), one can obtain

$$\begin{array}{cc}\hfill a_3\tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{\delta }_x^2{u}_i^n{\nabla }_t{u}_i^n& =a_3\tau \sum _{n=1}^N({\delta }_x^2{u}^n,{\nabla }_t{u}^n)\hfill \\ \hfill & =-a_3\sum _{n=1}^N\langle {\nabla }_x{u}^n,{\nabla }_x{u}^n-{\nabla }_x{u}^{n-1}\rangle \hfill \\ \hfill & \le -\frac{a_3}{2}\sum _{n=1}^N\left(\right||u^n{\left|\right|}_1^2-\left|\right|u^{n-1}{\left|\right|}_1^2)\hfill \\ \hfill & =\frac{a_3}{2}\left(\right||u^0{\left|\right|}_1^2-\left|\right|u^N{\left|\right|}_1^2).\hfill \end{array}$$

(58)

□

Using the inequality $xy\le \epsilon x^2+\frac{y^2}{4\epsilon }$, we can obtain

$$\begin{array}{cc}\hfill & \tau \sum _{n=1}^N\sum _{i=1}^{M-1}h{f}_i^n{\nabla }_t{u}_i^n\hfill \\ \hfill & \le \tau (\frac{a_1{T}^{1-\alpha }}{2\Gamma (2-\alpha )}+a_2)\sum _{n=1}^N\sum _{i=1}^{M-1}h{\left({\nabla }_t{u}_i^n\right)}^2+\frac{\tau }{4(\frac{a_1{T}^{1-\alpha }}{2\Gamma (2-\alpha )}+a_2)}\sum _{n=1}^N\sum _{i=1}^{M-1}h{f}_i^n\hfill \\ \hfill & =\tau (\frac{a_1{T}^{1-\alpha }}{2\Gamma (2-\alpha )}+a_2)\sum _{n=1}^N\left|\right|{\nabla }_t{u}^n{\left|\right|}_0^2+\frac{\tau \Gamma (2-\alpha )}{2(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\left|\right|f^n{\left|\right|}_0^2.\hfill \end{array}$$

(59)

Therefore, the right-hand side of Equation (56) is bounded by

$$\begin{array}{cc}\hfill R& \le \frac{a_3}{2}\left(\right||u^0{\left|\right|}_1^2-\left|\right|u^N{\left|\right|}_1^2)+\tau (\frac{a_1{T}^{1-\alpha }}{2\Gamma (2-\alpha )}+a_2)\sum _{n=1}^N\left|\right|{\nabla }_t{u}^n{\left|\right|}_0^2\hfill \\ \hfill & +\frac{\tau \Gamma (2-\alpha )}{2(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\left|\right|f^n{\left|\right|}_0^2.\hfill \end{array}$$

(60)

It follows from $L=R$ that

$$\begin{array}{c}\hfill \frac{a_4}{2}\left|\right|u^N{\left|\right|}_0^2+\frac{a_3}{2}\left|\right|u^N{\left|\right|}_1^2\le \frac{a_1{T}^{2-\alpha }}{2\Gamma (3-\alpha )}\left|\right|{\varphi }_2{\left|\right|}_0^2+(\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}+\frac{a_4}{2})\left|\right|u^0{\left|\right|}_0^2\\ \hfill +\frac{a_3}{2}\left|\right|u^0{\left|\right|}_1^2+\frac{\tau \Gamma (2-\alpha )}{2(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\left|\right|f^n{\left|\right|}_0^2.\end{array}$$

(61)

Then, this implies that

$$\begin{array}{cc}\hfill \left|\right|u^N{\left|\right|}_0^2& \le \frac{a_1{T}^{2-\alpha }}{a_4\Gamma (3-\alpha )}\left|\right|{\varphi }_2{\left|\right|}_0^2+(1+\frac{2a_5{\tau }^{-\beta }}{a_4\Gamma (2-\beta )})\left|\right|u^0{\left|\right|}_0^2+\frac{a_3}{a_4}\left|\right|u^0{\left|\right|}_1^2\hfill \\ \hfill & +\frac{\tau \Gamma (2-\alpha )}{a_4(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\left|\right|f^n{\left|\right|}_0^2,\hfill \end{array}$$

(62)

and

$$\begin{array}{cc}\hfill \left|\right|u^N{\left|\right|}_1^2& \le \frac{a_1{T}^{2-\alpha }}{a_3\Gamma (3-\alpha )}\left|\right|{\varphi }_2{\left|\right|}_0^2+(\frac{a_4}{a_3}+\frac{2a_5{\tau }^{-\beta }}{a_3\Gamma (2-\beta )})\left|\right|u^0{\left|\right|}_0^2+\left|\right|u^0{\left|\right|}_1^2\hfill \\ \hfill & +\frac{\tau \Gamma (2-\alpha )}{a_3(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\left|\right|f^n{\left|\right|}_0^2.\hfill \end{array}$$

(63)

It follows from Equation (63) and Lemma 4 that

$$\begin{array}{cc}\hfill \left|\right|u^N{\left|\right|}_{\infty }^2& \le \frac{L}{4}\left|\right|u^N{\left|\right|}_1^2\le \frac{L}{4}[\frac{a_1{T}^{2-\alpha }}{a_3\Gamma (3-\alpha )}\left|\right|{\varphi }_2{\left|\right|}_0^2+(\frac{a_4}{a_3}+\frac{2a_5{\tau }^{-\beta }}{a_3\Gamma (2-\beta )})\left|\right|u^0{\left|\right|}_0^2\hfill \\ & +\left|\right|u^0{\left|\right|}_1^2+\frac{\tau \Gamma (2-\alpha )}{a_3(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\left|\right|f^n{\left|\right|}_0^2].\hfill \end{array}$$

(64)

Hence, the inequalities in Equations (61)–(64) demonstrate that the difference scheme in Equation (30) is unconditionally stable.

4. Convergence Analysis of the Finite Difference Scheme

Let ${\tilde{u}}^n={[u(x_1,t_n),u(x_2,t_n),\cdots ,u(x_{M-1},t_n)]}^T$, and $u^n={[u_1^n,u_2^n,\cdots ,u_{M-1}^n]}^T$ be the exact solution and the numerical solution vectors, respectively. The theorem of convergence can be obtained as follows.

Theorem 3.

If the solution to Equation (30) satisfies $u(x,t)\in C_{x,t}^{4,3}(\Omega )$, then there are three positive constants $C_1$, $C_2$, and $C_3$, which are independent of τ and h such that

$$\left|\right|e^n{\left|\right|}_0\le C_1(\tau +h^2),$$

(65)

$$\left|\right|e^n{\left|\right|}_1\le C_2(\tau +h^2),$$

(66)

and

$$\left|\right|e^n{\left|\right|}_{\infty }\le C_3(\tau +h^2).$$

(67)

Proof. 

Denote the errors $e_i^n=u(x_i,t_n)-u_i^n$, $e^n={\tilde{u}}^n-u^n={[e_1^n,e_2^n,\cdots ,e_{M-1}^n]}^T$. By subtracting Equation (29) from Equation (28), we can obtain

$$\begin{array}{cc}\hfill & \frac{a_1{\tau }^{1-\alpha }}{\Gamma (3-\alpha )}\left[b_0^{\left(\alpha \right)}{\nabla }_t{e}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\alpha \right)}-b_{n-k}^{\left(\alpha \right)})e_i^k\right]+a_2{\nabla }_t{e}_i^n\hfill \\ & =a_3{\delta }_x^2{e}_i^n-a_4{e}_i^n-\frac{a_5{\tau }^{-\beta }}{\Gamma (2-\beta )}\left[b_0^{\left(\beta \right)}{e}_i^n-\sum _{k=1}^{n-1}(b_{n-k-1}^{\left(\beta \right)}-b_{n-k})e_i^k-b_{n-1}^{\left(\beta \right)}{e}_i^0\right]+R_i^n,\hfill \end{array}$$

(68)

with the initial conditions $e_i^0=0,i=1,2,\cdots ,M-1$ and $e_0^n=e_M^n=0$.

It follows from the proof of Theorem 2 that

$$\begin{array}{cc}\hfill \left|\right|e^n{\left|\right|}_0^2& \le \frac{\Gamma (2-\alpha )\tau h}{a_6(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\sum _{i=1}^{M-1}{\left(R_i^n\right)}^2\hfill \\ \hfill & \le \frac{\Gamma (2-\alpha )\tau h}{a_6(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}\sum _{n=1}^N\sum _{i=1}^{M-1}{C}^2{(\tau +h^2)}^2\hfill \\ \hfill & \le \frac{\Gamma (2-\alpha )C^{2}n\tau (M-1)h}{a_6(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}{(\tau +h^2)}^2\hfill \\ \hfill & \le \frac{\Gamma (2-\alpha )C^{2}TL}{a_6(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}{(\tau +h^2)}^2,\hfill \end{array}$$

(69)

which implies

$$\left|\right|e^n{\left|\right|}_0\le C_1(\tau +h^2),$$

(70)

where $C_1=C\sqrt{\frac{\Gamma (2-\alpha )TL}{a_6(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}}$. In the same way, we have

$$\left|\right|e^n{\left|\right|}_1\le C_2(\tau +h^2),$$

(71)

and

$$\left|\right|e^n{\left|\right|}_{\infty }\le C_3(\tau +h^2),$$

(72)

where $C_2=C\sqrt{\frac{\Gamma (2-\alpha )TL}{a_3(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}}$ and $C_3=\frac{CL}{2}\sqrt{\frac{\Gamma (2-\alpha )TL}{a_3(a_1{T}^{1-\alpha }+2a_2\Gamma (2-\alpha ))}}$. □

5. Numerical Simulation and Discussion

Example 1.

The non-Newtonian Maxwell fluid model with a multi-term fractional order is considered as follows:

$$\left\{\begin{array}{cc}\hfill & a_1{D}_t^{\alpha }u(x,t)+a_2\frac{\partial u(x,t)}{\partial t}=a_3\frac{{\partial }^{2}u(u,t)}{\partial x^2}-a_{4}u(x,t)-a_5{D}_t^{\beta }u(x,t)+f(x,t),\hfill \\ \hfill & u(x,0)=sin\pi x,u_t(x,0)=0,\hfill \\ \hfill & u(0,t)=0,u(1,t)=0,\hfill \end{array}\right.$$

(73)

where $(x,t)\in (0,1)\times (0,1]$, $0\le x\le 1$, $0\le t\le 1$, $1<\alpha <2,0<\beta <1$, and

$$f(x,t)=sin\pi x\left[a_1\frac{\Gamma \left(3\right)}{\Gamma (3-\alpha )}{t}^{2-\alpha }+2a_{2}t+a_3{\pi }^2(t^2+1)+a_4(t^2+1)+a_5\frac{\Gamma \left(3\right)}{\Gamma (3-\beta )}{t}^{2-\beta }\right].$$

(74)

The exact solution to the above fractional order formula (Equation (73)) is $u(x,t)=(t^2+1)sin\pi x$. In the following, we choose the parameters to be $a_1=a_2=a_3=a_4=a_5=1$.

The implicit finite difference scheme in Equation (30) is used to solve the fractional order formula (Equation (73)). We can obtain the numerical results, which are in

Table 1. The temporal error $L_2$, $L_{\infty }$ and convergence of the difference scheme I with different orders $\alpha$ and $\beta$ when $h=0.001$ and $t=1$.

$\alpha =1.5$, $\beta =0.5$	${\left|\right|E(h,\tau )\left|\right|}_0$	Order	${\left|\right|E(h,\tau )\left|\right|}_{\infty }$	Order
1/40	2.8621 × ${10}^{-3}$		4.0476 × ${10}^{-3}$
1/80	1.4087 × ${10}^{-3}$	1.023	1.9922 × ${10}^{-3}$	1.023
1/160	6.9717 × ${10}^{-4}$	1.015	9.8594 × ${10}^{-4}$	1.015
1/320	3.4645 × ${10}^{-4}$	1.009	4.8995 × ${10}^{-4}$	1.009
1/640	1.7278 × ${10}^{-4}$	1.004	2.4434 × ${10}^{-4}$	1.004
$\alpha =1.5$, $\beta =0.7$	${\left|\right|E(h,\tau )\left|\right|}_0$	Order	${\left|\right|E(h,\tau )\left|\right|}_{\infty }$	Order
1/40	3.1908 × ${10}^{-3}$		4.5125 × ${10}^{-3}$
1/80	1.5587 × ${10}^{-3}$	1.034	2.2044 × ${10}^{-3}$	1.034
1/160	7.6515 × ${10}^{-4}$	1.027	1.0821 × ${10}^{-3}$	1.027
1/320	3.7717 × ${10}^{-4}$	1.021	5.3340 × ${10}^{-4}$	1.021
1/640	1.8668 × ${10}^{-4}$	1.015	2.6400 × ${10}^{-4}$	1.015
$\alpha =1.8$, $\beta =0.7$	${\left|\right|E(h,\tau )\left|\right|}_0$	Order	${\left|\right|E(h,\tau )\left|\right|}_{\infty }$	Order
1/40	2.7492 × ${10}^{-3}$		3.8880 × ${10}^{-3}$
1/80	1.3171 × ${10}^{-3}$	1.061	1.8626 × ${10}^{-3}$	1.062
1/160	6.3821 × ${10}^{-4}$	1.045	9.0257 × ${10}^{-3}$	1.045
1/320	3.1190 × ${10}^{-4}$	1.033	4.4110 × ${10}^{-4}$	1.033
1/640	1.5350 × ${10}^{-4}$	1.023	2.1708 × ${10}^{-4}$	1.023

, and it yields the $L_2$ and $L_{\infty }$ error and a convergence order of τ with different $\alpha ,\beta$, and $h=1/1000$ at $t=1$. Meanwhile, we found that the results of the numerical simulation were in perfect agreement with the exact solution, which demonstrates that the convergence order reached the expected first order. The comparison of the exact solution and the numerical solution of the fractional order formula (Equation (73)) is given in Figure 1, which further demonstrates the effectiveness and accuracy of our methods. In addition, we used MATALAB R2019b for the numerical computations, and the computation time was 5.0183 s for $h=1/1000$ and $\tau =1/640$. The configuration of the Lenovo desktop was as follows: an Inter Core i7-6700HQ, 2.60 GHz with 16 GB of RAM.

Example 2.

The following non-Newtonian Maxwell fluid model with a multi-term fractional order is considered:

$$\left\{\begin{array}{cc}\hfill & a_1{D}_t^{\alpha }u(x,t)+a_2\frac{\partial u(x,t)}{\partial t}=a_3\frac{{\partial }^{2}u(u,t)}{\partial x^2}-a_{4}u(x,t)-a_5{D}_t^{\beta }u(x,t)+f(x,t),\hfill \\ \hfill & u(x,0)=2sin\pi x,u_t(x,0)=0,\hfill \\ \hfill & u(0,t)=0,u(1,t)=0,\hfill \end{array}\right.$$

(75)

where $(x,t)\in (0,1)\times (0,1]$, $0\le x\le 1$, $0\le t\le 1$, $1<\alpha <2,0<\beta <1$, and

$$f(x,t)=sin\pi x\left[a_1\frac{\Gamma \left(4\right)}{\Gamma (4-\alpha )}{t}^{3-\alpha }+3a_2{t}^2+a_3{\pi }^2(t^3+2)+a_4(t^3+2)+a_5\frac{\Gamma \left(4\right)}{\Gamma (4-\beta )}{t}^{3-\beta }\right].$$

(76)

The exact solution to the above fractional-order formula (Equation (75)) is $u(x,t)=(t^3+2)sin\pi x$. In the following, we choose the parameters to be $a_1=1,a_2=2,a_3=3,a_4=4$, and $a_5=5$.

The implicit finite difference scheme in Equation (30) is used to solve the fractional order formula in Equation (75). We can obtain the numerical results, which are in

Table 2. The temporal error $L_2$, $L_{\infty }$ and convergence of difference scheme I with different orders $\alpha$ and $\beta$ when $h=0.001$ and $t=1$.

$\alpha =1.5$, $\beta =0.5$	${\left|\right|E(h,\tau )\left|\right|}_0$	Order	${\left|\right|E(h,\tau )\left|\right|}_{\infty }$	Order
1/40	4.4875 × ${10}^{-3}$		6.3463 × ${10}^{-3}$
1/80	2.17887 × ${10}^{-3}$	1.042	3.0813 × ${10}^{-3}$	1.042
1/160	1.0650 × ${10}^{-3}$	1.033	1.5061 × ${10}^{-3}$	1.033
1/320	5.2380 × ${10}^{-4}$	1.024	7.4077 × ${10}^{-4}$	1.024
1/640	2.5912 × ${10}^{-4}$	1.015	3.6646 × ${10}^{-4}$	1.015
$\alpha =1.5$, $\beta =0.7$	${\left|\right|E(h,\tau )\left|\right|}_0$	Order	${\left|\right|E(h,\tau )\left|\right|}_{\infty }$	Order
1/40	5.3384 × ${10}^{-3}$		7.5497 × ${10}^{-3}$
1/80	2.5484 × ${10}^{-3}$	1.067	3.6040 × ${10}^{-3}$	1.067
1/160	1.2232 × ${10}^{-3}$	1.059	1.7300 × ${10}^{-3}$	1.059
1/320	5.9084 × ${10}^{-4}$	1.050	8.3558 × ${10}^{-4}$	1.050
1/640	2.8727 × ${10}^{-4}$	1.040	4.0627 × ${10}^{-4}$	1.040
$\alpha =1.8$, $\beta =0.7$	${\left|\right|E(h,\tau )\left|\right|}_0$	Order	${\left|\right|E(h,\tau )\left|\right|}_{\infty }$	Order
1/40	5.8816 × ${10}^{-3}$		8.3178 × ${10}^{-3}$
1/80	2.7960 × ${10}^{-3}$	1.073	3.9541 × ${10}^{-3}$	1.073
1/160	1.3360 × ${10}^{-3}$	1.066	1.8893 × ${10}^{-3}$	1.066
1/320	6.4204 × ${10}^{-4}$	1.057	9.0799 × ${10}^{-4}$	1.057
1/640	3.1052 × ${10}^{-4}$	1.048	4.3915 × ${10}^{-4}$	1.048

, and giving the $L_2$ and $L_{\infty }$ error and the convergence order of τ with different $\alpha ,\beta$, and $h=1/1000$ at $t=1$. Meanwhile, we found that the results of the numerical simulation were in perfect agreement with the exact solution, which demonstrates that the convergence order reached the expected first order. The comparison of the exact solution and the numerical solution of the fractional order formula in Equation (75) is given in Figure 2, which further demonstrates the effectiveness and accuracy of our methods. In addition, we used MATALAB R2019b to yield the numerical computations, and the computation time was 5.0879s for $h=1/1000$ and $\tau =1/640$. The configuration of the Lenovo desktop was as follows: an Inter Core i7-6700HQ, 2.60 GHz with 16 GB of RAM.

6. Conclusions

In our paper, we studied the application of the finite difference method to new unsteady, generalized Maxwell fluids with multi-term time-fractional derivatives. We proposed a new finite difference scheme to approximate the multi-term fractional-order Maxwell model. The stability and convergence analysis of this finite difference scheme was proposed based on the discrete $H^1$. Then, we proved that the accuracy of the new finite difference method was $O(\tau +h^2)$. Finally, a numerical example was used to illustrate the validity and rationality of the method. Our method is effective enough to be generalized to solve the fractional-order generalized Oldroyd-B fluid model, fractional-order generalized Burgers fluid model, and other non-Newtonian second-order fluid models. We will use our techniques and methods to simulate some new multi-term, fractional-order non-Newtonian fluid models, such as the Oldroyd-B fluid model and Burgers fluid model with high dimensions, in a future work.