Symmetrized two-scale finite element discretizations for partial differential equations with symmetric solutions

Pengyu Hou School of Mathematical Sciences, Beijing Normal University, Beijing 100875, China ([email protected]). Fang Liu School of Statistics and Mathematics, Central University of Finance and Economics, Beijing 102206, China ([email protected]). Aihui Zhou LSEC, Institute of Computational Mathematics and Scientific/Engineering Computing, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China ([email protected])

Abstract

In this paper, a symmetrized two-scale finite element method is proposed for a class of partial differential equations with symmetric solutions. With this method, the finite element approximation on a fine tensor product grid is reduced to the finite element approximations on a much coarse grid and a univariant fine grid. It is shown by both theory and numerics including electronic structure calculations that the resulting approximation still maintains an asymptotically optimal accuracy. Consequently the symmetrized two-scale finite element method reduces computational cost significantly.

Key words. symmetric, two-scale, finite element, partial differential equation.

AMS subject classifications. 65N15, 65N25, 65N30, 65N50.

1 Introduction

To improve the efficiency of solving multi-dimensional elliptic source and eigenvalue problems on tensor-product domains, the two-scale finite element method was proposed in [7, 12, 13]. With this method, low frequency parts of the finite element solution to an elliptic problem are captured on a coarse grid and high frequency components are handled by some univariant fine and coarse grids.

We understand that many quantities in science and engineering have symmetries such as wavefunctions, Hamiltonians, and interatomic potentials in quantum physics. It is shown in literature that efficient numerical methods can be designed to reduce the computational cost of solving problems with symmetric properties (see, e.g., [2, 6, 8] and references cited therein). For instance, a symmetry-based decomposition approach was proposed to solve eigenvalue problems with spatial symmetries in [6]. By this approach, the original differential eigenvalue problem is decomposed into some eigenvalue subproblems that can be carried out in parallel. In addition, each subproblem requires only a smaller number of eigenpairs compared with the original problem. As a result, the computational cost is reduced.

To reduce the computational cost of approximations of the elliptic source and eigenvalue problems with symmetric solutions, in this paper, we study some symmetrized two-scale finite element discretizations.

Let us give an informal description of the main ideas and results in this paper. Set $k$ be a positive integer. The symmetric group $\text{Sym}(k)$ is the set of bijections of $\{1,2,\cdots,k\}$ to itself. Specially, a permutation $\sigma\in\text{Sym}(k)$ which interchanges two letters $i$ and $j$ and leaves all the other letters unchanged is called a transposition. Denote the transposition $\sigma$ which interchanges $i$ and $j$ by $(i,j)$ [10].

Set $\Omega=(0,1)^{3}$ and $\mathbf{x}=(x_{1},x_{2},x_{3})$ . We say that $u$ is a symmetric function on $\bar{\Omega}$ [2] if

u(x_{\sigma(1)},x_{\sigma(2)},x_{\sigma(3)})=u(x_{1},x_{2},x_{3}),\quad\quad\forall\mathbf{x}\in\bar{\Omega},~{}\sigma\in\text{Sym}(3).

Let ${T}^{h_{1},h_{2},h_{3}}(\Omega)$ be a uniform tensor product grid on $\bar{\Omega}$ . $S^{h_{1},h_{2},h_{3}}(\Omega)$ is the associated tensor-product space of piecewise trilinear functions on $\bar{\Omega}$ . Let $u_{h_{1},h_{2},h_{3}}\in S^{h_{1},h_{2},h_{3}}(\Omega)$ be a finite element approximation to $u$ .

In computation, the grid points of $T^{h_{1},h_{2},h_{3}}(\Omega)$ are usually numbered along the $x_{1}$ -direction, the $x_{2}$ -direction, and the $x_{3}$ -direction, consecutively. Let $N_{i}=1/h_{i}(i=1,2,3)$ be a positive integer. Then the values of $u_{h_{1},h_{2},h_{3}}$ on these grid points are stored in a vector denoted by $U^{h_{1},h_{2},h_{3}}$ with the size $N_{dof}:=\prod_{i=1}^{3}(N_{i}+1)$ . Let $G:\mathbb{C}^{N_{dof}}\rightarrow S^{h_{1},h_{2},h_{3}}(\Omega)$ be a bijection satisfying

G(U^{h_{1},h_{2},h_{3}})=\sum_{i=1}^{N_{dof}}U^{h_{1},h_{2},h_{3}}_{i}\phi_{i},\quad\quad\forall U^{h_{1},h_{2},h_{3}}\in\mathbb{C}^{N_{dof}},

where $U^{h_{1},h_{2},h_{3}}_{i}$ is the $i$ -th component of $U^{h_{1},h_{2},h_{3}}$ and $\phi_{i}$ is the Lagrangian basis function of $S^{h_{1},h_{2},h_{3}}(\Omega)$ with $i=1,2,\cdots,N_{dof}$ . Obviously, we have $u_{h_{1},h_{2},h_{3}}=G(U^{h_{1},h_{2},h_{3}})$ .

For any $\sigma\in\text{Sym}(3)$ , there holds $u_{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}=G(U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}})$ . We define a vector transformation $\mathcal{T}_{\sigma}:\mathbb{C}^{N_{dof}}\rightarrow\mathbb{C}^{N_{dof}}$ satisfying

\mathcal{T}_{\sigma}(U^{h_{1},h_{2},h_{3}})=U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}},\quad\forall U^{h_{1},h_{2},h_{3}}\in\mathbb{C}^{N_{dof}}.

Then $U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}$ can be obtained from $U^{h_{1},h_{2},h_{3}}$ by a symmetrization operation (Algorithm 3.1). Let $R_{\sigma}:S^{h_{1},h_{2},h_{3}}(\Omega)\rightarrow S^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}(\Omega)$ be a bijection defined by

R_{\sigma}=G\circ\mathcal{T}_{\sigma}\circ G^{-1}.

(1.1)

We get

u_{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}=R_{\sigma}(u_{h_{1},h_{2},h_{3}}),\quad\forall u_{h_{1},h_{2},h_{3}}\in S^{h_{1},h_{2},h_{3}}(\Omega).

(1.2)

Note that the two-scale finite element solution [7, 12] of a linear elliptic source problem with the Dirichlet boundary condition is a linear combination of the standard finite element solutions on four different scale grids. That is,

B^{h}_{H,H,H}u=u_{h,H,H}+u_{H,h,H}+u_{H,H,h}-2u_{H,H,H}.

Hence, for a linear source problem with symmetric solution, we obtain from (1.2) that $u_{H,h,H}=R_{(1,2)}(u_{h,H,H})$ and $u_{H,H,h}=R_{(1,3)}(u_{h,H,H})$ . Consequently,

B^{h}_{H,H,H}u=u_{h,H,H}+R_{(1,2)}(u_{h,H,H})+R_{(1,3)}(u_{h,H,H})-2u_{H,H,H}.

It is seen from (1.1) that the amount of operations for getting $R_{(1,2)}(u_{h,H,H})$ or $R_{(1,3)}(u_{h,H,H})$ from $u_{h,H,H}$ is essentially that for calculating $\mathcal{T}_{(1,2)}(U^{h,H,H})$ or $\mathcal{T}_{(1,3)}(U^{h,H,H})$ from $U^{h,H,H}$ by the symmetrization algorithm (Algorithm 3.1), respectively. It is only $18N_{dof}$ with $N_{dof}:=(n-1)(N-1)(N-1)$ if $n=1/h$ and $N=1/H$ . While the amount of operations for computing $u_{H,h,H}$ or $u_{H,H,h}$ by the standard finite element method is far from $18N_{dof}$ , due to creating stiffness matrix and solving linear algebraic systems. Hence the symmetrization algorithm is quite efficient. Similarly, the symmetrization algorithm can be also combined with the postprocessed two-scale finite element method [11] and the two-scale postprocessed finite element method [15] respectively to construct new and more efficient algorithms.

This paper is structured as follows. In section 2, some preliminaries are presented and the existing related results that will be used in this paper are called. In section 3, a vector transformation operator associated with symmetric functions is constructed. A symmetrization algorithm is then proposed to perform this vector transformation operator. In section 4 and section 5, combining the symmetrization algorithm with the two-scale finite element method, some symmetrized two-scale finite element discretizations are developed for the linear source and eigenvalue problems with symmetric solutions, respectively. In section 6, several numerical results, including the applications to electronic structure calculations, are presented to illustrate the efficiency of our approaches. Finally, some concluding remarks are given in section 7.

2 Preliminaries

Let $\Omega=(0,1)^{d}$ with $d\geq 2$ . The standard notation for Sobolev spaces $W^{s,p}(\Omega)$ and their associated norms and seminorms will be used [4]. Set $H^{s}(\Omega)=W^{s,2}(\Omega)$ and $\lVert\ \cdot\ \rVert_{s,\Omega}=\lVert\ \cdot\ \rVert_{s,2,\Omega}$ . Let $H_{0}^{1}(\Omega)=\{v\in H^{1}(\Omega):v\mid_{\partial\Omega}=0\}$ , $H^{-1}(\Omega)$ be the dual of $H_{0}^{1}(\Omega)$ , and $(\cdot,\cdot)$ be the standard $L^{2}(\Omega)$ inner product. We shall use $A\lesssim B$ to mean that $A\leq CB$ for some constant $C$ which does not depend on any grid parameters.

Let $\mathbb{N}_{0}$ be the set of all nonnegative integers. For $w\in W^{s,p}(\Omega)$ , $\mathbf{x}=(x_{1},x_{2},\cdot\cdot\cdot,x_{d})\in\Omega$ , and $\mathbf{\alpha}=(\alpha_{1},\alpha_{2},\cdot\cdot\cdot,\alpha_{d})\in\mathbb{N}_{0}^{d}$ , set

(D^{\mathbf{\alpha}}w)(\mathbf{x})=(\frac{\partial^{\alpha_{1}}}{\partial x_{1}^{\alpha_{1}}}\cdot\cdot\cdot\frac{\partial^{\alpha_{d}}w}{\partial x_{d}^{\alpha_{d}}})(\mathbf{x})\quad\quad\mbox{and}\quad\quad\lvert\mathbf{\alpha}\rvert=\alpha_{1}+\alpha_{2}+\cdot\cdot\cdot+\alpha_{d}.

Furthermore, set $\textbf{0}=(0,\cdot\cdot\cdot,0)\in\mathbb{R}^{d}$ , $\textbf{e}=(1,\cdot\cdot\cdot,1)\in\mathbb{R}^{d}$ , and $\mathbb{Z}_{d}={1,2,\cdot\cdot\cdot,d}$ . For each $i\in\mathbb{Z}_{d}$ , let $\textbf{e}_{i}=(0,\cdot\cdot\cdot,0,1,0,\cdot\cdot\cdot,0)\in\mathbb{R}^{d}$ , which is the vector whose $i$ -th component is one and other components are zero. Write $\hat{\textbf{e}}_{i}=\textbf{e}-\textbf{e}_{i}$ .

We shall use the following Sobolev space

W^{G,4}(\Omega)=\{w\in H^{3}(\Omega):D^{\mathbf{\alpha}}w\in L^{2}(\Omega),\mathbf{0}\leq\mathbf{\alpha}\leq 3\mathbf{e},\lvert\mathbf{\alpha}\rvert=4\}

with its natural norm $\lVert\ \cdot\ \rVert_{W^{G,4}(\Omega)}$ [16].

2.1 A source problem

Consider the following elliptic source problem:

\displaystyle\left\{\begin{aligned} Lu&=f\ \ \mbox{in}\ \ \Omega,\\ u&=0\ \ \mbox{on}\ \ \partial\Omega.\end{aligned}\right.

(2.1)

Here $f\in H^{-1}(\Omega)$ and $L$ is defined by

Lu=-\sum_{i,j=1}^{d}\frac{\partial}{\partial x_{i}}\big{(}a_{ij}\frac{\partial u}{\partial x_{j}}\big{)}+\sum_{i=1}^{d}b_{i}\frac{\partial u}{\partial x_{i}}+Vu,

where $a_{ij}\in W^{1,\infty}(\Omega),\ b_{i},\ V\in L^{\infty}(\Omega)$ , and $\sum_{i,j=1}^{d}a_{ij}x_{i}x_{j}\geq c\sum_{k=1}^{d}x_{k}^{2}$ on $\bar{\Omega}$ for some constant $c>0$ .

A weak form of (2.1) is: find $u\in H_{0}^{1}(\Omega)$ such that

a(u,v)=(f,v)\ \ \ \ \forall v\in H_{0}^{1}(\Omega),

(2.2)

where

a(u,v)=\int_{\Omega}\Big{(}\sum_{i,j=1}^{d}a_{ij}\frac{\partial u}{\partial x_{i}}\frac{\partial v}{\partial x_{j}}+\sum_{i=1}^{d}b_{i}\frac{\partial u}{\partial x_{i}}v+Vuv\Big{)}\quad\mbox{and}\quad(f,v)=\int_{\Omega}fv.

(2.3)

For each $f\in H^{-1}(\Omega)$ , assume that there is a unique solution $u\in H_{0}^{1}(\Omega)$ of (2.2).

Set $N\in\mathbb{N}_{0}/\{0\}$ and $h=1/N\in[0,1]$ . Let $T^{h}[0,1]$ be the uniform grid with grid size $h$ . Let $h_{j}=1/N_{j}$ for $N_{j}\in\mathbb{N}_{0}/\{0\},\ j\in\mathbb{Z}_{d}$ . Set ${\bf h}=(h_{1},\ldots,h_{d})$ and $\displaystyle h=\max_{j\in\mathbb{Z}_{d}}h_{j}$ . Define a tensor-product grid on $\bar{\Omega}=[0,1]^{d}$ by

T^{\bf h}(\Omega):=T^{h_{1}}[0,1]\times\cdots\times T^{h_{d}}[0,1].

Let $S^{h}[0,1]\subset C[0,1]$ be the piecewise linear finite element space on $T^{h}[0,1]$ . The tensor product spaces of piecewise $d$ -linear functions on $\bar{\Omega}$ are

S^{\bf h}(\Omega):=S^{h_{1}}[0,1]\otimes\cdots\otimes S^{h_{d}}[0,1]\qquad\text{and}\qquad S_{0}^{\bf h}(\Omega):=S^{\bf h}(\Omega)\cap H_{0}^{1}(\Omega).

The standard Galerkin finite element method for solving (2.2) is: find $u_{\textbf{h}}\in S_{0}^{\textbf{h}}(\Omega)$ such that

a(u_{\textbf{h}},v)=(f,v)\ \ \ \forall\ v\in S_{0}^{\textbf{h}}(\Omega).

(2.4)

Define the Galerkin projector $P_{\bf h}:H^{1}_{0}(\Omega)\mapsto S^{\bf h}_{0}(\Omega)$ by

{a}(w-P_{\bf h}w,v)=0~{}\quad\forall v\in S^{\bf h}_{0}(\Omega),\quad\forall w\in H^{1}_{0}(\Omega).

(2.5)

Obviously $P_{\bf h}w=w_{\bf h},~{}\forall w\in H^{1}_{0}(\Omega)$ .

2.2 An eigenvalue problem

Assume that $a(\cdot,\cdot)$ is the same as (2.3) with $b_{i}=0,~{}i=1,2,\cdots,d$ and $(a_{ij})$ is symmetric. Consider the following eigenvalue problem

\displaystyle a(u,v)=\lambda(u,v)\quad\forall v\in H^{1}_{0}(\Omega),

(2.6)

where $\lambda$ is the eigenvalue of the symmetric bilinear form $a(\cdot,\cdot)$ relative to the $L^{2}(\Omega)$ inner product $(\cdot,\cdot)$ and $u$ is the corresponding eigenvector. The eigenvalues of (2.6) are real numbers satisfying $\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\cdots$ [1]. Assume that the corresponding eigenvectors $u_{1},u_{2},u_{3},\dots$ satisfy $(u_{i},u_{j})=\delta_{ij}$ for $i,j=1,2,\cdots$ with $\delta_{ij}$ the Dirac delta.

Consider the standard Galerkin finite element method for solving (2.6): find a pair $(\lambda_{\bf h},u_{\bf h})\in\mathbb{R}\times S^{\bf h}_{0}(\Omega)$ satisfying

\displaystyle a(u_{\bf h},v)=\lambda_{\bf h}(u_{\bf h},v)~{}~{}~{}\forall v\in S^{\bf h}_{0}(\Omega),

(2.7)

which has finite eigenvalues $\lambda_{1,{\bf h}}<\lambda_{2,{\bf h}}\leq\cdots\leq\lambda_{n_{\bf h},{\bf h}},$ where $n_{\bf h}=\mbox{dim}~{}S^{\bf h}_{0}(\Omega)$ . Assume that the corresponding eigenvectors $u_{1,{\bf h}},u_{2,{\bf h}},\dots,u_{n_{\bf h},{\bf h}}$ satisfy $(u_{i,{\bf h}},$ $u_{j,{\bf h}})=\delta_{ij}$ . The Rayleigh principle [1] implies that $\lambda_{i}\leq\lambda_{i,{\bf h}}\quad\text{for }i=1,2,\dots,n_{\bf h}.$

Suppose that $h=\displaystyle\max_{j\in\mathbb{Z}_{d}}h_{j}\ll 1,\lambda$ is simple, and $(\lambda_{\bf h},u_{\bf h})$ is an approximation to $(\lambda,u)$ which satisfies (2.7) and Lemma 3.2 in [11] with $\|u_{\bf h}\|_{0,\Omega}=1$ . There holds [1]

\lambda_{\bf h}-\lambda+\|u-u_{\bf h}\|_{0,\Omega}+h\|u-u_{\bf h}\|_{1,\Omega}\lesssim h^{2}.

(2.8)

Lemma 2.1.

[7] There holds

\lambda_{\bf h}-\lambda=\lambda(u,u-P_{\bf h}u)+\mathcal{O}(h^{4}).

(2.9)

3 A vector transformation operator and its implementation

In this section, we will introduce a vector transformation operator associated with symmetric functions. A so-called symmetrization algorithm will be proposed to implement this vector transformation. We will combine the symmetrization algorithm with the two-scale finite element method to design new and efficient algorithms in Section 4.

We call $u$ a symmetric function on $\bar{\Omega}$ if

u(x_{\sigma(1)},x_{\sigma(2)},\cdots,x_{\sigma(d)})=u(x_{1},x_{2},\cdots,x_{d}),\quad\quad\forall\mathbf{x}\in\bar{\Omega},~{}\sigma\in\text{Sym}(d).

(3.1)

Thus we have

	$\displaystyle u(i_{\sigma(1)}h_{\sigma(1)},i_{\sigma(2)}h_{\sigma(2)},\cdots,i_{\sigma(d)}h_{\sigma(d)})$	$\displaystyle=u(i_{1}h_{1},i_{2}h_{2},\cdots,i_{d}h_{d}),$
		$\displaystyle\quad\quad i_{k}=0,1,\cdots,N_{k};k=1,2,\cdots,d.$		(3.2)

Next we illustrate how to get the values of $u(\mathbf{x})$ on the grid points of $T^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}(\Omega)$ from those on the grid points of $T^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ in computation. The values of $u(\mathbf{x})$ on the grid points of $T^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ are stored in a vector $U^{h_{1},h_{2},\cdots,h_{d}}$ with the size $N_{dof}:=\prod_{i=1}^{d}(N_{i}+1)$ . Usually the grid points of $T^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ are numbered along the $x_{1}$ -direction, the $x_{2}$ -direction, $\cdots$ , and the $x_{d}$ -direction, consecutively. That is, the $\hat{I}$ -th component of $U^{h_{1},h_{2},\cdots,h_{d}}$ satisfies

	$\displaystyle U^{h_{1},h_{2},\cdots,h_{d}}_{\hat{I}}$	$\displaystyle=u(i_{1}h_{1},i_{2}h_{2},\cdots,i_{d}h_{d}),$
		$\displaystyle\quad\quad i_{k}=0,1,\cdots,N_{k};k=1,2,\cdots,d,$		(3.3)

where the subscript

\hat{I}:=I(i_{1},i_{2},\cdots,i_{d})=\sum_{l=2}^{d}\left(i_{l}\prod_{j=1}^{l-1}(N_{j}+1)\right)+i_{1}+1.

(3.4)

The values of $u(\mathbf{x})$ on the grid points of $T^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}(\Omega)$ are contained in the vector $U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}$ . We obtain from (3)-(3.4) that

	$\displaystyle U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}_{I(i_{\sigma(1)},i_{\sigma(2)},\cdots,i_{\sigma(d)})}$	$\displaystyle=U^{h_{1},h_{2},\cdots,h_{d}}_{I(i_{1},i_{2},\cdots,i_{d})},$
		$\displaystyle\quad\quad i_{k}=0,1,\cdots,N_{k};k=1,2,\cdots,d,$		(3.5)

with

I(i_{\sigma(1)},i_{\sigma(2)},\cdots,i_{\sigma(d)})=\sum_{l=2}^{d}\left(i_{\sigma(l)}\prod_{j=1}^{l-1}(N_{\sigma(j)}+1)\right)+i_{\sigma(1)}+1.

Hence for any $\sigma\in\text{Sym}(d)$ , one can get $U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}$ from $U^{h_{1},h_{2},\cdots,h_{d}}$ . We can define a vector transformation operator $\mathcal{T}_{\sigma}:\mathbb{C}^{N_{dof}}\rightarrow\mathbb{C}^{N_{dof}}$ satisfying

\mathcal{T}_{\sigma}(U^{h_{1},h_{2},\cdots,h_{d}})=U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}},\quad\quad\forall U^{h_{1},h_{2},\cdots,h_{d}}\in\mathbb{C}^{N_{dof}}.

(3.6)

This vector transformation operator can be implemented by Algorithm 3.1, which is called the symmetrization algorithm.

Algorithm 3.1.

Input: $h_{1}$ , $h_{2}$ , $\cdots$ , $h_{d}$ , $U^{h_{1},h_{2},\cdots,h_{d}}$ , $\sigma$ .
Output: $U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}$ .

N_{1}=1/h_{1};\ N_{2}=1/h_{2};\cdots;\ N_{d}=1/h_{d}

2:for

i_{1}=0:N_{1}

3: for

i_{2}=0:N_{2}

\cdots

5: for

i_{d}=0:N_{d}

U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}_{I(i_{\sigma(1)},i_{\sigma(2)},\cdots,i_{\sigma(d)})}=U^{h_{1},h_{2},\cdots,h_{d}}_{I(i_{1},i_{2},\cdots,i_{d})}

7:return

U^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}

Calculating $I(i_{1},i_{2},\cdots,i_{d})$ in (3.4) requires $d^{2}$ operations. Hence the amount of operations for Algorithm 3.1 is $2d^{2}N_{dof}$ . Next we illustrate the idea of Algorithm 3.1 further for $d=2,3$ , respectively.

3.1 Two dimensional case

Let $\Omega=(0,1)^{2}$ . We see that there is only one element $\sigma=(1,2)$ in $\text{Sym}(2)$ besides the identical operator. Assume that $u$ is a symmetric function on $\bar{\Omega}$ . Namely,

u(x_{2},x_{1})=u(x_{\sigma(1)},x_{\sigma(2)})=u(x_{1},x_{2}),\ \ \forall\mathbf{x}\in\bar{\Omega}.

(3.7)

Consequently for the grid points of $T^{h_{2},h_{1}}(\Omega)=T^{h_{\sigma(1)},h_{\sigma(2)}}(\Omega)$ and $T^{h_{1},h_{2}}(\Omega)$ , we have

	$\displaystyle u(i_{2}h_{2},i_{1}h_{1})=u(i_{\sigma(1)}h_{\sigma(1)},i_{\sigma(2)}h_{\sigma(2)})$	$\displaystyle=u(i_{1}h_{1},i_{2}h_{2}),$
		$\displaystyle i_{k}=0,1,\cdots,N_{k};k=1,2.$		(3.8)

In implementation, the values of $u(x_{1},x_{2})$ on the grid points of a uniform grid $T^{h_{1},h_{2}}(\Omega)$ are stored in a vector $U^{h_{1},h_{2}}$ with the size $N_{dof}:=(N_{1}+1)\times(N_{2}+1)$ . Usually the grid points of a uniform grid $T^{h_{1},h_{2}}(\Omega)$ are numbered along the $x_{1}$ -direction and the $x_{2}$ -direction, consecutively. That is, the $\hat{I}$ -th component of $U^{h_{1},h_{2}}$ satisfies

U^{h_{1},h_{2}}_{\hat{I}}=u(i_{1}h_{1},i_{2}h_{2}),\quad\quad i_{k}=0,1,\cdots,N_{k};k=1,2,

(3.9)

where the subscript

\hat{I}:=I(i_{1},i_{2})={i_{2}(N_{1}+1)+i_{1}+1}.

(3.10)

The values of $u(x_{2},x_{1})$ on the grid points of $T^{h_{2},h_{1}}(\Omega)$ are stored in the vector $U^{h_{2},h_{1}}$ . By (3.1) – (3.10), we have

U^{h_{2},h_{1}}_{I(i_{2},i_{1})}=U^{h_{1},h_{2}}_{I(i_{1},i_{2})},\quad\quad i_{k}=0,1,\cdots,N_{k};k=1,2,

(3.11)

with

I(i_{2},i_{1})={i_{1}(N_{2}+1)+i_{2}+1}.

That is, one can get $U^{h_{2},h_{1}}$ from $U^{h_{1},h_{2}}$ . Let $d=2$ in (3.6). It follows that $\mathcal{T}_{(1,2)}(U^{h_{1},h_{2}})=U^{h_{2},h_{1}}$ .

Refer to caption — Figure 1: The process of obtaining the vector $U^{1,\frac{1}{2}}$ from $U^{\frac{1}{2},1}$ . For example, $U^{1,\frac{1}{2}}_{2}=U^{1,\frac{1}{2}}_{I(1,0)}=u(1,0)=u(0,1)=U^{\frac{1}{2},1}_{I(0,1)}=U_{4}^{\frac{1}{2},1}$ by (3.9), (3.10), and (3.11).

For example, if $h_{1}=\frac{1}{2}$ and $h_{2}=1$ , it is shown in Figure 1 that the values of $u(x_{1},x_{2})$ on the grid points of the uniform grids $T^{\frac{1}{2},1}(\Omega)$ and $T^{1,\frac{1}{2}}(\Omega)$ are stored in the vectors $U^{\frac{1}{2},1}$ and $U^{1,\frac{1}{2}}$ , respectively. By Algorithm 3.1, $U^{1,\frac{1}{2}}=\mathcal{T}_{(1,2)}(U^{\frac{1}{2},1})$ can be obtained from $U^{\frac{1}{2},1}$ . Similarly, if $h_{1}=\frac{1}{4}$ and $h_{2}=\frac{1}{2}$ , The process of obtaining the vector $U^{\frac{1}{2},\frac{1}{4}}$ from $U^{\frac{1}{4},\frac{1}{2}}$ is illustrated in Figure 2, in which $\hat{U}:=U^{\frac{1}{4},\frac{1}{2}}$ and $\tilde{U}:=U^{\frac{1}{2},\frac{1}{4}}$ for simplicity.

The process of obtaining the vectors $U^{h_{\sigma(1)},h_{\sigma(2)}}=\mathcal{T}_{\sigma}(U^{h_{1},h_{2}})$ from $U^{h_{1},h_{2}}$ can be implemented by Algorithm 3.1 with $d=2$ , which demands $8N_{dof}$ operations.

3.2 Three dimensional case

Let $\Omega=(0,1)^{3}$ . Assume that $u$ is a symmetric function on $\bar{\Omega}$ . That is,

u(x_{\sigma(1)},x_{\sigma(2)},x_{\sigma(3)})=u(x_{1},x_{2},x_{3}),\quad\forall\mathbf{x}\in\bar{\Omega},~{}\sigma\in\text{Sym}(3).

Thus for the grid points of $T^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}(\Omega)$ and $T^{h_{1},h_{2},h_{3}}(\Omega)$ , we have

	$\displaystyle u(i_{\sigma(1)}h_{\sigma(1)},i_{\sigma(2)}h_{\sigma(2)},i_{\sigma(3)}h_{\sigma(3)})$	$\displaystyle=u(i_{1}h_{1},i_{2}h_{2},i_{3}h_{3}),$
		$\displaystyle i_{k}=0,1,\cdots,N_{k};k=1,2,3.$		(3.12)

In implementation, similarly, the values of $u(x_{1},x_{2},x_{3})$ on the grid points of $T^{h_{1},h_{2},h_{3}}(\Omega)$ are stored in a vector denoted by $U^{h_{1},h_{2},h_{3}}$ with the size $N_{dof}:=\prod_{i=1}^{3}(N_{i}+1)$ . The grid points of $T^{h_{1},h_{2},h_{3}}(\Omega)$ are usually numbered along the $x_{1}$ -direction, the $x_{2}$ -direction, and the $x_{3}$ -direction, consecutively. Thus the $\hat{I}$ -th component of $U^{h_{1},h_{2},h_{3}}$ satisfies

	$\displaystyle U^{h_{1},h_{2},h_{3}}_{\hat{I}}$	$\displaystyle=u(i_{1}h_{1},i_{2}h_{2},i_{3}h_{3}),$
		$\displaystyle i_{k}=0,1,\cdots,N_{k};k=1,2,3,$		(3.13)

where the subscript

\hat{I}:=I(i_{1},i_{2},i_{3})=i_{3}(N_{1}+1)(N_{2}+1)+i_{2}(N_{1}+1)+i_{1}+1.

(3.14)

The values of $u(\mathbf{x})$ on the grid points of $T^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}(\Omega)$ are contained in the vector $U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}$ . By (3.2), (3.2), and (3.14), we obtain

	$\displaystyle U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}_{I(i_{\sigma(1)},i_{\sigma(2)},i_{\sigma(3)})}$	$\displaystyle=U^{h_{1},h_{2},h_{3}}_{I(i_{1},i_{2},i_{3})},$
		$\displaystyle i_{k}=0,1,\cdots,N_{k};k=1,2,3,$		(3.15)

with

I(i_{\sigma(1)},i_{\sigma(2)},i_{\sigma(3)})=i_{\sigma(3)}(N_{\sigma(1)}+1)(N_{\sigma(2)}+1)+i_{\sigma(2)}(N_{\sigma(1)}+1)+i_{\sigma(1)}+1.

That is, the vector $U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}$ can be obtained from $U^{h_{1},h_{2},h_{3}}$ . Let $d=3$ in (3.6). It follows that $\mathcal{T}_{\sigma}(U^{h_{1},h_{2},h_{3}})=U^{h_{\sigma(1)},h_{\sigma(2)},h_{\sigma(3)}}$ .

For example, if $h_{1}=\frac{1}{4}$ , $h_{2}=\frac{1}{2}$ , and $h_{3}=\frac{1}{2}$ , the components of the vector $U^{\frac{1}{4},\frac{1}{2},\frac{1}{2}}$ are illustrated in Figures 3 and 4, where $\hat{U}:=U^{\frac{1}{4},\frac{1}{2},\frac{1}{2}}$ for simplicity.

The vectors $U^{\frac{1}{2},\frac{1}{4},\frac{1}{2}}=\mathcal{T}_{(1,2)}(U^{\frac{1}{4},\frac{1}{2},\frac{1}{2}})$ and $U^{\frac{1}{2},\frac{1}{2},\frac{1}{4}}=\mathcal{T}_{(1,3)}(U^{\frac{1}{4},\frac{1}{2},\frac{1}{2}})$ can be obtained from $U^{\frac{1}{4},\frac{1}{2},\frac{1}{2}}$ with $\sigma=(1,2)$ and $(1,3)$ in (3.2), respectively. It is illustrated in Figures 5 and 6.

The process of obtaining the vectors $U^{h_{2},h_{1},h_{3}}=U^{h_{(1,2)(1)},h_{(1,2)(2)},h_{(1,2)(3)}}$ or $U^{h_{3},h_{2},h_{1}}=U^{h_{(1,3)(1)},h_{(1,3)(2)},h_{(1,3)(3)}}$ from $U^{h_{1},h_{2},h_{3}}$ can be implemented by Algorithm 3.1 with $d=3$ , which requires $18N_{dof}$ operations.

4 The symmetrized two-scale finite element approximations

For the source problem (2.2) and eigenvalue problem (2.6), the two-scale finite element method has been developed [7, 12]. In this section, for (2.2) and (2.6) with symmetric solutions, Algorithm 3.1 will be combined with the two-scale finite element method.

Let $v$ be a symmetric function and $v_{h_{1},h_{2},\cdots,h_{d}}\in{S}^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ be a standard finite element approximation to $v$ . As stated in Section 3, usually the grid points of $T^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ are numbered along the $x_{1}$ -direction, the $x_{2}$ -direction, $\cdots$ , and the $x_{d}$ -direction, consecutively. The values of $v_{h_{1},h_{2},\cdots,h_{d}}$ on these nodes are stored in a vector denoted by $V^{h_{1},h_{2},\cdots,h_{d}}$ with the size $N_{dof}:=\prod_{i=1}^{d}(N_{i}+1)$ .

Define a bijection $G:\mathbb{C}^{N_{dof}}\rightarrow S^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ satisfying

G(V^{h_{1},h_{2},\cdots,h_{d}})=\sum_{i=1}^{N_{dof}}V^{h_{1},h_{2},\cdots,h_{d}}_{i}\phi_{i},\quad\quad\forall V^{h_{1},h_{2},\cdots,h_{d}}\in\mathbb{C}^{N_{dof}},

where $V^{h_{1},h_{2},\cdots,h_{d}}_{i}$ is the $i$ -th component of $V^{h_{1},h_{2},\cdots,h_{d}}$ and $\phi_{i}$ is the Lagrangian basis of $S^{h_{1},h_{2},\cdots,h_{d}}(\Omega)$ with $i=1,2,\cdots,N_{dof}$ . Obviously we have $v_{h_{1},h_{2},\cdots,h_{d}}=G(V^{h_{1},h_{2},\cdots,h_{d}})$ . For any $\sigma\in\text{Sym}(d)$ , there holds

v_{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}=G(V^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}).

Invoking (3.6) we obtain

\mathcal{T}_{\sigma}(V^{h_{1},h_{2},\cdots,h_{d}})=V^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}},\quad\forall V^{h_{1},h_{2},\cdots,h_{d}}\in\mathbb{C}^{N_{dof}}.

The process of obtaining $V^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}$ from $V^{h_{1},h_{2},\cdots,h_{d}}$ can be implemented by a symmetrization algorithm (Algorithm 3.1). Define a bijection $R_{\sigma}:S^{h_{1},h_{2},\cdots,h_{d}}(\Omega)\rightarrow S^{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}(\Omega)$ satisfying

R_{\sigma}=G\circ\mathcal{T}_{\sigma}\circ G^{-1}.

(4.1)

We have

v_{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}=R_{\sigma}(v_{h_{1},h_{2},\cdots,h_{d}}),\quad\forall v_{h_{1},h_{2},\cdots,h_{d}}\in S^{h_{1},h_{2},\cdots,h_{d}}(\Omega).

(4.2)

It is seen from (4.1) and (4.2) that the critical process of getting $R_{\sigma}(v_{h_{1},h_{2},\cdots,h_{d}})$ from $v_{h_{1},h_{2},\cdots,h_{d}}$ is to calculate $\mathcal{T}_{\sigma}(V^{h_{1},h_{2},\cdots,h_{d}})$ from $V^{h_{1},h_{2},\cdots,h_{d}}$ by Algorithm 3.1, which requires $2d^{2}N_{dof}$ operations with $N_{dof}:=\prod_{i=1}^{d}(N_{i}+1)$ . While computing $v_{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}$ by the standard finite element method usually demands $MN_{dof}$ operations with an amount $M\gg 2d^{2}$ , due to creating stiffness matrix and solving linear algebraic systems. Consequently, if $v_{h_{1},h_{2},\cdots,h_{d}}$ has been calculated, then obtaining $v_{h_{\sigma(1)},h_{\sigma(2)},\cdots,h_{\sigma(d)}}$ through Algorithm 3.1 is much more efficient than the standard finite element method. From this observation, some efficient algorithms for (2.2) and (2.6) with symmetric solutions can be established.

4.1 The source problem

For the source problem (2.2) with symmetric solution, we propose a symmetrized two-scale finite element method to reduce computational cost by combining Algorithm 3.1 with the two-scale finite element method in [7, 12]. Let $h,H\in(0,1)$ and assume that $H/h$ is a positive integer. Let $w_{h\alpha+H(\textbf{e}-\alpha)}\in S_{0}^{h\alpha+H(\textbf{e}-\alpha)}(\Omega)(\textbf{0}\leq\alpha\leq\textbf{e})$ , and define

B^{h}_{H\textbf{e}}w_{h\textbf{e}}=\sum_{i=1}^{d}w_{H\hat{\textbf{e}}_{i}+h\textbf{e}_{i}}-(d-1)w_{H\textbf{e}}.

For instance, $B^{h}_{H,H,H}w_{h,h,h}=w_{h,H,H}+w_{H,h,H}+w_{H,H,h}-2w_{H,H,H}$ , for $d=3$ . Invoking (4.2), a symmetrized two-scale finite element method for the linear source problem (2.2) with symmetric solution is described in Algorithm 4.2.

Algorithm 4.2.

1.

Solve (2.4) on a coarse grid: find $P_{H\textbf{e}}u\in S_{0}^{H\textbf{e}}(\Omega)$ such that

$a(P_{H\textbf{e}}u,v)=(f,v)\ \ \ \ \forall v\in S_{0}^{H\textbf{e}}(\Omega).$

Solve (2.4) on a grid which is fine in the first coordinate direction: find $P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u\in S_{0}^{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}(\Omega)$ such that

\displaystyle a(P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u,v)=(f,v)\ \ \ \ \forall v\in S_{0}^{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}(\Omega).

Get $R_{(1,i)}(P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u)\in S_{0}^{h\textbf{e}_{i}+H\hat{\textbf{e}}_{i}}(\Omega),\ i\in\mathbb{Z}_{d}/\{1\}$ from $P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u$ through Algorithm 3.1.

Set

\displaystyle B^{h}_{H\textbf{e}}P_{h\textbf{e}}u=P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u+\sum_{i=2}^{d}R_{(1,i)}(P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u)-(d-1)P_{H\textbf{e}}u.

For $d=2$ or $3$ , there holds the following error estimates for the symmetrized two-scale finite element solution.

Theorem 4.1.

Assume that $\partial_{x_{l}}a_{ij}\in W^{1,\infty}(\Omega)$ and $\partial_{x_{l}}b_{l}\in L^{\infty}(\Omega)$ for $i,j,l\in\mathbb{Z}_{d}$ . If $u\in H_{0}^{1}(\Omega)\cap W^{G,4}(\Omega)$ and $B^{h}_{H\textbf{e}}P_{h\textbf{e}}u$ is obtained from Algorithm 4.2, then

	$\displaystyle\lVert u-B^{h}_{H\textbf{e}}P_{h\textbf{e}}u\rVert_{1,\Omega}$	$\displaystyle\lesssim(h+H^{3})\lVert u\rVert_{W^{G,4}(\Omega)},$
	$\displaystyle\lVert u-B^{h}_{H\textbf{e}}P_{h\textbf{e}}u\rVert_{0,\Omega}$	$\displaystyle\lesssim(h^{2}+H^{4})\lVert u\rVert_{W^{G,4}(\Omega)}.$

Proof.

Invoking Theorem 3.3 and Lemma 4.1 in [15], we have

\lVert P_{h\textbf{e}}u-B^{h}_{H\textbf{e}}P_{h\textbf{e}}u\rVert_{1,\Omega}\lesssim(h^{2}+H^{3})\lVert u\rVert_{W^{G,4}(\Omega)}.

(4.3)

Then imitating the proof of Theorem 2.5 in [11], we complete the proof. ∎

Remark 4.1.

We may see from Algorithm 4.2 that the symmetrized two-scale finite element method is quite applicable for the source problem (2.1) with symmetric solution. $R_{(1,i)}(P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u)$ can be obtained from $P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u$ through Algorithm 3.1. The computational cost of Algorithm 3.1 is far smaller than computing $P_{h\textbf{e}_{i}+H\hat{\textbf{e}}_{i}}u$ by the standard finite element method. For example, when $d=3$ , the amount of operations for Algorithm 3.1 to get $P_{H,h,H}u=R_{(1,2)}(P_{h,H,H}u)$ from $P_{h,H,H}u$ is only $18N_{dof}$ with $N_{dof}:=(n-1)(N-1)(N-1)$ . While computing $P_{H,h,H}u$ by the standard finite element method usually requires $MN_{dof}$ operations with $M\gg 18$ , due to creating stiffness matrix and solving linear algebraic systems. Hence the computational cost is reduced significantly.

Remark 4.2.

Combining the two-scale finite element method with the postprocessing technique, the postprocessed two-scale finite element method and the two-scale postprocessed finite element method have been proposed in [11, 15], in which the interpolation postprocessing operator $\Pi_{2\textbf{h}}$ for uniform tensor product grids $T^{\bf h}(\Omega)$ is used. For the linear source problem (2.1) with symmetric solution, these methods can also be combined with Algorithm 3.1 similarly to reduce computational cost further.

4.2 The eigenvalue problem

Similar to Section 4.1, for the eigenvalue problem (2.6) with symmetric solution, we combine Algorithm 3.1 with the two-scale finite element method in [7, 12] to propose the symmetrized two-scale finite element method (Algorithm 4.3).

Algorithm 4.3.

1.

Solve (2.7) on a coarse grid: find $(\lambda_{H\textbf{e}},u_{H\textbf{e}})\in\mathbb{R}\times S_{0}^{H\textbf{e}}(\Omega)$ satisfying $\lVert u_{H\textbf{e}}\rVert_{0,\Omega}=1$ and

$a(u_{H\textbf{e}},v)=\lambda_{H\textbf{e}}(u_{H\textbf{e}},v)\ \ \ \ \forall v\in S_{0}^{H\textbf{e}}(\Omega).$

Solve (2.7) on a grid which is fine in the first coordinate direction:

find $(\lambda_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}},u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}})\in\mathbb{R}\times S_{0}^{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}(\Omega)$ satisfying $\lVert u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}\rVert_{0,\Omega}=1$ and

\displaystyle a(u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}},v)=\lambda_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}(u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}},v)\ \ \ \ \forall v\in S_{0}^{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}(\Omega).

For $i\in\mathbb{Z}_{d}/\{1\}$ , set $\lambda_{h\textbf{e}_{i}+H\hat{\textbf{e}}_{i}}=\lambda_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}$ and obtain $R_{(1,i)}(u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}})\in S_{0}^{h\textbf{e}_{i}+H\hat{\textbf{e}}_{i}}(\Omega)$ from $u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}\in S_{0}^{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}(\Omega)$ through Algorithm 3.1.

Set

	$\displaystyle B^{h}_{H\textbf{e}}\lambda_{h\textbf{e}}$	$\displaystyle=\sum_{i=1}^{d}\lambda_{h\textbf{e}_{i}+H\hat{\textbf{e}}_{i}}-(d-1)\lambda_{H\textbf{e}},$
	$\displaystyle B^{h}_{H\textbf{e}}u_{h\textbf{e}}$	$\displaystyle=u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}+\sum_{i=2}^{d}R_{(1,i)}(u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}})-(d-1)u_{H\textbf{e}}.$

From (4.2) we have $u_{h\textbf{e}_{i}+H\hat{\textbf{e}}_{i}}=R_{(1,i)}(u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}})$ . For $d=2$ or $3$ , there holds the following result, which is a refinement of the results in [7, 12].

Theorem 4.2.

Let $(\lambda_{H\textbf{e}},u_{H\textbf{e}})$ and $(\lambda_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}},u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}})$ be eigenpairs of (2.7) that approximate the same exact eigenpair $(\lambda,u)$ . If $\partial_{x_{l}}a_{ij}\in W^{1,\infty}(\Omega)$ for $i,j,l\in\mathbb{Z}_{d}$ and $u\in H_{0}^{1}(\Omega)\cap W^{G,4}(\Omega)$ , then

	$\displaystyle\lVert u-B^{h}_{H\textbf{e}}u_{h\textbf{e}}\rVert_{1,\Omega}$	$\displaystyle\lesssim h+H^{3},$		(4.4)
	$\displaystyle\lvert\lambda-B^{h}_{H\textbf{e}}\lambda_{h\textbf{e}}\rvert$	$\displaystyle\lesssim h^{2}+H^{4}.$		(4.5)

Proof.

Combining Lemma 5.2 with Theorem 3.3 in [15], we obtain

\lVert B^{h}_{H\textbf{e}}u_{h\textbf{e}}-u_{h\textbf{e}}\rVert_{1,\Omega}\lesssim h^{2}+H^{3}.

(4.6)

Thus we arrive at (4.4). By Lemma 2.1, we have

\displaystyle B^{h}_{H{\bf e}}\lambda_{h{\bf e}}-\lambda_{h{\bf e}}=\lambda(u,P_{h{\bf e}}u-B^{h}_{H{\bf e}}P_{h{\bf e}}u)+\mathcal{O}(H^{4}).

Applying Theorem 4.1, we complete the proof. ∎

Remark 4.3.

For the eigenvalue problem with symmetric eigenfunction, Algorithm 3.1 is combined with the two-scale finite element method to obtain new and more efficient algorithm (Algorithm 4.3). Algorithm 3.1 can also be combined with the postprocessed two-scale finite element method in [11] and the two-scale postprocessed finite element method in [15]. Besides, some two-scale finite element discretizations for the nonlinear eigenvalue problems have been proposed and analyzed in [9]. Consequently, for the nonlinear eigenvalue problems with symmetric eigenfunctions, the two-scale finite element methods in [9] can also be combined with Algorithm 3.1 to get new algorithms.

5 Numerical examples

In this section, several numerical examples, including the electronic structure calculations, are presented to illustrate the efficiency of our approaches. To optimize the cost of the computation, we choose $h=H^{2}$ approximately.

Example 1 Consider a source problem with a singular coefficient:

\displaystyle\vspace{-0.8cm}\left\{\begin{aligned} -\Delta u-\frac{1}{\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}}u+x_{1}x_{2}x_{3}u&=f\ \ in\ \ \Omega=(0,1)^{3},\\ u&=0\ \ on\ \ \partial\Omega\end{aligned}\right.\vspace{-0.8cm}

with $f$ chosen so that the exact solution is

u=2x_{1}x_{2}x_{3}(1-x_{1})(1-x_{2})(1-x_{3})\text{e}^{x_{1}+x_{2}+x_{3}}.

Let $L=1.0$ be the side length of $\Omega$ .

For this linear source problem with symmetric solution, the numerical results are presented in Figures 7-8 and Tables 1-2. Figure 7 supports the convergence results of Theorem 4.1. Time consumptions of different methods are presented in Figure 8. It is shown that the two-scale finite element method is more efficient than the standard finite element method. Algorithm 4.2 reduces the time consumption further compared with the two-scale finite element method, which is illustrated in more details in Tables 1-2. All of these results show that for this linear source problem with symmetric solution, considering both accuracy and computational cost, the symmetrized two-scale finite element method (Algorithm 4.2) is the preferred method.

Table 1: Example 1: The running time (s) by the two-scale FEM.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$P_{H\textbf{e}}u$	$P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u$	$P_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}}u$	$P_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}u$	$B_{H\textbf{e}}^{h}P_{h\textbf{e}}u$	Total
$16\times 4\times 4$	0.0030	0.0104	0.0110	0.0104	0.0006	0.0354
$36\times 6\times 6$	0.0106	0.0609	0.0624	0.0614	0.0052	0.2005
$64\times 8\times 8$	0.0364	0.2250	0.2220	0.2192	0.0281	0.7307
$100\times 10\times 10$	0.0754	0.5574	0.5484	0.5596	0.1104	1.8512
$144\times 12\times 12$	0.1005	1.2790	1.2835	1.2941	0.3264	4.2835

Table 2: Example 1: The running time (s) by Algorithm 4.2.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$P_{H\textbf{e}}u$	$P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u$	$P_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}}u,P_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}u$	$B_{H\textbf{e}}^{h}P_{h\textbf{e}}u$	Total
$16\times 4\times 4$	0.0030	0.0104	0.0001	0.0003	0.0138
$36\times 6\times 6$	0.0106	0.0609	0.0004	0.0026	0.0745
$64\times 8\times 8$	0.0364	0.2250	0.0027	0.0135	0.2776
$100\times 10\times 10$	0.0754	0.5574	0.0136	0.0581	0.7045
$144\times 12\times 12$	0.1005	1.2790	0.0461	0.1708	1.5964

Example 2 Consider a source problem:

\displaystyle\left\{\begin{aligned} -\sum\limits_{i,j=1}^{3}\frac{\partial}{\partial x_{j}}(a_{ij}\frac{\partial u}{\partial x_{i}})+\sum\limits_{i=1}^{3}b_{i}\frac{\partial u}{\partial x_{i}}+u&=f\ \ in\ \ \Omega=(1,2)^{3},\\ u&=0\ \ on\ \ \partial\Omega,\end{aligned}\right.

where

A=(a_{ij})=\begin{pmatrix}\text{e}^{x_{1}}&1&1\\ 1&\text{e}^{x_{2}}&1\\ 1&1&\text{e}^{x_{3}}\end{pmatrix},b_{1}=b_{2}=b_{3}=0.001,

and $f$ is chosen so that we have the symmetric exact solution as follows

u=(1-x_{1})(2-x_{1})(1-x_{2})(2-x_{2})(1-x_{3})(2-x_{3})(\sin\sqrt{x_{1}x_{2}x_{3}})\text{e}^{x_{1}+x_{2}+x_{3}}.

Let $L=1.0$ be the side length of $\Omega$ .

The numerical results are presented in Figures 9-10 and Tables 3-4 for the nonsymmetric source problem with symmetric solution. Similar to Example 1, Figure 9 also supports the convergence results of Theorem 4.1. It is shown that the symmetrized two-scale finite element method (Algorithm 4.2) is still better than the other methods.

Table 3: Example 2: The running time (s) by the two-scale FEM.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$P_{H\textbf{e}}u$	$P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u$	$P_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}}u$	$P_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}u$	$B_{H\textbf{e}}^{h}P_{h\textbf{e}}u$	Total
$16\times 4\times 4$	0.0062	0.0257	0.0243	0.0236	0.0008	0.0806
$36\times 6\times 6$	0.0292	0.1510	0.1539	0.1506	0.0052	0.4899
$64\times 8\times 8$	0.0677	0.5428	0.5385	0.5349	0.0275	1.7114
$100\times 10\times 10$	0.1424	1.4185	1.4174	1.4132	0.1095	4.5010
$144\times 12\times 12$	0.2577	3.1907	3.1658	3.1772	0.3328	10.1242

Table 4: Example 2: The running time (s) by Algorithm 4.2.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$P_{H\textbf{e}}u$	$P_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}u$	$P_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}}u,P_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}u$	$B_{H\textbf{e}}^{h}P_{h\textbf{e}}u$	Total
$16\times 4\times 4$	0.0062	0.0257	0.0001	0.0005	0.0325
$36\times 6\times 6$	0.0292	0.1510	0.0006	0.0028	0.1836
$64\times 8\times 8$	0.0677	0.5428	0.0031	0.0137	0.6273
$100\times 10\times 10$	0.1424	1.4185	0.0133	0.0570	1.6312
$144\times 12\times 12$	0.2577	3.1907	0.0478	0.1748	3.6710

Example 3 Consider a model for the quantum harmonic oscillator:

(-\frac{1}{2}\Delta+\frac{1}{2}r^{2})u=\lambda u\ \ \ \ in\ \mathbb{R}^{3}.

The minimum eigenvalue $\lambda=1.5$ and the associated eigenfunction $u=\pi^{-\frac{3}{4}}\text{e}^{-r^{2}/2}$ . Since the eigenfunctions decay exponentially, in our calculation, the following eigenvalue problem is considered.

\displaystyle\left\{\begin{aligned} (-\frac{1}{2}\Delta+\frac{1}{2}r^{2})u&=\lambda u\ \ \ \ in\ \Omega,\\ u&=0\ \ \ \ \ \ on\ \partial\Omega,\end{aligned}\right.

where $\Omega=[-5.0,5.0]^{3}$ and $r=(x^{2}+y^{2}+z^{2})^{\frac{1}{2}}$ . Let $L=10.0$ be the side length of $\Omega$ . The numerical results of the minimum eigenvalue and associated eigenfunction are presented in Figures 11-12 and Tables 5-6. Figure 11 supports the convergence results of Theorem 4.2. For this linear eigenvalue problem with symmetric solution, our results illustrate the efficiency of the symmetrized two-scale finite element method (Algorithm 4.3).

Table 5: Example 3: The running time (s) by the two-scale FEM.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$u_{H\textbf{e}}$	$u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}$	$u_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}}$	$u_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}$	$B_{H\textbf{e}}^{h}u_{h\textbf{e}}$	$B_{H\textbf{e}}^{h}\lambda_{h\textbf{e}}$	Total
$40\times 20\times 20$	0.7060	1.7362	1.6241	1.8173	0.1405	0.0001	6.0242
$62\times 25\times 25$	1.6529	5.3517	5.1389	6.1617	0.5373	0.0001	18.8426
$90\times 30\times 30$	3.5823	13.2750	13.4933	18.8052	1.6285	0.0001	50.7844
$122\times 35\times 35$	7.4273	35.3173	41.0032	40.0038	4.0751	0.0001	127.8268
$160\times 40\times 40$	13.1982	113.4475	105.7242	146.7036	9.8352	0.0016	388.9103

Table 6: Example 3: The running time (s) by Algorithm 4.3.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$u_{H\textbf{e}}$	$u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}$	$u_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}},u_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}$	$B_{H\textbf{e}}^{h}u_{h\textbf{e}}$	$B_{H\textbf{e}}^{h}\lambda_{h\textbf{e}}$	Total
$40\times 20\times 20$	0.7060	1.7362	0.0009	0.1338	0.0001	2.5770
$62\times 25\times 25$	1.6529	5.3517	0.0029	0.5145	0.0001	7.5221
$90\times 30\times 30$	3.5823	13.2750	0.0118	1.6162	0.0001	18.4854
$122\times 35\times 35$	7.4273	35.3173	0.0341	3.8900	0.0001	46.6688
$160\times 40\times 40$	13.1982	113.4475	0.0953	9.3284	0.0003	136.0697

Example 4 Consider the Kohn-Sham equation for the helium atom:

\left(-\frac{1}{2}\Delta-\frac{2}{|r|}+\int\frac{\rho(r^{\prime})}{|{r}-{r}^{\prime}|}d{r}^{\prime}+V_{xc}\right)u=\lambda u\ \ \ in\ \mathbb{R}^{3},

with $\|u\|_{0,\mathbb{R}^{3}}=1$ , where $|r|=(x^{2}+y^{2}+z^{2})^{\frac{1}{2}}$ and $\rho=2|u|^{2}$ . In our computation, we solve the following nonlinear eigenvalue problem [5]: find $(\lambda,u)\in\mathbb{R}\times H^{1}_{0}(\Omega)$ such that $\|u\|_{0,\Omega}=1$ and

\displaystyle\left\{\begin{aligned} \left(-\frac{1}{2}\Delta-\frac{2}{|r|}+\int_{\Omega}\frac{\rho(r^{\prime})}{|{r}-{r}^{\prime}|}d{r}^{\prime}+V_{xc}\right)u&=\lambda u\ \ \ in\ \Omega,\\ u&=0\ \ \ \ \ on\ \partial\Omega,\end{aligned}\right.

(5.1)

where $\Omega=[-5.0,5.0]^{3}$ . We choose $V_{xc}(\rho)=-\frac{3}{2}\alpha(\frac{3}{\pi}\rho)^{1/3}$ with $\alpha=0.77298$ . Let $L=10.0$ be the side length of $\Omega$ .

For this nonlinear eigenvalue problem, it has been shown in [9] that the two-scale finite element method is more economic than the standard finite element method. Because the first eigenfunction of (5.1) is symmetric, the two-scale finite element discretization proposed in [9] can be also combined with Algorithm 3.1. To illustrate this observation, numerical results are presented in Figure 13 and Table 7. One can see that by combining the idea of Algorithm 3.1 with the two-scale finite element method, the time consumption is reduced further.

Table 7: Example 4: The running time (s) by combining the two-scale FEM [9] with Algorithm 3.1.

$\frac{L}{h}\times\frac{L}{H}\times\frac{L}{H}$	$u_{H\textbf{e}}$	$u_{h\textbf{e}_{1}+H\hat{\textbf{e}}_{1}}$	$u_{h\textbf{e}_{2}+H\hat{\textbf{e}}_{2}},u_{h\textbf{e}_{3}+H\hat{\textbf{e}}_{3}}$	$u_{H\textbf{e}}^{h}$	$\lambda^{h}_{H\textbf{e}}$	Total
$32\times 18\times 18$	9.62	18.71	0.0006	0.16	7.55	36.10
$40\times 20\times 20$	15.11	34.77	0.0012	0.29	14.55	64.84
$48\times 22\times 22$	19.87	46.39	0.0015	0.54	26.05	93.00
$58\times 24\times 24$	27.12	79.33	0.0028	0.91	46.61	154.25
$68\times 26\times 26$	36.50	112.12	0.0042	1.47	76.90	227.41
$78\times 28\times 28$	49.10	187.78	0.0062	2.22	116.20	355.92

6 Conclusions

In this paper, a symmetrization algorithm is developed to perform a vector transformation associated with symmetric functions. By combining the symmetrization algorithm with the two-scale finite element methods, some symmetrized two-scale finite element methods are proposed for the elliptic source and eigenvalue problems with symmetric solutions. Numerical analyses and numerical results show that the symmetrized two-scale finite element methods save computational cost significantly compared with the corresponding two-scale finite element methods. The symmetrization algorithm can also be combined with the postprocessed two-scale finite element method to obtain efficient numerical methods [11, 15]. Besides, some two-scale finite element discretizations have been developed for a class of integral equation and integrodifferential equation [3, 14, 17], which can be also combined with the symmetrization algorithm to get more efficient algorithms. Moreover, we believe that our symmetrization algorithm is significant for solving higher dimensional problems on tensor product domains.

Acknowledgments

The authors thank Professor Huajie Chen for enlightening discussions. P. Hou was partially supported by the National Natural Science Foundation of China (grant 11971066). F. Liu was partially supported by the National Natural Science Foundation of China (grant 11771467) and the disciplinary funding of Central University of Finance and Economics. A. Zhou was partially supported by the National Key R & D Program of China under grants 2019YFA0709600 and 2019YFA0709601.

References

[1] I. Babuska and J.E. Osborn. Eigenvalue problems. In Handbook of Numerical Analysis, volume II, pages 641–787. North-Holland, Amsterdam, 1991.
[2] M. Bachmayr, G. Dusson, and C. Ortner. Polynomial approximation of symmetric functions. arXiv, 2109.14771, 2021.
[3] H. Chen, F. Liu, N. Reich, C. Winter, and A. Zhou. Two-scale finite element discretizations for integrodifferential equations. J. Integ. Equ. Appl., 23:351–381, 2011.
[4] P.G. Ciarlet. The Finite Element Method for Elliptic Problems, volume 4 of Studies in Mathematics and its Applications. North-Holland Publishing Co., Amsterdam, 1978.
[5] X. Dai and A. Zhou. Three-scale finite element discretizations for quantum eigenvalue problems. SIAM J. Numer. Anal., 46:295–324, 2008.
[6] J. Fang, X. Gao, and A. Zhou. A symmetry-based decomposition approach to eigenvalue problems. J. Sci. Comput., 57:638–669, 2013.
[7] X. Gao, F. Liu, and A. Zhou. Three-scale finite element eigenvalue discretizations. BIT, 48(3):533–562, 2008.
[8] J. Han, Y. Li, L. Lin, J. Lu, J. Zhang, and L. Zhang. Polynomial approximation of symmetric functions. arXiv, 1912.01765, 2022.
[9] P. Hou and F. Liu. Two-scale finite element discretizations for nonlinear eigenvalue problems in quantum physics. Adv. Comput. Math., 47:59, 2021.
[10] D. Joyner. Adventures in group theory: Rubik’s cube, Merlin’s machine, and other mathematical toys. Johns Hopkins University Press, second edition, 2008.
[11] F. Liu, M. Stynes, and A. Zhou. Postprocessed two-scale finite element discretizations, part I. SIAM J. Numer. Anal., 49:1947–1971, 2011.
[12] F. Liu and A. Zhou. Two-scale finite element discretizations for partial differential equations. J. Comput. Math., 24:373–392, 2006.
[13] F. Liu and A. Zhou. Localizations and parallelizations for two-scale finite element discretizations. Commun. Pure Appl. Anal., 6(3):757–773, 2007.
[14] F. Liu and A. Zhou. Two-scale Boolean Galerkin discretizations for Fredholm integral equations of the second kind. SIAM J. Numer. Anal., 45:296–312, 2007.
[15] F. Liu and J. Zhu. Two-scale sparse finite element approximations. Sci. China Math., 59(4):789–808, 2016.
[16] C. Pflaum and A. Zhou. Error analysis of the combination technique. Numer. Math., 84:327–350, 1999.
[17] Y. Xu and A. Zhou. Fast Boolean approximation methods for solving integral equations in high dimensions. J. Integ. Equ. Appl., 16:83–110, 2004.