The finite volume element method for the two-dimensional space-fractional convection–diffusion equation

We develop a fully discrete finite volume element scheme of the two-dimensional space-fractional convection–diffusion equation using the finite volume element method to discretize the space-fractional derivative and Crank–Nicholson scheme for time discretization. We also analyze and prove the stability and convergence of the given scheme. Finally, we validate our theoretical analysis by data from three examples.


Introduction
Fractional differential equations are generalizations of the integer-order differential equations; they contain noninteger-order derivatives and can effectively describe the memory and genetic properties of a variety of substances. Such equations play an increasingly important role in mathematical physics, mechanical applications, biological engineering, electronic science and technology, control theory, financial mathematics, and other fields [1][2][3][4][5].
Abnormal diffusion in physics originally developed from the random walk model [6][7][8]. Fractional convective diffusion equations are a powerful tool to simulate various abnormal diffusion phenomena. We consider the following one-dimensional random walk model. First, we denote by n the random walk step size and assume that the waiting time obeys an exponential distribution with mean τ n and the probability density function λ(x) of the jump step obeys a stable distribution. We define the probability density function W (x, t) of the particle at time t at the position x. Then where (t) = 1 -t 0 w(t ) dt is the probability that the particle does not jump within the time interval (0, t), η(x, t ) is the probability that the particle is at the position of x at time t multiplied by the probability that the particle jumps xx step in time tt , namely, the Markov supplement process in the random process η( where δ(x)δ(t) is the initial condition, and ψ(x, t) = w(t)λ(x). Then taking the Laplace transform of (t) and the Fourier and Laplace transforms of x and t in W (x, t) and η(x, t), respectively, we get .
On the basis of the above, we discuss the random movement of particles in a fluid with mean velocity V located in the porous medium. We assume that the particles are stuck in place while waiting for the next jump. Since the fluid is located in the porous medium, the particles may be confined to the pores, and so our hypothesis is reasonable. We denote by n the number of random walk steps, by α the ratio index, by σ α n the diffusion coefficient, and by ϑ the tilt angle of the transition probability density function. Then we introduce the similar variable ξ = x -Vt to obtain the corresponding probability density function φ(x, t) = ψ(x -Vt, t). In this case, the special function λ(κ) = exp -|κ| α exp i sign(κ)ϑπ/2 = 1 -|κ| α exp i sign(κ)ϑπ/2 n + O( 1 n ). After transformation and derivation, W (κ, s) can be rewritten as where K α = σ α τ , c + = c + (α, ϑ) = sin((α-ϑ)π /2) sin(απ) , and c -= c -(α, ϑ) = sin((α+ϑ)π /2) sin(απ) . As n → ∞, we get s W n (κ, s) = -K α c + (-iκ) α + c -(iκ) α W n (κ, s) -iV κ W n (κ, s).
The main work we have done in this paper is obtaining a fully discrete finite volume element scheme of the two-dimensional Riesz space fractional convection-diffusion equation (3) by using the CN-FVE method. Also, we prove the uniqueness, stability, and convergence of the scheme in the L 2 -norm. To save time, we consider the characteristics of Toplitz matrices, and use the FAST-BICGSTAB algorithm to solve the numerical examples, so as to verify the accuracy of our theoretical analysis.
The remainder of this paper is organized as follows. In Sect. 2, we give the fractional convection-diffusion equation (3) and its fully discrete finite volume element scheme. In Sect. 3, we analyze and prove the stability and convergence of the scheme. In Sect. 4, we give three numerical examples to verify the accuracy of the theoretical analysis. Finally in Sect. 5, we provide conclusions of this paper.
Based on the definition of the Riesz fractional derivative, we can write the following equivalent form of equation (3): Let M be a positive integer, let τ = T M , and let t m = mτ (m = 0, 1, . . . , M). We use the CN-scheme to discretize the time derivative in equation (3) and get Next, let N x and N y be positive integers, and let h x = L 1 N x +1 and h y = L 2 N y +1 . We divide the region into a uniform grid with grid nodes x i = ih x (0 ≤ i ≤ N x + 1) and y j = jh y (0 ≤ j ≤ N y + 1). Denote the element i, Besides, for i = 1, . . . , N x , j = 1, . . . , N y , m = 1, . . . , M, let u m i,j be the finite volume element approximation of u(x i , y j , t m ).
We take the trial function space S h ( ) ⊂ H 1 0 ( ) as the linear element space with respect to the space subdivision mentioned above: for u h ∈ S h ⊂ H 1 0 ( ), we have S h = {u h ∈ C( ) : u h | i,j ∈ P 1 ; u h | ∂ = 0}, where P 1 is the set of all linear polynomials on i,j . Then the approximate solution u h (x, y, t m ) of equation (3) can be expressed as follows: where φ x k (x) are linear basis functions in the direction of x, expressed as follows: e lsewher e, (6) and the definition of φ y l (y) in the y direction is similar.
]. Accordingly, we choose the test function space as the piecewise constant function space (4) and omit the time-truncated error term. Then we get the following variational form: We expand formula (8) to obtain the following fully discrete finite volume element scheme: To facilitate the calculation, we introduce the following lemma.
Let u m and F m be N := N x N y -dimensional vectors, defined as follows: wherē Then we define the following N x th-order mass matrix A x , N x th-order stiffness matrix B x , and N x th-order matrix C x : where T 1-β,N x is a Toeplitz matrix of the following form: Similarly, we can get three matrices A y , B y , and C y in the y direction defined as follows: We introduce the following definition and lemmas to give the matrix form of the fully discrete scheme (9).

Definition 1 ([41]) For
A ∈ R m×n and B ∈ R r×s , their Kronecker product is the partitioned mr × ns matrix

Lemma 4 ([41])
The Kronecker product satisfies the following bilinear and associative properties: where A, B, and C are matrices, and k is a constant.
Thus the finite volume element scheme (9) can be expressed in the following matrix form: where η β =

Stability and convergence analysis
In this section, to investigate the stability and convergence of scheme (22), we need the following lemmas.
Lemma 11 ([17]) Let A x and A y be defined by (16) and (19), respectively. Then A y ⊗ A x is symmetric and positive definite, and for all v ∈ R N , we have Lemma 12 For 0 < β, γ < 1, the stiffness matrices A y ⊗ B x and B y ⊗ A x are positive definite, that is, Proof According to Lemma 8, to prove that the matrices A y ⊗ B x and B y ⊗ A x are positive definite, we just have to prove that their symmetric parts According to (16) and (17), we find that the elements g i,j of G x satisfy By Lemma 9 we can infer that Thus the matrix G x is a symmetric strictly diagonally dominant matrix with positive diagonal elements. It follows that the matrix G x is positive definite. Similarly, we can prove that the matrix G y is symmetrically positive definite. The matrices A x and A y are also symmetrically positive definite. Therefore by Lemma 10 we can infer that the matrices G 1 and G 2 are positive definite. This completes the proof.

Lemma 13
Let C x and C y be defined by (16) and (19), respectively. Then the matrices

Because of the symmetry of
and since C x is an antisymmetric matrix, we obtain This completes the proof.
Next, according to Lemma 11, a two-dimensional weighted discrete norm is defined as and v m We can easily verify that · A and · l 2 are equivalent, and by Lemma 11 we have Theorem 1 For any 0 < β, γ < 1, the finite volume element scheme (22) is uniquely solvable.
Proof It is obvious from equation (9) that the scheme is solvable. Next, to prove that scheme (22) is uniquely solvable, we need to prove that the homogeneous system of equa- has only zero solution. Now left multiplying both sides of this system by h x h y (u m ) T , we organize it into the following form: By Lemmas 12 and 13 we know so u m 2 A ≤ 0, that is, u m 2 A = 0, and thus u m = 0. This completes the proof.
Theorem 2 For any 0 < β, γ < 1, the finite volume element scheme (22) is unconditionally stable in the sense of the discrete norm, that is, Proof First, we organize (22) into the following form: Then left multiplying both sides of equation (37) by h x h y (u m + u m-1 ) T , by Lemmas 11, 12, and 13 we get where in the last two steps we used the Cauchy-Schwarz inequality and formula (33), respectively. Subtracting u m + u m-1 A from both sides of this equation, we have By iteration we get Again by formula (33) we obtain This completes the proof. and Lemma 15 ([15]) Let W (x) ∈ C β+1 [0, L] for some 0 < β < 1. Then (45) Proof Applying the Lagrange mean value theorem, we get This completes this proof.

Corollary 2
Combining Lemmas 14 and 16, we obtain Proof By Lemma 14 we have and by Lemma 16 we get This completes this proof. (3), let u m i,j = u(x i , y j , t m ) be the numerical solution of the finite volume element scheme (22).

Then we have the error estimate
Proof By the previous analysis the exact solution U m i,j = u(x i , y j , t m ) of the form (22) has a local truncation error Let e m = U m -u m be the global truncation error of the initial value e 0 = 0. We can obtain the following error equation in the matrix form: where E = [1, 1, . . . , 1] T . Then according to the statement on stability in Theorem 2, we have This completes this proof.

Numerical examples
First, we give a fast algorithm for calculating the product of block Toeplitz matrices with vectors [21,43]. It is known that thematrices A y ⊗ A x , A y ⊗ B x , B y ⊗ A x , A y ⊗ C x , and C y ⊗ A x are block Toeplitz-Toeplitz block matrices. Taking the matrix A y ⊗ A x as an example, we first extend every N x × N x Toeplitz block matrix of A y ⊗ A x into a 2N x × 2N x cyclic matrix, so that the original matrix is expanded into an N y × N y block Toeplitzcyclic block matrix. Then we assemble the new block Toeplitz matrix into a 2N y × 2N y block cycle matrix, so that we generate a 2N y × 2N y block cycle matrix with each inner block of 2N x × 2N x block cycle, represented by C. Let F = F 2N x ⊗ F 2N y represent the twodimensional discrete Fourier transform matrix, and let c represent the first column of the newly assembled block cyclic-cyclic block matrix, so that we obtain the Fourier transform of the vector c: Next, we use the following property of the cyclic matrix: We can implement the fast matrix multiplication vector algorithm for two-dimensional problems, which reduces the computation amount of matrix vector multiplication from the traditional O(N 3 ) to O (N log N). We give several numerical examples to verify the validity of the finite volume element scheme. Just for simplicity of the rest of calculation, let h x = h y = h.
We define the L 2 -norm of the error between the exact solution U M and numerical solution u M as follows: Based on the above analysis, the L 2 -norm and convergence order of Example 1 calculated by MATLAB program are shown in the following tables. Tables 1 and 2 give the L 2 -norm and spatial convergence rate under the condition of τ = h and different β, γ , Tables 3 and 4 give the L 2 -norm and time convergence rate under the condition of h = 2 -8 and different β, γ . In agreement with the theoretical analysis given before, the time convergence rate and the space convergence rate are both of order 2. Tables 5 and 6 show the comparison of the calculation time of the three different calculation methods. We easily see that the fast bicgstab method greatly improves the calculation speed.
The L 2 -norm and convergence order of Example 2 calculated by MATLAB program are shown in Tables 7-12. As in Example 1, our final results are consistent with the theoretical analysis.
As in the previous two examples, we use the MATLAB program to produce the following data results (see Tables 13-14 and Fig. 5). We find that the calculated data are consistent with the theoretical results obtained from our previous analysis.

Conclusion
In this paper, we have successfully given the Crank-Nicolson finite volume element scheme for two-dimensional Riesz space-fractional convection-diffusion equations. We use the finite volume element scheme to discretize the space-fractional derivatives and the   CN-scheme to approximate the time derivatives. We show that the fully discrete scheme is stable and convergent. Finally, we verify the correctness and validity of the theoretical analysis through three examples.