A two-grid ADI finite element approximation for a nonlinear distributed-order fractional sub-diffusion equation

1.   Introduction

Distributed order differential equations ^{[4,10,11,12,13,14,15,16,38,39,40]} can be seen as a natural extension of single-term and multi-term fractional order differential equations. When the fractional order derivative term is sufficiently large, we take its limit state to obtain the distributed order derivative. Distributed order models are a broader class of models with a broader meaning. It can be used to describe processes that cannot be portrayed by a single-term or multi-term fractional differential equations, such as retarding sub-diffusion and accelerating superdiffusion ^[5,30]. The distribution order differential equations can be classified into time, space and space-time distributed order differential equations according to the location of the distributed order integral terms. Time distributed order differential equations are often used to describe some complex processes in which the diffusion index varies with time. It is now playing an important role in many fields and has become a popular research topic in the international academic community. However, the complexity and nonlocality of the distributed-order operator make it difficult to solve the exact solution of the distributed-order differential equations, so scholars have turned to its numerical solution and made important progress ^[13]. Among many algorithms, the finite element method is favored by scholars because of its strong regional adaptability, more flexible mesh parting, lower smoothness requirement, and strong generality ^[16].

In this paper, we consider the following nonlinear distributed-order fractional sub-diffusion equation

where $\Omega = I_x\times I_y = (0, L)\times(0, L)$, and the boundary $\partial\Omega$ is Lipschitz continuous. $J = (0, T]$ is the time interval, and the nonlinear reaction term $m(u)$ satisfies $|m(u)|\le C|u|$ with $|m^{'}(u)|\le C$, where $C$ is a positive constant. The term $f(x, y, t)$ is a given source function.

Define

where

and $\omega(\alpha)\geq0, \int_0^1\omega(\alpha)d\alpha = C_0 > 0$.

Inspired by the works ^{[19,23,24,34]}, in this article, we propose a two-grid ADI FE scheme with the WSGD approximation formula to solve nonlinear distributed-order fractional sub-diffusion equation. In what follows, we will introduce the advancements of the WSGD approximation formula, two-grid method, and ADI FE method.

In 2015, Tian et al. ^[31] proposed a WSGD approximation formula for the Riemann-Liouville space fractional derivative. Based on this operator, they got the second-order convergence, which is independent of the changed fractional parameters. The WSGD operator, with its outstanding advantages such as high-order approximation, has attracted the attention of many scholars and has been widely used, and many fruitful research results have been achieved. Wang and Vong in ^[32] discussed the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation by using the WSGD formula to approximate the Caputo fractional derivative of time. In 2016, Liu et al. ^[23] presented a two-grid FE scheme with the WSGD operator for solving the nonlinear fractional Cable equation. In 2017, Liu et al. ^[25] solved a Caputo time-fractional sub-diffusion equation by using a high-order local discontinuous Galerkin (LDG) method combined with the WSGD approximation. Wang et al. ^[33] discussed an $H^1$-Galerkin mixed finite element (MFE) method combined with the WSGD operator for solving nonlinear convection-diffusion equation with time fractional derivative. In 2020, Saffarian and Mohebbi ^[29] discussed the application of ADI Galerkin spectral element method with the WSGD operator in solving two-dimensional time fractional sub-diffusion equation.

In 1994, Xu ^[34] presented a new finite element discretization technique based on two (coarse and fine) subspaces for a semilinear elliptic boundary value problem. And then in 1996, Xu ^[35] did further study and explored about the two-grid method. This method has attracted many scholars in solving partial differential equations because of the advantage of saving computational time. Since its introduction, the method has been mainly applied to numerically solve integer-order PDEs by many computational scholars ^[21,27]. Until 2015, Liu et al. ^[26] presented the application of the two-grid finite element method to numerically solve the nonlinear fourth-order fractional differential equations, in which the fractional derivative is the Caputo type. And in 2016, Liu et al. ^[23] discussed the numerical solution of nonlinear fractional Cable equation with initial and boundary condition by using two-grid FE method with higher-order time approximate scheme. Chen et al. in ^[1] proposed a fully discrete two-grid modified method of characteristics (MMOC) scheme for solving nonlinear variable-order time-fractional advection-diffusion equations in two space dimensions. Zeng and Tan in ^[37] developed a two-grid FE methods for variable coefficient time fractional diffusion equations. However, there is still a lack of work on the numerical solution of fractional partial differential equations using the two-grid FE method, especially there are few studies on the solution of nonlinear time-distributed partial differential equations using the two-grid FE method, so the research on this method still needs to be focused. The ADI FE method is an efficient algorithm for solving multi-dimensional differential equations. It inherits the advantages of the ADI method of low computation and low storage, and also has the characteristics of high accuracy of the FE algorithm. Therefore, it has received a lot of attention from researchers in solving differential equations. In 1971, Douglas and Dupont ^[8] for the first time proposed the ADI FE method. Then Dendy ^[6,7] and Fernandes ^[9] and Zhang ^[36] have studied the method in depth and further extended the application of the method. In 2013, Li and Xu ^[19] discussed the Galerkin FE method in space and the ADI method in the time stepping for the two-dimensional fractional diffusion-wave equation. Li and Xu ^[18] in 2014 studied the two-dimensional time fractional evolution equation by ADI Galerkin FE method. In 2017, Chen and Li ^[3] proposed an efficient ADI Galerkin method for solving a time-fractional partial differential equation with damping. In the same year, Li and Huang ^[20] solved a class of 2D nonlinear fractional diffusion-wave equations with the Caputo-type temporal derivative and Riesz-type spatial derivative by using the ADI FE method. In 2020, Chen ^[2] used the ADI FE method to numerically solve two classes of Riesz space fractional partial differential equations. In the same year, Qiu et al. ^[28] presented the ADI Galerkin FE method for solving the distributed-order time-fractional mobile-immobile equation in two dimensions.

However, there is little research on solving distributed-order partial differential equations using the ADI FE method. In particular, we are not aware of any studies that apply the two-grid ADI FE method to numerically solve the time-distributed order reaction diffusion equation with nonlinear terms and time-integer order derivatives. Based on these considerations, we propose two-grid ADI FE method for solving this model (1.1), prove the stability, and derive the a priori error estimates. Finally, we verify the effectiveness and computational efficiency of the algorithm by carrying out numerical tests.

This paper is organized as follows. In Section 2, we give some preliminaries and lemmas. In Section 3, we present the ADI FE numerical approximation of the nonlinear distributed order equation. In Section 4, we analyze the stability and convergence of the fully discrete scheme. In Section 5, some numerical examples are presented to confirm the theoretical analysis. In Section 6, conclusions and future works are discussed.

2.   Preliminaries and lemmas

We insert the nodes $\alpha_k = k\Delta\alpha$, $k = 0, 1, 2\cdots, 2K$ in the interval $[0, 1]$, where $\Delta\alpha = \frac{1}{2K}$ and $0 = \alpha_0 < \alpha_1 < \alpha_2 < \cdots < \alpha_{2K} = 1$. $\Delta t = \frac{T}{N+1}$ is the step size, $t_n = n\Delta t$ $(n = 0, 1, 2, \cdots, N+1)$, and $0 = t_0 < t_1 < t_2 < \cdots < t_{N+1} = T$, where $N$ is positive integer, and for a smooth function $\psi$ on $[0, T]$, we denote $\psi^n = \psi(t_n)$, and $\delta_t\psi^{n+\frac{1}{2}} = \frac{\psi^{n+1}-\psi^{n}}{\Delta t}$, $\psi^{n+\frac{1}{2}} = \frac{\psi^{n+1}+\psi^{n}}{2}$.

Define $H^m(\Omega)$, $L^{\infty}(\Omega)$, $L^{2}(\Omega)$ and $\|\cdot\|_m$, $\|\cdot\|_\infty$, $\|\cdot\|$ be the usual Sobolev spaces and their norms, respectively. Meanwhile, the inner product of $L^2(\Omega)$ is denoted by $(\cdot, \cdot)$. For $\hbar > 0$ ($\hbar$ indicates fine grid size $h$ or coarse grid size $\widetilde{H}$), $r\geq2$, define $V_{\hbar}^r\subset H_0^1(\Omega)$ is a finite-dimensional subspace and $V_{\hbar}^r$ satisfies the following properties ^[9]

where

In order to give the numerical approximation of Eq (1.1), we use the following composite Trapezoid formula.

Lemma 2.1. Letting $s(\alpha)\in C^{2}[0, 1]$, we have

where

Lemma 2.2. ^[31,32] For $0 < \alpha < 1$, the following approximate formula holds

where

and

Lemma 2.3. ^[23] For series $\{\gamma_j^\alpha\}$ given by Lemma 2.2, we have

Lemma 2.4. ^[23,31,32] For series $\{q_\alpha(j)\}$ defined by formula (2.4). Then for any integer $n$ and any positive integer $N$ as well as real vector $(u^0, u^1, \ldots u^N)\in R^{N+1}$, we get

3.   Derivation of the two-grid ADI FE scheme

By using Lemma 2.1, we discretize the integral term of the distributed order equation. Suppose $s(\alpha)\in C^2[0, 1]$, and set $s(\alpha) = \omega(\alpha)_0^CD_t^\alpha u$, to arrive at

where $R_1 = \text{O}(\Delta\alpha^2)$.

By using Lemma 2.2 and formula (3.1), the time discretization scheme of Eq (1.1) can be written as

Denote

So we have

Finding $u_h^{n+1}\in V_h^r$, we arrive at the following FE scheme of formula (3.4)

To speed up the calculation, we create the following two-grid ADI finite element scheme.

Denote

Then, we get the following two-grid ADI finite element scheme.

Step 1: Letting $U_{\widetilde{H}}^{n+1}:[0, T]\longmapsto V_{\widetilde{H}}^r\subset V_h^r$ be the solution of the following nonlinear system which is based on the coarse grid $\mathcal{T}_{\widetilde{H}}$, we have

Step 2: Letting $\mathfrak{u}_h^{n+1}:[0, T]\longmapsto V_h^r$ be the solution of the following nonlinear system which is based on the fine grid $\mathcal{T}_h$, we have

where $h\ll {\widetilde{H}}.$

Next, we further discuss the numerical scheme (3.6). We rewrite it into a more familiar alternating direction finite element form. Assume that $V_{h}^r = V_{h, x}^r\otimes V_{h, y}^r$, where $V_{h, x}^r$ and $V_{h, y}^r$ are finite-dimensional subspaces of $H_0^1(\Omega)$. Let $\{\varphi_i\}_{i = 1}^{N_x-1}$ and $\{\chi_p\}_{p = 1}^{N_y-1}$ be bases of $V_{h, x}^r$ and $V_{h, y}^r$, respectively. So $\{\varphi_i\chi_p\}_{i = 1, p = 1}^{N_x-1, N_y-1}$ is the tensor product basis of $V_{h}^r$.

Let

where

Let $v_{\widetilde{H}} = \varphi_l\chi_m, l = 1, \cdots, N_x-1; m = 1, \cdots, N_y-1$, then the scheme (3.6) can be changed into

where

Denote

and let

So we obtain the matrix form of the ADI Galerkin scheme (3.10) as follows

where $\otimes$ denotes the matrix tensor product and $I_{N_x-1}$ and $I_{N_y-1}$ denote the identity matrices of order $N_x-1$ and $N_y-1$, respectively. By introducing the auxiliary vector $\hat{\beta}^{(n+1)}$, then the Eq (3.12) is equivalent to

Thus we determine $\beta^{(n+1)}$ by solving two sets of independent one-dimensional problems.

In $x$-direction, we calculate $\hat{\beta}_p^{(n+1)}$ by using the following equations

where

In $y$-direction, we can calculate $\beta_i^{(n+1)}$ by using the following equations

where

For the numerical scheme (3.7) we use a similar approach to the above process. In the next process, we consider the stability and convergence of the two-grid FE systems (3.6) and (3.7).

4.   Analysis of stability and convergence

4.1.   Stability

Theorem 4.1. For the systems (3.6) and (3.7), which is based on coarse grid $\mathcal{T}_{\widetilde{H}}$ and fine grid $\mathcal{T}_h$, the following stable inequality holds

Proof. Taking $v_h = 2\mathfrak{u}_h^{n+\frac12}$ in scheme (3.7), we obtain

By using Cauchy-Schwarz inequality and Young inequality, we can get

Multiply $\Delta t$ on both sides of inequality (4.3) and sum the resulting inequality (4.3) for $n$ from $0$ to $N$ to obtain

By using Lemma 2.4 and Gronwall lemma, we have

Next we discuss the term $\|U_{\widetilde{H}}^{n}\|^2$. Take $v_{\widetilde{H}} = 2U_{\widetilde{H}}^{n+\frac12}$ in scheme (3.6) and use the same process of the derivation for $\|\mathfrak{u}_h^N\|^2$ to arrive at

Use $\Delta t\sum\limits_{n = 0}^N\le T$ and substitute inequality (4.6) into inequality (4.5) to derive

Thus, the proof of stability is completed.

4.2.   Convergence

Define a Ritz projection operator $\Re_\hbar:H_0^1(\Omega)\rightarrow V_\hbar^r$ by

Then, we introduce the property of the projector $\Re_\hbar$.

Lemma 4.2. ^[9] If $\frac{\partial^lu}{\partial t^l}\in L^p(H^r)$, $l = 0, 1, 2$, $p = 2, \infty$, then there exists a constant $C$ that is independent of $\hbar$, such that

where, $k = 0, 1$, $1\le s\le r$ and $\hbar$ is coarse grid step length $\widetilde{H}$ or fine grid size $h$.

Lemma 4.3. ^[7] Let $D$ represent the operator $\frac{\partial}{\partial t}$ or $\frac{\partial^2}{\partial t^2}$. By using the triangle inequality and inequality (2.1), we have

To simplify the notations, we denote

Theorem 4.4. Let $u(t_n)$, $U_{\widetilde{H}}$ and $\mathfrak{u}_h$ be the solution of Eq (1.1), scheme (3.6) and scheme (3.7), respectively. Assume that $u(t_n)\in L^\infty(H^r)$, $\frac{\partial u}{\partial t}\in L^2((0, T];H^r)$, $\frac{\partial^3 u}{ \partial x\partial y\partial t}\in L^2((0, T];L^2)$ and $r\ge 2$, then we obtain the following error results

where $C$ is a positive constant independent of fine grid step length $h$, coarse grid step length ${\widetilde{H}}$ as well as time step parameter $\Delta t$.

Proof. Combining scheme (3.4) and scheme (3.7), the error equation is as follows

Taking $v_h = 2\eta_u^{n+\frac12}$, summing up Eq (4.11) for $n$ from $0$ to $N$, and multiplying $\Delta t$ on both sides of the above equation, we can get

In the following, we estimate the terms $E_i$, $i = 1, 2, \ldots, 8$.

From Lemma 4.2, we can get

Then we have

Use Cauchy-Schwarz inequality, Young inequality and Lemma 4.2 to arrive at

By using the Cauchy-Schwarz inequality and Young inequality, we have

By using the Young inequality and H${\rm\ddot{o}}$lder inequality, we have

From Lemma 4.2 and Lemma 4.3, we obtain

So we can get

Similarly, we have

Together, we have

Substituting inequalities (4.14)–(4.20) into Eq (4.12) and using Gronwall lemma, we have

For the next discussion, we give the estimate for the term $\Delta t\sum\limits_{n = 0}^N\|(u^{n+\frac12}-U_{\widetilde{H}}^{n+\frac12})^2\|^2$.

Subtract scheme (3.6) from scheme (3.4) to arrive at

Take $v_{\widetilde{H}} = 2\rho_u^{n+\frac12}$ in scheme (4.22) and use the similar process of proof to the estimate for $\| u(t_n)-\mathfrak{u}_h^{n}\|$ to obtain

By using the similar progress of ^[22], we have

By using the similar progress of ^[17], we can get

Substitute inequalities (4.24) and (4.25) into inequality (4.21) to have

By using triangle inequality, we finish the proof of the theorem.

5.   Numerical experiment

In this section, we carry out numerical experiments to illustrate our theoretical results.

Example 5.1. On the space-time domain $[0, 1]^2\times[0, \frac12]$, we choose $\omega(\alpha) = \Gamma(4-\alpha)$, nonlinear term $m(u) = \sin(u)$, source term

and the exact solution

In Table 1, with $\Delta\alpha = 1/500$, $\Delta t = 1/200$, ${\widetilde{H}}_x = {\widetilde{H}}_y = 1/2, 1/3, \, 1/4, \, 1/5$, $1/6$ and $h_x = h_y = 1/4, \, 1/9, \, 1/16, \, 1/25$, $1/36$, the error estimate result, second-order spatial convergence rates and computation time of $u$ are obtained. In Table 2, with $\Delta\alpha = \Delta t = h_x = h_y = {\widetilde{H}}_x^2 = {\widetilde{H}}_y^2 = 1/4, \, 1/16, \, 1/36, \, 1/64$ and $1/100$, we get the convergence in time and space. The data results show that the two-grid ADI finite element metnod can effectively solve the nonlinear time distributed-order reaction-diffusion equations.

In Figure 1, we give the surface for the exact solution $u$ at $t = 0.5$ with $h_x = h_y = 1/200$. Furthermore, in Figures 2–7, by taking $\Delta\alpha = 1/500$, $\Delta t = 1/200$, $h_x = h_y = {\widetilde{H}}_x^2 = {\widetilde{H}}_y^2 = 1/4, \, 1/9, \, 1/16, \, 1/25, \, 1/36$ and $1/49$, at $t = 0.5$, we show the surfaces for the numerical solutions $u_h$. In order to show the error behavior between the numerical solution and the exact solution, in Figures 8–13, we give the surfaces for the errors $u-u_h$ with $\Delta\alpha = 1/500$, $\Delta t = 1/200$, $h_x = h_y = {\widetilde{H}}_x^2 = {\widetilde{H}}_y^2 = 1/4, \, 1/9, \, 1/16, \, 1/25, \, 1/36$, $1/49$ at $t = 0.5$. It can be seen from the image display that the numerical method is effective in solving the nonlinear time distributed-order reaction-diffusion equation.

Example 5.2. On the space-time domain $[0, 1]^2\times[0, \frac12]$, we choose $\omega(\alpha) = \Gamma(5-\alpha)$, nonlinear term $m(u) = \arctan(u)$, source term

and the exact solution

In Table 3, by taking fixed temporal step length $\Delta t = 1/300$, fractional parameter $\Delta\alpha = 1/600$, and changed ${\widetilde{H}}_x = {\widetilde{H}}_y = 1/2, \, 1/3, \, 1/4, \, 1/5$, $1/6$, $h_x = h_y = 1/4, \, 1/9, \, 1/16, \, 1/25$, $1/36$, we show the error estimation results, the second-order spatial convergence rate and calculation time of $u$, from which one can see that we can obtain the similar results as that shown in Example 5.1.

Remark 5.1. As pointed out in some works, when the exact solution $u(x, y, t)$ has a strong singularity, the WSGD scheme maybe cannot preserve second-order approximation accuracy in time, so the initial correction technique presented in ^[44,45] can be considered.

6.   Conclusions and future works

In this paper, a two-dimensional nonlinear time distributed-order fractional sub-diffusion equation is solved by a proposed two-grid ADI FE method based on the WSGD operator. The unconditional stability and optimal error estimates with second-order convergence rate in spatial direction are obtained. The comparison of the numerical solution and the exact solution is made to demonstrate the efficiency of the numerical method. Compared with the traditional FE method, computing time can be saved and the storage space can be reduced in this method. Therefore, this method has further research value.

In the future, this method can be used to numerically solve nonlinear distributed-order time-space fractional sub-diffusion equations and nonlinear Volterra integro-differential equation with weakly singular kernel ^[42,43], and an ADI method with two-grid finite volume element method ^[41] for solving fractional models will be developed in another work.

7.   Acknowledgements

This work is supported by the Inner Mongolia University of Technology 2021 Research Start-up Grant Program (DC2200000908), the Natural Science Foundation of Inner Mongolia (2021MS01018), the Program for Innovative Research Team in Universities of Inner Mongolia Autonomous Region (NMGIRT2207), the Scientific research innovation project of postgraduate of Inner Mongolia University (11200-5223737).

Conflict of interest

The authors declare there is no conflict of interest.