A two-grid finite element method for the Schrödinger–Poisson system

Abstract

1. Introduction

The Schrödinger–Poisson (SP) system describes the thermodynamical and electrostatic equilibrium of electrons confined in narrow quantum wells and incorporates all relevant quantum phenomena. It is widely used in the simulation of quantum devices Ando et al. (1982); Vinter (1984); Beaumont & Sotomayor-Torres (1990). Much attention has been paid to the study of one- and two-dimensional problems Stern (1972); Ohkura (1990); Wu et al. (1996); Spinelli et al. (1998). To avoid numerical challenges of solving coupled quantum systems, physicists and engineers came up with approximate models instead of solving the microscopic Schrödinger model in Ancona & Tiersten (1987); Ancona & Iafrate (1989); Ben Abdallah & Unterreiter (1998); Pinnau (2002); de Falco et al. (2005); Mu & Zheng (2023). However, for extremely small devices like nanowires, numerical simulations based on quantum models become more and more preferable to fully account for quantum effects in such systems in Laux & Frank (1986); Landauer (1988).

The coupled SP system bears to be a popular model for semiconductor devices in the quantum regime. Based on the assumption that the bounded domain |$\varOmega \in \mathbb{R}^{d}$|  |$(d\leq 3)$| is of |${\cal C}^{2}$|-class or a convex polyhedron, the existence and uniqueness of the solution to the SP model are established in Nier (1993); Kaiser & Rehberg (1997) by using traces of differential operators. In the numerical aspect, Ben Abdallah, Cáceres, Carrillo and Vecil proposed a Newton–Raphson method for solving the SP problem based on the facts that their Schrödinger equation is one-dimensional (1D) and all eigenvalues are single in Ben Abdallah et al. (2009). Their numerical experiments show that the Newton–Raphson iteration is much superior to the decoupled Gummal iterative method. In Zhang et al. (2014), Zhang, Cao and Luo proposed a multi-scale asymptotic method for the stationary SP model with rapidly oscillating discontinuous coefficients. The finite element approximation to the homogenized model is studied. They proved the convergence of the numerical solution as the mesh size tends to zero. However, the convergence rate of numerical solutions is not obtained.

The objectives of this work are twofold. The first one is to propose a two-grid finite element method for solving the SP model on a bounded Lipschitz polyhedral domain |$\varOmega \subset{\mathbb{R}}^{3}$|⁠. Unlike Nier (1993); Kaiser & Rehberg (1997) where the solution is assumed to be |$H^{2}$|-regular, the domain here is not necessarily convex and the solution of the Poisson equation has a lower regularity, that is, |$V\in H^{1+s}(\varOmega )$| for some |$1/2<s\le 1$|⁠. The two-grid method is to solve the nonlinear Poisson equation on a fine mesh |${\cal T}_{h}$| of |$\varOmega $| and solve the Schrödinger eigenvalue problem on a coarse mesh |${\cal T}_{H}$|⁠. Provided that |$H\gg h$|⁠, the computational complexity for solving the eigenvalue problem is greatly reduced. Another issue concerns the electron-density operator |$n[V]$| in the Poisson equation. It is defined by an infinite series whose terms depend on eigenvalues and eigenfunctions of the Schrödinger equation. In practice, we must truncate the series into the sum of finitely many terms. To further reduce the computational complexity, we define the discrete electron-density operator |$n_{H}$| by the sum of |$L_{H}$| terms, where |$H$| is the mesh size of |${\cal T}_{H}$| and |$L_{H}=O(|\!\ln H|^{3/2})$|⁠. Therefore, only |$L_{H}$| eigenvalues and eigenfunctions of the discrete Schrödinger equation are computed on the coarse mesh.

A priori error estimate is obtained for linear Lagrangian finite element approximations to both the Poisson equation on |${\cal T}_{h}$| and the Schrödinger equation on |${\cal T}_{H}$|⁠. Provided that the mesh sizes satisfy |$h=O(H^{2})$|⁠, the |$H^{1}$|-error between the discrete electric potential |$V_{h}$| and the exact electric potential |$V$| is bounded by

where the constant |$C$| depends only on |$\left \|{V}\right \|_{L^{2}(\varOmega )}$| and given information of the system. If |$V\in H^{2}(\varOmega )$|⁠, optimal error estimate is obtained in the |$H^{1}$|-norm. Here for the SP system, an interesting observation is that |$\left \|{V}\right \|_{H^{1+s}(\varOmega )}$| is bounded by a constant which depends only on |$\left \|{V}\right \|_{L^{2}(\varOmega )}$|⁠.

The second objective is to propose an efficient method for solving the discrete nonlinear SP system. In Ben Abdallah et al. (2009), the Newton–Raphson method for the SP system is derived upon the feature that all eigenvalues of the 1D Schrödinger equation are single. The argument does not apply to our case since the Schrödinger equation is three-dimensional and multiple eigenvalues are unavoidable. Fortunately, we can follow Kaiser and Rehberg in Kaiser & Rehberg (1997) to derive an explicit expression of the Fréchet derivative |$n_{H}^{\prime}$| of the discrete electron-density operator by using a trace formula. This enables us to design the Newton–Raphson method for the discrete SP problem. Moreover, we prove that the convergence of the Newton–Raphson method is of second-order.

The Newton–Raphson method requires to solve a linearized problem of the nonlinear Poisson equation on |${\cal T}_{h}$|⁠. Although computing |$n_{H}$| only needs |$L_{H}$| eigenvalues and eigenfunctions of the Schrödinger equation, its Fréchet derivative |$n^{\prime}_{H}$| depends on all |$N_{H}$| eigenvalues and eigenfunctions. Building the stiffness matrix corresponding to |$n^{\prime}_{H}$| needs |$O\big (N_{h}N_{H}|\!\ln H|^{3/2}\big )=O\big (N_{h}^{3/2}|\!\ln H|^{3/2}\big )$| computations, provided that |$h=O(H^{2})$|⁠, where |$N_{h}$| and |$N_{H}$| are the numbers of DOFs (degrees of freedom) on |${\cal T}_{h}$| and |${\cal T}_{H}$|⁠, respectively. To reduce the computational complexity further, we propose an inexact Newton method by replacing |$n^{\prime}_{H}$| with an approximate derivative |$\tilde n_{H}^{\prime}$| which depends only on |$L_{H}$| eigenvalues and eigenfunctions. Treating the term concerning |$\tilde n_{H}^{\prime}$| only needs |$O(N_{h}|\!\ln H|^{3})$| operations. Numerical experiments show that the inexact Newton method is very efficient and superior to fixed-point iterations.

1.1 Outline

The paper is organized as follows. In Section 2, we introduce the SP model and present some preliminary results. In Section 3, we present a two-grid finite element discrete for the SP system and prove the well-posedness of the discrete problem. In Section 4, we establish a priori error estimate of numerical solution in |$H^{1}$|-norm regarding the regularity of solution. In Section 5, we propose two Newton-type methods for solving the nonlinear discrete SP problem. Convergence rate of the Newton–Raphson method is obtained. In Section 6, we present two numerical experiments to verify the convergence rate of the two-grid finite element method and demonstrate the efficiency of the inexact Newton method.

1.2 Notation

We introduce the notation to be used in this paper. Let |$\varOmega \subset{\mathbb{R}}^{3}$| be a bounded Lipschitz polyhedral domain and |$\partial \varOmega $| its boundary.

• The space of square integrable functions on |$\varOmega $| is denoted by and its inner product is denoted by |$(u,v)$|⁠.

• |$H^{m}(\varOmega )$| denotes the Sobolev space whose functions have square-integrable partial derivatives up to order |$m\ge 0$|⁠, and |$H^{1}_{0}(\varOmega )= \{v\in H^{1}(\varOmega ): v|_{\partial \varOmega }=0\}$|⁠.

• |$f$| denotes the Boltzmann or Fermi–Dirac distribution upon a given Fermi level |$\varepsilon _{F}$| (cf. Nier (1993))

  $$ \begin{align}& f(t)= f_{0}\,e^{-\mu (t-\varepsilon_{F})}, \qquad f(t)= f_{0}/(1+e^{\mu (t-\varepsilon_{F})}), \qquad f_{0}>0.\end{align} $$

  (1.1)

• |$V_{0}\in L^{2}(\varOmega )$| denotes an external electric potential.

• |$u\in L^{2}(\varOmega )$| denotes an arbitrary, but fixed function and |$\varTheta _{u}:=\left \|{u+V_{0}}\right \|_{L^{2}(\varOmega )}$|⁠.

• |$s\in (1/2,1]$| is the regularity constant such that, for any ⁠, the solution of the elliptic boundary value problem satisfies |$z\in H^{1+s}(\varOmega )$| (cf. Amrouche et al. (1998) and (Monk, 2003, Theorem 3.18)), where

  $$ \begin{align}& -\varDelta z = \xi \quad \hbox{in}\;\;\varOmega,\qquad z=0 \quad \hbox{on}\;\;\partial\varOmega.\end{align} $$

  (1.2)

• Suppose |$H$| is a separable Hilbert space and |$\{\varphi _{n}\}^{\infty }_{n=1}$| is an orthonormal basis of |$H$|⁠. For any |${\cal A}\in{\cal L}(H)$|⁠, |$\mathrm{tr}\,{\cal A}$| denotes the trace of |${\cal A}$| and is defined by |$\mathrm{tr}\,{\cal A} =\sum \limits ^\infty _{n=1}(\varphi _{n},{\cal A}\varphi _{n})_{H}$|⁠.

• An operator |${\cal A}\in{\cal L}(H)$| is of trace-class if and only if |$\mathrm{tr}\big ({\cal A}^{*}{\cal A}\big )^{1/2} <\infty $|⁠, where |${\cal A}^{*}$| denotes the adjoint of |${\cal A}$| with respect to |$(\cdot ,\cdot )_{H}$|⁠. The family of all trace-class operators is denoted by |$\mathscr{B}_{1}$|⁠.

• Some properties of the trace are listed below (cf. (Reed & Simon, 1975, Section VI.6)).

  • If |${\cal A}\in \mathscr{B}_{1}$|⁠, then |$\mathrm{tr}\,{\cal A}$| is independent of the choose of orthonormal basis.

  • |$\mathrm{tr}$| is a linear mapping from |$\mathscr{B}_{1}$| to |${\mathbb{R}}$|⁠.

  • If |${\cal A}\in \mathscr{B}_{1}$| and |${\cal B}\in{\cal L}(H)$|⁠, then |${\cal A}{\cal B},{\cal B}{\cal A}\in \mathscr{B}_{1}$| and |$\mathrm{tr}({\cal A}{\cal B})=\mathrm{tr}({\cal B}{\cal A})$|⁠.

2. The SP system

The SP model is the coupling of the Poisson equation for the electric potential and the eigenvalue problem of the Hamitonian. After proper nondimensionalization, the Poisson equation has the form

where |$V$| is the electrostatic potential, |$n_{D}\in L^{2}(\varOmega )$| is the doping profile and is the electron-density operator defined by

Here |$(\varepsilon _{l},\psi _{l})$|⁠, |$l\ge 1$|⁠, denote the eigenpairs of the Hamiltonian |${\cal H}[u] = -\varDelta +u+V_{0}$|⁠, namely

(2.3)

We write |$\varepsilon _{l}[u]$| to specify the dependency of |$\varepsilon _{l}$| on |$u$| and order them monotonically |$\varepsilon _{1}[u]<\varepsilon _{2}[u]\le \cdots $| by counting multiplicities.

From (Kaiser & Rehberg, 2000, Theorem 5.13), problem (2.1) has a unique solution ⁠. Moreover, the electron-density operator is Fréchet-differentiable and monotone, namely

(2.4)

Lemma 2.1.

There exists two positive constants |$C_{0},C_{1}$| depending only on |$\varOmega $| such that

(2.5)

Proof.

First we define the bilinear form by

(2.6)

The weak formulation of (2.3) is to find such that

(2.7)

Using Poincaré’s inequality and the Gagliardo–Nirenberg inequality, we have

(2.8)

This shows the right inequality of (2.5).

Similarly, we have for a constant |$\hat{C}>0$|⁠. Poincaré’s inequality implies ⁠. This yields the left inequality of (2.5).

Remark 2.2.

Let |$0<\lambda _{1}\le \lambda _{2}\le \cdots \le \lambda _{l}\le \cdots $| be the eigenvalues of ⁠. From (Kaiser & Rehberg, 1997, equation (2.1)) we know that |$\lim \limits _{l\to \infty }(l^{-2/3}\lambda _{l}) = C_\varOmega>0$|⁠. By arguments similar to (2.8), we can find two positive constants |$C_{\min }$| and |$C_{\max }$| depending only on |$\varOmega $| such that

$$ \begin{align} C_{\min}(l^{2/3}-\varTheta_{u}^{4})\le \varepsilon_{l}[u]\le C_{\max}(l^{2/3}-\varTheta_{u}^{4})\quad \hbox{as}\;\; l\to\infty.\end{align} $$

(2.9)

Theorem 2.3.

Suppose is the solution of problem (2.1). There exist two positive constants |$C$| and |$\alpha $| depending only on |$\varOmega $| such that

Proof.

Let |$\varepsilon _{l},\psi _{l}$| denote the eigenvalues and eigenfunctions of the Hamiltonian |${\cal H}[V]$|⁠. From Lemma 2.1, it is easy to see that

(2.10)

Using (1.2) and Remark 2.2, we get

for a constant |$\alpha>0$| depending only on |$\varOmega $|⁠. The proof is finished.

3. A two-grid finite element method

Let |${\cal T}_{H}$| be a coarse mesh of |$\varOmega $| and |${\cal T}_{h}$| a fine mesh obtained by refining |${\cal T}_{H}$|⁠. Both of them are quasi-uniform and consist of shape-regular tetrahedra. The mesh sizes of |${\cal T}_{H}$| and |${\cal T}_{h}$| are denoted by |$H$| and |$h$|⁠, respectively. The lowest-order continuous finite element spaces on |${\cal T}_{H}$| and |${\cal T}_{h}$| are defined by

where |$P_{1}(K)$| is the space of linear polynomials on |$K$|⁠. Clearly |$X_{H}\subset X_{h}$|⁠. Write |$N_{H}=\operatorname{dim} (X_{H})$| and |$N_{h}=\operatorname{dim} (X_{h})$| for convenience.

3.1 The discrete problem

To propose the discrete problem, we introduce a truncated distribution function |$f_{H}\in C^\infty ({\mathbb{R}})$| which satisfies |$f_{H}^{\prime}\le 0$| and

(3.1)

Since |$f$| deceases exponentially, we can assume that, with a constant |$C$| independent of |$H$|⁠,

The two-grid finite element method for solving (2.1) is to seek |$V_{h}\in X_{h}$| such that

where is the discrete electron-density operator defined as follows:

(3.4)

Here |$\varepsilon _{l,H}[u]\in{\mathbb{R}}$| and |$\psi _{l,H}\in X_{H}$| solve the discrete eigenvalue problem on the coarse mesh

(3.5)

The eigenvalues are ordered monotonically |$\varepsilon _{1,H}[u]\le \varepsilon _{2,H}[u] \le \cdots \le \varepsilon _{N_{H},H}[u]$| by counting multiplicities. The Galerkin approximation of the eigenvalue problem implies |$\varepsilon _{l}[u] \le \varepsilon _{l,H}[u]$|⁠.

3.2 Frechét derivative of |$n_{H}[u]$|

Following Kaiser and Rehberg Kaiser & Rehberg (1997), we define the iterated primitive functions of |$f_{H}$| by

Since |$f_{H}(t)=0$| for |$t\ge \varepsilon _{F}+2\mu ^{-1}|\!\ln H|+1$|⁠, we also have |$f_{H}^{(-j)}(t)=0$| for |$t\ge \varepsilon _{F}+2\mu ^{-1}|\!\ln H|+1$| for |$j=1,2$|⁠. For simplicity, we make the following assumptions which are rather mild for many practical distribution functions.

• (A1) |$f_{H}$| is the trace of some holomorphic function |$f_{H}(z)$| in a bounded domain |${\mathbb{G}}\subset{\mathbb{C}}$| which covers the interval |$[\varepsilon _{1,h},\varepsilon _{N_{h},h}]$| and is symmetric about the real axis of |${\mathbb{C}}$|⁠.

• (A2) There is a bounded contour |$\varGamma \subset{\mathbb{G}}$| and a constant |$\delta>0$| independent of |$H$| such that |$\varGamma $| encloses |$[\varepsilon _{1,h},\varepsilon _{N_{h},h}]$| and satisfies

  $$ \begin{align}& \mathrm{dist}(\varepsilon_{l,H},\varGamma)\geq\delta,\qquad 1\le l\le N_{H}.\end{align} $$

  (3.9)

  Moreover, |$\varGamma $| is oriented such that the interval |$[\varepsilon _{1,H}, \varepsilon _{N_{H},H}]$| lies on its left-hand side.
• (A3) |$f_{H}(z)$| has primitive functions |$f_{H}^{(-1)}(z)$| and |$f_{H}^{(-2)}(z)$| which are holomorphic in |${\mathbb{G}}$| and have real values on |${\mathbb{R}}$|⁠. There exists a constant |$C>0$| independent of |$H$| such that

  $$ \begin{align}& \int_{\varGamma}\left|{z}\right|{}^{k}\big|f_{H}^{(-j)}(z)\big|\,\mathrm{d} z \le C\quad \hbox{for}\;\; j=1,2\;\; \hbox{and}\;\; k=0, 1,2,3.\end{align} $$

  (3.10)

Lemma 3.4.

The discrete density operator |$n_{H}$| is Fréchet-differentiable at any |$u\in L^{2}(\varOmega )$| and the Fréchet derivative ⁠. Moreover, for any ⁠, we have

$$ \begin{align} n^{\prime}_{H}[u](\xi) =\sum^{N_{H}}_{i,j=1}\frac{f_{H}(\varepsilon_{i,H})-f_{H}(\varepsilon_{j,H})}{\varepsilon_{i,H}-\varepsilon_{j,H}} \Big(\int_{\varOmega}\xi\psi_{i,H}\psi_{j,H}\Big)\psi_{i,H}\psi_{j,H}.\end{align} $$

(3.11)

In case of |$\varepsilon _{i,H}=\varepsilon _{j,H}$|⁠, the coefficient |$\big [f_{H}(\varepsilon _{i,H})-f_{H}(\varepsilon _{j,H})\big ]/(\varepsilon _{i,H}-\varepsilon _{j,H})$| is replaced with |$f_{H}^{\prime}(\varepsilon _{i,H})$|⁠.

Proof.

First we define the discrete Schrödinger operator |${\cal H}_{H}[u]:X_{H}\to X_{H}$| as follows: for any |$\xi _{H}\in X_{H}$|⁠, |${\cal H}_{H}[u](\xi _{H})\in X_{H}$| is the unique solution to the discrete problem

$$ \begin{align} \big({\cal H}_{H}[u](\xi_{H}),\varphi_{H}\big) = a(u;\xi_{H},\varphi_{H}) \quad \forall\,\varphi_{H}\in X_{H}.\end{align} $$

(3.12)

Clearly the eigenpairs |$(\varepsilon _{l,H}, \psi _{l,H})$|⁠, |$1\le l\le N_{H}$|⁠, are also eigenpairs of |${\cal H}_{H}[u]$|⁠.

For any ⁠, let be the multiplying operator defined by

$$ \begin{align} ({\cal M}_{H}[\eta](v),\phi_{H}) =(\eta v,\phi_{H} ) \quad \forall\,\phi_{H}\in X_{H}.\end{align} $$

(3.13)

Since |$X_{h}\subset C(\overline \varOmega )$|⁠, the right-hand side is well-defined for any |$\phi _{H}\in X_{H}$|⁠. In fact, |${\cal M}_{H}[\eta ](v)$| is the -projection of |$\eta v$| onto |$X_{H}$|⁠.

For any |$z\in \varGamma $|⁠, we denote the resolvent of |${\cal H}_{H}[u]$| by |${\cal R}_{z}[u]=(z-{\cal H}_{H}[u])^{-1}$|⁠. For |$k=0,1,2$|⁠, we define the linear operators |${\cal U}_{H,k}:X_{H}\to X_{H}$| by

$$ \begin{align*} {\cal U}_{H,k}:=\frac{k!}{2\pi{\mathbf{i}}} {\cal M}_{H}[\eta]\int_{\varGamma}f^{(-k)}_{H}(z){\cal R}^{k+1}_{z}[u]\,\mathrm{d} z. \end{align*} $$

By the Residue theorem and the definition of trace in section 1.2, it is easy to see that

$$ \begin{align} \mathrm{tr}({\cal U}_{H,k})=\,& \frac{k!}{2\pi{\mathbf{i}}}\sum^{N_{H}}_{l=1} \int_{\varOmega} {\cal M}_{H}[\eta](\psi_{l,H}) \bigg(\int_{\varGamma}f_{H}^{(-k)}(z){\cal R}_{z}^{k+1}[u](\psi_{l,H})\,\mathrm{d} z\bigg) \notag \\ =\,&\frac{k!}{2\pi{\mathbf{i}}}\sum^{N_{H}}_{l=1} \int_{\varOmega} \eta\psi_{l,H} \bigg(\int_{\varGamma}\frac{f_{H}^{(-k)}(z)}{(z-\varepsilon_{l,H})^{k+1}}\psi_{l,H}\,\mathrm{d} z\bigg) \notag \\ =\,&\sum^{N_{H}}_{l=1}f(\varepsilon_{l,H})\int_\varOmega \eta\psi^{2}_{l,H} \notag \\ =\,& (n_{H}[u],\eta).\end{align} $$

(3.14)

First we write |${\cal N}_{H}[u](\xi ):= \sum \limits ^{N_{H}}_{i,j=1}\frac{f_{H}(\varepsilon _{i,H})-f_{H}(\varepsilon _{j,H})}{\varepsilon _{i,H}-\varepsilon _{j,H}} \big (\xi ,\psi _{i,H}\psi _{j,H}\big ) \psi _{i,H}\psi _{j,H}$|⁠. Since ⁠, we have ⁠, and by arguments similar to (3.14), find that

$$ \begin{align*} \big({\cal N}_{H}[u](\xi),\eta\big) =\sum_{j=1}^{2} \mathrm{tr}\bigg(\frac{{\cal M}_{H}[\eta]}{2\pi{\mathbf{i}}}\int_\varGamma f^{(-1)}_{H}(z) {\cal R}^{j}_{z}[u]{\cal M}_{H}[\xi]{\cal R}^{3-j}_{z}[u]\,\mathrm{d} z\bigg). \end{align*} $$

It is left to show |$n^{\prime}_{H}[u] ={\cal N}_{H}[u]$|⁠.

In view of (3.14), it is easy to see that

$$ \begin{align} &(n_{H}[u+\xi]-n_{H}[u]-{\cal N}_{H}[u](\xi),\eta) \\ &\quad =\, \mathrm{tr}\bigg(\frac{{\cal M}_{H}[\eta]}{2\pi{\mathbf{i}}}\int_\varGamma f^{(-1)}_{H}(z) \Big({\cal R}^{2}_{z}[u+\xi]-{\cal R}^{2}_{z}[u] -{\cal R}^{2}_{z}[u]{\cal M}_{H}[\xi]{\cal R}_{z}[u]-{\cal R}_{z}[u]{\cal M}_{H}[\xi]{\cal R}^{2}_{z}[u]\Big)\,\mathrm{d} z\bigg). \notag\end{align} $$

(3.15)

Since |$\mathrm{dist}(\varepsilon _{l,H}[u+\xi ],\varGamma )\geq \delta $| by Assumption (A2), we choose small enough such that |$\mathrm{dist}(\varepsilon _{l,H}[u+\xi ],\varGamma )\geq \delta /2$| for |$1\le l\le N_{H}$|⁠. It follows from (3.9) that

$$ \begin{align} \left\|{{\cal R}_{z}[u+\xi]}\right\|_{{\cal L}(X_{H})}\leq 2\delta^{-1},\qquad \left\|{{\cal R}_{z}[u]}\right\|_{{\cal L}(X_{H})}\leq \delta^{-1}.\end{align} $$

(3.16)

Since norms are equivalent on |$X_{H}$|⁠, there is a generic constant |$C_{H}$|⁠, depending on |$H$| such that

It is easy to see that |${\cal D}_{z}[\xi ]:={\cal R}_{z}[u+\xi ]-{\cal R}_{z}[u]$| satisfies

By the identity |${\cal R}_{z}[u+\xi ]-{\cal R}_{z}[u]={\cal R}_{z}[u+\xi ]{\cal M}_{H}[\xi ]{\cal R}_{z}[u]$| and proper arrangements, we get

Inserting the above estimate into the right-hand side of (3.15), we end up with

This shows

Then |$n_{H}$| is Fréchet-differentiable at and |$n^{\prime}_{H}[u]={\cal N}[u]$|⁠.

3.3 Well-posedness of the discrete problem

We end this section by proving the well-posedness of problem (3.3).

Theorem 3.8.

The discrete problem (3.3) has a unique solution. Moreover, there is a constant |$C>0$| depending only on |$\varOmega $| such that, with |$\varTheta _{0}=\left \|{V_{0}}\right \|_{L^{2}(\varOmega )}$|⁠,

(3.18)

Proof.

Corollary 3.5 implies

(3.19)

Classical theories for monotone operators show that problem (3.3) has a unique solution (cf. (Zeider, 1990, Section 25.4, pp. 503–506)).

Next we consider the discrete eigenvalue problem: find |$(\lambda _{l,H},\phi _{l,H})\in{\mathbb{R}}\times X_{H}$| which satisfy

Taking |$v_{H}=\phi _{l,H}$| and using Poincaré’s inequality, we can find a constant |$C$| depending only on |$\varOmega $| such that the estimates hold

There exists a generic constant |$C$| depending only on |$\varOmega $| such that

Since |$\lambda _{l,H}\ge \lambda _{l}$| and |$f$| decays exponentially, we deduce from (2.9) that

Finally, taking |$v_{h}=V_{h}$| in (3.5), we find that

Using (3.19) and Poincaré’s inequality, we obtain

The proof is finished.

4. Finite element error estimates

The purpose of this section is to establish the error estimate between the exact potential |$V$| and the discrete potential |$V_{h}$|⁠. First we estimate the error between |$n[u]$| and |$n_{H}[u]$|⁠. From Lemma 3.2, we have |$L_{H}=L_{H}[u]\le L=L[u]\ll N_{H}$|⁠.

4.1 Solution operators

By the Gagliardo–Nirenberg inequality, Poincaré’s inequality and Young’s inequality, there exists a constant |$C_\varOmega $| depending only on |$\varOmega $| such that, for ⁠,

Define |$\rho _{u}:=C_\varOmega \varTheta _{u}^{4}$|⁠. Clearly |$a(\rho _{u}+u;\boldsymbol{\cdot },\boldsymbol{\cdot })$| is continuous and coercive on ⁠, namely

(4.1)

where |$C_{2},C_{3}$| are two positive constants depending only on |$\varOmega $|⁠.

For any |$g\in H^{-1}(\varOmega )$|⁠, the elliptic problem

(4.2)

has a unique solution which satisfies ⁠. This defines a solution operator ⁠, satisfying |${\cal S}[\rho _{u}+u](g)=w$|⁠. In fact, |${\cal S}[\rho _{u}+u]$| is the inverse of and satisfies, upon a constant |$C$| depending only on |$\varOmega $|⁠,

(4.3)

Similarly, we define a discrete solution operator |${\cal S}_{H}[\rho _{u}+u]: H^{-1}(\varOmega ) \to X_{H}$| by the Galerkin approximation of problem (4.2), namely

Without specifications, we write |${\cal S}={\cal S}[u+\rho _{u}]$| and |${\cal S}_{H}={\cal S}_{H}[u+\rho _{u}]$| throughout this section.

Lemma 4.1.

Suppose and let |$s\in (1/2,1]$| be the regularity constant in (1.2). Then the solution of (4.2) admits ⁠. Moreover, there exists a constant |$C>0$| depending only on |$\varOmega $| such that

(4.5)

Proof.

Set |$p:=6/(3-2s)>3$| and |$q:= p/(p-1)<3/2$|⁠. For |${\boldsymbol{z}}\in{\boldsymbol{L}}^{q}(\varOmega )$|⁠, the elliptic problem

$$ \begin{eqnarray*} \varDelta\phi_{\boldsymbol{z}} = \operatorname{div}{\boldsymbol{z}} \quad \hbox{in}\;\;\varOmega,\qquad \phi_{\boldsymbol{z}} =0 \quad \hbox{on}\;\;\partial\varOmega, \end{eqnarray*} $$

has a unique solution that satisfies

$$ \begin{align} \left\|{\phi_{\boldsymbol{z}}}\right\|_{W^{1,q}(\varOmega)}\le C\left\|{\operatorname{div}{\boldsymbol{z}}}\right\|_{[W_{0}^{1,p}(\varOmega)]^{\prime}} \le C\left\|{{\boldsymbol{z}}}\right\|_{{\boldsymbol{L}}^{q}(\varOmega)}\!.\end{align} $$

(4.6)

The solution of (4.2) satisfies

$$ \begin{eqnarray*} \left\|{\nabla w}\right\|_{{\boldsymbol{L}}^{p}(\varOmega)} = \sup_{{\boldsymbol{z}}\in{\boldsymbol{L}}^{q}(\varOmega)} \frac{(\nabla w,{\boldsymbol{z}})}{\left\|{{\boldsymbol{z}}}\right\|_{{\boldsymbol{L}}^{q}(\varOmega)}} = \sup_{{\boldsymbol{z}}\in{\boldsymbol{L}}^{q}(\varOmega)} \frac{(\nabla w,\nabla\phi_{\boldsymbol{z}})}{\left\|{{\boldsymbol{z}}}\right\|_{{\boldsymbol{L}}^{q}(\varOmega)}} = \sup_{{\boldsymbol{z}}\in{\boldsymbol{L}}^{q}(\varOmega)} \frac{(g-(\rho_{u}+u+V_{0})w,\phi_{\boldsymbol{z}})}{\left\|{{\boldsymbol{z}}}\right\|_{{\boldsymbol{L}}^{q}(\varOmega)}}. \end{eqnarray*} $$

By (4.6) and the injection |$W^{1,q}(\varOmega )\hookrightarrow L^{3}(\varOmega )$|⁠, we have

In view of (4.3) and Poincaré’s inequality, it follows that

(4.7)

This yields ⁠. The error estimate (4.5) comes directly from classical finite element error estimates.

Let denote the Banach algebra of linear and bounded operators on ⁠, endowed with the operator norm. Suppose |${\cal Q}\in \mathscr{B}$| and only has point spectrum |$\{\nu _{l}\}_{l=1}^\infty $|⁠. For any |$p\in [1,+\infty )$|⁠, we define the Schatten |$p$|-norm of |${\cal Q}$| by

An operator |${\cal Q}$| is called |$p$|-summable if |$\left \|{{\cal Q}}\right \|_{p}<\infty $|⁠. The closed ideal of |$p$|-summable operators on is denoted by |$(\mathscr{B}_{p},\left \|{\boldsymbol{\cdot }}\right \|_{p})$|⁠.

Similar to (3.13), for |$\eta \in L^{q}(\varOmega )$| with |$q\ge 2$|⁠, we define a multiplying operator |${\cal M}[\eta ]$| as follows:

It is clear that |${\cal M}[\eta ]$| maps |$L^{\frac{qr}{q-r}}(\varOmega )$| to |$L^{r}(\varOmega )$| for |$1\le r\le q$| and maps to |$H^{-1}(\varOmega )$| if |$q\ge 3$|⁠. In view of (3.13), we find that |${\cal M}_{H}[\eta ](v)$| is the |$L^{2}$|-projection of |${\cal M}[\eta ](v)$| onto the finite element space |$X_{H}$|⁠.

Lemma 4.2.

Suppose and |$\eta \in L^{6}(\varOmega )$|⁠. There exists a constant |$C$| depending only on |$\varOmega $| and such that

Proof.

In view of (4.8), the multiplying operator |${\cal M}[\eta ]$| actually maps to |$L^{3/2}(\varOmega )$| and is bounded by |$\left \|{\eta }\right \|_{L^{6}(\varOmega )}$|⁠. Then (4.3) and the injection of |$L^{3/2}(\varOmega )\hookrightarrow H^{-1}(\varOmega )$| show that

Moreover, inequality (4.7) shows that ⁠. Then we have

The results for |${\cal S}_{H}$| can be proved similarly.

4.2 Estimate of |$n[u]-n_{H}[u]$|

First we split |$n[u]-n_{H}[u]$| into two terms

The first term is easily estimated with (3.1) and the Cauchy–Schwarz inequality. In fact,

From Lemma 2.1, equation (1.1) and Remark 2.2, there exists a constant |$\alpha>0$| such that

It is left to estimate ⁠. The key ingredient is to use trace representations of |$\hat{n}[u]$| and |$n_{H}[u]$|⁠. Note that |$\mu _{l}= (\varepsilon _{l}+\rho _{u})^{-1}$| and |$\mu _{l,H}= (\varepsilon _{l,H}+\rho _{u})^{-1}$| are eigenvalues of |${\cal S}$| and |${\cal S}_{H}$|⁠, respectively, and |$\psi _{l},\psi _{l,H}$| are the eigenfunctions belonging to |$\mu _{l}$| and |$\mu _{l,H}$|⁠, respectively. Then

Clearly |$\mu _{1},\cdots ,\mu _{L}$| and |$\mu _{1,H},\cdots ,\mu _{L_{H},H}$| are included in the interval |$[0,(\varepsilon _{1}+\rho _{u})^{-1}]$|⁠. By (3.1), we can extend |$\hat{f}_{H}$| by zero to |$(-\infty ,0]$| and the extension |$\hat{f}_{H}\in C^\infty ({\mathbb{R}})$|⁠. For |$k\ge 0$|⁠, the function |$\hat{f}_{k}(t):=t^{-k}\hat{f}_{H}(t)$| is smooth with the extended value |$\hat{f}_{k}(0):=0$| and vanishes on |$(-\infty ,\gamma _{H}]$|⁠, where |$\gamma _{H}:=(2\mu ^{-1}|\!\ln H|+\varepsilon _{F}+1+\rho _{u})^{-1}$|⁠.

Similar to section 3.2, we make two assumptions.

• (A4) For any integer |$k\ge 0$|⁠, |$\hat{f}_{k}$| is the trace of some holomorphic function |$\hat{f}_{k}(z)$| in |${\mathbb{G}}$|⁠.

• (A5) There is a bounded contour |$\hat \varGamma \subset{\mathbb{G}}$| which incloses |$[0,(\varepsilon _{1}+\rho _{u})^{-1}]$| and satisfies

  $$ \begin{align}& \mathrm{dist}(t,\hat\varGamma)\geq\delta\quad\;\; \forall\, t\in [0,(\varepsilon_{1}+\rho_{u})^{-1}].\end{align} $$

  (4.11)

  Moreover, |$\hat \varGamma $| is oriented such that the interval |$\big [0,(\varepsilon _{1}+\rho _{u})^{-1}\big ]$| lies on its left-hand side.

Since |$\hat{f}_{k}\in C^\infty ({\mathbb{R}})$|⁠, assumptions (A4) and (A5) are rather mild by similar analyses given in Remark 3.3. We do not elaborate on the discussions again here.

Lemma 4.3.

For and integer |$k\ge 1$|⁠, there hold

$$ \begin{align} \big(\hat{n}[u],\eta)=\,&\mathrm{tr}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{k}(z) (z-{\cal S})^{-1}{\cal S}^{k}\, \mathrm{d} z\bigg), \end{align} $$

(4.12)

 

$$ \begin{align} \big(n_{H}[u],\eta)=\,&\mathrm{tr}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{k}(z)(z-{\cal S}_{H})^{-1}{\cal S}_{H}^{k} \, \mathrm{d} z\bigg).\end{align} $$

(4.13)

Proof.

The proof is easy by using the Residue theorem. Since |$\hat{f}_{k}\equiv 0$| in |$(-\infty ,\gamma _{H})$|⁠, we have

$$ \begin{align*} \big(\hat{n}[u],\eta) =\,& \sum^{L}_{l=1}\frac{\mu_{l}^{k}}{2\pi{\mathbf{i}}} \int_{\varOmega}\bigg\{\eta\int_{\hat\varGamma} \frac{\hat{f}_{k}(z)}{z-\mu_{l}}\, \mathrm{d} z\bigg\}\psi^{2}_{l} = \sum^{\infty}_{l=1}\frac{\mu_{l}^{k}}{2\pi{\mathbf{i}}} \int_{\varOmega}\bigg\{\eta\int_{\hat\varGamma} \frac{\hat{f}_{k}(z)}{z-\mu_{l}}\, \mathrm{d} z\bigg\}\psi^{2}_{l} \\ =\,& \mathrm{tr}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma}\hat{f}_{k}(z)(z-{\cal S})^{-1}{\cal S}^{k}\, \mathrm{d} z\bigg). \end{align*} $$

The proof of (4.13) is similar. Since |$X_{H}$| is a subspace of ⁠, we have the orthogonal decomposition with respect to the inner product of ⁠. Note that |$\{\psi _{l,H}\}_{l=1}^{N_{H}}$| provides an orthonormal basis of |$X_{H}$|⁠. We choose an orthonormal basis of |$X^\perp _{H}$| by |$\big \{\psi ^\perp _{l}\big \}^{\infty }_{l=N_{H}+1}$|⁠. In view of (4.4), the operator |${\cal S}_{H} $| maps to |$X_{H}$|⁠, and maps |$X^\perp _{H}$| to |$\{0\}$|⁠. Using the Residue theorem and the fact that |$\hat{f}_{k}\equiv 0$| in |$(-\infty ,\gamma _{H})$|⁠, we have

$$ \begin{align*} \mathrm{tr}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{k}(z)(z-{\cal S}_{H})^{-1}{\cal S}_{H}^{k}\, \mathrm{d} z\bigg) =\,& \sum^{N_{H}}_{l=1}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{k}(z)(z-{\cal S}_{H})^{-1}{\cal S}_{H}^{k}\psi_{l,H}\, \mathrm{d} z,\psi_{l,H}\bigg)\\ \,&+\sum^{\infty}_{l=N_{H}+1}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{k}(z)(z-{\cal S}_{H})^{-1}{\cal S}_{H}^{k}\psi^\perp_{l}\, \mathrm{d} z,\psi^\perp_{l}\bigg)\\ =\,&\sum^{N_{H}}_{l=1}\bigg(\frac{{\cal M}[\eta]}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{k}(z)(z-{\cal S}_{H})^{-1}{\cal S}_{H}^{k}\psi_{l,H}\, \mathrm{d} z,\psi_{l,H}\bigg)\\ =\,&\sum^{L}_{l=1}\frac{\mu_{l,H}^{k}}{2\pi{\mathbf{i}}} \int_{\varOmega}\bigg\{\eta\int_{\hat\varGamma} \frac{\hat{f}_{k}(z)}{z-\mu_{l,H}}\, \mathrm{d} z\bigg\}\psi^{2}_{l,H} \\ =\,& \big(n_{H}[u],\eta). \end{align*} $$

The proof is finished.

Proof.

For any |$\eta \in L^{6}(\varOmega )$|⁠, using (4.12) and (4.13) with |$k_{1}=k_{2}=2$|⁠, we have

$$ \begin{align} \big(\hat{n}[u]-n_{H}[u],\eta)= \mathrm{tr}\bigg(\frac{1}{2\pi{\mathbf{i}}} \int_{\hat\varGamma} \hat{f}_{4}(z) {\cal M}[\eta] \big\{{\cal S}^{2}(z-{\cal S})^{-1}{\cal S}^{2} -{\cal S}_{H}^{2}(z-{\cal S}_{H})^{-1}{\cal S}_{H}^{2}\big\}\, \mathrm{d} z\bigg).\end{align} $$

(4.14)

Direct calculations show that

$$ \begin{align} &{\cal S}^{2}({z}-{\cal S})^{-1}{\cal S}^{2}-{\cal S}_{H}^{2}({z}-{\cal S}_{H})^{-1}{\cal S}_{H}^{2} \notag \\ &\quad =\,{\cal S}^{2}({z}-{\cal S})^{-1}\big\{{\cal S}({\cal S}-{\cal S}_{H})+({\cal S}-{\cal S}_{H}){\cal S}_{H} +({\cal S}-{\cal S}_{H})({z}-{\cal S}_{H})^{-1}{\cal S}_{H}^{2}\big\} \notag \\ &\qquad + {\cal S}({\cal S}-{\cal S}_{H})({z}-{\cal S}_{H})^{-1}{\cal S}_{H}^{2} +({\cal S}-{\cal S}_{H}){\cal S}_{H} ({z}-{\cal S}_{H})^{-1}{\cal S}_{H}^{2}.\end{align} $$

(4.15)

We shall estimate the three terms on the right-hand side.

For |$p>3/2$|⁠, by (2.9) and (4.1), there exists a constant |$C>0$| depending only on |$p$| and |$\varOmega $| such that |$\mu _{l}^{p} = (\varepsilon _{l}+\rho _{u})^{-p} \le C l^{-2p/3}$|⁠. Then

$$ \begin{align} \left\|{{\cal S}_{H}}\right\|_{p} =\sum^{N_{H}}_{l=1}\mu_{l,H}^{p} \le \left\|{{\cal S}}\right\|_{p} =\sum^{\infty}_{l=1}\mu_{l}^{p} \le C.\end{align} $$

(4.16)

Note that |$\|({z}-{\cal S})^{-1}\|\le \delta ^{-1}$| for any |$z\in \hat \varGamma $|⁠. By (4.5), (4.16) and Lemma 4.2, the first term on the right-hand side of (4.15) can be estimated as follows:

where the constant |$C>0$| depends only on |$\varOmega $| and |$\delta $|⁠. The second term on the right-hand side of (4.15) can be estimated similarly and is also bounded by ⁠. For the third term, we use the formula |$\mathrm{tr}({\cal A}{\cal B})=\mathrm{tr}({\cal B}{\cal A})$| and obtain

$$ \begin{align*} \mathrm{tr}\big({\cal M}[\eta]({\cal S}-{\cal S}_{H}){\cal S}_{H}({z}-{\cal S}_{H})^{-1}{\cal S}_{H}^{2}\big) =\,& \mathrm{tr}\big({\cal S}_{H}({z}-{\cal S}_{H})^{-1}{\cal S}_{H}^{2}{\cal M}[\eta]({\cal S}-{\cal S}_{H})\big) \\ \leq\, & \left\|{{\cal S}_{H}}\right\|^{2}_{2}\left\|{({z}-{\cal S}_{H})^{-1}}\right\| \left\|{{\cal S}_{H}{\cal M}[\eta]}\right\| \left\|{{\cal S}-{\cal S}_{H}}\right\|\\ \leq\,& C\rho_{u}^{4}\left\|{\eta}\right\|_{L^{6}(\varOmega)}H^{2s}. \end{align*} $$

Substituting these estimates first into (4.15) and then into (4.14), we end up with

$$ \begin{align*} \left|{\big(\hat{n}[u]-n_{H}[u],\eta)}\right| \le C \rho_{u}^{5} H^{2s}\left\|{\eta}\right\|_{L^{6}(\varOmega)} \int_{\hat\varGamma} \big|\hat{f}_{4}(z)\big| \, \mathrm{d} z \le C \rho_{u}^{5} H^{2s}\left\|{\eta}\right\|_{L^{6}(\varOmega)}\!. \end{align*} $$

The proof is finished.

4.3 Estimate of |$V-V_{h}$|

Now we present the main result of this section, namely the a priori error estimate of |$V_{h}$|⁠.

Theorem 4.5.

Let |$\alpha $| be the constant in Theorem 2.3 and define

There exist two positive constants |$C$| and |$\beta $| independent of |$H$| and |$h$| such that

Proof.

Let |$\pi _{h}$| be the nodal interpolation operator onto |$X_{h}$|⁠. By Theorem 2.3 and standard interpolation error estimates, we obtain

$$ \begin{align*} \left\|{V-\pi_{h} V}\right\|_{H^{m}(\varOmega)}\le h^{(2-m)s}\left\|{V}\right\|_{H^{1+s}(\varOmega)} \le CD_{V} h^{(2-m)s},\qquad m=0,1. \end{align*} $$

Take |$k_{1}=1,k_{2}=2$| and |$u=V_{h}$| in (4.13). By Theorem 3.8, Lemma 4.2 and arguments similar to the proof of Lemma 4.4, we obtain

where the constant |$C$| depends only on |$\varOmega $| and |$\delta $|⁠. This shows

(4.17)

Write |$e=V-V_{h}$|⁠, |$\xi _{h}=\pi _{h}V-V_{h}$| and |$\eta =V-\pi _{h}V$| for convenience. Using (3.3) and the equality |$(\nabla V, \nabla v_{h}) = (n[V]-n_{D}, v_{h})$| for all |$v_{h}\in X_{h}$|⁠, we have

(4.18)

Since |$(n_{H}[V]-n_{H}[V_{h}], V-V_{h})\le 0$|⁠, using (2.10) and (4.17), we deduce that

$$ \begin{align} (n_{H}[V]-n_{H}[V_{h}], \xi_{h}) \le (n_{H}[V_{h}]-n_{H}[V], \eta) \le C\big(D_{0}^{4}+\varLambda_{V}\big)D_{V} h^{2s}.\end{align} $$

(4.19)

Moreover, (4.9) and Lemma 4.4 indicate that

(4.20)

Inserting (4.19)–(4.20) into (4.18) yields

The proof is finished by using Poincaré’s inequality.

5. Newton-type methods

The purpose of this section is to study Newton-type methods for solving the nonlinear discrete problem (3.3). First we derive the Newton–Raphson method and prove that it has second-order convergence rate. Then we propose an inexact Newton method by truncating the sum on the right-hand side of (3.11). The inexact Newton method has a convergence speed comparable to the Newton–Raphson method, but is vastly superior in terms of computational complexity.

5.1 Stability of |$n_{H}^{\prime}[u]$|

To study the convergence rate, we need the stability of |$n^{\prime}_{H}[u]$| with respect to |$H$| for any |$u\in L^{2}(\varOmega )$|⁠. Recall the bounded contour |$\varGamma $| and the constant |$\delta>0$| in Assumption (A2) of section 3.2. For |$z\in \varGamma $|⁠, |${\cal R}_{z}[u](z)=(z-{\cal H}_{H}[u])^{-1}$| denotes the resolvent of |${\cal H}_{H}[u]$|⁠.

Proof.

From Remark 3.6, it is easy to see that

$$ \begin{align*} \big(n^{\prime}_{H}[u](\xi),\eta\big) =\,&\frac{1}{\pi{\mathbf{i}}}\int_\varGamma f^{(-2)}_{H}(z) \mathrm{tr}\Big( {\cal M}[\eta]{\cal R}_{z}^{3}[u]{\cal M}[\xi]{\cal R}_{z}[u] +{\cal M}[\eta]{\cal R}_{z}^{2}[u]{\cal M}[\xi]{\cal R}_{z}^{2}[u] \notag\\ &+ {\cal M}[\eta]{\cal R}_{z}[u]{\cal M}[\xi]{\cal R}_{z}^{3}[u]\Big)\, \mathrm{d} z. \end{align*} $$

By Lemma 5.1 and assumption (3.10), we have

Similarly, from Lemma 3.7, we have

$$ \begin{align*} &\left|{\big(n^{\prime}_{H}[u+w](\xi),\eta\big) -\big(n^{\prime}_{H}[u](\xi),\eta\big)}\right| \\ & \quad \le \, C\int_\varGamma \big|f^{(-2)}_{H}(z)\big| \left\|{{\cal M}[w] {\cal R}_{z}[u]}\right\| \left\|{{\cal M}[\eta]{\cal R}_{z}[u+w]}\right\| \big(\left\|{{\cal M}[\xi]{\cal R}_{z}[u+w]}\right\| + \left\|{{\cal M}[\xi]{\cal R}_{z}[u]}\right\| \big)\\ &\qquad \times \big(\left\|{{\cal R}_{z}[u+w]}\right\|_{2}^{2}+\left\|{{\cal R}_{z}[u]}\right\|_{2}^{2}\big)\, \mathrm{d} z. \end{align*} $$

Again we get (5.2) by using Lemma 5.1. The proof is finished.

5.2 The Newton–Raphson method

Next we define the discrete residual operator |${\cal P}_{h}: X_{h}\to X_{h}^{\prime}$| as follows:

Then (3.3) is equivalent to the nonlinear equation

From Lemma 3.4, |${\cal P}_{h}$| is Frechét-differentiable at |$w_{h}\in X_{h}$| and the Frechét derivative satisfies

Given an approximate solution |$V^{(k)}_{h}$| with |$k\ge 0$|⁠, the Newton–Raphson method for solving (5.4) is to update |$V^{(k+1)}_{h}=V^{(k)}_{h}-r_{k}$|⁠, where the correction |$r_{k}\in X_{h}$| solves the linearized problem

By Corollary 3.5, it is clear that problem (5.8) or (5.6) has a unique solution.

Now we present the Newton–Raphson method for solving (5.4), or (3.3) equivalently.

Theorem 5.3.

Let |$\delta>0$| be the constant in (3.9). There is a constant |$C_{\mathrm{NR}}\ge 1$| depending only on |$\varOmega $| and |$\delta $| such that for any |$V^{(1)}_{h}\in X_{h}$| satisfying ⁠, the iterative sequence generated by the Newton–Raphson method satisfies

Proof.

Define |$e_{k}=V^{(k)}_{h} -V_{h}$|⁠. From (5.5), we know that

Using Corollary 3.5 and inequality (5.2), we have

By Poincaré’s inequality, there exists a constant |$C_{\mathrm{NR}}>0$| depending only on |$\delta $| and |$\varOmega $| such that

Since ⁠, we conclude the result by induction.

Now we analyze the computational complexity of the stiffness matrix associated with the term |$(n^{\prime}_{H}(r_k), v_h)$| in Algorithm 5.3. Define

Here |$b_{1},\cdots ,b_{N_{h}}$| are nodal basis functions of |$X_{h}$|⁠. By (3.11), we find that |${\mathbb{N}} ={\mathbb{K}}^\top{\mathbb{F}}{\mathbb{K}}$|⁠, where |${\mathbb{K}}$| is a |$N_{H}^{2}\times N_{h}$| matrix and |${\mathbb{F}}$| is a |$N_{H}^{2}\times N_{H}^{2}$| diagonal matrix. Their entries are given by

where |$1\le i,j\le N_{H}$| and |$1\le l\le N_{h}$|⁠. Since |$f_{H}(\varepsilon _{i,H})=0$| for all |$i>L_{H}=O(|\!\ln H|^{3/2})$|⁠, |${\mathbb{F}}$| has only |$O(N_{H}|\!\ln H|^{3/2})$| nonzero diagonal entries. Then only |$O(N_{H}|\!\ln H|^{3/2}\times N_{h})$| entries of |${\mathbb{K}}$| are used in building |${\mathbb{N}}$|⁠. For any |$\mathbf{v}\in{\mathbb{R}}^{N_{h}}$|⁠, the computation of |${\mathbb{N}}{\mathbf{v}}$| consists of three successive steps |${\mathbf{v}}_{1}={\mathbb{K}}{\mathbf{v}}$|⁠, |${\mathbf{v}}_{2}={\mathbb{F}}{\mathbf{v}}_{1}$| and |${\mathbb{N}}{\mathbf{v}}={\mathbb{K}}^\top{\mathbf{v}}_{2}$|⁠. The total number of operations is |$O(N_{h}N_{H}|\!\ln H|^{3/2})$| too. Therefore, the computational complexity of dealing with |${\mathbb{N}}$| needs |$O\big (N_{h}^{3/2}|\!\ln h|^{3/2}\big )$| operations, provided that |$h=O(H^{2})$|⁠, and is very time-consuming.

5.3 Inexact Newton method

To further reduce the computational complexity, we truncate the sum |$n_{H}^{\prime}$| in (3.11) and define an approximate operator |$\tilde n_{H}^{\prime}\approx n_{H}^{\prime}$| as follows

From (3.6), we know that |$L_{H}:=L_{H}[\tilde V_{h}^{(k)}]\le C|\!\ln H|^{3/2}$|⁠. Usually, an efficient algorithm for computing the eigen-pairs |$(\varepsilon _{i,H},\psi _{i,H})$|⁠, |$1\le i\le L_{H}$|⁠, needs |$O(N_{H}L_{H}^{2})=O(N_{H}|\!\ln H|^{3})$| operations. Similar to |${\mathbb{N}}$|⁠, the stiffness matrix associated with |$\big (\tilde n^{\prime}_{H}[\tilde V_{h}^{(k)}](\tilde r_{k}), v_{h}\big )$| can be written as |$\tilde{\mathbb{N}} = \tilde{\mathbb{K}}^\top \tilde{\mathbb{F}}\tilde{\mathbb{K}}$|⁠, where |$\tilde{\mathbb{F}}$| is a diagonal matrix and |$\tilde{\mathbb{K}}$| is a low-rank matrix. Their nonzero entries are given by

where |$1\le i,j\le L_{H}$| and |$1\le l\le N_{H}$|⁠. The build of |$\tilde{\mathbb{K}}$| and |$\tilde{\mathbb{F}}$| needs |$O(N_{h}|\!\ln H|^{3})$| operations totally, and the computation of |$\tilde{\mathbb{N}}{\mathbf{v}}$| also needs |$O(N_{h}|\!\ln H|^{3})$| operations for any |$\mathbf{v}\in{\mathbb{R}}^{N_{h}}$|⁠. Therefore, the computational complexity of dealing with |$\tilde{\mathbb{N}}$| only needs |$O(N_{h}|\!\ln H|^{3})$| operations. Comparing |$\tilde{\mathbb{N}}$| with matrix |${\mathbb{N}}$| in the Newton–Raphson method, the computational complexity is reduced from |$O(N_{h}^{3/2}|\!\ln h|^{3/2})$| to |$O(N_{h}|\!\ln H|^{3})$|⁠. Tables 2 and 3 in the next section show the efficiencies of both the Newton–Raphson method and the inexact Newton method.

6. Numerical experiments

We report two numerical experiments to demonstrate the competitive behavior of the two-grid finite element method. The first one verifies the convergence rate of the two-grid method. The second one demonstrates the efficiency of the inexact Newton method.

Our code is based on the finite element toolbox PHG developed in the State Key Laboratory of Scientific and Engineering Computing (LSEC), Chinese Academy of Sciences. The numerical experiments are carried out on the super computer LSSC-IV of LSEC. The computational domain is set by |$\varOmega =(0,1)^{3}$| and the distribution function is chosen as |$f(t)=f_{0}\,e^{-t/10}$|⁠.

Throughout this section, we fix the number of computed eigenvalues by |$L_{H}=10$|⁠. Computing more eigenvalues does not improve the results due to the fact that |$f_{H}(\varepsilon _{l,H}[u])=0$| for |$l>L_{H}$|⁠. In practice, we can also compute the eigenvalues |$\varepsilon _{l,H}$| one by one from (3.5) until |$f_{H}(\varepsilon _{l,H}[u])=0$|⁠. Then the number |$L_{H}$| is set by |$l-1$|⁠.

Example 6.1.

The operator has a sequence of eigenvalues |$\lambda _{n_{x},n_{y},n_{z}}$| and corresponding functions |$\phi _{n_{x},n_{y},n_{z}}$|⁠, which are given by

$$ \begin{eqnarray*} \begin{cases} \lambda_{n_{x},n_{y},n_{z}}=(n^{2}_{x}+n^{2}_{y}+n^{2}_{z})\pi^{2},\\ \phi_{n_{x},n_{y},n_{z}}=8\pi^{-3}\sin(n_{x} \pi x)\sin(n_{y}\pi y)\sin(n_{z}\pi z), \end{cases} \quad n_{x},n_{y},n_{z}\in{\mathbb{Z}}_{+}. \end{eqnarray*} $$

We arrange the eigenvalues in ascending order and replace triple indices |$(n_{x},n_{y},n_{z})$| by single indices such that the eigenvalues are written as |$\lambda _{1}\leq \lambda _{2}\leq \lambda _{3}\leq \cdots $|⁠, and the corresponding eigenfunctions are denoted by |$\phi _{1}, \phi _{2}, \phi _{3},\cdots $|⁠.

Define |$V_{0}= \sin (\pi x)\sin (\pi y)\sin (\pi z)$|⁠, |$f_{0}=1$|⁠, |$\varepsilon _{F}=0$| and |$n_{D}=\sum \limits _{l=1}^\infty f(\lambda _{l}-\varepsilon _{F})\phi _{l}^{2} -\varDelta V_{0}$|⁠. The exact solution of the problem is given by |$V=V_{0}$|⁠.

To illustrate the convergence rates of the numerical solutions, we choose a tetrahedral mesh for the calculations. The fine mesh is obtained by uniform refinements of the coarse mesh. Table 1 shows that optimal convergence rates are obtained for the numerical potential when |$h\leq 2H^{2}$|⁠, that is, ⁠.

Table 1

Open in new tab

Convergence orders of |$V_{h}$|⁠. (Example 6.1)

|$N_{h}$| .	|$h$| .	|$N_{H}$| .	|$H$| .	.	order .
|$125$|	|$1/4$|	|$125$|	|$1/4$|	4.64e−01	—
|$729$|	|$1/8$|	|$125$|	|$1/4$|	1.89e−01	1.29
|$4913$|	|$1/16$|	|$729$|	|$1/8$|	8.25e−02	1.15
|$35937$|	|$1/32$|	|$729$|	|$1/8$|	3.99e−02	1.05
|$274625$|	|$1/64$|	|$4913$|	|$1/16$|	1.95e−02	1.03

|$N_{h}$| .	|$h$| .	|$N_{H}$| .	|$H$| .	.	order .
|$125$|	|$1/4$|	|$125$|	|$1/4$|	4.64e−01	—
|$729$|	|$1/8$|	|$125$|	|$1/4$|	1.89e−01	1.29
|$4913$|	|$1/16$|	|$729$|	|$1/8$|	8.25e−02	1.15
|$35937$|	|$1/32$|	|$729$|	|$1/8$|	3.99e−02	1.05
|$274625$|	|$1/64$|	|$4913$|	|$1/16$|	1.95e−02	1.03

Table 1

Open in new tab

Convergence orders of |$V_{h}$|⁠. (Example 6.1)

|$N_{h}$| .	|$h$| .	|$N_{H}$| .	|$H$| .	.	order .
|$125$|	|$1/4$|	|$125$|	|$1/4$|	4.64e−01	—
|$729$|	|$1/8$|	|$125$|	|$1/4$|	1.89e−01	1.29
|$4913$|	|$1/16$|	|$729$|	|$1/8$|	8.25e−02	1.15
|$35937$|	|$1/32$|	|$729$|	|$1/8$|	3.99e−02	1.05
|$274625$|	|$1/64$|	|$4913$|	|$1/16$|	1.95e−02	1.03

|$N_{h}$| .	|$h$| .	|$N_{H}$| .	|$H$| .	.	order .
|$125$|	|$1/4$|	|$125$|	|$1/4$|	4.64e−01	—
|$729$|	|$1/8$|	|$125$|	|$1/4$|	1.89e−01	1.29
|$4913$|	|$1/16$|	|$729$|	|$1/8$|	8.25e−02	1.15
|$35937$|	|$1/32$|	|$729$|	|$1/8$|	3.99e−02	1.05
|$274625$|	|$1/64$|	|$4913$|	|$1/16$|	1.95e−02	1.03

The information about both fine meshes and coarse meshes are listed in Table 1. As a comparison, we also present the numbers of fixed-point iterations for solving the nonlinear discrete Poisson problem. The tolerance for relative residuals is set by |$10^{-8}$|⁠.

Table 2 lists the numbers of iterations required for solving the coupled SP system with both the fixed-point iterative method and the inexact Newton method. Clearly the inexact Newton method needs much less iterations than the fixed-point iterative method and is robust as the number of DOFs increases.

In Fig. 1, we show the numbers of iterations for both the inexact Newton method and the fixed-point iterative method on a fine mesh with |$N_{h}=729$| and a coarse mesh with |$N_{H}=125$|⁠. We observe that the inexact Newton method has second-order convergence, while the fixed-point iterative method on a single mesh has first-order convergence.

Fig. 1.

Fixed-point iterative method versus inexact Newton method.

Open in new tabDownload slide

Table 3 lists the computing time for solving the problem by different methods. We find that the two-grid inexact Newton method takes less time than other methods. Moreover, the two Newton-type methods are much superior to the fixed-point method in terms of computing time.

Funding

Tao Cui was supported by National Key R & D Program of China 2019YFA0709600 and 2019YFA0709602. Weiying Zheng was supported in part by National Key R & D Program of China 2019YFA0709600 and 2019YFA0709602, and the National Science Fund for Distinguished Young Scholars 11725106.

References