On a Multigrid Method for Tempered Fractional Diffusion Equations

Abstract

In this paper, we develop a suitable multigrid iterative solution method for the numerical solution of second- and third-order discrete schemes for the tempered fractional diffusion equation. Our discretizations will be based on tempered weighted and shifted Grünwald difference (tempered-WSGD) operators in space and the Crank–Nicolson scheme in time. We will prove, and show numerically, that a classical multigrid method, based on direct coarse grid discretization and weighted Jacobi relaxation, performs highly satisfactory for this type of equation. We also employ the multigrid method to solve the second- and third-order discrete schemes for the tempered fractional Black–Scholes equation. Some numerical experiments are carried out to confirm accuracy and effectiveness of the proposed method.

1. Introduction

In this paper, we will develop a multigrid method to numerically solve, highly efficiently, the tempered fractional diffusion equation. Multigrid is known to be an efficient and powerful numerical technique, particularly, for solving elliptic partial differential equations (PDEs). The convergence rate is usually independent of the mesh size [1,2]. Different methodological enhancements have been considered to generalize the range of applicability of the multigrid method, such as nontrivial smoothing methods, for example, in [3,4], coarsening and interpolation, like in [5,6], convection dominated problems [7], and so on, each time enlarging the range of robust and efficient multigrid applications, see also [8].

Recently, multigrid methods have also been applied to solving fractional diffusion equations (FDE) in the literature [9,10,11,12,13]. Fractional diffusion equations are governed by their long range interactions, so that, after discretization, full matrices result. These full matrices may possess a favorable structure, like a Toeplitz matrix structure, which is beneficial regarding efficient matrix-vector multiplication. Pang and Sun [9], for example, developed multigrid methods where the coarse grid operator retained the Toeplitz-like structure, by means of the Meerschaet–Tadjeran method. Hamid et al. constructed multigrid methods for a two-dimensional FDE problem, which was discretized by means of a CN-WSGD scheme, and they confirmed that multigrid methods performed better than classical preconditioners based on multilevel circulant matrices, in [13]. Gu, et al. [14] reformulated the classical time-stepping schemes as a kind of parallel-in-time (PinT) methods for both one- and two-dimensional space fractional diffusion equations and the fast Krylov subspace method with tau preconditioners is used to solve the resulting discretized linear systems.

It is well-known that a fractional derivative can be employed to accurately describe memory properties and hereditary effects of materials and processes. Differential equations with fractional operators are nowadays commonly applied in different fields of science and engineering, for example, in physics [15,16,17,18], hydrology [19,20,21,22], biology [23,24], or even finance [25,26,27,28]. Fractional derivatives also have a natural application when studying anomalous diffusion (for an extensive review, we refer to [29]). In another setting, Lévy flight models are used to mathematically describe the super-diffusion phenomenon, whose jumps have infinite moments in complex systems. So-called tempered fractional operators were introduced to describe probability density functions related to the positions of particles, by applying an exponential tempering of the probability of large jumps occurring in these Lévy flights [30]. These tempered derivatives have applications in physics [31,32,33], ground water hydrology [34], and even in finance [35,36].

A recent significant research effort on discretization methods for differential equations with tempered fractional derivatives has resulted in accurate finite element techniques [37], finite differences [31,33,38], and also spectral methods [32,39,40]. For example, Cartea and Del-Castillo-Negrete [35] defined a finite difference scheme to price exotic options under Lévy processes. Zhang et al. [36] presented a second-order discretization for the tempered fractional Black–Scholes equation and analyzed the stability and convergence properties of it. Li and Deng [31] defined higher-order discretizations based on a weighted and shifted Grünwald type approximation for the tempered fractional derivative. They also provided stability and convergence results for a second-order discretization of the tempered fractional diffusion equation. Zhao et al. [41] designed the first-order fully implicit and semi-implicit schemes for the nonlinear tempered fractional diffusion equation with variable coefficients, where the stabilities and convergences of the two numerical schemes are proved under several assumptions. Then the PinT implementation of the fully implicit scheme is given and the resulting nonlinear system is solved by using the fast preconditioned iterative method.

In [42], we developed the third-order discretizations based on the weighted and shifted Grünwald type difference (WSGD) for the tempered fractional derivatives. We also analyzed the stability and convergence properties for the tempered fractional diffusion equation, and proved that the third-order accurate scheme is unconditionally stable for a large ranges of problem parameters. A third-order scheme for the tempered Black–Scholes equation is also proposed and tested numerically. In this paper, we focus on the multigrid solution method for the tempered fractional diffusion and the fractional Black–Scholes equation, discretized by means of the second and third order CN-WSGD schemes we proposed before. Numerical results confirm that the proposed method is accurate and efficient.

The paper is organized as follows. In Section 2, we will provide the discretization details for the tempered fractional diffusion equation. Section 3.1 and Section 3.2 describe the components of the multigrid method for the second-order and the third-order discrete schemes for the fractional diffusion equation. A contribution of this paper is the multigrid convergence analysis for these discrete schemes in this section. In Section 4, we then present some numerical results to confirm the accuracy and efficiency of the proposed methods. Moreover, we also solve the fractional Black–Scholes equation in this section. Finally, we summarize our findings in the last section.

2. Numerical Schemes for the Tempered Fractional Diffusion Equation

We consider the following tempered fractional diffusion equation

$$\left\{\begin{array}{l}\frac{\partial u(x,t)}{\partial t}=c_l(x,t)·{}_a{\mathbb{D}}_x^{\alpha ,{\lambda }_1}u(x,t)+c_r(x,t)·{}_x{\mathbb{D}}_b^{\alpha ,{\lambda }_2}u(x,t)+f(x,t),\\ (x,t)\in (a,b)\times (0,T),\\ u(a,t)=0,u(b,t)=0,t\in (0,T),\\ u(x,0)=S\left(x\right),x\in (a,b),\end{array}\right.$$

(1)

where $\alpha \in (1,2)$, $f(x,t)$ is the source term, $c_l(x,t),c_r(x,t)\ge 0$ with $c_l(x,t)+c_r(x,t)\ne 0$, ${}_a{\mathbb{D}}_x^{\alpha ,\lambda }u\left(x\right){=}_a{D}_x^{\alpha ,\lambda }u\left(x\right)-\alpha {\lambda }^{\alpha -1}{\partial }_{x}u\left(x\right)-{\lambda }^{\alpha }u\left(x\right),$ and ${}_x{\mathbb{D}}_b^{\alpha ,\lambda }u\left(x\right){=}_x{D}_b^{\alpha ,\lambda }u\left(x\right)+$$\alpha {\lambda }^{\alpha -1}{\partial }_{x}u\left(x\right)-{\lambda }^{\alpha }u\left(x\right).$

The Riemann–Liouville tempered fractional derivatives, that we encounter in this equation, are defined as follows.

Definition 1

(See [31]).For $\alpha \in (n-1,n)$, let $u\left(x\right)$ be $(n-1)$-times continuously differentiable on $(a,b)$ with its nth derivative integrable on any subinterval of $[a,b]$, and $\lambda \ge 0.$ Then, the left Riemann–Liouville tempered fractional derivative of order α is defined as

$${}_a{D}_x^{\alpha ,\lambda }u\left(x\right)=\left(e^{-\lambda x}{}_a{D}_x^{\alpha }{e}^{\lambda x}\right)u\left(x\right)=\frac{e^{-\lambda x}}{\mathsf{\Gamma }(n-\alpha )}\frac{d^n}{d{x}^n}{\int }_a^x\frac{e^{\lambda \xi }u\left(\xi \right)}{{(x-\xi )}^{\alpha -n+1}}d\xi ;$$

(2)

the right Riemann–Liouville tempered fractional derivative of order α is defined as

$${}_x{D}_b^{\alpha ,\lambda }u\left(x\right)=\left(e^{\lambda x}{}_x{D}_b^{\alpha }{e}^{-\lambda x}\right)u\left(x\right)=\frac{{(-1)}^n{e}^{\lambda x}}{\mathsf{\Gamma }(n-\alpha )}\frac{d^n}{d{x}^n}{\int }_x^b\frac{e^{-\lambda \xi }u\left(\xi \right)}{{(\xi -x)}^{\alpha -n+1}}d\xi ,$$

(3)

where ’a’ and ’b’ can be extended to $-\infty$ and ∞, respectively.

We will construct a high-order scheme based on the tempered-WSGD operators for the tempered fractional derivative in space. The following results are developed for the tempered fractional operators in [31,42].

Remark 1.

In this paper, we consider a well-defined function $u\left(x\right)$ on a bounded interval $[a,b]$, and the function $u\left(x\right)$ will be zero extended for $x<a$ or $x>b$, so that $u\left(x\right)\in L^1\left(\mathbb{R}\right)$, and ${}_a{D}_x^{\alpha +l,\lambda }u\left(x\right)$, ${}_x{D}_b^{\alpha +l,\lambda }u\left(x\right)$ and their Fourier transforms belong to $L^1\left(\mathbb{R}\right)$. The α-th order left and right Riemann–Liouville tempered fractional derivatives of $u\left(x\right)$ at grid point x can then be approximated by tempered-WSGD operators ${}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }u$ and ${}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }u$, as follows

$$\begin{array}{ll}{}_a{D}_x^{\alpha ,\lambda }u\left(x\right)-{\lambda }^{\alpha }u\left(x\right)& =\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\left[\frac{x-a}{h}\right]+p}}{g}_{l,\lambda }^{(k,\alpha )}u(x-(l-p)h)-\frac{1}{h^{\alpha }}{\displaystyle \sum _{j=1}^m}{\gamma }_j{e}^{p_{j}h\lambda }{(1-e^{-h\lambda })}^{\alpha }u\left(x\right)+O\left(h^k\right)\\ & ={}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }u\left(x\right)+O\left(h^k\right),\end{array}$$

$$\begin{array}{ll}{}_x{D}_b^{\alpha ,\lambda }u\left(x\right)-{\lambda }^{\alpha }u\left(x\right)& =\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\left[\frac{b-x}{h}\right]+p}}{g}_{l,\lambda }^{(k,\alpha )}u(x+(l-p)h)-\frac{1}{h^{\alpha }}{\displaystyle \sum _{j=1}^m}{\gamma }_j{e}^{p_{j}h\lambda }{(1-e^{-h\lambda })}^{\alpha }u\left(x\right)+O\left(h^k\right)\\ & ={}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }u\left(x\right)+O\left(h^k\right),\end{array}$$

see [31,42] for details. The second- and third-order operators are given in Section 2.1 and Section 2.2.

Let the equidistant time partition, $t_j=j\tau$$(0\le t_j\le T,j=0,\cdots ,N)$, and spatial grid, $x_i=a+ih$$(a\le x_i\le b,i=0,\cdots ,M)$, be defined, where $\tau =T/N$ and $h=(b-a)/M$. Using the high-order tempered-WSGD operators ${}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }u$ and ${}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }u$ (as explained in Remark 1), high-order scheme for the first-order spatial derivative with ${\delta }_{k,x}u={\partial }_{x}u+O\left(h^k\right)$, and a Crank–Nicolson discretization in time, the numerical scheme for (1) reads

$$\begin{array}{l}\frac{u_i^{j+1}-u_i^j}{\tau }=c_{l,i}^{j+\frac{1}{2}}·\left({}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{u}_i^{j+\frac{1}{2}}-\alpha {\lambda }^{\alpha -1}{\delta }_{k,x}{u}_i^{j+\frac{1}{2}}\right)\\ +c_{r,i}^{j+\frac{1}{2}}·\left({}_R{\mathcal{D}}_{h,l}^{\alpha ,\lambda }{u}_i^{j+\frac{1}{2}}+\alpha {\lambda }^{\alpha -1}{\delta }_{k,x}{u}_i^{j+\frac{1}{2}}\right)+f_i^{j+\frac{1}{2}}+O({\tau }^2+h^k),\end{array}$$

(4)

where $u_i^j$ represents the solution of (1) at the point $(x_i,t_j),$$c_{l,i}^j=c_l(x_i,t_j),c_{r,i}^j=c_r(x_i,t_j)$ and $f_i^{j+\frac{1}{2}}=\frac{1}{2}(f(x_i,t_j)+f(x_i,t_{j+1})).$ Rewriting gives us

$$\begin{array}{l}{u}_i^{j+1}-\frac{\tau }{2}\left[c_{l,i}^{j+1}\left({}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{u}_i^{j+1}\right)+c_{r,i}^{j+1}\left({}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{u}_i^{j+1}\right)-\alpha {\lambda }^{\alpha -1}\left(c_{l,i}^{j+1}-c_{r,i}^{j+1}\right){\delta }_{k,x}{u}_i^{j+1}\right]\\ =u_i^j+\frac{\tau }{2}\left[c_{l,i}^j\left({}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{u}_i^j\right)+c_{r,i}^j\left({}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{u}_i^j\right)-\alpha {\lambda }^{\alpha -1}\left(c_{l,i}^{j+1}-c_{r,i}^{j+1}\right){\delta }_{k,x}{u}_i^j\right]\\ +\tau f_i^{j+\frac{1}{2}}+O({\tau }^3+\tau h^k).\end{array}$$

(5)

Denote by $U_i^j$ the solution of the numerical scheme for (1) at point $(x_i,t_j)$. The numerical scheme can now be written as

$$\begin{array}{l}{U}_i^{j+1}-\frac{\tau }{2}\left[c_{l,i}^{j+1}\left({}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{U}_i^{j+1}\right)+c_{r,i}^{j+1}\left({}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{U}_i^{j+1}\right)-\alpha {\lambda }^{\alpha -1}\left(c_{l,i}^{j+1}-c_{r,i}^{j+1}\right){\delta }_{k,x}{U}_i^{j+1}\right]\\ =U_i^j+\frac{\tau }{2}\left[c_{l,i}^j\left({}_L{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{U}_i^j\right)+c_{r,i}^j\left({}_R{\mathcal{D}}_{h,k}^{\alpha ,\lambda }{U}_i^j\right)-\alpha {\lambda }^{\alpha -1}\left(c_{l,i}^{j+1}-c_{r,i}^{j+1}\right){\delta }_{k,x}{U}_i^j\right]+\tau f_i^{j+\frac{1}{2}}.\end{array}$$

(6)

We will use the following notations, for vector $U^n={(U_1^n,U_2^n,\dots ,U_{M-1}^n)}^T$. Further, ${\varphi }^{l,m}\left(\lambda \right)={\sum }_{j=1}^m{\gamma }_j{e}^{p_j\lambda }{(1-e^{-h\lambda })}^{\alpha }$, $C_l^j=diag(c_{l,1}^j,c_{l,2}^j,\cdots ,c_{l,M-1}^j)$, $C_l^j=diag(c_{r,1}^j,c_{r,2}^j,$$\cdots ,c_{r,M-1}^j)$, and

$$B_{k,\lambda }=\left(\begin{array}{ccccc}{g}_{1,\lambda }^{(k,\alpha )}-{\varphi }^{k,m}\left(\lambda \right)& g_{0,\lambda }^{(k,\alpha )}& & & \\ g_{2,\lambda }^{(k,\alpha )}& g_{1,\lambda }^{(k,\alpha )}-{\varphi }^{k,m}\left(\lambda \right)& g_{0,\lambda }^{(k,\alpha )}& & \\ ⋮& g_{2,\lambda }^{(k,\alpha )}& g_{1,\lambda }^{(k,\alpha )}-{\varphi }^{k,m}\left(\lambda \right)& \ddots & \\ g_{n-2,\lambda }^{(k,\alpha )}& \cdots & \ddots & \ddots & g_{0,\lambda }^{(k,\alpha )}\\ g_{n-1,\lambda }^{(k,\alpha )}& g_{n-2,\lambda }^{(k,\alpha )}& \cdots & g_{2,\lambda }^{(k,\alpha )}& g_{1,\lambda }^{(k,\alpha )}-{\varphi }^{k,m}\left(\lambda \right)\end{array}\right).$$

(7)

The corresponding matrix form of (6) then reads

$${\widehat{A}}_{h,k}^{j+1}{U}^{j+1}+\frac{\alpha {\lambda }^{\alpha -1}\tau }{2}(C_l^{j+1}+C_r^{j+1}){\delta }_{k,x}{U}^{j+1}=f_{h,k}^{j+1},$$

(8)

where

$${\widehat{A}}_{h,k}^{j+1}=\left(I-\frac{\tau }{2h^{\alpha }}(C_l^{j+1}{B}_{k,\lambda }-C_r^{j+1}{B}_{k,\lambda }^T)\right),$$

(9)

$$f_{h,k}^{j+1}=\left(I+\frac{\tau }{2h^{\alpha }}(C_l^j{B}_{k,\lambda }+C_r^j{B}_{k,\lambda }^T)\right)U^j-\frac{\alpha {\lambda }^{\alpha -1}\tau }{2}(C_l^j+C_r^j){\delta }_{k,x}{U}^j+\tau {\widehat{F}}_i^{j+\frac{1}{2}},$$

(10)

and

$${\widehat{F}}^{n+\frac{1}{2}}=\left(\begin{array}{c}{f}_1^{n+\frac{1}{2}}\\ f_2^{n+\frac{1}{2}}\\ ⋮\\ f_{M-2}^{n+\frac{1}{2}}\\ f_{M-1}^{n+\frac{1}{2}}\end{array}\right)+\frac{U_0^{n+\frac{1}{2}}}{h^{\alpha }}\left(\begin{array}{c}{c}_{l,1}^{n+\frac{1}{2}}{g}_{2,\lambda }^{(k,\alpha )}+c_{r,1}^{n+\frac{1}{2}}{g}_{0,\lambda }^{(k,\alpha )}\\ c_{l,2}^{n+\frac{1}{2}}{g}_{3,\lambda }^{(k,\alpha )}\\ ⋮\\ c_{l,M-2}^{n+\frac{1}{2}}{g}_{M-1,\lambda }^{(k,\alpha )}\\ c_{l,M-1}^{n+\frac{1}{2}}{g}_{M,\lambda }^{(k,\alpha )}\end{array}\right)+\frac{U_M^{n+\frac{1}{2}}}{h^{\alpha }}\left(\begin{array}{c}{c}_{r,1}^{n+\frac{1}{2}}{g}_{M,\lambda }^{(k,\alpha )}\\ c_{r,2}^{n+\frac{1}{2}}{g}_{M-1,\lambda }^{(k,\alpha )}\\ ⋮\\ c_{r,M-2}^{n+\frac{1}{2}}{g}_{3,\lambda }^{(k,\alpha )}\\ c_{r,M-1}^{n+\frac{1}{2}}{g}_{2,\lambda }^{(k,\alpha )}+c_{l,1}^{n+\frac{1}{2}}{g}_{0,\lambda }^{(k,\alpha )}\end{array}\right).$$

(11)

2.1. Second-Order Discrete Scheme for the Tempered Fractional Diffusion Equation

We first present the second-order scheme for the tempered fractional diffusion Equation (1). Here, the second-order operators are defined as follows,

$$\begin{array}{cc}\hfill {}_L{\mathcal{D}}_{h,2}^{\alpha ,\lambda }u\left(x_j\right)=& \hfill \frac{1}{h^{\alpha }}\sum _{k=0}^{j+1}{g}_{k,\lambda }^{(2,\alpha )}u\left(x_{j-k+1}\right)-\frac{1}{h^{\alpha }}{\gamma }_2(h,\lambda ){(1-e^{-h\lambda })}^{\alpha }u\left(x_j\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {}_R{\mathcal{D}}_{h,2}^{\alpha ,\lambda }u\left(x_j\right)=& \hfill \frac{1}{h^{\alpha }}\sum _{k=0}^{N-j+1}{g}_{k,\lambda }^{(2,\alpha )}u\left(x_{j+k-1}\right)-\frac{1}{h^{\alpha }}{\gamma }_2(h,\lambda ){(1-e^{-h\lambda })}^{\alpha }u\left(x_j\right),\hfill \end{array}$$

(13)

where

$${\gamma }_2(h,\lambda )={\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda },$$

with ${\gamma }_j$ satisfying the following conditions,

$$\left\{\begin{array}{l}{\gamma }_1=\frac{\alpha }{2}+{\gamma }_3,\\ {\gamma }_2=\frac{2-\alpha }{2}-2{\gamma }_3,\end{array}\right.$$

(14)

and the weights $g_{k,\lambda }^{(2,\alpha )},k=0,\cdots ,j+1,$ are given by

$$\begin{array}{l}{g}_{0,\lambda }^{(2,\alpha )}={\gamma }_1{\omega }_0^{\left(\alpha \right)}{e}^{h\lambda },g_{1,\lambda }^{(2,\alpha )}={\gamma }_1{\omega }_1^{\left(\alpha \right)}+{\gamma }_2{\omega }_0^{\left(\alpha \right)},\\ g_{k,\lambda }^{(2,\alpha )}=\left({\gamma }_1{\omega }_k^{\left(\alpha \right)}+{\gamma }_2{\omega }_{k-1}^{\left(\alpha \right)}+{\gamma }_3{\omega }_{k-2}^{\left(\alpha \right)}\right)e^{-(k-1)h\lambda },k\ge 2.\end{array}$$

(15)

We will present the second-order scheme for the tempered fractional diffusion equation in detail here. Using the tempered-WSGD operators, ${}_L{\mathcal{D}}_{h,2}^{\alpha ,\lambda }{=}_L{\mathcal{D}}_{h,-1,0,1}^{\alpha ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_4}$ and ${}_R{\mathcal{D}}_{h,2}^{\alpha ,\lambda }{=}_R{\mathcal{D}}_{h,-1,0,1}^{\alpha ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_4}$ for the tempered fractional derivatives, and the second-order scheme for the first-order spatial derivative, the numerical scheme can now be written as,

$$\begin{array}{cc}\hfill & \hfill U_i^{j+1}-{\displaystyle \frac{\tau }{2}}\left[c_{l,i}^{j+1}·\left({}_L{\mathcal{D}}_{h,2}^{\alpha ,\lambda }{U}_i^{j+1}\right)+c_{r,i}^{j+1}·\left({}_R{\mathcal{D}}_{h,2}^{\alpha ,\lambda }{U}_i^{j+1}\right)\right]+{\displaystyle \frac{\tau \alpha {\lambda }^{\alpha -1}}{4h}}\left(c_{l,i}^{j+1}-c_{r,i}^{j+1}\right)\left(U_{i+1}^{j+1}-U_{i-1}^{j+1}\right)\hfill \\ \hfill & \hfill =U_i^j+{\displaystyle \frac{\tau }{2}}\left[c_{l,i}^j·\left({}_L{\mathcal{D}}_{h,2}^{\alpha ,\lambda }{U}_i^j\right)+c_{r,i}^j·\left({}_R{\mathcal{D}}_{h,2}^{\alpha ,\lambda }{U}_i^j\right)\right]-{\displaystyle \frac{\tau \alpha {\lambda }^{\alpha -1}}{4h}}\left(c_{l,i}^j-c_{r,i}^j\right)\left(U_{i+1}^j-U_{i-1}^j\right)+f_i^{j+\frac{1}{2}}.\hfill \end{array}$$

(16)

The corresponding matrix form of (16) then reads

$$A_{h,2}^{j+1}{U}^{j+1}=f_{h,2}^{j+1},$$

(17)

with ${\widehat{A}}_{h,2}^{j+1}$, ${\widehat{F}}^{n+\frac{1}{2}}$ as defined in (9), (11) when $k=2$, $H_2=tridiag\{-1,0,1\}$,

$$A_{h,2}^{j+1}={\widehat{A}}_{h,2}^{j+1}+\frac{\tau \alpha {\lambda }^{\alpha -1}}{4h}(C_l^{j+1}-C_r^{j+1})H_2,$$

(18)

and

$$f_{h,2}^{j+1}=\left(I+\frac{\tau }{2h^{\alpha }}(C_l^j{B}_{2,\lambda }+C_r^j{B}_{2,\lambda }^T)-\frac{\tau \alpha {\lambda }^{\alpha -1}}{4h}(C_l^{j+1}-C_r^{j+1})H_2\right)U^j+\tau {\widehat{F}}_i^{j+\frac{1}{2}}.$$

(19)

The stability and convergence of the second-order scheme for the tempered fractional diffusion Equation (1), when $c_l(x,t)$ and $c_r(x,t)$ are constants have already been presented in [31]. In a similar way, we can derive and prove the following theorem, based on the lemma below.

Lemma 1

(From [31]).For $1<\alpha <2$ and $\lambda \ge 0$, if

$$max\left\{\frac{(2-\alpha )({\alpha }^2+\alpha -8)}{2({\alpha }^2+3\alpha +2)},\frac{(1-\alpha )({\alpha }^2+2\alpha )}{2({\alpha }^2+3\alpha +4)}\right\}<{\gamma }_3<\frac{(2-\alpha )({\alpha }^2+2\alpha -3)}{2({\alpha }^2+3\alpha +2)},$$

then the weight coefficients ${\omega }_k^{\left(\alpha \right)}$ and $g_{k,\lambda }^{(2,\alpha )}$ satisfy the following properties,

1. 

${\omega }_0^{\left(\alpha \right)}=1$, ${\omega }_1^{\left(\alpha \right)}=-\alpha ,$$0\le \dots \le {\omega }_3^{\left(\alpha \right)}\le {\omega }_2^{\left(\alpha \right)}\le 1$, ${\displaystyle \sum _{k=0}^{\infty }{\omega }_k^{\left(\alpha \right)}=0}$,

2. 

${\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }=1+{\gamma }_1\left(e^{h\lambda }+e^{-h\lambda }-2\right)+\frac{{\displaystyle \alpha }}{{\displaystyle 2}}\left(1-e^{-h\lambda }\right)>1$,

3. 

$g_{1,\lambda }^{(2,\alpha )}\le 0$, $g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(2,\alpha )}\ge 0$, $g_{k,\lambda }^{(2,\alpha )}\ge 0(k\ge 3).$

Theorem 1.

For $1<\alpha <2$, and $\lambda \ge 0$, if $a_1<{\gamma }_3<a_2$, the numerical scheme (16) is stable for

$$a_1=max\{\frac{(2-\alpha )({\alpha }^2+\alpha -8)}{2({\alpha }^2+3\alpha +2)},\frac{(1-\alpha )({\alpha }^2+2\alpha )}{2({\alpha }^2+3\alpha +4)}\}$$

and

$$a_2=\frac{(2-\alpha )({\alpha }^2+2\alpha -3)}{2({\alpha }^2+3\alpha +2)}.$$

Denoting $e_i^j=u_i^j-U_i^j,i=1,2,\dots ,M-1$ and $E^j={(e_1^j,e_2^j,\dots ,e_{M-1}^j)}^T,j=1,2,\dots ,N$, moreover, it is found that

$$E_h^j\le c({\tau }^2+h^2),1\le j\le N-1.$$

(20)

2.2. Third-Order Discrete Scheme for the Tempered Fractional Diffusion Equation

In this work, we also consider the third-order operators. They are defined as

$$\begin{array}{c}\hfill {}_L{\mathcal{D}}_{h,3}^{\alpha ,\lambda }u\left(x_j\right)=\frac{1}{h^{\alpha }}\left(\sum _{k=0}^{\left[\frac{x-a}{h}\right]+1}{g}_{k,\lambda }^{(3,\alpha )}u\left(x_{j-k+1}\right)-{\gamma }_3(h,\lambda ){(1-e^{-h\lambda })}^{\alpha }u\left(x_j\right)\right),\end{array}$$

(21)

and

$$\begin{array}{c}\hfill {}_R{\mathcal{D}}_{h,3}^{\alpha ,\lambda }u\left(x_j\right)=\frac{1}{h^{\alpha }}\left(\sum _{k=0}^{\left[\frac{b-x}{h}\right]+1}{g}_{k,\lambda }^{(3,\alpha )}u\left(x_{j+k-1}\right)-{\gamma }_3(h,\lambda ){(1-e^{-h\lambda })}^{\alpha }u\left(x_j\right)\right),\end{array}$$

(22)

where

$${\gamma }_3(h,\lambda )={\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }+{\gamma }_4{e}^{-2h\lambda },$$

with ${\gamma }_j$ satisfying the following conditions

$$\left\{\begin{array}{l}{\gamma }_1=\frac{{\alpha }^2}{8}+\frac{5}{24}\alpha -{\gamma }_4,\\ {\gamma }_2=-\frac{{\alpha }^2}{4}+\frac{1}{12}\alpha +1+3{\gamma }_4,\\ {\gamma }_3=\frac{{\alpha }^2}{8}-\frac{7}{24}\alpha -3{\gamma }_4.\end{array}\right.$$

(23)

The weights, $g_{k,\lambda }^{(3,\alpha )},k=0,\cdots ,j+1,$ are found to be

$$\begin{array}{l}{g}_{0,\lambda }^{(3,\alpha )}={\gamma }_1{\omega }_0^{\left(\alpha \right)}{e}^{h\lambda },g_{1,\lambda }^{(3,\alpha )}={\gamma }_1{\omega }_1^{\left(\alpha \right)}+{\gamma }_2{\omega }_0^{\left(\alpha \right)},\\ g_{2,\lambda }^{(3,\alpha )}=\left({\gamma }_1{\omega }_2^{\left(\alpha \right)}+{\gamma }_2{\omega }_1^{\left(\alpha \right)}+{\gamma }_3{\omega }_0^{\left(\alpha \right)}\right)e^{-h\lambda },\\ g_{k,\lambda }^{(3,\alpha )}=\left({\gamma }_1{\omega }_k^{\left(\alpha \right)}+{\gamma }_2{\omega }_{k-1}^{\left(\alpha \right)}+{\gamma }_3{\omega }_{k-2}^{\left(\alpha \right)}+{\gamma }_4{\omega }_{k-3}^{\left(\alpha \right)}\right)e^{-(k-1)h\lambda },k\ge 3.\end{array}$$

(24)

Using the tempered-WSGD operators, ${}_L{\mathcal{D}}_{h,3}^{\alpha ,\lambda }{=}_L{\mathcal{D}}_{h,-1,0,1,2}^{\alpha ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_4}$ and ${}_R{\mathcal{D}}_{h,3}^{\alpha ,\lambda }{=}_R{\mathcal{D}}_{h,-1,0,1,2}^{\alpha ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_4}$, for the tempered fractional derivatives, and the fourth-order scheme for the first-order spatial derivative, we find the following numerical discretization for (1)

$$\begin{array}{l}{U}_i^{j+1}-{\displaystyle \frac{\tau }{2}}\left[c_{l,i}^{j+1}·\left({}_L{\mathcal{D}}_{h,3}^{\alpha ,\lambda }{U}_i^{j+1}\right)+c_{r,i}^{j+1}·\left({}_R{\mathcal{D}}_{h,3}^{\alpha ,\lambda }{U}_i^{j+1}\right)\right]\\ +\frac{\tau \alpha {\lambda }^{\alpha -1}}{24h}\left(c_{l,i}^{j+1}-c_{r,i}^{j+1}\right)\left(8(U_{i+1}^{j+1}-U_{i-1}^{j+1})-(U_{i+2}^{j+1}-U_{i-2}^{j+1})\right)\\ =U_i^j+{\displaystyle \frac{\tau }{2}}\left[c_{l,i}^j·\left({}_L{\mathcal{D}}_{h,3}^{\alpha ,\lambda }{U}_i^j\right)V_{i+1}^j)+c_{r,i}^j·\left({}_R{\mathcal{D}}_{h,3}^{\alpha ,\lambda }{U}_i^j\right)\right]\\ -\frac{\tau \alpha {\lambda }^{\alpha -1}}{24h}\left(c_{l,i}^j-c_{r,i}^j\right)\left(8(U_{i+1}^j-U_{i-1}^j)-(U_{i+2}^j-U_{i-2}^j)\right)+f_i^{j+\frac{1}{2}}.\end{array}$$

(25)

The matrix form now looks as follows

$$A_{h,3}^{j+1}{U}^{j+1}=f_{h,3}^{j+1},$$

(26)

with ${\widehat{A}}_{h,3}^{j+1},$${\widehat{F}}^{n+\frac{1}{2}}$ as defined in (9), (11) with $k=3$,

$$f_{h,3}^{j+1}=\left(I+\frac{\tau }{2h^{\alpha }}(C_l^j{B}_{3,\lambda }+C_r^j{B}_{3,\lambda }^T)-\frac{\tau \alpha {\lambda }^{\alpha -1}}{24h}(C_l^{j+1}-C_r^{j+1})H_3\right)U^j+\tau {\widehat{F}}_i^{j+\frac{1}{2}},$$

(27)

and

$$A_{h,3}^{j+1}={\widehat{A}}_{h,3}^{j+1}+\frac{\tau \alpha {\lambda }^{\alpha -1}}{24h}(C_l^{j+1}-C_r^{j+1})H_3,$$

(28)

with

$$H_3=\left(\begin{array}{cccccc}0& 8& -1& & & \\ -8& 0& 8& -1& & \\ 1& -8& 0& 8& -1& \\ ⋮& 1& -8& 0& \ddots & \\ 0& \cdots & \ddots & \ddots & 0& 8\\ 0& \cdots & \cdots & 1& -8& 0\end{array}\right).$$

(29)

We have already discussed the stability and convergence of the third-order scheme for the tempered fractional diffusion Equation (1) when $c_l(x,t)$ and $c_r(x,t)$ are constants in the paper [42]. Before we introduce the stability and convergence of the third-order scheme (25), we define the functions

$$f_B\left(x\right)=\sum _{k=0}^{N-1}{g}_{k,\lambda }^{(3,\alpha )}{e}^{i(k-1)x}-\varphi \left(\lambda \right),f_{B^T}\left(x\right)=\sum _{k=0}^{N-1}{g}_{k,\lambda }^{(3,\alpha )}{e}^{-i(k-1)x}-\varphi \left(\lambda \right),$$

(30)

and the generating function,

$$f(\alpha ,\lambda ;x)=\frac{f_B\left(x\right)+f_{B^T}\left(x\right)}{2}.$$

(31)

We obtain the stability for the numerical scheme (25), based on the theorem below.

Lemma 2

(From [42]).Let the matrices $B_{3,\lambda }$ and $B_{3,\lambda }^T$ be given via the numerical scheme (7). For $\lambda \ge 0,h>0$ and $\alpha \in [1,2]$, if we can find (analytically, or with the help of numerical techniques) values of ${\gamma }_j$ for which the generating functions $f(\alpha ,\lambda ;x)$ of $B_{\lambda }$ are negative, then the eigenvalues of the matrix $B_{3,\lambda }$ are negative too.

In a similar way as in [42], we have the following theorem.

Theorem 2.

If, for $1<\alpha <2$, the generating functions $f(\alpha ,\lambda ;x)$ given in (31), are negative, the numerical scheme (25) is stable.

Theorem 3.

Assuming function $u(x,t)$ to be the solution of Equation (1) on a bounded interval $(a,b)\times (0,T)$, which can be zero extended for $x<a$ or $x>b$, so that $u\in L^1(0,T;\mathbb{R})$, and ${}_a{D}_x^{\alpha +3,\lambda }u$ and its Fourier transform also belong to $L^1(0,T;\mathbb{R})$. Let’s denote by $e_i^j=u_i^j-U_i^j,i=1,2,\dots ,M-1,$ and $E^j={(e_1^j,e_2^j,\dots ,e_{M-1}^j)}^T,j=1,2,\dots ,N$. With solutions $u_i^j$ and $U_i^j$ of Equations (5) when $k=3$ and (25), respectively, we have, for $1<\alpha <2$, if $f(\alpha ,{\lambda }_i;x)<0$, $i=1,2$,

$$E_h^j\le c({\tau }^2+h^3),1\le j\le N-1.$$

(32)

It is our aim in this paper to solve the resulting discrete equations by means of a multigrid technique. The challenge here is, of course, the occurance of the nonlocality of these discretization schemes for the tempered fractional derivatives.

3. Multigrid Method for Tempered Fractional Diffusion Equation

In this subsection, we provide a multigrid method (see, for example in [8]) to solve the presented linear systems originating from the discretized fractional diffusion equations. Actually, the classical multigrid setting will be employed here, based on the direct coarse grid discretization. The corresponding two-grid algorithmic description is the following:

1.

Pre-smoothing:

• Compute ${\tilde{U}}_h^{j+1,m}$ by applying ${\nu }_1(\ge 0)$ steps of a smoothing procedure to $U_h^{j+1,m}$

  $${\tilde{U}}_h^{j+1,m}={\mathcal{S}}^{{\nu }_1}(U_h^{j+1,m},A_{h,k},f_{h,k}^{j+1})$$

2.

Coarse-grid correction:

• Define the residual: $r_h=f_{h,k}^{j+1}-A_{h,k}{\tilde{U}}_h^{m,j+1}$,
• Restrict the residual (fine-to-coarse grid transfer): $r_H=I_h^H{r}_h$,
• Solve $A_{H,k}{\widehat{v}}_H=r_H$,
• Interpolate the correction (coarse-to-fine grid transfer): ${\widehat{v}}_h=I_H^h{\widehat{v}}_H$,
• Compute a new approximation: ${\breve{U}}_h^{j+1,m}={\tilde{U}}_h^{j+1,m}+{\widehat{v}}_h$.

3.

Post-smoothing:

• Compute $U_h^{j+1,m+1}$ by applying ${\nu }_2(\ge 0)$ steps of a smoothing procedure to ${\breve{U}}_h^{j+1,m}$

  $$U_h^{j+1,m+1}={\mathcal{S}}^{{\nu }_2}({\breve{U}}_h^{m,j+1},A_{h,k},f_{h,k}^{j+1}).$$

For the above description, the notation is as follows:

• ${\nu }_1,{\nu }_2$ denote the number of smoothing steps. We will use ${\nu }_i=0,1,2$.
• The classical fine-to-coarse restriction operator $I_h^H$ is employed,

  $$I_h^H=\frac{1}{4}\left(\begin{array}{ccccccc}1& 2& 1& & & & \\ & 1& 2& 1& & & \\ & & & \ddots & & & \\ & & & & 1& 2& 1\end{array}\right).$$

  (33)
• The scaled transpose of the restriction is the coarse-to-fine interpolation operator, i.e., $I_H^h=2{\left(I_h^H\right)}^{\prime }$.
• The fine grid operator $A_{h,k}$ and the coarse grid operator $A_{H,k}$ are defined as in Equations (18) or (28). Obviously, the coefficient matrix possess the Toeplitz-like structure.
• The recursive generalization of this classical two-grid scheme towards multiple grids is well-known.

For the tempered fractional diffusion equation, we will be using the damped Jacobi iteration as the smoother. Here we use

$$diag\left(A_{h,k}\right)z^{m+1}=\left(diag\left(A_{h,k}\right)-A_{h,k}\right)u^m+f^m.$$

(34)

Then

$$\begin{array}{ll}{u}^{m+1}& =u^m+\omega \left(z^{m+1}-u^m\right)\\ & =u^m+\omega \left(diag{\left(A_{h,k}\right)}^{-1}\left(\left(diag\left(A_{h,k}\right)-A_{h,k}\right)u^m+f^m\right)-u^m\right)\\ & =\left(I-\omega ·diag{\left(A_{h,k}\right)}^{-1}·A_{h,k}\right)u^m+\omega ·diag{\left(A_{h,k}\right)}^{-1}·f^m.\end{array}$$

(35)

We can easily generalize the classical Jacobi iteration by introducing a relaxation parameter $\omega$, in the standard way, i.e.,

$$S_{h,\omega }=I-\omega ·diag{\left(A_{h,k}\right)}^{-1}·A_{h,k},$$

(36)

where $diag\left(A_{h,k}\right)$ represents the main diagonal of matrix $A_{h,k}$. The smoother then becomes,

$${\mathcal{S}}_{\omega }(u_0,A_{h,k},f_{h,k})=S_{h,\omega }{u}_0+\omega ·\mathrm{diag}{\left(A_{h,k}\right)}^{-1}·f_{h,k}^{j+1}.$$

(37)

Remark 2.

For general full matrices, a matrix-vector multiplication is an expensive task. However, in the present context with the fractional diffusion equations, we can benefit from the choices made within the discretization and regarding the multigrid components. The overall computational complexity of the multigrid method here is therefore O($MlogM$) at each time step, despite the fact that we’re dealing with a full matrix. In this paper, the resulting coefficient matrix which contains three Toeplitz matrices possesses a Toeplitz-like structure. It is however nontrivial to use fast Toeplitz solvers directly when the coefficients $c_l,c_r$ would depend on the spatial position, i.e., $c=c(x,t)$.

3.1. Multigrid Convergence Analysis for Second-Order Tempered Fractional Diffusion Scheme

Here, we will analyze the convergence of the multigrid method for the second-order discrete scheme. To simplify the analysis, we assume that $c_r=c_l=c>0$. Then we have $A_{h,k}={\widehat{A}}_{h,k}$. We denote by

$$\begin{array}{lll}{B}_h& =& -\frac{c\tau }{2h^{\alpha }}(B_{2,\lambda }+B_{2,\lambda }^T),\\ a_0& =& 1-2d\left(g_{1,\lambda }^{(2,\alpha )}-\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }\right),\\ a_1& =& a_{-1}=-d\left(g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(k,\alpha )}\right)\\ a_j& =& a_{-j}=-d{g}_{j+1,\lambda }^{(2,\alpha )},\end{array}$$

with $j=2,3,4,\dots$. Then, we find that $A_{h,2}=I+B_h$ is a symmetric Toeplitz matrix, of the following form,

$$A_{h,2}=I-d(B_{2,\lambda }+B_{2,\lambda }^T)=\left(\begin{array}{ccccc}{a}_0& a_1& a_2& \cdots & a_{N-2}\\ a_{-1}& a_0& a_1& \cdots & a_{N-3}\\ a_{-2}& a_{-1}& a_0& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & a_2\\ a_{-(N-3)}& \ddots & \ddots & \ddots & a_1\\ a_{-(N-2)}& a_{-(N-3)}& \cdots & a_{-1}& a_0\end{array}\right),$$

(38)

where $d=\frac{{\displaystyle c\tau }}{{\displaystyle 2h^{\alpha }}}$.

We need the following lemmas regarding the properties of our matrix, in order to prove multigrid convergence.

Lemma 3

(From [31]).Using the notation,

$$b_1=max\left\{\frac{(2-\alpha )({\alpha }^2+\alpha -8)}{2({\alpha }^2+3\alpha +2)},\frac{(1-\alpha )({\alpha }^2+2\alpha )}{2({\alpha }^2+3\alpha +4)}\right\}and{b}_2=\frac{(2-\alpha )({\alpha }^2+2\alpha -3)}{2({\alpha }^2+3\alpha +2)}.$$

For $1<\alpha <2$, and $\lambda \ge 0$, if $b_1<{\gamma }_3<b_2$, then the matrix,

$$B_h=\frac{B+B^T}{2},$$

is diagonally dominant and all eigenvalues of $B_h$ are negative.

Moreover, if a matrix is real-valued, symmetric, strictly diagonally dominant or irreducibly diagonally dominant, with positive diagonal entries, then it is positive [31].

We have the following lemma regarding our matrix $A_{h,2}$.

Lemma 4.

For $1<\alpha <2$ and $\lambda \ge 0$, if $b_1<{\gamma }_3<b_2$, the matrix $A_{h,2}$ defined in (18) is diagonally dominant and all eigenvalues of $A_{h,2}$ are positive.

Proof. 

From Lemma 1, we obtain,

$$\begin{array}{l}{a}_0>1,\\ a_j=a_{-j}<0,j=1,2,3,\cdots .\end{array}$$

(39)

Therefore, we have the following result for the ith row of the matrix $A_{h,2}$

$$\begin{array}{ll}{\displaystyle \sum _{j=-(i-1)}^{N-i-1}}{a}_j\ge & {\displaystyle \sum _{j=-\infty }^{\infty }}{a}_j\\ & =|a_0|-2{\displaystyle \sum _{j=1}^{\infty }}\left|a_j\right|\\ & =1-2d\left({\displaystyle \sum _{j=0}^{\infty }}{g}_{j,\lambda }^{(2,\alpha )}-\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }\right)\\ & =1.\end{array}$$

(40)

It can seen that the matrix $A_{h,2}$ is strictly diagonally dominant for $1<\alpha <2$. From Lemma (3), we then conclude that $A_{h,2}$ is a symmetric, positive definite matrix. □

We will use the following inner products,

$${\langle u,v\rangle }_0=\langle diag\left(A_{h,2}\right)u,v\rangle ,{\langle u,v\rangle }_1=\langle A_{h,2}u,v\rangle ,{\langle u,v\rangle }_2=\langle diag{\left(A_{h,2}\right)}^{-1}{A}_{h,2}u,A_{h,2}v\rangle .$$

Here $\langle ·,·\rangle$ is the Euclidean inner product.

Theorem 4

(From [43]).For a symmetric, positive definite matrix $A_{h,k}$, suppose that the damping parameter ω, in the damped Jacobi smoother, in (36), is properly chosen, to fulfill

$$1/\omega \ge \rho \left(diag{\left(A_{h,k}\right)}^{-1}{A}_{h,k}\right),$$

(41)

where $\rho (·)$ denotes the spectral radius of the matrix. Then, $S_{h,\omega }$ in (36) satisfies

$${\parallel S_{h,\omega }{e}_h\parallel }_1^2\le {\parallel e_h\parallel }_1^2-\omega {\parallel e_h\parallel }_2^2,\forall e_h\in {\mathbb{R}}^{N-1}.$$

(42)

The inequality (42) is the well-known smoothing property [2]. We find that $\parallel S_{h,\omega }\parallel \le 1$, when $\omega$ satisfies (41). As $A_{h,2}$ is symmetric, positive definite and diagonally dominant, we have

$$\rho \left(diag{\left(A_{h,2}\right)}^{-1}[A_{h,2}-diag\left(A_{h,2}\right)]\right)\le 1,$$

hence

$$\rho \left(diag{\left(A_{h,2}\right)}^{-1}{A}_{h,2}\right)\le 2.$$

(43)

Here we choose $0<\omega \le 1/2$ which satisfies (41).

For the two-grid method (TGM), the correction operator is given by

$$T_{TGM}=I-I_H^h{\left(A_{H,2}\right)}^{-1}{I}_h^H{A}_{h,2}.$$

Therefore, the convergence factor of the TGM reads $\parallel {\left(S_{h,\omega }\right)}^{{\nu }_2}{T}_{TGM}{\left(S_{h,\omega }\right)}^{{\nu }_1}\parallel$. For convenience, we consider here the case that ${\nu }_1=0,$ and ${\nu }_2=1$.

Theorem 5

(From [43]).Let $A_h$ be symmetric and positive definite and let $\omega >0$ be chosen such that $S_{h,\omega }$ satisfies the smoothing condition (42), i.e.,

$${\parallel S_{h,\omega }{e}_h\parallel }_1^2\le {\parallel e_h\parallel }_1^2-\omega {\parallel e_h\parallel }_2^2,\forall e_h\in {\mathbb{R}}^{N-1}.$$

(44)

Suppose that $I_H^h$ has full rank and that there exists a scalar $\beta >0$, such that

$$\underset{e_H\in {\mathbb{R}}^{N/2-1}}{\min }{\parallel e_h-I_H^h{e}_H\parallel }_0^2\le \beta {\parallel e_h\parallel }_1^2,\forall e_h\in {\mathbb{R}}^{N-1}.$$

(45)

Then, $\beta \ge \omega$ and the convergence factor of the TGM satisfies

$${\parallel S_{h,\omega }·T_{TGM}\parallel }_1\le \sqrt{1-\omega /\beta }.$$

(46)

In other words, we just need to find a suitable $\beta$-value to satisfy (46) and then we find that the convergence of TGM is independent of N. Let $L_{N-1}=tridiag(-1,2,-1)$ be the $(N-1)\times (N-1)$ one-dimensional discrete Laplacian matrix. Then $L_{N-1}$ is also a symmetric positive definite Toeplitz matrix.

Lemma 5.

For $1<\alpha <2$ and $\lambda \ge 0$, if $b_1<{\gamma }_3<b_2$, we denote $B_{rest}=B_h+a_1{L}_{N-1}$. Then $B_{rest}$ is symmetric positive definite.

Proof. 

Since both $B_h$ and $L_{N-1}$ are symmetric Toeplitz, $B_{rest}$ is also symmetric Toeplitz. We have

$$\begin{array}{ll}{B}_{rest}\hfill & =\left(\begin{array}{ccccc}{\tilde{b}}_0& {\tilde{b}}_1& {\tilde{b}}_2& \cdots & {\tilde{b}}_{N-2}\\ {\tilde{b}}_{-1}& {\tilde{b}}_0& {\tilde{b}}_1& \cdots & {\tilde{b}}_{N-3}\\ {\tilde{b}}_{-2}& {\tilde{b}}_1& {\tilde{b}}_0& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & {\tilde{b}}_2\\ {\tilde{b}}_{-(N-3)}& \ddots & \ddots & \ddots & {\tilde{b}}_1\\ {\tilde{b}}_{-(N-2)}& {\tilde{b}}_{-(N-3)}& \cdots & {\tilde{b}}_{-1}& {\tilde{b}}_0\end{array}\right)\hfill \\ & =\left(\begin{array}{ccccc}{a}_0-1+2a_1& 0& a_2& \cdots & a_{N-2}\\ 0& a_0-1+2a_1& 0& \cdots & a_{N-3}\\ a_{-2}& 0& a_0-1+2a_1& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & a_2\\ a_{-(N-3)}& \ddots & \ddots & \ddots & 0\\ a_{-(N-2)}& a_{-(N-3)}& \cdots & 0& a_0-1+2a_1\end{array}\right),\hfill \end{array}$$

(47)

where ${\tilde{b}}_0=-2d(g_0+g_2+g_1-\varphi \left(x\right))>0$, ${\tilde{b}}_1={\tilde{b}}_{-1}=0$ and

$${\tilde{b}}_j={\tilde{b}}_{-j}=-d{g}_{j+1}<0,j=2,3,\cdots ,N-1.$$

For the k-th row, we then obtain that

$$\begin{array}{ll}{\tilde{b}}_0-{\displaystyle \sum _{j=1-k,j\ne 0}^{N-(k+1)}}|\tilde{b}|\hfill & =-2d\left(g_0+g_2+g_1-\varphi \left(x\right)\right)-d({\displaystyle \sum _{j=3}^k}{g}_j+{\displaystyle \sum _{j=3}^{N-k}}{g}_j)\hfill \\ \hfill & \ge -2d\left(g_0+g_2+g_1-\varphi \left(x\right)\right)-d({\displaystyle \sum _{j=3}^{N-1}}{g}_j+{\displaystyle \sum _{j=3}^{N-1}}{g}_j)\hfill \\ \hfill & =-2d{\displaystyle \sum _{j=0}^{N-1}}{g}_j>0.\hfill \end{array}$$

(48)

Hence, $B_{rest}$ is strictly diagonally dominant, and we know that $B_{rest}$ is positive since ${\tilde{b}}_0>0$. □

From Lemma 5, it’s easy to see that

$$\langle A_{h,2}{u}_h,u_h\rangle =\langle (I+B_{rest}-a_1{L}_{N-1})u_h,u_h\rangle \ge \langle (I-a_1{L}_{N-1})u_h,u_h\rangle ,\forall u_h\in {\mathbb{R}}^{N-1}.$$

(49)

Now we are ready to provide the proof for the TGM convergence, in the following theorem.

Theorem 6.

Suppose that $A_{h,2}$ is defined in (18) and $\omega \le 1/2$ such that $S_{h,\omega }$ satisfies the smoothing condition (42). The convergence factor of the TGM satisfies

$${\parallel S_{h,\omega }·T_{TGM}\parallel }_1<1.$$

(50)

Proof. 

We denote

$$u_h=(u_1,u_2,\cdots ,u_{N-1})\in {\mathbb{R}}^{N-1}$$

and

$$u_H=({\tilde{u}}_1,{\tilde{u}}_2,\cdots ,{\tilde{u}}_{N/2-1})\in {\mathbb{R}}^{N/2-1},$$

where

${\tilde{u}}_j=u_{2j},$$1\le j\le N/2-1$.

Let $u_0=u_N=0$. Then we have,

$$\begin{array}{ll}{\parallel u_h-I_H^h{u}_H\parallel }_0^2\hfill & =a_0{\displaystyle \sum _{j=0}^{N/2-1}}{\left(u_{2j+1}-\frac{1}{2}{u}_{2j}-\frac{1}{2}{u}_{2j+2}\right)}^2\hfill \\ & =a_0{\displaystyle \sum _{j=0}^{N/2-1}}\left(u_{2j+1}^2+\frac{1}{4}{u}_{2j}^2+\frac{1}{4}{u}_{2j+2}^2-u_{2j+1}{u}_{2j}-u_{2j+1}{u}_{2j+2}+\frac{1}{2}{u}_{2j}{u}_{2j+2}\right)\\ & \le a_0{\displaystyle \sum _{j=0}^{N/2-1}}\left(u_{2j+1}^2+\frac{1}{2}{u}_{2j}^2+\frac{1}{2}{u}_{2j+2}^2-u_{2j+1}{u}_{2j}-u_{2j+1}{u}_{2j+2}\right)\\ & =a_0{\displaystyle \sum _{j=0}^{N-1}}\left(u_j^2-u_j{u}_{j+1}\right).\end{array}$$

(51)

It suggests that

$$\sum _{j=1}^{N-1}{u}_j^2=\sum _{j=0}^{N-1}\frac{1}{2}(u_j^2+u_{j+1}^2)\ge \sum _{j=1}^{N-1}{u}_j{u}_{j+1}.$$

(52)

From Lemma 5, we obtain

$$\begin{array}{ll}{\parallel u_h\parallel }_1^2=\langle A_{h,2}{u}_h,u_h\rangle \hfill & \ge \langle (I-a_1{L}_{N-1})u_h,u_h\rangle \hfill \\ & ={\displaystyle \sum _{j=1}^{N-1}}\left[(1-2a_1)u_j^2+2a_1{u}_j{u}_{j+1}\right]\\ & \ge -2a_1{\displaystyle \sum _{j=1}^{N-1}}\left(u_j^2-u_j{u}_{j+1}\right).\end{array}$$

(53)

To satisfy (45), we here take

$$\begin{array}{ll}\beta & =-\frac{a_0}{2a_1}=\frac{1-2d\left(g_{1,\lambda }^{(2,\alpha )}-\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }\right)}{2d\left(g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(2,\alpha )}\right)}\\ & =\frac{-g_{1,\lambda }^{(2,\alpha )}+\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }}{g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(2,\alpha )}}+\frac{1}{2d\left(g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(2,\alpha )}\right)}.\end{array}$$

(54)

From Lemma 3, we have

$$\frac{g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(k,\alpha )}}{2}=\frac{|g_{0,\lambda }^{(2,\alpha )}+g_{2,\lambda }^{(k,\alpha )}|}{2}<-g_{1,\lambda }^{(2,\alpha )}+\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }$$

(55)

Combining (54) and (55) with $\omega \le \frac{1}{2}$, we have $\beta >\frac{1}{2}$, and

$$1-\frac{\omega }{\beta }>0.$$

(56)

Based on Theorem 5, we obtain

$${\parallel S_{h,\omega }·T_{TGM}\parallel }_1\le \sqrt{1-\omega /\beta }<1.$$

(57)

□

Remark 3.

It can be seen that the TGM converges linearly from Theorem 6 when $A_{h,2}$ is defined in (18) and $\omega \le 1/2$. In fact, the TGM will also be stable when $\omega >1/2$ but satisfies ${\parallel S_{h,\omega }·T_{TGM}\parallel }_1<1$. The numerical examples in Section 4 also show cases like this.

3.2. Multigrid Convergence for the Third-Order Tempered Fractional Diffusion Discretization

In this subsection, we will repeat the analysis of the multigrid convergence, but now for third-order accurate schemes for the tempered fractional diffusion equation.

To simplify the multigrid analysis of the third-order scheme, we assume $c_r=c_l=c>0$, and we denote by

$$\begin{array}{lll}{\tilde{B}}_h& =& -\frac{c\tau }{2h^{\alpha }}(B_{3,\lambda }+B_{3,\lambda }^T),\\ p_0& =& 1-2d\left(g_{1,\lambda }^{(3,\alpha )}-\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }+{\gamma }_4{e}^{-2h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }\right),\\ p_1& =& p_{-1}=-d\left(g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}\right),and{p}_j=p_{-j}=-d{g}_{j+1,\lambda }^{(3,\alpha )},\end{array}$$

with $j=2,3,4,\dots$. We then find that $A_{h,3}=P_h=I+{\tilde{B}}_h$ is a symmetric Toeplitz matrix of the following form,

$$P_h=I+{\tilde{B}}_h=I-d(B_{3,\lambda }+B_{3,\lambda }^T)=\left(\begin{array}{ccccc}{p}_0& p_1& p_2& \cdots & p_{N-2}\\ p_{-1}& p_0& p_1& \cdots & p_{N-3}\\ p_{-2}& p_{-1}& p_0& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & p_2\\ p_{-(N-3)}& \ddots & \ddots & \ddots & p_1\\ p_{-(N-2)}& p_{-(N-3)}& \cdots & p_{-1}& p_0\end{array}\right),$$

(58)

where $d=\frac{{\displaystyle c\tau }}{{\displaystyle 2h^{\alpha }}}$.

For $\alpha \in (1,2)$, we denote

$$q_1=max\left\{\frac{{\displaystyle \frac{{\alpha }^5}{8}+\frac{7}{12}{\alpha }^4-\frac{5}{8}{\alpha }^3-\frac{49}{12}{\alpha }^2+3\alpha }}{{\displaystyle {\alpha }^3+6{\alpha }^2+11\alpha +6}},\frac{{\displaystyle \frac{{\alpha }^5}{8}+\frac{{\alpha }^4}{3}-\frac{67}{24}{\alpha }^3-\frac{23}{6}{\alpha }^2+\frac{175}{6}\alpha -30}}{{\displaystyle {\alpha }^3+6{\alpha }^2+11\alpha +6}}\right\},$$

(59)

and

$$q_2=min\left\{\frac{{\displaystyle \frac{1}{8}{\alpha }^4+\frac{7}{12}{\alpha }^3+\frac{1}{8}{\alpha }^2-\frac{13}{6}\alpha }}{{\displaystyle {\alpha }^2+5\alpha +8}},\frac{{\displaystyle \frac{{\alpha }^5}{8}+\frac{11}{24}{\alpha }^4-\frac{41}{24}{\alpha }^3-\frac{107}{24}{\alpha }^2+\frac{163}{12}\alpha -8}}{{\displaystyle {\alpha }^3+6{\alpha }^2+11\alpha +6}}\right\},$$

(60)

where $\alpha \in (1,2)$.

The impact of varying $\alpha$ on $q_1$ and $q_2$ is graphically illustrated in Figure 1. Particularly, it can be observed that when $\alpha \in (1.26,1.71)$, $q_1<q_2$ and $(q_1,q_2)\ne \varnothing$.

Again, we will be looking into the matrix properties for the third-order discretization. For this, we will be using the following lemmas and theorems.

Theorem 7

(From [42]).For $\alpha \in (1.26,1.71)$, $\lambda \ge 0$ and $q_1\le {\gamma }_4\le q_2$, the following properties are satisfied, $g_{1,\lambda }^{(3,\alpha )}\le 0,g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}\ge 0,g_{k,\lambda }^{(3,\alpha )}\ge 0(k\ge 3).$

Lemma 6

(From [42]).Let the matrices $B_{3,\lambda }$ and $B_{3,\lambda }^T$ be given by (7) when $k=3$. For $\lambda \ge 0,h>0$ and $\alpha \in (1.26,1.71)$, let $f(\alpha ,\lambda ;x)$ be the generating function of $H=\frac{B_{3,\lambda }+B_{3,\lambda }^T}{2}$. If ${\gamma }_4\in (q_1,q_2)$, we have $f(\alpha ,\lambda ;x)<0$ and $B_{3,\lambda }$ is negative.

For $\alpha \in (1.26,1.71)$, we obtain the following result, which is similar to Lemma 1.

Lemma 7.

For $1.26<\alpha <1.71$, $\lambda \ge 0$, and $q_1\le {\gamma }_4\le q_2$, the matrix $P_h=I-d(B_{3,\lambda }+B_{3,\lambda }^T)$ is diagonally dominant and all eigenvalues of $P_h$ are positive.

Proof. 

From Theorem 7, we obtain

$$\begin{array}{l}{p}_0>1,\\ p_j=p_{-j}<0,j=1,2,3,\dots .\end{array}$$

(61)

Therefore, we find the following result for the i-th row of matrix $P_h$

$$\begin{array}{ll}\sum _{j=-(i-1)}^{N-i-1}{p}_j\ge & {\displaystyle \sum _{j=-\infty }^{\infty }}{p}_j\\ & =|p_0|-2{\displaystyle \sum _{j=1}^{\infty }}\left|p_j\right|\\ & =1-2d\left({\displaystyle \sum _{j=0}^{\infty }}{g}_{j,\lambda }^{(3,\alpha )}-\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }+{\gamma }_4{e}^{-2h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }\right)\\ & =1.\end{array}$$

(62)

It can seen that matrix $P_h$ is strictly diagonally dominant, for $1.26<\alpha <1.71$. From Lemma 3, we conclude that $P_h$ is a symmetric, positive definite matrix. □

Using the same inner products as for the second-order case, i.e.,

$${\langle u,v\rangle }_0=\langle diag\left(P_h\right)u,v\rangle ,{\langle u,v\rangle }_1=\langle P_{h}u,v\rangle ,{\langle u,v\rangle }_2=\langle diag{\left(P_h\right)}^{-1}{P}_{h}u,P_{h}v\rangle .$$

with $\langle ·,·\rangle$ the Euclidean inner product, and also using $\omega \le 1/2$, which satisfies (41) for the third-order scheme (25), we have, similar to Lemma 5, the following lemma

Lemma 8.

For $\alpha \in (1.26,1.71)$, we denote ${\tilde{B}}_{rest}={\tilde{B}}_h+p_1{L}_{N-1}$. Then ${\tilde{B}}_{rest}$ is symmetric, positive definite, when $q_1\le {\gamma }_4\le q_2$.

We thus obtain the following theorem regarding the TGM convergence.

Theorem 8.

Suppose that $P_h$ is defined as in (58) and $\omega \le 1/2$ such that ${\tilde{S}}_{h,\omega }$ satisfies the smoothing condition (42). The convergence factor of the TGM then satisfies

$${\parallel {\tilde{S}}_{h,\omega }·T_{TGM}\parallel }_1<1.$$

(63)

Proof. 

The definition of $u_h$ and $u_H$ are the same as in Theorem 6, with $u_0=u_N=0$. We have the following result, which is similar to (51)

$$\begin{array}{ll}{\parallel u_h-I_H^h{u}_H\parallel }_0^2& =p_0{\displaystyle \sum _{j=0}^{N/2-1}}{\left(u_{2j+1}-\frac{1}{2}{u}_{2j}-\frac{1}{2}{u}_{2j+2}\right)}^2\\ & \le p_0{\displaystyle \sum _{j=0}^{N-1}}\left(u_j^2-u_j{u}_{j+1}\right).\end{array}$$

(64)

This result suggests that

$$\sum _{j=1}^{N-1}{u}_j^2\ge \sum _{j=1}^{N-1}{u}_j{u}_{j+1}.$$

(65)

From Lemma 8, we obtain

$$\begin{array}{ll}{\parallel u_h\parallel }_1^2=\langle P_h{u}_h,u_h\rangle & \ge \langle (I-p_1{L}_{N-1})u_h,u_h\rangle \\ & ={\displaystyle \sum _{j=1}^{N-1}}\left[(1-2q_1)u_j^2+2q_1{u}_j{u}_{j+1}\right]\\ & \ge -2q_1{\displaystyle \sum _{j=1}^{N-1}}\left(u_j^2-u_j{u}_{j+1}\right).\end{array}$$

(66)

To satisfy (45), we will here use

$$\begin{array}{ll}\beta & =-\frac{q_0}{2q_1}=\frac{1-2d\left(g_{1,\lambda }^{(3,\alpha )}-\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }+{\gamma }_4{e}^{-2h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }\right)}{2d\left(g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}\right)}\\ & =\frac{-g_{1,\lambda }^{(3,\alpha )}+\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }+{\gamma }_4{e}^{-2h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }}{g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}}+\frac{1}{2d\left(g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}\right)}.\end{array}$$

(67)

With Lemma 6, we find

$$\frac{g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}}{2}=\frac{|g_{0,\lambda }^{(3,\alpha )}+g_{2,\lambda }^{(3,\alpha )}|}{2}<-g_{1,\lambda }^{(3,\alpha )}+\left({\gamma }_1{e}^{h\lambda }+{\gamma }_2+{\gamma }_3{e}^{-h\lambda }+{\gamma }_4{e}^{-2h\lambda }\right){(1-e^{-h\lambda })}^{\alpha }.$$

(68)

Combining (67) and (68), with $\omega \le \frac{1}{2}$, gives us $\beta >\frac{1}{2}$, and,

$$1-\frac{\omega }{\beta }>0.$$

(69)

Based on Theorem 5, we thus obtain

$${\parallel S_{h,\omega }·T_{TGM}\parallel }_1\le \sqrt{1-\omega /\beta }<1.$$

(70)

□

4. Numerical Example

In this section, we use the V-cycle and provide some numerical results for the tempered fractional diffusion equation, and for the tempered fractional Black–Scholes equation to verify the theoretical multigrid results. So, we will analyze the practical multigrid convergence with a classical multigrid scheme for a number of test cases with tempered fractional derivatives. Here we use the the stopping criterion as follows

$$\frac{\parallel r_h^{\left(k\right)}\parallel }{\parallel r_h^{\left(0\right)}\parallel }<{10}^{-7},$$

(71)

where $r_h^{\left(k\right)}$ is the residual vector after k iterations.

4.1. The Tempered Fractional Diffusion Equation

Example 1.

We first consider the tempered fractional diffusion equation, which is defined as follows

$$\left\{\begin{array}{l}\frac{\partial u(x,t)}{\partial t}=\left({}_0{\mathbb{D}}_x^{\alpha ,\lambda }u(x,t)\right)+\left({}_x{\mathbb{D}}_1^{\alpha ,\lambda }u(x,t)\right)+p(x,t),\\ (x,t)\in (0,1)\times (0,1)\\ u(0,t)=0,u(1,t)=0,t\in (0,1)\\ u(x,0)=e^{-\lambda x}{x}^{1+\alpha },x\in [0,1].\end{array}\right.$$

(72)

In this example, the exact solution for (72) is given by $u(x,t)=e^{-t}{x}^3{(1-x)}^3$, and the source term

$$\begin{array}{ll}p(x,t)=& -e^{-t}{x}^3{(1-x)}^3-e^{-\lambda x-t}\left({}_0{D}_x^{\alpha }\left[e^{\lambda x}(x^3-3x^4+3x^5-x^6)\right]\right)\\ & -e^{\lambda x-t}\left({}_x{D}_1^{\alpha }\left[e^{-\lambda x}({(1-x)}^3-3{(1-x)}^4+3{(1-x)}^5-{(1-x)}^6)\right]\right),\end{array}$$

is prescribed accordingly.

We compute $p(x,t)$ by using the following formulae

$$\begin{array}{c}\hfill {}_0{D}_x^{\alpha }\left(e^{\lambda x}{x}^m\right){=}_0{D}_x^{\alpha }\left(\sum _{n=0}^{\infty }\frac{{\lambda }^n}{n!}{x}^{n+m}\right)=\sum _{n=0}^{\infty }\frac{{\lambda }^n\mathsf{\Gamma }(n+m+1)}{n!\mathsf{\Gamma }(n+m-\alpha +1)}{x}^{n+m-\alpha },\end{array}$$

(73)

and

$$\begin{array}{ll}{}_x{D}_1^{\alpha }\left(e^{-\lambda x}{(1-x)}^m\right)& ={}_x{D}_1^{\alpha }\left({\displaystyle \sum _{n=0}^{\infty }}\frac{{\lambda }^n{e}^{-\lambda }}{n!}{(1-x)}^{n+m}\right)\\ & =e^{-\lambda }{\displaystyle \sum _{n=0}^{\infty }}\frac{{\lambda }^n\mathsf{\Gamma }(n+m+1)}{n!\mathsf{\Gamma }(n+m-\alpha +1)}{(1-x)}^{n+m-\alpha }.\end{array}$$

(74)

In the numerical experiment, we take $\alpha =1.5\in (1.26,1.71)$ (which is in the interval for which we have proven multigrid convergence) and $\lambda =0.5$ for the tempered fractional diffusion Equation (72). We use ${\gamma }_3=0.01\in (a_1,a_2)$ for the second-order scheme, and ${\gamma }_4=-0.03$ for the third-order scheme (36) with $N=M$ when $k=2$ and $k=3$.

Table 1. $L^2$ errors for (72) with different $\omega$ by the second-order scheme.

M	Error	Order	$\mathit{\omega }=0.4$	$\mathit{\omega }=0.5$	$\mathit{\omega }=0.6$	$\mathit{\omega }=0.7$
			Iter	Iter	Iter	Iter
$2^5$	$1.18\times {10}^{-5}$		8	7	6	5
$2^6$	$3.04\times {10}^{-6}$	1.96	8	6	5	4
$2^7$	$7.72\times {10}^{-7}$	1.98	7	6	5	4
$2^8$	$1.95\times {10}^{-7}$	1.99	7	5	4	4
$2^9$	$4.89\times {10}^{-8}$	1.99	6	5	4	3

and

Table 2. $L^2$ errors for (72) with $\omega =0.5$ by the second-order scheme.

M	Error	V(0,1)	V(1,0)	V(1,1)
		Iter	Iter	Iter
$2^5$	$1.18\times {10}^{-5}$	13	12	7
$2^6$	$3.04\times {10}^{-6}$	10	12	6
$2^7$	$7.72\times {10}^{-7}$	8	11	6
$2^8$	$1.95\times {10}^{-7}$	6	10	5
$2^9$	$4.89\times {10}^{-8}$	6	8	5

present the corresponding $L^2$ discretization errors and the number of multigrid iterations based on the second-order scheme to reach the tolerance.

Table 3. $L^2$ errors for (72) with different $\omega$ by the third-order scheme.

M	Error	Order	$\mathit{\omega }=0.4$	$\mathit{\omega }=0.5$	$\mathit{\omega }=0.6$	$\mathit{\omega }=0.7$
			Iter	Iter	Iter	Iter
$2^5$	$8.05\times {10}^{-6}$		12	9	7	5
$2^6$	$1.14\times {10}^{-6}$	2.82	11	8	6	5
$2^7$	$1.57\times {10}^{-7}$	2.86	10	7	6	4
$2^8$	$2.16\times {10}^{-8}$	2.87	9	7	5	4
$2^9$	$3.02\times {10}^{-9}$	2.84	8	6	5	4

and

Table 4. $L^2$ errors for (72) with $\omega =0.5$ by the third-order scheme.

M	Error	V(0,1)	V(1,0)	V(1,1)
		Iter	Iter	Iter
$2^5$	$8.05\times {10}^{-6}$	16	17	9
$2^6$	$1.14\times {10}^{-6}$	15	15	8
$2^7$	$1.57\times {10}^{-7}$	13	14	7
$2^8$	$2.16\times {10}^{-8}$	12	13	7
$2^9$	$3.02\times {10}^{-9}$	10	11	6

show the corresponding $L^2$ discretization errors and the number of multigrid iterations for the third-order scheme. Here we use $V({\nu }_1,{\nu }_2)$ to denote the multigrid V-cycle, where ${\nu }_1$ denotes the number of pre-smoothing steps and ${\nu }_2$ the number of post-smoothing steps. From the result, we clearly see the h-independent convergence of multigrid for these involved tempered fractional operators.

4.2. The Tempered Fractional Black–Scholes Equations

In this subsection, we consider the following tempered fractional Black–Scholes equation

$$\begin{array}{cc}\hfill & \hfill \frac{\partial u(x,t)}{\partial t}+a·\frac{\partial u(x,t)}{\partial x}+b{·}_{B_d}{D}_x^{\alpha ,{\lambda }_1}u(x,t)+d{·}_x{D}_{B_u}^{\alpha ,{\lambda }_2}u(x,t)=c·u(x,t)+b{\lambda }_1^{\alpha }u(x,t)+d{\lambda }_2^{\alpha }u(x,t),\hfill \end{array}$$

(75)

where $\alpha \in (1,2)$, the parameters $b,d,c,{\lambda }_1$ and ${\lambda }_2$ are all non-negative. Here, we show, experimentally, that the proposed schemes are robust and accurate, without any proof of stability/convergence.

We consider the following problem, with a source term $p(x,t)$, which was added to test the numerical scheme, as follows,

$$\left\{\begin{array}{l}\frac{\partial u(x,t)}{\partial t}+a\frac{\partial u(x,t)}{\partial x}+b_{B_d}{D}_x^{\alpha ,{\lambda }_1}u(x,t)+d_x{D}_{B_u}^{\alpha ,{\lambda }_2}u(x,t)\\ =c·u(x,t)++b{\lambda }_1^{\alpha }u(x,t)+d{\lambda }_2^{\alpha }u(x,t)+p(x,t),(x,t)\in (B_d,B_u)\times (0,T)\\ u(B_d,t)=0,u(B_u,t)=0,t\in (0,T)\\ u(x,T)=S\left(x\right),x\in (B_d,B_u),x\in (B_d,B_u).\end{array}\right.$$

(76)

4.2.1. Multigrid Results with the Second-Order Scheme

For this test case, we take $t_j=(N-j)\tau$, $0\le t_j\le T,j=0,\cdots ,N$ and $x_i=B_d+ih,$$B_d\le x_i\le B_u,$ where $\tau =T/N$ and $h=(B_u-B_d)/M$. We use the tempered-WSGD operators

$${}_L{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_1}{=}_L{\mathcal{D}}_{h,-1,0,,1}^{\alpha ,{\gamma }_1,{\gamma }_2,{\gamma }_3}and{}_R{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_2}{=}_R{\mathcal{D}}_{h,-1,0,1}^{\alpha ,{\gamma }_1,{\gamma }_2,{\gamma }_3},$$

for the tempered fractional derivatives, and the second-order scheme for the first-order spatial derivative, so that the discretization for (76) reads

$$\begin{array}{c}\hfill \frac{\partial u_i}{\partial t}+a\frac{u_{i+1}-u_{i-1}}{2h}+b\left({}_L{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_1}{u}_i\right)+d\left({}_R{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_2}{u}_i\right)=c·u_i+p_i+O\left(h^2\right),\end{array}$$

(77)

where $u_i$ is the solution of (76) when $x=x_i$, and $p_i=p(x_i,t)$.

The discretization in time is based on the Crank–Nicolson scheme, which, for (77), reads,

$$\begin{array}{l}-\frac{u_i^{j+1}-u_i^j}{\tau }+a\frac{u_{i+1}^{j+\frac{1}{2}}-u_{i-1}^{j+\frac{1}{2}}}{2h}+b\left({}_L{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_1}{u}_i^{j+\frac{1}{2}}\right)+d\left({}_R{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_2}{u}_i^{j+\frac{1}{2}}\right)\\ =c·u_i^{j+\frac{1}{2}}+p_i^{j+\frac{1}{2}}+O({\tau }^2+h^2),\end{array}$$

(78)

with $u_i^j$ the solution of (76) at point $(x_i,t_j)$, and $p_i^{j+\frac{1}{2}}=p(x_i,t_{j+\frac{1}{2}})$.

The numerical discretization in space and time can be written as follows,

$$\begin{array}{l}(1+\frac{\tau }{2}c)U_i^{j+1}-\frac{\tau }{2}\left[a\frac{U_{i+1}^{j+1}-U_{i-1}^{j+1}}{2h}+b_L{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_1}{U}_i^{j+1}+d_R{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_2}{U}_i^{j+1}\right]\\ =(1-\frac{\tau }{2}c)U_i^j+\frac{\tau }{2}\left[a\frac{U_{i+1}^j-U_{i-1}^j}{2h}+b_L{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_1}{U}_i^j+d_R{\mathcal{D}}_{h,2}^{\alpha ,{\lambda }_2}{U}_i^j\right]-P_i^{j+\frac{1}{2}},\end{array}$$

(79)

where $U_i^j$ is the numerical solution for (76) at point $(x_i,t_j)$, and $P_i^{j+\frac{1}{2}}=\frac{\tau }{2}(p_i^j+p_i^{j+1})$.

The matrix equation for (79) is given by

$$A{U}^{j+1}={\widehat{P}}^{j+1},$$

(80)

where

$$A=\left(\begin{array}{ccccc}1+\frac{\tau }{2}c& -\frac{a\tau }{4h}& & & \\ \frac{a\tau }{4h}& 1+\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ ⋮& \frac{a\tau }{4h}& 1+\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ \cdots & \ddots & \ddots & \ddots & -\frac{a\tau }{4h}\\ \cdots & \cdots & \cdots & \frac{a\tau }{4h}& 1+\frac{\tau }{2}c\end{array}\right)-\frac{\tau }{2}\left(b{B}_{2,{\lambda }_1}+d{B}_{2,{\lambda }_2}^{\mathrm{T}}\right),$$

$$Q=\left(\begin{array}{ccccc}1-\frac{\tau }{2}c& -\frac{a\tau }{4h}& & & \\ \frac{a\tau }{4h}& 1-\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ ⋮& \frac{a\tau }{4h}& 1-\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ \cdots & \ddots & \ddots & \ddots & -\frac{a\tau }{4h}\\ \cdots & \cdots & \cdots & \frac{a\tau }{4h}& 1-\frac{\tau }{2}c\end{array}\right)+\frac{\tau }{2}\left(b{B}_{2,{\lambda }_1}+d{B}_{2,{\lambda }_2}^{\mathrm{T}}\right),$$

and

$${\widehat{P}}^{j+1}=Q{U}^i+{\left(P_1^{j+\frac{1}{2}},P_2^{j+\frac{1}{2}},\dots ,P_{M-1}^{j+\frac{1}{2}}\right)}^T.$$

4.2.2. Multigrid Results for the Third-Order Scheme

Using the tempered-WSGD operators, ${}_L{\mathcal{D}}_{h,3}^{\alpha ,{\lambda }_1}{=}_L{\mathcal{D}}_{h,-1,0,\alpha -1,1}^{\alpha ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_4}$ and ${}_R{\mathcal{D}}_{h,3}^{\alpha ,{\lambda }_2}=$${}_R{\mathcal{D}}_{h,-1,0,\alpha -1,1}^{\alpha ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_4}$, for the tempered fractional derivatives, and the fourth-order scheme for the first-order spatial derivative, we obtain the following space discretization for (76),

$$\begin{array}{l}\frac{\partial u_i}{\partial t}+a\frac{8(u_{i+1}-u_{i-1})-(u_{i+2}-u_{i-2})}{12h}+b\left({}_L{\mathcal{D}}_h^{\alpha ,{\lambda }_1}{u}_i\right)+d\left({}_R{\mathcal{D}}_h^{\alpha ,{\lambda }_2}{u}_i\right)\\ =c·u_i+p_i+O\left(h^3\right).\end{array}$$

(81)

For (81), the Crank–Nicolson time discretization can now be written as

$$\begin{array}{l}(1+\frac{\tau }{2}c)U_i^{j+1}-\frac{\tau }{2}\left[a\frac{8(U_{i+1}^{j+1}-U_{i-1}^{j+1})-(U_{i+2}^{j+1}-U_{i-2}^{j+1})}{12h}+b_L{\mathcal{D}}_{h,3}^{\alpha ,{\lambda }_1}{U}_i^{j+1}+d_R{\mathcal{D}}_{h,3}^{\alpha ,{\lambda }_2}{U}_i^{j+1}\right]\\ =(1-\frac{\tau }{2}c)U_i^j+\frac{\tau }{2}\left[a\frac{8(U_{i+1}^j-U_{i-1}^j)-(U_{i+2}^j-U_{i-2}^j)}{12h}+b_L{\mathcal{D}}_{h,3}^{\alpha ,{\lambda }_1}{U}_i^j+d_R{\mathcal{D}}_{h,3}^{\alpha ,{\lambda }_2}{U}_i^j\right]-P_i^{j+\frac{1}{2}}.\end{array}$$

(82)

The matrix form for (82) is

$$\widehat{A}{U}^{j+1}={\tilde{P}}^{j+1},$$

(83)

where

$$\widehat{A}=\left(\begin{array}{ccccc}1+\frac{\tau }{2}c& -\frac{a\tau }{4h}& & & \\ \frac{a\tau }{4h}& 1+\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ ⋮& \frac{a\tau }{4h}& 1+\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ \cdots & \cdots & \ddots & \ddots & -\frac{a\tau }{4h}\\ \cdots & \cdots & \cdots & \frac{a\tau }{4h}& 1+\frac{\tau }{2}c\end{array}\right)-\frac{\tau }{2}\left(b{B}_{3,{\lambda }_1}+d{B}_{3,{\lambda }_2}^{\mathrm{T}}\right),$$

$$\widehat{Q}=\left(\begin{array}{ccccc}1-\frac{\tau }{2}c& -\frac{a\tau }{4h}& & & \\ \frac{a\tau }{4h}& 1-\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ ⋮& \frac{a\tau }{4h}& 1-\frac{\tau }{2}c& -\frac{a\tau }{4h}& \\ \cdots & \ddots & \ddots & \ddots & -\frac{a\tau }{4h}\\ \cdots & \cdots & \cdots & \frac{a\tau }{4h}& 1-\frac{\tau }{2}c\end{array}\right)+\frac{\tau }{2}\left(b{B}_{3,{\lambda }_1}+d{B}_{3,{\lambda }_2}^{\mathrm{T}}\right),$$

and

$${\tilde{P}}^{j+1}=\widehat{Q}{U}^i+{\left(P_1^{j+\frac{1}{2}},P_2^{j+\frac{1}{2}},\dots ,P_{M-1}^{j+\frac{1}{2}}\right)}^T.$$

Example 2.

We finally consider the following tempered fractional model

$$\left\{\begin{array}{l}\frac{\partial u(x,t)}{\partial t}+a\frac{\partial u(x,t)}{\partial x}+b\left({}_0{\mathfrak{D}}_x^{\alpha ,{\lambda }_1}u(x,t)\right)+d\left({}_0{\mathfrak{D}}_x^{\alpha ,{\lambda }_1}u(x,t)\right)=c·u(x,t)+p(x,t),\\ (x,t)\in (0,1)\times (0,T)\\ u(0,t)=0,u(1,t)=0,t\in (0,T)\\ u(x,T)=S\left(x\right),x\in (0,1),\end{array}\right.$$

(84)

where

$$\begin{array}{ll}p(x,t)=& -(1+a{\lambda }_1^{\alpha }+b{\lambda }_2^{\alpha }+p)u(x,t)+3a{e}^{T-t}{x}^2{(1-x)}^2(1-2x)\\ & +b{e}^{-{\lambda }_{1}x+(T-t)}\left({}_0{D}_x^{\alpha }{e}^{{\lambda }_{1}x}u(x,t)\right)+c{e}^{{\lambda }_{2}x+(T-t)}\left({}_x{D}_1^{\alpha }{e}^{-{\lambda }_{2}x}u(x,t)\right).\end{array}$$

The exact solution of the above equation is given by $u(x,t)=e^{-\lambda x+(T-t)}{x}^3(1-x)$.

We will use the following parameters in the numerical tests, $b=c=d=1$ and $a=-0.5$. We choose $\alpha =1.8,{\lambda }_1=0.5,{\lambda }_2=1,$ and $\tau ={10}^{-4}$ in this case.

Table 5. $L^2$ errors for (84) with different $\omega$ by the second-order scheme.

M	Error	Order	$\mathit{\omega }=0.4$	$\mathit{\omega }=0.5$	$\mathit{\omega }=0.6$	$\mathit{\omega }=0.7$
			Iter	Iter	Iter	Iter
$2^6$	$1.63\times {10}^{-5}$		9	7	6	6
$2^7$	$4.10\times {10}^{-6}$	1.99	8	7	6	6
$2^8$	$1.03\times {10}^{-6}$	2.00	7	6	5	6
$2^9$	$2.57\times {10}^{-7}$	2.00	6	5	4	6
$2^{10}$	$6.44\times {10}^{-8}$	2.00	5	4	4	5

and

Table 6. $L^2$ errors for (84) with $\omega =0.5$ by the second-order scheme.

M	Error	V(0,1)	V(1,0)	V(1,1)
		Iter	Iter	Iter
$2^6$	$1.63\times {10}^{-5}$	12	13	7
$2^7$	$4.10\times {10}^{-6}$	10	12	7
$2^8$	$1.03\times {10}^{-6}$	9	11	6
$2^9$	$2.57\times {10}^{-7}$	7	10	5
$2^{10}$	$6.44\times {10}^{-8}$	6	8	4

present the corresponding $L^2$ discretization errors and the number of multigrid iterations for the second-order scheme.

Table 7. $L^2$ errors for (84) with different $\omega$ by the third-order scheme.

M	Error	Order	$\mathit{\omega }=0.4$	$\mathit{\omega }=0.5$	$\mathit{\omega }=0.6$	$\mathit{\omega }=0.7$
			Iter	Iter	Iter	Iter
$2^6$	$3.55\times {10}^{-6}$		9	7	6	5
$2^7$	$4.69\times {10}^{-7}$	2.92	9	7	5	5
$2^8$	$6.04\times {10}^{-8}$	2.96	9	7	5	4
$2^9$	$7.63\times {10}^{-9}$	2.99	8	6	5	4
$2^{10}$	$9.38\times {10}^{-10}$	3.02	8	6	5	4

and

Table 8. $L^2$ errors for (84) with $\omega =0.5$, by the third-order scheme.

M	Error	V(0,1)	V(1,0)	V(1,1)
		Iter	Iter	Iter
$2^6$	$3.55\times {10}^{-6}$	17	18	9
$2^7$	$4.69\times {10}^{-7}$	17	18	9
$2^8$	$6.04\times {10}^{-8}$	16	17	9
$2^9$	$7.63\times {10}^{-9}$	15	16	8
$2^{10}$	$9.38\times {10}^{-10}$	14	15	8

show the corresponding $L^2$ errors and the number of multigrid iterations based on the third-order scheme.

5. Conclusions

In this paper, we analyzed a classical multigrid method for second- and third-order numerical schemes for the tempered fractional diffusion equation. We have detailed the classical multigrid components, like the damped Jacobi smoothing iteration, and the direct coarse grid approximation, which is based on the second- and third-order discrete schemes. A focus of this paper was the multigrid convergence analysis, which was based on the properties of the occurring discretization matrices.

Moreover, we have also shown that the multigrid method converged very well for the tempered fractional Black–Scholes equation.

Obviously, the numerical schemes presented in this paper are computationally highly accurate and efficient.