A positivity preserving iterative method for finding the ground states of saturable nonlinear Schr\"odinger equations

TL;DR

This paper introduces a globally convergent iterative method for computing positive ground states of saturable nonlinear Schrödinger equations, ensuring positivity and quadratic convergence.

Contribution

The paper presents a novel iterative approach with a halving procedure for parameter selection, guaranteeing positivity and convergence for solving the nonlinear eigenvalue problem.

Findings

Method converges globally with quadratic rate

Ensures positivity of solutions during iteration

Numerical experiments validate theoretical results

Abstract

In this paper, we propose an iterative method to compute the positive ground states of saturable nonlinear Schr\"odinger equations. A discretization of the saturable nonlinear Schr\"odinger equation leads to a nonlinear algebraic eigenvalue problem (NAEP). For any initial positive vector, we prove that this method converges globally with a locally quadratic convergence rate to a positive solution of NAEP. During the iteration process, the method requires the selection of a positive parameter $\theta_k$ in the $k$th iteration, and generates a positive vector sequence approximating the eigenvector of NAEP and a scalar sequence approximating the corresponding eigenvalue. We also present a halving procedure to determine the parameters $\theta_k$, starting with $\theta_k=1$ for each iteration, such that the scalar sequence is strictly monotonic increasing. This method can thus be used to illustrate the existence of positive ground states of saturable nonlinear Schr\"odinger equations. Numerical experiments are provided to support the theoretical results.

Tables

Table 1. Table 1: Numerical results for NNI

Parameters			NNI
$n$	$𝐚$		Iter	Residual
$2500$	$𝐚\ge 1$		$6$	2.52e-16
$10000$	$𝐚\ge 1$		$6$	2.38e-16
$40000$	$𝐚\ge 1$		$6$	2.79e-16
$2500$	$1>𝐚>0$		$13$	2.15e-16
$10000$	$1>𝐚>0$		$16$	2.25e-16
$40000$	$1>𝐚>0$		$23$	2.50e-16
$2500$	$𝐚>0$		$13$	1.98e-16
$10000$	$𝐚>0$		$15$	1.34e-15
$40000$	$𝐚>0$		$21$	7.04e-16

Equations

$$i\frac{\partial \varphi }{\partial z}+△\varphi +\Gamma f(|\varphi {|}^2)\varphi =0\text{ for some constant }\Gamma \in \mathbb{R},$$

$$i\frac{\partial \varphi }{\partial z}=-△\varphi +\Gamma \left(1-\frac{1}{a(x)+|\varphi {|}^2}\right)\varphi ,\text{ for }\Gamma >0,$$

$$m(\Gamma ,I)=\inf \{H(u)|u\in H^1({\mathbb{R}}^2),{\int }_{{\mathbb{R}}^2}{u}^2=1\},$$

$$H(u)={\int }_{{\mathbb{R}}^2}|\nabla u{|}^2+\Gamma \left[u^2-\ln \left(1+\frac{u^2}{a(x)}\right)\right]dx.$$

$$-\Delta u+\Gamma \left(1-\frac{1}{a(x)+u^2}\right)u=\lambda u,$$

$$A\mathbf{u}+\Gamma \mathrm{diag}\left(\mathbf{e}-\frac{\mathbf{e}}{\mathbf{a}+{\mathbf{u}}^{[2]}}\right)\mathbf{u}=\lambda \mathbf{u},{\mathbf{u}}^T\mathbf{u}=1,$$

$$A_{ij}=0,\forall i\in S,\forall j\notin S.$$

$$\max \left(\frac{\mathbf{w}}{\mathbf{v}}\right)=\underset{i}{\max }\left(\frac{{\mathbf{w}}^{(i)}}{{\mathbf{v}}^{(i)}}\right),\min \left(\frac{\mathbf{w}}{\mathbf{v}}\right)=\underset{i}{\min }\left(\frac{{\mathbf{w}}^{(i)}}{{\mathbf{v}}^{(i)}}\right).$$

$$\mathcal{A}(\mathbf{u})\mathbf{u}=\lambda \mathbf{u},$$

$$\mathcal{A}(\mathbf{u})=A+\Gamma \mathrm{diag}\left(\mathbf{e}-\frac{\mathbf{e}}{\mathbf{a}+{\mathbf{u}}^{[2]}}\right)$$

$$\mathbf{r}(\mathbf{u,}\lambda)=\mathcal{A}(\mathbf{u})\mathbf{u}-\lambda\mathbf{u},\quad\mathbf{F}(\mathbf{u},\lambda)=\left[\begin{array}[]{c}-\mathbf{r}(\mathbf{u,}\lambda)\\ \frac{1}{2}\left(1-\mathbf{u}^{T}\mathbf{u}\right)\end{array}\right].$$

$$F^{\prime}(\mathbf{u},\lambda)=\left[\begin{array}[]{cc}J({\mathbf{u}})&-\mathbf{u}\\ -\mathbf{u}^{T}&0\end{array}\right],$$

$$J(\mathbf{u})=A+(\Gamma -\lambda )I-\Gamma \mathrm{diag}\left(\frac{\mathbf{a}-{\mathbf{u}}^{[2]}}{(\mathbf{a}+{\mathbf{u}}^{[2]}{)}^{[2]}}\right).$$

$$\displaystyle\left[\begin{array}[]{cc}J(\mathbf{u}_{k})&-\mathbf{u}_{k}\\ -\mathbf{u}_{k}^{T}&0\end{array}\right]\left[\begin{array}[]{c}\Delta_{k}\\ \delta_{k}\end{array}\right]$$

$${\mathbf{u}}_{k+1}$$

$${\hat{\lambda }}_{k+1}$$

$$J({\mathbf{u}}_k)({\Delta }_k+{\mathbf{u}}_k)$$

$${\mathbf{u}}_{k+1}=J({\mathbf{u}}_k{)}^{-1}\left({\delta }_k{\mathbf{u}}_k+2\Gamma \mathrm{diag}\left(\frac{{\mathbf{u}}_k^{[2]}}{(\mathbf{a}+{\mathbf{u}}_k^{[2]}{)}^{[2]}}\right){\mathbf{u}}_k\right).$$

$${\delta }_k>0,J({\mathbf{u}}_k{)}^{-1}\ge 0.$$

$${\lambda }_k=\min \left(\frac{\mathcal{A}({\mathbf{u}}_k){\mathbf{u}}_k}{{\mathbf{u}}_k}\right),$$

$$J({\mathbf{u}}_k)-(\mathcal{A}({\mathbf{u}}_k)-{\lambda }_{k}I)=2\Gamma \mathrm{diag}\left(\frac{{\mathbf{u}}_k^{[2]}}{(\mathbf{a}+{\mathbf{u}}_k^{[2]}{)}^{[2]}}\right),$$

$$F^{\prime}(\mathbf{u},\lambda)=\left[\begin{array}[]{cc}I&0\\ -\mathbf{u}^{T}(J(\mathbf{u}))^{-1}&1\end{array}\right]\left[\begin{array}[]{cc}J({\mathbf{u}})&-\mathbf{u}\\ 0&-\mathbf{u}^{T}(J(\mathbf{u}))^{-1}\mathbf{u}\end{array}\right]$$

$$\displaystyle(F^{\prime}(\mathbf{u},\lambda))^{-1}=\left[\begin{array}[]{cc}(J(\mathbf{u}))^{-1}-\frac{(J(\mathbf{u}))^{-1}\mathbf{u}\mathbf{u}^{T}(J(\mathbf{u}))^{-1}}{\mathbf{u}^{T}(J(\mathbf{u}))^{-1}\mathbf{u}}&-\frac{(J(\mathbf{u}))^{-1}\mathbf{u}}{\mathbf{u}^{T}(J(\mathbf{u}))^{-1}\mathbf{u}}\\ -\frac{\mathbf{u}^{T}(J(\mathbf{u}))^{-1}}{\mathbf{u}^{T}(J(\mathbf{u}))^{-1}\mathbf{u}}&-\frac{1}{\mathbf{u}^{T}(J(\mathbf{u}))^{-1}\mathbf{u}}\end{array}\right].$$

$$\frac{1}{2}(1+{\mathbf{u}}_k^T{\mathbf{u}}_k)-\Gamma {\mathbf{u}}_k^T(J({\mathbf{u}}_k){)}^{-1}\mathrm{diag}\left(\frac{2{\mathbf{u}}_k^{[2]}}{(\mathbf{a}+{\mathbf{u}}_k^{[2]}{)}^{[2]}}\right){\mathbf{u}}_k⩾0.$$

$$\Gamma (J({\mathbf{u}}_k){)}^{-1}\mathrm{diag}\left(\frac{2{\mathbf{u}}_k^{[2]}}{(\mathbf{a}+{\mathbf{u}}_k^{[2]}{)}^{[2]}}\right){\mathbf{u}}_k=\left(I-(J({\mathbf{u}}_k){)}^{-1}(\mathcal{A}({\mathbf{u}}_k)-{\lambda }_{k}I)\right){\mathbf{u}}_k.$$

$$\frac{1}{2}$$

$$=\frac{1}{2}(1+{\mathbf{u}}_k^T{\mathbf{u}}_k)-{\mathbf{u}}_k^T{\mathbf{u}}_k+{\mathbf{u}}_k^T(J({\mathbf{u}}_k){)}^{-1}(\mathcal{A}({\mathbf{u}}_k){\mathbf{u}}_k-{\lambda }_k{\mathbf{u}}_k)$$

$$=\frac{1}{2}(1-{\mathbf{u}}_k^T{\mathbf{u}}_k)+{\mathbf{u}}_k^T(J({\mathbf{u}}_k){)}^{-1}\mathbf{r}({\mathbf{u}}_k,{\lambda }_k).$$

$$\displaystyle F^{\prime}(\mathbf{u_{k}},\lambda_{k})\left[\begin{array}[]{c}{\bf u}_{k}+\Delta_{k}\\ \delta_{k}\end{array}\right]$$

$$\displaystyle=\left[\begin{array}[]{c}\Gamma{\rm diag}\left(\frac{2{\bf u}_{k}^{[2]}}{({\bf a}+{\bf u}_{k}^{[2]})^{[2]}}\right){\bf u}_{k}\\ -\frac{1}{2}\left(1+{\bf u}_{k}^{T}{\bf u}_{k}\right)\end{array}\right].$$

$${\mathbf{u}}_k+{\Delta }_k=$$

$$+\frac{1+{\mathbf{u}}_k^T{\mathbf{u}}_k}{2{\mathbf{u}}_k^T(J({\mathbf{u}}_k){)}^{-1}{\mathbf{u}}_k}(J({\mathbf{u}}_k){)}^{-1}{\mathbf{u}}_k$$

$$=$$

$$+\frac{\frac{1}{2}(1+{\mathbf{u}}_k^T{\mathbf{u}}_k)-\Gamma {\mathbf{u}}_k^T(J({\mathbf{u}}_k){)}^{-1}\mathrm{diag}\left(\frac{2{\mathbf{u}}_k^{[2]}}{(\mathbf{a}+{\mathbf{u}}_k^{[2]}{)}^{[2]}}\right){\mathbf{u}}_k}{{\mathbf{u}}_k^T(J({\mathbf{u}}_k){)}^{-1}{\mathbf{u}}_k}(J({\mathbf{u}}_k){)}^{-1}{\mathbf{u}}_k.$$

Full text

1 Introduction

The nonlinear Schrödinger (NLS) equation [14] is a nonlinear variation of the Schrödinger equation and is a general model in nonlinear science and mathematics. Such an equation can be expressed as follows:

[TABLE]

where $\phi=\phi(x,z):\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{C}$, the function $f$ denotes the nonlinearity and $i$ is the imaginary unit. A NLS equation is called a saturable NLS equation [3, 9] if the nonlinear function $f(s)=1-1/(a+s^{2})$, that is,

[TABLE]

where $a(x)>0$ is a bounded function. A saturable NLS equation is of interest in several applications [5, 7, 8, 10, 11, 12], and has been extensively studied in the past thirty years. In many application areas, one is interested in finding the ground state vector of equation (2). The ground state of equation (2) is defined as the minimizer of the energy function, which is determined by the following constrained optimization problem [3, 9]:

[TABLE]

where

[TABLE]

Therefore, the associated Euler-Lagrange equation of (3) is as follows:

[TABLE]

where $a(x)）>0,\int_{-\infty}^{\infty}u^{2}(x)dx=1$, $(\lambda,u)$ is the eigenpair. In general, the eigenfunction $u(x)$ describes the probability distribution of finding a particle in a particular region in space. Therefore, the existence of positive solutions $u(x)$ [9] and the problem of computing these solutions has attracted much attention in recent years. Here we consider the finite-difference discretization of the nonlinear eigenvalue problem (4) with Dirichlet boundary conditions, and the discretization gives a nonlinear algebraic eigenvalue problem (NAEP)

[TABLE]

where $\mathbf{a}>0,\Gamma>0$, $\mathbf{u}=[u_{1},u_{2},\ldots,u_{n}]^{T}\in\mathbb{R}^{n},$ $\mathbf{u}^{[2]}=[u_{1}^{2},u_{2}^{2},\ldots,u_{n}^{2}]^{T},$ $A$ is an irreducible nonsingular $M$-matrix and $\mathbf{e}=[1,\ldots,1]^{T}$. We aim to provide a structure-preserving algorithm with fast convergence rate for computing positive eigenvectors $\mathbf{u}_{\ast}$ and eigenvalues $\lambda_{\ast}$of NAEP (5) and giving a detailed convergence analysis.

In many applications, the positivity structure of the approximate solutions is important; if the approximations lose positivity structure, then they may be meaningless and unexplained. Therefore, in this paper, we propose a positivity preserving iteration for nonlinear algebraic eigenvalue problems (5) by combining the idea of Newton’s method with the idea of the Noda iteration [13], called the Newton-Noda iteration (NNI). NNI is a Newton iterative method with a new type of full Newton steps, it has the advantage that no line-searches are needed, and naturally preserves the strict positivity of the target eigenvector $\mathbf{u}_{\ast}$ in its approximations at all iterations. We also present a halving procedure to determine the parameters $\theta_{k}$, starting with $\theta_{k}=1$ for each iteration, such that the sequence approximating target eigenvalue $\lambda_{\ast}$ is strictly monotonic increasing and bounded, and thus its global convergence is guaranteed. Another advantage of NNI is that it converges quadratically and computes the desired eigenpair correctly for any positive initial vector.

The rest of this paper is organized as follows. In Section 2, we present a Newton-Noda iteration. In Section 3, we prove some basic properties for Newton-Noda iteration. Section 4 addresses the global convergence and the local convergence rate of NNI. In Section 5, we provide numerical examples to verify the theoretical results and the performance of NNI. Some concluding remarks are given in the last section.

Throughout this paper, we use the bold face letters to denote a vector and use the $2$-norm for vectors and matrices. The superscript $T$ denotes the transpose of a vector or matrix, and we use $\mathbf{v}^{(i)}$ to represent the $i$th element of a vector $\mathbf{v}$. $\mathbf{v}^{[m]}$ denotes element-by-element powers, i.e., $\mathbf{v}^{[m]}=[v_{1}^{m},v_{2}^{m},\ldots,v_{n}^{m}]^{T}.$ A real matrix $A=\left[A_{ij}\right]\in\mathbb{R}^{n\times k}$ is called nonnegative (positive) if $A_{ij}\geq 0$ $(A_{ij}>0)$. For real matrices $A$ and $B$ of the same size, we write $A\geq B$ ($A>B$) if $A-B$ is nonnegative (positive). A real square matrix $A$ is called a $Z$-matrix if all its off-diagonal elements are nonpositive. A matrix $A$ is called a M-matrix if it is a Z-matrix with $A^{-1}\geq 0$. A matrix $A$ is called reducible [2, 6] if there exists a nonempty proper index subset $S\subset\left\{1,2,\ldots,n\right\}$ such that

[TABLE]

If $A$ is not reducible, then we call $A$ irreducible. For a pair of positive vectors $\mathbf{v}$ and $\mathbf{w}$, define

[TABLE]

2 The Newton-Noda iteration

In this section, we will present a Newton-Noda iteration (NNI) for computing a positive eigenvector $\mathbf{u}_{\ast}$ of NAEP (5), and then we prove some basic properties of NNI in Section 3, which will be used to establish its convergence theory in Section 4.

First, NAEP (5) can be simplified as follows:

[TABLE]

where

[TABLE]

and $\mathrm{diag}\left(\ast\right)$ returns a square diagonal matrix with the elements of vector $\ast$ on the main diagonal. We define two vector-valued functions $\mathbf{r}:$ $\mathbb{R}_{+}^{n+1}\mathbb{\rightarrow R}^{n}$ and $F:$ $\mathbb{R}_{+}^{n+1}\mathbb{\rightarrow R}^{n+1}$ as follows:

[TABLE]

The Fréchet derivative of $F$ is given by

[TABLE]

where

[TABLE]

Next, we consider using Newton’s method to solve the equation $\mathbf{F}(\mathbf{u,}\lambda)=0$. Given an approximation $(\mathbf{u}_{k},\widehat{\lambda}_{k})$, Newton’s method produces the next approximation $(\mathbf{u}_{k+1},\widehat{\lambda}_{k+1})$ as follows:

[TABLE]

From the first equation of (14), we have

[TABLE]

Hence,

[TABLE]

Since $\mathbf{u}_{k}$ is going to approximate the positive eigenvector of NAEP, we will also require $\mathbf{u}_{k}>0$. However, we cannot guarantee $\mathbf{u}_{k+1}>0$ in (17) unless we have

[TABLE]

What is needed here is that $J(\mathbf{u}_{k})$ is a nonsingular M-matrix. For $\mathbf{u}_{k}>0$, we suggests taking

[TABLE]

which is precisely the idea of the Noda iteration [13]. This implies that the $Z$-matrix $\mathcal{A}(\mathbf{u}_{k})-\lambda_{k}I$ is such that $(\mathcal{A}(\mathbf{u}_{k})-\lambda_{k}I)\mathbf{u}_{k}\geq 0$. Thus $\mathcal{A}(\mathbf{u}_{k})-\lambda_{k}I$ is a nonsingular $M$-matrix when $(\mathbf{u}_{k},\lambda_{k})$ is not an eigenpair, and is a singular $M$-matrix when $(\mathbf{u}_{k},\lambda_{k})$ is an eigenpair. Since

[TABLE]

we have $J(\mathbf{u}_{k})\mathbf{u}_{k}>0$. Thus $J(\mathbf{u}_{k})$ is a nonsingular $M$-matrix. Based on (14), (15) and (18), we can present NNI as Algorithm 2.1.

In what follows, we will prove the positivity of $\mathbf{u}_{k}$ and give a strategy for choosing $\theta_{k}$. These results will show that Algorithm 2.1 is a positivity preserving algorithm.

2.1 Positivity of $\mathbf{u}_{k}$

Suppose that $\left\{\mathbf{u}_{k},{\lambda}_{k}\right\}$ is generated by Algorithm 2.1. We now prove that the parameter $\theta_{k}\in(0,1]$ in Algorithm 2.1 naturally preserves the strict positivity of $\mathbf{u}_{k}$ at all iterations.

For any vector $\mathbf{u}>0$, from (7), it follows that

[TABLE]

is nonsingular and

[TABLE]

Lemma 1**.**

Given $\mathbf{a},\Gamma>0$. Suppose that $\lambda_{k}\in\mathbb{R}$ and $\mathbf{u}_{k}>0$ with $\|\mathbf{u}_{k}\|=1$ such that $\mathbf{r}(\mathbf{u}_{k},\lambda_{k})\geqslant 0$. Then

[TABLE]

Moreover, the equality holds if and only if $\mathbf{r}(\mathbf{u}_{k},\lambda_{k})=0$.

Proof.

Since $\mathcal{A}({\bf u}_{k}){\bf u}_{k}-\lambda_{k}{\bf u}_{k}\geqslant 0$, $J({\bf u}_{k})=(\mathcal{A}({\bf u}_{k})-\lambda_{k}I)+\Gamma{\rm diag}\left(\frac{2{\bf u}_{k}^{[2]}}{({\bf a}+{\bf u}_{k}^{[2]})^{[2]}}\right)$ is nonsingular $M$-matrix. Then we have

[TABLE]

Hence,

[TABLE]

Since $\|{\bf u}_{k}\|=1$ and $\mathbf{r}(\mathbf{u}_{k},\lambda_{k})\geqslant 0$, (24) holds by using ${\bf u}_{k}^{T}(J({\bf u}_{k}))^{-1}>0$. It is easily seen that the equality of (24) holds if and only if $\mathbf{r}(\mathbf{u}_{k},\lambda_{k})=0$. ∎

Theorem 2**.**

Given $\mathbf{a},\Gamma>0$. Assume $\left\{\mathbf{u}_{k},{\lambda}_{k}\right\}$ is generated by Algorithm 2.1. If $\theta_{k}\in(0,1]$, then $\mathbf{u}_{k}>0$ for all $k$.

Proof.

Since ${\bf u}_{0}>0$, by mathematical induction, it suffices to show that if ${\bf u}_{k}>0$ then ${\bf u}_{k+1}>0$. Suppose that ${\bf u}_{k}>0$, it follows from the step 3 of Algorithm 2.1 that

[TABLE]

By (23), we have

[TABLE]

Since $\Gamma>0$, ${\bf u}_{k}>0$, $J({\bf u}_{k})$ is a nonsingular $M$-matrix and $\mathbf{u}_{k}^{T}(J({\bf u}_{k}))^{-1}{\bf u}_{k}>0$, it follows from Lemma 1 that ${\bf u}_{k}+\Delta_{k}>0$. Therefore, ${\bf w}_{k+1}={\bf u}_{k}+\theta_{k}\Delta_{k}>0$ if $0<\theta_{k}\leq 1$, and hence, $\mathbf{u}_{k+1}={\mathbf{w}}_{k+1}/\|{\mathbf{w}}_{k+1}\|>0$. ∎

Remark 1**.**

---

1. (i)

$\mathbf{u}_{k}^{T}\Delta_{k}=0$: From (2.1), it is easily seen that $\mathbf{u}_{k}^{T}\Delta_{k}=\frac{1}{2}(1-\mathbf{u}_{k}^{T}\mathbf{u}_{k})=0.$ 2. (ii)

$\delta_{k}\geqslant 0$: From the step 3 of Algorithm 2.1 and using (23), we have

[TABLE]

Lemma 3**.**

If $\delta_{k}$, $\Delta_{k}$ and $\mathbf{u}_{k}$ are generated by Algorithm 2.1, then the following statements are equivalent:

[TABLE]

Proof.

From the step 3 of Algorithm 2.1, we have

[TABLE]

(i)$\Rightarrow$(ii): From (ii) of Remark 1, we get $\delta_{k}=0$ if and only if $\mathbf{r}(\mathbf{u}_{k},\lambda_{k})=0$.

(ii)$\Rightarrow$(iii): Since $\mathbf{r}(\mathbf{u}_{k},\lambda_{k})=0$ and $F^{\prime}(\mathbf{u_{k}},\lambda_{k})$ is a nonsingular matrix, we have $\Delta_{k}=0$ and $\delta_{k}=0$.

(iii) $\Rightarrow$ (i): If $\Delta_{k}=0$, then

[TABLE]

and it follows

[TABLE]

which implies

[TABLE]

Then

[TABLE]

which means $\delta_{k}=0$. ∎

2.2 The strategy for choosing $\theta_{k}$

In this section, we would like to choose $\theta_{k}\in(0,1]$ such that the sequence $\left\{\lambda_{k}\right\}$ is strictly increasing and bounded above.

Lemma 4**.**

Given a unit vector $\mathbf{u}_{k}>0$ and $\theta_{k}\in(0,1]$, then

[TABLE]

where $\mathbf{h}_{k}(\theta_{k})=\mathbf{r}(\mathbf{u}_{k+1},\lambda_{k})$. Moreover, $\mathbf{h}_{k}(\theta_{k})$ can be also expressed in the form

[TABLE]

where $\|\mathbf{R}(\theta_{k}\Delta_{k})\|\leq M\|\theta_{k}\Delta_{k}\|^{2}$.

Proof.

By Theorem 2, we know that $\mathbf{u}_{k}>0$ for all $k$.

[TABLE]

where $\mathbf{h}_{k}(\theta_{k})=\mathbf{r}(\mathbf{u}_{k+1},\lambda_{k})$. By Taylor’s theorem, we have

[TABLE]

where $\|\mathbf{E}_{k}\|\leq M_{1}\|\mathbf{u}_{k+1}-\mathbf{u}_{k}\|^{2}$.

Since $\mathbf{u}_{k}^{T}\Delta_{k}=0$ from Remark 1, we have $\|{\mathbf{w}}_{k+1}\|=\sqrt{1+\|\theta_{k}\Delta_{k}\|^{2}}.$ Hence, the third term in the right-hand side of (28) is bounded by

[TABLE]

and the upper bound of $\|\mathbf{E}_{k}\|$ can be re-estimated as follows:

[TABLE]

From the above relation (28)-(30), we have

[TABLE]

where

[TABLE]

with $\|\mathbf{R}(\theta_{k}\Delta_{k})\|\leq M\|\theta_{k}\Delta_{k}\|^{2}$ and $M=\frac{\Gamma}{2\min(\mathbf{a})}+\frac{M_{1}}{2}+M_{1}$. ∎

We next show that $\lambda_{k}$ is strictly increasing and bounded above for suitable $\theta_{k}$, unless $\mathbf{u}_{k}$ is an eigenvector of NAEP for some $k$, in which case NNI terminates with $\lambda_{k}$.

Theorem 5**.**

Suppose $A$ be an irreducible M-matrix and $\eta>0$ be a fixed constant. Given a unit vector $\mathbf{u}_{k}>0$, suppose $\mathbf{u}_{k}\not=\mathbf{u}_{\ast}$ and $\theta_{k}$ in Algorithm 2.1 satisfies

[TABLE]

where for each $k$ with $\mathbf{h}_{k}(1)<\frac{\delta_{k}\mathbf{u}_{k}}{(1+\eta)\|{\mathbf{w}}_{k+1}\|}$,

[TABLE]

Then $0<\eta_{k}<1$ whenever it is defined, and

[TABLE]

Proof.

By Lemma 4, we have

[TABLE]

We need to prove $\mathbf{h}_{k}(\theta_{k})>0.$

From (27) and $\|\mathbf{R}(\theta_{k}\Delta_{k})\|\leq M\|\theta_{k}\Delta_{k}\|^{2}$, we have

[TABLE]

If $\eta_{k}\geq 1$, then

[TABLE]

and it follows

[TABLE]

Thus

[TABLE]

If $\eta_{k}<1$, we have

[TABLE]

which ensures the inequality

[TABLE]

Substituting (36) into (33), we obtain

[TABLE]

Therefore,

[TABLE]

Next, we prove that the sequence $\left\{\lambda_{k}\right\}$ is bounded above. Suppose that $\left\{\lambda_{k}\right\}$ is unbounded. This implies that $\lambda_{k}\geq N>0$ for $k$ large enough. Since $A(\mathbf{u}_{k}){\bf u}_{k}\geq{\lambda}_{k}\mathbf{u}_{k}$, we then have

[TABLE]

which is a contradiction. ∎

From (31), we know that the inequality $\mathbf{h}_{k}(1)\geq\frac{\delta_{k}\mathbf{u}_{k}}{(1+\eta)\|{\mathbf{w}}_{k+1}\|}$ depends on the parameter $\eta$. Therefore, if $\eta$ large enough, then we can choose $\theta_{k}=1$ for which $\mathbf{h}_{k}(1)>0$ holds. By Theorem 5, we can indeed choose $\theta_{k}\in(0,1]$ in NNI such that the sequence $\{\lambda_{k}\}$ is strictly increasing. However, in practice, it is difficult to determine $\eta_{k}$. Therefore, we can determine $\theta_{k}$ by repeated halving technique. More precisely, for each $k$, we can take $\theta_{k}=1$ first and check whether $\mathbf{h}_{k}(1)>0$ holds. If not, then we update $\theta_{k}$ using $\theta_{k}\leftarrow\theta_{k}/2$ and check again until we get $\theta_{k}$ for which $\mathbf{h}_{k}(1)>0$ holds. This process of repeatedly halving will be referred to as the halving procedure. As long as $\theta_{k}$ is bounded below by a positive constant, which will be mentioned in the next section.

3 Some basic properties of Newton-Noda iteration

In this section, we prove a number of basic properties of NNI, which will be used to establish its convergence theory in Section 4.

Lemma 6**.**

Let $A$ be an irreducible M-matrix. Assume that the sequence $\left\{\lambda_{k},\mathbf{u}_{k},\mathbf{w}_{k}\right\}$ is generated by Algorithm 2.1. For any subsequence $\left\{\mathbf{u}_{k_{j}}\right\}\subseteq\left\{\mathbf{u}_{k}\right\},$ we have the following results:

1. (i)

If $\mathbf{u}_{k_{j}}\rightarrow\mathbf{v}$ as $j\rightarrow\infty,$ then $\mathbf{v}>0.$ 2. (ii)

$\min(\mathbf{u}_{k})\geq m$* for some positive constant $m$.* 3. (iii)

$\|\mathbf{w}_{k}\|\leq\frac{1}{m}$.

Proof.

(i). If $\lim_{j\rightarrow\infty}\mathbf{u}_{k_{j}}=\mathbf{v,}$ then $\mathbf{v}\geq 0$. Let $S$ be the set of all indices $i$ such that $\lim_{j\to\infty}\mathbf{u}_{k_{j}}^{(i)}=\mathbf{v}^{(i)}=0$. Since $\left\|\mathbf{u}_{k_{j}}\right\|=1$, $S$ is a proper subset of $\{1,2,\ldots,n\}$. Suppose $S$ is nonempty. Then by the definition of $\lambda_{k},$

[TABLE]

Since $\lim_{j\rightarrow\infty}\mathbf{u}_{k_{j}}^{(i)}=0$ for $i\in S$, it holds that $\lim_{j\rightarrow\infty}\left(\mathcal{A}(\mathbf{u}_{k_{j}}){\bf u}_{k_{j}}\right)^{(i)}=\left(\mathcal{A}(\mathbf{v}){\bf v}\right)^{(i)}=0$ for $i\in S$. Thus, $\mathcal{A}(\mathbf{v})_{i,j}=0$ for all $i\in S$ and for all $j\notin S$, which contradicts the irreducibility of $\mathcal{A}(\mathbf{v})$. Therefore, $S$ is empty and thus $\mathbf{v}>0$.

(ii). Suppose $\min(\mathbf{u}_{k})$ is not bounded below by a positive constant. Then there exists a subsequence $\{k_{j}\}$ such that $\lim_{j\rightarrow\infty}\min(\mathbf{u}_{k_{j}})=0$. Since $\|\mathbf{u}_{k_{j}}\|=1$, we may assume that $\lim_{j\to\infty}\mathbf{u}_{k_{j}}=\mathbf{v}$ exists. Then $\lim_{j\rightarrow\infty}\min(\mathbf{u}_{k_{j}})=\min(\mathbf{v})=0$. This is a contradiction since $\mathbf{v}>0$ by (i). Therefore, $\min(\mathbf{u}_{k})$ is bounded below by a positive constant. That is $\min(\mathbf{u}_{k})\geq m$ for some positive constant $m$.

(iii). From Remark 1, we have $\mathbf{u}_{k}^{T}\mathbf{w}_{k+1}=1$ and then

[TABLE]

Since $\mathbf{u}_{k}>0$ and $\mathbf{u}_{k+1}>0$ with $\|\mathbf{u}_{k}\|=\|\mathbf{u}_{k+1}\|=1$, we have

[TABLE]

where $\|\cdot\|_{1}$ is the vector 1-norm. Form (ii),

[TABLE]

∎

Lemma 7**.**

Assume that the sequence $\left\{\Delta_{k},\delta_{k},\theta_{k}\right\}$ is generated by Algorithm 2.1. We have the following results:

1. (i)

There exists a constant $\beta>0$ such that $\beta\|\Delta_{k}\|\leq\delta_{k}$. 2. (ii)

$\theta_{k}=1$* if $\|\Delta_{k}\|\leq\frac{\eta\beta}{(1+\eta)M}$ where $(\eta,M)$ is as in Theorem 5.* 3. (iii)

$\theta_{k}\geq\xi$* for some positive constant $\xi$.*

Proof.

(i). From the step 3 of Algorithm 2.1, we have

[TABLE]

Since $\|J(\mathbf{u})^{-1}\mathbf{u}\|$ a continuous function achieves its extreme values in a compact set, it follows

[TABLE]

Therefore, $\|J(\mathbf{u})^{-1}\mathbf{u}\|\leq M_{2}$ for some constant $M_{2}$.

On the other hand, from (ii) of Remark 1, we have

[TABLE]

Since $\mathbf{u}_{k}>0$ and $J(\mathbf{u}_{k})^{-1}\mathbf{r}(\mathbf{u}_{k},\lambda_{k})>0$, by using the same proving technique of (iii) of Lemma 6, we have

[TABLE]

which implies

[TABLE]

From (39) and the above inequality,

[TABLE]

(ii). If $\|\Delta_{k}\|\leq\frac{\eta\beta}{(1+\eta)M}$, then

[TABLE]

From the proof of Theorem 5, $\theta_{k}=1$ when $\eta_{k}\geq 1.$

(iii). From (31), we recall that

[TABLE]

where $\eta_{k}=\frac{\eta\delta_{k}\min\left(\mathbf{u}_{k}\right)}{(1+\eta)M\|{\mathbf{w}}_{k+1}\|\left\|\Delta_{k}\right\|^{2}}.$ Suppose $\theta_{k}$ is not bounded below by $\xi>0$. Since $\mathbf{u}_{k}$ is bounded, we can find a subsequence $\{k_{j}\}$ such that

[TABLE]

Note that $\mathbf{v}>0$ by Lemma 6.

From (20), $F^{\prime}(\mathbf{u}_{k_{j}},\lambda_{k_{j}})$ is a nonsingular matrix, and the vector $\left[\Delta_{k_{j}}^{T},\delta_{k_{j}}\right]^{T}$ satisfies

[TABLE]

Since the sequence $\left\{\lambda_{k}\right\}$ is monotonically increasing and bounded above, we have $\lim_{k\rightarrow\infty}\lambda_{k}=\alpha$. Therefore,

[TABLE]

which means $\left\|\Delta_{k_{j}}\right\|$ is bounded. If $\eta_{k}$ is defined only on a finite subset of $\{k_{j}\}$, then $\theta_{k_{j}}=1$ except for a finite number of $j$ values, contradicting $\lim_{j\rightarrow\infty}\theta_{k_{j}}=0$. If $\eta_{k}$ is defined on an infinite subset $\{k_{j_{i}}\}$ of $\{k_{j}\}$, then

[TABLE]

It follows that $\lim_{i\rightarrow\infty}\|\Delta_{k_{j_{i}}}\|=\infty.$ This is contradictory to $\|\Delta_{k_{j}}\|<\infty$. ∎

4 Convergence analysis

In this section, we prove that the convergence of NNI is global and quadratic, assuming that $\mathbf{u}_{k}\neq\mathbf{u}_{\ast}$ for each $k$.

4.1 Global convergence of NNI

Theorem 5 shows that the sequence $\left\{\lambda_{k}\right\}$ is strictly increasing and bounded above by a constant and hence converges. We now show that the limit of $\lambda_{k}$ is precisely the eigenvalue $\lambda_{\ast}$ of NAEP (5).

Theorem 8**.**

Let $A$ be an irreducible M-matrix and the sequence $\left\{\lambda_{k}\right\}$ is generated by Algorithm 2.1. If $\mathbf{a},\Gamma>0$, then the NAEP (5) has a positive eigenvecor.

Proof.

From (26), (37) and Lemma 7, we have

[TABLE]

From (iii) of Lemma 6, we have $\|\mathbf{w}_{k+1}\|\leq\frac{1}{m}<\infty$. It follows from (42) that $\lim_{k\rightarrow\infty}\delta_{k}\min(\mathbf{u}_{k})=0$. From (ii) of Lemma 6, $\min(\mathbf{u}_{k})$ is bounded below by a positive constant, and thus $\lim_{k\rightarrow\infty}\delta_{k}=0.$

Let $\mathbf{v}$ be any limit point of $\{\mathbf{u}_{k}\}$, with $\lim_{j\to\infty}\mathbf{u}_{k_{j}}=\mathbf{v}>0$. From Lemma 3, we then have $\lim_{j\to\infty}\delta_{k_{j}}=0$ if and only if $\lim_{j\to\infty}\left(\mathcal{A}(\mathbf{u}_{k_{j}})\mathbf{u}_{k_{j}}-\lambda_{k_{j}}\mathbf{u}_{k_{j}}\right)=0$, which means $\mathcal{A}(\mathbf{v})\mathbf{v}=\lambda\mathbf{v}$. Therefore, $\mathbf{v}$ is a positive eigenvector of NAEP and $\lambda=\min\left(\frac{\mathcal{A}(\mathbf{v})\mathbf{v}}{\mathbf{v}}\right)$ is the corresponding eigenvalue, i.e., $\mathbf{u}_{\ast}=\mathbf{v}$ and $\lambda_{\ast}=\lambda.$ ∎

The above theorem guarantees the global convergence of NNI and also proves the existence of positive eigenvectors of NAEP.

4.2 Quadratic convergence of NNI

In the previous section, we discussed the global convergence of NNI. In the following section, we will establish a convergence rate analysis by exploiting a connection between NNI and Newton’s method. So we start with the following result about Newton’s method.

Lemma 9**.**

Suppose that $\left(\mathbf{u}_{k},\lambda_{k}\right)$ form Algorithm 2.1 is sufficiently close to an eigenpair $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$ with $\mathbf{u}_{\ast}>0$ and $\|\mathbf{u}_{\ast}\|=1$. Let $\left\{\allowbreak\widehat{\mathbf{u}}_{k},\widehat{\lambda}_{k}\right\}$ be obtained by Newton’s method from $\left(\mathbf{u}_{k},\lambda_{k}\right)$, i.e.,

[TABLE]

where $\Delta_{k}$ and $\delta_{k}$ as in Algorithm 2.1. Then there is a constant $\beta$ such that for all $\left(\mathbf{u}_{k},\lambda_{k}\right)$ sufficiently close to $\left(\mathbf{u}_{\ast},\lambda_{\ast})\right)$

[TABLE]

Proof.

We already know that $F^{\prime}\left(\mathbf{u}_{k},\lambda_{k}\right)$ is nonsingular. It is also clear that $F^{\prime}\left(\mathbf{u},\lambda\right)$ satisfies a Lipschitz condition at $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$ since its Fréchet derivative is continuous in a neighborhood of $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$. The inequality (43) is then a basic result of Newton’s method. ∎

Lemma 10**.**

Let $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$ be an eigenpair with $\mathbf{u}_{\ast}>0$ and $\|\mathbf{u}_{\ast}\|=1$. Let $\left\{\mathbf{u}_{k},\lambda_{k}\right\}$ be generated by NNI. Then there are constants $c_{2}>0$ such that $\left|\lambda_{k}-\lambda_{\ast}\right|\leq c_{2}\|\mathbf{u}_{k}-\mathbf{u}_{\ast}\|$ for all $\mathbf{u}_{k}$ sufficiently close to $\mathbf{u}_{\ast}$.

Proof.

Since

[TABLE]

Since the Fréchet derivative of $\frac{\mathcal{A}(\mathbf{u})\mathbf{u}}{\mathbf{u}}$ is continuous in a neighborhood of $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$, we have $\left|\lambda_{k}-\lambda_{\ast}\right|\leq c_{2}\|\mathbf{u}_{k}-\mathbf{u}_{\ast}\|$ for all $\left(\mathbf{u}_{k},\lambda_{k}\right)$ sufficiently close to $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$. ∎

We now prove the local quadratic convergence of Algorithm 2.1.

Theorem 11**.**

Assume $\left\{\mathbf{u}_{k},\lambda_{k}\right\}$ be generated by NNI. Suppose that $\left(\mathbf{u}_{k_{0}},\lambda_{k_{0}}\right)$ is sufficiently close to an eigenpair $\left(\mathbf{u}_{\ast},\lambda_{\ast}\right)$ with $\mathbf{u}_{\ast}>0$ and $\|\mathbf{u}_{\ast}\|=1$. Then $\lambda_{k}$ converges to $\lambda_{\ast}$ quadratically and $\mathbf{u}_{k}$ converges to $\mathbf{u}_{\ast}$ quadratically.

Proof.

For some $\delta\in(0,\min\mathbf{u}_{\ast})$, there are positive constants $c_{1}$, $c_{2}$ and $c_{3}$ such that

[TABLE]

whenever \left\|\left[\begin{array}[]{c}{\mathbf{u}}_{k}\\ {\lambda}_{k}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{u}_{\ast}\\ \lambda_{\ast}\end{array}\right]\right\|<\delta,

[TABLE]

whenever $\left\|{\mathbf{u}}_{k}-\mathbf{u}_{\ast}\right\|<\delta$, and

[TABLE]

whenever \left\|\left[\begin{array}[]{c}\widehat{\mathbf{u}}_{k+1}\\ \widehat{\lambda}_{k+1}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{u}_{\ast}\\ \lambda_{\ast}\end{array}\right]\right\|<\delta. Note that $\mathbf{u}_{k}>0$ is guaranteed. Now for all $\epsilon>0$ sufficiently small,assume that \left\|\left[\begin{array}[]{c}{\mathbf{u}}_{k}\\ {\lambda}_{k}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{u}_{\ast}\\ \lambda_{\ast}\end{array}\right]\right\|<\epsilon for $k=k_{0}$. By (44) and (45) we have (with $\epsilon\leq\delta$)

[TABLE]

By (46) and (44) we have (with $\epsilon\leq\delta$, $c\epsilon^{2}\leq\delta$)

[TABLE]

Then $\|\widehat{\mathbf{u}}_{k+1}\|\geq\frac{1}{2}$ (with $c_{3}c(1+c_{2})^{2}\epsilon^{2}<\frac{3}{8}$). Now

[TABLE]

Then

[TABLE]

For $\epsilon$ with $(1+\frac{4}{3}c_{3})c(1+c_{2})^{2}\epsilon\leq\frac{1}{1+c_{2}}$, we have $\|{\mathbf{u}}_{k+1}-\mathbf{u}_{\ast}\|<\frac{1}{1+c_{2}}\epsilon$ and thus $\left|{{\lambda}}_{k+1}-\lambda_{\ast}\right|\leq c_{2}\|\mathbf{u}_{k+1}-\mathbf{u}_{\ast}\|<\frac{c_{2}}{1+c_{2}}\epsilon$. Therefore, \left\|\left[\begin{array}[]{c}{\mathbf{u}}_{k+1}\\ {\lambda}_{k+1}\end{array}\right]-\left[\begin{array}[]{c}\mathbf{u}_{\ast}\\ \lambda_{\ast}\end{array}\right]\right\|=\|{\mathbf{u}}_{k+1}-\mathbf{u}_{\ast}\|+\left|{\lambda}_{k+1}-\lambda_{\ast}\right|<\epsilon. We can then repeat the above process to get $\|{\mathbf{u}}_{k+1}-\mathbf{u}_{\ast}\|\leq d\|{\mathbf{u}}_{k}-\mathbf{u}_{\ast}\|^{2}$ for all $k\geq k_{0}$ and $d=\left(1+\frac{4}{3}c_{3}\right)c(1+c_{2})^{2}$. Thus $\mathbf{u}_{k}$ converges to $\mathbf{u}_{\ast}$ quadratically and then ${\lambda}_{k}$ converges to $\lambda_{\ast}$ quadratically by (45).

∎

5 Numerical experiments

In this section, we present some numerical results to support our theory for NNI and illustrate its effectiveness. All numerical tests were performed on 4.2GHz quad-core Intel Core i7 with $32$ GB memory using Matlab R$2018b$ with machine precision $\varepsilon=2.22\times 10^{-16}$ under the macOS High Sierra. Throughout the experiments, the initial vector is $\mathbf{u}_{0}=\frac{1}{\sqrt{n}}[1,\ldots,1]^{T}\in\mathbb{R}^{n}$. In the experiments, the stopping criterion for NNI is the relative residual

[TABLE]

where we use the cheaply computable $(\|\mathcal{A}(\mathbf{u}_{k})\|_{1}\|\mathcal{A}(\mathbf{u}_{k})\|_{\infty})^{1/2}$ to estimate the 2-norm $\|\mathcal{A}(\mathbf{u}_{k})\|$, which is more reasonable than the individual $\|\mathcal{A}(\mathbf{u}_{k})\|_{1}$ or $\|\mathcal{A}(\mathbf{u}_{k})\|_{\infty}$ with $\|\cdot\|_{\infty}$ the infinity norm of a matrix.

Example 1**.**

Consider the finite-difference discretization of the nonlinear eigenvalue problem (4) with Dirichlet boundary conditions on $[0,1]\times[0,1]$, i.e.,

[TABLE]

where $A\in\mathbb{R}^{n\times n}$ is a negative 2D Laplacian matrix with Dirichlet boundary conditions.

For Example 1, Figure 1 depicts how the relative residual evolves versus the number of iterations for NNI. It shows that NNI uses $8$ iterations to achieve the required accuracy, clearly indicating its quadratic convergence.

Figure 2 shows that the magnitude of parameter $\Gamma$ affects the total number of iterations to achieve convergence. As we see, NNI requires more iterations to achieve convergence for lager parameter $\Gamma$.

Table 1 reports the results obtained by NNI. In the table, $n$ specifies the dimension, $\mathbf{a}$ is a parameter to adjust the diagonal elements of $\mathrm{diag}\left(\mathbf{e}-\frac{\mathbf{e}}{\mathbf{a}+\mathbf{u}^{[2]}}\right)$. “$\mathbf{a}\geq 1$”denotes that each element of $\mathbf{a}$ is larger than $1$, “$1>\mathbf{a}>0$”denotes that each element of $\mathbf{a}$ is between [math] and $1$, “$\mathbf{a}>0$”denotes that each element of $\mathbf{a}$ is larger than [math]. “Iter”denotes the number of iterations to achieve convergence, “Residual”denotes the relative residual when NNI is terminated. From the table, we see that the number of iterations for NNI is at most $23$, clearly indicating its rapid convergence. For this example, $\mathbf{h}_{k}(\theta_{k})>0$ holds with $\theta_{k}=1$ for each iteration of NNI and the halving procedure is not used. These results indicate that our theory of NNI can be conservative.

6 Conclusion

In this paper, we are concerned with the nonlinear algebraic eigenvalue problem (NAEP) generated by the discretization of the saturable nonlinear Schrödinger equation. Based on the idea of Noda’s iteration and Newton’s method, we have proposed an effective method for computing the positive eigenvectors of NAEP, called Newton–Noda iteration. It involves the selection of a positive parameter $\theta_{k}$ in the $k$th iteration. We have presented a halving procedure to determine the parameters $\theta_{k}$, starting with $\theta_{k}=1$ for each iteration, such that the sequence approximating target eigenvalue $\lambda_{\ast}$ is strictly monotonic increasing and bounded, and thus its global convergence is guaranteed for any initial positive unit vector. Additionally, another advantage of the presented method is its local convergence speed. We have shown that the parameter $\theta_{k}$ is chosen eventually equal to $1$ and locally quadratic convergence is achieved. The numerical experiments have indicated that the halving procedure will often return $\theta_{k}=1$ (i.e., no halving is actually used) for each $k$, and near convergence the halving procedure will always return $\theta_{k}=1$. These results confirm our theory and demonstrate that our theoretical results can be realistic and pronounced.

This iterative method has several nice features: Structure Preserving–It maintains positivity in its computation of positive ground state vectors, and its convergence is global and quadratic. Easy-to-implement –The structure of the new algorithm is still very simple, although its convergence analysis is rather involved for nonlinear algebraic eigenvalue problems. On the other hand, it gives an alternative approach to approximate the solution of the nonlinear Schrödinger equation by constructing a sequence. This is precisely the way we use to prove the existence of solutions of the discrete nonlinear Schrödinger equation.

Acknowledgements

I am very grateful to the Ministry of Science and Technology in Taiwan for funding this research, and I would like to thank Chun-Hua Guo, Wen-Wei Lin and Tai-Chia Lin for their valuable comments on this paper.