# Homotopy Methods for Eigenvector-Dependent Nonlinear Eigenvalue Problems
††thanks: The research was supported in part by the National Natural Science
Foundation of China (11971092)

Xuping Zhang School of Mathematical Sciences, Dalian University of Technology,
Dalian, Liaoning 116025, P. R. China (zhangxp@dlut.edu.cn). Haimei Huo School
of Mathematical Sciences, Dalian University of Technology, Dalian, Liaoning
116025, P. R. China (ab1234@mail.dlut.edu.cn).

###### Abstract

Eigenvector-dependent nonlinear eigenvalue problems are considered which arise
from the finite difference discretizations of the Gross-Pitaevskii equation.
Existence and uniqueness of positive eigenvector for both one and two
dimensional cases and existence of antisymmetric eigenvector for one
dimensional case are proved. In order to compute eigenpairs corresponding to
excited states as well as ground state, homotopies for both one and two
dimensional problems are constructed respectively and the homotopy paths are
proved to be regular and bounded. Numerical results are presented to verify
the theories derived for both one and two dimensional problems.

Key Words eigenvector-dependent nonlinear eigenvalue problem, Gross-Pitaevskii
equation, homotopy continuation method

Subject Classification(AMS):65H17, 65H20, 65N06, 65N25

##  1 Introduction

In this paper, we are concerned with the eigenvector-dependent nonlinear
eigenvalue problems resulting from the finite difference discretizations of
the Gross-Pitaevskii equation (GPE) describing Bose-Einstein condensates
(BEC). BEC are clouds of ultracold alkali-metal atoms or molecules that occupy
a single quantum state [1, 2]. The properties of a BEC at temperature $T$ much
smaller than the critical condensation temperature $T_{c}$ are usually
described by the nonlinear Schrödinger equation (NLS) for the macroscopic wave
function known as the Gross-Pitaevskii equation

$\begin{array}[]{lcl}{i\psi_{t}=-\frac{1}{2}\Delta\psi+V(x)\psi+\beta\
|\psi|^{2}\psi},&&{t>0,~{}x\in\Omega},\\\ {\psi(x,t)=0},&&{t\geq
0,~{}x\in\partial\Omega},\end{array}\\\ $ (1)

where $\psi=\psi(x,t)$ is the macroscopic wave function of the BEC,
$V(x)=\frac{1}{2}(x_{1}^{2}+x_{2}^{2}+\cdots+x_{N}^{2})$ is a typical trapping
potential, $\Omega$ is a bounded domain in $\mathbb{R}^{N}$, $N\leq 3$, and
$\beta$ positive or negative corresponds to the defocusing or focusing NLS.
Two important invariants of GPE are the normalization of the wave function

$\begin{array}[]{lcl}N(\psi)=\int_{\Omega}|\psi(x,t)|^{2}dx=1,&&t\geq
0,\end{array}$ (2)

and the energy

$\displaystyle
E(\psi(x,t))=\int_{\Omega}\left[\frac{1}{2}|\nabla\psi|^{2}+V(x)|\psi|^{2}+\frac{\beta}{2}|\psi|^{4}\right]dx=E(\psi(x,0)).$
(3)

To find stationary solution of (1), we substitute the formula
$\psi(x,t)=e^{-i\lambda t}\phi(x)$ into (1) and (2) and obtain the time-
independent Schrödinger equation with Dirichlet boundary condition and the
normalized condition

$\displaystyle\lambda\phi(x)$
$\displaystyle=-\frac{1}{2}\Delta\phi(x)+V(x)\phi(x)+\beta\phi^{3}(x),\quad{x\in\Omega},$
(4) $\displaystyle\phi(x)$ $\displaystyle=0,\quad x\in\partial\Omega,$ (5)
$\displaystyle\int_{\Omega}|\phi(x)|^{2}dx$ $\displaystyle=1,$ (6)

where $\lambda$ is the chemical potential of the condensate and $\phi(x)$ is a
real function independent of $t$ [3]. (4)-(6) is a nonlinear eigenvalue
problem. The eigenfunction corresponding to the minimum energy is called
ground state and other eigenfunctions corresponding to larger energy are
called excited states in the literature.

There have been many theoretical studies as well as numerical studies for the
time-independent Schrödinger equation. Bao and Cai [4] pointed out when
$\beta>0$, the positive ground state is unique, and if $V(x)$ is radially
symmetric in 2D, the positive ground state must be radially symmetric. Bao and
Tang [1] proposed methods by directly minimizing the energy functional via
finite element approximation to obtain the ground state and by continuation
method to obtain excited states. Edwards and Burnett [5] presented a Runge-
Kutta type method and employed it to solve the spherically symmetric time-
independent GPE. Adhikari [6] used this approach to get the ground state
solution of GPE in 2D with radial symmetry. Chang and Chien [7] and Chang,
Chien and Jeng [8] investigated stationary state solutions of $(\ref{sec 1:
transformed nonlinear problem})$ using numerical continuation method, where
$\lambda$ was treated as a continuation parameter. The solution curves
branching from the first few bifurcation points of $(\ref{sec 1: transformed
nonlinear problem})$ were numerically traced using continuation method under
the normalization condition $(\ref{sec 1:transformed constraint})$.

Since nonlinearity rather than discretization method is our main concern and
finite difference discretization will lead to a simpler nonlinear structure,
finite difference discretization is adopted in this paper. The finite
difference discretization of (4)-(6) is the following eigenvector-dependent
nonlinear eigenvalue problem,

$\begin{array}[]{c}D\varphi+\beta\varphi^{3}=\lambda\varphi,\\\
h\varphi^{\mathrm{T}}\varphi-1=0,\end{array}$ (7)

where $D=\frac{1}{2}D_{1}+V$, $D_{1}$ is the coefficient matrix corresponding
to $-\Delta$, $V$ is the diagonal matrix corresponding to the potential
$V(x)$, $\lambda$ and $\varphi$ are the unknowns, and $h$ is a constant
related to mesh size. $\varphi^{3}$ represents the vector with elements being
the corresponding elements of $\varphi$ to the power 3. This convention will
be used throughout this paper.

With respect to the theoretical aspects of eigenvector-dependent nonlinear
eigenvalue problem, [9] and [10] studied the following general nonlinear
eigen-value problem

$Ax+F(x)=\lambda x,$ (8)

where $A$ is an $n\times n$ irreducible Stieltjes matrix, i.e, an irreducible
symmetric positive definite matrix with off-diagonal entries nonpositive,
$F(x)=(f_{1}(x_{1}),\ldots,f_{n}(x_{n}))^{\mathrm{T}}$ and
$x=(x_{1},\cdots,x_{n})^{\mathrm{T}}$. The functions $f_{i}(x_{i})$ are
assumed to have the property that $f_{i}(x_{i})>0$, when $x_{i}>0$,
$i=1,\ldots,n$. It is shown that under certain conditions on $F(x)$, there
exists a positive eigenvector $x(\lambda)$ if and only if $\lambda>\mu$, where
$\mu$ is the smallest eigenvalue of $A$, and for every $\lambda>\mu$, the
positive eigenvector is unique. Moreover, such a solution is a monotone
increasing function of $\lambda$. The most popular numerical method to the
eigenvector-dependent nonlinear eigenvalue problems is the self-consistent
field (SCF) iteration, which is suitable for computing the ground state; for
instance, see [11, 12] and the references therein. In [13], inverse iteration
method was applied to solve eigenvector-dependent nonlinear eigenvalue
problems. Most of the above papers concentrate on the ground state and the
first excited state. As far as we know, there are only a few numerical works
on other excited states, such as [1, 14, 15, 16]. The main purpose of this
paper is to design algorithms for computing excited states of high energy.

Homotopy method is one of the effective methods for solving eigenvalue
problems. A great advantage of the homotopy method is that it is to a large
degree parallel, in the sense that each eigenpath is traced independently of
the others. There are several works on homotopy methods for linear eigenvalue
problems. Remarkable numerical results have been obtained by using homotopy
algorithm on eigenvalue problems of tridiagonal symmetric matrices [17, 18].
Solving eigenvalue problems of real nonsymmetric matrices with real homotopy
was developed in [19, 20]. The homotopy method is also used to solve the
generalized eigenvalue problem [21]. For eigenvalue-dependent nonlinear eigen-
problems such as $\lambda$-matrix problems, a homotopy was given by Chu, Li
and Sauer [22].

The major part of this paper is the construction of homotopy for computing
many eigenpairs of the eigenvector-dependent nonlinear eigen-problem. Key
issues encountered in constructing the homotopy are the selection of the
homotopy parameter and that of an appropriate initial eigenvalue problem so
that the homotopy paths determined by the homotopy equation are regular and
the numerical work in following these paths is at reasonable cost. The
parameter $\beta$ in the original problem seems to be a natural choice for the
homotopy parameter. However, it seems difficult to prove that 0 is a regular
value for such homotopy. In fact, 0 is probably not a regular value of the
natural homotopy with parameter $\beta$. Instead, an artificial parameter $t$
is chosen as the homotopy parameter to connect a constructed initial
eigenvalue problem and the target one. As for the selection of an initial
eigenvalue problem, random matrix with certain sparse structure is designed,
which guarantees that 0 is a regular value of the homopoty with probability
one and which renders the initial problem and the target problem possess
similar structures.

The rest of this paper is organized as follows. In Section 2, the time-
independent GPE Dirichlet problem (4)-(6) is discretized by finite difference
method and existence of certain types of solution of the discretized problems
is derived. In Section 3, homotopies for (7) are constructed with
$\Omega\subset\mathbb{R}$ and $\Omega\subset\mathbb{R}^{2}$ respectively, and
regularity and boundedness of the homotopy paths are proved. In Section 4,
numerical results are presented to verify the theoretical results derived for
$\Omega\subset\mathbb{R}$ and $\Omega\subset\mathbb{R}^{2}$ respectively.
Conclusions are drawn in the last section.

##  2 Discretizations of the nonlinear eigenvalue problem

###  2.1 Finite difference discretizations

For one dimensional problem (4)-(6) with $\Omega=[a,b]\subset\mathbb{R}$, the
grid points are $x_{j}=a+jh$, $j=0,\ldots,n+1$, where $n\in\mathbb{N}^{+}$ and
$h=\frac{b-a}{n+1}$ is the mesh size. The finite difference discretization of
the differential equation and a simple quadrature of the normalization
condition lead to the following system of algebraic equations,

$\begin{array}[]{c}D\varphi+\beta\varphi^{3}-\lambda\varphi=0,\\\
\frac{1}{2}\left(\frac{1}{h}-\varphi^{\mathrm{T}}\varphi\right)=0,\end{array}$
(9)

where $\varphi=(\varphi_{1},\cdots,\varphi_{n})^{\mathrm{T}}$, $\varphi_{j}$
are the approximations of $\phi(x_{j})$, $v_{j}=V(x_{j})$, $j=1,\ldots,n$, and
$D=\frac{1}{2}D_{1}+V$ with

$\displaystyle D_{1}=\frac{1}{h^{2}}\left(\begin{array}[]{cccc}2&-1&\\\
-1&2&\ddots&\\\ &\ddots&\ddots&-1\\\ &&-1&2\end{array}\right),\quad
V=\left(\begin{array}[]{ccc}v_{1}&&\\\ &\ddots&\\\ &&v_{n}\end{array}\right).$
(17)

The discretization of the normalized condition is rewritten so that the
Jacobian matrix of the nonlinear mapping with respect to $(\varphi,\lambda)$
is symmetric, as will be seen below.

For two dimensional problem (4)-(6) with $\Omega=[a,b]\times[c,d]$, the domain
is divided into a $(m+1)\times(n+1)$ mesh with step size
$h_{1}=\frac{b-a}{m+1}$ in $x$-direction, $h_{2}=\frac{d-c}{n+1}$ in
$y$-direction. The grid points $(x_{i},y_{j})$ are $x_{i}=a+ih_{1}$,
$i=0,\ldots,m+1$, and $y_{j}=c+jh_{2}$, $j=0,\ldots,n+1$. Using central
difference, we get

$\begin{array}[]{c}D\varphi+\beta\varphi^{3}-\lambda\varphi=0,\\\
\frac{1}{2}\left(\frac{1}{h_{1}h_{2}}-\varphi^{\mathrm{T}}\varphi\right)=0,\end{array}$
(18)

where
$\varphi=(\varphi_{11},\cdots,\varphi_{1n},\cdots,\varphi_{m1},\cdots,\varphi_{mn})^{\mathrm{T}}$,
$\varphi_{ij}$ are the approximations of $\phi(x_{i},y_{j})$,
$v_{ij}=V(x_{i},y_{j})$, $i=1,\ldots,m$, $j=1,\ldots,n$ and D is a block
tridiagonal matrix

$D=\frac{1}{2}\left(\begin{array}[]{cccc}D_{11}&D_{12}&\\\
D_{21}&D_{22}&\ddots&\\\ &\ddots&\ddots&D_{{m-1},m}\\\
&&D_{m,{m-1}}&D_{mm}\end{array}\right)+\left(\begin{array}[]{cccc}V_{1}&&\\\
&V_{2}&&\\\ &&\ddots&\\\ &&&V_{m}\end{array}\right),$

where

$\displaystyle
D_{ii}=\left(\begin{array}[]{cccc}\frac{2}{h_{1}^{2}}+\frac{2}{h_{2}^{2}}&-\frac{1}{h_{2}^{2}}&\\\
-\frac{1}{h_{2}^{2}}&\frac{2}{h_{1}^{2}}+\frac{2}{h_{2}^{2}}&\ddots&\\\
&\ddots&\ddots&-\frac{1}{h_{2}^{2}}\\\
&&-\frac{1}{h_{2}^{2}}&\frac{2}{h_{1}^{2}}+\frac{2}{h_{2}^{2}}\end{array}\right)\in\mathbb{R}^{n\times
n},$ (23)

$\displaystyle D_{i-1,i}=\left(\begin{array}[]{cccc}-\frac{1}{h_{1}^{2}}&&&\\\
&-\frac{1}{h_{1}^{2}}&&\\\ &&\ddots&\\\
&&&-\frac{1}{h_{1}^{2}}\end{array}\right)\in\mathbb{R}^{n\times n},\quad
D_{{i-1},i}=D_{i,{i-1}},$ (28)

$\displaystyle V_{i}=\left(\begin{array}[]{cccc}v_{i1}&&&\\\ &v_{i2}&&\\\
&&\ddots&\\\ &&&v_{in}\end{array}\right),\quad i=1,\ldots,m.$ (33)

###### Remark 2.1

For both one and two dimensional cases, D is an irreducible symmetric diagonal
dominant matrix and the diagonal entries of D are all positive. Therefore D is
positive definite.

###  2.2 Existence of certain types of solution

In this subsection, we will study the existence of solution for the
discretized nonlinear eigenvalue problem. From [4], we know for (4)-(6), when
$\beta>0$, the positive ground state is unique, and if $V(x)$ is radially
symmetric in 2D, the positive ground state must be radially symmetric. We will
prove the existence of positive solution and the existence of antisymmetric
solution for discretized nonlinear eigenvalue problem. For convenient reading,
two underlying theorems from [9] are quoted as underlying lemmas.

###### Lemma 2.2

([9]) Let A be an irreducible Stieltjes matrix and $\mu$ be the smallest
positive eigenvalue of A. Let $\lambda>\mu$ and let

$F(x)=\left(\begin{array}[]{c}f_{1}(x_{1})\\\ \vdots\\\
f_{n}(x_{n})\end{array}\right),$ (34)

where for $i=1,\ldots,n$, $f_{i}(x):[0,\infty)\rightarrow[0,\infty)$ are
$C^{1}$ functions satisfying the conditions:

$\lim\limits_{t\rightarrow
0}\frac{f_{i}(t)}{t}=0,\qquad\lim\limits_{t\rightarrow\infty}\frac{f_{i}(t)}{t}=\infty.$
(35)

Then $Ax+F(x)=\lambda x$ has a positive solution. If, in addition, for
$i=1,\ldots,n$,

$\frac{f_{i}(s)}{s}<\frac{f_{i}(t)}{t}$ (36)

whenever $0<s<t$, then the solution is unique.

###### Lemma 2.3

([9]) Let the conditions $(\ref{sec 2: limit conditon})$ and $(\ref{sec 2:
inequality condition})$ of Lemma 2.2 be satisfied and let $x(\lambda)$ denote
the unique positive eigenvector corresponding to $\lambda\in(\mu,\infty)$.
Then:

  1. (i)

$x(\lambda_{1})<x(\lambda_{2})$, if $\mu<\lambda_{1}<\lambda_{2}<\infty$;

  2. (ii)

$x(\lambda)$ is continuous on $(\mu,\infty)$;

  3. (iii)

$\lim\limits_{\lambda\rightarrow\infty}{x_{i}(\lambda)}=\infty$,
$i=1,\cdots,n$;

  4. (iv)

$\lim\limits_{\lambda\rightarrow\mu^{+}}{x_{i}(\lambda)}=0$, $i=1,\cdots,n$.

###### Remark 2.4

Lemma 2.3 indicates for any given normalization $r>0$, there exist a unique
$\lambda>\mu$ and unique positive $x(\lambda)$ such that $||x(\lambda)||=r$.

###### Theorem 2.5

If $\beta>0$, there exist unique positive eigenvectors for problem $(\ref{sec
2: finite difference - one dimen})$ and $(\ref{sec 2: finite difference - two
dimen})$ respectively.

Proof. From Remark 2.1, we know $D$ in both $(\ref{sec 2: finite difference -
one dimen})$ and $(\ref{sec 2: finite difference - two dimen})$ is an
irreducible Stieltjes matrix. In addition it can be verified that
$\beta\varphi^{3}$ in both $(\ref{sec 2: finite difference - one dimen})$ and
$(\ref{sec 2: finite difference - two dimen})$ satisfies the conditions of
$F(x)$ in Lemma 2.2. From Remark 2.4, the claim is proved.

###### Theorem 2.6

Let $\beta>0$ and $\Omega=[-a,a]$ for problem $(\ref{sec 2: finite difference
- one dimen})$.

  1. (i)

When $n$ is odd and the grid points $x_{i}$, $i=1,\ldots,n$, satisfy

$x_{1}=-x_{n},~{}x_{2}=-x_{n-1},~{}\cdots,~{}x_{\frac{n-1}{2}}=-x_{\frac{n+3}{2}},~{}x_{\frac{n+1}{2}}=0,$

there exists a unique solution
$\varphi=(\varphi_{1},\varphi_{2},\cdots,\varphi_{n})^{\mathrm{T}}$ with
$\varphi_{1}=-\varphi_{n}$, $\varphi_{2}=-\varphi_{n-1},\cdots$,
$\varphi_{\frac{n-1}{2}}=-\varphi_{\frac{n+3}{2}}$,
$\varphi_{\frac{n+1}{2}}=0$, $\varphi_{j}>0$, $j=1,\ldots,n$.

  2. (ii)

When $n$ is even and the grid points $x_{i}$, $i=1,\ldots,n$, satisfy

$x_{1}=-x_{n},~{}x_{2}=-x_{n-1},~{}\cdots,~{}x_{\frac{n}{2}}=-x_{\frac{n}{2}+1},$

there exists a unique solution
$\varphi=(\varphi_{1},\varphi_{2},\cdots,\varphi_{n})^{\mathrm{T}}$ with
$\varphi_{1}=-\varphi_{n}$, $\varphi_{2}=-\varphi_{n-1}$, $\cdots$,
$\varphi_{\frac{n}{2}}=-\varphi_{\frac{n}{2}+1}$, $\varphi_{j}>0$,
$j=1,\ldots,n$.

Proof. (i). When $n$ is odd, consider the following equations

$\begin{array}[]{c}D_{2}\varphi+\beta\varphi^{3}=\lambda\varphi,\\\
\varphi^{\mathrm{T}}\varphi-\frac{1}{2h}=0,\end{array}$ (37)

where

$\displaystyle
D_{2}=\left(\begin{array}[]{cccc}\frac{1}{h^{2}}+\frac{1}{2}x_{1}^{2}&-\frac{1}{2h^{2}}&&\\\
-\frac{1}{2h^{2}}&\frac{1}{h^{2}}+\frac{1}{2}x_{2}^{2}&\ddots&\\\
&\ddots&\ddots&-\frac{1}{2h^{2}}\\\
&&-\frac{1}{2h^{2}}&\frac{1}{h^{2}}+\frac{1}{2}x_{\frac{n-1}{2}}^{2}\end{array}\right)\in\mathbb{R}^{\frac{n-1}{2}\times\frac{n-1}{2}}.$
(42)

Note that $D_{2}$ is an irreducible Stieltjes matrix and $\beta\varphi^{3}$
satisfies the conditions of $F(x)$ in Lemma 2.2. Therefore there exists a
unique positive solution
$\varphi=(\varphi_{1},\varphi_{2},\ldots,\varphi_{\frac{n-1}{2}})^{\mathrm{T}}$
for (37). Due to the relations
$x_{1}=-x_{n},~{}x_{2}=-x_{n-1},~{}\cdots,~{}x_{\frac{n-1}{2}}=-x_{\frac{n+3}{2}},~{}x_{\frac{n+1}{2}}=0$,
set $\varphi_{n}=-\varphi_{1}$, $\varphi_{n-1}=-\varphi_{2}$, $\cdots$,
$\varphi_{\frac{n+3}{2}}=-\varphi_{\frac{n-1}{2}}$,
$\varphi_{\frac{n+1}{2}}=0$. Then
$\varphi=(\varphi_{1},\varphi_{2},\ldots,\varphi_{n})^{\mathrm{T}}$ is a
solution of $(\ref{sec 2: finite difference - one dimen})$.

(ii). When $n$ is even, consider the following equations

$\begin{array}[]{c}D_{2}\varphi+\beta\varphi^{3}=\lambda\varphi,\\\
\varphi^{\mathrm{T}}\varphi-\frac{1}{2h}=0,\end{array}$ (43)

where

$\displaystyle
D_{2}=\left(\begin{array}[]{cccc}\frac{1}{h^{2}}+\frac{1}{2}x_{1}^{2}&-\frac{1}{2h^{2}}&&\\\
-\frac{1}{2h^{2}}&\frac{1}{h^{2}}+\frac{1}{2}x_{2}^{2}&\ddots&\\\
&\ddots&\ddots&-\frac{1}{2h^{2}}\\\
&&-\frac{1}{2h^{2}}&\frac{3}{2h^{2}}+\frac{1}{2}x_{\frac{n}{2}}^{2}\end{array}\right)\in\mathbb{R}^{\frac{n}{2}\times\frac{n}{2}}.$
(48)

Similarly there exists a unique positive solution
$\varphi=(\varphi_{1},\varphi_{2},\ldots,\varphi_{\frac{n}{2}})^{\mathrm{T}}$
for (43). Due to the relations
$x_{1}=-x_{n},~{}x_{2}=-x_{n-1},~{}\cdots,~{}x_{\frac{n}{2}}=-x_{\frac{n}{2}+1}$,
set $\varphi_{n}=-\varphi_{1}$, $\varphi_{n-1}=-\varphi_{2}$, $\cdots$,
$\varphi_{\frac{n}{2}+1}=-\varphi_{\frac{n}{2}}$. Then
$\varphi=(\varphi_{1},\varphi_{2},\ldots,\varphi_{n})^{\mathrm{T}}$ is a
solution of $(\ref{sec 2: finite difference - one dimen})$.

##  3 Homotopy methods

In this section, in order to compute many eigenpairs, we construct homotopy
equations for 1D discretized problem $(\ref{sec 2: finite difference - one
dimen})$ and 2D discretized problem $(\ref{sec 2: finite difference - two
dimen})$ respectively. We shall prove the regularity and boundedness of the
homotopy paths. The regularity of homotopy paths can be usually obtained by
random perturbations of appropriate parameters, so the most important feature
of our construction is the choice of appropriate parameters. In addition, if
the initial matrix can be chosen as close to the matrix $D$ as possible, then
most of the homotopy paths are close to straight lines and will be easy to
follow [19].

###  3.1 One dimensional case

For $(\ref{sec 2: finite difference - one dimen})$, a homotopy
$H:\mathbb{R}^{n}\times\mathbb{R}\times[0,1]\rightarrow\mathbb{R}^{n}\times\mathbb{R}$
is defined as

$H(\varphi,\lambda,t)=\left(\begin{array}[]{c}(1-t)A(K)\varphi+D\varphi+t\beta\varphi^{3}-\lambda\varphi\\\
\frac{1}{2}\left(\frac{1}{h}-\varphi^{\mathrm{T}}\varphi\right)\end{array}\right)=0,$
(49)

where $A(K)=\mbox{diag}(a_{1},\cdots,a_{n})$ is a random diagonal matrix with
$K=(a_{1},\cdots,a_{n})^{\mathrm{T}}\in\mathbb{R}^{n}$. At $t=0$,
$H(\varphi,\lambda,0)$ corresponds to the linear eigenvalue problem

$H(\varphi,\lambda,0)=\left(\begin{array}[]{c}A(K)\varphi+D\varphi-\lambda\varphi\\\
\frac{1}{2}\left(\frac{1}{h}-\varphi^{\mathrm{T}}\varphi\right)\end{array}\right)=0,$
(50)

while at $t=1$, $H(\varphi,\lambda,1)=0$ corresponds to the problem (9).
Assuming that the eigenpairs of $H(\varphi,\lambda,0)=0$ are
$(\varphi^{(i)},\lambda_{i})$, $i=1,\ldots,n$,
$\lambda_{1}\leq\cdots\leq\lambda_{n}$, we shall use these $n$ points
$(\varphi^{(i)},\lambda_{i},0)$ as our initial points when tracing the
homotopy curves leading to the desired solutions of $H(\varphi,\lambda,1)=0$.

The choice of the initial matrix $A(K)+D$ provides some advantages. First, it
makes sure that $\forall t\in[0,1)$, what we need to solve is a sparse
nonlinear eigenvalue problem with constraint. Second, since $\forall
K\in\mathbb{R}^{n}$, $\forall t\in[0,1)$ and $\varphi\in\mathbb{R}^{n}$,
$(1-t)A(K)+D+t\beta\mbox{diag}(\varphi^{2})$ is a real symmetric matrix, the
solution curves starting from the initial points at $t=0$ are real. Finally,
since $A(K)+D$ is a tridiagonal matrix with all the subdiagonal and
supdiagonal entries nonzero, all the eigenvalues of $A(K)+D$ are simple and
the Jacobian matrix of $H$ at $(\varphi_{0},\lambda_{0},0)$ is nonsingular for
$\forall(\varphi_{0},\lambda_{0})$ such that $H(\varphi_{0},\lambda_{0},0)=0$.
Thus locally a unique curve around $(\varphi_{0},\lambda_{0},0)$ is
guaranteed.

The effectiveness of the homotopy is based on the following Parametrized
Sard’s Theorem.

###### Theorem 3.1

(Parametrized Sard’s Theorem) Let $f:M\times
P\subset\mathbb{R}^{m}\times\mathbb{R}^{q}\to\mathbb{R}^{n}$ be a $C^{k}$
mapping with $k>max(0,m-n)$, where $M$ and $P$ are open sets in
$\mathbb{R}^{m}$ and $\mathbb{R}^{q}$ respectively. If $y$ is a regular value
of f, then y is also a regular value of $f(\cdot,p)$ for almost all $p\in P$.

In the rest of this paper, we will denote the $i$-th row of a matrix $M$ as
$M(i,:)$ and the $j$-th column of $M$ as $M(:,j)$. If $I_{1}$ is a row index
set, $M(I_{1},:)$ will be the submatrix formed by the $I_{1}$ rows of $M$.
$M(i:end,:)$ will be the submatrix formed by the rows from the $i$-th row to
the last. Similarly $M(:,J_{1})$ will be the submatrix formed by the $J_{1}$
columns of $M$. If $I_{1}$ and $I_{2}$ are two index sets,
$[M(I_{1},:);M(I_{2},:)]$ denotes the submatrix formed by the $I_{1}$ rows and
the $I_{2}$ rows of $M$. In Theorem 3.2, we prove regularity and boundedness
of the homotopy paths determined by the homotopy equation (49).

###### Theorem 3.2

For the homotopy
$H:\mathbb{R}^{n}\times\mathbb{R}\times[0,1)\rightarrow\mathbb{R}^{n}\times\mathbb{R}$
in (49), for almost all $K\in\mathbb{R}^{n}$,

  1. (i)

0 is a regular value of $H$ and therefore the homotopy curves corresponding to
different initial points do not intersect each other for $t\in[0,1)$;

  2. (ii)

Every homotopy path $(\varphi(s),\lambda(s),t(s))\subset H^{-1}(0)$ is
bounded.

Proof. (i) Define a mapping
$\tilde{H}:\mathbb{R}^{n}\times\mathbb{R}\times[0,1]\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\times\mathbb{R}$
related to $H$ as follows

$\tilde{H}(\varphi,\lambda,t,K)=\left(\begin{array}[]{c}(1-t)A(K)\varphi+D\varphi+t\beta\varphi^{3}-\lambda\varphi\\\
\frac{1}{2}\left(\frac{1}{h}-\varphi^{\mathrm{T}}\varphi\right)\end{array}\right)=0,$
(51)

such that $H(\varphi,\lambda,t)=\tilde{H}(\varphi,\lambda,t,K)$. The Jacobian
matrix of $\tilde{H}$,
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$, is

$\left(\begin{array}[]{cccc}(1-t)A(K)+D+3t\beta\mbox{diag}(\varphi^{2})-\lambda
I&-\varphi&\beta\varphi^{3}-A(K)\varphi&(1-t)\mbox{diag}(\varphi)\\\
-\varphi^{\mathrm{T}}&0&0&0\end{array}\right).$

Divide $\\{1,\cdots,n\\}$ into two parts $C^{0}$ and $C^{*}$, where $C^{0}$
denotes the indices $i$ such that $\varphi_{i}=0$ and $C^{*}$ denotes the
indices $i$ such that $\varphi_{i}\neq 0$. Since $\forall i\in C^{0}$, both
the columns $3t\beta\mbox{diag}(\varphi^{2})(:,i)$ and
$t\beta\mbox{diag}(\varphi^{2})(:,i)$ are equal to zero, it holds that,

$\displaystyle\left((1-t)A(K)+D+3t\beta\mbox{diag}(\varphi^{2})-\lambda
I\right)(:,C^{0})$ $\displaystyle=$
$\displaystyle\left((1-t)A(K)+D+t\beta\mbox{diag}(\varphi^{2})-\lambda
I\right)(:,C^{0}).$ (52)

For the columns in $C^{*}$, the diagonal matrix $\mbox{diag}(\varphi)$ is
nonzero. By elementary column transformations, the Jacobian matrix
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$ is transformed to
the following,

$\left(\begin{array}[]{cccc}(1-t)A(K)+D+t\beta\mbox{diag}(\varphi^{2})-\lambda
I&-\varphi&\beta\varphi^{3}-A(K)\varphi&(1-t)\mbox{diag}(\varphi)\\\
-\varphi^{\mathrm{T}}&0&0&0\end{array}\right).$

Define

$\displaystyle
F_{1}=\left(\begin{array}[]{cc}(1-t)A(K)+D+t\beta\mbox{diag}(\varphi^{2})-\lambda
I&-\varphi\\\ -\varphi^{\mathrm{T}}&0\end{array}\right).$ (55)

For $\forall K\in\mathbb{R}^{n}$,
$\forall(\varphi,\lambda,t)\in\mathbb{R}^{n}\times R\times[0,1)$ satisfying
$(\ref{sec3: mapping-related-to-homotopy-1D})$, since the subdiagonals and
supdiagonals of the matrix $(1-t)A(K)+D+t\beta\mbox{diag}(\varphi^{2})$ are
nonzero, $\lambda$ is a simple eigenvalue of
$(1-t)A(K)+D+t\beta\mbox{diag}(\varphi^{2})$. Therefore $F_{1}$ is nonsingular
and $\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$ is row full
rank. As a result, 0 is a regular value of the mapping $\tilde{H}$. From
Theorem 3.1, for almost all $K\in\mathbb{R}^{n}$, 0 is a regular value of the
restricted mapping $\tilde{H}(\cdot,\cdot,\cdot,K)$, i.e., $H$.

(ii) From $(\ref{sec 3: homotopy one dimen})$, for fixed $K\in\mathbb{R}^{n}$,
$\forall(\varphi,\lambda,t)\in\mathbb{R}^{n}\times\mathbb{R}\times[0,1)$
satisfying $H(\varphi,\lambda,t)$$=$$0$, $\|\varphi\|=1/\sqrt{h}$ and
$\lambda=h\left((1-t)\varphi^{\mathrm{T}}A(K)\varphi+\varphi^{\mathrm{T}}D\varphi+t\beta\varphi^{\mathrm{T}}\varphi^{3}\right)$.
Then

$\displaystyle|\lambda|$ $\displaystyle=$ $\displaystyle
h|(1-t)\varphi^{\mathrm{T}}A(K)\varphi+\varphi^{\mathrm{T}}D\varphi+t\beta\varphi^{\mathrm{T}}\varphi^{3}|$
$\displaystyle\leq$ $\displaystyle(1-t)\rho(A(K))+\rho(D)+\frac{t\beta}{h}$
$\displaystyle\leq$ $\displaystyle\rho(A(K))+\rho(D)+\frac{\beta}{h},$

where $\rho(M)$ denotes the spectral radius of $M$.

###  3.2 Two dimensional case

####  3.2.1 The homotopy with random tridiagonal matrix

For two dimensional case, the homotopy
$H:\mathbb{R}^{mn}\times\mathbb{R}\times[0,1]\rightarrow\mathbb{R}^{mn}\times\mathbb{R}$
is constructed as follows

$H(\varphi,\lambda,t)=\left(\begin{array}[]{c}(1-t)A(K)\varphi+D\varphi+t\beta\varphi^{3}-\lambda\varphi\\\
\frac{1}{2}\left(\frac{1}{h_{1}h_{2}}-\varphi^{\mathrm{T}}\varphi\right)\end{array}\right)=0,$
(56)

where $A(K)\in\mathbb{R}^{mn\times mn}$ is a random block diagonal matrix with
tridiagonal blocks, namely,

$A(K)=\left(\begin{array}[]{cccc}A_{1}&&&\\\ &A_{2}&&\\\ &&\ddots&\\\
&&&A_{m}\end{array}\right)\mbox{with}~{}A_{i}=\left(\begin{array}[]{ccccc}a_{11}^{(i)}&a_{12}^{(i)}&&&\\\
a_{12}^{(i)}&a_{22}^{(i)}&a_{23}^{(i)}&&\\\
&a_{23}^{(i)}&a_{33}^{(i)}&\ddots&\\\ &&\ddots&\ddots&a_{n-1,n}^{(i)}\\\
&&&a_{n-1,n}^{(i)}&a_{nn}^{(i)}\end{array}\right),~{}i=1,\ldots,m,$

and $K=\left(K_{1},K_{2},\ldots,K_{m}\right)^{\mathrm{T}}$ with
$K_{i}=\left(a_{11}^{(i)},a_{12}^{(i)},a_{22}^{(i)},a_{23}^{(i)},\ldots,a_{n-1,n-1}^{(i)},a_{n-1,n}^{(i)},a_{nn}^{(i)}\right)$.
$H(\varphi,\lambda,0)$ corresponds to the linear eigenvalue problem

$H(\varphi,\lambda,0)=\left(\begin{array}[]{c}A(K)\varphi+D\varphi-\lambda\varphi\\\
\frac{1}{2}\left(\frac{1}{h_{1}h_{2}}-\varphi^{\mathrm{T}}\varphi\right)\end{array}\right)=0,$

while $H(\varphi,\lambda,1)$ corresponds to the problem $(\ref{sec 2: finite
difference - two dimen})$.

In order to show the effectiveness of this homotopy $H$ by the Parametrized
Sard’s Theorem, define a mapping
$\tilde{H}:\mathbb{R}^{mn}\times\mathbb{R}\times[0,1]\times\mathbb{R}^{(2nm-m)}\rightarrow\mathbb{R}^{mn}\times\mathbb{R}$
related to $H$ as follows

$\tilde{H}(\varphi,\lambda,t,K)=\left(\begin{array}[]{c}(1-t)A(K)\varphi+D\varphi+t\beta\varphi^{3}-\lambda\varphi\\\
\frac{1}{2}\left(\frac{1}{h_{1}h_{2}}-\varphi^{\mathrm{T}}\varphi\right)\end{array}\right)=0,$
(57)

such that $H(\varphi,\lambda,t)=\tilde{H}(\varphi,\lambda,t,K)$. The Jacobian
matrix of $\tilde{H}$,
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$, is

$\left(\begin{array}[]{cccc}(1-t)A(K)+D+3t\beta\mbox{diag}(\varphi^{2})-\lambda
I&-\varphi&\beta\varphi^{3}-A(K)\varphi&(1-t)B\\\
-\varphi^{\mathrm{T}}&0&0&0\end{array}\right),$ (58)

where $B=\frac{\partial(A(K)\varphi)}{\partial
K}\in\mathbb{R}^{{(mn)}\times{(2nm-m)}}$. Denote
$\varphi_{i}=\left(\varphi_{i1},\varphi_{i2},\ldots,\varphi_{in}\right)^{\mathrm{T}}$.
It can be verified that

$B=\left(\begin{array}[]{cccc}B_{1}&&&\\\ &B_{2}&&\\\ &&\ddots&\\\
&&&B_{m}\end{array}\right)$ (59)

with

$B_{i}=\frac{\partial(A_{i}\varphi_{i})}{\partial
K_{i}}=\left(\begin{array}[]{ccccccccccc}\varphi_{i1}&\varphi_{i2}&&&&&&&&&\\\
&\varphi_{i1}&\varphi_{i2}&\varphi_{i3}&&&&&&&\\\
&&&\varphi_{i2}&\varphi_{i3}&\varphi_{i4}&&&&&\\\ &&&&&&\ddots&&&&\\\
&&&&&&&\varphi_{i,n-2}&\varphi_{i,n-1}&\varphi_{in}&\\\
&&&&&&&&&\varphi_{i,n-1}&\varphi_{in}\end{array}\right).$ (60)

For example, when $m=2,n=3$, we have

$B=\begin{pmatrix}\varphi_{11}&\varphi_{12}&{0}&{0}&{0}&{0}&{0}&{0}&{0}&{0}\\\
{0}&\varphi_{11}&\varphi_{12}&\varphi_{13}&{0}&{0}&{0}&{0}&{0}&{0}\\\
{0}&{0}&{0}&\varphi_{12}&\varphi_{13}&{0}&{0}&{0}&{0}&{0}\\\
{0}&{0}&{0}&{0}&{0}&\varphi_{21}&\varphi_{22}&{0}&{0}&{0}\\\
{0}&{0}&{0}&{0}&{0}&{0}&\varphi_{21}&\varphi_{22}&\varphi_{23}&{0}\\\
{0}&{0}&{0}&{0}&{0}&{0}&{0}&{0}&\varphi_{22}&\varphi_{23}\end{pmatrix}.$

Denote by $G_{a}$ all the inner grid points and by $R_{a}$ the ordering of the
grid points in $G_{a}$, i.e.,

$\displaystyle G_{a}=\\{(i,j):i=1,\ldots,m,j=1,\ldots,n\\},$ (61)
$\displaystyle R_{a}=\\{j+(i-1)*n:(i,j)\in G_{a}\\},$ (62)

and by $G^{0}$ the grid points with function value being zero,

$\displaystyle G^{0}=\\{(i,j)\in G_{a}:\varphi_{ij}=0\\}.$ (63)

Note that $G_{a}$ and $R_{a}$ have a one to one correspondence. Define such
correspondence as a mapping $\Gamma:G_{a}\rightarrow R_{a}$,

$\displaystyle\Gamma(i,j)=j+(i-1)*n.$ (64)

Denote by $R^{0}$ the indices of rows in which $B$ is zero, by $R^{*}$ the
indices of rows, in which $B$ is not zero, and $S_{i}^{0}$ and $S_{i}^{*}$
with similar meanings for $B_{i}$,

$\displaystyle R^{0}=\\{r:B(r,:)=0\\},\quad R^{*}=\\{r:B(r,:)\neq 0\\},$ (65)
$\displaystyle S^{0}_{i}=\\{r:B_{i}(r,:)=0\\},\quad
S_{i}^{*}=\\{r:B_{i}(r,:)\neq 0\\}.$ (66)

It can be verified that

$\displaystyle R^{0}=\bigcup_{i=1}^{m}S^{0}_{i},\quad
R^{*}=\bigcup_{i=1}^{m}S_{i}^{*}.$ (67)

###### Lemma 3.3

Let $B$ be the matrix defined in (59). Then the nonzero rows of $B$ are
linearly independent, i.e., the submatrix $B(R^{*},:)$ is row full rank.

Proof. Note that $B(R^{*},:)=[B(S_{1}^{*},:);\ldots;B(S_{m}^{*},:)]$. Due to
the block structure of $B$, it suffices to prove that for any $i$, $1\leq
i\leq m$, $B_{i}(S_{i}^{*},:)$ is row full rank if $S_{i}^{*}$ is not empty.

Now suppose that $S_{i}^{*}$ is not empty, that is $\varphi_{i}\neq 0$. Let
the nonzero components of $\varphi_{i}$ be
$\varphi_{i,i_{1}},\varphi_{i,i_{2}},\ldots,\varphi_{i,i_{r}}$, with
$i_{1}<i_{2}<\ldots<i_{r}$. Denote by
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ the rows of $B_{i}$ containing
$\varphi_{i,i_{1}},\varphi_{i,i_{2}},\ldots,\varphi_{i,i_{k}}$, $1\leq k\leq
r$. Then $B_{i}(S_{i}^{*},:)=Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{r}})$.
The claim will be proved by successively adding rows with nonzero component of
$\varphi_{i}$ to $Z$.

(1) Prove $Z(\varphi_{i,i_{1}})$ is row full rank. If $i_{1}=1$,
$Z(\varphi_{i,i_{1}})$ is the first two rows of $B_{i}$ and it is row full
rank. If $1<i_{1}<n$, $Z(\varphi_{i,i_{1}})$ is three successive rows of
$B_{i}$ which contains the following submatrix involving $\varphi_{i,i_{1}}$,

$\left(\begin{array}[]{ccc}\varphi_{i,i_{1}}&0&0\\\ &\varphi_{i,i_{1}}&*\\\
0&0&\varphi_{i,i_{1}}\end{array}\right)$

where $*$ stands for an element which may be zero or nonzero. Therefore
$Z(\varphi_{i,i_{1}})$ is row full rank. If $i_{1}=n$, $Z(\varphi_{i,i_{1}})$
is the last two rows of $B_{i}$ and is row full rank too.

(2) Prove that when $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ is row
full rank, $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$ is also row full
rank, $1\leq k<r$. If $i_{k}=n-1$, then $i_{k+1}=n$ and
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})=Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$.
Next suppose $i_{k}<n-1$. If $i_{k+1}=i_{k}+1$, i.e., $\varphi_{i,i_{k}}$ and
$\varphi_{i,i_{k+1}}$ are successive, then
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$ consists of
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ and a new row with
$\varphi_{i,i_{k+1}}$ as its first nonzero element and with
$\varphi_{i,i_{k+1}}$ located in a column different from those of
$\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}}$, the connecting submatrix
illustrated in the following,

$\left(\begin{array}[]{cccccc}*&\varphi_{i,i_{k}}&\varphi_{i,i_{k+1}}&0&0&0\\\
0&0&\varphi_{i,i_{k}}&\varphi_{i,i_{k+1}}&*&0\\\
0&0&0&0&\varphi_{i,i_{k+1}}&*\end{array}\right)$

Thus $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$ is row full rank. If
$i_{k+1}=i_{k}+2<n$, then $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$
consists of $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ and two new rows
and similarly these two new rows are linearly independent with
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$. If $i_{k+1}\geq i_{k}+3$ and
$i_{k+1}<n$, then $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$ consists
of $Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ and three new rows. If
$i_{k+1}=n$, $i_{k+1}=i_{k}+2$, then
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$ consists of
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ and a new row. If $i_{k+1}=n$,
$i_{k+1}=i_{k}+s$ and $s\geq 3$, then
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k+1}})$ consists of
$Z(\varphi_{i,i_{1}},\ldots,\varphi_{i,i_{k}})$ and two new rows. All the
latter cases can be similarly proved.

In the following, first we prove there exists a zero measure set
$U_{1}\subset\mathbb{R}^{(2nm-m)}$, such that if
$K\in\mathbb{R}^{(2nm-m)}\setminus U_{1}$, the eigenvalues of $A(K)+D$ are
simple. Then we prove that 0 is a regular value of $H(\varphi,\lambda,t)$ for
almost all $K\in\mathbb{R}^{(2nm-m)}\setminus(U_{1}\cup U_{2})$, where
$U_{2}=(\mathbb{R}^{+})^{(2nm-m)}$. The removal of $U_{2}$ is to make the
elements in the subdiagonal and supdiagonal of the matrix $(1-t)A(K)+D$
negative for $t\in[0,1)$.

###### Lemma 3.4

The eigenvalues of $A(K)+D$ are simple for $K$ almost everywhere in
$\mathbb{R}^{(2nm-m)}$ except on a subset of real codimension 1.

Proof. Let $f(\lambda)=det(A(K)+D-\lambda I)$. The polynomial $f(\lambda)$ has
no multiple roots if and only if its discriminant $R(K)$ is nonzero [22]. It
is obvious that $R(K)$ is not identically zero. Furthermore, since $R(K)$ is a
polynomial in the elements of vector $K$, it can vanish only on a hypersurface
of real codimension 1. The hypersurface is

$U_{1}=\\{K|R(K)=0\\}.$

In the following Lemma 3.5, we prove that for any matrix $F$ consisting of
several submatrices by row, if every submatrix of $F$ is row full rank and the
index sets of nonzero columns do not intersect for any two submatrices, then
$F$ is row full rank.

###### Lemma 3.5

Let $F\in\mathbb{R}^{p\times q}$ be a matrix. Suppose
$\\{1,\ldots,p\\}=\bigcup_{i=1}^{s}I_{i}$, with $I_{i}\bigcap
I_{j}=\emptyset$, if $i\neq j$. Denote

$E_{i}=\\{k\in\\{1,\ldots,q\\}:F(I_{i},k)\neq 0\\}.$

Suppose that for any $1\leq i\leq s$, $F(I_{i},:)$ is row full rank and
$E_{i}\bigcap E_{j}=\emptyset$, if $i\neq j$. Then $F$ is row full rank.

Proof. Since for any $1\leq i\leq s$, $F(I_{i},:)$ is row full rank, there
exist a column index set $J_{i}\subset E_{i}$ such that $F(I_{i},J_{i})$ is a
nonsingular submatrix. Correspondingly, the matrix
$[F(I_{1},:);F(I_{2},:);\ldots;F(I_{s},:)]$ has a nonsingular submatrix as
follows,

$\left(\begin{array}[]{cccc}F(I_{1},J_{1})&&&\\\ &F(I_{2},J_{2})&&\\\
&&\ddots&\\\ &&&F(I_{s},J_{s})\end{array}\right).$

As a result, the matrix $[F(I_{1},:);F(I_{2},:);\ldots;F(I_{s},:)]$ is row
full rank and so is $F$.

To prove that $0$ is a regular value of
$\tilde{H}(\varphi,\lambda,t,K):\mathbb{R}^{mn}\times
R\times[0,1)\times\mathbb{R}^{(2nm-m)}\setminus(U_{1}\cup
U_{2})\to\mathbb{R}^{mn+1}$, we need to prove
$\forall(\varphi,\lambda,t,K)\in\mathbb{R}^{mn}\times
R\times[0,1)\times\mathbb{R}^{(2nm-m)}\setminus(U_{1}\cup U_{2})$ satisfying
$\tilde{H}(\varphi,\lambda,t,K)=0$, the Jacobian matrix of
$\tilde{H}(\varphi,\lambda,t,K)$ is row full rank.
$\forall(\varphi,\lambda,t,K)$ satisfying $\tilde{H}(\varphi,\lambda,t,K)=0$,
for $B$ defined in (59), we have $B(R^{0},:)=0$. Correspondingly for the
Jacobian matrix defined in (58), we have

$\displaystyle((1-t)A+D+3t\beta\mbox{diag}(\varphi)^{2}-\lambda
I)(R^{0},:)=((1-t)A+D-\lambda I)(R^{0},:).$

Through row permutations, the Jacobian matrix of
$\tilde{H}(\varphi,\lambda,t,K)$,
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$, can be rewritten as

$\left(\begin{array}[]{cccc}((1-t)A+D+3t\beta\mbox{diag}(\varphi)^{2}-\lambda
I)(R^{*},:)&-\varphi(R^{*})&(\beta\varphi^{3}-A\varphi)(R^{*})&(1-t)B(R^{*},:)\\\
((1-t)A+D-\lambda
I)(R^{0},:)&-\varphi(R^{0})&(\beta\varphi^{3}-A\varphi)(R^{0})&0\\\
-\varphi^{\mathrm{T}}&0&0&0\end{array}\right).$

$\forall(\varphi,\lambda,t,K)$ satisfying $\tilde{H}(\varphi,\lambda,t,K)=0$,
from Lemma 3.3, we know $B(R^{*},:)$ is row full rank. Therefore, if we can
prove

$\left(\begin{array}[]{c}((1-t)A+D-\lambda I)(R^{0},:)\\\
-\varphi^{\mathrm{T}}\end{array}\right)$ (68)

is row full rank, then
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$ is row full rank.
From $\tilde{H}=0$, we have $((1-t)A+D-\lambda I)(R^{0},:)\varphi=0$, i.e.,
$\varphi$ is orthogonal to the rows of $((1-t)A+D-\lambda I)(R^{0},:)$.
Therefore, if we can prove $((1-t)A+D-\lambda I)(R^{0},:)$ is row full rank,
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,K)}$ is row full rank.
Now the problem is turned into proving that $((1-t)A+D-\lambda I)(R^{0},:)$ is
row full rank.

For easy exposition, some concepts concerning the topology of the grid points
with zero function value are introduced. In addition, $(i,j)$ will be
considered as grid point in the rest of this section, representing
$(x_{i},y_{j})$.

###### Definition 3.6

In the grid $G_{a}$, a zero valued node $(i,j)$ is a grid point with
$\varphi_{i,j}=0$.

###### Definition 3.7

Two zero valued nodes $(i_{1},j_{1})$ and $(i_{p},j_{p})$ are said to be zero
valued connected, if there exist a sequence of zero valued nodes,

$\displaystyle(i_{1},j_{1}),(i_{2},j_{2}),\ldots,(i_{p},j_{p}),$ (69)

such that for any two successive nodes $(i_{k},j_{k})$ and $(i_{k+1},j_{k+1})$
of the sequence, the distance of these two nodes is 1 in the sense that

$\displaystyle|i_{k}-i_{k+1}|+|j_{k}-j_{k+1}|=1.$ (70)

###### Definition 3.8

A set $S$ consisting of zero valued nodes is called a zero valued connected
set if any two nodes of $S$ are zero valued connected.

###### Definition 3.9

A set $S$ consisting of zero valued nodes is called a zero valued connected
component if $S$ is connected and S is the largest zero connected set
containing $S$.

From (56), the discretization of the differential equation in a stencil is
written explicitly,

$\displaystyle\alpha_{1}^{(ij)}\varphi_{ij}+\alpha_{2}^{(ij)}\varphi_{i-1,j}+\alpha_{3}^{(ij)}\varphi_{i+1,j}+\alpha_{4}^{(ij)}\varphi_{i,j+1}+\alpha_{5}^{(ij)}\varphi_{i,j-1}=0,$
(71) $\displaystyle\varphi_{0j}=\varphi_{m+1,j}=0,$ (72)
$\displaystyle\varphi_{i0}=\varphi_{i,n+1}=0,$ (73)
$\displaystyle\varphi^{\mathrm{T}}\varphi-\frac{1}{h_{1}h_{2}}=0,$ (74)

where
$\alpha_{1}^{(ij)}=(1-t)a_{jj}^{(i)}+\frac{1}{h_{1}^{2}}+\frac{1}{h_{2}^{2}}+v_{ij}+t\beta\varphi_{ij}^{2}-\lambda$,
$\alpha_{2}^{(ij)}=\alpha_{3}^{(ij)}=(1-t)a_{jj+1}^{(i)}-\frac{1}{2h_{2}^{2}}$,
$\alpha_{4}^{(ij)}=\alpha_{5}^{(ij)}=-\frac{1}{2h_{1}^{2}}$, and
$i=1,\ldots,m$, $j=1,\ldots,n$. When $K\in R^{(2nm-m)}\setminus(U_{1}\bigcup
U_{2})$, the sign relationships among the components of $\varphi$ are given in
the following remark.

###### Remark 3.10

Let $(i,j)$ be an inner zero valued node. Assume that all except two of its
neighbouring points are known to be zero valued nodes. If one of the rest two
points is a zero valued node, then so is the other; If one of the rest two
points is not a zero valued node, then neither is the other.

For any set $G$ of zero valued nodes, denote by $R_{G}$ the set of indices of
rows corresponding to $G$, in which the matrix $B$ is zero, i.e.,

$\displaystyle R_{G}=\\{r:B(r,:)=0,r=j+(i-1)*n,(i,j)\in G\\}.$ (75)

Note that if a zero valued node $(i,j)\in G$ is such that $B(s,:)=0$ with
$s=j+(i-1)*n$, then the neighbouring inner grid points in $y$-direction should
be zero valued nodes, that is, both $(i,j-1)$ and $(i,j+1)$ are zero valued
nodes if $1<j<n$, or $(i,j+1)$ is a zero valued node if $j=1$, or $(i,j-1)$ is
a zero valued node if $j=n$. If $(i,j)\in G$ and the upper point $(i,j+1)$ or
the lower point $(i,j-1)$ is not a zero valued node, then the corresponding
row index of $s=\Gamma(i,j)$ will not be in $R_{G}$. Therefore
$\Gamma^{-1}(R_{G})\subset G$. Note that $R_{G}$ may be empty even if $G$ is
not empty.

Let $M=(1-t)A+D-\lambda I$. Denote by $C_{G}$ the indices of columns
corresponding to the zero valued connected set $G$, in which $M(R_{G},:)$ is
not zero, i.e.,

$\displaystyle C_{G}=\\{c\in R_{a}:M(R_{G},c)\neq 0\\}.$ (76)

If $R_{G}=\emptyset$, define $C_{G}=\emptyset$. For any $s\in R_{G}$, let
$(i,j)=\Gamma^{-1}(s)$. Then the $\Gamma$ images of $(i,j)$ and its
neighbouring inner grid points are possibly included in $C_{G}$.

Denote by $g_{i}$ the $i$-th column of grid points in $x$-direction, i.e.,

$\displaystyle g_{i}=\\{(i,j):1\leq j\leq n\\},$ (77)

and by $O_{1}$ ($O_{m}$) the ordering of all the inner grid points except the
first (last) column in $y$-direction,

$\displaystyle O_{1}=\\{r=j+(i-1)*n:(i,j)\in G_{a}\backslash g_{1}\\},$ (78)
$\displaystyle O_{m}=\\{r=j+(i-1)*n:(i,j)\in G_{a}\backslash g_{m}\\}.$ (79)

###### Lemma 3.11

Both $M(O_{1},:)$ and $M(O_{m},:)$ are row full rank.

Proof. It is obvious that $M(O_{1},:)$ has the following form:

$\left(\begin{array}[]{ccccc}*&\ldots&\ldots&\ldots&\ldots\\\
&*&\ldots&\ldots&\ldots\\\ &&\ddots&\ldots&\ldots\\\
&&&*&\ldots\end{array}\right),$

where $*$ represents nonzero element. Therefore $M(O_{1},:)$ is row full rank.

$M(O_{m},:)$ has the following form:

$\left(\begin{array}[]{ccccc}\ldots&*&&&\\\ \ldots&\ldots&*&&\\\
\ldots&\ldots&\ldots&\ddots&\\\
\ldots&\ldots&\ldots&\ldots&*\end{array}\right).$

Therefore $M(O_{m},:)$ is row full rank.

Note that if $G^{0}\bigcap g_{1}=\emptyset$, then $R^{0}\subset O_{1}$, from
Lemma 3.11, $M(R^{0},:)$ is row full rank. In the following, we consider the
case $G^{0}\bigcap g_{1}\not=\emptyset$. The set $G^{0}\bigcap g_{1}$ can have
its own zero valued connected components.

###### Lemma 3.12

For $m\geq n\geq 6$, suppose that there is a zero valued connected component
$\gamma$ of $G^{0}\bigcap g_{1}$ with $s$ points, $s\geq 2$. Then

  1. (i)

$s<n$;

  2. (ii)

If $\gamma$ contains the point $(1,1)$ or contains the point $(1,n)$, then the
zero valued connected component $G$ of $G^{0}$ containing $\gamma$ will be the
set of zero valued nodes starting from $\gamma$ and ending on column $s$ with
$s+1-i$ zero valued nodes on column $i$, $1\leq i\leq s$;

  3. (iii)

Suppose that $\gamma$ is located in the inner part of $g_{1}$, i.e., the grid
point $(1,j)$ of $\gamma$ is such that $2\leq j\leq n-1$. If $s$ is even, then
the zero valued connected component $G$ of $G^{0}$ containing $\gamma$ will be
the set of zero valued nodes starting from $\gamma$ and ending on column $s/2$
with $s+2-2i$ zero valued nodes on column $i$, $1\leq i\leq s/2$;

  4. (iv)

Suppose that $\gamma$ is located in the inner part of $g_{1}$. If $s$ is odd,
then there is a zero valued connected set $G$ of $G^{0}$ containing $\gamma$
and arriving at a single zero valued node on column $(s+1)/2$ with $s+2-2i$
zero valued nodes on column $i$, $1\leq i\leq(s+1)/2$;

  5. (v)

For the cases (2), (3), (4), if $s\geq 3$, the index set of nonzero columns of
$M$ corresponding to $G$ satisfies $C_{G}=\Gamma(G)$.

Proof.

  1. (i)

If $s=n$, we have $\varphi_{11}=\varphi_{12}=\ldots=\varphi_{1n}=0$. From the
sign relationships among the components of $\varphi$ as stated in Remark 3.10,
the grid points in the column 2 are all zero valued nodes, namely,
$\varphi_{21}=\varphi_{22}=\ldots=\varphi_{2n}=0$. By induction, all the grid
points are zero valued nodes, namely,
$\varphi_{11}=\ldots=\varphi_{1n}=\ldots=\varphi_{mn}=0$, a contradiction with
the condition that $\varphi^{T}\varphi\neq 0$. Therefore $s<n$.

  2. (ii)

If $\gamma$ starts from the point $(1,1)$, from the sign relationships among
the components of $\varphi$, the zero valued nodes connecting $\gamma$ on the
column 2 of grid points are $(2,j)$, $1\leq j\leq s-1$. By induction, it can
be seen that the zero valued connected component $G$ of $G^{0}$ containing
$\gamma$ will end on column $s$. The zero valued nodes connecting $\gamma$ are
illustrated in Fig 1 and form a zero valued component, where black dot
represents zero and white dot represents nonzero. Similarly, the case $\gamma$
ends at the point $(1,n)$ can be proved.

Figure 1: Illustration of a
flag$(1,1)$$(1,2)$$(1,s-1)$$(1,s)$$(1,s+1)$$(1,n)$$(2,1)$$(s-1,1)$$(s,1)$$(s+1,1)$$(m,1)$

  3. (iii)

Assume that $\gamma$ starts from the point $(1,j_{1})$ and ends at the point
$(1,j_{2})$ with $j_{1}\geq 2$ and $j_{2}\leq n-1$. From the sign
relationships among the components of $\varphi$, the zero valued nodes
connecting $\gamma$ on the column 2 of grid points are $(2,j)$, $3\leq j\leq
s-2$. The zero valued nodes connecting $\gamma$ are illustrated in Fig 2. The
zero valued nodes connecting $\gamma$ will end on column $s/2$ with two zero
valued nodes and form a zero valued component $G$.

Figure 2: Illustration of a
flag$(1,1)$$(1,t-1)$$(1,t)$$(1,t+1)$$(1,t+\frac{s}{2}-2)$$(1,t+\frac{s}{2}-1)$$(1,t+\frac{s}{2})$$(1,t+\frac{s}{2}+1)$$(1,t+s-2)$$(1,t+s-1)$$(1,t+s)$$(1,n)$$(2,1)$$(\frac{s}{2}-1,1)$$(\frac{s}{2},1)$$(m,1)$

  4. (iv)

Similar to the case (3), the zero valued nodes connecting $\gamma$ are
illustrated in Fig 3. However, since $s$ is odd, the zero valued nodes
connecting $\gamma$ arrives at one single point on column $(s+1)/2$. All these
zero valued nodes connecting $\gamma$ from column 1 to column $(s+1)/2$ form a
zero valued connected set $G$. It may be or may not be a zero valued connected
component.

Figure 3: Illustration of a
flag$(1,10)$$(1,t-1)$$(1,t)$$(1,t+1)$$(1,t+\frac{s-3}{2})$$(1,t+\frac{s-1}{2})$$(1,t+\frac{s+1}{2})$$(1,t+s-2)$$(1,t+s-1)$$(1,t+s)$$(1,n)$$(2,1)$$(\frac{s-1}{2},1)$$(\frac{s+1}{2},1)$$(m,1)$

  5. (v)

From the figures for cases (2), (3), (4), it can be seen that
$\Gamma^{-1}(R_{G})$ is the set of zero valued nodes by shrinking $G$ one
layer from its nonzero boundary. $\Gamma^{-1}(C_{G})$ is the set of zero
valued nodes by extending $\Gamma^{-1}(R_{G})$ one layer towards the nonzero
boundary, which is $G$ itself.

###### Definition 3.13

Let $\gamma$ be a zero valued connected component of $G^{0}\bigcap g_{1}$ with
$s$ points, $s\geq 2$. The zero valued connected set $G$ mentioned in (2),
(3), (4) in Lemma 3.12 is called a flag of $\gamma$, denoted as $F_{\gamma}$.

###### Lemma 3.14

Let $\gamma$ be a zero valued connected component of $G^{0}\bigcap g_{1}$ of
$s$ points with $s>2$ an odd number and $F_{\gamma}$ be its flag. Let $G$ be
another set of zero valued nodes. Let $I_{1}=R_{F_{\gamma}}$ and $I_{2}=R_{G}$
be the row index sets as defined in (75) such that $M(I_{2},:)$ is row full
rank. Denote

$\displaystyle J_{1}=\\{j\in R_{a}:M(I_{1},j)\neq 0\\},\quad J_{2}=\\{j\in
R_{a}:M(I_{2},j)\neq 0\\},$ (80) $\displaystyle J_{12}=J_{1}\bigcap J_{2}.$
(81)

Suppose that $J_{12}$ has only one element $k$ and that both the column
$M(I_{1},k)$ and the column $M(I_{2},k)$ have only one nonzero element,
denoted as $M(s_{1},k)\neq 0$ and $M(s_{2},k)\neq 0$ respectively. Let
$(i_{1},j_{1})$ and $(i_{2},j_{2})$ be the grid point corresponding to $s_{1}$
and $s_{2}$ respectively. Suppose $i_{2}=i_{1}+2$ and $j_{1}=j_{2}$. Then
$[M(I_{1},:);M(I_{2},:)]$ is row full rank.

Proof. Note that since $J_{12}$ is not empty, the zero valued connected
component $\gamma$ corresponds to the case (4) of Lemma 3.12. Geometrically,
for that $k$, the two grid points $(i_{1},j_{1})$ and $(i_{2},j_{2})$
corresponding to $M(s_{1},k)\neq 0$ and $M(s_{2},k)\neq 0$ lie in the same row
of grid points with distance 2 in the sense that
$|i_{2}-i_{1}|+|j_{2}-j_{1}|=2$, as illustrated locally in Figure 4.

Figure 4: Local illustration of connection with a
flag$(i_{1},j_{1})$$(i_{2},j_{2})$

Note that by Lemma 3.11, $M(I_{1},:)$ is row full rank. By row permutations
and column permutations, $[M(I_{1},:);M(I_{2},:)]$ is transformed to the
following form, denoted as $W$,

$\begin{array}[]{l@{\hspace{-5pt}}llll}\begin{array}[]{@{\hspace{10pt}}l@{\hspace{4pt}}l@{\hspace{2pt}}l@{\hspace{18pt}}l}\hskip
10.0pt\lx@intercol\overbrace{\hphantom{\begin{array}[]{ccccccc}\ast&\ast&\ast&\ddots&\ddots&\ddots&\ast\end{array}}}^{\displaystyle
s}\hfil\hskip
4.0&\overbrace{\hphantom{\begin{array}[]{lcccc}\ast&\ddots&\ddots&\ast&\ast\end{array}}}^{\displaystyle
s-2}\hfil\hskip
2.0&\overbrace{\hphantom{\begin{array}[]{lcc}\ast&\ddots&\ast\end{array}}}^{\displaystyle
s-4}\hfil\hskip
18.0&\overbrace{\hphantom{\begin{array}[]{lcc}\ast&\ast&\ast\end{array}}}^{\displaystyle
3}\par\end{array}&\\\
\left(\begin{array}[]{ccccccc|ccccc|ccc|c|ccc|c|cc}\ast&\ast&\ast&&&&&\ast&&&&&&&&&&&&&&\\\
&\ast&\ast&\ast&&&&&\ast&&&&&&&&&&&&&\\\
&&\ast&\ast&\ast&&&&&\ast&&&&&&&&&&&&\\\
&&&\ddots&\ddots&\ddots&&&&&\ddots&&&&&&&&&&&\\\
&&&&\ast&\ast&\ast&&&&&\ast&&&&&&&&&\\\
&&\ast&&&&&\ast&\ast&\ast&&&\ast&&&&&&&&&\\\
&&&\ddots&&&&&\ddots&\ddots&\ddots&&&\ddots&&&&&&&&\\\
&&&&\ast&&&&&&\ast&\ast&&&\ast&&&&&&&\\\
&&&&&&&&&\ddots&&&&\ddots&&\ddots&&&&&&\\\ &&&&&&&&&&&&&&&&\ast&&&&&\\\
&&&&&&&&&&&&&&&&&\ast&&&&\\\ &&&&&&&&&&&&&&&&&&\ast&&&\\\
&&&&&&&&&&&&&&&\cdots&\ast&\ast&\ast&\ast&&\\\
\hline\cr&&&&&&&&&&&&&&&&&&&\ast&\ast&\cdots\\\ &&&&&&&&&&&&&&&&&&&&\ast&\\\
&&&&&&&&&&&&&&&&&&&&&\ddots\end{array}\right)&\begin{array}[]{l}\vspace{35pt}\left.\vphantom{\begin{array}[]{c}\ast\\\
\ast\\\ \ast\\\ \ddots\\\ \ast\\\ \ast\\\ \ddots\\\ \ast\\\ \ddots\\\ \ast\\\
\ast\\\ \ast\\\
\ast\end{array}}\right\\}\frac{(s-1)^{2}}{4}\end{array}\end{array}.$

Note that the diagonal entries of the submatrix
$W(1:\frac{(s-1)^{2}}{4},s+1:\frac{(s+1)^{2}}{4})$ are all nonzero. Therefore,
by elementary matrix transformations on the first $\frac{(s+1)^{2}}{4}$
columns, the above matrix $W$ is transformed to the following, denoted as
$\overline{W}$,

$\begin{array}[]{l@{\hspace{-5pt}}llll}\begin{array}[]{@{\hspace{1pt}}l@{\hspace{2pt}}l@{\hspace{2pt}}l@{\hspace{15pt}}l}\hskip
1.0pt\lx@intercol\overbrace{\hphantom{\begin{array}[]{ccccccc}\ast&\ast&\ast&\ddots&\ddots&\ddots&\ast\end{array}}}^{\displaystyle
s}\hfil\hskip
2.0&\overbrace{\hphantom{\begin{array}[]{lcccc}\ast&\ddots&\ddots&\ast&\ast\end{array}}}^{\displaystyle
s-2}\hfil\hskip
2.0&\overbrace{\hphantom{\begin{array}[]{lcc}\ast&\ddots&\ast\end{array}}}^{\displaystyle
s-4}\hfil\hskip
15.0&\overbrace{\hphantom{\begin{array}[]{lcc}\ast&\ast&\ast\end{array}}}^{\displaystyle
3}\par\end{array}&\\\
\left(\begin{array}[]{ccccccc|ccccc|ccc|c|ccc|c|cc}&&&&&&&\ast&&&&&&&&&&&&&&\\\
&&&&&&&&\ast&&&&&&&&&&&&&\\\ &&&&&&&&&\ast&&&&&&&&&&&&\\\
&&&&&&&&&&\ddots&&&&&&&&&&&\\\ &&&&&&&&&&&\ast&&&&&&&&&\\\
&&&&&&&&&&&&\ast&&&&&&&&&\\\ &&&&&&&&&&&&&\ddots&&&&&&&&\\\
&&&&&&&&&&&&&&\ast&&&&&&&\\\ &&&&&&&&&&&&&&&\ddots&&&&&&\\\
&&&&&&&&&&&&&&&&\ast&&&&&\\\ &&&&&&&&&&&&&&&&&\ast&&&&\\\
&&&&&&&&&&&&&&&&&&\ast&&&\\\ &&&&&&&&&&&&&&&&&&&\ast&&\\\
\hline\cr\times&\ast&\ast&\cdots&\ast&\ast&\ast&\ast&\ast&\cdots&\ast&\ast&\ast&\cdots&\ast&\cdots&\ast&\ast&\ast&\ast&\ast&\cdots\\\
&&&&&&&&&&&&&&&&&&&&\ast&\\\
&&&&&&&&&&&&&&&&&&&&&\ddots\end{array}\right)&\begin{array}[]{l}\vspace{35pt}\left.\vphantom{\begin{array}[]{c}\ast\\\
\ast\\\ \ast\\\ \ddots\\\ \ast\\\ \ast\\\ \ddots\\\ \ast\\\ \ddots\\\ \ast\\\
\ast\\\ \ast\\\
\ast\end{array}}\right\\}\frac{(s-1)^{2}}{4}\end{array}\end{array}.$

By elementary matrix transformations on the first $\frac{(s-1)^{2}}{4}+1$
rows, the above matrix $\overline{W}$ is further transformed to the following,
denoted as $\overline{\overline{W}}$,

$\begin{array}[]{l@{\hspace{-5pt}}llll}\begin{array}[]{@{\hspace{0.1pt}}l@{\hspace{1pt}}l@{\hspace{1pt}}l@{\hspace{9pt}}l}\hskip
0.1pt\lx@intercol\overbrace{\hphantom{\begin{array}[]{ccccccc}\ast&\ast&\ast&\ddots&\ddots&\ddots&\ast\end{array}}}^{\displaystyle
s}\hfil\hskip
1.0&\overbrace{\hphantom{\begin{array}[]{lcccc}\ast&\ddots&\ddots&\ast&\ast\end{array}}}^{\displaystyle
s-2}\hfil\hskip
1.0&\overbrace{\hphantom{\begin{array}[]{lcc}\ast&\ddots&\ast\end{array}}}^{\displaystyle
s-4}\hfil\hskip
9.0&\overbrace{\hphantom{\begin{array}[]{lcc}\ast&\ast&\ast\end{array}}}^{\displaystyle
3}\par\end{array}&\\\
\left(\begin{array}[]{ccccccc|ccccc|ccc|c|ccc|c|cc}&&&&&&&\ast&&&&&&&&&&&&&&\\\
&&&&&&&&\ast&&&&&&&&&&&&&\\\ &&&&&&&&&\ast&&&&&&&&&&&&\\\
&&&&&&&&&&\ddots&&&&&&&&&&&\\\ &&&&&&&&&&&\ast&&&&&&&&&\\\
&&&&&&&&&&&&\ast&&&&&&&&&\\\ &&&&&&&&&&&&&\ddots&&&&&&&&\\\
&&&&&&&&&&&&&&\ast&&&&&&&\\\ &&&&&&&&&&&&&&&\ddots&&&&&&\\\
&&&&&&&&&&&&&&&&\ast&&&&&\\\ &&&&&&&&&&&&&&&&&\ast&&&&\\\
&&&&&&&&&&&&&&&&&&\ast&&&\\\ &&&&&&&&&&&&&&&&&&&\ast&&\\\
\hline\cr\times&\ast&\ast&\cdots&\ast&\ast&\ast&&&&&&&&&&&&&&\ast&\cdots\\\
&&&&&&&&&&&&&&&&&&&&\ast&\\\
&&&&&&&&&&&&&&&&&&&&&\ddots\end{array}\right)&\begin{array}[]{l}\vspace{35pt}\left.\vphantom{\begin{array}[]{c}\ast\\\
\ast\\\ \ast\\\ \ddots\\\ \ast\\\ \ast\\\ \ddots\\\ \ast\\\ \ddots\\\ \ast\\\
\ast\\\ \ast\\\
\ast\end{array}}\right\\}\frac{(s-1)^{2}}{4}\end{array}\end{array}.$

Note that in the row $\frac{(s-1)^{2}}{4}+1$, at least one element
$\overline{\overline{W}}(\frac{(s-1)^{2}}{4}+1,1)$, denoted as
${}^{\prime}\times^{\prime}$, is not zero. Therefore, the lower submatrix
$\overline{\overline{W}}(\frac{(s-1)^{2}}{4}+1:end,:)$ is still row full rank.
Now Lemma 3.5 can be applied to conclude that $\overline{\overline{W}}$ is row
full rank. Therefore $[M(I_{1},:);M(I_{2},:)]$ is row full rank.

###### Theorem 3.15

For the homotopy
$H:\mathbb{R}^{mn}\times\mathbb{R}\times[0,1)\rightarrow\mathbb{R}^{mn}\times\mathbb{R}$
in (56), $\forall n\in\mathbb{N}^{+}$, $m\geq n\geq 6$, for almost all
$K\in\mathbb{R}^{(2nm-m)}\setminus(U_{1}\cup U_{2})$,

  1. (i)

0 is a regular value of $H$ defined in $(\ref{sec3:eqn-homotopy-3diagonal})$
and therefore the homotopy paths corresponding to different initial points do
not intersect each other for $t\in[0,1)$;

  2. (ii)

Every homotopy path $(\varphi(s),\lambda(s),t(s))\subset H^{-1}(0)$ is
bounded.

Proof. (i). It suffices to prove that
$\forall(\varphi,\lambda,t,K)\in\mathbb{R}^{mn}\times\mathbb{R}\times[0,1)\times\mathbb{R}^{(2nm-m)}\setminus(U_{1}\cup
U_{2})$ satisfying $\tilde{H}(\varphi,\lambda,t,K)=0$, $((1-t)A+D-\lambda
I)(R^{0},:)$ is row full rank, or in short hand notation $M(R^{0},:)$ is row
full rank. Note that $R^{0}$ corresponds to the zero valued nodes $G^{0}$. If
$G^{0}\bigcap g_{1}=\emptyset$, then $R^{0}\subset O_{1}$. By Lemma 3.11
$M(R^{0},:)$ is row full rank. Next $G^{0}\bigcap g_{1}\neq\emptyset$ is
assumed. If the zero valued connected components of the set $G^{0}\bigcap
g_{1}$ are all single point sets, $R^{0}$ is a subset of $O_{1}$, and again by
Lemma 3.11, $M(R^{0},:)$ is row full rank.

Assume that there are $q$ zero valued connected components of the set
$G^{0}\bigcap g_{1}$, each of which has more than one point. These connected
components of the set $G^{0}\bigcap g_{1}$ are denoted as $\gamma_{i}$, with
corresponding flags $G_{i}=F_{\gamma_{i}}$, $i=1,\cdots,q$. Denote

$\displaystyle\overline{G}=G_{a}\setminus(\bigcup_{i=1}^{q}G_{i}),$ (82)
$\displaystyle\overline{R}=R_{\overline{G}},\quad\overline{C}=C_{\overline{G}},\quad
R_{i}=R_{G_{i}},\quad C_{i}=C_{G_{i}},\quad\quad\forall 1\leq i\leq q.$ (83)

By Lemma 3.12, it can be verified that

$\displaystyle R^{0}=\overline{R}\bigcup\left(\bigcup_{i=1}^{q}R_{i}\right),$
(84) $\displaystyle R_{i}\bigcap\overline{R}=\emptyset,\quad\forall 1\leq
i\leq q,$ (85) $\displaystyle R_{i}\bigcap R_{j}=\emptyset,\quad C_{i}\bigcap
C_{j}=\emptyset,\quad\forall i\neq j.$ (86)

Note that if a zero valued connected component $\gamma$ of $G^{0}\bigcap
g_{1}$ has only two nodes, then the flag of $\gamma$ is the $\gamma$ itself,
and $R_{\gamma}=\emptyset$. Assume that all the zero valued connected
components $\gamma_{i}$ have more than two nodes. Since

$\displaystyle\overline{R}\subset O_{1},\quad R_{i}\subset O_{m},\quad\forall
1\leq i\leq q,$ (87)

by Lemma 3.11, all of $M(\overline{R},:)$ and $M(R_{i},:)$, $1\leq i\leq q$,
are row full rank. It is possible that for some $\gamma_{i}$,
$C_{i}\bigcap\overline{C}\neq\emptyset$. If so, $\gamma_{i}$ is the case
described in case $(4)$ in Lemma 3.12, $C_{i}\bigcap\overline{C}$ has only one
element and the corresponding grid points are illustrated in Figure 5.

Figure 5: Local illustration of a connection$G_{i}$$\overline{G}$

Without loss of generality, suppose that all the $\gamma_{i}$ satisfying such
property are the first $p$ $\gamma_{i}$, i.e.,

$\displaystyle C_{i}\bigcap\overline{C}\neq\emptyset,\quad 1\leq i\leq p,$
(88) $\displaystyle C_{i}\bigcap\overline{C}=\emptyset,\quad p+1\leq i\leq q.$
(89)

Denote

$\displaystyle G$
$\displaystyle=\overline{G}\bigcup\left(\bigcup_{i=p+1}^{q}G_{i}\right),$ (90)
$\displaystyle Z(G_{k},\ldots,G_{1},G)$
$\displaystyle=[M(R_{k},:);\ldots;M(R_{1},:);M(R_{G},:)],\quad\forall 1\leq
k\leq p.$ (91)

Then $M(R^{0},:)=Z(G_{p},\ldots,G_{1},G)$. The claim that $M(R^{0},:)$ is row
full rank will be proved by recursion.

Firstly, prove $Z(G_{1},G)$ is row full rank. Take $I_{1}=R_{1}$ and
$I_{2}=R_{G}$. The conditions of Lemma 3.14 are satisfied. Thus
$[M(R_{1},:);M(R_{G},:)]$ is row full rank.

Secondly, prove $Z(G_{k+1},\ldots,G_{1},G)$ is row full rank if
$Z(G_{k},\ldots,G_{1},G)$ is, $\forall 1\leq k<p$. Take $I_{1}=R_{k+1}$ and
$I_{2}=R_{k}\cup\cdots\cup R_{1}\cup R_{G}$. Since $C_{k+1}\cap
C_{i}=\emptyset$, $1\leq i\leq k$, and $C_{k+1}\bigcap C_{G}$ has only one
element, $I_{1}$ and $I_{2}$ satisfy the conditions of Lemma 3.14. Thus
$[M(R_{k},:);\ldots;M(R_{1},:);M(R_{G},:)]$ is row full rank.

(ii). Similar to the one dimensional case, that $H^{-1}(0)$ is bounded can be
proved also.

####  3.2.2 The homotopy with random pentadiagonal matrix

A homotopy with random pentadiagonal matrix is also possible. Specifically,
for $\tilde{H}$ defined in (57), replacing $A(K)$ with
$\overline{A}(\overline{K})$, where
$\overline{A}(\overline{K})\in\mathbb{R}^{mn\times mn}$ is a random
pentadiagonal matrix with the same sparse structure as $D$, namely,

$\overline{A}(\overline{K})=\left(\begin{array}[]{cccc}A_{1}&\overline{A}_{1}&&\\\
\overline{A}_{1}&A_{2}&\overline{A}_{2}&\\\ &&\ddots&\overline{A}_{m-1}\\\
&&\overline{A}_{m-1}&A_{m}\end{array}\right),$

where

$A_{i}=\left(\begin{array}[]{ccccc}a_{11}^{(i)}&a_{12}^{(i)}&&&\\\
a_{12}^{(i)}&a_{22}^{(i)}&a_{23}^{(i)}&&\\\
&a_{23}^{(i)}&a_{33}^{(i)}&\ddots&\\\ &&\ddots&\ddots&a_{n-1,n}^{(i)}\\\
&&&a_{n-1,n}^{(i)}&a_{nn}^{(i)}\end{array}\right),\quad\overline{A}_{i}=\left(\begin{array}[]{ccccc}b_{1}^{(i)}&&&&\\\
&b_{2}^{(i)}&&&\\\ &&b_{3}^{(i)}&&\\\ &&&\ddots&\\\
&&&&b_{n}^{(i)}\end{array}\right),$

$\displaystyle\overline{K}=\left(K_{1},K_{2},\ldots,K_{m},\overline{K}_{1},\ldots,\overline{K}_{m-1}\right)^{\mathrm{T}},$
$\displaystyle
K_{i}=\left(a_{11}^{(i)},a_{12}^{(i)},a_{22}^{(i)},a_{23}^{(i)},\ldots,a_{n-1,n-1}^{(i)},a_{n-1,n}^{(i)},a_{nn}^{(i)}\right),~{}i=1,\ldots,m,$
$\displaystyle\overline{K}_{i}=\left(b_{1}^{(i)},b_{2}^{(i)},\ldots,b_{n}^{(i)}\right),~{}i=1,\ldots,m-1.$

The Jacobian matrix of $\tilde{H}(\varphi,\lambda,t,\overline{K})$,
$\frac{\partial\tilde{H}}{\partial(\varphi,\lambda,t,\overline{K})}$, is :

$\left(\begin{array}[]{cccc}(1-t)\overline{A}(\overline{K})+D+3t\beta\mbox{diag}(\varphi^{2})-\lambda
I&-\varphi&\beta\varphi^{3}-\overline{A}(\overline{K})\varphi&(1-t)\overline{B}\\\
-\varphi^{\mathrm{T}}&0&0&0\end{array}\right),$

where
$\overline{B}=\frac{\partial(\overline{A}(\overline{K})\varphi)}{\partial\overline{K}}\in\mathbb{R}^{{(mn)}\times{(3nm-
m-n)}}$. Recall
$\varphi_{i}=\left(\varphi_{i1},\varphi_{i2},\ldots,\varphi_{in}\right)^{\mathrm{T}}$.
It can be verified that

$\displaystyle\overline{B}$
$\displaystyle=\left(\begin{array}[]{ccccc|cccc}B_{1}&&&&&\overline{B}_{2}&&&\\\
&B_{2}&&&&\overline{B}_{1}&\overline{B}_{3}&&\\\
&&\ddots&&&&\overline{B}_{2}&\overline{B}_{4}&\\\
&&&\ddots&&&&\ddots&\ddots\\\
&&&&B_{m}&&&&\overline{B}_{m-1}\end{array}\right)$ (97)
$\displaystyle=\left(\begin{array}[]{ccccc|cccc}&&&&&\overline{B}_{2}&&&\\\
&&&&&\overline{B}_{1}&\overline{B}_{3}&&\\\
&&B&&&&\overline{B}_{2}&\overline{B}_{4}&\\\ &&&&&&&\ddots&\ddots\\\
&&&&&&&&\overline{B}_{m-1}\end{array}\right),$ (103)

with

$\displaystyle\overline{B}_{i}=\left(\begin{array}[]{cccc}\varphi_{i1}&&&\\\
&\varphi_{i2}&&\\\ &&\ddots&\\\ &&&\varphi_{in}\end{array}\right).$ (108)

Note that the left part $B$ of $\overline{B}$ is nothing but the matrix in
(59).

###### Lemma 3.16

The eigenvalues of $\overline{A}(\overline{K})+D$ are simple for
$\overline{K}$ almost everywhere except on a subset of real codimension 1.

Proof. Similar to the proof of Lemma 3.4.

For our discussion, $\overline{U}_{2}=(\mathbb{R}^{+})^{(3nm-m-n)}$ is
removed. To prove that
$\forall(\varphi,\lambda,t,\overline{K})\in\mathbb{R}^{mn}\times
R\times[0,1)\times\mathbb{R}^{(3nm-
m-n)}\setminus(\overline{U}_{1}\cup\overline{U}_{2})$ satisfying
$\tilde{H}(\varphi,\lambda,t,\overline{K})=0$, the Jacobian matrix of
$\tilde{H}(\varphi,\lambda,t,\overline{K})$ in (3.2.2) is row full rank, it
suffices to prove that the following submatrix of the Jacobian matrix in
(3.2.2)

$\left(\begin{array}[]{cccc}(1-t)\overline{A}(\overline{K})+D+3t\beta\mbox{diag}(\varphi^{2})-\lambda
I&-\varphi&\beta\varphi^{3}-\overline{A}(\overline{K})\varphi&(1-t)B\\\
-\varphi^{\mathrm{T}}&0&0&0\end{array}\right)$ (109)

is row full rank. With the notations $R^{0}$ and $R^{*}$ as in (65) and
similar arguments as in Subsection 3.2.1, it can be proved that
$((1-t)\overline{A}+D-\lambda I))(R^{0},:)$ is row full rank. Therefore the
matrix (109) is row full rank. That is $0$ is a regular value of the homotopy
$H$ with random pentadiagonal matrix for almost all
$\overline{K}\in\mathbb{R}^{(3nm-m-n)}\setminus(U_{1}\cup U_{2})$ with $m\geq
n\geq 6$.

##  4 Algorithm and numerical results

###  4.1 Algorithm

Thanks to Theorems 3.2 and 3.15, since 0 is a regular value of the homotopies
constructed, the homotopy paths determined by the homotopy equations (49) and
(56) have no bifurcation points with probability one. Therefore the usual path
following algorithm, i.e., the predictor-corrector method as in [23, 24], can
be adapted to trace the homotopy paths of the homotopy equations (49) and
(56). The adapted algorithm is stated in Algorithm 1. For notation
convenience, in Algorithm 1, we denote the iterate at step $k$ as
$x_{k}=(\varphi^{(k)},\lambda_{k})$, $k=1,\ldots$.

Initialization: Set $(x_{0},t_{0})=(x_{0},0)$, $k=0$, the minimum step size
$ds_{min}$, the initial step size $ds$. Compute the tangent vector
$(\dot{x}_{0},\dot{t}_{0})$ such that $\dot{t}_{0}>0$ and record the
orientation $ori$ ${H_{x}\dot{x}_{0}+H_{t}\dot{t}_{0}=0},\qquad
ori=\mbox{sign}\left(\left|\begin{array}[]{cc}H_{x}&H_{t}\\\
\dot{x}_{0}&\dot{t}_{0}\end{array}\right|\right).$ (110) while  _$t_{k} <1$_
do

Predictor:
$(\bar{x}_{k+1},\bar{t}_{k+1})=(x_{k},t_{k})+ds(\bar{x}_{k},\bar{t}_{k})$;

if  _$\bar{t}_{k+1} >1$_ then

change ds such that $\bar{t}_{k+1}=1$;

end if

Corrector: if  _$\bar{t}_{k+1}=1$_ then

$(v,\tau)=(0,1)$;

else

$(v,\tau)=(\dot{x}_{k},\dot{t}_{k})$;

end if

Employ Newton method to solve the following nonlinear equations:
$\begin{array}[]{ccc}\left(\begin{array}[]{c}H(x,t)\\\
v^{T}(x-\bar{x}_{k+1})+\tau(t-\bar{t}_{k+1})\end{array}\right)&=&0.\end{array}$
(111)

Judgement: if  _the above iteration converges to $(x_{k+1},t_{k+1})$_ then

compute the tangent vector $(\dot{x}_{k+1},\dot{t}_{k+1})$ satisfying
${H_{x}\dot{x}_{k+1}+H_{t}\dot{t}_{k+1}=0},\qquad\mbox{sign}\left(\left|\begin{array}[]{cc}H_{x}&H_{t}\\\
\dot{x}_{k+1}&\dot{t}_{k+1}\end{array}\right|\right)=ori$ (112)

else

$ds=\frac{1}{2}ds$;

go to Predictor;

end if

Compute angle $\theta$ of $(\bar{x}_{k},\bar{t}_{k})$ and
$(\bar{x}_{k+1},\bar{t}_{k+1})$;

if  _$\theta >18^{0}$_ then

$ds=\frac{1}{2}ds,$ go to Predictor;

end if

Accept iterates:
$(x_{k},t_{k})=(x_{k+1},t_{k+1})$,$(\dot{x}_{k},\dot{t}_{k})=(\dot{x}_{k+1},\dot{t}_{k+1})$;

if  _$\theta <6^{0}$_ then

$ds=2ds$;

end if

if  _$ds <ds_{min}$_ then

stop the algorithm;

end if

end while

Algorithm 1 Predictor-Corrector method

###  4.2 Numerical results

The first numerical example is the 1D discretized problem (9) with $\beta=20$,
$V(x)=\frac{1}{2}x^{2}$, $\Omega=[-2,2]$ and $n=999$. Some eigenvectors, i.e.,
approximate eigenfunctions are plotted in Figures 6(a) to 6(i). The
eigenvector in Figure 6(a) is the unique positive solution of (9) as stated in
Theorem 2.5, which corresponds to the unique positive ground state, as proved
in [4] for the continuous nonlinear eigenvalue problem (4)-(6) when $\beta>0$.
The approximate eigenfunction in Figure 6(b) is antisymmetric as described in
Theorem 2.6 and is an approximate first excited state. Others are approximate
excited states corresponding to higher energy. The order preserving property
of the eigenvalue curves as stated in [17] was observed, that is, if
$\lambda(0)$ is the $k$th smallest eigenvalue of the initial problem, then
$\lambda(t)$ is the $k$th smallest eigenvalue of the intermediate problem for
each $t\in[0,1)$. However, we are not able to prove such property for the
eigenvector-dependent nonlinear eigen-problem yet.

(a) $\lambda=6.76$

(b) $\lambda=8.39$

(c) $\lambda=10.35$

(d) $\lambda=12.73$

(e) $\lambda=15.62$

(f) $\lambda=19.08$

(g) $\lambda=23.13$

(h) $\lambda=27.79$

(i) $\lambda=33.05$

The second numerical example is the 2D discretized problem $(\ref{sec 2:
finite difference - two dimen})$ with $\beta=20$,
$V(x)=\frac{1}{2}(x_{1}^{2}+x_{2}^{2})$, $\Omega=[0,1]\times[0,1]$ and
$m=n=29$. Some approximate eigenfunctions are collected in Figures 6(j) to
6(u). The approximate eigenfunction in Figure 6(j) corresponds to the unique
positive ground state. Others are approximate excited states corresponding to
higher energy. The order preserving property of the eigenvalue curves is also
observed for the 2D case.

(j) $\lambda=43.36$

(k) $\lambda=60.89$

(l) $\lambda=65.31$

(m) $\lambda=78.81$

(n) $\lambda=86.28$

(o) $\lambda=94.06$

(p) $\lambda=104.4$

(q) $\lambda=108.82$

(r) $\lambda=120.47$

(s) $\lambda=129.67$

(t) $\lambda=133.81$

(u) $\lambda=138.68$

##  5 Conclusion

Solutions to the discretized problem with the finite difference disretization
for the GPE inherit certain properties of the solutions to the continuous
problem, such as the existence and uniqueness of positive eigenvector
(eigenfunction). The designed homotopy continuation methods are suitable for
computing eigenpairs corresponding to excited states of high energy as well as
the ground state and the first excited state. In order to make sure that the
homotopy paths are regular and that the path following is efficient,
artificial homotopy parameter and random matrices with certain structures in
the homotopies seem indispensable.

## References

  * [1] W. Bao and W. Tang, Ground state solution of Bose-Einstein condensate by directly minimizing the energy functional, J. Comput. Phys., 187 (2003), 230–254.
  * [2] S.-L. Chang, H.-S. Chen, B.-W. Jeng, and C.-S. Chien, A spectral-Galerkin continuation method for numerical solutions of the Gross-Pitaevskii equation, J. Comput. Appl. Math., 254 (2013), 2–16.
  * [3] W. Bao and Q. Du, Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, SIAM J. Sci. Comput., 25(5) (2003), 1674–1697.
  * [4] W. Bao and Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, J. Kinetic & Related Models., 6 (2013), 1–135.
  * [5] M. Edwards and K. Burnett, Numerical solution of the nonlinear schrödinger equation for small samples of trapped neutral atoms, Phys. Rev. A., 51 (1995), 1382–1386.
  * [6] S. K. Adhikari, Numerical solution of the two-dimensional Gross-Pitaevskii equation for trapped interacting atoms, Phys. Lett. A., 265 (2000), 91–96.
  * [7] S.-L. Chang and C.-S. Chien, Numerical continuation for nonlinear schrödinger equations, Int. J. Bifurcat. Chaos., 17 (2007), 641–656.
  * [8] S.-L. Chang, C.-S. Chien, and B.-W. Jeng, Computing wave functions of nonlinear Schrödinger equations: a time-independent approach, J. Comput. Phys., 226 (2007), 104–130.
  * [9] Y. S. Choi, I. Koltracht, P. J. McKenna, and N. Savytska, Global monotone convergence of Newton iteration for a nonlinear eigen-problem, Linear Algebra Appl., 357 (2002), 217–228.
  * [10] Y. S. Choi, I. Koltracht, and P. J. McKenna, A generalizaion of the Perron-Frobenius theorem for nonlinear perturbations of stiltjes matrices, Contemp. Math. AMS., 281 (2001), 325–330.
  * [11] Y. Cai, L.-H. Zhang, Z. Bai, and R.-C. Li, On an eigenvector-dependent nonlinear eigenvalue problem, SIAM J. Matrix Anal. Appl., 39 (2018), 1360–1382.
  * [12] Y. Saad, J. R. Chelikowsky, and S. M. Shontz, Numerical methods for electronic structure calculations of materials, SIAM Rev., 52 (2010), 3–54.
  * [13] E. Jarlebring, S. Kvaal, and W.Michiels, An inverse iteration method for eigenvalue problems with eigenvector nonlinearities, SIAM J. Sci. Comput., 36 (2014), 1978–2001.
  * [14] X. D. Yao and J. X. Zhou, Numerical methods for computing nonlinear eigenpairs: part I. Iso-homogeneous cases, SIAM J. Sci. Comput., 29 (2007), 135–1374.
  * [15] Q. Z. Yang, P. F. Huang and Y. J. Liu, Numerical examples for solving a class of nonlinear eigenvalue problems, Journal on Numerical Methods and Computer Applications, 40 (2019), 130–142.
  * [16] H. H. Xie, A multigrid method for nonlinear eigenvalue problems (in Chinese), Sci. Sin. Math., 45 (2015), 1193–1204.
  * [17] T. Y. Li and N. H. Rhee, Homotopy algorithm for symmetric eigenvalue problems, Numer. Math., 55 (1989), 265–280.
  * [18] T. Y. Li, H. Zhang, and X. H. Sun, Parallel homotopy algorithm for the symmetric tridiagonal eigenvalue problems, SIAM J. Sci. Stat. Comput., 12 (2006), 469–487.
  * [19] T. Y. Li, Z. Zeng, and L. Cong, Solving eigenvalue problems of real nonsymmetric matrices with real homotopies, SIAM J. Numer. Anal., 29 (1992), 229–248.
  * [20] S. H. Lui, H. B. Keller, and T. W. Kwok, Homotopy method for the large, sparse, real nonsymmetric eigenvalue problem, SIAM J. Matrix Anal. Appl. 18 (1997), 312–333.
  * [21] T. Y. Li And T. Sauer, Homotopy method for generalized eigenvalue problems Ax=$\lambda$Bx, Linear Algebra. Appl., 91 (1987), 65–74.
  * [22] M. T. Chu, T. Y. Li, and T. Sauer, Homotopy method for general $\lambda$-matrix problems, SIAM J. Matrix Anal. Appl., 9 (1988), 528–536.
  * [23] E. L. Allgower and K. Georg, Introduction to Numerical Continuation Methods, Classics in Applied Mathematics, 45, SIAM, Philadelphia, 2003.
  * [24] X. D. Huang, Z. G. Zeng and Y. N. Ma, Theories and Methods of Nonlinear Numerical Analysis (in Chinese), Wuhan University Press, Wuhan, 2004.