Two-grid $H^1$-Galerkin mixed finite elements combined with $L1$ scheme for nonlinear time fractional parabolic equations

1.   Introduction

We focus on the following nonlinear time fractional partial differential equation (TFPDE):

where $J = (0, T]$ and $\Omega\subset \mathbb{R}^2$ is a convex polygonal with boundary $\partial\Omega$, the nonlinear function $f(y)$ fulfills

The Caputo derivative ^[1] $\mathcal{D}_t^\alpha$ defined as follows:

Here, $\Gamma(z) = \int_0^{+\infty}\kappa^{z-1}e^{-\kappa}d\kappa$ is the usual Gamma function.

Due to the excellent performance in describing memory and genetic properties, TFPDEs have been widely used to model anomalous diffusion, electromagnetic wave, hydrodynamics, signal processing, material science etc. Unfortunately, the analytical solutions of most TFPDEs are almost impossible to calculate. For the past twenty years, many researchers proposed extensive numerical methods of TFPDEs. Related books can be found in ^[2,3,4]. In the current literatures, we can see finite difference ^[5,6,7,8], spectral method ^[9,10], finite volume ^[11], finite element (FE) ^[12,13,14], weak FE ^[15], mixed finite element (MFE) ^{[16,17,18,19]}, virtual element method ^[20], fast algorithm ^[21,22,23], higher-order schemes ^{[24,25,26,27]} and so on, which mostly focused on the linear fractional differential equation case. Recently, Zheng et al. ^[28] established the stability and error estimates of fully discrete FE for a hidden-memory variable-order time-fractional optimal control model.

In recent years, scholars pay much attention to nonlinear TFPDEs. In ^[29], Li et al. considered a numerical approximation of nonlinear TFPDEs. In ^[30], the authors introduced a Newton linearized Galerkin FE approximation for nonlinear TFPDEs with non-smooth solutions. The two-grid method was first proposed by Xu ^[31] for solving nonsymmetric and nonlinear partial differential equations (PDEs) under finite element approximation. As an effective numerical method, it has been expanded to solve various nonlinear fractional PDEs ^{[32,33,34,35]}.

For the classic MFE, its discrete space pairs must satisfy the Ladyženskaja-Babuška-Brezzi (LBB) condition. To avoid the LBB constraints, Pehlivanov et al. ^[36] developed a least-squares MFE, Pani ^[37] presented an $H^1$-Galerkin mixed finite element (HMFE) approximation by requiring extra regularity on the solution and Yang ^[38] proposed a splitting positive definite MFE by separating the pressure equation from flux equation. Recently, Liu and Hou investigated the HMFE approximation of one-dimensional TFPDEs ^[39,40] and semilinear parabolic integro-differential equations ^[41], respectively. Our aim of this article is to provide a two-grid algorithm (TGA) for the nonlinear problem (1.1) by HMFE combined with $L1$ scheme discretization and analyze its stability and convergence.

The layout of the paper is as follows. We give a fully discrete HMFE approximation scheme (FDHAS) of the problem (1.1) in Section 2. The stability and convergence of the discrete scheme is analyzed in Section 3. The TGA of (1.1) and its stability and convergence are discussed in Section 4. In Section 5, we give several examples to verify our theoretical findings.

2.   FDHAS of nonlinear TFPEDs

The FDHAS of the nonlinear TFPDE (1.1) will be constructed in this section. Throughout this article, we use standard notations for Sobolev spaces $W^{m, q}(\Omega)$ with a semi-norm $|w|_{m, q}$ and a norm $\|w\|_{m, q}$. For $q = 2$, let $H^m(\Omega) = W^{m, 2}(\Omega), \ H_0^m(\Omega) = \{w\in H^m(\Omega): w|_{\partial \Omega} = 0\}$, $\|\cdot\| = \|\cdot\|_{0, 2}, \, \|\cdot\|_{m} = \|\cdot\|_{m, 2}$. $L^p(J; W^{m, q}(\Omega))$ be all $L^p$ integrable function space from $J$ into $W^{m, q}(\Omega)$ with norm $\|w\|_{L^p(J; W^{m, q}(\Omega))} = \Big(\int_0^T||w||_{W^{m, q}(\Omega)}^pdt\Big)^{1/p}\, \text{for}\, p\in [1, \infty),$ and the standard modification for $p = \infty$. In addition, $c > 0$ and $C > 0$ are generic constant.

Let us assume that the Eq (1.1) has a solution $y(t, x)$ such that

Set $W = H_0^1(\Omega)$ and $\pmb{V} = H(\mbox{div}, \Omega)$ with

Like in ^[42], an $H^1$-Galerkin mixed variational formulation of (1.1) is: Find $\{y, \pmb{Y}\} : (0, T] \rightarrow W\times\pmb{V}$, for any $\pmb{v}\in\pmb{V}$ and $w\in W$, such that

Let $\mathcal{T}_h$ be a uniform rectangular partition of $\Omega$, where $h = \max\limits_{e\in \mathcal{T}_h}\{h_e\}$ denotes the spatial mesh size and $h_e$ deotes the diameter of element $e$. Associated with the rectangular partition $\mathcal{T}_h$ of $\Omega$, we define subspaces $W_h\times\pmb{V}_h\subset W\times\pmb{V}$ as follows ^[43]:

where $Q_{m, n}(e)$ be the space of polynomials with degree no more than $m$ and $n$ in the $x_1$ direction and the $x_2$ direction, respectively.

Then a semidiscrete HMFE approximation of (1.1) is: Find $\{y_h, \pmb{Y}_h\} : (0, T] \rightarrow W_h\times\pmb{V}_h$, for any $\pmb{v}_h\in\pmb{V}_h$ and $w_h\in W_h$, such that

where ^[41] $P_h:W\rightarrow W_h$, which fulfills: For any $w_h\in W_h$ and $\varphi\in W$

We introduce ^[43,44] $\Pi_h:\pmb{V}\rightarrow \pmb{V}_h$, which fulfills: For any $\pmb{v}_h\in\pmb{V}_h$ and $\pmb{\psi}\in\pmb{V}$

Since the solution may be weakly singular at $t = 0$, we use the graded mesh for time discretization. For $n = 0, 1, \cdots, N$ with $N\in\mathbb{Z}^+$, let $t_n = T(n/N)^\gamma$, where $\gamma\geq 1$ will be adapted to the strength of the singularity and chosen by the user and the temporal mesh size $\tau = \max\{\tau_n\}_{n = 1}^N$ with $\tau_n = t_{n}-t_{n-1}$. We set $\varphi^n = \varphi(t_n, x)$ for $n = 0, 1, \cdots, N$. The $L1$ approximation scheme is given by ^[2]:

where

Lemma 2.1. (^[22]) If $|\varphi^{\prime\prime}(t)|\leq Ct^{\alpha-2}, 0 < t\leq T$, there exists a constant $C$ independent of $\tau$, such that

where $r = \min\{2-\alpha, \gamma\alpha\}$.

Then the FDHAS of (1.1) is as follows: Find $\left(y_h^n, \pmb{Y}_h^n\right)\in W_h\times\pmb{V}_h, n = 1, 2, \cdots, N$ for any $\pmb{v}_h\in\pmb{V}_h$ and $w_h\in W_h$, such that

3.   Stability and convergence analysis of the FDHAS

The stability and convergence of the FDHAS (2.13)–(2.15) will be analyzed in this section. The discrete Grönwall lemma will be used in the following analysis.

Lemma 3.1. (^[26,45]) Let $\{\xi^n\}_{n = 1}^{N}$, $\{g^n\}_{n = 1}^{N}$ and $\{\lambda_n\}_{n = 0}^{N-1}$ be given nonnegative sequences. If there is a constant $\Lambda$ independent of $\tau$ satisfies

Then, for any nonnegative sequence $\{V^k\}_{k = 0}^N$ and $1\leq n\leq N$ such that

it holds that

where $E_{\alpha}(z) = \sum\limits_{k = 0}^\infty\frac{z^k}{\Gamma(1+k\alpha)}$ is the Mittag-Leffler function.

3.1.   Stability analysis

We first show the stability of the solution to (2.13)–(2.15).

Theorem 3.1. Let $\left(y_h, \pmb{Y}_h\right)$ be the solution of (2.13)–(2.15) and all the conditions in Lemma 3.1 are valid. Then

Proof. Choosing $\pmb{v}_h = 2\pmb{Y}_h^n$ in (2.13) leads to

where we have used the Taylor's formula, $\varepsilon$-Cauchy inequality and the condition (1.2).

From the definition of $\mathcal{D}_N^\alpha\pmb{Y}_h^n$ and Hölder inequality, we can derive

Taking $w_h = y_h^n$ in (2.14) and using Hölder's inequality, we have

Employing the Poincaré inequality and (3.5)–(3.7), we obtain

Then (3.3) follows from Lemma 3.1 and (3.8). From (3.3) and (3.7), we can easily arrive at (3.4). □

3.2.   Convergence analysis

From (2.1)–(2.3), (2.11) and (2.13)–(2.15), we obtain error equations

For convenience, we set

Theorem 3.2. Let $\left(y, \pmb{Y}\right)$ and $\left(y_h, \pmb{Y}_h\right)$ be the solutions of (2.1)–(2.3) and (2.13)–(2.15), respectively. Suppose that $y\in L^\infty(J; H^2(\Omega))\cap H^2(J; L^2(\Omega))$, $\pmb{Y}\in (L^\infty(J; H^1(\Omega))^2\cap (H^2(J; H^1(\Omega)^2$ and all the conditions in Lemma 2.1 and Theorem 3.1 are valid. Then for $n = 1, 2, \cdots, N$, there hold

Proof. From the definitions of the projection operators $P_h$ and $\Pi_h$ and error Eqs (3.9)–(3.10). For any $\pmb{v}_h\in\pmb{V}_h$ and $w_h\in W_h$, we have

Taking $\pmb{v}_h = \pmb{\vartheta}^n$ in (3.13), then there yields

According to the mean-value theorem, the conditions (1.2) and (2.7), we obtain

It follows from Hölder's inequality and (2.12) that

From (2.10) and Poincaré inequality, we have

Using (3.15)–(3.18), Hölder's inequality and $\varepsilon$-Cauchy inequality, we get

Noting that $\left(\mathcal{D}_N^\alpha\pmb{\vartheta}^n, \pmb{\vartheta}^n\right)\geq\frac{1}{2}\mathcal{D}_N^\alpha\|\pmb{\vartheta}^n\|^2$, from (3.16)–(3.19), we have

Choosing $w_h = \eta^n$ in (3.14), we get

Using (3.21), Cauchy-Schwartz inequality, (2.9) and Poincaré inequality, we derive

According to (3.22), $\pmb{Y}\in (L^\infty(J; H^2(\Omega))^2$ and Poincaré inequality, we get

Combining (3.20) and (3.23) yields

By (3.24) and Lemma 3.1, we obtain

Noting that $\pmb{Y}\in (L^\infty(J; H^2(\Omega))^2$ and using triangle inequality, (2.8) and (3.25), we derive

From (3.23) and (3.25), we have

It follows from (2.7), triangle inequality and (3.27) that

Thus, the proof is completed. □

4.   TGA and error estimates

A two-grid HMFE algorithm to solve the nonlinear TFPDE (1.1) will be proposed in this section. Let ${\mathcal {T}}_H$ and ${\mathcal {T}}_h$ be two uniform rectangular partitions of $\Omega$ with different size $H$ and $h \; (h\ll H)$. The fine mesh ${\mathcal {T}}_h$ is obtained by uniformly refined the coarse mesh ${\mathcal {T}}_H$. Associated with ${\mathcal {T}}_H$ and ${\mathcal {T}}_h$ are HMFE spaces $W_H\times\pmb{V}_H$ and $W_h\times\pmb{V}_h$, respectively. It is obvious that $W_H\times\pmb{V}_H \subset W_h\times\pmb{V}_h$. We show the TGA as follows:

Two-grid algorithm (TGA).

Step 1. On ${\mathcal {T}}_H$: Find $\left(y_H^n, \pmb{Y}_H^n\right)\in W_H\times\pmb{V}_H$ for $n = 0, 1, \cdots, N$ and any $\pmb{v}_H\in\pmb{V}_H$, $w_H\in W_H$, such that

Step 2. On ${\mathcal {T}}_h$: Given $y_H^n$ for $n = 1, 2, \cdots, N$, find $\left(y_h^{*, n}, \pmb{Y}_h^{*, n}\right)\in W_h\times \pmb{V}_h$ for $n = 0, 1, \cdots, N$ and any $\pmb{v}_h\in\pmb{V}_h$, $w_h\in W_h$, such that

4.1.   Stability analysis

Now, we analyze the stability of the TGA.

Theorem 4.1. Let $\left(y_h^{*}, \pmb{Y}_h^{*}\right)$ be the solution of the two-grid algorithm (4.1)–(4.6) and all the conditions in Lemma 3.1 are valid. Then we have

Proof. Setting $\pmb{v}_h = 2\pmb{Y}_h^{*, n}$ in (4.4), we get

where we have used the Taylor's formula, $\varepsilon$-Cauchy inequality and the condition (1.2).

Choosing $w_h = y_h^{*, n}$ in (4.5) and using Hölder's inequality, we obtain

Noting $2\left(\mathcal{D}_N^\alpha\pmb{Y}_h^{*, n}, \pmb{Y}_h^{*, n}\right)\geq \mathcal{D}_N^\alpha\left\|\pmb{Y}_h^{*, n}\right\|^2$ for $n = 1, 2, \cdots, N$. Applying the Poincaré inequality, (3.7), (4.9) and (4.10), we derive

Then (4.7) follows from (4.11) and Lemma 3.1. According to (4.7) and (4.10), it is easy to get (4.8).□

4.2.   Error estimates

Subtracting (4.4) and (4.5) from (2.1) and (2.2), for any $\pmb{v}_h\in\pmb{V}_h$ and $w_h\in W_h$, we can obtain equations

In order to conveniently derive error estimation results, we set

Theorem 4.2. Let $\left(y, \pmb{Y}\right)$ and $\left(y_h^{*}, \pmb{Y}_h^{*}\right)$ be the solutions of (2.1)–(2.3) and (4.1)–(4.6), respectively. Suppose that $y\in L^\infty(J; H^2(\Omega))\cap H^2(J; L^2(\Omega))$, $\pmb{Y}\in (L^\infty(J; H^2(\Omega))^2\cap (H^2(J; L^2(\Omega)^2$ and all the conditions in Theorem 3.2 and Theorem 4.1 are valid. Then for $n = 1, 2, \cdots, N$, there hold

Proof. Subtracting (4.2) from (2.2) and utilizing the definition of $P_H$, we have

Choosing $w_H = P_Hy^n-y_H^n$ in (4.16), we derive

From Poincaré inequality, (4.17) and Theorem 3.2, we obtain

According to triangle inequality, (2.7) and (4.18), we arrive at

Using (4.12) and (4.13) and the definitions of $P_h$ and $\Pi_h$, for any $\pmb{v}_h\in\pmb{V}_h$ and $w_h\in W_h$, we get

From the Taylor expansion formula, we have

where $\delta$ between $y^n$ and $y_H^n$.

Substituting (4.22) into (4.20) and choosing $\pmb{v}_h = \pmb{\rho}^n$ yields

We can estimate the right-hand side of (4.23) as follows:

and

and

where we have used (2.7), (2.10), Lemma 2.1, Hölder inequality, embedding theorem, Theorem 3.2 and Poincaré inequality.

Taking $w_h = \xi^n$ in (4.21), we get

From (2.9), (4.27) and Cauchy-Schwartz inequality, we derive

By (4.28) and Poincaré inequality, we have

Combining (4.23)–(4.29) and using $(a+b)^4\leq 8(a^4+b^4)$ with $a, b > 0$, we yield

It follows from Lemma 3.1 and (4.30) that

Using triangle inequality, (2.8) and (4.31), we arrive at

From (4.29) and (4.31), we have

It follows from (2.7), triangle inequality and (4.33) that

We complete the proof of Theorem 4.2. □

5.   Numerical experiment

Several examples are provided to demonstrate our theoretical findings in this section. All numerical examples will be solved by the FDHAS described as (2.13)–(2.15) and TGA described as (4.1)–(4.6), where the program codes are based on AFEPack ^[46]. Let $J = (0, 1]$ and $\Omega = (0, 1)^2$. We numerically solve the following nonlinear time fractional equation:

For simplicity, we define $|||\psi||| = \max\limits_{0\leq n\leq N}\{\|\psi^n\|\}$. The convergence rate are computed by $Rate = \frac{\ln\left(|||\psi_{k+1}|||\right)-\ln\left(|||\psi_k|||\right)}{\ln(s_{k+1})-\ln(s_k)}$, where $|||\psi_{k+1}|||\left(|||\psi_k|||\right)$ is the error with the spatial or temporal mesh step $s_{k+1}\left(s_k\right)$.

Example 1. The initial condition and the right function $g(t, x)$ are suitably chosen such that $y(t, x) = t^2\sin(\pi x_1)\sin (\pi x_2)$ and the nonlinear term $f(y) = y(1-y)$.

When the spatial step $h = \frac{1}{100}$ and $H = \sqrt{h} = \frac{1}{10}$ are fixed, in Table 1, we present the errors $|||\pmb{Y}_h-\pmb{Y}|||$ and $|||\pmb{Y}_h^*-\pmb{Y}|||$, the temporal convergence orders and the CPU time for $\alpha = 0.4, 0.5$ and $0.8$ by using the general FDHAS (2.13)–(2.15) and TGA (4.1)–(4.6), respectively.

The numerical results in Table 1 show that the TGA (4.1)–(4.6) can save significant computational costs compared with FDHAS (2.13)–(2.15) without losing accuracy. Fixed $N = 1000$, $\gamma = \frac{2-\alpha}{\alpha}$ and $h = H^2$, the results in Table 2 reflect $|||\pmb{Y}_h-\pmb{Y}||| = \mathcal{O}\left(h\right)$, $|||y_h-y||| = \mathcal{O}\left(h^2\right)$ and $|||\pmb{Y}_h^*-\pmb{Y}||| = \mathcal{O}\left(h\right)$. The convergence rate results in time and space direction are consistent with our theoretical results.

Example 2. The initial condition and the right function $g(t, x)$ are suitably chosen such that $y(t, x) = \left(E_\alpha(-t^\alpha)+t^3\right)x_1(1-x_1)x_2(1- x_2)$ and the nonlinear term $f(y) = y^3$.

In Table 3, we also give the numerical results with the fixed spatial step $h = \frac{1}{100}$ and $H = \sqrt{h} = \frac{1}{10}$ for $\alpha = 0.4, 0.5$ and $0.8$ by utilizing the general FDHAS (2.13)–(2.15) and TGA (4.1)–(4.6), respectively. For fixed $N = 1000$, $\gamma = \frac{2-\alpha}{\alpha}$ and $h = H^2$, we show the numerical results in Table 4. The numerical results demonstrate that the TGA is more efficient than the FDHAS. It is agreement with our theoretical analysis.

6.   Conclusions

In this paper, we proposed the TGA for the nonlinear TFPDEs (1.1) discretized by $H^1$-Galerkin mixed finite element on spatial rectangular mesh combined with $L1$ scheme on temporal graded mesh. The stability and optimal convergence of the TGA are rigorously proved. Our theoretical results seem to be new in the literature. Numerical experimental results show that the TGA (4.1)–(4.6) can save a lot of computing cost compared with FDHAS (2.13)–(2.15) without losing accuracy. Although our TGA in this paper focuses on a two-dimensional case, it can be directly applied to three-dimensional problems. Future work includes the developments of two-grid finite element methods combined with some higher-order schemes or fast algorithm for nonlinear TFPDEs.

Use of AI tools declaration

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

Acknowledgments

This work is supported by the Scientific Research Foundation of Hunan Provincial Department of Education (20A211) and the Natural Science Foundation of Hunan Province (2020JJ4323).

Conflict of interest

The authors declare that there are no conflicts of interest.