On a Multigrid Method for Tempered Fractional Diffusion Equations

: 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 conﬁrm accuracy and effectiveness of the proposed method.


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. Sections 3.1 and 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.

Numerical Schemes for the Tempered Fractional Diffusion Equation
We consider the following tempered fractional diffusion equation The Riemann-Liouville tempered fractional derivatives, that we encounter in this equation, are defined as follows.
Definition 1 (See [31]). For α ∈ (n − 1, n), let u(x) be (n − 1)-times continuously differentiable on (a, b) with its nth derivative integrable on any subinterval of [a, b], and λ ≥ 0. Then, the left Riemann-Liouville tempered fractional derivative of order α is defined as (2) the right Riemann-Liouville tempered fractional derivative of order α is defined as where 'a' and 'b' can be extended to −∞ 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(x) on a bounded interval [a, b], and the function u(x) will be zero extended for x < a or x > b, so that u(x) ∈ L 1 (R), and a D α+l,λ x u(x), x D α+l,λ b u(x) and their Fourier transforms belong to L 1 (R). The α-th order left and right Riemann-Liouville tempered fractional derivatives of u(x) at grid point x can then be approximated by tempered-WSGD operators L D α,λ h,k u and R D α,λ h,k u, as follows [31,42] for details. The second-and third-order operators are given in Sections 2.1 and 2.2.
Let the equidistant time partition, t j = jτ (0 ≤ t j ≤ T, j = 0, . . . , N), and spatial grid, Using the high-order tempered-WSGD operators L D α,λ h,k u and R D α,λ h,k u (as explained in Remark 1), high-order scheme for the first-order spatial derivative with δ k,x u = ∂ x u + O(h k ), and a Crank-Nicolson discretization in time, the numerical scheme for (1) reads where u j i represents the solution of (1) at the point ( Denote by U j i the solution of the numerical scheme for (1) at point (x i , t j ). The numerical scheme can now be written as We will use the following notations, for vector U n = (U n 1 , U n 2 , . . . , The corresponding matrix form of (6) then readŝ and

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, where with γ j satisfying the following conditions, and the weights g k,λ , k = 0, · · · , j + 1, are given by We will present the second-order scheme for the tempered fractional diffusion equation in detail here. Using the tempered-WSGD operators, for the tempered fractional derivatives, and the second-order scheme for the first-order spatial derivative, the numerical scheme can now be written as, The corresponding matrix form of (16) then reads withÂ j+1 h,2 ,F n+ 1 2 as defined in (9), (11) when k = 2, H 2 = tridiag{−1, 0, 1}, and 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.

Third-Order Discrete Scheme for the Tempered Fractional Diffusion Equation
In this work, we also consider the third-order operators. They are defined as and where with γ j satisfying the following conditions The weights, g (3,α) k,λ , k = 0, · · · , j + 1, are found to be Using the tempered-WSGD operators, L D α,λ h,3 = L D α,γ 1 ,γ 2 ,...,γ 4 h,−1,0,1,2 and R D α,λ h,3 = R D α,γ 1 ,γ 2 ,...,γ 4 h,−1,0,1,2 , for the tempered fractional derivatives, and the fourth-order scheme for the first-order spatial derivative, we find the following numerical discretization for (1) The matrix form now looks as follows withÂ j+1 h,3 ,F n+ 1 2 as defined in (9), (11) and with 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 and the generating function, We obtain the stability for the numerical scheme (25), based on the theorem below.
Lemma 2 (From [42]). Let the matrices B 3,λ and B T 3,λ be given via the numerical scheme (7). For λ ≥ 0, h > 0 and α ∈ [1, 2], if we can find (analytically, or with the help of numerical techniques) values of γ j for which the generating functions f (α, λ; x) of B λ are negative, then the eigenvalues of the matrix B 3,λ are negative too.
In a similar way as in [42], we have the following theorem.
Theorem 2. If, for 1 < α < 2, the generating functions f (α, λ; 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) × (0, T), which can be zero extended for x < a or x > b, so that u ∈ L 1 (0, T; R), and a D α+3,λ x u and its Fourier transform also belong to L 1 (0, T; R).
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.

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Ũ For the above description, the notation is as follows: • ν 1 , ν 2 denote the number of smoothing steps. We will use ν i = 0, 1, 2. • The classical fine-to-coarse restriction operator I H h is employed, • The scaled transpose of the restriction is the coarse-to-fine interpolation operator, i.e., 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 Then We can easily generalize the classical Jacobi iteration by introducing a relaxation parameter ω, in the standard way, i.e., where diag(A h,k ) represents the main diagonal of matrix A h,k . The smoother then becomes, 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(M log M) 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).

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 =Â h,k . We denote by j+1,λ , with j = 2, 3, 4, . . .. Then, we find that A h,2 = I + B h is a symmetric Toeplitz matrix, of the following form, where d = cτ 2h α . We need the following lemmas regarding the properties of our matrix, in order to prove multigrid convergence. (From [31]). Using the notation,
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 . (18) is diagonally dominant and all eigenvalues of A h,2 are positive.
Proof. From Lemma 1, we obtain, Therefore, we have the following result for the ith row of the matrix A h,2 It can seen that the matrix A h,2 is strictly diagonally dominant for 1 < α < 2. From Lemma (3), we then conclude that A h,2 is a symmetric, positive definite matrix.
We will use the following inner products, Here ·, · is the Euclidean inner product.
Then, β ≥ ω and the convergence factor of the TGM satisfies In other words, we just need to find a suitable β-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) × (N − 1) one-dimensional discrete Laplacian matrix. Then L N−1 is also a symmetric positive definite Toeplitz matrix.
Proof. Since both B h and L N−1 are symmetric Toeplitz, B rest is also symmetric Toeplitz. We have For the k-th row, we then obtain that Hence, B rest is strictly diagonally dominant, and we know that B rest is positive sincẽ b 0 > 0.

From Lemma 5, it's easy to see that
Now we are ready to provide the proof for the TGM convergence, in the following theorem.
The numerical examples in Section 4 also show cases like this.

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 bỹ , and p j = p −j = −dg (3,α) j+1,λ , with j = 2, 3, 4, . . .. We then find that A h,3 = P h = I +B h is a symmetric Toeplitz matrix of the following form, where d = cτ 2h α .
For α ∈ (1, 2), we denote and where α ∈ (1, 2). The impact of varying α on q 1 and q 2 is graphically illustrated in Figure 1. Particularly, it can be observed that when α ∈ (1.26, 1.71), q 1 < q 2 and (q 1 , q 2 ) = ∅. 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.
Therefore, we find the following result for the i-th row of matrix P h It can seen that matrix P h is strictly diagonally dominant, for 1.26 < α < 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., with ·, · the Euclidean inner product, and also using ω ≤ 1/2, which satisfies (41) for the third-order scheme (25), we have, similar to Lemma 5, the following lemma We thus obtain the following theorem regarding the TGM convergence.
Theorem 8. Suppose that P h is defined as in (58) and ω ≤ 1/2 such thatS h,ω satisfies the smoothing condition (42). The convergence factor of the TGM then satisfies S h,ω · T TGM 1 < 1.
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) This result suggests that From Lemma 8, we obtain To satisfy (45), we will here use With Lemma 6, we find Combining (67) and (68), with ω ≤ 1 2 , gives us β > 1 2 , and, Based on Theorem 5, we thus obtain

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 where r (k) h is the residual vector after k iterations.

The Tempered Fractional Diffusion Equation
Example 1. We first consider the tempered fractional diffusion equation, which is defined as follows 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 is prescribed accordingly.
We compute p(x, t) by using the following formulae In the numerical experiment, we take α = 1.5 ∈ (1.26, 1.71) (which is in the interval for which we have proven multigrid convergence) and λ = 0.5 for the tempered fractional diffusion Equation (72). We use γ 3 = 0.01 ∈ (a 1 , a 2 ) for the second-order scheme, and γ 4 = −0.03 for the third-order scheme (36) with N = M when k = 2 and k = 3. Tables 1 and 2 present the corresponding L 2 discretization errors and the number of multigrid iterations based on the second-order scheme to reach the tolerance. Tables 3 and 4 show the corresponding L 2 discretization errors and the number of multigrid iterations for the third-order scheme. Here we use V(ν 1 , ν 2 ) to denote the multigrid V-cycle, where ν 1 denotes the number of pre-smoothing steps and ν 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.   3.02 × 10 −9 10 11 6

The Tempered Fractional Black-Scholes Equations
In this subsection, we consider the following tempered fractional Black-Scholes equation where α ∈ (1, 2), the parameters b, d, c, λ 1 and λ 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,