Two-grid finite element method with an H2N2 interpolation for two-dimensional nonlinear fractional multi-term mixed sub-diffusion and diffusion wave equation

1.   Introduction

In recent decades, fractional partial differential equations (FPDEs) have become popular for modeling anomalous transport processes. In virtue of the difficulty for looking for the exact solutions of FPDEs, the construction of numerical methods for FPDEs have been studied extensively by many scholars which mainly cover the finite difference method, finite element method, spectral method and so on. Unlike general PDEs, it has been found that fractional derivatives are nonlocal, with long memory and weak singular kernels, and it takes more time to solve fractional differential equations when using some well-known L-type discretization formulas. In order to increase computational efficiency, some fast algorithms have been presented, such as the sum-of-exponentials technique ^[1,2] and the two-grid technique ^{[3,4,5,6,7,8,9]}.

There have been many studies on the approximation of time fractional order equations. Schumer ^[10] first developed the fractional-order mobile/immobile model in which the time drift term $\partial u/\partial t$ is added to describe the motion time, thus helps to distinguish the status of particles conveniently. Zhang and his coauthors ^[11] proposed a difference scheme combining the compact difference approach for spatial discretization and the alternating direction implicit (ADI) method in the time stepping for the two-dimensional time fractional diffusion-wave equation. In ^[12], the authors established a fully-discrete approximate scheme for the 2D multi-term time-fractional mixed diffusion and diffusion-wave equations with a spatial variable coefficient by using the linear triangle finite element method in space and classical L1 time-stepping method combined with the Crank–Nicolson scheme in time. Sun et al. ^[13] gave an L-type difference scheme for time-fractional mixed sub-diffusion and diffusion wave equations, the new analytical technique with a min$\{2-\alpha, 3-\beta\}$ order accuracy in the discrete $H^1$ norm and second order accuracy in space, where $\alpha\in (0, 1), \; \beta\in (1, 2)$. In recent years, Shen et al. ^[14] derived the H2N2 method to develop a known numerical differential formula of the Caputo fractional derivative to obtain a $(3-\alpha)$ order temporal convergence rate with $\alpha\in(1, 2)$. Moreover, Shen and his coauthors ^[15] also gave the finite difference methods based on an H2N2 interpolation for two-dimensional time fractional mixed sub-diffusion and diffusion-wave equations with a $(3-\beta)$ order temporal convergence rate and second order accuracy in space.

With further research, it is found that nonlinear partial differential equations can describe some model problems more effectively compared to linear fractional differential equations. As we know, numerical methods for nonlinear parabolic fractional partial differential problems have been extensively studied, especially for the two-grid approximation of nonlinear reaction terms. In ^[4,5], Xu utilized the two-grid method to discretize asymmetric, indefinite, and nonlinear partial differential equations. Afterward, many authors devoted themselves to the research of the two-grid method. For example, Chen and Chen ^[6] investigated a scheme for nonlinear reaction-diffusion equations by using the mixed finite element methods, in which the two-grid method (TGM) studied provided a new approach to take the advantage of some nice properties hidden in a complex problem. Based on the finite difference approximation in the time direction and the finite element method in the spatial direction, a class of time fractional fourth order reaction diffusion problems with nonlinear reaction terms was solved by Liu et al. ^[16]. Almost simultaneously, Liu et al. ^[17] considered a nonlinear fractional Cable equation by a two-grid algorithm combined arriving at a second-order convergence rate independent of fractional parameters $\alpha\; (0 < \alpha < 1)$ and $\beta\; (0 < \beta < 1)$ with finite element method. In addition, Bi et al. ^[18] established a discontinuous Galerkin finite element scheme for the second-order nonlinear elliptic problem, using piecewise polynomials of degree $r\ge2$ based on pointwise error estimation, and provided an error estimate of the energy norm of a two-grid mesh algorithm. Li and Rui ^[8] introduced and analyzed a two-grid block-centered finite difference scheme for the nonlinear time-fractional parabolic equation with the Neumann condition, and the error estimates established on a non-uniform rectangular grid with the discrete $L^{\infty}(L^2)$ and $L^2(H^1)$ errors were $O(t^{2-\alpha}+h^2+H^3)$. Jin et al. ^[19] demonstrated a general criterion for showing the fractional discrete Gr$\ddot{{\mathrm{u}}}$onwall inequality and verified it for the L1 scheme and convolution quadrature generated by backward difference formulas for time-fractional nonlinear parabolic partial differential equations. In 2019, Li et al. ^[20] presented and developed a Newton linearized Galerkin finite element method to solve nonlinear time fractional parabolic problems with non-smooth solutions in time direction. Recently, Chen et al. ^[7] constructed and studied a two-grid finite element method for 2D nonlinear time fractional two-term mixed sub-diffusion and diffusion wave equations and got min(2-$\alpha, 3-\beta)$ order in the time direction where $\alpha\in(0, 1), \; \beta\in(1, 2)$.

However, to our knowledge, there is no two-grid finite element method with an H2N2 interpolation for two-dimensional nonlinear fractional multi-term mixed sub-diffusion and diffusion wave equations, and our aim is to present the corresponding algorithm in this paper. We consider the numerical solution of the following nonlinear time-fractional multi-term mixed sub-diffusion and diffusion wave equations:

where $M_{1}, M_{2}$ are positive integers, $p_{i}, q_{j}$ are nonnegative constants, $\Omega\subset\mathbb{R}^2$ is a bounded convex polygonal region with boundary $\partial\Omega, {\boldsymbol{x}} = [x, y]$ and $g(\cdot)$ is twice continuously differentiable. $\Delta$ is a Laplace operator, and $u_{0}({\boldsymbol{x}}), u_{t}^0({\boldsymbol{x}})$ are given as sufficiently smooth functions. The Caputo fractional derivative $\partial^{\alpha_{i}}_{t}u({\boldsymbol{x}}, t), \partial^{\beta_{j}}_{t}u({\boldsymbol{x}}, t)$ is defined by

with $0 < \alpha_{0} < \dots < \alpha_{i} < \dots < \alpha_{M_{1}} < 1\; (i = 1, 2, \dots, M_{1})$ and $1 < \beta_{0} < \dots < \beta_{j} < \dots < \beta_{M_{2}} < \; 2\; (j = \; 1, 2, \dots, M_{2})$. For convenience, we define $\alpha_{M_{1}} = \alpha, \beta_{M_{2}} = \beta.$ Throughout this paper, the notation $C$ denotes a generic constant, which may vary at different occurrences, but it is always independent of the coarse grid mesh size $H$, fine grid mesh size $h$ and time step size $\tau$.

The remaining outline of the article is constructed as follows. In section two, we propose some preliminary knowledge of fractional derivatives and some lemmas on time approximations. The unconditional stability and the error estimates of the two-grid finite element method (FEM) are discussed in section three. In section four, numerical examples are presented to demonstrate the efficiency of the theoretical results. Finally, a brief conclusion of the article was provided.

2.   Construction of fully-discrete scheme

Let $W_p^m(\Omega)$ be the standard Sobolev space with a norm $\left\|\cdot\right\|_{m, p}$ given by $\left\|v\right\|_{m, p}^p = { \sum_{\left|\alpha\right|\le m}}\left\|D^\alpha v\right\|^p_{L^P(\Omega)}$. For $p = 2$, we denote $H^m(\Omega) = W_2^m(\Omega)$ and $H^1_0(\Omega)$ to be the subspace of $H^1(\Omega)$ consisting of functions with a vanishing trace on $\partial\Omega$. $\left\|\cdot\right\|_m = \left\|\cdot\right\|_{m, 2}$ and $\left\|\cdot\right\| = \left\|\cdot\right\|_{0, 2}$. Take an integer $N$ and denote $\tau = \frac{T}{N}, \; t_n = n\tau, \; t_{n-\frac{1}{2}} = \frac{1}{2}(t_n+t_{n-1})$. For a sequence of smooth funtions $\{u(t)\}_{n = 0}^N$ on $[0, T]$, we denote

For any function $u$ defined on the interval $[0, t_1]$, we consider the Hermite quadratic interpolation polynomial $H_{2, 0}(t)$. On the interval $[t_{k-1}, t_{k+1}]\; (1\le k\le N-1)$, the Newton quadratic interpolation polynomial $N_{2, k}(t)$ will be used. With the Caputo-fractional derivatives (1.1) on the time grid points of the form $\{t_\frac{1}{2}, t_\frac{3}{2}, \dots, t_{N-\frac{1}{2}}\}$ ^[14,15], we can obtain

where

The H2N2 interpolation for the Caputo derivative $\mathcal{D}_{t}^{\beta_j}u(t_{n-\frac{1}{2}})$ can be written as follows:

Here

Lemma 2.1. For $A_0^{(n, \alpha_i)}, a_k^{(n, \alpha_i)}, b_k^{(n, \beta_j)}, \; 0\le k\le n-1$ defined by (2.2)–(2.4), (2.6) and (2.7), we gain

Lemma 2.2. ^[14,15] Suppose $u(t)\in C^3[t_0, t_n]$, and denote $R_{\alpha_i}^{n-\frac{1}{2}} = \sum_{i = 1}^{M_{1}}p_{i}\partial_t^{\alpha_i}u(t_{n-\frac{1}{2}})-\sum_{i = 1}^{M_{1}}p_{i}\mathcal{D}_t^{\alpha_i}u(t_{n-\frac{1}{2}})$ then we can get

If $u(t)\in C^3[t_0, t_n]$, denote $R_{\beta_j}^{n-\frac{1}{2}} = \sum_{j = 1}^{M_{2}}q_{j}\partial_t^{\beta_j}u(t_{n-\frac{1}{2}})-\sum_{j = 1}^{M_{2}}q_{j}\mathcal{D}_t^{\beta_j}u(t_{n-\frac{1}{2}})$, then we have

The basic mechanism in this algorithm is the construction of two shape-regular triangulations of $\Omega$, which we denote by $\mathcal{T}_H$ and $\mathcal{T}_h$, with different mesh sizes $H$ and $h\; (h\ll H)$. This can be accomplished by successively refining the triangulation $\mathcal{T}_H$ to obtain the fine mesh $\mathcal{T}_h$. The corresponding finite element spaces are $V_H$ and $V_h$, which will be called coarse and fine space, respectively. We note that by such a construction we have $V_H\subset V_h$. Let $V_h$ be the two-dimensional subspace of $H^1_0(\Omega)$, which consists of continuous piecewise polynomials of degree $r(r\ge1)$ on $\mathcal{T}_h$ and $V_h^0 = \{v\in V_h, v|\partial\Omega = 0\}$, $V_H^0 = \{v\in V_H, v|\partial\Omega = 0\}$. {Note} that the corresponding weak formulation of (1.1) is to find $u({\boldsymbol{x}}, t):(0, T]\to H^1_0(\Omega)$, then

First, we derive the two-grid fully discrete scheme for problem (1.1). The process is divided into two steps. Step 1: On the coarse grid $\mathcal{T}_H$, for any $v_H\in V_H^0,$ find $u_H^n\in V_H^n$ for the following nonlinear system, such that

Step 2: On the fine grid $\mathcal{T}_h$, find $u_h^n\in V_h^n$ for the following linear system, then

Next, the two-grid fully discrete scheme based on the Newton iteration idea will be analyzed. In the two-grid algorithm, we solve the nonlinear factional equation on the coarse grid $\mathcal{T}_H$ to produce a rough approximation, and then use the rough approximation as the initial guess to slove the linearized equation on the fine grid $\mathcal{T}_h$.

3.   Unconditional stability and error estimate

Below, the stability and convergence of the two-grid fully discrete scheme will be derived. For the demand of analysis, we introduce some lemmas as follows:

Lemma 3.1. ^[21] Let $R_{\tilde{h}}:H^1_0(\Omega)\to V^0_{\tilde{h}}$, be the Rize-projection operator satisfying

For $u\in H^1_0(\Omega)\cap H^{r+1}(\Omega)$, it holds that

where $\tilde{h}$ is coarse grid step length $H$ or fine grid size $h$.

Lemma 3.2. ^[22] If $A_m, B_s$ are nonnegative real sequences and the sequence $Z_m$ satisfies

where $A_0\ge0$, then the sequence ${Z_m}$ satisfies

3.1.   Stability

In this subsection, we shall present the main algorithm of the paper. The energy method is applied to estabilish the stability of the two-grid fully discrete scheme. First, we derive the stability of the coarse grid system.

Theorem 3.1. For $u_H^n\in V_H^0$, the two-grid finite element analysis system (2.10) can obtain the following inequality

Proof. Suppose $\{u_H^{n-\frac{1}{2}}\}$ is the solution of a semi discrete scheme (2.10), the $v_H = \delta_tu_H^{n-\frac{1}{2}}\in V_H^0$, then the first item of (2.10) holds that

Using Eq (2.1), the second item of (2.10) follows that

By substituting (2.5) into (2.10), we can organize and obtain the third item of (2.10)

Substituting (3.3)–(3.5) into (2.10), and multiplying $2\tau$ on both sides of this equatity, we obtain

Denoting $E^0 = \left\|\nabla u_{H}^{0}\right\|^2$ and $E^n = \left\|\nabla u_{H}^{n}\right\|^2+\sum_{i = 1}^{M_{1}}\frac{\tau p_i}{\Gamma(2-\alpha_i)}\sum_{k = 1}^{n}a_{n-k}^{(n, \alpha_i)}\left\|\delta_{t}u_H^{k-\frac{1}{2}}\right\|^2+\sum_{j = 1}^{M_{2}}\frac{\tau q_j}{\Gamma(2-\beta_j)}\sum_{k = 1}^{n}b_{n-k}^{(n, \beta_j)}\cdot$ $\left\|\delta_{t}u_H^{k-\frac{1}{2}}\right\|^2$, we can calculate that

which, by induction, gives

Using the Cauchy-Schwarz inequality, the Young's inequality and the nonlinear property of $g(\cdot)$, it is shown that

where

Combing with the Lemma 2.1, above analysis uses

Substituting the estimates (3.8)–(3.10) into (3.5), and using the Poincar${\mathrm{\acute{e}}}$ inequality, we have

By the Lemma 3.2, we have the desired result.

Theorem 3.2. For the fine grid system (2.11), the following stable inequality for $u_h^n\in V_h^0$

Proof. Considering (2.11), we can deduce that

where $C, C_1, C_2$ are positive constants.

Let $v_h = \delta_{t}u_h^{n-\frac{1}{2}}, \; F^n = \left\|\nabla u_{h}^{n}\right\|^2+\sum_{i = 1}^{M_{1}}\frac{\tau p_i}{\Gamma(2-\alpha_i)}\sum_{k = 1}^{n}a_{n-k}^{(n, \alpha_i)}\left\|\delta_{t}u_h^{k-\frac{1}{2}}\right\|^2$ $+\sum_{j = 1}^{M_{2}}\frac{\tau q_j}{\Gamma(2-\beta_j)}\sum_{k = 1}^{n}b_{n-k}^{(n, \beta_j)}\left\|\delta_{t}u_h^{k-\frac{1}{2}}\right\|^2$. Using the recurrence relation, we have

It follows from (3.10) that

then combing with (3.8) and substituting (3.14) into (3.13) yields

Using the Poincar${\mathrm{\acute{e}}}$ inequality and the Lemma 3.2, we obtain the proof.

3.2.   Covergence

The convergence of the FEM system is considered by the energy method. First, we shall derive the convergence of the coarse grid system as follows:

Theorem 3.3. Let $u(t_n)\in H_0^1(\Omega)\cap H^{r+1}(\Omega), \; u_t\in H^2(\Omega), \; u_{tt}\in L^2(\Omega), \; u_{ttt}\in L^2(\Omega), \; u_H^n\in V_H^n$ and $u_H^0 = R_hu_0(x)$, then we have

Proof. We combine (1.1) with (2.10) to get

Denoting $\eta^n = u^n-R_hu^n, \; \theta^n = R_hu^n-u_H^n$ and choosing $v_h = \delta _t\theta^{n-\frac{1}{2}}$ in (3.15), we obtain

From the definitions of $\sum_{i = 1}^{M_1}p_i\mathcal{D}_t^{\alpha_i}u_H^{n-\frac{1}{2}}, \; \sum_{j = 1}^{M_2}q_j\mathcal{D}_t^{\beta_j} u_H^{n-\frac{1}{2}}$ and $\delta_tu_H^{n-\frac{1}{2}}$, we can have that

and

Since $g(\cdot)$ is twice continuously differentiable, we have

Substituting the above results into (3.16) and multiplying $2\tau$ on both sides of the resulting identity, it holds that

Denoting $G^0 = \left\|\nabla\theta^{0}\right\|^2$ and $G^n = \left\|\nabla\theta^{n}\right\|^2+\sum_{i = 1}^{M_1}\sum_{k = 1}^{n}\frac{\tau p_ia_{n-k}^{(n, \alpha_i)}}{\Gamma(2-\alpha_i)}\left\|\delta_{t}\theta^{k-\frac{1}{2}}\right\|^2+\sum_{j = 1}^{M_2}\sum_{k = 1}^{n}\frac{\tau q_jb_{n-k}^{(n, \beta_j)}}{\Gamma(2-\beta_j)}\left\|\delta_{t}\theta^{k-\frac{1}{2}}\right\|^2$, we have that

which, by induction, gives

By the Cauchy-Schwarz inequality, the Young's inequality and Lemma 2.2, we obtain

The definition of $\Xi$ is the same as (3.9). By Lemma 2.1, we have

Based on the above estimates and Lemma 3.1, we have

Combining (3.22)–(3.37) can give

By the Poincar${\mathrm{\acute{e}}}$ inequality, Lemma 3.1 and the Lemma 3.2, it follows that

The proof is complete.

Theorem 3.4. Let $u(t_n)\in H_0^1(\Omega)\cap H^{r+1}(\Omega), \; u_t\in H^2(\Omega), \; u_{tt}\in L^2(\Omega), \; u_{ttt}\in L^2(\Omega), \; u^n_h\in V_h^0$, and $u^0_h = R_hu_0({\boldsymbol{x}})$, then we arrive at the following

Proof. Subtracting (1.1) from (2.10), it yields

By the Taylor expansion, we can have the estimate that

then it follows from (3.40) that

Let $u^{n-\frac{1}{2}}-R_hu^{n-\frac{1}{2}} = \eta^{n-\frac{1}{2}}, \; R_hu^{n-\frac{1}{2}}-u_h^{n-\frac{1}{2}} = \theta^{n-\frac{1}{2}}$, and $v_h = \delta_t\theta^{n-\frac{1}{2}}$. Hence, (3.42) becomes

Defining $H^0 = \nabla\theta^{0}, \; H^n = \left\|\nabla\theta^{n}\right\|^2+\sum_{i = 1}^{M_1}\sum_{k = 1}^{n}\frac{\tau p_i}{\Gamma(2-\alpha_i)}a_{n-k}^{(n, \alpha_i)}\left\|\delta_{t}\theta^{k-\frac{1}{2}}\right\|^2$ $+\sum_{j = 1}^{M_2}\sum_{k = 1}^{n}\frac{\tau q_j}{\Gamma(2-\beta_j)}b_{n-k}^{(n, \beta_j)}\left\|\delta_{t}\theta^{k-\frac{1}{2}}\right\|^2$, we have following results:

We use the similar process of proof to the estimate (3.37). By the Poincar${\mathrm{\acute{e}}}$ inequality, Lemma 3.1, Lemma 3.2 and the estimate of (3.37), it follows that

which leads to

This concludes the proof.

Remark 4.1. In the estimate (3.44), we observe that the TGM algorithm can achieve the convergence rate $h^{r+1}$ as long as the mesh sizes satisfy $H = O(h^{\frac{1}{2}})$.

4.   Numerical examples

In this section, we present a numerical example to demonstrate the theoretical analysis and illustrate the efficiency of the algorithm discussed in section three. To investigate the spatial and temporal convergence order, we use a bilinear finite element approximation and the computation is performed by using Matlab.

Example 4.1. The following equation has an exact solution $u(x, y, t) = t^{3+\alpha+\beta}{\mathrm{sin}}\pi x{\mathrm{sin}}\pi y$:

where $\Omega$ is the unit square $(0, 1)\times (0, 1)$, $T = 1$ and $g$ is a known function. The domain $\Omega$ is divided into families $\mathcal{T}_H$ and $\mathcal{T}_h$ of quadrilaterals, and $V_H, V_h\subset H_0^1(\Omega)$ are linear spaces of piecewise continuous bilinear functions defined on $\mathcal{T}_H$ and $\mathcal{T}_h$, respectively.

Let $H_x = H_y = H, \; h_x = h_y = h$ and $h = H^2$. Tables 1–3 show that the spatial convergence rates in the $L^2$-norm of FEM and Algorithm 3.1 are both equivalent to two. The convergence results are in good agreement with the results $O(h^{r+1})$ of the theoretical analysis. Moreover, combining Tables 2–4, it can be seen that the TGM scheme will save more time than the general FEM scheme as $M$ increases.

Example 4.2. The following equation has the exact solution $u(x, y, t) = (1+t^{3+\alpha+\beta}){\mathrm{sin}}\pi x{\mathrm{sin}}\pi y$:

where $\Omega$ is the unit square $(0, 1)\times (0, 1)$, $T = 1$ and $g$ is a known function.

5.   Conclusions

In this paper, we proposed a new H2N2 formula to approximate the multi-term fractional derivative $\sum_{i = 1}^{M_{1}}p_{i}\partial_{t}^{\alpha_{i}}u({\boldsymbol{x}}, t), \; \sum_{j = 1}^{M_{2}}q_{j}\partial_{t}^{\beta_{j}}u({\boldsymbol{x}}, t), \; \alpha_{i}\in(0, 1), \; \beta_{j}\in(1, 2)$. Based on the H2N2 approximation in time and the finite element method for the spatial discretization, we have presented a fully discrete TGM scheme for two-dimensional nonlinear time fractional multi-term mixed sub-diffusion and diffusion wave equations and proved they are of second-order convergence in space and can reach the optimal convergence order $3-\beta$ in time, which is not related to $\alpha$. In future work, we will consider the L2-1$_\sigma$ formulation, which can reach second-order accuracy in time.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the State Key Program of the National Natural Science Foundation of China (11931003), and by the Grants (41974133).

Conflict of interest

There is no conflicts of interest between all authors.