Iterative two-level algorithm for nonsymmetric or indefinite elliptic problems

We consider the following Dirichlet boundary value problems for general second-order elliptic partial differential equations

where the coefficient $\boldsymbol{\alpha}(x)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ is smooth functions on $\bar{\Omega}$, and satisfies consistent ellipticity, namely, there are minimum and maximum eigenvalues $\alpha_{*}$ and $\alpha^{*}$ respectively, satisfying $\alpha_{*}|\xi|^{2}\leqslant\mathbf{\alpha}(x)\xi\cdot\xi\leqslant\alpha^{*}|\xi|^{2},~{}~{}\forall\xi\in\mathbb{R}^{d}$. Both $\boldsymbol{\beta}(x)\in\mathbb{R}^{d}$ and $\gamma(x)\in\mathbb{R}^{1}$ are smooth functions on $\bar{\Omega}$, and $f\in L^{2}(\Omega)$ is a given function. Noting that when $\boldsymbol{\beta}(x)\neq\boldsymbol{0}$, the continuous variational problems (CVP) of (1.1)-(1.2) are nonsymmetric, and when $\gamma(x)<0$, the CVP of (1.1)-(1.2) may be indefinite.

The two-grid (TG) algorithm was first introduced by Xu [6] for solving nonasymmetric or indefinite problems. The main idea of the TG algorithm is to solve the original problems on the coarse mesh to obtain an approximate finite element solution, and then use the approximate solution to solve the corresponding linear symmetric positive definite problems on the fine mesh. Over the last two decades, the TG algorithm has been widely used to solve many problems, such as nonlinear elliptic problems [1, 10], semilinear parabolic equations [5], nonlinear parabolic equations [3, 2], Poisson-Nernst-Planck problems [7], and Maxwell equations [8, 9].

The main objective of this paper is to propose a new iterative two-level algorithm. Compared with the traditional iterative two-grid algorithm, our algorithm needs only one mesh and takes less CPU time to achieve the same accuracy.

In this work, we first present elliptic problems discretized by the finite element method. Then, we propose an iterative two-level algorithm. Finally, some numerical results are presented to illustrate the efficiency of the proposed algorithms.

Given $S\subset\mathbb{R}^{d}(d=2,3)$, we denote $W^{1,2}(S)$ as the standard Sobolev space with norm $\|\cdot\|_{1,2,S}$. For simplicity of notation, we denote $\|\cdot\|_{1}=\|\cdot\|_{1,2,\Omega}$, and $H^{1}_{0}(\Omega):=\{u\in H^{1}(\Omega):u|_{\partial\Omega}=0\}$ in the sense of trace. We define the following two bilinear forms

Then, we obtain the CVP of (1.1)-(1.2): Find $u\in H_{0}^{1}(\Omega)$, such that

We assume that $\Omega$ is partitioned by a quasi-uniform division $\mathcal{T}_{h}=\{\tau_{i}\}$. By this we mean that $\tau_{i}$’s are simplexes of the size $h$ with $h\in(0,1)$ and $\bar{\Omega}=\cup_{i}\tau_{i}$. For the given quasi-uniform division $\mathcal{T}_{h}$, the conforming finite element space is defined as follows

where $\mathbb{P}_{l}(\tau)$ is the space of polynomial of degree not greater than a positive integer $l$ on the subdivision element $\tau$. The discrete variational problems of (2.3) is to find $u_{h}^{l}\in V_{h}^{l}$, such that

The basic mechanism of the classical iterative two-grid algorithm is two quasi-uniform tetrahedral nested meshes of $\Omega$, namely the fine space $V_{h}^{k,l}$ and the coarse space $V_{H}^{k,l}$, with two different meshes sizes $h$ and $H(h<H)$. Furthermore, in the application given in the succeeding text, we shall always assume that $H=O\left(h^{\lambda}\right)\text{ for some }0<\lambda<1.$

Assume that $u_{h}^{l,k+1}\in V_{h}^{l}$ is the solution obtained by Algorithm 3.1 with $k\geq 1$, then we have (see Theorem 4.4 of [6])

which means that the two-grid solution given by Algorithm 3.1 can effectively approximate the finite element solution $u_{h}^{l}$ of (2.4).

Next, we consider replacing $V^{l}_{h}$ with $V_{H}^{s}(s\geq l+1)$ in the second step of Algorithm 3.1, and obtain the following algorithm.

In this section, numerical experiments are carried out to verify the effectiveness of the iterative two-level algorithm. We performed all experiments for our iterative two-level algorithm with the help of the software package Fenics [4].

We first partition the $x-$ axis and $y-$ axis of the domain $\Omega$ into equally distributed $M$ subintervals, then divide each square into two triangles by using the line with slope $-1$. Hence, we obtain a sequence of nested and structured grids and the corresponding meshes as $\mathcal{T}_{H}$ with $H=1/M$, where $M\geq 2$ is an integer, see Figure 1. We choose the piecewise conform $l$ order finite element spaces $V_{H}^{l}$ based on the meshes $\mathcal{T}_{H}$. For Algorithm 3.1, we choose $h=H^{2}$.

From Table 2, we can observe that both our algorithm and the traditional iterative mesh can reach the optimal convergence order. And to achieve the same accuracy, our algorithm uses less CPU time.

From the Tables 2- 3, it can be found that when the degree of the coarse space polynomial degree $l=3$ is fixed, as the degree of the fine space polynomial increases once, the convergence order of the error estimates in $H^{1}$-norm for the solution of the algorithm 3.2 increase by one order. It can be seen that the errors in $H^{1}$-norm for the solution of the algorithm 3.2 depend on the value of $s$.

The same conclusions can be observed in Tables 4 and 5 as in Tables 2 and 3.

The first, second and fourth authors are supported by the National Natural Science Foundation of China (No. 12071160). The second author is also supported by the National Natural Science Foundation of China (No. 11901212). The third author is supported by the National Natural Science Foundation of China (No. 12161026), Guangxi Natural Science Foundation (No. 2020GXNSFAA159098).