, , ,

## A Minimal Energy Tracking Method for Non-radially

## Symmetric Solutions of Coupled Nonlinear Schr¨

## odinger

## Equations

Yueh-Cheng Kuoa, Wen-Wei Linb, Shih-Feng Shiehc, Weichung Wangd

a_{Department of Applied Mathematics, National University of Kaohsiung, Kaohsiung 811, Taiwan}

(yckuo@nuk.edu.tw).

b_{Department of Applied Mathematics, National Chiao-Tung University, Hsinchu 300, Taiwan}

(wwlin@math.nctu.edu.tw).

c_{Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan}

(sfshieh@math.ntnu.edu.tw).

d_{Department of Mathematics, National Taiwan University, Taipei 106, Taiwan}

(wwang@math.ntu.edu.tw).

Abstract

We aim at developing methods to track minimal energy solutions of time-independent m-coupled discrete nonlinear Schr¨odinger (DNLS) equations. We first propose a method to find energy minimizers of the 1-component decoupled DNLS equation and use it as the initial point of the m-coupled DNLS equations in a continuation scheme. We then show that the change of local optimality occurs only at the bifurcation points. The fact leads to a minimal energy tracking method that guides the choice of bifurcation branch corresponding to the minimal energy solution curve. By combining all these techniques with a parameter-switching scheme, we successfully compute a non-radially symmetric energy minimizer that can not be computed by existing numerical schemes straightforwardly.

Key words: Coupled nonlinear Schr¨odinger equations, continuation method, ground states, minimal energy, non-radially symmetric solutions.

1. Introduction

The main theme of this article is to compute the energy minimizers of m-coupled discrete nonlinear Schr¨odinger (NLS) equations. We consider the time-dependent m-coupled NLS equations that are defined as

− ι∂
∂tΦj= [4 − V (z)]Φj+ µj|Φj|
2_{Φ}
j+
m
X
i6=j,i=1
βij|Φi|2Φj, (1a)
Φj= Φj(t, z) ∈ C, z ∈ Ω ⊆ Rn, j = 1, . . . , m, (1b)
Φj(t, z) = 0, as z ∈ ∂Ω, (1c)

where the µj are positive constants, n ≤ 3, V (z) > 0 and βij = βji (i 6= j) are coupling

coefficients. If Ω = Rn_{, the boundary condition (1c) becomes}

Φj(t, z) → 0, as |z| → ∞, t > 0.

To obtain the solitary wave solutions, we set Φj(t, z) = e−ιλjtφj(z) and transform (1)

into the time-independent m-coupled NLS equations

[4 − V (z)]φj− λjφj+ µj|φj|2φj+ m X i6=j,i=1 βij|φi|2φj= 0, in Rn, φj > 0 in Ω ⊆ Rn, j = 1, . . . , m, φj(z) = 0, as z ∈ ∂Ω.

To solve Equations (2) numerically, we consider the corresponding m-coupled discrete nonlinear Schr¨odinger (DNLS) equations

Auj− λjuj+ µju 2 j ◦ uj+P m i6=j,i=1βiju 2 i ◦ uj= 0, uj > 0, uj ∈ RN, for j = 1, . . . , m, (3)

where uj∈ RN denotes the approximation of φj(z), for j = 1, . . . , m. Here A ∈ RN ×N

is the standard central finite difference discretization matrix of the operator [4 − V (z)] with the homogeneous Dirichlet boundary conditions on a finite domain Ω ⊆ Rn, with uj|∂Ω = 0. In addition, A is an irreducible and symmetric negative definite matrix.

The size of N depends on the approximation domain and grid sizes. For example, if a uniform grid size h is applied over a two-dimensional finite domain [−d, d] × [−d, d] and (2d

h − 1) is a positive integer, we have a N = ( 2d

h − 1)

2_{. For u = (u}

1, . . . , uN)>,

v = (v1, . . . , vN)> ∈ RN, u ◦ v = (u1v1, . . . , uNvN)> denotes the Hadamard product of

u and v and u
r _{= u ◦ · · · ◦ u denotes the r-time Hadamard product of u. Note that}

in the presence of strong periodic trapping potentials, the NLS equations (2) can be approximated by the DNLS equations.

The NLS equations (1) model a physical phenomenon in nonlinear optics [1], where the solution Φj denotes the j-th component of the beam in Kerr-like photorefractive

media. The positive constant µj measures the self-focusing in the j-th component of the

beam; and λj is the chemical potential. The coupling coefficient βij is the interaction

between the i-th and j-th components of the beam. For βij > 0, the interaction is

attractive; otherwise, the interaction is repulsive. Furthermore, Equations (3) describe a large class of discrete nonlinear systems such as optical fibers [11, 12], small molecules such as benzene [13], and, more recently, dilute Bose-Einstein condensates trapped in a multiwell periodic potential [3, 6, 23, 24].

The problem described in (2) has been considered in several cases. For n = 1, i.e. in one spatial dimension, the system (2) is integrable. Many analytical and numerical results on solitary wave solutions of m-coupled NLS equations are well-studied in, e.g., [10, 14, 15, 16]. For n = 2 and m = 1, physical experiments in [21] 2-dimensional photorefractive screening solutions and a 2-dimensional self-trapped beam. It is natural to believe that there are 2-dimensional m-component (m ≥ 2) solitons and self-trapped beams. A general theorem on the existence of high dimensional m-component solitons was first proved in [20]. In that paper, the authors show that the signs of the coupling

**Algorithm 1 ** **Algorithm 2 ** **Algorithm 3 **
Set*m*=3, ( ) 0,*V***z** = and

*ij* *ij*

*b* =*bd* in Eq. (3).

Set* b = and use Algorithm 1 *0
to solve the 1-comp. DNLS Eqs

**Use Algorithm 2 to track **
solution curves of the 1-comp.

DNLS Eqs with different
continuation parametersb &*d ij*

Output computed solutions No

Yes Input initial guess

Compute the next iterate

Convergence

Output solution _{Yes }

No Input initial solution

Track solution until bifurcates

Keep tracking

Output computed solutions Detect the next minimal

energy solution curve

Figure 1: Flowcharts of the algorithms proposed in this article. (a) Algorithm 1 (Fixed Point Iteration) is used to find energy minimizers of the 1-component decoupled DNLS equation. (b) Algorithm 2 (Minimal Energy Tracking Continuation Method) is used to track solution curves associated with minimal energy. (c) Algorithm 3 (Continuation Parameter Switching Method) is used to find the non-radially symmetric ground state predicted by [20].

coefficients βij’s are crucial for the existence of ground state solutions. For m = 3, when

all the βij’s are positive, there exists a ground state solution, which is radially symmetric.

The motivations of this study are twofold. First, while the search of minimal energy solutions of the DNLS equations is of interest, non-linearity of the equations of the solutions prevents existing numerical schemes from being complete. In particular, while continuation and homotopy type schemes provide feasible ways for computing various solutions of the DNLS equations; it remains an challenge to assert the minimality of the computed solutions. Another motivation arises in a theoretical prediction reported recently. Lin and Wei [20] consider the case V (z) = 0 and m = 3 with one repulsive and two attractive interactions. They show that if the coupling coefficients |βij| 1 and

the ground state solution exists, then the ground state cannot be radially symmetric. Computation of such non-radially symmetric energy minimizer can hardly be achieved by using continuation methods straightforwardly. The solution is actually “hidden” in the set of solutions and satiable routes leading to the target solution need to be identified. In this article, we develop numerical schemes for computing minimal energy solutions of the m-coupled DNLS equations. The main contributions of this article are highlighted in the following items and the flowcharts in Figure 1.

• We develop and prove a globally convergent fixed point iteration method (Figure 1a) for computing the energy minimizers of the 1-component decoupled DNLS equation. In the case that the decoupled DNLS system has multiple minimizers, Algorithm 1 converges to a local minimizer that depends on the initial solution. However, if the decoupled DNLS system has a unique (and thus global) minimizer, which is true

when V (z) = 0 [19], Algorithm 1 is expected to converge to the global minimizer. For the 3-coupled DNLS equations with V (z) = 0, these computed global energy minimizers are also used as the initial point of a continuation method for tracking other minimal energy solutions (Figure 1c).

• We propose a scheme for tracking minimal energy solutions by a continuation method (Figure 1b). To detect solution curves with minimal energy, we first show that the local minimal energy properties of the solutions remain unchanged be-tween solution curve bifurcation points under mild assumptions. Consequently, whenever we detect a bifurcation point, we can probe each of the available bifur-cation branches by checking the solution properties one step ahead to determine the target branch associated with the minimal energy solutions. By following this chosen branch, we can then track minimal energy solutions along a suitable solution curve.

• By integrating the above two schemes, we further develop a parameter-switching scheme (Figure 1c) to qualitatively find the non-radially symmetric ground state in a 3-component DNLS equations system with two attractive interactions and one repulsive interaction and V (z) = 0. The solution is predicted by [20] in an asymptotic setting. The proposed scheme confirms the existence of the solution numerically and visually.

It is worth noting that, in some other models, the DNLS equations is equipped with an extra normalization constraint on uj. The induced problem consequently becomes an

eigenvalue problem. Numerical schemes such as [4, 5, 7, 8, 9] have been developed for solving such eigenvalue problems. However, they are not suitable for the target problem considered here.

To define the ground state solution to m-component cases, we first set the Nehari’s
manifold
Nm= {φ = (φ1, ..., φm) ∈ (H1(Ω))m| φj≥ 0, φj≡/ 0 and
Z
Ω
|∇φj|2+
Z
Ω
V (z)φ2_{j}+ λj
Z
Ω
φ2_{j} = µj
Z
Ω
φ4_{j} +
m
X
i6=j,i=1
βij
Z
Ω
φ2_{i}φ2_{j}, j = 1, ..., m },

and the energy functional

E(φ) =
m
X
j=1
R
Ω|∇φj|
2_{+}R
ΩV (z)φ
2
j+ λj
R
Ωφ
2
j+
µj
2
R
Ωφ
4
j
2 +
!
−1
4
m
X
i,j=1,i6=j
βij
Z
Ω
φ2_{i}φ2_{j},
(4)
where φ = (φ1, . . . , φm) ∈ (H1(Ω))m. Then, we consider the minimization problem

inf

φ∈Nm

E(φ). If φ = (φ1, . . . , φm) ∈ Nm has the following properties:

(i) φj > 0 for all j and φ satisfy (2);

(ii) E(φ) ≤ E(ψ) for any other solution ψ of (2),

then φ is called a ground state solution of (2). In the case of m = 1 (1-component), the ground state solution of (2) can be obtained from the minimization problem

inf
φ ≥ 0
φ ∈ H1 (Ω)
R
Ω|∇φ|
2_{+}R
ΩV (z)φ
2_{+ λ}R
Ωφ
2
R
Ωφ
41/2 , (5)
up to a suitable normalization.

We define the corresponding energy functional for the m-coupled DNLS equations (3) as follows. The energy functional E(φ) in (4) becomes

E(u) = m X j=1 −u> jAuj+ λju>juj− µj 2u 2> j u 2 j 2 ! −1 4 m X i,j=1,i6=j βiju 2> i u 2 j , (6)

where the vector u = (u>_{1}, . . . , u>_{m})>∈ RN m _{is in}

Nm= {(u>1, . . . , u>m)> ∈ R N m | uj ≥ 0, uj 6= 0 and − u>jAuj+ λju>j uj = µju j2>u 2 j + m X i6=j,i=1 βiju i2>u 2 j , j = 1, . . . , m }. (7)

Throughout this paper, we use bold face letters or symbols to denote matrices or vectors. For u = (u1, . . . , uN)>, [[u]] := diag(u) denotes the diagonal matrix of u and

kuk4 = (u 2>u 2)1/4. For A ∈ RN ×N, A > 0 (≥ 0) denotes a positive (nonnegative)

matrix with positive (nonnegative) entries, A 0 (with A> = A) denotes a symmetric positive definite matrix, σ(A) denotes the spectrum of A, and N (A) and R(A) denote the null and range spaces of A, respectively.

This paper is organized as follows. In Section 2, we develop an iterative method to compute the positive ground state solution of 1-component DNLS equations and show that the iterative method is globally convergent. In Section 3, we discuss how we can track bifurcation branches in a continuation method to obtain local minimal energy solutions of the m-coupled DNLS equations. In Section 4, we propose a parameter-switching scheme to compute the non-radially symmetric ground state for the 3-component DNLS. Conclusion of the paper is given in Section 5.

2. Local and global energy minimizers of decoupled DNLS equations

We begin our discussion of numerical schemes and the corresponding analysis for solving 1-component decoupled DNLS equations. The solution of the decoupled equation is used as the initial solution of the m-coupled DNLS equations that m > 1. A fixed point iteration method is developed and analyzed for solving this 1-component DNLS equation. The computed solutions are local or global minimizers of the corresponding energy functional according to circumstances.

The 1-component DNLS equation is given by

Au − λu + µu
2 _{◦ u = 0,}

u > 0, u ∈ RN, (8)

where λ and µ are positive constants. The variational problem corresponding to (5) can be formulated as

inf

u≥0E(u),b (9a)

where

b

E(u) = −u

>_{Au + λu}>_{u}

(u
2>_{u}
2_{)}1/2 . (9b)

Formulation (9) is equivalent to minimize (6) subject to the constraint (7) with m = 1.
Consequently, any solution u ∈ RN _{of (6) is a local minimum or a saddle point on the}

Nehari’s manifold (7) with m = 1.

Next, we develop a fixed point iteration method for finding minimizers of (9) based on the following observations. The matrix A in (8) is generically diagonally dominant with nonnegative off-diagonal entries. That is, −A is an irreducible M-matrix. Let

¯

A = λI − A. (10)

Then ¯A is an irreducible M-matrix because λ > 0. It follows that ¯A−1is positive definite with positive entries (i.e., ¯A−1 0 and ¯A−1> 0). We define the set

M =_{u ∈ R}N_{| kuk}

4= 1, u ≥ 0 , (11)

It is easy to verify that if u ∈ M, then ¯

A−1u = (λI − A)−1u > 0. (12)

We now define a map f : M → M by
f (u) =
¯
A−1u
3
k ¯A−1_{u}
3_{k}
4
. (13)

Since the map f is well-defined by (12) and (13), we can use f to define the fixed point iteration ui+1= f (ui), as in the following algorithm.

Algorithm 1. [Fixed Point Iteration]

(i) Let ¯_{A ∈ R}N ×N, u0> 0 with ku0k4= 1, and i = 0;

(ii) Solve the linear system

¯

Aui+1 = u 3 i .

Compute ui+1= ui+1/kui+1k4.

(iii) If convergent, then u∗← ui+1, stop; else i ← i + 1, go to (ii).

Now, we analyze the convergence behavior of Algorithm 1. Detailed proofs of the theorems are given in the appendix. First, we show that the 1-component DNLS equation (8) has a solution ¯u(µ) and the solution can be computed by using the fixed point of f .

Theorem 1. The map f : M → M given in (13) has a fixed point u∗ in

◦

M. Further-more, the vector

¯
u(µ) = 1
µ1/2k ¯A
−1_{u}∗
3_{k}−1/2
4 u
∗_{∈ N}
1 (14)

solves the 1-component DNLS equation (8).

Theorem 1 suggests solving the 1-component DNLS equation (8) by Algorithm 1. The following theorems further discuss how the solution sequence generated by Algorithm 1 converges to a fixed point of f . In Theorem 2 below, we first show that the energy sequence corresponding to the iterates is decreasing, and therefore a subsequence of the iterates converges to a fixed point in M of f . In Theorem 5 below, by making a mild assumption, we further show that the whole sequence {ui}∞i=0 generated by Algorithm 1

converges to u∗∈ M globally.

Theorem 2. (i) If u ∈ M and v = f (u), then bE(v) ≤ bE(u), where bE(·) is defined in (9b), and the equality holds if and only if u is a fixed point of f : M → M, i.e., f (u) = u. (ii) For a sequence {ui}∞i=0 generated by Algorithm 1, there exists a subsequence {uni}

∞
i=0
such that
lim
i→∞uni = u
∗_{,} _{(15)}

where u∗∈ M is a fixed point of the function f defined in (13).

The following corollary can be easily obtained by applying Theorem 2.

Corollary 3. If the minimization problem (9) has a unique global minimizer u∗ ∈ M, then there exists a neighborhood Ru∗ of u∗ such that the fixed point iteration converges

to u∗ for any initial vector u0 ∈ Ru∗. In addition, ¯u(µ), defined in (14), is a global

minimizer of (9).

Now, we discuss how the entire sequence generated by Algorithm 1 converges. To
do so, we first define the ¯A-norm of u by kuk_{A}¯ =

√

u>_{Au and introduce the following}¯
lemma. Note that the definition of ¯A-norm is well-defined, as ¯A is positive definite.
Lemma 4. Let {ui}∞i=0 be the sequence generated by Algorithm 1. We have

lim

i→∞kui+1− uikA¯ = 0. (16)

Theorem 5 (Existence of a globally convergent sequence). If u∗given in (15) is a strictly local minimum of (9), then the sequence {ui}∞i=0 generated by Algorithm 1

converges to u∗∈ M.

In Theorem 5, we have shown that if a limit point u∗ of {ui}∞i=0 is a strictly local

minimum of (9), then the sequence {ui}∞i=0generated by Algorithm 1 converges globally

to u∗∈M and u◦ ∗ _{satisfies}

(λI − A)u∗= τ u∗
3 _{where τ = k ¯}_{A}−1_{u}∗
3_{k}
4.

We can then compute ¯u(µ) by (14) to find the minimizers of the 1-component DNLS equation (8).

Remark 1. The 1-component NLS (i.e. Equation (2) with m = 1) may associate with multiple local minimizers. In such cases, finding the global minimizer is a challeng-ing task. However, in the case that V (z) = 0, the 1-component NLS equation has a unique global solution [19]. Consequently, the computed solution ¯u(µ) is expected to be the ground state solution of the 1-component DNLS equation (8), though in the absence of a rigorous proof. This computed ground state solution is then used in Section 4 as the initial of a continuation method for tracking the non-radially symmetric energy minimizer predicted in [20] that also assumes V (z) = 0.

Remark 2. Algorithm 1 is similar to the continuous normalized gradient flow method (CNGFM) [5] in the sense that both methods generate energy decreasing sequences with mass conservation. Actually, Algorithm 1 can be formulated as a special case of CNGFM as shown below. However, a difference between the two methods does exist. The corresponding energy in each iteration of Algorithm 1 is decreased with step size equals 1 and the energy diminishing property holds for any initial vector u0∈ M. In contrast, CNGFM does not guarantee the energy diminishing property,

except the virtual time step ∆t is sufficiently small. Such a requirement may slow down convergence rate in the scheme.

Now, we verify the claim given above. We first note that we consider problem (9) for the 1-component DNLS equations. Since bE(u) = bE(cu) for c > 0, it is natural to find the minimizers on a “unit sphere”. Here, we choose M, which is defined in (11), as the restriction set. Problem (9) is thus transformed to be the form:

min

u∈Mu

>_{Au,}_{¯} _{(17)}

where ¯A = λI − A. Let eE(u) = 1 2u

>_{Au −}_{¯} 1
2

√

u
2>_{u}
2 _{and consider the }

minimiza-tion problem

min

u∈ME(u).e (18)

One can easily verify that the problems (9), (17) and (18) have the same minimizer. In the following, we shall focus on problem (18). Taking gradient of eE(u), we have

∇ eE(u) = ¯Au − (u
2>_{u}
2_{)}−1/2_{u}
3_{.}

By applying the CNGFM, we have (

ut(t) = −∇ eE(u(t)), tn < t < tn+1= tn+ ∆t, n > 0,

u(tn+1) := u(t+n+1) = u(t −

n+1)/ku(t − n+1)k4.

(19)

By further discretizing (19), we have the following semi-implicit time discretization
scheme
(_{u}_{˜}
n+1−un
∆t = −( ¯A − I)˜un+1− un+ (u
2>
n u
n2)−1/2u
n3,
un+1 := ˜un+1/k˜un+1k4.
(20)
8

Taking ∆t = 1 and using the fact that kunk4= 1, (20) can be reduced to

( ˜

un+1 = ¯A−1u n3,

un+1 := ˜un+1/k˜un+1k4.

This is exactly the fixed point iteration (Algorithm 1) proposed in this section for solving minimizers of the 1-component DNLS equations. Consequently, we claim that the fixed point iteration is a special case of the CNGFM with a large time step ∆t = 1.

3. The Minimal Energy Tracking Continuation Method

Now, we focus on how to track minimal energy solutions in the framework of con-tinuation methods. After a brief introduction of concon-tinuation methods, we discuss the technique for tracking minimal energy solutions. At the end of this section, we integrate all the proposed ideas into a continuation method algorithm for tracking minimal energy solutions.

3.1. General Framework of Continuation Methods

We briefly introduce a general framework of a continuation method for the m-coupled DNLS equation (3). For a detailed discussion of the continuation scheme, see [2, 8, 17, 18], for example.

Denote the continuation parameter by β ≥ 0, and rewrite the m-coupled DNLS equation (3) as

G(x, β) = 0, (21)

where x = (u>1, . . . , u>m)> ∈ RmN and G = (G1, . . . , Gm) : RmN × R → RmN is a

smooth mapping with

Gj(x, β) = Auj− λj(β)uj+ µj(β)u 2 j ◦ uj+ m X i6=j,i=1 βij(β)u 2 i ◦ uj, (22)

for j = 1, . . . , m. The parameters λj, µj, and βij in (22) may depend on β. One example

is to fix λj and µj and set βij(β) = ββij. Furthermore, we define the solution curve of

(21) as

C =y(s) = (x(s)>_{, β(s))}>

| G(y(s)) = 0, s ∈ R , (23) assuming that a parametrization via arc-length s is available.

Two main components of a continuation method are to follow the solution curve and to test bifurcation points. To follow the solution curve, we use the prediction-correction process. Suppose yi(s) = (x>i (s), βi(s))> ∈ RmN +1 is a solution lying (approximately)

on the solution curve C. Starting from the point yi(s), standard continuation methods

usually take the tangent vector of the solution curve at yi(s) as the prediction vector.

In particular, the tangent vector can be computed by solving the linear system DG(y(s)) ˙y(s) = 0,

which is obtained by differentiating equation (21) with respect to s. Here ˙y(s) = ( ˙x(s)>, ˙β(s))> is a tangent vector of C at y(s), and

DG(y(s)) = [Gx(y(s)), Gβ(y(s))] ∈ RmN ×(mN +1) (24)

denotes the Jacobian matrix of G at y(s). To track the solution curve described in (23), we first find the Euler predictor

yi+1,1= yi+ hiy˙i, (25)

where hi> 0 is the step length and ˙yi is the normalized tangent vector at yi. Newton’s

method is then used to find the corrector to improve the accuracy of yi+1,1. More

precisely, for the correction vector δl, the iteration

yi+1,l+1= yi+1,l+ δl (26)

is computed for l = 1, 2, . . . until a convergence criterion is satisfied for l = l∞ and we

take

yi+1 = yi+1,l∞ (27)

as a new approximate solution on the solution curve C.

To test bifurcations, we rely on Theorem 6. The theorem gives conditions under which a point on the solution curve is a bifurcation point. The theoretical and numerical details for detecting bifurcation points of the solution curve C and for tracing the bifurcation branches can be found in [8, 18]. Note that from Theorem 6, we see that a bifurcation point occurs when the matrix Gx is singular.

Theorem 6. (Bifurcation Test [17]) Let C in (23) be a smooth curve of (21) parameter-ized by s. Suppose det(Gx(y(s))) changes sign at s∗. Then y(s∗) is a bifurcation point

of (21).

We have discussed how we follow the solution curve and detect bifurcation points in the continuation method. In the next section, we focus on how we may determine the bifurcation branches associated with minimal energy solutions.

3.2. Minimal Energy Tracking

In this subsection, we discuss how to determine whether a solution to the DNLS equation (3) is a local minimum of E(x) on Nehari’s manifold Nm. The main result is

that each bifurcation point of the solution curve C coincides with a bifurcation of critical points for E(x) on Nm, as will be shown in Theorem 10.

First, we define the necessary notation and discuss how to verify local minimum of E(x) on Nmby applying standard optimization techniques to the optimization problem

inf

x∈Nm

E(x), (28)

where the energy functional E(x) and Nehari’s manifold Nm are defined by (6) and (7),

respectively. We define the Lagrangian function of the optimization problem (28) as L(x, ν) = E(x) − m X j=1 νjgj(x), (29) 10

where gj(x) = u>jAuj−λjuj>uj+µju j2>u 2 j + Pm i6=j,i=1βiju 2> i u 2 j , and ν = (ν1, . . . , νm)> ∈

Rmare the Lagrange multipliers. Furthermore, the discretized Nehari’s manifold Nm

de-fined in (7) can be written as
Nm=(u>1, . . . , u
>
m)
>
∈ RmN_{|u}
j ≥ 0, uj 6= 0 and gj(x) = 0, j = 1, . . . , m .

The total derivative of the function g = (g1, . . . , gm)>

∇g(x) =
∇g1(x)
..
.
∇gm(x)
∈ R
m×mN_{,} _{(30)}

where ∇gj(x) is the row vector given by

(∇gj(x))_{i}=
2βijui◦ u
j2, i 6= j,
2(Auj− λjuj+ 2µju
j2 ◦ uj+
m
X
i6=j,i=1
βiju
i2 ◦ uj), i = j. (31)

The following theorem gives a test to determine whether a point x ∈ RmN _{is a local}

minimum of E(x) on Nmor not.

Theorem 7. (Nocedal-Wright [22]) Let x be a point in RmN_{. Suppose that the }

Karush-Kuhn-Tucker (KKT) conditions

∇xL(x, ν) = 0 and x ∈ Nm (32)

are satisfied for a certain ν ∈ Rm_{. Suppose also that}

w>∇2

xxL(x, ν)w > 0, for all w ∈ C(x, ν), w 6= 0, (33)

where C(x, ν) is the null space of ∇g(x) in (30). Then x is a strict local minimum solution of E(x) on Nm.

Furthermore, we derive the relations between the DNLS equations (21) and the opti-mization problem (28) in the following remarks.

1. All solution points lying on the solution curve C of DNLS equation satisfy the KKT conditions (32) with ν = 0.

From (29) and Theorem 7, a point x that satisfies the KKT conditions is the solution of

∇xE(x)>− ∇g(x)>ν = 0 and x ∈ Nm,

where ∇g(x) is given in (30) and ∇xE(x) = (∇u1E(x), . . . , ∇umE(x)) ∈ R

mN _{with}
∇ujE(x)
>_{= Au}
j− λjuj+ µju
j2 ◦ uj+
m
X
i6=j,i=1
βiju
i2 ◦ uj.

Since the equation ∇xE(x) = 0 is actually the DNLS equation (21), we see that if

y = (x>, β)>is a solution of the DNLS equation (21), then x ∈ Nmand ∇xL(x, 0) =

0. That is, all points on the solution curve C of the DNLS equation satisfy the KKT conditions (32) with ν = 0.

2. We define the projected Hessian matrix H(y(s∗)), to be used later for testing positive definiteness. Let y = (x>, β)> be a point on C. We define the projected Hessian matrix at y to be

H(y) ≡ P(x)>∇2xxL(x, 0)P(x) = P(x)>Gx(y)P(x), (34)

where Gx(y) is given in (24). Here P(x) ∈ RmN ×` is the matrix whose columns

form an orthonormal basis of the null space of ∇g(x). In other words,

[∇g(x)]P(x) = 0, (35)

and if ∇g(x) is of full row rank, then ` = mN − m.

3. Applying the above two remarks and considering the path y(s) that describes the solution curve C, observe that if y(s) is a local minimum of (28) for s < s∗ and a saddle point for s > s∗, then from (33) in Theorem 7, the projected Hessian matrix H(y(s∗)) is singular.

We have presented the relations between the DNLS equations (21) and the optimiza-tion problem (28). Now, we develop the relaoptimiza-tions between the soluoptimiza-tion curve bifurcaoptimiza-tion and the critical point bifurcation by applying Theorem 6 (for the solution curve bifur-cation test) and Theorem 7 (for the optimality test). Specifically, in Lemma 8, we show that for each bifurcation point y∗= y(s∗) ∈ C, i.e., where the matrix Gx(y∗) is singular,

the projected Hessian matrix H(y∗) is singular. In Lemma 9, we show that the converse of Lemma 8 is true whenever the auxiliary matrix Σ(y) is invertible. Here

Σ(y) = (2βiju 2> i u 2 j ), (36) and we set βii = µi.

Lemma 8. Let y∗ be a bifurcation point of the DNLS equation (21) along the solution curve C. Then det(H(y∗)) = 0.

Lemma 9. Let y∗_{∈ C. If det(H(y}∗_{)) = 0 and Σ(y}∗_{) is invertible, then det(G}

x(y∗)) =

0.

The above two lemmas suggest the following result. Let y∗ ∈ C and Σ(y∗_{) be an}

invertible matrix, then det(H(y∗)) = 0 if and only if det(Gx(y∗)) = 0. In other words,

the bifurcation of the solution curves C and the bifurcation of the critical points for E(x) on Nmoccur at exactly the same points.

Combining Lemma 8 and Lemma 9, we can see that solutions lying on a so-called “solution segment” have the same critical point characteristics. Before giving the rigorous statement of the result in Theorem 10, we first define the solution segment. Let y0 =

y(s0) be any regular point (i.e., det(Gx(y0)) 6= 0) of the solution curve C. Let Γseg(y0) ⊆

C be a maximal connected set without bifurcation point and containing y0. That is, a

solution segment

Γseg(y0) = {y(s) ∈ C | det(Gx(y(¯s))) 6= 0 for ¯s between s0 and s}. (37)

Now, we state the theorem that can be used to choose suitable bifurcation branches that lead to minimal energy solutions while bifurcations occur.

Theorem 10. Let y0= (x>0, β0)> ∈ C with det(Gx(y0)) 6= 0 and suppose that Σ(y) is

invertible for each y ∈ Γseg(y0). If x0is a strict local minimum solution of E(x) on Nm,

then for each y = (x>, β)>∈ Γseg(y0), x is a strict local minimum solution of E(x) on

Nm.

In summary, Theorem 10 suggests that along the solution curve C tracked by a con-tinuation method, if the initial point is a local minimum of E(x) on Nm and the local

minimum becomes a saddle point somewhere along the solution curve, then it must meet a bifurcation point under some mild assumptions. Consequently, whenever we detect a bifurcation point, we can track each of the available bifurcation branches one step ahead and test the positivity of the corresponding projected Hessian matrices. According to Theorem 10, we can then determine which one is the local minimum energy branch. 3.3. The Overall Algorithm

Finally, we conclude this section by proposing the minimal energy tracking continu-ation method (METCM) in Algorithm 2.

Algorithm 2. [Minimal Energy Tracking Continuation Method (METCM)]. (i) [Initialization]

Determine an initial solution y(0) by letting β(0) = 0 (e.g., use Algorithm 1 to solve the resulting decoupled DNLS equations).

(ii) [Solution curve following] Iterate until bifurcation occurs.

(a) Compute the next solution y(s) by the predictor-corrector scheme described in (25), (26), and (27).

(b) Detect the bifurcation point by Theorem 6 and techniques described in [8, 18]. (iii) [Minimal energy solution curve detection]

(a) Track one step ahead for each of the available bifurcation branches. (b) Test the positiveness of the corresponding projected Hessian matrices.

(c) Pick one bifurcation branch with positive projected Hessian matrices as the next minimal energy solution curve to be followed.

(iv) Go to Step (ii) to follow the next solution curve or stop.

Note that in Step (i) of the algorithm, we may set βij = 0 and solve the decoupled

system DNLS equation (3) by Algorithm 1 and (14) to obtain the initial solution y(0) =
((¯u, ¯u, ¯u)>, 0)>for m = 3. In Step (iii-b), it is possible to identify more than one minimal
energy branch. We can track all these branches simultaneously in parallel computations.
To conclude this section, we give a simple example to illustrate the key components
discussed in this section. The main goal is to show how we can track energy minimizers
of E(x) = 0.5x4_{+ (1 − β)x}2

by the METCM. Here, x ∈ R and β is the continuation
parameter. By setting G(x, β) := E0(x) = 2x2_{+ 2(1 − β)x, we can see that (i) the}

solutions of G(x, β) = 0 are actually the critical points of E(x), (ii) the solution curve C of G(x, β) = 0 is the union of the trivial curve C0 = {(0, β)| β ∈ R} and the parabolic

curve Cp = {(±

√

β − 1, β)| β > 1}, and (iii) y∗ = (0, 1) ∈ C0 is a bifurcation point, as

C0 is a smooth curve with parameter β and Gx(0, β) = 2(1 − β) changes sign at β∗ = 1

(Theorem 6).

Starting from y0= (0, 0) ∈ C, we track the solution curve until we meet the

bifurca-tion point y∗ = (0, 1). As the point y0 = (0, 0) is a strict local minimum of E(x) with

det(Gx(y0)) 6= 0 and det(Gx(y∗)) = 0, Theorem 10 suggests that all the solutions on

the solution segment, as defined in (37), Γseg(y0) = {(0, β)| β < 1} are all strict local

minimizers. To detect minimal energy solution curve as suggested in Step (iii) of Algo-rithm 2, we take one step ahead on the three available bifurcation branches and test the points (√β − 1, β)), (−√β − 1, β)) and (0, β) for β = 1+ to determine the next minimal energy curve to be followed. By doing so, we can track the local minimum energy curve Cp rather than {(0, β)| β > 1}. It is worth noting that, a straightforward calculation can

verify that the local minimum curve of E(x) is Γseg(y0) ∪ Cp, which matches the results

obtained by the METCM exactly.

4. Minimal Energy Tracking for Non-Radially Symmetric Solutions

In this section, we demonstrate the capabilities of the METCM by finding non-trivial minimal energy solutions of a 3-coupled DNLS problem that has one repulsive and two attractive interactions and assumes V (z) = 0. Our main tools are threefold: the compu-tation of the ground states of the decoupled DNLS equations in Section 2, the minimal energy tracking continuation method in Section 3, and a parameter-switching scheme to be discussed below. By combining these techniques, we can find a 3-component non-radially symmetric energy minimizer while βij approaches zero. The existence of such

non-radially symmetric solution has been predicted by Lin and Wei theoretically [20] for m = 3 and V (z) = 0 in Equations (2). They show that with one repulsive and two attractive interactions, if the coupling coefficients |βij| 1, and the ground state

solution exists, then the ground state must be non-radially symmetric. Furthermore, the corresponding energy is smaller than the energy of the positive radially symmetric solution.

Using notations similar to those in [20], we consider the following 3-coupled DNLS equations by assuming λ1= λ2 = λ3 = µ1= µ2 = µ3= 1. We also rewrite βij = βδij

and assume δ12= δ21, δ13= δ31, and δ23= δ32,

Au1− u1+ u 3 1 + βδ21u 2 2 u1+ βδ31u 2 3 u1= 0, Au2− u2+ u 3 2 + βδ12u 2 1 u2+ βδ32u 2 3 u2= 0, Au3− u3+ u 33 + βδ13u1 2u3+ βδ23u 22u3= 0.

Note that the matrix A in Equations (38) is the discretization matrix of the Laplacian only as V (z) = 0.

To the best of our knowledge, such non-radially symmetric solutions have not been computed and visualized numerically. Simple straightforward numerical methods cannot lead to non-radially symmetric solutions. For example, only radially symmetric solutions are found for small βij’s if we simply follow the solution curve

Cβ=(x>, β)>| G(x, β) = 0 is given in (38) with

δ12= δ13= 1 and δ23= −1, for β ∈ R+} .

by starting from β = 0 [18].

### 1

### 1

### −

1###

###

2 3###

*3N*

###

### ( , , )

**u u u**

23
### δ

0.2### β

Figure 2: Illustration for the curves C1, C2 and C3in Steps (ii)-(iv) of Algorithm 3.

Now, we describe how we can find the non-radially symmetric solutions for the case of one repulsive and two attractive interactions in the 3-coupled DNLS Equations (38). We first propose a continuation parameter-switching scheme with an auxiliary illustration in Figure 2. Then we describe the motivations behind the scheme.

Algorithm 3. [Continuation Parameter Switching Method]

(i) Compute the initial solution (¯u, ¯u, ¯u) by Algorithm 1, where ¯u is the solution of the 1-component DNLS (or the decoupled DNLS for β = 0).

(ii) Let β be the continuation parameter. Use METCM to track the solution curve C1=(x>, β)>| G(x, β) = 0 is given in (38) with

δ12= δ13= δ23= 1, for 0 ≤ β ≤ 0.2}

from β = 0 to β = 0.2.

(iii) Let δ23 be the continuation parameter and fix β = 0.2. Use METCM to track the

solution curve

C2=(x>, δ23)>| G(x, δ23) = 0 is given in (38) with

β = 0.2, δ12= δ13= 1, for − 1 ≤ δ23≤ 1} .

from δ23= 1 to δ23= −1.

(iv) Let β be the continuation parameter. Use METCM to track the solution curve C3=(x>, β)>| G(x, β) = 0 is given in (38) with

δ12= δ13= 1 and δ23= −1, for β ∈ R} ,

from β = 0.2 to β ≈ 0+.

We implement Algorithm 3 on a square domain [−5, 5]×[−5, 5] with grid size h = 0.2. We plot the conceptual solution curves (with bifurcation points) and the corresponding

energy curves of C2 in Figures 3 and 4, respectively. We use the same curve styles in

these two figures to indicate the corresponding solution curves. Similarly, the bifurcation diagram and energy curves of C3 are shown in Figures 5 and 6, respectively. In the

solution curves (i.e., Figures 3 and 5) the corresponding nodal domains of three positive bound state solutions of certain segments of the solution curves are attached in triples near the solution curves. In each of the nodal domain triples, the left, middle and right figures are the density plots of u1, u2 and u3, respectively. As the solution formats

remain similar unless bifurcation occurs, only one representative nodal domain triple is shown for each of the solution curve segments. In Figures 4 and 6, the plots of the nodal domains (in the form of squared sums) overlap to show their relative positions. The corresponding triples and overlapping nodal domains are labeled by the same capital letters in different figures, e.g. Figures 3 and 4.

Now, we justify the continuation parameter-switching scheme and note some obser-vations from the numerical results shown in the figures.

1. In Step (ii) of Algorithm 3, we follow the solution curve C1 by increasing β from 0

to 0.2. As the initial solution (¯u, ¯u, ¯u) is the global (thus also local) minimal energy solution, and there is no bifurcation found in the interval 0 ≤ β ≤ 0.2, the states of the solutions are thus unchanged by Theorem 10. That is, all the intermediate solutions are all local minimal solutions corresponding to each of the β’s.

Furthermore, we anticipate that all these solutions corresponding to each of the β’s are ground state solutions, due to the following observation. In this setting, all the interactions are attractive. Consequently, the three components tend to gather together and concentrate at the center of the domain to achieve minimal energy. Such solution profiles are similar to the solution profile of the initial solution (¯u, ¯u, ¯u) for β = 0. It is thus reasonable that (¯u, ¯u, ¯u) is a good initial guess for the global minimal solution of a DNLS with a small positive β. By using the continuation method, we thus can track the global minimal solutions in C1.

2. The curve C2 acts as a “bridge” connecting the two settings in Step (ii) and Step

(iv). In particular, δ23 is changed from 1 (three attractive interactions) to −1 (one

repulsive and two attractive interactions). Figure 3 shows that there is only one bifurcation point in C2, where δ23 = −0.314. At this bifurcation point, METCM

suggests tracking either the upper or the lower bifurcation branch (the red curves) to retain minimal energy solutions. By tracking either one of the branches, non-radially symmetric solutions are observed. Furthermore, as shown in Figure 4, the solutions of these bifurcation branches have lower energies than the ones on the primal stalk of C2 (the blue curve).

3. In Step (iv) of Algorithm 3, we switch to the target setting, in which one repulsive and two attractive interactions are assumed. We use the non-radially symmetric solution obtained in the terminal point of C2, in which δ23= −1, as the initial guess

for tracking the curve C3. If we decrease β from 0.2 to 0+, the target non-radially

symmetric solutions are obtained while β approaches zero. As shown in Figure 5, no bifurcation occurs for 0 < β ≤ 0.2, so all the computed solutions remain minimal energy solutions. Consequently, we find the non-radially symmetric solution for one repulsive and two attractive interactions with small β.

In short, we have shown how the parameter-switching scheme can lead to the non-radially symmetric solution predicted in [20]. These non-non-radially symmetric positive

1 0 1
− −0.314
23
δ
*3N*
A
B
C
D

Figure 3: A bifurcation diagram of the solution curve for C2= {(x>, δ23)>| G(x, δ23) = 0 is given in

(38) with β = 0.2, δ12= δ13= 1, for −1 ≤ δ23≤ 1}.
1
− 0 1
17
16.5
16
15.5
15
14.5
23
δ
(
)
*E*
**x**
0.314
−
A
C
B, D

Figure 4: Energy curve of C2.

solutions are expected to be the ground state solutions of (38). This conjecture is based on the following observations. We start from the ground state solution for β = 0 and then track C1 until β = 0.2. We then track the path with lower energy solutions in C2,

when the only bifurcation occurs, and obtain non-radially symmetric solutions. As there is no other bifurcation point in C3, the value for β changes from 0.2 to 0+. The ground

state character of the initial solutions is preserved in the tracked non-radially symmetric solutions.

Finally, we make the following remarks on the numerical experiments.

Remark 3. We can also track the primal stalk of C2 in Step (iii) of Algorithm 3 until

δ23= −1 (i.e., the blue curve in Figure 3), and then track Cβby decreasing β = 0.2

to 0+. However, we would only find radially symmetric solutions, whose energies are higher than those of the non-radially symmetric solutions we have found. Remark 4. In Figure 7, we compare the energy curves of C3and Cβ. The energy curve of

C3is lower than that of Cβ. This result is obviously consistent with the consequence

reported in [20].

Remark 5. Another type of non-radially symmetric positive solution can be found by following the solution curve C3 and increasing the values of β from 0.2 to 1−. A

representative of such solutions is shown in Figure 5 for β ≈ 0.969. 17

0.969
β
*3N*
A
B
C

Figure 5: A bifurcation diagram of the solution curve for C3= {(x>, β)>| G(x, β) = 0 is given in (38)

with δ12= δ13= 1 and δ23= −1, for β ∈ R}.

0
11
12
14
16
18
0.969
β
(
)
*E*
**x** A _{B, C}

Figure 6: Energy curve of C3.

0.969
11
13
15
17
β
(
)
*E*
**x**
3
0
β

Figure 7: Compared with Cβ, the solution curve C3has lower energies.

5. Conclusion

This article focused on the use of continuation methods to solve the time-independent m-coupled discrete nonlinear Schr¨odinger equation. In particular, we propose a new algorithm that is capable of tracking the local minimal energy solutions along the solution curves. We have also shown how we may compute the ground states of the decoupled discrete nonlinear Schr¨odinger equation. By combining these two techniques with a parameter-switching scheme, we find non-radially symmetric minimal energy solutions for the case with one repulsive and two attractive interactions and small coupling coefficients. We believe our minimal energy tracking continuation method can be applied to other coupled elliptic partial differential equations, probably with suitable modifications. The method thus acts as a useful tool for exploring various steady-state solutions of the differential equations.

Acknowledgments

The authors are grateful to the anonymous referees for their valuable and helpful com-ments and suggestions. This work is partially supported by the National Science Council, the National Center for Theoretical Sciences, and the Taida Institute of Mathematical Sciences in Taiwan.

References

[1] N. Akhmediev and A. Ankiewicz. Partially coherent solitons on a finite background. Phys. Rev. Lett., 82:2661–2664, 1999.

[2] E. L. Allgother and K. Gerog. Numerical path following. In P. G. Ciarlet and J. L. Lions, editors, Handbook of Nemerical Analysis, volume V. Elsevier, 1997.

[3] D. Bambusi and A. Sacchetti. Exponential times in the one-dimensional Gross-Pitaevskii equation with multiple well potential. Communications in Mathematical Physics, 275(1):1–36, October 2007. [4] W.-Z. Bao. Ground states and dynamics of multi-component Bose-Einstein condensates. SIAM

Multiscale Modeling and Simulation, 2(2):210–236, 2004.

[5] W.-Z. Bao and Q. Du. Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. SIAM J. Sci. Comput., 25(5):1674–1697, 2004.

[6] F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio. Josephson junction arrays with Bose-Einstein condensates. Science, 293:843, 2001. [7] S.-L. Chang, C.-S. Chien, and B.-W. Jeng. Liapunov-schmidt reduction and continuation for

non-linear Schr¨odinger equations. SIAM J. Sci. Comput., 27:729–755, 2007.

[8] S.-M. Chang, Y.-C. Kuo, W.-W. Lin, and S.-F. Shieh. A continuation BSOR-Lanczos-Galerkin method for positive bound states of a multi-component Bose-Einstein condensate. J. Comp. Phys., 210:439–458, 2005.

[9] S.-M. Chang, W.-W. Lin, and S.-F. Shieh. Gauss-Seidel-type methods for energy states of a multi-component Bose-Einstein condensate. J. Comp. Phy., 202:367–390, 2005.

[10] K.W. Chow. Periodic solutions for a system of four coupled nonlinear Schr¨odinger equations. Phys. Lett. A, 285:319–326, 2001.

[11] D. N. Christodoulides and R. I. Joseph. Discrete self-focusing in nonlinear arrays of coupled waveg-uides. Opt. Lett., 13(9):794, 1988.

[12] D. N. Christodoulides, F. Lederer, and Y. Silberberg. Discretizing light behaviour in linear and nonlinear waveguide lattices. Nature, 424:817–823, 2003.

[13] A. S. Davydov. The theory of contraction of proteins under their excitation. Journal of Theoretical Biology, 38(3):559–569, March 1973.

[14] F. T. Hioe. Solitary waves for n coupled nonlinear Schr¨odinger equations. Phys. Rev. Lett., 82:1152– 1155, 1999.

[15] F. T. Hioe and T. S. Salter. Special set and solutions of coupled nonlinear Schr¨odinger equations. J. Phys. A: Math. Gen., 35:8913–8928, 2002.

[16] T. Kanna and M. Lakshmanan. Exact soliton solutions, shape changing collisions, and partially coherent solitons in coupled nonlinear Schr¨odinger equations. Phys. Rev. Lett., 86:5043–5046, 2001. [17] H. B. Keller. Lectures on Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin,

1987.

[18] Y.-C. Kuo, W.-W. Lin, S.-F. Shieh, and W. Wang. A hyperplane-constrained continuation method for near singularity in coupled nonlinear Schr¨odinger equations. Preprint, 2008.

[19] M. K. Kwong. Uniqueness of positive solutions of ∆u − u + up_{= 0 in R}n_{. Arch. Rational Mech.}

Anal, 105(3):243–266, 1989.

[20] T.-C. Lin and J. Wei. Ground state of n coupled nonlinear Schr¨odinger equations in Rn_{, n ≤ 3.}

Commun. Math. Phys., 255:629–653, 2005.

[21] M. Mitchell, Z. Chen, M. Shih, and M. Segev. Self-trapping of partially spatially incoherent light. Phys. Rev. Lett., 77:490–493, 1996.

[22] J. Nocedal and S. Wright. Numerical Optimization. Springer, 2nd edition, 2006.

[23] A. Trombettoni and A. Smerzi. Discrete solitons and breathers with dilute Bose-Einstein conden-sates. Phys. Rev. Lett., 86:2353–2356, 2001.

[24] A. Trombettoni, A. Smerzi, and A. R. Bishop. Discrete nonlinear Schr¨odinger equation with defects. Phys. Rev. E, 67(1):016607–, January 2003.

Appendix

Proof of Theorem 1. Equations (12) and (13) imply that f is continuous on M. By the definition of (11), M is homeomorphic to an (N − 1)-dimensional standard simplex, which is convex and compact. Applying the Schauder fixed point theorem to f , we can see that there is a point u∗∈ M satisfying

f (u∗) = u∗. (39)

The existence of the fixed point u∗∈M follows from the fact that the function f in (13)◦ maps M into

◦

M. From (10), (13) and (39), we have
k ¯A−1u∗
3_{k}−1
4 u
∗
3 _{= (λI − A)u}∗_{.} _{(40)}
Multiplying (40) by 1
µ1/2k ¯A−1u∗
3_{k}−1/2

4 from the left and setting

¯
u(µ) = 1
µ1/2k ¯A
−1_{u}∗
3
k−1/2_{4} u∗,
we obtain

A¯u(µ) − λ¯u(µ) + µ¯u(µ)
3 _{= 0.}

It is easy to verify that ¯u(µ) belongs to N1. This completes the proof.

Proof of Theorem 2. (i) Since u, v ∈ M, we have kuk4= 1, kvk4= 1 and

b

E(v) = v>(λI − A)v. (41)

Substituting v = f (u) =
¯
A−1u
3
k ¯A−1_{u}
3_{k}
4
into (41), we get
b

E(v) = v>(λI − A)v = 1
k ¯A−1_{u}
3_{k}

4

v>u
3_{.}

Letting c = 1

k ¯A−1_{u}3
_{k}_{4} and applying the H¨older inequality (|x

>_{y| ≤ kxk}

pkykq where 1

p + 1

q = 1 and p > 1) with p = 4, q = 4/3, x = v and y = u

3_{, we obtain}

v>(λI − A)v ≤ ckvk4kuk34= c. (42)

Since kuk4= 1, we have

c = u

>_{A ¯}¯_{A}−1_{u}
3

k ¯A−1_{u}
3_{k}
4

= u>(λI − A)v. (43)

Since λI − A is positive definite, it has the Cholesky factorization (λI − A) = L>L. Applying the Cauchy-Schwarz inequality to (43), we obtain

c = u>L>Lv ≤ √

u>_{L}>_{Lu ·}√_{v}>_{L}>_{Lv}

= q

u>_{(λI − A)u ·}q_{v}>_{(λI − A)v.} _{(44)}

From (42) and (44), it follows that v>(λI − A)v ≤ c ≤

q

u>_{(λI − A)u ·}q_{v}>_{(λI − A)v} _{(45)}

and therefore

q

v>_{(λI − A)v ≤}q_{u}>_{(λI − A)u.}

Using the fact that kvk4= kuk4= 1, we have

b
E(v) =v
>_{(λI − A)v}
kvk2
4
≤ u
>_{(λI − A)u}
kuk2
4
= bE(u). (46)

The equality in (46) holds if and only if the inequalities in (42) and (44) become equalities. Furthermore, both inequalities in (42) and (44) hold if and only if the vectors v and u are linearly dependent, i.e., v = au for some a ∈ R. Since v > 0, u > 0 with kvk4= kuk4= 1, we have v = u. Hence, the equality in (46) holds if and only if u is a

fixed point of f .

(ii) Since the sequence {ui}∞i=0⊂ M is bounded, there exists a convergent subsequence

{uni} ∞

i=0 and a point u∗∈ M such that

lim

i→∞uni = u
∗_{.}

Consequently, we have lim

i→∞E(ub ni) = bE(u

∗_{) and lim}

i→∞E(f (ub ni)) = bE(f (u

∗_{)),} _{(47)}

as f and bE are continuous. Furthermore, since the cost function bE(·) in (9b) is continuous
on the compact set M, the function b_{E(·) : M → R}+ attains its minimum value on M.

From part (i) of this theorem, it can be easily seen that the sequence { bE(ui)}∞i=1converges

to a certain positive number bE∗. That is, lim

i→∞E(ub i) = bE

∗_{.} _{(48)}

By Equations (47), (48), and the fact that { bE(uni)} ∞

i=0 and { bE(f (uni))} ∞

i=0 are

sub-sequences of { bE(ui)}∞i=0, we see that the three sequences converge to the same value.

Consequently, we have

b

E(f (u∗)) = bE(u∗).
By part (i) of this theorem, we conclude that f (u∗) = u∗_{. }

Proof of Lemma 4. By definition,

kui+1− uik2_{A}¯ = (ui+1− ui)>A(u¯ i+1− ui)

= ui+1> Au¯ i+1+ u>i Au¯ i− 2u>i+1Au¯ i

= bE(ui+1) + bE(ui) − 2u>i+1Au¯ i. (49)

From (42), (43) and (45), we have

b

E(ui+1) ≤ u>i+1Au¯ i≤

q b E(ui+1) q b E(ui). (50)

Furthermore, by (49) and (50), it follows that
kui+1− uik_{A}¯ ≤

q b

E(ui) − bE(ui+1), (51)

or equivalently, limi→∞kui+1− uik_{A}¯ = 0.

Proof of Theorem 5. Since u∗ is a strictly local minimum of the optimization problem
(9), the Hessian matrix H(u∗) of bE(u) is positive definite. Therefore, there is a δ > 0
such that H(u) is positive definite, i.e., bE(u) is convex, for u ∈ M and ku − u∗_{k}_{¯}

A< δ.

For any positive number 0 < ε < δ/2, we let
b
Eε= min
ku−u∗_{k}_{¯}
A=ε
b
E(u) > bE∗, (52)

where bE∗is given by (48), and define
B(u∗, bEε) =
n
u ∈ M
ku − u
∗_{k}
¯
A< ε, bE(u) < bEε
o
. (53)

From (15) and (16), there exists N0∈ N such that

unj ∈ B(u
∗_{, b}_{E}

ε) and kui+1− uikA¯ < ε for nj, i > N0. (54)

Since 2ε < δ, if ui ∈ B(u∗, bEε) and kui+1− uik_{A}¯ < ε, then kui+1− u∗k_{A}¯ < 2ε. On

the other hand, using the fact that bE(ui+1) ≤ bE(ui) < bEε and bE(u) is convex on

ku − u∗_{k}
¯

A< δ it holds that ui+1∈ B(u∗, bEε). Thus, we have

kui− u∗k_{A}¯ < ε for all i > N0.

This completes the proof.

Proof of Lemma 8. Since y∗ = (x∗>, β∗)> is a bifurcation point, the matrix Gx(y∗)

is singular. Thus there exists a nonzero vector z = (z>1, . . . , z>m)> ∈ RmN such that

Gx(y∗)z = 0. (55)

Now, we claim that ∇g(x∗)z = 0. Since y∗ is a solution of DNLS equation (21) and
x∗= (u∗>_{1} , . . . , u∗>_{m})>, the vector ∇gj(x∗) in (31) can be written as

∇gj(x∗)>=
2β1ju∗1◦ u
∗
2
j
..
.
2β(j−1)ju∗j−1◦ u
∗
2
j
2µju∗j ◦ u
∗
2
j
2β(j+1)ju∗j+1◦ u
∗
2
j
..
.
2βmju∗m◦ u
∗
2
j
∈ RmN, (56)
for j = 1, 2, . . . , m. Let
U(x∗) =
u∗1 0 · · · 0
0 u∗
2 . .. 0
..
. ... . .. 0
0 · · · u∗_{m}
∈ RmN ×m. (57)
A calculation leads to
U(x∗)>Gx(y∗) = ∇g(x∗). (58)

From (55) and (58), it is easy to verify that ∇g(x∗_{)z = 0, that is, there exists a nonzero}

vector ¯_{z ∈ R}` _{such that z = P(x}∗_{)¯}_{z. Hence, from (34), we obtain H(y}∗_{)¯}_{z = 0, i.e.,}

det(H(y∗_{)) = 0. }

Proof of Lemma 9. First, we note that from (36), (57), (58), and (59), we have

U(x)>Gx(y)U(x) = Σ(y), (59)

for each y = (x>, β)>∈ C.

Since y∗ = (x∗>, β∗)> ∈ C and det(H(y∗_{)) = 0, there exists a nonzero vector ¯}_{z ∈}

RmN −m such that

H(y∗)¯z = 0. (60)

Let z = P(x∗)¯_{z ∈ R}mN_{. Note that z is a nonzero vector, since P(x}∗_{) is of full column}

rank. Now we claim that Gx(y∗)z = 0. From (59) and (58), it follows that

U(x∗)>Gx(y∗)z = ∇g(x∗)z = 0. (61)

On the other hand, from (34) and (60) we have

P(x∗)>Gx(y∗)z = 0. (62)

Combing (61) and (62) gives
P(x∗)>
U(x∗)>
Gx(y∗)z = 0. (63)
We calculate that
P(x∗)>
∇g(x∗_{)}
[P(x∗), U(x∗)] =
Im(N −1) P(x∗)>U(x∗)
0 Σ(x∗)
. (64)

Here, the (1, 1) entry in (64) follows from the fact that the column vectors of P(x∗) are orthonormal, the (2, 1) entry follows from equation (35) and the (2, 2) entry follows from equations (58) and (59). Since Σ(y∗) is invertible, it follows that the mN × mN matrix [P(x∗), U(x∗)] is invertible. From (63), we obtain Gx(y∗)z = 0, that is, det(Gx(y∗)) =

0.