A continuation method for positive bound states of
m-coupled nonlinear Schr¨odinger equations
Yuen-Cheng Kuo
∗Wen-Wei Lin
†Shih-Feng Shieh
‡Weichung Wang
∗Abstract
We develop a stable continuation method for the computation of positive bound states of time-independent, m-coupled nonlinear algebra equation (DNLS) which describe a m-coupled nonlinear Schr¨odinger equation. The solution curve with respect to some parameter of the DNLS is then followed by the proposed method. For a one-component DNLS, we develop a iterative method for compu-tation of positive ground state solution. For a m-coupled DNLS, we prove that m radially symmetric positive bound states will bifurcate into m different positive bound states at a finite repulsive inter-component scattering length. Numerical results show that various positive bound states of a three-component DNLS are solved efficiently and reliably by the continuation method.
1
Introduction
We consider the m-coupled nonlinear Schr¨odinger equations
△φj − λjφj+ µj|φj|2φj+ X i6=j βij|φi|2φj = 0, in Rn, (1.1a) φj > 0 in Rn, j = 1, . . . , m, (1.1b) φj(x)→ 0, as |x| → ∞, (1.1c)
∗Department of Applied Mathematics, National University of Kaohsiung, Kaohsiung, 811, Taiwan
([email protected], [email protected])
†Department of Mathematics, National Tsinghua University, Hsinchu, 300, Taiwan (
‡Department of Mathematics, National Taiwan Normal University, Taipei, 11677, Taiwan
where λj, µj > 0 are positive constants, n ≤ 3, and βij = βji for i 6= j are coupling
constants. As βij = βji > 0, the interaction between the i-th and j-th component of
the beam is attractive, otherwise the interaction is repulsive.
In the one-component case (m = 1), the solution of (1.1) can be obtained through the following minimization
inf φ ≥ 0 φ ∈ H1(Rn ) R Rn|∇φ|2 + λ R Rnφ2 R Rnφ4 1/2 . (1.2)
An equivalent formulation is to consider the following minimization problem inf φ∈N1 E(φ) (1.3a) where N1 = φ∈ H1(Rn)| φ ≥ 0, φ ≡/ 0, Z Rn|∇φ| 2+ λ Z Rn φ2 = µ Z Rn φ4 (1.3b) and E(φ) = 1 2 Z Rn|∇φ| 2+λ 2 Z Rn φ2− µ 4 Z Rn φ4. (1.3c)
It is easy to see that (1.2) and (1.3) are equivalent. If φ satisfies (1.3) then φ is called a ground state solution. Hereafter, we extend the definition of ground state solutions to m-component case. We define
Nm = φ= (φ1, φ2, . . . , φm)∈ (H1(Rn))m| φj ≥ 0, φj ≡/ 0, Z Rn|∇φj| 2+ λ j Z Rn φ2j = µj Z Rn φ4j + X i6=j βij Z Rn φ2iφ2j, j = 1, . . . , m ) (1.4)
and the energy functional
E(φ) = m X j=1 1 2 Z Rn|∇φj| 2+λj 2 Z Rn φ2j − µj 4 Z Rn φ4j − 1 4 m X i, j = 1 i6= j βij Z Rn φ2iφ2j, (1.5)
where φ = (φ1, . . . , φm)∈ (H1(Rn))m. Then, consider the minimization problem
inf
φ∈Nm
E(φ). (1.6)
1. φj > 0 for all j and φ satisfy (1.1).
2. E(φ)≤ E(ψ) for any other solution ψ of (1.1). then φ is called a ground state solution of (1.1).
The main propose of this paper is first to discretize the m-coupled nonlinear Schr¨odinger equations (1.1) to an m-coupled discrete nonlinear Schr¨odinger equations (DNLS) and to develop a global convergence method for the computation of positive ground state solution of one-component case (m = 1). Secondly, we can use the positive ground state solution of one-component nonlinear Schr¨odinger equation to describe the primal stalk of solution curve for this m-component DNLS and we develop a continuation method for the computation of all possible positive bound states of a m-coupled non-linear Schr¨odinger equations (1.1). Thirdly, for the three-component DNLS with n = 2, µ1 = µ2 = µ3 = λ1 = λ2 = λ3 = 1 and β =−β12 =−β13 = β23> 0, we can prove that
the primal stalk of the solution curve of DNLS will occur at least N − p bifurcation points at finite values of coupling constants β, where N is the number of grid points and p is the number of nonnegative eigenvalues of A− I + 3[[u 2
∗ ]] (see Theorem 4.2).
This paper is organized as follows. In Section 2, we develop a iterative method for the computation of positive ground state solution of one-component DNLS. We also prove the iterative method is a global convergence method. In Section 3, we develop a stability continuation method for positive solutions of m-coupled DNLS. In Section 4, we prove the existence of the bifurcation of a 3-coupled DNLS at finite value of coupling constant β. Numerical results of positive bound states of 3-component nonlinear Schr¨odinger equation by solving the DNLS is presented in Section 5. Finally a conclusion is given in Section 6.
Throughout this paper, we use the bold face letters or symbols to denote a matrix or a vector. For u = (u1, . . . , uN)⊤, v = (v1, . . . , vN)⊤∈ RN, u◦v = (u1v1, . . . , uNvN)⊤
denotes the Hadamard product of u and v, u r = u◦· · ·◦u denotes the r-time Hadamard
product of u, [[u]] := diag(u) denotes the diagonal matrix of u. 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 and σ(A)
denotes the spectrum of A and N (A) and R(A) denote the null and range spaces of A, respectively.
2
Global convergence method for one-component
DNLS
To compute the nonlinear Schr¨odinger equations (1.1) numerically by continuation (e.g. [1, 7]), it is natural to first discretize the differential equations in (1.1) by a
finite difference method or a finite element method. Here, the domain is restricted to [−d, d] × [−d, d] where d is a sufficiently large positive real number. Suppose that the Laplacian operator △ in (1.1) is discretized by the central difference approximation with the grid size h. Then the discretization matrix A ∈ RN×N corresponding to
the operator△ with the Dirichlet boundary condition is an irreducible and symmetric negative definite matrix. ujis defined by the a approximation of φj(x) for j = 1, . . . , m.
Then the discretization of (1.1) is called m-coupled discrete nonlinear Schr¨odinger equation (DNLS), can be formulated as follows
Auj− λjuj + µju j2 ◦ uj+Pmi6=j,i=1βiju i2 ◦ uj = 0,
uj > 0 for j = 1, . . . , m,
(2.1)
where λj, µj are positive constants and βij = βji, i 6= j are coupling constants. The
energy functional E(φ) in (1.5) becomes
E(x) = m X j=1 −12u⊤j Auj+ λj 2u ⊤ j uj − µj 4u 2⊤ j u 2 j − 14 m X i6=j,i=1 βiju i2⊤u 2 j , (2.2)
where the vector x = (u⊤1, . . . , u⊤m)⊤ ∈ RN m is in
Nm = (u⊤1, . . . , u⊤m)⊤ ∈ RN m| uj ≥ 0, x ≡/ 0 and −u⊤j Auj + λju⊤j uj = µju j2⊤u 2 j + m X i6=j,i=1 βiju i2⊤u 2 j , j = 1, . . . , m ) . (2.3)
We first study the ground state of one-component nonlinear Schr¨odinger equation (1.1) described by the one-component DNLS
Au− λu + µu 2 ◦ u = 0,
u > 0, (2.4)
where λ and µ are positive real numbers. The discretization of minimization problem (1.2), can be formulated as follows
inf
u≥0E(u),b (2.5a)
where
b
E(u) = −u⊤Au + λu⊤u
The equivalent formulation is to consider the equation (1.3). We see that the associated minimization problem of (2.5) becomes
inf u∈N1E(u), (2.6a) where E(u) = −1 2 u ⊤Au + λ 2u ⊤u−µ 4u 2⊤u 2 (2.6b) and N1 =
u∈ RN| u ≥ 0, u ≡/ 0, −u⊤Au + λu⊤u = µu 2⊤
u 2
. (2.6c) Any solution u ∈ RN of (2.4) solves a local minimum or a saddle point of (2.5) or
(2.6).
In the following, we shall develop a numerical algorithm for finding the global minimum of (2.5) or (2.6), i.e., the ground state solution of (1.2) or (1.3). From (2.1), the property of A is diagonal dominant with positive off-diagonal entries, i.e., −A is an irreducible M-matrix. Let
¯
A = λI− A. (2.7)
Since λ in (2.4) is a positive constant, it is easily seen that ¯A is an irreducible M-matrix and ¯A−1 is positive definite matrix with positive entries ( ¯A−1 ≻ 0 and ¯A−1 > 0).
Define the set
M =u∈ RN| kuk
4 = 1, u≥ 0
, M = interior of M,◦ (2.8) where kuk4 = (u 2⊤u 2 )1/4. It is easy to verify that if u∈ M then
¯
A−1u = (λI− A)−1u > 0. (2.9) We now define a mapping f :M → M by
f (u) = A¯
−1u 3
k ¯A−1u 3 k
4
. (2.10)
Then f is well-defined and maps M to M by (2.9). The function f in (2.10) can be used to define the fixed point iteration:
ui+1 = f (ui).
The algorithm is given in the following: Algorithm 2.1 Fixed Point Iteration.
(i) Given ¯A∈ RN×N and u
0 > 0 with ku0k4 = 1, let i = 0;
(ii) Solve the linear system
¯
Aui+1= u i3 .
Compute ui+1 = ui+1/kui+1k4.
(iii) If converges, then u∗ ← u
i+1, stop; else i← i + 1, go to Repeat (ii).
Theorem 2.1. The function f :M → M given in (2.10) has a fixed point u∗ in M,◦
and the point
¯
u(µ) = 1 µ1/2k ¯A
−1u∗ 3
k−1/24 u∗ ∈ N1, (2.11)
which solve the DNLS (2.4).
Proof. From (2.9) and (2.10), it is easily seen that f is continuous on M and M is homeomorphic to (N−1)−dimensional standard simplex which is convex and compact. Apply Schauder fixed point theorem to f , there is a point u∗ ∈ M such that
f (u∗) = u∗. (2.12)
The fixed point u∗ ∈ M follows from the fact that the function f in (2.10) maps M◦ to M. From (2.7), (2.10) and (2.12), we have◦
k ¯A−1u∗ 3 k−14 u∗ 3 = (λI− A)u∗. (2.13) Multiplying (2.13) by µ1/21 k ¯A−1u∗ 3
k−1/24 from the left and set
¯
u(µ) = 1 µ1/2k ¯A
−1u∗ 3
k−1/24 u∗, (2.14)
then we can obtain that
A¯u(µ)− λ¯u(µ) + µ¯u(µ) 3
= 0. and it is easy to verify that ¯u(µ) belong toN1.
By Theorem 2.1, the DNLS (2.4) has a solution ¯u(µ) which can be obtained by the fixed point of f . The question is how to find the fixed point of the mapping f . The following theorems show that the iteration method proposed by Algorithm 2.1 converges globally to a fixed point u∗ ∈M of f.◦
Theorem 2.2. If u ∈ M and v = f(u) then bE(v) ≤ bE(u), where bE(·) is defined in (2.5b). The equality holds if and only if u is a fixed point of f : M → M, i.e., f (u) = u.
Proof. Since u, v∈ M, then we have kuk4 = 1 and
b
E(v) = v⊤(λI− A)v. (2.15)
Substituting v = f (u) = A¯−1u 3 k ¯A−1u 3k 4 into (2.15), we get v⊤(λI− A)v = 1 k ¯A−1u 3 k 4 v⊤u 3 .
Let c = k ¯A−11u 3k4 and apply H¨older inequality (|x⊤y| ≤ kxkpkykq 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. (2.16)
Using the fact kuk4 = 1, we have
c = u⊤A ¯¯A−1u
3
k ¯A−1u 3 k
4
= u⊤(λI− A)v. (2.17)
Since λI− A is a positive definite matrix, it has the Cholesky factorization (λI − A) = L⊤L. Applying Cauchy-Schwarz inequality to (2.17), we obtain that
c = u⊤L⊤Lv≤√u⊤L⊤Lu·√v⊤L⊤Lv
=pu⊤(λI− A)u ·pv⊤(λI− A)v. (2.18)
From (2.16) and (2.18), we have
v⊤(λI− A)v ≤ c ≤pu⊤(λI− A)u ·pv⊤(λI− A)v. (2.19)
Multiplying (v⊤(λI− A)v)−1/2 to each side of the inequality, we obtain
p
v⊤(λI− A)v ≤pu⊤(λI− A)u.
Use the fact kvk4 =kuk4 = 1, then we have
b
E(v) = v⊤(λI− A)v kvk2 4 ≤ u⊤(λI− A)u kuk2 4 = bE(u). (2.20)
The equality in (2.20) holds if and only if the equality of H¨older inequality (2.16) and the Cauchy-Schwarz inequality (2.18) hold, if and only if the vectors v and u are dependent, i.e., v = au for some a ∈ R. Since v > 0, u > 0 and kvk4 =kuk4 = 1, we
obtain v = u. Hence, the equality in (2.20) holds if and only if u is a fixed point of f .
In fact, bE(ku∗) = bE(u∗) for any positive real number k. If we can compute a
global minimizer u∗ ∈ M of (2.5), then the positive solution ¯u defined by (2.14) which
is global minimizer of (2.5), (2.6) and is the solution of DNLS (2.4). The following consequence can be easily obtained from Theorem 2.2:
Corollary 2.3. If the minimization problem (2.5) has a unique global minimizer u∗ ∈
M, then there exist a neighborhood Ru∗ of u∗ such that the fixed point iteration
con-verges to u∗ for any initial vector u0 ∈ Ru∗. In addition, ¯u(µ), defined in (2.11) is a
global minimizer of (2.5), (2.6).
Since the cost function bE(·) in (2.5b) is continuous on the compact set M, then b
E(·) : M → R+ attains a minimum value. From Theorem 2.2, it is easily seen that
the sequence { bE(ui)}∞i=1 converges to some positive number bE∗, that is
lim
i→∞E(ub i) = bE
∗. (2.21)
The sequence{ui}∞i=0⊂ M is bound, so there exists a convergent subsequence {uni}∞i=0
and u∗ ∈ M such that
lim
i→∞uni = u
∗. (2.22)
By Theorem 2.2, we have f (u∗) = u∗. This implies bE(f (u∗)) = bE(u∗).
Since ¯A is positive definite, we define the ¯A-norm of u by kukA¯ =
√
u⊤Au. Then¯
kui+1− uik2A¯ = (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, (2.23)
and from (2.17) and (2.19) we have
b
E(ui+1)≤ u⊤i+1Au¯ i ≤
q b E(ui+1) q b E(ui). (2.24)
Thus by (2.23) and (2.24) we see that
kui+1− uikA¯ ≤
q b
that is,
lim
i→∞kui+1− uikA¯ = 0. (2.26)
Theorem 2.4. If u∗ given in (2.22) is strictly local minimum of (2.5) or (2.6), then
the sequence {ui}∞i=0 generated by Algorithm 2.1 converges to u∗ ∈ M.
Proof. Since u∗ is strictly local minimum of the optimization problems (2.5) or (2.6), then the Hessian matrix H(u∗) of bE(u) is positive definite. Then there exists δ
1 > 0
such that H(u) is positive definite, i.e., bE(u) is convex, for u∈ M and ku−u∗k¯ A < δ1. Let b E1 = min ku−u∗k¯ A=δ1 b E(u) > bE∗,
then there exists δ2 > 0 such that
k bE(u)− bE∗k < Eb1− bE∗
2 for ku − u
∗k¯
A < δ2. (2.27)
From (2.22) and (2.26), for any positive number ε with ε < min{δ1/2, δ2}, there
exists N0 ∈ N such that
kunj − u
∗k ¯
A < ε and kui+1− uikA¯ < ε for nj, i > N0. (2.28)
Since ε < min{δ1/2, δ2}, using the fact bE(ui) is decreasing in i, it is easily seen that if
kui − u∗kA¯ < ε and kui+1− uikA¯ < ε thenkui+1− u∗kA¯ < ε. Thus
kui− u∗kA¯ < ε, for all i > N0.
This complete the proof.
From above theorem, we know that if one of limit points of {ui}∞i=0, u∗, which
is strictly local minimum of (2.5) or (2.6), then {ui}∞i=0 generated by Algorithm 2.1
converges globally to a fixed point u∗ ∈M and u◦ ∗ satisfies
(λI− A)u∗ = τ u∗ 3
where τ =k ¯A−1u∗ 3
k4. (2.29)
From Theorem 2.1, we can compute ¯u(µ) which solve the DNLS (2.4).
Remark 2.1. Numerical experience shows that for any arbitrary initial positive vector u0 withku0k4 = 1, the fixed point iteration converges to the global minimizer of (2.5).
3
Stable continuation method for DNLS
In this section, we develop a stable continuation method to solve the m-coupled DNLS (2.1). To study the phase separation of m-coupled DNLS, we assume that
βij = δijβ, i6= j, i, j = 1, . . . , m, (3.1)
where δij = δji is a constant for i6= j and β ≥ 0 is a parameter. As δij = δji > 0, the
interaction between the i-th and j-th component of the beam is attractive, otherwise the interaction is repulsive.
Let x = (u⊤
1, . . . , u⊤m)⊤ ∈ RN m. The m-coupled DNLS can be rewritten by
G(x, β) = 0, (3.2a) where G = (G1, . . . , Gm) : RN m× R → RN m is a smooth mapping with
Gj(x, β) = Auj− λjuj + µju j2 ◦ uj+ β m
X
i6=j,i=1
δiju i2 ◦ uj, j = 1, . . . , m. (3.2b)
We denote the Jacobian of G by
DG = [Gx, Gβ]∈ RM×(M+1)
with M = N m, and the solution curve C of (2.1) by
C =y(s) = (x(s)⊤, β(s))⊤| G(y(s)) = 0, s ∈ R . (3.3) Here we assume a parametrization via arc lengths is available on C. By differentiating the equation (3.2) with respect to s we obtain
DG(y(s)) ˙y(s) = 0, (3.4) where ˙y(s) = ( ˙x(s)⊤, ˙β(s))⊤ denotes a tangent vector toC at y(s), that is, the tangent vector to C at y(s) is the nontrivial solution of homogenous system DG(y(s))w = 0, where DG(y(s)) is M × (M + 1) matrix.
Remark 3.1. We consider 2-coupled nonlinear Schr¨odinger equation (1.1) with n = 1, that is, the domain of φj is R. If we differentiate both side of (1.1a) then we have
L1 2β12φ1φ2 2β12φ1φ2 L2 φ′1 φ′2 = 0, (3.5) where L1 = d 2 dx2− λ1+ 3µ1φ21+ β12φ22 and L2 = d 2 dx2− λ2+ 3µ2φ22+ β12φ21. From (3.5) we
see that the matrix
L1 2β12φ1φ2
2β12φ1φ2 L2
is singular. It easily seen that the solution set of (1.1) is one dimensional. Actually, the solution set of (1.1) is a one dimensional solution curve, {φr(x)|φr(x) = (φ1(x− r), φ2(x− r)), r ∈ R}.
Now, we consider the m-coupled DNLS (3.1) with domain [−d, d] × [−d, d]. If the size of the square domain d is sufficiently large and the grid size h is sufficiently small. From Remark 3.1, the numerical rank of Gx(x(s), β(s)) should be less than M and the
numerical null space is spanned by
K0 = span{Dxx(s), Dyx(s)}, (3.6)
where
Dx= diag{Dx, . . . , Dx}, Dy = diag{Dy, . . . , Dy} ∈ RN m×Nm, (3.7)
Dx, Dy ∈ RN×N are discretization matrices of the differential operators ∂x∂ , ∂y∂ ,
respec-tively, with discretized by the central difference approximation of the grid size h. If we use vectors [(Dxx(s))⊤, 0]⊤, [(Dyx(s))⊤, 0]⊤ as our prediction directions then,
nu-merically, there are solutions which are translation by x-axis and y-axis from x(s), respectively. Those solutions have very small residual (kG(x, β)k < ε), so those solu-tions are also called “ε-solusolu-tions”. Since the square domain is bounded, the residual is getting large if we keep shifting x(s) along the x-axis or y-axis, where (x(s), β(s)) is solution of G(y(s)) = 0. In the prediction-correction scheme (of continuation method), there are three serious problems if the ε-solutions exist.
1. In prediction, we can not compute good prediction direction by solving (3.4). 2. In correction, the Newton’s method is chosen as a corrector. Since the Jacobian
matrix Gxis near singular, the Newton’s method will lose quadratic convergence.
3. Since there is a 2-dimensional ε-solution set containing the solution curveC (3.3) and we can not get the “good” prediction direction, then the computed solution path by continuation method may appear randomly in the ε-solution set. That is, we can not trace the solution curve C (3.3) by continuation method if the ε-solutions exist.
Let yi = (x⊤i , βi)⊤ ∈ RM +1 be a point that has been accepted as an approximating
point for the solution curve C. A “good” prediction direction ˙yi = ( ˙x⊤i , ˙βi)⊤ should be
satisfy (3.5) and the vector ˙x⊤
i should be in K0⊥, where K0 is given in (3.6). So the
good prediction direction ˙yi ∈ RM +1 satisfies the linear bordered system
Gx Gβ a⊤x 0 a⊤ y 0 c⊤ i ci ˙yi = 0 0 0 1 , (3.8)
where ax = Dxxi and ay = Dyxi and some suitable constant vector (c⊤i , ci)⊤ ∈ RM +1.
Suppose that the Euler predictor,
yi+1,1 = yi+ hi˙yi
is used to predict a new point yi+1,1, where hi > 0 is the step length and ˙yi is the
unit tangent vector at yi which is obtained by solving the linear bordered system (3.8).
The accuracy of the approximation yi+1,1 to the solution curve C can be improved
by a corection process. Typically, Newton’s method is chosen as a corector, i.e., the following linear bordered system
Gx(yi+1,l) Gβ(yi+1,l) a⊤ x 0 a⊤y 0 ˙x⊤ i β˙i δl= −G(yi+1,l) ρx,l ρy,l ρl , l = 1, 2, . . . , (3.9)
with ρl = ˙yi⊤(yi+1,l− yi+1,1), ρx,l = a⊤x(xi+1,l− xi+1,1) and ρy,l = a⊤y(xi+1,l− xi+1,1), is
solved by setting yi+1,l+1 = yi+1,l+ δl, l = 1, 2, . . . . If {yi+1,l} converges until l = l∞,
then we accept yi+1 = yi+1,l∞ as a new approximation to the solution curveC.
Those two linear systems (3.8) and (3.9) are overdetermined systems. It is easily seen that if we consider the least square problem on those two linear systems (3.8) and (3.9) then the minimum residual is very small. In fact, linear systems (3.8) and (3.9) can be rewritten in the form
B f g⊤ γ x σ = q ρ , x(M + 1, 1) = x(M + 2, 1) = 0, (3.10) where B = Gx ax ay a⊤ x 0 0 a⊤ y 0 0 ∈ R(M +2)×(M+2) (3.11)
is symmetric and f , g, q ∈ R(M +2). The linear system (3.10) can be easily solved by the
well-known block elimination (BE) algorithm (see e.g., [7]) when B is well-conditioned. However, near turning points or branch points, B in (3.10) becomes nearly singular, i.e., B is ill-conditioned. Then the linear system should be solved by the deflated block elimination (DBE) algorithm by Chan [3], or the more efficient, backward stable, mixed block elimination (BEM) algorithm proposed by Govaerts [6, 5].
(i) Solve ξ⊤B = g⊤, (ii) Compute δ1 = γ− ξ⊤f , σ = (ρ− ξ⊤q)/δ1, (iii) Solve Bv = f , (iv) Compute δ = γ− g⊤v, q 1 = q− fσ, ρ1 = ρ− γσ, (v) Solve Bw = q1, (vi) Compute σ1 = (ρ1− ξ⊤w)/δ, (vii) Compute x = w− vσ1, σ = σ + σ1.
From Algorithm 3.1, we see that the main step in (3.8) or in (3.9) is to solve a linear system of the form Bξ = g.
Remark 3.2. Since those linear systems (3.8) and (3.9) have very small minimum resid-ual, then the equation (3.10) is near consistent. Thus, the solution x solved by Algo-rithm 3.1 satisfies x(M + 1, 1)≈ x(M + 2, 1) ≈ 0 automatically.
Testing for Bifurcation.
LetC be the curve defined in (3.3), y(s) ∈ C and J(s) = Gx(y(s)) Gβ(y(s)) a⊤ x 0 a⊤ y 0 ∈ R(M +2)×(M+1), (3.12a) J(s) = Gx(y(s)) a⊤x a⊤y ∈ R(M +2)×M. (3.12b)
As was described in [1, 4, 7] a point y(s) ∈ C is said to be a regular point if rank(J(s)) = M (dimN (J(s)) = 1), and is a singular point if rank(J(s)) ≤ M −1 (dim N (J(s)) ≥ 2). For a regular point y(s), the tangent vector ˙y(s) is uniquely determined by the linear system (3.8).
Now our task is to design algorithms to detect singular points of the solution curve C and to compute tangent vectors if y(s) is singular point. For simplicity, here we only consider the case
[Gβ(y(s))⊤, 0, 0]⊤ ∈ R(J(s)) for each singular points y(s) ∈ C, (3.13)
i.e., dimN (J(s)) = dim N (J(s)) − 1 and dim N (J(s)) ≥ 2.
Remark 3.3. We only consider the condition (3.13) holds at singular points y(s) ∈ C not all of points inC. This condition say that the tangent vectors at singular point have nonzero component at ˙β(s) and the condition (3.13) con be expected in a situation of the DNLS (3.2).
For convenience, hereafter we assume that B(s) is the matrix B given by (3.10) at the point y(s)∈ C. The following theorem can be easily obtained.
Theorem 3.1. If the condition (3.13) holds, then the following statements are equiv-alent
(1) rank(J(s))≤ M − 1, (2) N (J(s)) 6= {0},
(3) B(s) is singular and there is a nonzero vector z ∈ RM such that [z⊤, 0, 0]⊤ ∈
N (B(s)).
We propose the following algorithm to find the tangent vectors at singularity.
Algorithm 3.2 Tangent Vectors at Singularity. (I) For dimN (J(s∗)) = 1 [7, p. 88-99]:
(i) Compute the unit right null vector φ = ( ¯φ⊤, 0, 0)⊤ of B(s∗), and solve
J(s∗) ¯φ
0 = −[Gβ(s∗)⊤, 0, 0]⊤ with ¯φ⊤φ¯0 = 0, by using sparse SVDPACK
[2]; (ii) Form φ1 = ¯φ 0 and φ2 = ¯φ0 1 ;
(iii) Solve the real vector roots{(ˆµk, ˆνk)}2k=1 of a11µ2+ 2a12µν + a22ν2 with
a11 = ¯φ⊤Gxx(s∗) ¯φ ¯φ, a12 = ¯φ⊤[Gxx(s∗) ¯φ0+ Gxβ] ¯φ,
a22= ¯φ⊤[Gxx(s∗) ¯φ0φ¯0+ 2Gxβ(s∗) ¯φ0+ Gββ(s∗)];
(iv) Form tangent vectors ˙yk(s∗) = ˆµkφ1+ ˆνkφ2, k = 1, 2.
(II) For dimN (J(s∗)) = ℓ≥ 2:
(i) Compute the unit right null vectors φ(1), . . . , φ(ℓ) of B(s∗) with φ(k) =
( ¯φ(k)⊤, 0, 0)⊤, k = 1, . . . , ℓ and solve J(s∗) ¯φ0 = −[Gβ(s∗)⊤, 0, 0]⊤ with
¯
φ(k)⊤φ¯0 = 0, k = 1, . . . , ℓ, by using sparse SVDPACK [2]; (ii) Form φk= ¯ φ(k) 0 , k = 1, . . . , ℓ and φℓ+1 = ¯φ0 1 ;
From Theorem 3.1, to detect singular points of the solution curveC is equivalent to find the null space of symmetric matrix B defined in (3.10). We use the inverse power method to to compute the smallest eigenvalues in modulus of B.
Algorithm 3.2 Inverse Power Method.
(i) Given a unit vector ζ0 ∈ R(M +2), and let l = 1,
(ii) Repeat l: until convergence,
Solve Bbζl= ζl−1, where B is given in (3.11). Set ζl = bζl/kbζlk2, µ(l) = ζ⊤l Bζl;
(iii) If converges, then µ(s)← µ(l); else l← l + 1, Goto Repeat (ii).
Let µ(s1) and µ(s2) be the smallest eigenvalues in modulus of B(s1) and B(s2),
re-spectively, where s1 < s2 are two consecutive parameters. If µ(s1) > 0 and µ(s2) < 0,
then there is a s∗ ∈ (s1, s2) such that B(s∗) is singular. We propose the following
algorithm to detect the singular point of C.
Algorithm 3.3 Detection of Singularity of C.
(i) Given µ(si) the smallest eigenvalue in modulus of B(si), i = 1, 2, where µ(s1) > 0,
µ(s2) < 0, e.g., |µ(s1)| ≈ |µ(s2)| ≈ 10−4.
(ii) Do Secant Method: until convergence,
(a) Compute y1(s0) := y(s0) = y(s1) +
tµ(s1)
µ(s2)− µ(s1)
, where t = y(s1)− y(s2),
(b) Perform Newton Correction (3.9): until convergence (i.e., ℓ = ℓ∞),
Solve Gx(yℓ(s0)) Gβ(yℓ(s0)) a⊤x 0 a⊤ y 0 t⊤ δℓ = −G(yℓ(s0)) ρx,ℓ ρy,ℓ ρℓ with ρℓ = t⊤(yℓ(s0)− y1(s0)), set yℓ+1(s0) = yℓ(s0) + δℓ, ℓ← ℓ + 1, Goto (b).
(c) Compute µ(s0) of B(s0) with y(s0) = yℓ∞(s0) using Algorithm 3.2,
(e) If µ(s0) > 0, s1 ← s0, else s2 ← s0, Goto (ii);
(iii) Call Algorithm 3.1 to compute the desired tangent vectors with y(s∗) = y
ℓ∞(s0).
4
3-coupled nonlinear Schr¨
odinger equations
In this section, we focus on the special case
m = 3, λ1 = λ2 = λ3 = µ1 = µ2 = µ3 = 1. (4.1)
Then the 3-coupled DNLS G(x, β, δ) = 0 in (3.1), where δ = (δ12, δ13, δ23) and β > 0,
can be rewritten by
Au1− u1+ u 13 + βδ12u2 2u1+ βδ13u 32u1 = 0, (4.2a)
Au2− u2+ u 23 + βδ12u1 2u2+ βδ23u 32u2 = 0, (4.2b)
Au3− u3+ u 33 + βδ13u1 2u3+ βδ23u 22u3 = 0. (4.2c)
Define the matrix Σ = 1 |βδ12| |βδ13| |βδ12| 1 |βδ23| |βδ13| |βδ23| 1 . It is shown in [8] that 1. If δ12< 0, δ13< 0 and δ23< 0, then the ground state solution doesn’t exist.
2. If δ12 > 0, δ13 > 0, δ23 > 0 and Σ is positive definite then the ground state
solution exists.
3. If δ12 < 0, δ13 < 0, δ23 > 0 and Σ is positive definite then the ground state
solution doesn’t exist.
4. If δ12 > 0, δ13 > 0, δ23 < 0, β ≪ 1 and the ground state solution exists then it
must be non-radially symmetric.
To simplify our computations, we shall focus on the following two cases
δ12 = δ13 =−1, δ23 = 1, (4.3a)
δ12 = δ13 = 1, δ23 =−1. (4.3b)
For the case (4.3a), the DNLS G(x, β) = 0 in (4.2) can be rewritten by Au1− u1+ u 13 − βu 2 2 u1− βu 32 u1 = 0, (4.4a) Au2− u2+ u 23 − βu 2 1 u2+ βu 32 u2 = 0, (4.4b) Au3− u3+ u 33 − βu 2 1 u3+ βu 22 u3 = 0, (4.4c)
where β > 0. It is easily seen that if we set β := −β then (4.4) will become to case (4.3b). So, to investigate these two cases, we only need to consider the DNLS (4.4) for β ∈ R. By letting
u2 = u3 = κu1. (4.5)
It follows equation (4.4b) and (4.4c) are identical. Thus, the equation (4.4) can be reduced to
Au1− u1+ (1− 2βκ2)u 13 = 0,
Au1− u1+ (κ2− β + βκ2)u 13 = 0.
(4.6)
From Theorem 2.1, (4.6) has positive solution u1 if 1− 2βκ2 = κ2− β + βκ2, then
κ =± s
1 + β 1 + 3β. Since u2 = u3 = κu1 > 0, this implies
κ = s
1 + β
1 + 3β. (4.7)
Then the positive solution u1 of equation (4.6) satisfy
Au1− u1+
1 + β− 2β2
1 + 3β u
3
1 = 0. (4.8)
From Theorem 2.1 the solution of (4.8) is
u1 =
s
1 + 3β
1 + β− 2β2u∗ = ηu∗, (4.9)
where u∗ is a solution of Au− u + u 3 = 0. From Theorem 2.1,
u∗ = ¯u(1) =k ¯A−1u∗ 3
k−1/24 u∗, (4.10)
where u∗ is the fixed point of f :M → M with λ = 1.
• If β → 1− then κ→q1
2 (by (4.7)) and from (4.5) and (4.9) we have
• If β → −1 3 −
then u1 → 0 (by (4.9)) and from (4.4b), (4.4c) and Theorem 2.2, we
have
u2 = u3 →
r 3 2u∗.
Theorem 4.1. The primal stalk of 3-coupled DNLS (4.4) can be described by u1 = q 1+3β 1+β−2β2u∗, u2 = u3 = q 1+β 1+β−2β2u∗, −1 3 ≤ β < 1 (4.11) where u∗ is given in (4.10).
The solution curve C is defined as in (3.3) coresponding to the DNLS (4.4) can be rewritten by
C = {y(s) = (x⊤(s), β(s))⊤|G(y(s)) = 0 is given in (4.4) and s ∈ R} (4.12) where x(s) = (u1(s), u2(s), u3(s))⊤.
Theorem 4.2. The primal stalk of solution curveC in (4.12) given by (4.11), undergoes at least N − p bifurcation points at finite values 0 < β = β∗
q < 1, q = 1, . . . , N − p,
where p is the number of nonnegative eigenvalues of A− I + 3[[u 2
∗ ]].
Proof. Since (4.4) has positive solution curve u2(β) = u3(β) = κu1(β) for 0 < β < 1
where κ is defined in (4.7), the Jacobian matrix of (4.4) with respect to u is of the form Gu(y(β)) = B1 E1 E1 E1 B2 E2 E1 E2 B2 (4.13) where B1 = A− I + [[3u 12 − 2βu 2 2 ]], (4.14a) B2 = A− I + [[(3 + β)u 22 − βu 2 1 ]], (4.14b) E1 =−2β[[u1◦ u2]], (4.14c) E2 = 2β[[u 22]]. (4.14d)
From (4.5), (4.7) and (4.9), we have
u2 = u3 = s 1 + β 1 + 3βu1 = s 1 + β 1 + β − 2β2u∗ (4.15)
where u∗ is given in (4.10). Substituting (4.9) and (4.15) into (4.14), we get B1 = A− I + 3 + 7β− 2β2 1 + β− 2β2 [[u 2 ∗ ]], (4.16a) B2 = A− I + 3 + 3β− 2β2 1 + β− 2β2 [[u 2 ∗ ]], (4.16b) E1 =−2β p (1 + β)(1 + 3β) 1 + β− 2β2 [[u 2 ∗ ]], (4.16c) E2 = 2β + 2β2 1 + β− 2β2[[u 2 ∗ ]]. (4.16d) Let Q = I 0 0 0 I I 0 0 I , then QGu(y(β))Q−1 = B1 E1 0 2E1 B2+ E2 0 E1 E2 B2− E2 . (4.17) Hence σ(Gu(y(β))) = σ B1 E1 2E1 B2+ E2 ∪ σ(B2− E2). (4.18)
On the other hand, B1 E1 2E1 B2+ E2 = B1 0 0 B1 − 2β 0 κη2[[u 2 ∗ ]] 2κη2[[u 2 ∗ ]] 2β+11 [[u 2 ∗ ]] = 1 0 0 1 ⊗ B1− 0 1 2 (2β+1)κη1 2 ⊗ 2βκη2[[u 2 ∗ ]] = 1 0 0 1 ⊗ B1− 0 1 2 a ⊗ 2βκη2[[u 2 ∗ ]] (4.19)
where κ and η are given by (4.7) and (4.9), respectively, and
a = 1 (2β + 1)κη2 = 1− β p (1 + 3β)(1 + β). (4.20) Since a+√a2+8 2 and a− √ a2+8 2 are eigenvalues of 0 1 2 a
, then from (4.16), (4.18) and (4.19), it follows
where Λ1(β) = σ A− I + 4β + 3 2β + 1[[u 2 ∗ ]] , (4.22a) Λ2(β) = σ A− I + 3 + 7β− 2β2 1 + β− 2β2 − (a + √ a2+ 8)βκη2 [[u 2 ∗ ]] , (4.22b) Λ3(β) = σ A− I + 3 + 7β− 2β2 1 + β− 2β2 − (a − √ a2+ 8)βκη2 [[u 2 ∗ ]] , (4.22c) Since (a +√a2+ 8)βκη2 = (a + √ a2+ 8)βp(1 + β)(1 + 3β) 1 + β − 2β2 = β(1− β) + β p (1− β)2+ 8(1 + 3β)(1 + β) 1 + β − 2β2 = 4β 2+ 4β 1 + β − 2β2 =−3 + 3 + 7β− 2β 2 1 + β− 2β2 . then Λ2(β) = σ(A− I + 3[[u ∗2]]). (4.23) Hence σ(Gu(y(β)))|β=0= Λ1(0)∪ Λ2(0)∪ Λ3(0) = σ(A− I + 3[[u 2
∗ ]])∪ σ(A − I + 3[[u ∗2]])∪ σ(A − I + 3[[u ∗2]]).
(4.24) When β → 1−, then Λ1(β)→ σ(A − I + 7 3[[u 2 ∗ ]]), (4.25) and 3+7β1+β−2β−2β22 − (a − √
a2+ 8)βκη2 → ∞, so there exists a β∗, 0 < β∗ < 1, such that
Λ3(β)⊂ R+, for β > β∗. (4.26)
If the number of nonnegative eigenvalues of Λ3(0) = σ(A− I + 3[[u ∗2 ]]) is p then
from (4.23), (4.24), (4.25) and (4.26), we see that the primal stalk of solution curve C in (4.12), undergoes at least N − p bifurcation points at finite values 0 < β∗
q < 1,
Remark 4.1. It is proved in [8, Lemma 1] that the number of nonnegative eigenvalues of △φ − φ + 3ω2 ∗φ = λφ, φ∈ H2(Rn),
is n + 1, where ω∗ is the unique solution of △φ − φ + φ3 = 0, φ > 0 in Rn, ω(x)→ 0 as |x| → ∞.
In the square domain (n = 2), it seems that the number of nonnegative eigenvalues of Λ1(0) = σ(A − I + 3[[u ∗2 ]]), where u∗ is given in (4.10), is 3. In numerical test,
we compute the eigenvalues of Λ1(0), and the number of nonnegative eigenvalues (the
eigenvalue bigger than −10−3) is 3.
5
Numerical examples
We use the fixed point iteration algorithm which is developed in Section 2 to compute the ground state solution of one-component nonlinear schr¨odinger equation, i.e., the positive solution of DNLS (2.4). In Section 3, we developed a stability continuation method which can compute positive bound states of m-coupled nonlinear schr¨odinger equations. In the numerical simulation, we consider the DNLS (2.1) with m = 3 and λ1 = λ2 = λ3 = µ1 = µ2 = µ3 = 1.
Simulation 1.
We first compute the solution curve of positive solutions of DNLS (4.2) with δij is
given in (4.3a) on a square domain [−5, 5] × [−5, 5] with grid size h = 1/3, i.e., the solution curve of
C+ =(x⊤, β)⊤| G(x, β) = 0 is given in (4.4) for β ∈ R +
. (5.1) In Figure 1 and 2, we plot the bifurcation diagram of the solution curve and energy curves versus β ∈ [0, 1.2]. The nodal domains of positive bound states are attached near the solution curve, where the left, middle and right figures are the level set of u1,
u2 and u3, respectively.
Simulation 2.
We compute the solution curve of positive solutions of DNLS (4.2) with δij is given
0.237 0.253 0.385 0.387 1 R 3 N β Figure 1. Bifurcation curves of DNLS.
E
(x
)
β Figure 2. Energy curves of DNLS.
curve of
C−=(x⊤, β)⊤| G(x, β) = 0 is given in (4.2) with (5.2) δ12 = δ13 = 1 and δ23=−1, for β ∈ R+} .
In Figure 3 and 4, we plot the bifurcation diagram of the solution curve and energy curves versus β ∈ [0, 0.5]. The nodal domains of positive bound states are attached near the solution curve, where the left, middle and right figures are the level set of u1,
u2 and u3, respectively.
In [8], it is shown that in the case of Simulation 2, when β sufficiently small, there exit a non-radially positive solution and its energy is lower than the energy of the positive solution given in (4.11). But in Figure 3, when β sufficiently small, we can only compute radially solution. For computing the non-radially positive solution described in [8], we can follow the following steps to compute the non-radially positive solution.
0.373 0.333 0.319 R 3 N β Figure 3. Bifurcation curves of DNLS.
E
(x
)
β Figure 4. Energy curves of DNLS.
Step 1. Let C1 =
(x⊤, β)⊤| G(x, β) = 0 is given in (4.2) with (5.3) δ12= δ13= δ23 = 1, for 0≤ β ≤ 0.2} .
We trace the solution curve with 0 ≤ β ≤ 0.2 by continuation method. Step 2. Fix β = 0.2, then trace the solution curve
C2 =
(x⊤, δ23)⊤| G(x, δ23) = 0 is given in (4.2) with (5.4)
β = 0.2, δ12 = δ13= 1, for − 1 ≤ δ23≤ 1}
from δ23= 1 to δ23 =−1 by continuation method.
Step 3. Fix δ23=−1, then trace the solution curve
C3 =
(x⊤, β)⊤| G(x, β) = 0 is given in (4.2) with (5.5) δ12 = δ13 = 1 and δ23=−1, for β ∈ R}
from β = 0.2 by continuation method. Simulation 3.
We compute the solution curve of positive solution of DNLS (4.2) following the previous steps on a square domain [−5, 5] × [−5, 5] with grid size h = 1/3. In Figure
5 and 6, we plot the bifurcation diagram of the solution curve C2 and energy curves
versus δ23 ∈ [−1, 1]. The nodal domains of positive bound states are attached near the
solution curve, where the left, middle and right figures are the level set of u1, u2 and
u3, respectively. In Figure 5, there is only one bifurcation point in the solution curve
C2 and if we trace the primal stalk ofC2 then the terminal point attain the primal stalk
of C−. In Figure 6, it show that solution curve of bifurcation branch has lower energy.
d
R
3
N
Figure 5. Energy curves of DNLS.
d
E
(x
)
Figure 6. Energy curves of DNLS.
In Figure 7 and 8, we plot the bifurcation diagram of the solution curve C3 and
energy curves versus β ∈ [0, 1] and the nodal domains of positive bound states are attached near the solution curve. In Figure 7, it show that the non-radially positive solution which is described in [8] for β sufficiently small. In Figure 6, it show that the energy of solution curve C3 is lower than solution curve C−.
References
[1] E.L. Allgother and K. Gerog. Numerical path following. Handbook of Nemerical Analysis, Edited by P. G. , 5, 1996.
0.985
β
R
3
N
Figure 7. Energy curves of DNLS.
β
E
(x
)
Figure 8. Energy curves of DNLS.
[2] M. Berry. Large scale singular value computations. Internat. J. Supercomputer Appl., 6(1):13–49, 1992.
[3] T. F. Chan. Deflation techniques and block-elimination algorithm for solving bor-dered sungular systems. SIAM J. Sci. Stat. Compu., 5:121–134, 1984.
[4] Willy J. F. Govaerts. Numerical Methods for Bifurcations of Dynamical Equilibria Computations. SIAM, Philadelphia, 2000.
[5] W. Govaerts and J.D. Pryce. Mixed block elimination for linear systems with wider borders. IMA J. Numer. Anal., 13:161–180, 1993.
[6] W. Govaerts. Stable solvers and block elimination for bordered systems. SIAM J. Matrix Anal. Appl., 12:469–483, 1991.
[7] H. B. Keller. Lectures on Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin, 1987.
[8] Tai-Chia Lin and Juncheng Wei. Ground state of n coupled nonlinear Schr¨odinger equations in Rn, n≤ 3. Commun. Math. Phys., 255:629–653, 2005.
β
E
(x
)
