High-order compact finite volume scheme for the 2D multi-term time fractional sub-diffusion equation

*Correspondence: ziwjiang@126.com 1School of Mathematics and Statistics, Shandong Normal University, Jinan, Shandong, China Abstract Based on an L1 interpolation operator, a new high-order compact finite volume scheme is derived for the 2D multi-term time fractional sub-diffusion equation. It is shown that the difference scheme is unconditionally convergent and stable in L∞-norm. The convergence order is O(τ 2–α + h1 + h2), where τ is the temporal step size and h1 is the spatial step size in one direction, h2 is the spatial step size in another direction. Two numerical examples are implemented, testifying to their efficiency and confirming their convergence order.


Introduction
The study of fractional partial differential equations has attracted many scholars' attention in recent decades [27][28][29]32]. Fractional differential equations are invoked in various scientific and engineering applications now. Physical phenomena in these fields, such as semiconductor, nuclear magnetic resonance, and signal processing, are successfully described by differential equations involving derivatives of fractional order [5,6,24,25]. Moreover, some physical phenomena cannot be described by single-term time fractional partial differential equations, and they have to be described by multi-term time fractional partial differential equations [4]. For example, when modeling various types of viscoelastic damping one must use multi-term time fractional partial differential equations. Because the fractional integrals and derivatives satisfy the nonlocal properties, fractional-order partial differential equations are different from classical partial differential equations.
Many theoretical analyses about the multi-term time fractional partial differential equations have been carried out, but the analytic solutions cannot be obtained exactly in most fractional partial differential equations [10,23]. Other scholars have discussed the numerical solution of fractional partial differential equations by using different methods, such as finite difference method, compact finite difference method, finite element method, alternating direction implicit schemes, and so on [1-3, 7-9, 11-13, 15-17, 19-22, 26, 30, 34-36, 38-41, 43-45, 47, 48, 50, 51, 53, 54]. Agarwal et al. discussed the solvability a linear loaded integro-differential equation in an infinite three-dimensional domain [1]. Then, Agarwal et al. analyzed special functions of differential equations [2]. Alderremy et al. structured new models of the multi-space fractional Gardner equation [3]. Wang et al. considered some tractable cases of limit values in which either a difference of two ingredients or a difference equation is used coupled with the relevant functional equations to give rise to unexpected results [41]. El-Sayed et al. introduced a numerical technique for solving a class of multi-term variable-order fractional differential equations [11]. Chen and Liu structured a finite difference method for two-dimensional anomalous sub-diffusion equation, and they analyzed the stability and convergence of the scheme [7]. Yuste and Acedo found a finite difference method which can solve the time fractional diffusion equation by using a forward Euler scheme, and they discussed the stability and convergence of the scheme [48]. Ren and Sun proposed some efficient and stable numerical methods for multi-term time fractional sub-diffusion equations, and they proved that these methods are stable convergent in the sense of the maximum norm [30]. Liu has made a great contribution to fractional partial differential equations by using a finite difference method and discussed the stability and convergence by using a new energy method [19,21]. Tadjeran et al. structured a second-order accurate numerical approximation for the fractional diffusion equation [36]. Sun et al. introduced a linearized second-order difference scheme for the nonlinear time fractional fourth-order reaction-diffusion equation [34]. Cui obtained finite difference schemes for the variable coefficients single-and multi-term time fractional diffusion equations [8]. To obtain higher order accuracy, a compact finite difference scheme is put forward. Wang et al. proposed compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation [44]. Cui structured a compact difference scheme for time fractional fourth-order equation with first Dirichlet boundary condition and discussed the stability and convergence [9]. Gao and Liu proposed a compact difference scheme for fourth-order temporal multi-term fractional wave equations and obtained maximum error estimates [13]. Later, Ji and Sun provided a high-order compact finite difference scheme for the fractional subdiffusion equation [16]. Zhuang et al. structured a new solution of the implicit numerical methods for the anomalous sub-diffusion equation and a numerical method of the modified anomalous sub-diffusion equation with a nonlinear source term [20,53]. Liu et al. provided finite element approximation for a modified anomalous sub-diffusion equation and a finite element method for two-dimensional multi-term time fractional diffusion wave equations [22,45]. Langlands et al. introduced an implicit solution method for the fractional diffusion equation and proved the accuracy and stability of the scheme [17]. Wang et al. proposed a high-order compact ADI method for two-dimensional fractional convection sub-diffusion equations with Neumann boundary conditions and proved its unconditional stability and global convergence to be O(τ 2-α + h 4 ) in the discrete L2 norm [43]. Zhang and Sun obtained alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation, and they proved that schemes are unconditionally stable and the numerical solution converges in the maximum norm [50]. Yang et al. provided a variational iteration method for solving multi-order fractional differential equations [47]. Furthermore, Zhuang et al. structured an implicit numerical method for the nonlinear fractional reaction sub-diffusion process, and they proved that the scheme was conditionally stable using the Fourier method [54]. Later, Zhang et al. proposed a Crank-Nicolson type difference scheme for the sub-diffusion equation which has a variable time and proved its unconditional stability in the maximum norm [51]. Tuan et al. analyzed a backward problem for fractional diffusion equation with Riemann-Liouville derivative [40]. Tuan et al. solved final value problem for nonlinear time fractional reaction-diffusion equation with discrete data [38]. Tuan et al. handled a terminal value problem for a generalization of fractional diffusion equation with hyper-Bessel operator [39]. Nguyen et al. settled a terminal value problem for time fractional diffusion equation by using regularization [26]. In addition, more researchers have found other numerical schemes such as the finite element method and others [12,15,18,33,35,37,46,52]. Zhou et al. structured a weak Galerkin finite element method for multi-term time fractional diffusion equations [52]. Lin et al. obtained a separable freconditioner for time-space fractional Caputo-Riesz diffusion equations [18]. Sheng and Shen proposed a space-time Petrov-Galerkin spectral method for time fractional diffusion equation [33]. Wu and Zhou provided parareal algorithms with local time integrators for time fractional differential equations [46]. Tang et al. introduced energy dissipation theory and numerical stability for time fractional phase field equations [37]. Obviously, the study of the time fractional sub-diffusion equation, which included only one item time fractional order, is thorough. In fractional physics, especially diffusion movement, the concept of Brown movement is extended because of the generalization of Gauss probability function. The scope of nuclear magnetic resonance is expanding by added resolving power, so one item time fractional order cannot explain this kind of problem. Multi-term time fractional sub-diffusion equations play an important role in actual production and life.
The finite volume method is also called the control volume method. This method is an important approximation method. We usually mesh the space and there are nonrepetitive control volumes near each grid point. Then, we integrate the equations separately on each control volume. Finally, we approximate the first-order partial derivatives with the function values of nodes. Zhang and Jiang provided a compact finite volume method for one-dimensional Schrödinger equation [49]. Wang et al. proposed a compact finite volume scheme for one-dimensional parabolic equations with third boundary conditions, and proved that the scheme is unconditionally stable [42]. The finite volume method has obvious advantages of integral conservation [31].
For the numerical approximation, take three positive integers M 1 , M 2 , N , and let h 1 = For every grid function u, we define the following notations: By using an L1 interpolation operator and a compact operator, we can construct a highorder compact finite scheme and establish the corresponding error estimate. The remainder of the article is arranged as follows. In Sect. 2, the compact finite volume scheme is derived. In Sect. 3, the existence and uniqueness, stability and convergence of the finite volume scheme are proved. In Sect. 4, two numerical examples are given to demonstrate the theoretical results. Finally, we obtain a brief conclusion.

The derivation of the compact finite volume scheme
To construct the compact finite volume scheme, we introduce some important lemmas. Especially, , we obtain that Especially, Proof We only prove the first integral approximation. The proof process of the second integral approximation is similar to the first one. We use Lagrange interpolation to approximate g(x) at the points ( . We can obtain that and Integrating g(x) by parts, we have To obtain the approximation order, we use substitution to simplify integral terms. For example, let ξ = x-ih 1 h 1 , so the integral terms can be simplified into and and Following the proof process of the first conclusion, we easily obtain the second conclusion.
Define a grid function space U h 1 = {g|g = (g 0 , g 1 , . . . , g M 1 )}. For every g ∈ U h 1 , define the integral operator as follows: Similarly, for every w ∈ U h 2 , define the integral operator as follows: Especially, Similarly, if w(y) ∈ C 5 [0, L 2 ] and y j+ 1 Especially, Proof We also only prove the first equation. Based on the Taylor expansion, we obtain the following equations: Similar to Lemma 1, to obtain the approximation order, we use substitution to simplify integral remainder. For example, let s = (i -1 2 )h + λ 2 h, so the integral remainder of (2) can be simplified into So, we can simplify the integral remainder of (2)-(5) by using the same way, then Similarly, we easily draw the second conclusion. ∈ U h 1 , define a compact operator as follows: .
Similarly, define a fractional index grid function space For every w j-1 2 ∈ U h 2 , define a compact operator as follows: .
Define an L1 approximation operator D α τ as follows: where the definition of a α k is the same as in Lemma 3.
Let us now construct a compact finite volume scheme for problem (1). On hτ , we now define the grid functions Suppose u(x, y, t) ∈ C 5,5,2 x,y,t ( × [0, T]). We consider equation (1) at the point (x, y, t n ), and we can have Integrating (6) on dual elements, we have We obtain that By using Lemma 2, we get that And then, we apply Lemma 1 to deal with spatial integral and use Lemma 3 to discretize the time fractional derivative. To optimize formula layout, we use the following notations: We obtain where r n ij = C 1 (τ 2-α ) + C 2 (h 4 1 + h 4 2 ). Notice that u(x, y, t) = ψ(x, y, t) and u(x, y, 0) = ϕ(x, y), so we have U n 0,j = ψ(0, y j , t n ), U n L 1 ,j = ψ(L 1 , y j , t n ), So, we leave out the infinitesimal and replace exact solution U n i,j with numerical solution u n i,j , then we obtain that u n 0,j = ψ(0, y j , t n ), u n L 1 ,j = ψ(L 1 , y j , t n ), These above equations are the compact finite volume scheme of question (1).

Analysis of the compact finite volume scheme
In this section, we prove the existence and uniqueness, stability and convergence of the compact finite volume scheme. First, we introduce the norms in the space U h . Let g = (g 00 , g 01 , . . . , Some important lemmas, which can be used in the following verifications, will be given.

Theorem 1 The solution of the compact finite volume scheme (13)-(15) is existent and unique.
Proof Let u n = (u n 00 , u n 01 , . . . , u n M 1 M 2 ) ∈ U h , and the numerical solution of u 0 can be obtained by (14). If the numerical solutions of u 0 , u 1 , . . . , u n-1 are existing and unique, we can obtain nonhomogeneous linear equations about u n from (13) and (15). If homogeneous linear equations only have a zero solution can be proved, we can prove that the existence and uniqueness of u n can be proved. Define S α = τ α (2α), S β = τ β (2β).
What we should do next is to prove that (16) and (17) only have a zero solution. Equation (16) can be rewritten into By using Lemma 5, equation mentioned above also can be transformed into If the equation holds on, we can obtain that u n ∞ = 0. Furthermore, u n = 0. According to the inductive principle, equations (13)-(15) have a unique solution.

Theorem 3
The solution of the compact finite volume scheme (13)-(15) is convergent.
Denote the maximum norm errors Tables 1 and 2 show the maximum norm errors and the convergence order of the finite volume scheme. In order to simplify the calculation, we let h = h 1 = h 2 . Suppose If τ is small enough, then E ∞ (h, τ ) ≈ O(h p ). Consequently,   So, p is the convergence order with respect to the spatial step size. Similarly, we can obtain If under different meshes Mass h keep about the same size, we can explain the integral conservation of the compact finite volume scheme.
We compare the exact solution and the numerical solution under different meshes. In Tables 1 and 2, temporal and spatial convergence orders are shown. We fix sufficiently small spatial step size h = 1 1000 and vary the temporal step sizes. Table 1 lists the numerical results for different temporal step sizes. Similarly, we fix sufficiently small temporal step  Figure 3 The red line represents temporal error order, the blue line represents spatial error order sizes τ = 1 1000 and vary the spatial step sizes. Also, Table 2 lists the numerical results for different spatial step sizes.
In order to observe error orders more intuitively, we plot Fig. 3 about error orders, where slope represents the error order. We can observe that the temporal error order is about 2α and the spatial error order is about h 4 . It is also not difficult to obtain the corresponding source term f (x, y, t) and the initial and boundary conditions with α = 0.8, β = 0.3, which are respectively f (x, y, t) = (4) (4α) xx 2 yy 2 t 3-α + (4) (4β) xx 2 yy 2 t 3-β + 2t 3 yy 2 + 2t 3 xx 2 and ψ(x, y, t) = 0, ϕ(x, y) = 0.
Like the previous example, we compute the maximum norm errors of the numerical solution to test the convergence of the compact finite volume scheme, which is defined the same as in Example 1. Denote In Table 3, we compute the problem for a longer time by fixing N = 8, 16, 32, 64, 128, and still choosing h = 1 1000 . And then, we compute the problem for a longer space by fixing M = 8, 16, 32, 64, 128, and still choosing τ = 1 1000 in Table 4. As in Example 1, we plot Fig. 6 about error orders, where slope represents the error order. We can observe that the temporal error order is about 2α and the spatial error order is about h 4 .     Figure 6 The red line represents temporal error order, the blue line represents spatial error order

Conclusion
In this article, we constructed a compact finite volume scheme for the multi-term time fractional sub-diffusion equation. We proved the existence and uniqueness, stability and L ∞ convergence of our scheme. Two numerical results show that the scheme is conserved and convergent with the order O(τ 2-α + h 4 ).