# The topological complexity of the ordered configuration space of disks in a
strip

Nicholas Wawrykow

###### Abstract

How hard is it to program $n$ robots to move about a long narrow aisle such
that only $w$ of them can fit across the width of the aisle? In this paper, we
answer that question by calculating the topological complexity of
$\text{conf}(n,w)$, the ordered configuration space of open unit-diameter
disks in the infinite strip of width $w$. By studying its cohomology ring, we
prove that, as long as $n$ is greater than $w$, the topological complexity of
$\text{conf}(n,w)$ is $2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$, providing a
lower bound for the minimum number of cases such a program must consider.

##  1 Introduction

How hard is it to program a fleet of $n$ robots to move about a long narrow
aisle such that only $w$ of the robots can fit across the width of the aisle?
If we are given the initial and final configurations of the robots, it is
rather easy to tell the robots how to move without bumping into each other or
the walls. However, if we seek a universal program that can guide the robots
between every possible combination of initial and final configurations and
avoids collisions, the programing task becomes far more daunting. We could
write a program that considers every possible relative configuration of the
robots, but as $n$ grows the number of cases quickly gets out of hand. Instead
of dealing with such a task, we would like to develop a program that breaks
this problem into as few cases as possible.

In this _motion planning problem_ we seek to minimize the number of different
cases of starting and ending positions of the fleet of robots for which we can
find a continuous rule, or _motion planner_. In the case of $n$ robots in a
long aisle such that at most $w$ of the robots fit across the width of the
aisle, this problem is equivalent to determining the _topological complexity,_
denoted TC, of the _ordered configuration space of $n$ open unit-diameter
disks in the infinite strip of width $w$_, see Figure 1. This space,
abbreviated $\text{conf}(n,w)$, is the space of ways of embedding $n$ open
unit-diameter disks in the infinite strip of width $w$, i.e., it is the
following subspace of $\mathbb{R}^{2n}$

$\text{conf}(n,w):=\big{\\{}(x_{1},y_{1},\dots,x_{n},y_{n})\in\mathbb{R}^{2n}|(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}\geq
1\text{ for }i\neq j\text{ and }\frac{1}{2}\leq y_{i}\leq w-\frac{1}{2}\text{
for all }i\big{\\}},$

and its topology as has been studied in [AKM21, AM24, BBK14, Waw22, Waw24]
among others. It is a natural restriction of the ordered configuration space
of $n$ points in the plane, $F_{n}(\mathbb{R}^{2})$, which Farber and
Yuzvinsky proved has topological complexity $2n-2$ if $n>1$ [FY02, Theorem 1].

In this paper, we determine $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$ by
studying the cohomology ring of $\text{conf}(n,w)$. This gives us a lower
bound for $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$. This, combined with an
identical upper bound arising from a space homotopy equivalent to
$\text{conf}(n,w)$, yields the following result.

######  Theorem 1.1.

$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}=\begin{cases}1&\mbox{if }n=1,\\\
2n-2&\mbox{if }1<n\leq w,\\\
2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1&\mbox{if }n>w.\end{cases}$

It follows that, as long as we have sufficiently many robots, only $w$ of
which fit across the aisle, we could never write a program that uses fewer
than $2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$ cases for the pairs of
initial and final configurations. Moreover, if we are very clever, we could
write a program with exactly $2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$
rules.

Figure 1: A point in $\text{conf}(5,2)$. Note that two disks can be aligned
vertically, but three cannot.

###  1.1 Acknowledgements

The author would like to thank Benson Farb for comments on an earlier draft of
this paper.

##  2 Topological complexity

Given a space $X$, there is a natural motion planning problem: one would like
a motion planner that takes in pair of points $x_{1},x_{2}\in X$ and produces
a continuous path from $x_{1}$ to $x_{2}$. For example, given two points in
$\text{conf}(n,w)$, i.e., two configurations of ordered disks, we would like a
rule that tells us how to move from one point to another, i.e., how to move
the disks in the first configuration to get the second configuration, keeping
the disks in the strip and avoiding collisions. Moreover, small perturbations
to $x_{1}$ or $x_{2}$ should not change the path drastically, i.e., this
motion planner should be continuous. In the case of $\text{conf}(n,w)$, this
means that configurations near to the initial and final ones should lead to
similar movements of the disks.

We formalize the general problem by letting $X^{I}$ denote the space of all
continuous paths

$\gamma:[0,1]\to X,$

with the compact-open topology, and write

$\displaystyle\pi:X^{I}$ $\displaystyle\to X\times X$
$\displaystyle\pi(\gamma)$ $\displaystyle=\big{(}\gamma(0),\gamma(1)\big{)}$

for the map $\pi$ that associates a path $\gamma\in X^{I}$ to its endpoints.
In this setting the problem of finding a continuous motion planner is
equivalent to finding a continuous section of $\pi$, i.e., a map

$s:X\times X\to X^{I},$

such that $\pi\circ s=\text{id}_{X\times X}$. If $X$ is convex, such a motion
planner exists, namely $s$ takes the points $x_{1}$ and $x_{2}$ to the
constant speed straight-line path between them, i.e.,

$s(x_{1},x_{2})(t)=tx_{1}+(1-t)x_{2}.$

Unfortunately, continuous motion planners $s$ only exist for slightly more
general $X$, as Farber proved that such an $s$ exists if and only if $X$ is
contractible [Far03, Theorem 1].

Since only the simplest spaces have continuous motion planners, we change our
question and ask how bad can the motion planning problem really be? Instead of
seeking a single motion planner $s$ for all of $X\times X$, we seek to find
the minimal number $k$ of open subsets $U_{i}\subset X\times X$ covering $X$,
such that there are continuous sections

$s_{i}:U_{i}\to X^{I},$

such that $\pi\circ s_{i}=\text{id}_{X\times X}|_{U_{i}}$. This number $k$ is
the _topological complexity of $X$_, and we say $\textbf{TC}(X)=k$. If no such
$k$ exists, then we set $\textbf{TC}(X)=\infty$.

######  Example 2.1.

The configuration space of a robot arm consisting of $n$ bars
$L_{1},\dots,L_{n}$, such that $L_{i}$ and $L_{i+1}$ are connected a flexible
joint can be viewed as the $n$-fold product of spheres, $S^{1}$s if the joints
are hinges with one dimension of freedom and $S^{2}$s if the joints are balls
and sockets with two dimensions of freedom. In the hinge case the topology
complexity is $n+1$, whereas in the ball and socket case the topological
complexity is $2n+1$ [Far03, Theorem 12].

######  Example 2.2.

If we shrink our robots from disks to points, we move from $\text{conf}(n,w)$
to $F_{n}(\mathbb{R}^{2})$, the ordered configuration space of $n$ points in
the plane. Farber and Yuzvinksy proved that
$\textbf{TC}\big{(}F_{n}(\mathbb{R}^{2})\big{)}=2n-2$ [FY02, Theorem 1], and
Farber and Grant extended this to ordered configurations in arbitrary
$\mathbb{R}^{m}$, proving that
$\textbf{TC}\big{(}F_{n}(\mathbb{R}^{m})\big{)}=2n-1$ is $m$ is even and
$\textbf{TC}\big{(}F_{n}(\mathbb{R}^{m})\big{)}=2n-2$ if $m$ is odd [FG09,
Theorem 1].

Often the space $X$ is hard to work with directly, and this makes computing
$\textbf{TC}(X)$ a challenge. Fortunately, Farber proved that topological
complexity is a homotopy invariant of $X$, [Far03, Theorem 3]. This result
will allow us to replace $\text{conf}(n,w)$ with a finite cellular complex
called $\text{cell}(n,w)$ that is significantly easier to manipulate.

######  Proposition 2.3.

(Farber [Far03, Theorem 3]) Topological complexity is a homotopy invariant.

Even though we can often replace $X$ with a simpler space $Y$, finding open
sets $\\{V_{1},\cdots,V_{k}\\}$ of $Y\times Y$ with continuous sections
$s_{i}:V_{i}\to Y^{I}$ with $k$ relatively small is a challenge. Proving that
such a $k$ is minimal is even harder. As such we seek upper and lower bounds
for topological complexity. Farber gave such bounds in [Far03, Sections 3 and
4], and we will show that they are sharp in the case of $\text{conf}(n,w)$
when $n>w$.

###  2.1 An upper bound for topological complexity

We seek an upper bound for the topological complexity of $X$. As long as $X$
is path connected and paracompact such a bound arises from the dimension of
$X$.

######  Proposition 2.4.

(Farber [Far03, Theorem 4]) For any path-connected paracompact space $X$,

$\textbf{TC}(X)\leq 2\dim X+1.$

We may naturally view $\text{conf}(n,w)$ a subspace of $\mathbb{R}^{2n}$,
proving that

$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}\leq 4n+1.$

This bound is sub-optimal. As we will see in the next section,
$\text{conf}(n,w)$ is homotopy equivalent to an
$\big{(}n-\big{\lceil}\frac{n}{w}\big{\rceil}\big{)}$-dimensional cellular
complex. It follows that Proposition 2.3 proves that
$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$ is at most
$2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$. We will show that this bound is
sharp when $w<n$.

###  2.2 A lower bound for topological complexity

Finding meaningful lower bounds for the topological complexity of a space,
$X$, is a harder problem, as this corresponds to showing that _any_ motion
planning program must consider a certain number of cases. Fortunately, the cup
product structure on the cohomology ring of $X$ leads to a reasonable lower
bound that will be sharp in our case.

Given a space $X$, let $\Delta_{X}\subset X\times X$ denote the diagonal,
i.e., the set of points $\big{\\{}(x,x)\in X\times X\big{\\}}$. A cohomology
class

$\phi\in H^{*}(X\times X;R)$

is called a _zero-divisor_ if its restriction to the diagonal is trivial,
i.e.,

$\phi|_{\Delta_{X}}=0\in H^{*}(X;R).$

We say that the _zero-divisor-cup-length_ of $H^{*}(X;R)$ is the length of the
longest non-trivial cup product $\big{(}$in $H^{*}(X\times X;R)\big{)}$ of
zero-divisors. Farber proved that the zero-divisor-cup-length serves as a
lower bound for the topological complexity of $X$.

######  Proposition 2.5.

(Farber [Far03, Theorem 7]) The zero-divisor-cup-length of $H^{*}(X;R)$ is
less than the topological complexity of $X$.

Farber also noted that if $\alpha\in H^{*}(X;R)$, then

$\overline{\alpha}:=1\otimes\alpha-\alpha\otimes 1\in H^{*}(X;R)\otimes
H^{*}(X;R)\cong H^{*}(X\times X;R)$

is a zero-divisor, where the isomorphism $H^{*}(X;R)\otimes H^{*}(X;R)\cong
H^{*}(X\times X;R)$ is given by the Künneth theorem. We will use this
construction to find a family of $2n-2\big{\lceil}\frac{n}{w}\big{\rceil}$
zero-divisors in $H^{*}\big{(}\text{conf}(n,w)\times\text{conf}(n,w)\big{)}$
that have non-trivial cup product, yielding a lower bound for the topological
complexity of $\text{conf}(n,w)$.

##  3 Bounds for $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$

In this section we determine sharp upper and lower bounds for the topological
complexity of $\text{conf}(n,w)$. We do this by recalling the definitions of
two spaces homotopy equivalent to $\text{conf}(n,w)$ that will allow us to use
Propositions 2.4 and 2.5 to bound the topological complexity of
$\text{conf}(n,w)$. The first space is the cellular complex
$\text{cell}(n,w)$, whose dimension can easily be determined, yielding an
upper bound for $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$. The second space,
denoted $M(n,w)$, will allow us to compute certain cup products in
$H^{*}\big{(}\text{conf}(n,w)\big{)}$, which we will use to find a lower bound
for $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$. We begin by stating several
facts about $\text{cell}(n,w)$. In the next subsection we will recall Alpert
and Manin’s construction of a discrete gradient vector field on
$\text{cell}(n,w)$, the critical cells of which correspond to a basis for the
cohomology of $M(n,w)$.

###  3.1 $\text{cell}(n,w)$ and an upper bound for
$\textbf{TC}\big{(}\text{conf}(n,2)\big{)}$

In this subsection we recall the definition of the cellular complex
$\text{cell}(n,w)$, which was introduced in [BZ14], and is homotopy equivalent
to $\text{conf}(n,w)$.

Let _$\text{cell}(n)$_ denote the cellular complex whose cells are represented
by _symbols_ consisting of an ordering of the numbers $1,\dots,n$ separated by
vertical bars into _blocks_ such that no block is empty. A cell
$f\in\text{cell}(n)$ is a codimensional-$1$ _face_ of a cell
$g\in\text{cell}(n)$ if $g$ can be obtained by deleting a bar in $f$ and
shuffling the resulting block. We say that $g$ is a _co-face_ of $f$. It
follows that cells represented by symbols with no bars are the top-dimensional
cells of $\text{cell}(n)$. It follows that we can view $\text{cell}(n)$ as an
$(n-1)$-dimensional complex, and each bar in a symbol lowers the dimension of
the corresponding cell by $1$.

By restricting how big a block can be in $\text{cell}(n)$, one gets the
cellular complex $\text{cell}(n,w)$. The cellular complex _$\text{cell}(n,w)$_
is the subcomplex of $\text{cell}(n)$ consisting of the cells represented by
symbols whose blocks have at most $w$ elements; see Figure 2. Alpert, Kahle,
and MacPherson proved that $\text{cell}(n,w)$ is homotopy equivalent to
$\text{conf}(n,w)$ [AKM21, Theorem 3.1], allowing us to study it in place of
$\text{conf}(n,w)$.

######  Proposition 3.1.

(Alpert–Kahle–MacPherson [AKM21, Theorem 3.1]) There is an $S_{n}$-equivariant
homotopy equivalence $\text{conf}(n,w)\simeq\text{cell}(n,w)$.

It follows from Proposition 2.3 that if we can determine the dimension of
$\text{cell}(n,w)$ we will get an upper bound for the topological complexity
of $\text{conf}(n,w)$.

######  Proposition 3.2.

The dimension of $\text{cell}(n,w)$ is
$n-\big{\lceil}\frac{n}{w}\big{\rceil}$.

###### Proof.

We seek to find the cells represented by the symbols with the least number of
vertical bars. The symbol

$1\,\cdots\,w\Big{|}w+1\,\cdots\,2w\Big{|}\cdots\Big{|}\Big{\lceil}\frac{n}{w}\Big{\rceil}-n+1\,\cdots\,n$

corresponds to a top dimensional cell of $\text{cell}(n,w)$ as it has the
fewest possible vertical bars, as any fewer bars would result in a block of
size greater than $w$ and every block, perhaps other than the last, has $w$
elements. This cell is
$\big{(}n-\big{\lceil}\frac{n}{w}\big{\rceil}\big{)}$-dimensional as it has
$\big{\lceil}\frac{n}{w}\big{\rceil}-1$ bars. ∎

This gives improves our upper bound on the topological complexity of
$\text{conf}(n,w)$.

######  Proposition 3.3.

$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}\leq
2n-2\Big{\lceil}\frac{n}{w}\Big{\rceil}+1.$

###### Proof.

This follows immediately from Propositions 2.4 and 3.2. ∎

Next, we find a lower bound the topological complexity of $\text{conf}(n,w)$.
Our bound will coincide with that of Proposition 3.3, proving that they are
sharp.

###  3.2 $M(n,w)$ and a lower bound for
$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$

Having found an upper bound for $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$ we
seek a suitable lower bound, i.e., the minimal number of cases any program
controlling the movement of the $n$ robots about the aisle must contain. In
order to use Proposition 2.5 we need to understand the cohomology ring of
$\text{conf}(n,w)$. Alpert and Manin construct a discrete gradient vector
field on $\text{cell}(n,w)$. The critical $j$-cells of this vector field will
correspond to a basis for $H^{j}\big{(}\text{conf}(n,w)\big{)}$, via
Poincaré–Lefschetz duality applied to another space called $M(n,w)$ that is
homotopy equivalent to $\text{conf}(n,w)$.

We begin by describing the critical cells of Alpert and Manin’s discrete
gradient vector field as they will correspond to elements in the cohomology of
$\text{conf}(n,w)$; see [AM24, Sections 3 and 4] for further details. Given a
cell $f\in\text{cell}(n,w)$, we break $f$ into its blocks. Each block is
further decomposed into _wheels_ , where each entry of a block is the _axle_
of a wheel if it is the largest entry of the block up to that point, and the
wheel consists of the axle and all the following smaller entries before the
next axle. A block $f_{i}$ of $f$ is a _unicycle_ if it consists of a single
wheel, i.e., the largest element in $f_{i}$ is also the first element. We say
that $f_{i+1}$ is a _follower_ if $f_{i}$ is a unicycle (its _leader_), which
is not itself a follower and whose wheel is smaller than any of the follower’s
wheels. The cell $f$ is critical if every block $f_{i}$ is either

  * •

a unicycle which is not a follower, or

  * •

a follower such that it and its leader have at least $w+1$ elements in total.

See Figure 2 for an example.

Figure 2: The cellular complex $\text{cell}(3,2)$ with Alpert and Manin’s
discrete gradient vector field. The critical cells of this vector field are in
red. This vector field arises from a total ordering of the cells of
$\text{cell}(3,2)$, and a vector $[f,g]$ is included if and only if $f$ is the
greatest face of $g$ and $g$ is the least coface of $f$. For more on Alpert
and Manin’s ordering see [AM24, Sections 3 and 4].

We relate these critical cells to the cohomology ring of $\text{conf}(n,w)$ by
studying another space, $M(n,w)$, that will allow us to take advantage of
Poincaré–Lefschetz duality. The _follower-free_ critical cells of
$\text{cell}(n,w)$, i.e., those such that no block is a follower, will be of
special interest to us, as the corresponding basis elements behave especially
nicely under taking cup products.

Let $\emph{M(n,w)}\subset\mathbb{R}^{2n}$ be the configuration space of open
disks of radius $1$ in a strip of any finite length $N>n$ and width
$w+\epsilon$ for any $0<\epsilon<1$. This space is a manifold (with boundary)
of dimension $2n$, and it homotopy equivalent to $\text{conf}(n,w)$. It
follows from Poincaré–Lefschetz duality that

$H^{i}\big{(}\text{conf}(n,w)\big{)}\cong H^{i}\big{(}M(n,w)\big{)}\cong
H_{2n-i}\big{(}M(n,w),\partial M(n,w)\big{)},$

and the cup product between classes in $H^{i}\big{(}M(n,w)\big{)}$ and
$H^{j}\big{(}M(n,w)\big{)}$ is given by the transverse intersection map

$\pitchfork H_{2n-i}\big{(}M(n,w),\partial M(n,w)\big{)}\otimes
H_{2n-j}\big{(}M(n,w),\partial M(n,w)\big{)}\to
H_{2n-i-j}\big{(}M(n,w),\partial M(n,w)\big{)}.$

Alpert and Manin describe a basis for $H_{2n-j}\big{(}M(n,w),\partial
M(n,w)\big{)}$, hence a basis for $H^{j}\big{(}\text{conf}(n,w)\big{)}$, by
associating basis elements to critical $j$-cells of $\text{cell}(n,w)$. Let
$f$ be the label of a critical cell, then the submanifold $V(f)\subset M(n,w)$
is the set of disk configurations such that

  1. 1.

The disks in each block of $f$ are lined up vertically in order.

  2. 2.

If a block $b_{1}$ comes before $b_{2}$ in $f$, they have a least $w+1$
elements combined, and one of them is a follower, then the column of disks
labeled by elements of $b_{1}$ is to the left of the column of disks labeled
by elements of $b_{2}$.

See Figure 3 for an example.

Figure 3: A point in $V(14\,10\,9\,8|13\,7\,6\,5|12\,4\,3\,2|11\,1)$. Note
that the disks may wiggle up and down while remaining aligned. The columns of
the disks are allowed to move left and right, and in this submanifold the
order of the columns can be permuted, so this point is one of the twenty four
components of $V(14\,10\,9\,8|13\,7\,6\,5|12\,4\,3\,2|11\,1)$.

Clearly this is a $(2n-j)$-dimensional submanifold of $M(n,w)$, and one can
check that $\partial V(f)=V(f)\cap\partial M(n,w)$, proving that it represents
a class in $H_{2n-j}\big{(}M(n,w),\partial M(n,w)\big{)}$. Lemma 6.1 of Alpert
and Manin proves that the $V(f)$ are in fact a basis for
$H_{2n-j}\big{(}M(n,w),\partial M(n,w)\big{)}$ [AM24, Lemma 6.1], and we will
write $\nu(f)$ for the Poincaré–Lefschetz dual of $V(f)$. Additionally, note
that the follower-free critical cells correspond to submanifolds in which the
columns of disks can be arranged in any order.

Given Proposition 2.5, we wish to understand the cup product structure of
$H^{*}\big{(}\text{conf}(n,w)\big{)}$. First we recall a result of Alpert and
Manin that will simplify our computations.

######  Proposition 3.4.

(Alpert–Manin [AM24, Proposition 6.2]) If there is a pair of labels $i$ and
$j$ which are contained in the same block of critical cells $f$ and $g$, then
$\nu(f)\cup\nu(g)=0$.

Next, we note a fact that will greatly reduce the number of cup products we
will need to take to determine whether a cup product is nonzero. This result,
along with the previous one, will allow us to use the identification
$H^{2n-2\big{\lceil}\frac{n}{w}\big{\rceil}}\big{(}\text{conf}(n,w)\times\text{conf}(n,w)\big{)}\cong
H^{n-\big{\lceil}\frac{n}{w}\big{\rceil}}\big{(}\text{conf}(n,w)\big{)}\otimes
H^{n-\big{\lceil}\frac{n}{w}\big{\rceil}}\big{(}\text{conf}(n,w)\big{)}$.

######  Proposition 3.5.

For $k\geq w$ and $i>j_{l}$ where $1\leq l\leq k$, let

$a_{l}=i\,j_{l}|n|\cdots|\widehat{i}|\cdots|\widehat{j_{l}}|\cdots|1$

be a follower-free critical cell of $\text{cell}(n,w)$. Then,

$\prod\nu(a_{l})=0.$

###### Proof.

Consider the intersection of the $V(a_{l})$. In such an intersection the disk
labeled $i$ is aligned above the disks labeled $j_{1},\dots,j_{k}$. Since
$k\geq w$, this means that $k+1>w$ disks must be aligned vertically in the
intersection, which is impossible. Therefore, the intersection is empty and
the cup product is $0$. ∎

The following two propositions will allow us to determine what terms might
appear when we take the cup product of $n-\big{\lceil}\frac{n}{w}\big{\rceil}$
classes in $H^{1}\big{(}\text{conf}(n,w)\big{)}$.

######  Proposition 3.6.

Let

$f=i\,i_{1}\,\cdots\,i_{k}|n|\cdots|\widehat{i}|\cdots|\widehat{i_{l}}|\cdots|1$

be a follower-free critical $k$-cell of $\text{cell}(n,w)$ such that $k+1<w$,
and let

$a=i\,j|n|\cdots|\widehat{i}|\cdots|\widehat{j}|\cdots|1$

be a follower-free critical $1$-cell of $\text{cell}(n,w)$ such that $j\neq
i_{1},\dots,i_{k}$. Then,

$\nu(f)\cup\nu(a)\sum\nu(e_{l}),$

where for $0\leq l\leq k$,

$e_{l}=i\,i_{1}\,\cdots\,i_{l}\,j\,i_{l+1}\,\cdots\,i_{k}|n|\cdots|\widehat{i}|\cdots|\widehat{i_{l}}|\cdots|1.$

is a follower-free critical $(k+1)$-cell of $\text{cell}(n,w)$.

###### Proof.

The intersection of $V(f)$ and $V(a)$ is transverse. It consists of all
configurations such that the disk labeled $i$ is aligned above the disks
labeled $j,i_{1},\dots,i_{k}$, and the disk labeled $i_{l}$ is above the disk
labeled $i_{l+1}$. This exactly the submanifold $\sqcup V(e_{l})$, completing
the proof. ∎

######  Proposition 3.7.

Let

$f=i_{1,1}\,\cdots\,i_{1,n_{1}}|\cdots|i_{j,1}\,\cdots\,i_{j,n_{j}}|i_{n_{1}+\cdots+n_{j}+1}|\cdots|i_{n}$

be a follower-free critical cell of $\text{cell}(n,w)$, and let

$g=i_{j+1,1}\,\cdots\,i_{j+1,n_{j+1}}|i^{\prime}_{n_{j+1}}|\cdots|i^{\prime}_{n}$

be a follower-free critical cell of $\text{cell}(n,w)$, such that
$i_{j+1,k}\in\\{i_{n_{1}+\cdots+n_{j}+1},\dots,i_{n}\\}$ for $1\leq k\leq
n_{j+1}$. Then

$\nu(f)\cup\nu(g)=\nu(e),$

where

$e=i_{1,1}\,\cdots\,i_{1,n_{1}}|\cdots|i_{j+1,1}\,\cdots\,i_{j+1,n_{j+1}}|\cdots|i_{j,1}\,\cdots\,i_{j,n_{j}}|i_{n_{1}+\cdots+n_{j}+1}|\cdots|\widehat{i_{j+1,1}}|\cdots|\widehat{i_{j+1,n_{j+1}}}|\cdots|i_{n}$

is a follower-free critical cell of $\text{cell}(n,w)$, i.e., there is an
$1\leq l\leq j$ such that $n_{l}\leq n_{j+1}\leq n_{l+1}$ and if
$n_{l}=n_{j+1}$, then $i_{1,l}>i_{1,j+1}$, and if $n_{j+1}=n_{l+1}$, then
$i_{1,j+1}>i_{1,l+1}$.

###### Proof.

It is clear that $V(f)$ and $V(g)$ have transverse intersection. Moreover,
their intersection is the submanifold of $V(f)$ such that the disks
$i_{j+1,1},\dots,i_{j+1,n_{j+1}}$ are aligned vertically in that order. It
follows that this column of disks can be anywhere among the columns of the
other disks of $V(f)$. This is exactly $V(e)$. ∎

Now we are ready to calculate the zero-divisor-cup-length of
$H^{*}\big{(}\text{conf}(n,w)\big{)}$ for $n>w$. Our proof relies on the
observation that the submanifold of $M(n,w)$ corresponding to $V(f)$ where $f$
is a follower-free top dimensional critical cell of $\text{cell}(n,w)$ such
that the $\big{\lceil}\frac{n}{w}\big{\rceil}$ blocks of $f$ have first
element in
$\big{\\{}n,\dots,n-\big{\lceil}\frac{n}{w}\big{\rceil}+1\big{\\}}$, see
Figure 3, is a submanifold contained in the intersection of
$n-\big{\lceil}\frac{n}{w}\big{\rceil}$ submanifolds arising from follower-
free critical $1$-cells of $\text{cell}(n,w)$ such that the first element of
the block with two elements is in
$\big{\\{}n,\dots,n-\big{\lceil}\frac{n}{w}\big{\rceil}+1\big{\\}}$, and the
second element is in
$\big{\\{}n-\big{\lceil}\frac{n}{w}\big{\rceil},\dots,1\big{\\}}$. If we
carefully choose two disjoint sets $A$ and $B$ of
$n-\big{\lceil}\frac{n}{w}\big{\rceil}$ such follower-free critical $1$-cells,
we will get a non-trivial product in
$H^{*}\big{(}\text{conf}(n,w)\otimes\text{conf}(n,w)\big{)}$.

######  Lemma 3.8.

For $n>w$ the zero-divisor-cup-length of $H^{*}\big{(}\text{conf}(n,w)\big{)}$
is $2n-2\big{\lceil}\frac{n}{w}\big{\rceil}$.

###### Proof.

Set $m:=\big{\lceil}\frac{n}{w}\big{\rceil}$. For $0\leq i\leq m-2$ and $0\leq
j\leq w-2$, let

$a_{i,j}=(n-i)\,\big{(}n-m-(w-1)i-j\big{)}|n|\cdots|\widehat{(n-i)}|\cdots|\widehat{\big{(}n-m-(w-1)i-j\big{)}}|\cdots|1$

and, for $i=m-1$ and $0\leq j\leq(n\text{ mod}w)-2$, let

$a_{i,j}=(n-m+1)\,\big{(}n-m-(w-1)(m-1)-j\big{)}|n|\cdots|\widehat{(n-m+1)}|\cdots|\widehat{\big{(}n-m-(w-1)(m-1)-j\big{)}}|\cdots|1$

be critical $1$-cells of $\text{cell}(n,w)$. Similarly, for $0\leq i\leq m-2$
and $0\leq j\leq w-2$, let

$b_{i,j}=(n-i+1)\,\big{(}n-m-(w-1)i-j\big{)}|n|\cdots|\widehat{(n-i)}|\cdots|\widehat{\big{(}n-m-(w-1)i-j\big{)}}|\cdots|1$

and, for $i=m-1$ and $0\leq j\leq(n\text{ mod}w)-2$,

$b_{i,j}=(n)\,\big{(}n-m-(w-1)(m-1)-j\big{)}|n|\cdots|\widehat{(n-m+1)}|\cdots|\widehat{\big{(}n-m-(w-1)(m-1)-j\big{)}}|\cdots|1$

be critical $1$-cells of $\text{cell}(n,w)$. There are precisely
$n-\big{\lceil}\frac{n}{w}\big{\rceil}$ of each of the $a_{i,j}$s and
$b_{i,j}$s.

Set

$\alpha_{i,j}=\nu(a_{i,j})\indent\text{and}\indent\beta_{i,j}=\nu(b_{i,j}),$

and

$\overline{\alpha_{i,j}}=\alpha_{i,j}\otimes
1-1\otimes\alpha_{i,j}\indent\text{and}\indent\overline{\beta_{i,j}}=\beta_{i,j}\otimes
1-1\otimes\beta_{i,j}.$

We claim that

$\prod\overline{\alpha_{i,j}}\cup\overline{\beta_{i,j}}\neq 0.$

To see this note that any nonzero term in the cup product in
$H^{2n-2\big{\lceil}\frac{n}{w}\big{\rceil}}\big{(}\text{conf}(n,w)\times\text{conf}(n,w)\big{)}$
must be in $H^{n-\lceil\frac{n}{w}\rceil}\big{(}\text{conf}(n,w)\big{)}\otimes
H^{n-\lceil\frac{n}{w}\rceil}\big{(}\text{conf}(n,w)\big{)}$ as
$H^{k}\big{(}\text{conf}(n,w)\big{)}=0$ for
$k>n-\big{\lceil}\frac{n}{w}\big{\rceil}$.

Consider any cup product of a total of $n-\big{\lceil}\frac{n}{w}\big{\rceil}$
of the $\alpha_{i,j}$s and $\beta_{i,j}$s. Since the cohomology ring is graded
commutative, it follows that we may first take the product of all the terms
such that $i=0$, take the product of all the terms such that $i=1$, etc, and
then take the product of the resulting classes. Note that by Proposition 3.5,
if for any $0\leq i\leq m-1$ there are more than $w-1$ total terms of the form
$\alpha_{i,j}$ or $\beta_{i,j}$, then the cup product is $0$. By Proposition
3.6 for each $0\leq i\leq m-1$ we may rewrite the product of the
$\alpha_{i,j}$s and $\beta_{i,j}$s as a sum of products of $\nu(e_{i})$, where
the $e_{i}$ are follower-free critical cells of the form

$e_{i}=n-i\,i_{i,2}\,\cdots,i_{i,l}|n|\widehat{n-i}|\cdots|\widehat{i_{i,2}}|\cdots|\widehat{i_{i,l}}|\cdots|1,$

where the $i_{2},\dots,i_{l}$ are permutations of the second elements in the
block with two elements of the critical cells of the form $a_{i,j}$ and
$b_{i,j}$ corresponding to the $\alpha_{i,j}$s and $\beta_{i,j}$s. Moreover,
if some element appears as the second element in the block with two elements
in two of the $a_{i,j}$s and $b_{i,j}$s, then the cup product is trivial. This
follows from the fact the first element of the two blocks must be different,
and since we may assume that no first entry $n-i$ appears more than $w-1$
times in product, any two first elements of a block must appear more than
$w-1$ times in the product. It follows that every non-trivial product of
$n-\big{\lceil}\frac{n}{w}\big{\rceil}$ of the $\alpha_{i,j}$s and
$\beta_{i,j}$s can be written as sum of products the $\nu(e_{i})$s where for
each product of the $\alpha_{i,j}$s and $\beta_{i,j}$s and at least $1$ of
$i$, the labels in the first block of the $e_{i}$s is different than in any
other product. By Proposition 3.7 every such product can be written as a sum
of follower-free critical
$\big{(}n-\big{\lceil}\frac{n}{w}\big{\rceil}\big{)}$-cells of
$\text{cell}(n,w)$ whose $m$ blocks are the first blocks of the $e_{i}$s for
$0\leq i\leq m-1$. Additionally, these are all basis elements. Moreover, every
nontrivial product of the $\alpha_{i,j}$s and $\beta_{i,j}$s is such that the
elements in at least one of the $m$ blocks are unique to that set of
$n-\big{\lceil}\frac{n}{w}\big{\rceil}$ of the $\alpha_{i,j}$s and
$\beta_{i,j}$s. Therefore,
$\prod\overline{\alpha_{i,j}}\cup\overline{\beta_{i,j}}\neq 0$. ∎

Combining this lemma and Proposition 2.5 yields the following lower bound for
$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}$.

######  Corollary 3.9.

$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}>2n-2\Big{\lceil}\frac{n}{w}\Big{\rceil}.$

It follows that we have a proof of Theorem 1.1, which we restate for
convenience.

######  Theorem 1.1.

$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}=\begin{cases}1&\mbox{if }n=1,\\\
2n-2&\mbox{if }1<n\leq w,\\\
2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1&\mbox{if }n>w.\end{cases}$

###### Proof.

If $n=1$, then $\text{conf}(n,w)$ is contractible, so [Far03, Theorem 1]
proves that $\textbf{TC}\big{(}\text{conf}(1,w)\big{)}=1$. If $1<n\leq w$,
then $\text{conf}(n,w)$ is homotopy equivalent to $F_{n}(\mathbb{R}^{2})$, so
[FY02, Theorem 1] proves that
$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}=2n-2$. Finally, if $n>w$, then
Corollary 3.9 proves that
$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}>2n-2\big{\lceil}\frac{n}{w}\big{\rceil}$
and Proposition 3.3 proves that $\textbf{TC}\big{(}\text{conf}(n,w)\big{)}\leq
2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$, forcing

$\textbf{TC}\big{(}\text{conf}(n,w)\big{)}=2n-2\Big{\lceil}\frac{n}{w}\Big{\rceil}+1.$

∎

It follows that any program guiding $n$ robots in the width $w$ aisle must
consider at least $2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$ cases. Moreover,
if we are very clever, we could develop a program that considers exactly
$2n-2\big{\lceil}\frac{n}{w}\big{\rceil}+1$ cases.

## References

  * [AKM21] Hannah Alpert, Matthew Kahle, and Robert MacPherson, _Configuration spaces of disks in an infinite strip_ , Journal of Applied and Computational Topology (2021), 1–34.
  * [AM24] Hannah Alpert and Fedor Manin, _Configuration spaces of disks in a strip, twisted algebras, persistence, and other stories_ , Geometry & Topology (2024).
  * [BBK14] Yuliy Baryshnikov, Peter Bubenik, and Matthew Kahle, _Min-type Morse theory for configuration spaces of hard spheres_ , International Mathematics Research Notices 2014 (2014), no. 9, 2577–2592.
  * [BZ14] Pavle V.M. Blagojević and Günter M. Ziegler, _Convex equipartitions via equivariant obstruction theory_ , Israel Journal of Mathematics 200 (2014), no. 1, 49–77.
  * [Far03] Michael Farber, _Topological complexity of motion planning_ , Discrete & Computational Geometry 29 (2003), 211–221.
  * [FG09] Michael Farber and Mark Grant, _Topological complexity of configuration spaces_ , Proceedings of the American Mathematical Society 137 (2009), no. 5, 1841–1847.
  * [FY02] Michael Farber and Sergey Yuzvinsky, _Topological robotics: subspace arrangements and collision free motion planning_ , arXiv preprint math/0210115 (2002).
  * [Waw22] Nicholas Wawrykow, _On the symmetric group action on rigid disks in a strip_ , Journal of Applied and Computational Topology 6 (2022).
  * [Waw24] Nicholas Wawrykow, _Representation stability for disks in a strip_ , Journal of Topology and Analysis (2024).


∎

11institutetext:  F. Xu 22institutetext: Institute of Computational
Mathematics, Department of Mathematics, Faculty of Science, Beijing University
of Technology, Beijing, 100124, China.
22email<EMAIL_ADDRESS>33institutetext: M. Xie 44institutetext: Center
for Applied Mathematics, Tianjin University, Tianjin 300072, China.
44email<EMAIL_ADDRESS>55institutetext: M. Yue 66institutetext: School of
Mathematics and Statistics, Beijing Technology and Business University,
Beijing 100048, China. 66email<EMAIL_ADDRESS>
# Multigrid method for nonlinear eigenvalue problems based on Newton
iteration††thanks: This work was supported by General projects of science and
technology plan of Beijing Municipal Education Commission (Grant No.
KM202110005011), National Natural Science Foundation of China (Grant Nos.
11801021).

Fei Xu  Manting Xie  Meiling Yue

(Received: date / Accepted: date)

###### Abstract

In this paper, a novel multigrid method based on Newton iteration is proposed
to solve nonlinear eigenvalue problems. Instead of handling the eigenvalue
$\lambda$ and eigenfunction $u$ separately, we treat the eigenpair
$(\lambda,u)$ as one element in a product space $\mathbb{R}\times
H_{0}^{1}(\Omega)$. Then in the presented multigrid method, only one discrete
linear boundary value problem needs to be solved for each level of the
multigrid sequence. Because we avoid solving large-scale nonlinear eigenvalue
problems directly, the overall efficiency is significantly improved. The
optimal error estimate and linear computational complexity can be derived
simultaneously. In addition, we also provide an improved multigrid method
coupled with a mixing scheme to further guarantee the convergence and
stability of the iteration scheme. More importantly, we prove convergence for
the residuals after each iteration step. For nonlinear eigenvalue problems,
such theoretical analysis is missing from the existing literatures on the
mixing iteration scheme.

###### Keywords:

Multigrid method nonlinear eigenvalue problems Newton iteration.

††journal: Journal of Scientific Computing

##  1 Introduction

Solving large-scale nonlinear eigenvalue problems is a basic and challenging
task in the field of scientific and engineering computing. Various practical
problems reduce to solve nonlinear eigenvalue problems. For example, the
Schrödinger-Newton equation which has been widely used to model the quantum
state reduction harri , the Gross-Pitaevskii equation which describes Bose-
Einstein condensates baowei , and the Kohn-Sham equation which is used to
describe the electronic structure in quantum chemistry CancesChakirMaday ;
Canceschakir . In this paper, we will study a class of nonlinear eigenvalue
models with convex energy functional. It is quite difficult to solve such
nonlinear eigenvalue problems whose computational work grows exponentially as
the problem size increases. Additionally, there are far fewer available
numerical methods for nonlinear eigenvalue problems than those for the
boundary value problems.

In this study, we resort to the multigrid method to improve the solving
efficiency for nonlinear eigenvalue problems. The multigrid method was first
proposed by Fedorenko fedore  in 1961. Later, Brandt brandt  demonstrated the
efficiency of the multigrid method, which was then a widely examined subject
in the computational mathematics field. Subsequently, Hackbusch, Xu, and many
others began to use the tool of functional analysis to analyze this algorithm
Hackbusch_Book ; Xu ; Shaidurov , which made the multigrid method experience
rapid development. As we know, the multigrid method is able to derive an
approximate solution possessing the optimal error estimates with the linear
computational complexity. For eigenvalue problems, the study on the multigrid
method is relatively limited. Xu and Zhou XuZhou_2001  first designed a two-
grid method for the linear eigenvalue problem. This algorithm requires solving
a small-scale eigenvalue problem on a coarse mesh and a large-scale boundary
value problem on a fine mesh. When the mesh sizes of the coarse mesh ($H$) and
fine mesh ($h$) have an appropriate proportional relation ($H=\sqrt{h}$), the
optimal estimate for the approximate solution can be derived. Later, Chen and
Liu et al. chenliu  extended the two-grid method to solve nonlinear eigenvalue
problems using specially selected coarse and fine meshes. However, owing to
the strict constrains on the ratio (i.e., $H=\sqrt{h}$), the two-grid method
only performs on two mesh levels and cannot be used in the multigrid sequence.

To solve the nonlinear eigenvalue problems via the multigrid technique, the
Newton iteration is adopted. In this paper, the eigenvalue $\lambda$ and
eigenfunction $u$ are treated as one element in a product space
$\mathbb{R}\times H_{0}^{1}(\Omega)$. Then, the nonlinear eigenvalue equation
can be viewed as a special nonlinear equation defined in $\mathbb{R}\times
H_{0}^{1}(\Omega)$. Next, based on the Newton iteration technique, we simply
solve a linear boundary value problem in each layer of the multigrid sequence.
Solving the large-scale nonlinear eigenvalue problem is avoided in our
algorithm. Besides, the involved linearized boundary value problems can be
solved quickly by many mature numerical algorithms, such as the classical
multigrid method. Hence, the presented algorithm can significantly improve the
solving efficiency for nonlinear eigenvalue problems. The well-posedness for
these linearized boundary value problems and the convergence for our entire
algorithm have been rigorously analyzed in this paper.

For the above-mentioned multigrid method, we find that the Newton iteration
scheme may fail to converge for some complicated nonlinear eigenvalue models.
In fact, this problem exists widely during the self-consistent field iteration
for nonlinear eigenvalue problems. To derive a convergence result, the mixing
scheme (or damping method) is introduced. The most widely used mixing scheme
is Anderson acceleration, which is used to improve the convergence rate for
fixed-point iterations. Anderson acceleration was first designed by D. G.
Anderson in 1965 to solve the integral equation andersion . Subsequently, this
technique was used to solve various models, such as, the self-consistent field
iteration in electronic structure computations fang , flow problems vc1 ; vc2
, molecular interaction stasiak , etc. Anderson acceleration adopts the mixing
scheme to make the fixed-point iterations converge for general equations.
Despite its widespread applications, the first mathematical convergence result
for linear and nonlinear problems did not appear until 2015 in Toth . For
nonlinear eigenvalue problems, there is still no strictly theoretical
analysis. In general, the mixing scheme is used in fixed-point iterations to
accelerate the convergence rate. Thus, we use the mixing scheme in the Newton
iteration scheme to design a more efficient multigrid method. By treating the
nonlinear eigenvalue problem as a special nonlinear equation, we just solve a
series of linearized boundary value problems derived from the Newton iteration
scheme defined on a multigrid space sequence. Next, we mix the approximate
solution in each iteration step to generate a more accurate approximation.
More importantly, we can rigorously prove convergence for the residuals of the
nonlinear eigenvalue problem after a one-time iteration step. This may provide
some inspiration to prove the theoretical conclusions of the self-consistent
field iteration for nonlinear eigenvalue problems.

The remainder of this paper is organized as follows. In Section 2, some
preliminaries about the nonlinear eigenvalue problems and the corresponding
finite element method are presented. In Section 3, we use the Newton iteration
to solve the nonlinear eigenvalue problems and provide several useful error
estimates. In Section 4, the multigrid method based on Newton iteration is
introduced. In Section 5, we adopt the mixing scheme to provide an improved
multigrid algorithm. Numerical experiments are presented in Section 6 to
verify the theoretical results derived in this paper. Some concluding remarks
are presented in the last section.

##  2 Preliminaries of nonlinear eigenvalue problems

In this section, we first introduce some standard notations for Sobolev space
$W^{s,p}(\Omega)$ and the associated norm $\|\cdot\|_{s,p,\Omega}$ on a
bounded domain $\Omega\in\mathbb{R}^{d}\ (d\geq 1)$. For $p=2$, we denote
$\|\cdot\|_{s,\Omega}=\|\cdot\|_{s,2,\Omega}$, $H^{s}(\Omega)=W^{s,2}(\Omega)$
and $H_{0}^{1}(\Omega)=\\{v\in H^{1}(\Omega):v|_{\partial\Omega}=0\\}$.
Hereafter, we use $\lesssim$, $\gtrsim$, $\approx$ to denote $\leq$, $\geq$,
$=$ with a mesh-independent coefficient for simplicity.

In this paper, we focus on the nonlinear eigenvalue problems arising from the
following variational problem:

$\displaystyle I=\inf\\{E(v),v\in
H_{0}^{1}(\Omega),\int_{\Omega}v^{2}d\Omega=1\\},$ (1)

where the energy functional $E$ is defined by

$\displaystyle
E(v)=\dfrac{1}{2}a(v,v)+\dfrac{1}{2}\int_{\Omega}F(v^{2})d\Omega,$

with

$\displaystyle a(u,v)=\int_{\Omega}(\mathcal{A}\nabla u\cdot\nabla
v+Vuv)d\Omega.$

To simplify the notation, let us define $f(t)=F^{\prime}(t)$. Then, for all
$v\in H_{0}^{1}(\Omega)$, $E^{\prime}(v)=A_{v}v$ holds with

$\displaystyle A_{v}=-{\rm div}(\mathcal{A}\nabla\cdot)+V+f(v^{2}).$ (2)

Next, applying the Lagrange multiplier method to (1) yields the following
nonlinear eigenvalue problem:

$\left\\{\begin{array}[]{rcl}A_{u}u&=&\lambda u,\ \ \text{ in }\Omega,\\\
u&=&0,\ \ \ \ \text{ on }\partial\Omega,\\\
\int_{\Omega}u^{2}d\Omega&=&1,\end{array}\right.$ (3)

where $\lambda$ comes from the Lagrange multiplier of the constraint
$\|u\|_{0}^{2}=1$.

For future theoretical analysis, we assume that the following conditions hold
true (see CancesChakirMaday ; ChenZhou ):

  * •

$\mathcal{A}$ is a symmetric positive-definite matrix;

  * •

$V\in L^{p}(\Omega)$ for some $p>\max(1,d/2)$;

  * •

$F\in C^{1}([0,\infty),\mathbb{R})\cap C^{2}((0,\infty),\mathbb{R})$ and
$F^{\prime\prime}>0$ on $(0,\infty)$;

  * •

$\exists\ 0\leq m<2\ \text{and}\ C\in\mathbb{R}_{+},\ \text{s.t.}\ \forall
t\geq 0,\ |F^{\prime}(t)|\leq C(1+t^{m})$;

  * •

$F^{\prime\prime}(t)t$ is locally bounded in $[0,\infty)$;

  * •

$\exists\ 1<r\leq 2$ and $0\leq s\leq 5-r$ such that $\forall a>0,\exists\
C\in\mathbb{R}_{+},\ \text{s.t.}\ \forall\ 0<t_{1}\leq a,\ \forall\
t_{2}\in\mathbb{R},$

$\displaystyle|f(t_{2}^{2})t_{2}-f(t_{1}^{2})t_{1}-2f^{\prime}(t_{1}^{2})t_{1}^{2}(t_{2}-t_{1})|\leq
C(1+|t_{2}|^{s})|t_{2}-t_{1}|^{r};$

  * •

$\exists\ 0\leq q<1\ \text{and}\ C\in\mathbb{R}_{+},\ \text{s.t.}\ \forall
t\geq 0,\ |f^{\prime}(t)|+|f^{\prime\prime}(t)t|\leq C(1+t^{q})$.

In the remainder of this section, we present two lemmas for the nonlinear
eigenvalue problem (3), and the detailed proof can be found in
CancesChakirMaday .

###### Lemma 1

There exist $M\in\mathbb{R}_{+}$ and $\beta>0$ such that for all $v\in
H_{0}^{1}(\Omega)$, there holds

$\displaystyle 0\leq((A_{u}-\lambda)v,v)\leq M\|v\|_{1}^{2},$ (4)

and

$\displaystyle\beta\|v\|_{1}^{2}\leq((E^{\prime\prime}(u)-\lambda)v,v)\leq
M\|v\|_{1}^{2},$ (5)

where $E^{\prime\prime}(u)$ has the following form

$\displaystyle(E^{\prime\prime}(u)w,v)=(A_{u}w,v)+2(f^{\prime}(u^{2})u^{2}w,v).$
(6)

Now, we introduce the finite element method for the nonlinear eigenvalue
problem (3). To use the finite element method, we define the variational form
for (3) as follows: Find $(\lambda,u)\in\mathbb{R}\times H_{0}^{1}(\Omega)$
such that $\|u\|_{0}=1$ and

$\displaystyle(\mathcal{A}\nabla u,\nabla
v)+(Vu+f(u^{2})u,v)=\lambda(u,v),\quad\forall v\in H_{0}^{1}(\Omega).$ (7)

Next, we define the finite element space $V_{h}$ on a triangulation
$\mathcal{T}_{h}$. The finite element space $V_{h}$ is composed of piecewise
polynomials such that $V_{h}\subset H_{0}^{1}(\Omega)$ and

$\lim_{h\rightarrow 0}\inf_{v_{h}\in V_{h}}\|w-v_{h}\|_{1}=0,\ \ \ \forall
w\in H_{0}^{1}(\Omega).$ (8)

Based on the finite element space $V_{h}$, we can derive the finite element
solution for (7) by solving the discrete nonlinear eigenvalue problem as
follows: Find $(\lambda_{h},u_{h})\in\mathbb{R}\times V_{h}$ such that
$\|u_{h}\|_{0}=1$ and

$(\mathcal{A}\nabla u_{h},\nabla
v_{h})+(Vu_{h}+f(u_{h}^{2})u_{h},v_{h})=\lambda_{h}(u_{h},v_{h}),\quad\forall
v_{h}\in V_{h}.$ (9)

The standard error estimates for the finite element approximate eigenpair
$(\lambda_{h},u_{h})$ are described in the following lemma.

###### Lemma 2

There exists $h_{0}>0$, such that for all $0<h<h_{0}$, the following error
estimates hold true

$\displaystyle\|u-u_{h}\|_{1}\lesssim\delta_{h}(u),$ (10)

and

$\displaystyle|\lambda-\lambda_{h}|\lesssim\|u-u_{h}\|_{1}^{2}+\|u-u_{h}\|_{L^{6/(5-2m)}},$
(11)

where

$\delta_{h}(u):=\inf_{v_{h}\in V_{h}}\|u-v_{h}\|_{1}.$ (12)

##  3 Newton iteration for nonlinear eigenvalue problems

To introduce Newton iteration for nonlinear eigenvalue problem (7), we first
give some symbols to simplify the description. Let us define

$\displaystyle\langle\mathcal{F}(\lambda,u),v\rangle=(\mathcal{A}\nabla
u,\nabla v)+(Vu+f(u^{2})u-\lambda u,v),\quad\forall v\in H_{0}^{1}(\Omega).$
(13)

The Fréchet derivation of $\mathcal{F}$ with respect to $u$ at $(\lambda,u)$
is denoted by

$\displaystyle\langle\mathcal{F}^{\prime}_{u}(\lambda,u)w,v\rangle$
$\displaystyle=$ $\displaystyle(\mathcal{A}\nabla w,\nabla
v)+((V+f(u^{2})-\lambda)w,v)+2(f^{\prime}(u^{2})u^{2}w,v)$ (14)
$\displaystyle=$ $\displaystyle(E^{\prime\prime}(u)w,v)-\lambda(w,v).$

Next, we denote the product space $\mathbb{R}\times H_{0}^{1}(\Omega)$ by $X$.
Let us define the norms in $X$ as

$\displaystyle\|(\lambda,u)\|_{X}=|\lambda|+\|u\|_{1}\quad\text{and}\quad\|(\lambda,u)\|_{0}=|\lambda|+\|u\|_{0},\quad\forall(\lambda,u)\in
X.$ (15)

Let us define

$\displaystyle\langle\mathcal{G}(\lambda,u),(\mu,v)\rangle=\langle\mathcal{F}(\lambda,u),v\rangle+\frac{1}{2}\mu(1-\int_{\Omega}u^{2}d\Omega),\quad\forall(\mu,v)\in
X.$ (16)

Obviously, the nonlinear eigenvalue problem (7) is equivalent to the following
nonlinear equation defined in the product space $X$: Find $(\lambda,u)\in X$
such that

$\displaystyle\langle\mathcal{G}(\lambda,u),(\mu,v)\rangle=0,\quad\forall(\mu,v)\in
X.$ (17)

Similarly, to solve (17) using the finite element method, we also introduce a
finite dimensional space $X_{h}:=\mathbb{R}\times V_{h}$ as the approximation
for $X$. Then based on $X_{h}$, the finite element equation (9) is equivalent
to the following nonlinear equation defined in the tensor product space
$X_{h}$: Find $(\lambda_{h},u_{h})\in X_{h}$ such that

$\displaystyle\langle\mathcal{G}(\lambda_{h},u_{h}),(\mu,v_{h})\rangle=0,\quad\forall(\mu,v_{h})\in
X_{h}.$ (18)

The Fréchet derivation of $\mathcal{G}$ at $(\lambda,u)$ is denoted by:

$\displaystyle\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle=\langle\mathcal{F}_{u}^{\prime}(\lambda,u)w,v\rangle-\gamma(u,v)-\mu(u,w),\quad\forall(\mu,v)\in
X.$ (19)

Next, we begin to introduce the Newton iteration for (18). Given an initial
value $(\lambda^{(0)},u^{(0)})\in X_{h}$, the one step of the Newton iteration
is defined as follows, which will generate a more accurate approximate
solution: Find $(\lambda^{(1)},u^{(1)})\in X_{h}$ such that

$\displaystyle\big{\langle}\mathcal{G}^{\prime}(\lambda^{(0)},u^{(0)})(\lambda^{(1)}-\lambda^{(0)},u^{(1)}-u^{(0)}),(\mu,v_{h})\big{\rangle}=-\big{\langle}\mathcal{G}(\lambda^{(0)},u^{(0)}),(\mu,v_{h})\big{\rangle},\forall(\mu,v_{h})\in
X_{h}.$ (20)

Based on the definitions (16) and (19), the linearized boundary value problem
(20) can be rewritten as follows: Find
$(\lambda^{(1)},u^{(1)})\in\mathbb{R}\times V_{h}$ such that for any
$(\mu,v_{h})\in\mathbb{R}\times V_{h}$, there holds

$\displaystyle\left\\{\begin{array}[]{lll}a(\lambda^{(0)},u^{(0)};u^{(1)},v_{h})+b(u^{(0)};v_{h},\lambda^{(1)})&=&\langle\mathcal{F}^{\prime}_{u}(\lambda^{(0)},u^{(0)})u^{(0)},v_{h}\rangle\\\
&&-\langle\mathcal{F}(\lambda^{(0)},u^{(0)}),v_{h}\rangle-\lambda^{(0)}(u^{(0)},v_{h}),\\\
b(u^{(0)};u^{(1)},\mu)&=&-\mu/2-\mu(u^{(0)},u^{(0)})/2,\end{array}\right.$
(24)

where

$a(\lambda^{(0)},u^{(0)};u^{(1)},v_{h})=\langle\mathcal{F}^{\prime}_{u}(\lambda^{(0)},u^{(0)})u^{(1)},v_{h}\rangle,\
\ \ \ b(u^{(0)};u^{(1)},\mu)=-\mu(u^{(0)},u^{(1)}).$

To guarantee the well-posedness of (20) and (24), an isomorphism property is
proved in Theorem 3.1:

###### Theorem 3.1

Let $(\lambda,u)$ be the exact solution of (17). When the mesh size $h$ is
sufficiently small, we can derive the following isomorphism properties for the
operator $\mathcal{G}^{\prime}$:

$\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\gtrsim$
$\displaystyle\|(\gamma,w)\|_{X},\quad\forall(\gamma,w)\in X,$ (25)

and

$\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X_{h}}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\gtrsim$
$\displaystyle\|(\gamma,w)\|_{X},\quad\forall(\gamma,w)\in X_{h}.$ (26)

For any $(\widehat{\lambda},\widehat{u})\in X$ such that
$\|(\widehat{\lambda}-\lambda,\widehat{u}-u)\|_{X}$ is sufficiently small, we
can also derive

$\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X_{h}}\dfrac{\langle\mathcal{G}^{\prime}(\widehat{\lambda},\widehat{u})(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\gtrsim$
$\displaystyle\|(\gamma,w)\|_{X},\quad\forall(\gamma,w)\in X_{h}.$ (27)

###### Proof

Firstly, proving (25) is equivalent to prove that for any
$(\tau,f)\in\mathbb{R}\times H^{-1}(\Omega)$, the following equation is
uniquely solvable: Find $(\gamma,w)\in\mathbb{R}\times H_{0}^{1}(\Omega)$ such
that

$\displaystyle\left\\{\begin{array}[]{ll}a(\lambda,u;w,v)+b(u;v,\gamma)=(f,v),&\forall
v\in H_{0}^{1}(\Omega),\\\
b(u;w,\mu)=\mu\tau,&\forall\mu\in\mathbb{R}.\end{array}\right.$ (30)

According to the Brezzi theory Fortin , we just need to prove the following
two conclusions:
(1) The following equation

$\displaystyle a(\lambda,u;w,v)=(f,v),\quad\forall v\in V_{0},$ (31)

is uniquely solvable in $V_{0}:=\\{v\in H_{0}^{1}(\Omega):b(u;v,\mu)=0,\
\forall\mu\in\mathbb{R}\\}$.
(2) The following inf-sup condition holds true

$\displaystyle\inf_{\mu\in\mathbb{R}}\sup_{v\in
H_{0}^{1}(\Omega)}\dfrac{b(u;v,\mu)}{\|v\|_{1}|\mu|}\geq k_{b},$ (32)

for some constant $k_{b}>0$.

Firstly, from (5), we can easily find that (31) is uniquely solvable. Next,
since $\|u\|_{0}=1$, taking $v=-\mu u$ leads to

$\displaystyle\inf_{\mu\in\mathbb{R}}\sup_{v\in
H_{0}^{1}(\Omega)}\dfrac{b(u;v,\mu)}{\|v\|_{1}|\mu|}=\frac{1}{\|u\|_{1}}=:k_{b}.$

Then we derive the inf-sup condition (32). Thus, we can derive the well-
posedness of (30) based on Brezzi theory, which further yields (25).

Next, we begin to prove the second formula (26). Based on (5), let us define
the projection $P_{h}:H_{0}^{1}(\Omega)\rightarrow V_{h}$ by

$\displaystyle a(\lambda,u;w,v-P_{h}v)=0,\quad\forall w\in V_{h},\ \forall
v\in H_{0}^{1}(\Omega).$ (33)

Using (5) and (33), we can derive

$\displaystyle\|P_{h}v\|_{1}\lesssim\|v\|_{1},\quad\forall v\in
H_{0}^{1}(\Omega),$ (34)

and

$\displaystyle\|v-P_{h}v\|_{0}\lesssim h\|v-P_{h}v\|_{1}\lesssim
h\|v\|_{1},\quad\forall v\in H_{0}^{1}(\Omega).$ (35)

Combining (25), (33), (34) and (35) leads to

$\displaystyle\|(\gamma,w)\|_{X}$ $\displaystyle\lesssim$
$\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}}$
(36) $\displaystyle=$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(0,v-P_{h}v)\rangle+\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle=$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{a(\lambda,u;w,v-P_{h}v)-\gamma(u,v-P_{h}v)+\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle=$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{-\gamma(u,v-P_{h}v)+\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\lesssim$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\gamma\|u\|_{0}\|v-P_{h}v\|_{0}+\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\lesssim$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\gamma
h\|v\|_{1}+\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\lesssim$ $\displaystyle
h\|(\gamma,w)\|_{X}+\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}.$

Thus, when the mesh size $h$ is small enough, we can obtain

$\displaystyle\|(\gamma,w)\|_{X}$ $\displaystyle\lesssim$
$\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,v)\|_{X}}$
(37) $\displaystyle\lesssim$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,P_{h}v)\rangle}{\|(\mu,P_{h}v)\|_{X}}$
$\displaystyle\lesssim$ $\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X_{h}}\dfrac{\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}},$

which is just the second formula (26).

For the last formula (27), we first prove a Lipschitz continuity for
$\mathcal{G}^{\prime}$ based on the assumptions presented in Section 2, Hölder
inequality and the imbedding theorem (see Adams ):

$\displaystyle\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle-\langle\mathcal{G}^{\prime}(\widehat{\lambda},\widehat{u})(\gamma,w),(\mu,v)\rangle$
$\displaystyle=$
$\displaystyle\big{(}(f(u^{2})-f(\widehat{u}^{2}))w,v\big{)}+((\widehat{\lambda}-\lambda)w,v)$
$\displaystyle+2\big{(}(f^{\prime}(u^{2})u^{2}-f^{\prime}(\widehat{u}^{2})\widehat{u}^{2})w,v\big{)}-\gamma(\widehat{u}-u,v)-\mu(\widehat{u}-u,w)$
$\displaystyle=$ $\displaystyle
2\big{(}f^{\prime}(\xi^{2})\xi(u-\widehat{u})w,v\big{)}+((\widehat{\lambda}-\lambda)w,v)$
$\displaystyle+4\big{(}(f^{\prime\prime}(\xi^{2})\xi^{3}+f^{\prime}(\xi^{2})\xi)(u-\widehat{u})w,v\big{)}-\gamma(\widehat{u}-u,v)-\mu(\widehat{u}-u,w)$
$\displaystyle\lesssim$
$\displaystyle\big{(}(1+\xi^{2q})\xi(u-\widehat{u})w,v\big{)}+((\widehat{\lambda}-\lambda)w,v)$
$\displaystyle+\big{(}(1+\xi^{2q})\xi(u-\widehat{u})w,v\big{)}-\gamma(\widehat{u}-u,v)-\mu(\widehat{u}-u,w)$
$\displaystyle\lesssim$
$\displaystyle\|\xi\|_{0,4}\|u-\widehat{u}\|_{0,4}\|w\|_{0,4}\|v\|_{0,4}+\|\xi^{2q}\|_{0,3/q}\|\xi\|_{0,6}\|u-\widehat{u}\|_{0,6}\|w\|_{0,6}\|v\|_{0,6/(3-2q)}$
$\displaystyle+|\widehat{\lambda}-\lambda|\|w\|_{0}\|v\|_{0}+|\gamma|\|\widehat{u}-u\|_{0}\|v\|_{0}+|\mu|\|\widehat{u}-u\|_{0}\|w\|_{0}$
$\displaystyle\lesssim$
$\displaystyle\|\xi\|_{1}\|u-\widehat{u}\|_{1}\|w\|_{1}\|v\|_{1}+\|\xi\|_{1}^{2q}\|\xi\|_{1}\|u-\widehat{u}\|_{1}\|w\|_{1}\|v\|_{1}$
$\displaystyle+|\widehat{\lambda}-\lambda|\|w\|_{0}\|v\|_{0}+|\gamma|\|\widehat{u}-u\|_{0}\|v\|_{0}+|\mu|\|\widehat{u}-u\|_{0}\|w\|_{0},$

where $\xi=\tau u+(1-\tau)\widehat{u}$ with $\tau\in(0,1)$.

Thus, we can derive

$\displaystyle\langle\mathcal{G}^{\prime}(\lambda,u)(\gamma,w),(\mu,v)\rangle-\langle\mathcal{G}^{\prime}(\widehat{\lambda},\widehat{u})(\gamma,w),(\mu,v)\rangle$
(38) $\displaystyle\lesssim$
$\displaystyle\|(\lambda-\widehat{\lambda},u-\widehat{u})\|_{X}\|(\gamma,w)\|_{X}\|(\mu,v)\|_{X}.$

Let $\varepsilon$ be a sufficiently small constant and
$\|(\lambda-\widehat{\lambda},u-\widehat{u})\|_{X}\leq\varepsilon$. Based on
(26) and (38), for any $(\gamma,w)\in X_{h}$, we have

$\displaystyle\|(\gamma,w)\|_{X}$ $\displaystyle\lesssim$
$\displaystyle\mathop{\sup}\limits_{(\mu,v)\in
X_{h}}\dfrac{\varepsilon\|(\gamma,w)\|_{X}\|(\mu,v)\|_{X}+\langle\mathcal{G}^{\prime}(\widehat{\lambda},\widehat{u})(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}}$
$\displaystyle\lesssim$
$\displaystyle\varepsilon\|(\gamma,w)\|_{X}+\mathop{\sup}\limits_{(\mu,v)\in
X_{h}}\dfrac{\langle\mathcal{G}^{\prime}(\widehat{\lambda},\widehat{u})(\gamma,w),(\mu,v)\rangle}{\|(\mu,v)\|_{X}}.$

Then $(\ref{coercive of g3})$ can be derived using the same way as the proof
process in (37). This completes the proof.

Next, in Theorem 3.2, we provide the estimate for the residual of the Newton
iteration scheme (20).

###### Theorem 3.2

For any $(\lambda^{(0)},u^{(0)})\in X,(\lambda^{(1)},u^{(1)})\in X$, we can
derive the following expansion

$\displaystyle\langle\mathcal{G}(\lambda^{(1)},u^{(1)}),(\sigma,\eta)\rangle$
$\displaystyle=$
$\displaystyle\langle\mathcal{G}(\lambda^{(0)},u^{(0)}),(\sigma,\eta)\rangle+\langle\mathcal{G}^{\prime}(\lambda^{(0)},u^{(0)})(\lambda^{(1)}-\lambda^{(0)},u^{(1)}-u^{(0)}),(\sigma,\eta)\rangle$
(39)
$\displaystyle+~{}R\big{(}(\lambda^{(0)},u^{(0)}),(\lambda^{(1)},u^{(1)}),(\sigma,\eta)\big{)},\quad\quad\quad\quad\quad\quad\forall(\sigma,\eta)\in
X,$

where the residual
$R\big{(}(\lambda^{(0)},u^{(0)}),(\lambda^{(1)},u^{(1)}),(\sigma,\eta)\big{)}$
satisfies the following estimate

$\displaystyle\big{|}R\big{(}(\lambda^{(0)},u^{(0)}),(\lambda^{(1)},u^{(1)}),(\sigma,\eta)\big{)}\big{|}\lesssim\|(\lambda^{(0)}-\lambda^{(1)},u^{(0)}-u^{(1)})\|_{X}^{2}\|(\sigma,\eta)\|_{X}.$

###### Proof

We prove this theorem through the Taylor expansion. Let us define

$\displaystyle\eta(t)=\big{\langle}\mathcal{G}\big{(}(\lambda^{(0)},u^{(0)})+t\\{(\lambda^{(1)},u^{(1)})-(\lambda^{(0)},u^{(0)})\\}\big{)},(\sigma,\eta)\big{\rangle},\
\ \ t\in[0,1].$ (40)

Obviously,

$\displaystyle\eta(0)=\langle\mathcal{G}(\lambda^{(0)},u^{(0)}),(\sigma,\eta)\rangle\
\ \ \text{and}\ \ \
\eta(1)=\langle\mathcal{G}(\lambda^{(1)},u^{(1)}),(\sigma,\eta)\rangle.$ (41)

The first order derivative and the second order derivative of $\eta(t)$ with
respect to $t$ are listed as follows:

$\displaystyle\eta^{\prime}(t)$ $\displaystyle=$
$\displaystyle\big{\langle}\mathcal{G}^{\prime}((\lambda^{(0)},u^{(0)})+t\\{(\lambda^{(1)},u^{(1)})-(\lambda^{(0)},u^{(0)})\\})((\lambda^{(1)},u^{(1)})-(\lambda^{(0)},u^{(0)})),(\sigma,\eta)\big{\rangle}$
(42) $\displaystyle=$
$\displaystyle\big{\langle}\mathcal{F}_{u}^{\prime}\big{(}(\lambda^{(0)},u^{(0)})+t\\{(\lambda^{(1)},u^{(1)})-(\lambda^{(0)},u^{(0)})\\}\big{)}(u^{(1)}-u^{(0)}),\eta\big{\rangle}$
$\displaystyle-(\lambda^{(1)}-\lambda^{(0)})(u^{(0)}+t(u^{(1)}-u^{(0)}),\eta)-\sigma(u^{(0)}+t(u^{(1)}-u^{(0)}),u^{(1)}-u^{(0)})$
$\displaystyle=$
$\displaystyle(\mathcal{A}\nabla(u^{(1)}-u^{(0)}),\nabla\eta)+\big{(}(u^{(0)}+f((u^{(0)}+t(u^{(1)}-u^{(0)}))^{2}))(u^{(1)}-u^{(0)}),\eta\big{)}$
$\displaystyle-\big{(}(\lambda^{(0)}+t(\lambda^{(1)}-\lambda^{(0)}))(u^{(1)}-u^{(0)}),\eta\big{)}$
$\displaystyle+2\big{(}f^{\prime}((u^{(0)}+t(u^{(1)}-u^{(0)}))^{2})(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2}(u^{(1)}-u^{(0)}),\eta\big{)}$
$\displaystyle-(\lambda^{(1)}-\lambda^{(0)})(u^{(0)}+t(u^{(1)}-u^{(0)}),\eta)-\sigma(u^{(0)}+t(u^{(1)}-u^{(0)}),u^{(1)}-u^{(0)}),$

and

$\displaystyle\eta^{\prime\prime}(t)$ $\displaystyle=$
$\displaystyle-2(\lambda^{(1)}-\lambda^{(0)})(u^{(1)}-u^{(0)},\eta)-\sigma(u^{(1)}-u^{(0)},u^{(1)}-u^{(0)})$
(43)
$\displaystyle+6\big{(}f^{\prime}\big{(}(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2}\big{)}\big{(}u^{(0)}+t(u^{(1)}-u^{(0)})\big{)}(u^{(1)}-u^{(0)})^{2},\eta\big{)}$
$\displaystyle+4\big{\langle}f^{\prime\prime}\big{(}(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2}\big{)}\big{(}u^{(0)}+t(u^{(1)}-u^{(0)})\big{)}^{3}(u^{(1)}-u^{(0)})^{2},\eta\big{\rangle}.$

For the nonlinear terms involved in (43), we have the following estimates

$\displaystyle\big{|}\big{(}f^{\prime}\big{(}(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2}\big{)}\big{(}u^{(0)}+t(u^{(1)}-u^{(0)})\big{)}(u^{(1)}-u^{(0)})^{2},\eta\big{)}\big{|}$
(44) $\displaystyle\lesssim$
$\displaystyle\big{|}\big{(}\big{(}1+(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2q}\big{)}\big{(}u^{(0)}+t(u^{(1)}-u^{(0)})\big{)}(u^{(1)}-u^{(0)})^{2},\eta\big{)}\big{|}$
$\displaystyle\lesssim$
$\displaystyle\big{|}\big{(}(u^{(0)}+t(u^{(1)}-u^{(0)}))(u^{(1)}-u^{(0)})^{2},\eta\big{)}$
$\displaystyle+\big{(}(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2q}(u^{(0)}+t(u^{(1)}-u^{(0)}))(u^{(1)}-u^{(0)})^{2},\eta\big{)}\big{|}$
$\displaystyle\lesssim$
$\displaystyle\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{0,4}\|u^{(1)}-u^{(0)}\|_{0,4}^{2}\|\eta\|_{0,4}$
$\displaystyle+\|(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2q}\|_{0,3/q}\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{0,6}\|u^{(1)}-u^{(0)}\|_{0,6}^{2}\|\eta\|_{0,6/(3-2q)}$
$\displaystyle\lesssim$
$\displaystyle\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{0,4}\|u^{(1)}-u^{(0)}\|_{0,4}^{2}\|\eta\|_{0,4}$
$\displaystyle+\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{0,6}^{2q}\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{0,6}\|u^{(1)}-u^{(0)}\|_{0,6}^{2}\|\eta\|_{0,6/(3-2q)}$
$\displaystyle\lesssim$
$\displaystyle\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{1}\|u^{(1)}-u^{(0)}\|_{1}^{2}\|\eta\|_{1}$
$\displaystyle+\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{1}^{2q}\|u^{(0)}+t(u^{(1)}-u^{(0)})\|_{1}\|u^{(1)}-u^{(0)}\|_{1}^{2}\|\eta\|_{1}$
$\displaystyle\lesssim$
$\displaystyle\|u^{(1)}-u^{(0)}\|_{1}^{2}\|\eta\|_{1},$

and

$\displaystyle\big{|}\big{\langle}f^{\prime\prime}\big{(}(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2}\big{)}\big{(}u^{(0)}+t(u^{(1)}-u^{(0)})\big{)}^{3}(u^{(1)}-u^{(0)})^{2},\eta\big{\rangle}\big{|}$
(45) $\displaystyle\lesssim$
$\displaystyle\big{|}\big{(}(1+(u^{(0)}+t(u^{(1)}-u^{(0)}))^{2q}\big{)}\big{(}u^{(0)}+t(u^{(1)}-u^{(0)})\big{)}(u^{(1)}-u^{(0)})^{2},\eta)\big{|}$
$\displaystyle\lesssim$
$\displaystyle\|u^{(1)}-u^{(0)}\|_{1}^{2}\|\eta\|_{1},$

where the Hölder inequality and the imbedding theorem are used.

Based on (43), (44) and (45), we can derive

$\displaystyle|\eta^{\prime\prime}(t)|$ $\displaystyle\lesssim$
$\displaystyle|\lambda^{(1)}-\lambda^{(0)}|\|u^{(1)}-u^{(0)}\|_{0}\|\eta\|_{0}+|\sigma|\|u^{(1)}-u^{(0)}\|_{0}^{2}+\|u^{(1)}-u^{(0)}\|_{1}^{2}\|\eta\|_{1}$
(46) $\displaystyle\lesssim$
$\displaystyle\|(\lambda^{(0)}-\lambda^{(1)},u^{(0)}-u^{(1)})\|_{X}^{2}\|(\sigma,\eta)\|_{X}.$

Thus, using (46) and the following Taylor expansion:

$\eta(1)=\eta(0)+\eta^{\prime}(0)+\int_{0}^{1}\eta^{\prime\prime}(t)(1-t)dt,$

we can derive

$\displaystyle\langle\mathcal{G}(\lambda^{(1)},u^{(1)}),(\sigma,\eta)\rangle$
$\displaystyle=$
$\displaystyle\langle\mathcal{G}(\lambda^{(0)},u^{(0)}),(\sigma,\eta)\rangle+\langle\mathcal{G}^{\prime}(\lambda^{(0)},u^{(0)})((\lambda^{(1)},u^{(1)})-(\lambda^{(0)},u^{(0)})),(\sigma,\eta)\rangle$
$\displaystyle+R((\lambda^{(0)},u^{(0)}),(\lambda^{(1)},u^{(1)}),(\sigma,\eta))$

with the residual

$\displaystyle|R((\lambda^{(0)},u^{(0)}),(\lambda^{(1)},u^{(1)}),(\sigma,\eta))|\lesssim\|(\lambda^{(0)}-\lambda^{(1)},u^{(0)}-u^{(1)})\|_{X}^{2}\|(\sigma,\eta)\|_{X}.$

Then the proof is completed.

##  4 Multigrid method for nonlinear eigenvalue problems based on Newton
iteration

In this section, a novel multigrid method is proposed for solving the
nonlinear eigenvalue problems based on the Newton iteration. To design the
multigrid method, let us construct a nested  sequence of finite element spaces
based on a nested multigrid sequence in the first step. The sequence of meshes
is denoted by $\\{\mathcal{T}_{h_{k}}\\}_{k=1}^{n}$. For any two consecutive
meshes $\mathcal{T}_{h_{k}}$ and $\mathcal{T}_{h_{k-1}}(k\geq 2)$,
$\mathcal{T}_{h_{k}}$ is produced through a one-time uniform refinement form
$\mathcal{T}_{h_{k-1}}$. The corresponding mesh sizes satisfy the following
conditions:

$\displaystyle h_{k}\approx\dfrac{1}{\beta}h_{k-1},\quad k=2,\cdots,n,$ (47)

where the refinement index $\beta>1$. Meanwhile, the following approximate
relationship holds true:

$\displaystyle\delta_{h_{k}}(u)\approx\dfrac{1}{\beta}\delta_{h_{k-1}}(u),\quad
k=2,\cdots,n.$ (48)

Based on the multigrid sequence $\\{\mathcal{T}_{h_{k}}\\}_{k=1}^{n}$, we
denote the corresponding  sequence of finite element spaces by:

$\displaystyle V_{h_{1}}\subset V_{h_{2}}\subset\cdots\subset V_{h_{n}}\subset
H_{0}^{1}(\Omega),$ (49)

and denote the sequence of tensor product finite element spaces by:

$\displaystyle X_{h_{1}}\subset X_{h_{2}}\subset\cdots\subset X_{h_{n}}\subset
X.$ (50)

###  4.1 One step of the multigrid method

To describe the multigrid method more clearly, we first introduce how to
perform one step of the Newton iteration.

For two consecutive finite element spaces $X_{h_{k}}$ and $X_{h_{k+1}}$,
assume that we have obtained an approximate solution
$(\lambda^{h_{k}},u^{h_{k}})\in X_{h_{k}}$, Algorithm 4.1 shows the procedure
for obtaining a new approximate solution $(\lambda^{h_{k+1}},u^{h_{k+1}})\in
X_{h_{k+1}}$.

###### Algorithm 4.1

One step of the multigrid method

  1. 1.

Assume that we have obtained an approximate solution
$(\lambda^{h_{k}},u^{h_{k}})\in X_{h_{k}}$.

  2. 2.

Solve the linearized equation derived from Newton iteration in the finite
element space $X_{h_{k+1}}$: Find $(\lambda^{h_{k+1}},u^{h_{k+1}})\in
X_{h_{k+1}}$ such that for any $(\mu,v_{h_{k+1}})\in X_{h_{k+1}}$, there holds

$\displaystyle\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle$
(51) $\displaystyle=$
$\displaystyle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle+\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle.$

The equation (51) can be solved by the mixed finite element method in the
following form: Find $(\lambda^{h_{k+1}},u^{h_{k+1}})\in\mathbb{R}\times
V_{h_{k+1}}$ such that for any $(\mu,v_{h_{k+1}})\in\mathbb{R}\times
V_{h_{k+1}}$, there holds

$\displaystyle\left\\{\begin{array}[]{lll}a(\lambda^{h_{k}},u^{h_{k}};u^{h_{k+1}},v_{h_{k+1}})+b(u^{h_{k}};v_{h_{k+1}},\lambda^{h_{k+1}})&=&\langle\mathcal{F}_{u}^{\prime}(\lambda^{h_{k}},u^{h_{k}})u^{h_{k}},v_{h_{k+1}}\rangle\\\
&&-\langle\mathcal{F}(\lambda^{h_{k}},u^{h_{k}}),v_{h_{k+1}}\rangle-\lambda^{h_{k}}(u^{h_{k}},v_{h_{k+1}}),\\\
b(u^{h_{k}};u^{h_{k+1}},\mu)&=&-\mu/2-\mu(u^{h_{k}},u^{h_{k}})/2.\end{array}\right.$

In order to simplify the notation, we use the following symbol to denote the
above two solving steps:

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=Newton_{-}Iteration(\lambda^{h_{k}},u^{h_{k}},X_{h_{k+1}}).$

Next, we can prove that the new approximate solution
$(\lambda^{h_{k+1}},u^{h_{k+1}})\in X_{h_{k+1}}$ derived by Algorithm 4.1 has
a better accuracy than $(\lambda^{h_{k}},u^{h_{k}})\in X_{h_{k}}$. The
detailed conclusion is presented in Theorem 4.1.

###### Theorem 4.1

After implementing Algorithm 4.1, the new approximate solution
$(\lambda^{h_{k+1}},u^{h_{k+1}})$ has the following error estimate

$\displaystyle\|(\lambda-\lambda^{h_{k+1}},u-u^{h_{k+1}})\|_{X}\lesssim\delta_{h_{k+1}}(u)+\|(\lambda-\lambda^{h_{k}},u-u^{h_{k}})\|_{X}^{2}.$

###### Proof

We use $(\lambda_{h_{k+1}},u_{h_{k+1}})$ to denote the standard finite element
solution of (9), that is $(\lambda_{h_{k+1}},u_{h_{k+1}})$ satisfies

$\displaystyle(\mathcal{A}\nabla u_{h_{k+1}},\nabla
v_{h_{k+1}})+(Vu_{h_{k+1}},v_{h_{k+1}})+(f(u_{h_{k+1}}^{2})u_{h_{k+1}},v_{h_{k+1}})$
$\displaystyle=$
$\displaystyle\lambda_{h_{k+1}}(u_{h_{k+1}},v_{h_{k+1}}),\quad\forall
v_{h_{k+1}}\in V_{h_{k+1}},$

or equivalently

$\displaystyle\langle\mathcal{G}(\lambda_{h_{k+1}},u_{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle=0,\quad\forall(\mu,v_{h_{k+1}})\in
X_{h_{k+1}}.$ (53)

Combining (51) and (53) leads to

$\displaystyle\big{\langle}\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda_{h_{k+1}}-\lambda^{h_{k+1}},u_{h_{k+1}}-u^{h_{k+1}}),(\mu,v_{h_{k+1}})\big{\rangle}$
(54) $\displaystyle=$
$\displaystyle\big{\langle}\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\big{\rangle}+\big{\langle}\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda_{h_{k+1}}-\lambda^{h_{k}},u_{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\big{\rangle}$
$\displaystyle=$
$\displaystyle-\big{\langle}\mathcal{G}(\lambda_{h_{k+1}},u_{h_{k+1}}),(\mu,v_{h_{k+1}})\big{\rangle}+\big{\langle}\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\big{\rangle}$
$\displaystyle+\big{\langle}\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda_{h_{k+1}}-\lambda^{h_{k}},u_{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\big{\rangle}$
$\displaystyle=$
$\displaystyle-R((\lambda^{h_{k}},u^{h_{k}}),(\lambda_{h_{k+1}},u_{h_{k+1}}),(\mu,v_{h_{k+1}})).$

Combining (27), Theorem 3.2 and (54), we can derive

$\displaystyle\|(\lambda_{h_{k+1}}-\lambda^{h_{k+1}},u_{h_{k+1}}-u^{h_{k+1}})\|_{X}\lesssim\|(\lambda_{h_{k+1}}-\lambda^{h_{k}},u_{h_{k+1}}-u^{h_{k}})\|_{X}^{2}.$
(55)

Using the standard error estimates in Lemma 2 and (55), we can obtain

$\displaystyle\|(\lambda-\lambda^{h_{k+1}},u-u^{h_{k+1}})\|_{X}$
$\displaystyle\lesssim$
$\displaystyle\|(\lambda-\lambda_{h_{k+1}},u-u_{h_{k+1}})\|_{X}+\|(\lambda_{h_{k+1}}-\lambda^{h_{k+1}},u_{h_{k+1}}-u^{h_{k+1}})\|_{X}$
$\displaystyle\lesssim$
$\displaystyle\delta_{h_{k+1}}(u)+\|(\lambda-\lambda^{h_{k}},u-u^{h_{k}})\|_{X}^{2}.$

Then we complete the proof.

###  4.2 Multigrid method based on Newton iteration

In this subsection, we introduce the multigrid method to solve the nonlinear
eigenvalue problems. We use the following nested  sequence of finite element
spaces:

$\displaystyle V_{h_{1}}\subset V_{h_{2}}\subset\cdots\subset V_{h_{n}}\subset
H_{0}^{1}(\Omega),$ (56)

and

$\displaystyle X_{h_{1}}\subset X_{h_{2}}\subset\cdots\subset X_{h_{n}}\subset
X.$ (57)

The multigrid method is presented in Algorithm 4.2. The new algorithm requires
solving the nonlinear eigenvalue problem directly in the initial space to
obtain an initial value; then performing Algorithm 4.1 in the subsequent
spaces. The approximate solution derived from the last space is used as the
initial iterative value in the current space.

###### Algorithm 4.2

Multigrid method based on Newton iteration

  1. 1.

Solve the following nonlinear eigenvalue problem directly in the initial space
$V_{h_{1}}$: Find $(\lambda^{h_{1}},u^{h_{1}})\in\mathbb{R}\times V_{h_{1}}$
such that $\|u^{h_{1}}\|_{0}=1$ and

$(\mathcal{A}\nabla u^{h_{1}},\nabla
v_{h_{1}})+(Vu^{h_{1}},v_{h_{1}})+(f((u^{h_{1}})^{2})u^{h_{1}},v_{h_{1}})=\lambda^{h_{1}}(u^{h_{1}},v_{h_{1}}),\quad\forall
v_{h_{1}}\in V_{h_{1}}.$

  2. 2.

For $k=1,\cdots,n-1$, obtain a new approximate solution
$(\lambda^{h_{k+1}},u^{h_{k+1}})\in X_{h_{k+1}}$ through

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=Newton_{-}Iteration(\lambda^{h_{k}},u^{h_{k}},X_{h_{k+1}}),$

End for.

  3. 3.

Finally, we obtain the approximate solution $(\lambda^{h_{n}},u^{h_{n}})\in
X_{h_{n}}$.

The error estimate for the final approximate solution
$(\lambda^{h_{n}},u^{h_{n}})\in X_{h_{n}}$ derived by Algorithm 4.2 is
presented in Theorem 4.2.

###### Theorem 4.2

After performing Algorithm 4.2, the final approximate solution
$(\lambda^{h_{n}},u^{h_{n}})\in X_{h_{n}}$ satisfies the following error
estimate

$\displaystyle\|(\lambda-\lambda^{h_{n}},u-u^{h_{n}})\|_{X}\lesssim\delta_{h_{n}}(u).$
(58)

###### Proof

We prove this theorem by the mathematical induction method. Since we solve the
nonlinear eigenvalue problem directly in the initial space, using Lemma 2
leads to

$\displaystyle\|(\lambda-{\lambda}^{h_{1}},u-{u}^{h_{1}})\|_{X}\lesssim\delta_{h_{1}}(u).$

Assume that (58) holds true for $(\lambda^{h_{n-1}},{u}^{h_{n-1}})\in
X_{h_{n-1}}$. According to Theorem 4.1, we can derive

$\displaystyle\|(\lambda-{\lambda}^{h_{n}},u-{u}^{h_{n}})\|_{X}$
$\displaystyle\lesssim$
$\displaystyle\delta_{h_{n}}(u)+\|(\lambda-\lambda^{h_{n-1}},u-u^{h_{n-1}})\|_{X}^{2}$
$\displaystyle\lesssim$
$\displaystyle\delta_{h_{n}}(u)+\delta_{h_{n-1}}^{2}(u)$
$\displaystyle\lesssim$ $\displaystyle\delta_{h_{n}}(u).$

So (58) also holds true for $(\lambda^{h_{n}},{u}^{h_{n}})\in X_{h_{n}}$. Then
the proof is completed.

###  4.3 Work estimate of Algorithm 4.2

In this subsection, we end Section 4 by estimating the computational work of
Algorithm 4.2 briefly. Let us use $N_{k}$ to denote the dimensions of
$V_{h_{k}}$. The dimensions of the  sequence of finite element spaces satisfy

$N_{k}\approx\beta^{d(k-n)}N_{n},\ \ \ k=1,2,\cdots,n.$

Then the computational work of Algorithm 4.2 is presented in Theorem 4.3,
which shows that Algorithm 4.2 has linear computational complexity.

###### Theorem 4.3

Assume that solving the nonlinear eigenvalue problem directly in the initial
space $V_{h_{1}}$ needs work $\mathcal{O}(M_{h_{1}})$, and solving the
linearized boundary value problem (51) in $V_{h_{k}}$ needs work
$\mathcal{O}(N_{k})$ for $k=2,3,\cdots,n$. Then the total computational work
of Algorithm 4.2 is $\mathcal{O}(N_{n}+M_{h_{1}})$. Furthermore, the linear
computational complexity $\mathcal{O}(N_{n})$ can be derived provided
$M_{h_{1}}\leq N_{n}$.

###### Proof

Let us use $W_{k}$ to denote the computational work of $V_{h_{k}}$. Then the
computational work $W$ of Algorithm 4.2 can be estimated by:

$\displaystyle W$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{n}W_{k}=\mathcal{O}(M_{h_{1}}+\sum_{k=2}^{n}N_{k})$
$\displaystyle=$
$\displaystyle\mathcal{O}(\sum_{k=2}^{n}\beta^{d(k-n)}N_{n}+M_{h_{1}})$
$\displaystyle=$ $\displaystyle\mathcal{O}(N_{n}+M_{h_{1}}).$

Then we get the desired conclusion $\mathcal{O}(N_{n}+M_{h_{1}})$, and when
$M_{h_{1}}\leq N_{n}$, the total computational work changes to
$\mathcal{O}(N_{n})$.

###### Remark 1

The linearized boundary value problem (51) can be solved efficiently by the
multigrid method with linear computational complexity $\mathcal{O}(N_{k})$
(see e.g., Saddle-book ; Shaidurov ). Because the dimension of the initial
space $V_{h_{1}}$ is small and independent of the final finite element space,
it is easy to get $M_{h_{1}}\leq N_{n}$. Thus, the linear computational
complexity for Algorithm 4.2 can be obtained.

##  5 Multigrid method for nonlinear eigenvalue problems based on mixing
scheme

In this section, an improved multigrid method is designed. The motivation is
that when solving the nonlinear equation (18) by Newton iteration, we may
encounter the non-convergence issue for some complicated models. This issue is
the same as that of the classical self-consistent field iteration for the
nonlinear eigenvalue problems. In order to overcome this difficulty, the
mixing theme (damping method) is introduced to solve the nonlinear eigenvalue
problems. The most widely used mixing scheme is Anderson acceleration
andersion , which is used to improve the convergence rate for fixed-point
iterations. Although Anderson acceleration was widely used for solving various
models, the first mathematical convergence result for linear and nonlinear
problems did not appear until 2015 in Toth . For nonlinear eigenvalue
problems, there is still no strictly theoretical analysis. In general, the
mixing scheme is used in fixed-point iterations to accelerate the convergence
rate. In this section, we use the mixing scheme in Newton iteration to design
an improved multigrid method for nonlinear eigenvalue problems. Most
importantly, we can prove that the norm of the residual decreases in each step
of the Newton iteration. This may provide some inspiration to prove the
theoretical conclusions for the more general mixing scheme for nonlinear
eigenvalue problems.

Assume that we have obtained an approximate solution
$(\lambda^{h_{k}},u^{h_{k}})\in X_{h_{k}}$, we introduce a novel iteration
step in Algorithm 5.1 to obtain a better approximation
$(\lambda^{h_{k+1}},u^{h_{k+1}})\in X_{h_{k+1}}$ on the basis of Newton
iteration and the mixing scheme.

###### Algorithm 5.1

One step of the mixing scheme

  1. 1.

Given a parameter $\theta_{k+1}\in(0,1)$ and an approximate solution
$(\lambda^{h_{k}},u^{h_{k}})\in X_{h_{k}}$.

  2. 2.

Solve the linearized equation derived from Newton iteration in the finite
element space $X_{h_{k+1}}$: Find
$(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}})\in X_{h_{k+1}}$ such that
for any $(\mu,v_{h_{k+1}})\in X_{h_{k+1}}$, there holds

$\displaystyle\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle$
(59) $\displaystyle=$
$\displaystyle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle+\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle.$

The equation (59) can be solved by the mixed finite element method in the
following form: Find
$(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}})\in\mathbb{R}\times
V_{h_{k+1}}$ such that for any $(\mu,v_{h_{k+1}})\in\mathbb{R}\times
V_{h_{k+1}}$, there holds

$\displaystyle\left\\{\begin{array}[]{lll}a(\lambda^{h_{k}},u^{h_{k}};\widehat{u}^{h_{k+1}},v_{h_{k+1}})+b(u^{h_{k}};v_{h_{k+1}},\widehat{\lambda}^{h_{k+1}})&=&\langle\mathcal{F}_{u}^{\prime}(\lambda^{h_{k}},u^{h_{k}})u^{h_{k}},v_{h_{k+1}}\rangle\\\
&&-\langle\mathcal{F}(\lambda^{h_{k}},u^{h_{k}}),v_{h_{k+1}}\rangle-\lambda^{h_{k}}(u^{h_{k}},v_{h_{k+1}}),\\\
b(u^{h_{k}};\widehat{u}^{h_{k+1}},\mu)&=&-\mu/2-\mu(u^{h_{k}},u^{h_{k}})/2.\end{array}\right.$

  3. 3.

Set

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=(1-\theta_{k+1})(\lambda^{h_{k}},u^{h_{k}})+\theta_{k+1}(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}}).$

Next, we can prove that the value of the residual
$|\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle|$
monotonically decreases after performing Algorithm 5.1.

###### Theorem 5.1

Assume that we have obtained an approximate solution
$(\lambda^{h_{k}},u^{h_{k}})$, then, there exists a sufficiently small
$\theta_{k+1}$ such that the new approximation
$(\lambda^{h_{k+1}},u^{h_{k+1}})$ derived by Algorithm 5.1 has the following
convergence

$\displaystyle|\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle|\leq(1-\frac{\theta_{k+1}}{2})|\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle|,\
\ \forall(\mu,v_{h_{k+1}})\in X_{h_{k+1}}.$

###### Proof

From the following mixing scheme

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=(1-\theta_{k+1})(\lambda^{h_{k}},u^{h_{k}})+\theta_{k+1}(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}}),$

used in Algorithm 5.1, we can derive

$\displaystyle(\lambda^{h_{k+1}}-\lambda^{h_{k}},u^{h_{k+1}}-u^{h_{k}})=\theta_{k+1}(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}).$
(61)

Combining (61) and the mean value theorem leads to

$\displaystyle\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
(62)
$\displaystyle-\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k+1}}-\lambda^{h_{k}},u^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle=$
$\displaystyle\langle\mathcal{G}^{\prime}(\widetilde{\lambda}^{h_{k}},\widetilde{u}^{h_{k}})(\lambda^{h_{k+1}}-\lambda^{h_{k}},u^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle-\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k+1}}-\lambda^{h_{k}},u^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle=$
$\displaystyle\theta_{k+1}\langle\mathcal{G}^{\prime}(\widetilde{\lambda}^{h_{k}},\widetilde{u}^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle-\theta_{k+1}\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle,$

where
$(\widetilde{\lambda}^{h_{k}},\widetilde{u}^{h_{k}})=(\alpha\lambda^{h_{k+1}}+(1-\alpha)\lambda^{h_{k}},\alpha
u^{h_{k+1}}+(1-\alpha)u^{h_{k}})$ with $\alpha\in(0,1)$.

From (61), we know that when $\theta_{k+1}\rightarrow 0$, there holds

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})\rightarrow(\lambda^{h_{k}},u^{h_{k}}),$

and thus

$\displaystyle(\widetilde{\lambda}^{h_{k}},\widetilde{u}^{h_{k}})\rightarrow(\lambda^{h_{k}},u^{h_{k}}),$

which also indicates

$\displaystyle\big{\langle}\mathcal{G}^{\prime}(\widetilde{\lambda}^{h_{k}},\widetilde{u}^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\big{\rangle}$
$\displaystyle\rightarrow$
$\displaystyle\big{\langle}\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\big{\rangle}.$

Hence, we can choose a sufficiently small $\theta_{k+1}$, such that

$\displaystyle\big{|}\langle\mathcal{G}^{\prime}(\widetilde{\lambda}^{h_{k}},\widetilde{u}^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
(63)
$\displaystyle-\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle\leq$
$\displaystyle\frac{1}{2}\big{|}\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{|}\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}.$

Based on (62) and (63), we can derive

$\displaystyle\big{|}\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
(64)
$\displaystyle-\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k+1}}-\lambda^{h_{k}},u^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle\leq$
$\displaystyle\frac{\theta_{k+1}}{2}|\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}.$

Besides, using (59) and (61), we can also derive

$\displaystyle\big{|}\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
(65)
$\displaystyle-\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k+1}}-\lambda^{h_{k}},u^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle=$
$\displaystyle\big{|}\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle-\theta_{k+1}\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}}-\lambda^{h_{k}},\widehat{u}^{h_{k+1}}-u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle=$
$\displaystyle\big{|}\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle+\theta_{k+1}\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle=$
$\displaystyle\big{|}\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-(1-\theta_{k+1})\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}.$

Combining (64) and (65) leads to

$\displaystyle\big{|}\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle-(1-\theta_{k+1})\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle\big{|}$
$\displaystyle\leq$
$\displaystyle\frac{\theta_{k+1}}{2}|\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle|.$

That is

$\displaystyle|\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle|$
$\displaystyle\leq$
$\displaystyle(1-\frac{\theta_{k+1}}{2})|\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle|.$

Then we complete the proof.

In order to give a practical algorithm, there are two remaining problems. The
first one is to choose an appropriate parameter $\theta_{k+1}$. The second one
is to compute the residual after each iteration step.

For the first problem, we provide an adaptive strategy for choosing
$\theta_{k+1}$ in Algorithm 5.2, which gradually reduces the value of
$\theta_{k+1}$. For the second problem, in order to compare the following two
values
$|\langle\mathcal{G}(\lambda^{h_{k+1}},u^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle|$
and $|\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle|$
for any $(\mu,v_{h_{k+1}})\in X_{h_{k+1}}$, we define the residual for the
nonlinear eigenvalue problem (3) directly in the following way:

$\displaystyle
Resi(\lambda^{h_{k}},u^{h_{k}})=\|A_{u^{h_{k}}}u^{h_{k}}-\lambda^{h_{k}}u^{h_{k}}\|_{1}+\frac{1}{2}|1-\|u^{h_{k}}\|_{0}^{2}|.$

Then next, we use $Resi(\lambda^{h_{k}},u^{h_{k}})$ to judge the error
reduction.

The modified mixing scheme and the corresponding multigrid method are
presented in Algorithms 5.2 and 5.3, respectively.

###### Algorithm 5.2

One step of the modified mixing scheme

  1. 1.

Assume that we have obtained an approximate solution
$(\lambda^{h_{k}},u^{h_{k}})\in X_{h_{k}}$.

  2. 2.

Set $\theta_{k+1}=1$.

  3. 3.

Solve the linearized equation derived from Newton iteration in the finite
element space $X_{h_{k+1}}$: Find
$(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}})\in X_{h_{k+1}}$ such that
for any $(\mu,v_{h_{k+1}})\in X_{h_{k+1}}$, there holds

$\displaystyle\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle=$
$\displaystyle-\langle\mathcal{G}(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle$
$\displaystyle+\langle\mathcal{G}^{\prime}(\lambda^{h_{k}},u^{h_{k}})(\lambda^{h_{k}},u^{h_{k}}),(\mu,v_{h_{k+1}})\rangle.$

The above equation can be solved by the mixed finite element method in the
following form: Find
$(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}})\in\mathbb{R}\times
V_{h_{k+1}}$ such that for any $(\mu,v_{h_{k+1}})\in\mathbb{R}\times
V_{h_{k+1}}$, there holds

$\displaystyle\left\\{\begin{array}[]{lll}a(\lambda^{h_{k}},u^{h_{k}};\widehat{u}^{h_{k+1}},v_{h_{k+1}})+b(u^{h_{k}};v_{h_{k+1}},\widehat{\lambda}^{h_{k+1}})&=&\langle\mathcal{F}_{u}^{\prime}(\lambda^{h_{k}},u^{h_{k}})u^{h_{k}},v_{h_{k+1}}\rangle\\\
&&-\langle\mathcal{F}(\lambda^{h_{k}},u^{h_{k}}),v_{h_{k+1}}\rangle-\lambda^{h_{k}}(u^{h_{k}},v_{h_{k+1}}),\\\
b(u^{h_{k}};\widehat{u}^{h_{k+1}},\mu)&=&-\mu/2-\mu(u^{h_{k}},u^{h_{k}})/2.\end{array}\right.$

  4. 4.

Set

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=(1-\theta_{k+1})(\lambda^{h_{k}},u^{h_{k}})+\theta_{k+1}(\widehat{\lambda}^{h_{k+1}},\widehat{u}^{h_{k+1}}).$

  5. 5.

If $Resi(\lambda^{h_{k+1}},u^{h_{k+1}})\leq Resi(\lambda^{h_{k}},u^{h_{k}})$,
stop. Else set $\theta_{k+1}=\theta_{k+1}/2$ and goto Step 4.

In order to simplify the notation, we use the following symbol to denote the
above solving steps:

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=Mixing_{-}Iteration(\lambda^{h_{k}},u^{h_{k}},\theta_{k+1},X_{h_{k+1}}).$

###### Algorithm 5.3

Multigrid method based on mixing scheme

  1. 1.

Solve the following nonlinear eigenvalue problem directly in the initial space
$V_{h_{1}}$: Find $(\lambda^{h_{1}},u^{h_{1}})\in\mathbb{R}\times V_{h_{1}}$
such that $\|u^{h_{1}}\|_{0}=1$ and

$(\mathcal{A}\nabla u^{h_{1}},\nabla
v_{h_{1}})+(Vu^{h_{1}},v_{h_{1}})+(f((u^{h_{1}})^{2})u^{h_{1}},v_{h_{1}})=\lambda^{h_{1}}(u^{h_{1}},v_{h_{1}}),\quad\forall
v_{h_{1}}\in V_{h_{1}}.$

  2. 2.

For $k=1,\cdots,n-1$, obtain a new approximate solution
$(\lambda^{h_{k+1}},u^{h_{k+1}})\in X_{h_{k+1}}$ through

$\displaystyle(\lambda^{h_{k+1}},u^{h_{k+1}})=Mixing_{-}Iteration(\lambda^{h_{k}},u^{h_{k}},\theta_{k+1},X_{h_{k+1}}),$

End for.

  3. 3.

Finally, we obtain the approximate solution $(\lambda^{h_{n}},u^{h_{n}})\in
X_{h_{n}}$.

##  6 Numerical results

In this section, we propose some numerical experiments to support our
theoretical conclusions and demonstrate the efficiency of the presented
multigrid methods. For the linearized boundary value problems derived by
Newton iteration, we used the V-cycle multigrid method to obtain the numerical
solutions. The multigrid method includes two pre-smoothing steps and two post-
smoothing steps. The adopted smoother for the pre-smoothing and post-smoothing
steps is the distributive Gauss-Seidel (DGS) iteration Saddle-book .

###  6.1 Example 1

In the first example, we use Algorithm 4.2 to solve the Gross-Pitaevskii
equation: Find $(\lambda,u)$ such that

$\left\\{\begin{array}[]{rcl}-\triangle u+Vu+\zeta|u|^{2}u&=&\lambda
u,\quad{\rm in}\ \Omega,\\\ u&=&0,\ \ \quad{\rm on}\ \partial\Omega,\\\
\int_{\Omega}u^{2}d\Omega&=&1,\end{array}\right.$ (67)

where $\Omega=[0,1]^{3}$, $V=x_{1}^{2}+2x_{2}^{2}+4x_{3}^{2}$ and $\zeta=1$.

The quadratic finite element space was adopted in this example. The sequence
of meshes was produced through a one-time uniform refinement. Thus, the
refinement index $\beta$ between two consecutive meshes equals $2$. The
initial mesh is depicted in Figure 1.

Figure 1:  The initial mesh for Example 1.

Because the exact solution of (67) is unknown, we select an adequately
accurate approximation on a sufficiently fine mesh as the exact solution
($\lambda=34.819449$). The error estimates of the approximate solutions
derived by Algorithm 4.2 are presented in Figure 2. The results show that
Algorithm 4.2 is able to derive the approximate solutions with the optimal
error estimates. In order to intuitively illustrate the efficiency of
Algorithm 4.2, the computational time of Algorithm 4.2 is presented in Table
1. Besides, the computational time of the direct finite element method (solve
the nonlinear eigenvalue problem directly in the final finite element space)
is also presented. We can see that the linear computational complexity of
Algorithm 4.2 can be obtained with the refinement of mesh. In addition,
Algorithm 4.2 has a great advantage over the direct finite element method.

Figure 2: The errors of approximate solutions derived by Algorithm 4.2 for Example 1. Table 1: The CPU time of Algorithm 4.2 and the direct finite element method for Example 1.  Mesh level | Number of degrees of freedom | Time of Algorithm 4.2 | Time of direct FEM
---|---|---|---
1 | 729 | 0.1643 | 0.1643
2 | 4913 | 0.2642 | 1.1434
3 | 35937 | 1.1011 | 12.2292
4 | 274625 | 8.3957 | 151.3532
5 | 2146689 | 66.9435 | 3947.5318
6 | 16974593 | 540.2690 | -

###  6.2 Example 2

In the second example, we use Algorithm 5.3 to solve the Gross-Pitaevskii
equation (67) in $\Omega=[0,1]^{3}$, where
$W=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+sin^{2}(2\pi x_{1})+sin^{2}(2\pi
x_{2})+sin^{2}(2\pi x_{3})$ and $\zeta=100$.

In this example, due to the strong nonlinearity, Algorithm 4.2 does not
converge. Thus, we used Algorithm 5.3 to solve the nonlinear eigenvalue
problem and the convergent results are obtained.

The quadratic finite element space was adopted in this example. The sequence
of meshes was produced through a one-time uniform refinement. Thus, the
refinement index $\beta$ between two consecutive meshes equals $2$. The
initial mesh is depicted in Figure 3.

Figure 3:  The initial mesh for Example 2. Figure 4: The errors of approximate
solutions derived by Algorithm 5.3 for Example 2.

Because the exact solution of (67) is unknown, we select an adequately
accurate approximation on a sufficiently fine mesh as the exact solution
($\lambda=205.112532$). The error estimates of the approximate solutions
derived by Algorithm 5.3 are shown in Figure 4. The results show that
Algorithm 5.3 is able to derive the approximate solutions with the optimal
error estimates. In order to intuitively illustrate the efficiency of
Algorithm 5.3, the computational time of Algorithm 5.3 and the direct finite
element method are presented in Table 2. We can see that the linear
computational complexity can be obtained with the refinement of mesh. In
addition, Algorithm 5.3 has a great advantage over the direct finite element
method.

In Table 3, we also present the value of $\theta_{k}$ in each finite element
space $X_{h_{k}}$. The residuals corresponding to different $\theta_{k}$ are
also presented. Based on the adaptive strategy, we can see that the value of
the residual monotonically decreases with the refinement of mesh, and the
results are consistent with Theorem 5.1.

Table 2: The CPU time of Algorithm 5.3 and the direct finite element method for Example 2. Mesh level | Number of degrees of freedom | Time of Algorithm 5.3 | Time of direct FEM
---|---|---|---
1 | 4913 | 2.1944 | 2.1944
2 | 35937 | 3.0003 | 24.5324
3 | 274625 | 10.7252 | 289.6608
4 | 2146689 | 80.1640 | 7560.2502
5 | 16974593 | 676.3930 | -
Table 3: The value of $\theta_{k}$ in Algorithm 5.3 for Example 2.  Mesh level | The value of $\theta_{k}$ | The value of $Resi$
---|---|---
1 | 0.5 | 17.8814
2 | 0.5 | 8.9953
3 | 0.5 | 4.4768
4 | 0.5 | 2.2358
5 | 0.5 | 1.1355

##  7 Concluding remarks

In this study, we designed a novel multigrid method to solve the nonlinear
eigenvalue problem on the basis of the Newton iteration. The novel scheme
transforms the nonlinear eigenvalue problem into a series of linearized
boundary value problems defined in a sequence of product finite element
spaces. Because of avoiding solving large-scale nonlinear eigenvalue problems
directly, the overall solving efficiency is significantly improved. Besides,
the optimal error estimate and linear complexity can be derived
simultaneously. In addition, an improved multigrid method coupled with the
mixing scheme is also introduced. The improved scheme can make the iteration
scheme converge for more complicated models. More importantly, a convergence
result is derived which is missing in the existing literature on the mixing
scheme for nonlinear eigenvalue problems.

## Declarations

Conflict of interest: The authors declare that they have no conflicts of
interest/competing interests.
Code Availability: The custom codes generated during the current study are
available from the corresponding author on reasonable request.

## References

  * (1) R. A. Adams, Sobolev Spaces, Academic Press, New York, 1975.
  * (2) D. G. Anderson, Iterative procedures for nonlinear integral equations, J. Assoc. Comput. Mach., 12 (1965), 547–560.
  * (3) G. Bao, G. Hu and D. Liu, An h-adaptive finite element solver for the calculations of the electronic structures, J. Comput. Phys., 231(14) (2012), 4967–4979.
  * (4) W. Bao, Q. Du, Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, SIAM J. Sci. Comput., 25 (2004), 1674–1697.
  * (5) M. Benzi, G. H. Golub and J. Liesen, Numerical solution of saddle point problems, Acta Numer., 14 (2005), 1–137.
  * (6) A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comput., 31 (1977), 333–390.
  * (7) F. Brezzi and F. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991.
  * (8) S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994.
  * (9) Y. Cai, L. Zhang, Z. Bai and R. Li, On an eigenvector-dependent nonlinear eigenvalue problem, SIAM J. Matrix Anal. Appl., 39(3) (2018), 1360–1382.
  * (10) E. Cancès, R. Chakir and Y. Maday, Numerical analysis of nonlinear eigenvalue problems, J. Sci. Comput., 45(1-3) (2010), 90–117.
  * (11) E. Cancès, R. Chakir and Y. Maday, Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models, ESAIM Math. Model. Numer. Anal., 46 (2012), 341–388.
  * (12) H. Chen, L. He and A. Zhou, Finite element approximations of nonlinear eigenvalue problems in quantum physics, Comput. Methods Appl. Mech. Engrg., 200(21) (2011), 1846–1865.
  * (13) H. Chen, F. Liu and A. Zhou, A two-scale higher-order finite element discretization for Schrödinger equation, J. Comput. Math., 27 (2009), 315–337.
  * (14) H. Chen, H. Xie and F. Xu, A full multigrid method for eigenvalue problems, J. Comput. Phys., 322 (2016), 747–759.
  * (15) P. G. Ciarlet and J. L. Lions (eds.): Finite Element Methods, Handbook of Numerical Analysis, vol. II, North Holland, Amsterdam, 1991.
  * (16) H. Fang and Y. Saad, Two classes of multisecant methods for nonlinear acceleration, Numer. Linear Algebra Appl., 16(3) (2009), 197–221.
  * (17) R. P. Fedorenko, A relaxation method for solving elliptic difference equations, USSR Computational Mathematics and Mathematical Physics., 1(4) (1961), 1092–1096.
  * (18) W. Hackbusch, Multi-grid Methods and Applications, Springer-Verlag, Berlin, 1985
  * (19) R. Harrison, I. Moroz and K.P. Tod, A numerical study of the Schrödinger-Newton equations, Nonlinearity, 16 (2003), 101–122.
  * (20) G. Hu, H. Xie and F. Xu, A multilevel correction adaptive finite element method for Kohn-Sham equation, J. Comput. Phys., 355 (2018), 436–449.
  * (21) S. Jia, H. Xie, M. Xie and F. Xu, A full multigrid method for nonlinear eigenvalue problems, Sci. China Math., 59 (2016), 2037–2048.
  * (22) W. Kohn and L. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. A, 140 (1965), 1133–1138.
  * (23) E. H. Lieb, Thomas-Fermi and related theories of atoms and molecules, Rev. Mod. Phys., 53 (1981), 603–641.
  * (24) L. Lin and C. Yang, Elliptic preconditioner for accelerating the self-consistent field iteration in Kohn-Sham density functional theory, SIAM J. Sci. Comput., 35(5) (2013), S277–S298.
  * (25) P. A. Lott, H. F. Walker, C. S. Woodward, and U. M. Yang, An accelerated Picard method for nonlinear systems related to variably saturated flow, Adv. Water Resour., 38 (2012), 92–101.
  * (26) P. Motamarri, M. Iyer, J. Knap, V. Gavini, Higher-order adaptive finite-element methods for orbital-free density functional theory, J. Comput. Phys., 231 (2012), 6596–6621.
  * (27) S. Pollock, L. Rebholz, and M. Xiao, Anderson-accelerated convergence of picard iterations for incompressible Navier-Stokes equations, SIAM J. Numer. Anal., 57(2) (2019), 615–637.
  * (28) P. Stasiak and M.W. Matsen, Efficiency of pseudo-spectral algorithms with anderson mixing for the SCFT of periodic block-copolymer phases, Eur. Phys. J. E, 34(110) (2011), 1–9.
  * (29) V. V. Shaidurov, Multigrid Methods for Finite Element, Kluwer Academic Publics, Netherlands, 1995.
  * (30) A. Toth and C. T. Kelley, Convergence analysis for Anderson acceleration, SIAM J. Numer. Anal., 53(2) (2015), 805–819.
  * (31) J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Rev., 34 (1992), 581–613.
  * (32) J. Xu and A. Zhou, A two-grid discretization scheme for eigenvalue problems, Math. Comp., 70 (2001), 17–25.
  * (33) H. Yserentant, On the regularity of the electronic Schrödinger equation in Hilbert spaces of mixed derivatives, Numer. Math., 98(4) (2004), 731–759.
  * (34) A. Zhou, An analysis of finite-dimensional approximations for the ground state solution of Bose-Einstein condensates, Nonlinearity, 17 (2004), 541–550.
  * (35) A. Zhou, Finite dimensional approximations for the electronic ground state solution of a molecular system, Math. Methods Appl. Sci., 30 (2007), 429–447.