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

The study of fractional partial differential equations has attracted many scholars’ attention in recent decades [27–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 sub-diffusion 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(\tau ^{2-\alpha }+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].

Consider the following two-dimensional multi-term time fractional sub-diffusion equation with homogeneous source term on the computational domain \(\Omega \times [0,T]\):

where \(\Omega =(0,L_{1})\times (0,L_{2})\), ∂Ω is the boundary, \(0<\beta <\alpha <1\), and the operator \({}^{C}_{0}D^{\alpha }_{t}\) denotes the Caputo fractional derivative of order α defined by

\(\Gamma (\cdot )\) is the gamma function, and \(\psi (x,y,t)\), \(\varphi (x,y)\), and \(f(x,y,t)\) are the known smooth functions, \(\psi (x,y,0)=\varphi (x,y)\).

For the numerical approximation, take three positive integers \(M_{1}\), \(M_{2}\), N, and let \(h_{1}=\frac{L_{1}}{M_{1}}\), \(h_{2}=\frac{L_{2}}{M_{2}}\), \(\tau =\frac{T}{N}\). Define \(x_{i}=ih_{1}\) (\(0\leq i\leq M_{1}\)), \(y_{j}=jh_{2}\) (\(0\leq i\leq M_{2}\)), \(t_{n}=n\tau\) (\(0\leq n\leq N\)), \(\Omega _{h_{1}}=\{x_{i}\mid 0\leq i\leq M_{1}\}\), \(\Omega _{h_{2}}=\{y_{j}\mid 0\leq i\leq M_{2}\}\), \(\Omega _{\tau }=\{t_{n}\mid 0\leq i\leq N\}\). Define \(\overline{\Omega }_{h}=\Omega _{h_{1}}\times \Omega _{h_{2}}\), \(\Omega _{h}=\overline{\Omega }_{h}\cup \partial \Omega \), \(\Omega _{h\tau }=\Omega _{h}\times \Omega _{\tau }\), then the computational domain \(\Omega \times [0,T]\) is covered by \(\Omega _{h\tau }\). Define primal partition \(I_{h}\), grid element \(I_{ij}=[x_{i},x_{i+1}]\times [y_{j},y_{j+1}]\) (\(0\leq i\leq M_{1}-1\), \(0 \leq j\leq M_{2}-1\)). Define dual partition \(I_{h}^{\ast }\), \(I_{0j}^{\ast }=[0,x_{\frac{1}{2}}]\times [y_{j-\frac{1}{2}},y_{j+ \frac{1}{2}}]\), \(I_{0i}^{\ast }=[0,y_{\frac{1}{2}}]\times [x_{i-\frac{1}{2}},x_{i+ \frac{1}{2}}]\), \(I_{ijp}^{\ast }=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times [y_{j- \frac{1}{2}},y_{j+\frac{1}{2}}]\) (\(1\leq i\leq M_{1}-1\), \(1\leq j\leq M_{2}-1\)), \(I_{jM_{1}}^{\ast }=[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]\times [x_{M_{1}- \frac{1}{2}},L_{1}]\), \(I_{iM_{2}}^{\ast }=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times [x_{M_{2}- \frac{1}{2}},L_{2}]\). Suppose that \(u=\{u_{i,j}^{n}\mid 0\leq i\leq M_{1},0\leq j\leq M_{2},0\leq n\leq N \}\) is a grid function on \(\Omega _{h\tau }\). For every grid function u, we define the following notations:

By using an L1 interpolation operator and a compact operator, we can construct a high-order 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.

To construct the compact finite volume scheme, we introduce some important lemmas.

Lemma 1

If \(g(x)\in C^{3}[0,L_{1}]\) and \(x_{i+\frac{1}{2}}=(i+\frac{1}{2})h_{1}\), \(0\leq i\leq M_{1}-1\), then

where

Especially,

and

Similarly, if \(w(y)\in C^{3}[0,L_{2}]\), we obtain that

where

Especially,

and

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 \((x_{i-1},g(x_{i-1}))\), \((x_{i},g(x_{i}))\), and \((x_{i+1},g(x_{i+1}))\). We can obtain that

where

and

Integrating \(g(x)\) by parts, we have

To obtain the approximation order, we use substitution to simplify integral terms. For example, let \(\xi =\frac{x-ih_{1}}{h_{1}}\), so the integral terms can be simplified into

and

Then

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}, \ldots , g_{M_{1}})\}\). For every \(g\in U_{h_{1}}\), define the integral operator as follows:

with \((S_{x}g)_{0}=\frac{h_{1}}{3}g_{0}+\frac{5h_{1}}{24}g_{1}- \frac{h_{1}}{24}g_{2}\), \((S_{x}g)_{M_{1}}=-\frac{h_{1}}{24}g_{M_{1}-2}+ \frac{5h_{1}}{24}g_{M_{1}-1}+\frac{h_{1}}{3}g_{M_{1}}\).

Similarly, for every \(w\in U_{h_{2}}\), define the integral operator as follows:

with \((S_{y}w)_{0}=\frac{h_{2}}{3}w_{0}+\frac{5h_{2}}{24}w_{1}- \frac{h_{2}}{24}w_{2}\), \((S_{y}w)_{M_{2}}=-\frac{h_{2}}{24}w_{M_{2}-2}+ \frac{5h_{2}}{24}w_{M_{2}-1}+\frac{h_{2}}{3}w_{M_{2}}\).

Lemma 2

If \(g(x)\in C^{5}[0,L_{1}]\) and \(x_{i+\frac{1}{2}}=(i+\frac{1}{2})h_{1}\), \(0\leq i\leq M_{1}-1\), then

where

Especially,

and

Similarly, if \(w(y)\in C^{5}[0,L_{2}]\) and \(y_{j+\frac{1}{2}}=(j+\frac{1}{2})h_{2}\), \(0\leq i\leq M_{2}-1\), then

where

Especially,

and

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-\frac{1}{2})h+\frac{\lambda }{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. □

Define a fractional index grid function space \(\overline{U}_{h_{1}}=\{g|g=(g_{\frac{1}{2}},g_{\frac{3}{2}}, \ldots , g_{M_{1}-\frac{1}{2}})\}\). For every \(g_{i-\frac{1}{2}}\in \overline{U}_{h_{1}}\), define a compact operator as follows:

with

Similarly, define a fractional index grid function space \(\overline{U}_{h_{2}}=\{w|w=(w_{\frac{1}{2}},w_{\frac{3}{2}}, \ldots , w_{M_{2}-\frac{1}{2}})\}\). For every \(w_{j-\frac{1}{2}}\in \overline{U}_{h_{2}}\), define a compact operator as follows:

with \((A_{y}w)_{\frac{1}{2}}=3w(y_{\frac{1}{2}})-5w(y_{\frac{3}{2}})+4w(y_{ \frac{5}{2}})-w(y_{\frac{7}{2}})\), \((A_{y}w)_{M_{2}-\frac{1}{2}}=-w(y_{M_{2}- \frac{7}{2}})+4w(y_{M_{2}-\frac{5}{2}})-5w(y_{M_{2}-\frac{3}{2}})+3w(y_{M_{2}- \frac{1}{2}})\).

Lemma 3

([14])

If \(0<\alpha <1\), \(g\in C^{2}[0,t_{n}]\), then

and

where \(a_{k}^{\alpha }=(k+1)^{1-\alpha }-k^{1-\alpha }\).

Define an L1 approximation operator \(D_{\tau }^{\alpha }\) as follows:

where the definition of \(a^{\alpha }_{k}\) is the same as in Lemma 3.

Let us now construct a compact finite volume scheme for problem (1). On \(\Omega _{h\tau }\), we now define the grid functions

Suppose \(u(x,y,t)\in C^{5,5,2}_{x,y,t}(\Omega \times [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_{ij}^{n}=C_{1}(\tau ^{2-\alpha })+C_{2}(h_{1}^{4}+h_{2}^{4})\).

Notice that \(u(x,y,t)=\psi (x,y,t)\) and \(u(x,y,0)=\varphi (x,y)\), so we have

So, we leave out the infinitesimal and replace exact solution \(U_{i,j}^{n}\) with numerical solution \(u_{i,j}^{n}\), then we obtain that

These above equations are the compact finite volume scheme of question (1).

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}, \ldots , g_{M_{1}M_{2}})\in U_{h}\), where \(U_{h}=U_{h_{1}}\times U_{h_{2}}\). Denote

Some important lemmas, which can be used in the following verifications, will be given.

Lemma 4

([14])

Suppose \(\alpha \in (0,1)\), and the definition of \(a_{l}^{\alpha }\) is the same as in Lemma (3), \(l=0,1,\ldots \) , then

1. (I)

   \(1=a_{0}^{\alpha }>a_{1}^{\alpha }>a_{2}^{ \alpha }>\cdots >a_{l}^{\alpha }; a_{l}^{\alpha }\rightarrow 0\), when \(l\rightarrow \infty \);

2. (II)

   \((1-\alpha )l^{-\alpha }< a_{l-1}^{\alpha }<(1- \alpha )(l-1)^{-\alpha }\), \(l\geq 1\).

Lemma 5

Let \(p=(p_{00}, p_{01}, \ldots , p_{M_{1}M_{2}})\in U_{h}\), \(q=(q_{00}, q_{01}, \ldots ,q_{M_{1}M_{2}})\in U_{h}\), then

Theorem 1

The solution of the compact finite volume scheme (13)–(15) is existent and unique.

Proof

Let \(u^{n}=(u_{00}^{n}, u_{01}^{n}, \ldots , u_{M_{1}M_{2}}^{n})\in U_{h}\), and the numerical solution of \(u^{0}\) can be obtained by (14). If the numerical solutions of \(u^{0}, u^{1}, \ldots , 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_{\alpha }=\tau ^{\alpha }\Gamma (2-\alpha )\), \(S_{\beta }=\tau ^{\beta } \Gamma (2-\beta )\).

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 \(\Vert u^{n}\Vert _{\infty }=0\). Furthermore, \(u^{n}=0\).

According to the inductive principle, equations (13)–(15) have a unique solution. □

Theorem 2

Suppose that the compact finite volume scheme (13)–(15) has a solution. We record it as \(\{v^{n}_{ij}\vert 0\leq i\leq M_{1}, 0\leq j\leq M_{2}, 0 \leq n\leq N\}\), then

where \(c=\frac{1}{2}\max \{T^{\alpha }\Gamma (1-\alpha ), T^{\beta }\Gamma (1- \beta )\}\), \(\Vert f^{m}\Vert _{\infty }=\max_{1\leq i\leq M_{1}-1,1\leq j \leq M_{2}-1}\vert f^{m}_{ij}\vert \).

Proof

To optimize formula layout, we use the following notations:

Equation (13) can be written into

By using Lemma 5, we can obtain

By using Lemma 4, we can obtain

Bring (20) and (21) into (19), then

where \(c=\frac{1}{2}\max \{T^{\alpha }\Gamma (1-\alpha ), T^{\beta }\Gamma (1- \beta )\}\).

We prove the establishment of (18) now.

According to (22), we can obtain that when \(n=1\),

It also can be written as

So, when \(k=1\), equation (18) is valid.

Suppose that equation (22) also satisfies conditions when \(k=1, 2, \ldots , n-1\). According to (22), we can obtain

It is equivalent to

So, when \(k=n\), equation (18) satisfies conditions.

According to the inductive principle, Theorem 2 holds. □

Theorem 3

The solution of the compact finite volume scheme (13)–(15) is convergent. Suppose that \(\{U^{n}_{ij}\vert 0\leq i\leq M_{1}, 0\leq j\leq M_{2}, 0 \leq n\leq N\}\) is the solution of equation (1) and \(\{u^{n}_{ij}\vert 0\leq i\leq M_{1}, 0\leq j\leq M_{2}, 0 \leq n\leq N\}\) is the solution of the compact finite volume scheme (13)–(15). Define

then

Proof

By subtracting equations (13)–(15) from the compact finite volume scheme (10)–(12), we can obtain error equations:

By using Theorem 2, we can obtain

where \(c=\frac{1}{2}\max \{T^{\alpha }\Gamma (1-\alpha ), T^{\beta }\Gamma (1- \beta )\}\).

Furthermore, it is easy to obtain that

So, the result

is proved. □

We now present two numerical examples using the numerical methods as discussed above, where the errors involved are measured by comparing the numerical solutions with exact solutions. To simplify the calculation, let \(h=h_{1}=h_{2}\).

Example 1

Let \(\Omega =[0,1]\times [0.1]\), \(T=1\). In order to test the convergence of the finite volume scheme, we refer to the exact solution of problem (1) as

We can obtain the corresponding source term \(f(x,y,t)\) and the initial and boundary conditions with \(\alpha =0.4\), \(\beta =0.2\), which are respectively

and

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_{\infty }(h,\tau )\approx O(h^{p})\). Consequently,

So, p is the convergence order with respect to the spatial step size. Similarly, we can obtain

for small enough h. Denote

Define \(\mathit{Mass}_{h}=\sum_{i=0}^{M_{1}}\sum_{j=0}^{M_{2}}u_{ij}^{N}h^{2}\). If under different meshes \(\mathit{Mass}_{h}\) keep about the same size, we can explain the integral conservation of the compact finite volume scheme.

The effect of numerical solution and exact solution at fixed \(h=\frac{1}{1000}\) and \(\tau =\frac{1}{64}\)

Full size image

The effect of numerical solution and exact solution at fixed \(\tau =\frac{1}{1000}\) and \(h=\frac{1}{64}\)

Full size image

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=\frac{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 sizes \(\tau =\frac{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-\alpha \) and the spatial error order is about \(h^{4}\).

The red line represents temporal error order, the blue line represents spatial error order

Full size image

The effect of numerical solution and exact solution at fixed \(h=\frac{1}{1000}\) and \(\tau =\frac{1}{64}\)

Full size image

The effect of numerical solution and exact solution at fixed \(h=\frac{1}{1000}\) and \(\tau =\frac{1}{64}\)

Full size image

Example 2

Let \(\Omega =[0,1]\times [0.1]\), \(T=1\). As before, we refer to the exact solution of problem (1) as

It is also not difficult to obtain the corresponding source term \(f(x,y,t)\) and the initial and boundary conditions with \(\alpha =0.8\), \(\beta =0.3\), which are respectively

and

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=\frac{1}{1000}\). And then, we compute the problem for a longer space by fixing \(M=8, 16, 32, 64, 128\), and still choosing \(\tau =\frac{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-\alpha \) and the spatial error order is about \(h^{4}\).

The red line represents temporal error order, the blue line represents spatial error order

Full size image

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_{\infty }\) convergence of our scheme. Two numerical results show that the scheme is conserved and convergent with the order \(O(\tau ^{2-\alpha }+h^{4})\).