A Hyperplane-Constrained Continuation
Method for Near Singularity in Coupled
Nonlinear Schr¨
odinger Equations
Yuen-Cheng Kuo
∗Wen-Wei Lin
†Shih-Feng Shieh
‡Weichung Wang
§May 31, 2010
Abstract
The continuation method is a useful numerical tool for solving dif-ferential equations to obtain multiform solutions and allow bifurca-tion analysis. However, when a standard continuabifurca-tion method is used to solve a type of time-independent m-coupled nonlinear Schr¨odinger (NLS) equations that can be used to model nonlinear optics, nearly singular systems arise in the computations of prediction and correc-tion search direccorrec-tions and deteccorrec-tions of bifurcacorrec-tions. To overcome the stability and efficiency problems that exist in standard continua-tion methods, we propose a new hyperplane-constrained continuacontinua-tion method by adding additional constraints to prevent the singularities while tracking the solution curves. Aimed at the 3-coupled DNLS equations, we conduct theoretical analysis to the solutions and bifur-cations on the primal stalk solution curve. The proposed algorithms have been implemented successfully to demonstrate numerical solu-tion profiles, energies, and bifurcasolu-tion diagrams in various settings. Keywords. Hyperplane-constrained continuation method, coupled nonlinear Schr¨odinger equations, numerical solutions, primal stalk solutions, bifurcation analysis.
∗Department of Applied Mathematics, National University of Kaohsiung, Kaohsiung
811, Taiwan ([email protected]).
†Department of Mathematics, National Tsing Hua University, Hsinchu 300, Taiwan
‡Department of Mathematics, National Taiwan Normal University, Taipei 116,
Tai-wan ([email protected]).
§Department of Mathematics, National Taiwan University, Taipei 106, Taiwan
1
Introduction
We intend to study bound state solutions of the following time-independent m-coupled nonlinear Schr¨odinger (NLS) equations numerically,
4φj− λjφj+ µj|φj|2φj+ m X i6=j,i=1 βij|φi|2φj = 0, in Rn, (1a) φj > 0 in Rn, j = 1, . . . , m, (1b) φj(z) → 0, as |z| → ∞, (1c)
where λj, µj > 0 are positive constants, n ≤ 3, and βij = βji (i 6= j) are
coupling coefficients. The NLS equations may be used to model a physical phenomenon in nonlinear optics [1]. In such a case, the solution φj denotes
the j-th component of the beam in Kerr-like photorefractive media. The positive constant µj is for self-focusing in the j-th component of the beam.
λj is referred to as the chemical potential. The coupling constant βij is the
interaction between the i-th and j-th component of the beam. For βij> 0,
the interaction is attractive; otherwise, the interaction is repulsive. To solve the NLS equations (1) numerically, we consider the correspond-ing 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, (2)
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 Laplacian operator with the homogeneous Dirichlet bound-ary conditions. Additionally, it can be seen that A is an irreducible and symmetric negative definite matrix. The size of N depends on the com-putational domain and grid sizes. For example, if a uniform grid size h is applied on a square finite domain [−d, d] × [−d, d] for n = 2, we have N = (2dh − 1)2. For u = (u
1, . . . , uN)> and 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.
While continuation methods may act as an effective technique for find-ing multiple solutions of the NLS equations, one particular characteristic of Equation (1) prevents standard continuation methods from being practical. Since the solution domain is unbounded in (1), a shift in any solution of the system remains a solution of the system. Consequently, a small shift in a numerical solution of an m-coupled DNLS equations (on a bounded computational domain) can also be an approximate solution to another m-coupled DNLS equations with a small residual. This phenomenon prevents classical continuation methods from being applicable or efficient for solving
the target problem. First, improper prediction directions may be obtained as the prediction directions can not be uniquely determined by numerics. Second, the Jacobian matrix may be nearly singular and the corresponding Newton’s correction process becomes inaccurate and inefficient. Third, the numerical singularity also makes detections of bifurcation points difficult. Finally, these improper search directions may result in undesired solution curves in continuation methods.
Aiming at the m-coupled DNLS equations, we make the following con-tributions to both numerical scheme development and numerical analysis. • To circumvent the obstacles mentioned above, we propose a novel hyperplane-constrained continuation method to compute all possi-ble positive bound states of m-coupled DNLS equations by adding additional hyperplane constraints. By doing so, the prediction direc-tions can be uniquely determined and the correction direcdirec-tions can be computed efficiently.
• Aiming at the 3-coupled DNLS equations, we characterize the primal stalk solutions and show that the primal stalk solution curve has at least N − p bifurcation points at finite values of coupling constants β. Here, N is the number of grid points and p is the number of nonnegative eigenvalues of a specified matrix.
• By using the proposed hyperplane-constrained continuation method, we conduct numerical experiments to explore the versatility of numer-ical solutions by presenting solution profiles, bifurcation diagrams, and the corresponding energies for various settings.
The target problem has been considered in several cases. For n = 1, i.e. the spatial dimension is one, the system (1) is integrable. Many analytical and numerical results on solitary wave solutions of m-coupled NLS equations have been well-studied [9, 13, 14, 15]. For n = 2 and m = 1, physical experiments in [18] are conducted to observe 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) solutions and self-trapped beams. A general theorem for the existence of higher dimensional m-component solutions was first proven in [17]. The sign of coupling constants βij’s is crucial for the existence of ground state
solutions. For m = 3, when all βij’s are positive, there exists a ground state
solution that is radially symmetric. Furthermore, a positive bound state solution that is non-radially symmetric is also found. See [17] for details.
It is worth mentioning that if the NLS equations are equipped with trap potentials, the difficulties resulting from solution shifting no longer exist and following numerical methods can be used to solve the equations. Bao
proposed a normalized gradient flow method [3, 4] and a time-splitting sine-spectral method [3]. For the time-independent case, a Gauss-Seidel-type iteration has been proposed in [8]. Furthermore, a continuation BSOR-Lanczos-Galerkin method has been developed in [7, ?]. More recently, the technique of Liapunov-Schmidt reduction and continuation method has been developed in [?].
This paper is organized as follows. In Section 2, we develop a hyperplane-constrained continuation method to compute positive bound state solutions of m-coupled DNLS equations. In Section 3, we prove the existence of the bifurcation of a 3-coupled DNLS equations at finite values of the coupling constant β. Numerical results of positive bound states of some 3-coupled DNLS equations are presented in Section 4. We conclude the paper in Section 5.
Throughout this paper, we use boldfaced letters or symbols to denote a matrix or a vector. For u = (u1, . . . , uN)>, [[u]] := diag(u) denotes the
di-agonal 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) en-tries, 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. Finally, we define 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 i 6= j βij Z Rn φ2iφ 2 j, (3)
where φ = (φ1, . . . , φm) ∈ (H1(Rn))m. The corresponding energy
func-tional for the m-coupled DNLS equations (2) is defined as
E(x) = m X j=1 −1 2u > j Auj+ λj 2 u > juj− µj 4 u 2> j u 2 j −1 4 m X i, j = 1 i 6= j βiju 2> i u 2 j , (4)
where the vector x = (u>1, . . . , u>m)>∈ RN m
2
The hyperplane-constrained continuation
method
In this section, we develop a novel hyperplane-constrained continuation method to solve the m-coupled DNLS equations (2). We assume that
the variable βij is changeable and then introduce the continuation method
parameter β ≥ 0 into Equation (2) by rewriting
βij= ζijβ (5)
for i, j = 1, . . . , m and i 6= j. Here ζij’s are nonzero constants and ζij = ζji.
If ζij > 0, the interaction between the i-th and the j-th components is
attractive; if ζij < 0, the interaction is repulsive. Furthermore, to fit
the framework of a continuation method better, we rewrite the m-coupled DNLS equations (2) as
G(x, β) = 0, (6a)
where x = (u>1, . . . , u>m)> ∈ RN m and 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. (6b) We let DG denote the Jacobian matrix of G; in particular,
DG = [Gx, Gβ] ∈ RM ×(M +1),
where M = N m. We define the solution curve of (2) as
C =y(s) = (x(s)>, β(s))>| G(y(s)) = 0, s ∈ R . (7)
Here we assume a parametrization via arc-length s is available on C. By differentiating Equation (6) with respect to s, we obtain
DG(y(s)) ˙y(s) = 0, (8)
where ˙y(s) = ( ˙x(s)>, ˙β(s))> is a tangent vector to C at y(s).
Equation (8) suggests that the tangent vector to C at y(s) is the natural nontrivial solution of the M × (M + 1) homogeneous system DG(y(s))w = 0, when DG(y(s)) is of full row rank. However, if DG(y(s)) is not of full rank at a certain s, then the the continuation algorithm is not well-defined. In this case, typical continuation methods may not follow the solution curve successfully.
The motivation for our hyperplane-constrained continuation method can be illustrated in a straightforward manner by observing the NLS equa-tion for n = 1. Based on the observaequa-tion, we then return to the DNLS equations for n = 2. Afterwards, we discuss how we may circumvent the obstacle for general cases in Section 2.1.
0 0.2 0.4 0.6 0.8 1 β Singular values 0 0.2 0.4 β 0.6 0.8 1 Singular values 10−12 10−10 10−8 10−6 10−4 10−2 100 0 0.2 0.4 β 0.6 0.8 1 Singular values ( )a ( )b ( )c 10−12 10−10 10−8 10−6 10−4 10−2 100 10−12 10−10 10−8 10−6 10−4 10−2 100
Figure 1: The smallest two singular values of Gx vs β with domains (a)
[−5, 5]; (b) [−10, 10] and (c) [−15, 15].
Considering the 2-coupled NLS equations (1) with n = 1 and a fixed β12, we differentiate both sides of (1a) with respect to the spatial variable
z defined in (1c) to obtain Gx(y(s)) φ01 φ02 = L1 2β12φ1φ2 2β12φ1φ2 L2 φ01 φ02 = 0, (9) where L1= d 2 dz2− λ1+ 3µ1φ21+ β12φ22and L2= d 2 dz2− λ2+ 3µ2φ22+ β12φ21.
Equation (9) implies that the matrix Gx(y(s)) is actually singular with
the corresponding singular vector (φ01, φ02)>. Consequently, we may have a one-dimensional solution set with tangent vector (φ01, φ02)> of (1) for a fixed β12
{φr(x)|φr(x) = (φ1(x − r), φ2(x − r)), r ∈ R} (10)
that contains all the translation invariant solutions. To specify this phe-nomena of singularity arising in DNLS equations. We use the case λ1 =
λ2= 1 and β12= β21= β in (2) to illustrate figures regarding the smallest
two singular values of Gxvs β in Figure 1 for the domains [−5, 5], [−10, 10]
and [−15, 15] with grid size h = 0.2. We can see that the second singular value is of O(10−1) which is almost independent of domain size, and the
smallest singular value of Gxdecades very fast as the domain size increases.
This is the reason why we call it “near singularity” in this paper. Denote D the discretization matrix of d
dx. To verify the approximation of the singular
vector of Gxby [(Du1)>, (Du2)>]>, we further illustrate the difference of
the normalized vectors, dx, of [(Du1)>, (Du2)>]> and the singular vector,
v1, of the Jacobian matrix Gx in various domains in Figure 2.
Similarly, singularity can be observed in the discrete version of the m-coupled NLS equations (2) with n = 2. We consider a bounded domain
2.40 2.40 2.40 2.40 2.40 2.40 2.40 2.41 2.41 x 10−4 0 0.2 0.4 0.6 0.8 1 1 2 x d v − ‖‖ ‖ 3.21 0 0.2 0.4 0.6 0.8 1 3.21 3.21 3.21 3.21 3.21 3.21 1 2 x d v − ‖‖ ‖ x 10−2 2.36 2.36 2.36 2.36 2.36 2.36 2.36x 10 −4 0 0.2 0.4 0.6 0.8 1 1 2 x d v − ‖‖ ‖ β β β ( )a ( )b ( )c
Figure 2: The difference of dx and the singular vector of Gx in various
domains: (a) [−5, 5]; (b) [−10, 10] and (c) [−15, 15].
[−d, d] × [−d, d], where the size d is sufficiently large and the uniform grid size h is sufficiently small. Let
D = 1 2h 0 1 0 −1 . .. . .. . .. . .. 1 0 −1 0 ∈ R √ N ×√N
be the central difference operator. We define
b
Dx= D>⊗ I√N, Dby = I√N⊗ D ∈ R N ×N
as the discretization matrices of the differential operators ∂x∂ and ∂y∂ , re-spectively. Then one can check that, similar to the operator dxd in (9),
Ks= span{Dxx(s), Dyx(s)} (11)
forms a numerical null space of Gx(x(s), β(s)), where
Dx= diag{ bDx, . . . , bDx}, Dy= diag{ bDy, . . . , bDy} ∈ RN m×N m. (12)
In other words, Gx(x(s), β(s)) is nearly singular and therefore may result
in numerical difficulties in practice.
If we use p ≡ a[(Dxx(s))>, 0]>+ b[(Dyx(s))>, 0]> (a2+ b2= 1) as the
prediction direction, then, numerically, the “new solution” x(˜s) (˜s ≈ s) computed by the continuation method is a small shift of x(s) in the p-direction. These new solutions generally have very small residual
( )
a
( )
b
Figure 3: Conceptual illustration of (a) translation invariant solutions and (b) ε-solutions.
kG(x(s), β(s))k < ε and are called “ε-solutions”. Note that as the domain of (2) is bounded, the residual gets larger if we keep shifting x(s) along the p-direction, where (x(s), β(s)) is the solution of G(y(s)) = 0.
The above arguments pose the following two difficulties that may be encountered by a standard continuation method. First, the near singular-ity may cause accuracy and efficiency problems while solving the resulting linear system. Second, a continuation method may be trapped, or become hard to keep moving ahead correctly, due to the near singularity. Figure 3 further conceptually illustrates the effects of ε-solutions. The solid curves in Part (a) of the figure represent translation invariant solutions of the NLS equations on a unbounded domain. One example is defined in (10). In con-trast, no such translation invariant solution exists in exact arithmetics on a bounded computational domain. However, in a standard continuation method, any of the ε-solutions (dashed curves in part (b)) may be consid-ered to be a reasonable approximation of the target solution (solid curve) and a particular ε-solution would be accepted as a new solution. Here, the target solution is the one that is supposed to be located on a certain solution curve being tracked. Unfortunately, the chosen solution is actu-ally the exact solution of another solution curve. Therefore, the computed solutions may simply jump between different solution curves, rather than follow a certain solution curve.
We now summarize that such ε-solutions result in the following chal-lenges to the prediction-correction scheme of a standard continuation method. C1. In the prediction step, we cannot compute a unique prediction
direc-tion by solving (8).
C2. In the correction step, since the Jacobian matrix Gxis nearly singular,
Newton’s correction will lose quadratic convergence and accuracy. C3. The bifurcation points are difficult to detect due to the numerical
singularity of the Jacobian matrix.
C4. Since the solution manifold C in (7) contains a 2-dimensional ε-solution set and the “good” prediction direction can not be uniquely com-puted, the solutions computed by the continuation method may ap-pear to be random or trapped in the ε-solution set. That is, we can not follow the desired solution curve C in (7) efficiently by a standard continuation method.
To overcome these difficulties, we develop a hyperplane-constrained con-tinuation method for solving (6) in the following subsections.
2.1
Prediction and correction with hyperplane
con-straints
Let yi = (x>i , βi)> ∈ RM +1 be a point that has been accepted as an
ap-proximation point for the solution curve C. A “good” prediction direction ˙
yi = ( ˙x>i , ˙βi)> should satisfy (8) and the vector ˙xi at the arc length
pa-rameter s = si should be in Ks⊥i, where Ksi is given in (11). It follows
that the prediction direction ˙yi∈ RM +1should satisfy the bordered linear
system Gx Gβ a>x 0 a>y 0 c>i ci ˙ yi = 0 0 0 1 , (13)
where ax= Dxxi, ay= Dyxi, and (c>i , ci)>∈ RM +1is a suitable constant
vector. In other words, we first use the Euler predictor yi+1,1= yi+ hiy˙i
to predict a new point yi+1,1, where hi> 0 is the step length and ˙yiis the
unit tangent vector at yithat is obtained by normalizing the solution of the
bordered linear system (13). Note that Equation (13) can be interpreted geometrically as follows. The next solution yi+1 must pass through the
two hyperplanes whose normal vectors are ax and ay. Furthermore, we
accordingly coin the name “hyperplane-constrained continuation method” due to these two additional hyperplane constraints.
Now the solution curve C is determined by the underlying system of equations G(x, β) = 0, a>xx = 0, a>yx = 0, ˙x>i x + ˙βiβ = 0.
Starting from the predictor, the accuracy of the approximation yi+1,1 to
the solution curve C can be improved by a correction process. Typically, Newton’s method is chosen as a corrector. By setting yi+1,l+1= yi+1,l+ δl
for l = 1, 2, . . ., we solve the bordered linear 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 , (14)
with ρx,l = a>xxi+1,l, ρy,l = a>yxi+1,l, and ρl = ˙y>i yi+1,l. If {yi+1,l}
con-verges until l = l∞, then we accept yi+1= yi+1,l∞ as a new approximation
to the solution curve C.
The two linear systems (13) and (14) are overdetermined systems and can be solved by the least squares method with very small minimal residual. Another efficient way to solve (13) and (14) is to rewrite them in the form
B f g> γ x σ = q ρ , x(M + 1) = x(M + 2) = 0, (15) where B = Gx ax ay a> x 0 0 a> y 0 0 ∈ R(M +2)×(M +2) (16)
is symmetric and f , g, q ∈ R(M +2). The linear system (15) can be
eas-ily solved with the well-known block elimination (BE) algorithm (see e.g. [16]) when B is well-conditioned. However, near turning points or branch points, B in (15) becomes nearly singular, i.e., B is ill-conditioned. In this case, the linear system should be solved by the deflated block elimination (DBE) algorithm by Chan [6], or the more efficient and backward stable, mixed block elimination (MBE) algorithm proposed by Govaerts [11, 10] as shown below.
Algorithm 1 [Mixed Block Elimination]. (i) Solve ξ>B = g>,
(ii) Compute δ1= γ − ξ>f , σ = (ρ − ξ>q)/δ1,
(iii) Solve Bv = f ,
(v) Solve Bw = q1,
(vi) Compute σ1= (ρ1− ξ>w)/δ,
(vii) Compute x = w − vσ1, σ = σ + σ1.
We finally make the following remarks regarding Algorithm 1. (i) The main step in Algorithm 1 is to solve the linear system of the form Bξ = g. (ii) Since linear systems (13) and (14) have very small minimal residuals, Equation (15) is nearly consistent. Thus, the solution x solved by Algo-rithm 1 satisfies x(M + 1) ≈ x(M + 2) ≈ 0 automatically. In the next section, we discuss how we deal with the numerical singularity of the Ja-cobian matrix to find the bifurcation points.
2.2
Testing for Bifurcation
Let C be the solution curve defined in (7), y(s) ∈ C and
J(s) = Gx(y(s)) Gβ(y(s)) a> x 0 a>y 0 ∈ R(M +2)×(M +1), (17a) J(s) = Gx(y(s)) a>x a>y ∈ R (M +2)×M. (17b)
As described in [2, 12, 16], a point y(s) ∈ C is said to be a regular point if rank(J(s)) = M (i.e., dim N (J(s)) = 1) and is a singular point if rank(J(s)) ≤ M − 1 (i.e., dim N (J(s)) ≥ 2). For a regular point y(s), the tangent vector ˙y(s) is uniquely determined by the linear system (13).
Now, our task is to design an algorithm to detect singular points of the solution curve C and to compute tangent vectors if y(s) is a singular point. For simplicity, we here only consider the case
(Gβ(y(s))>, 0, 0)>∈ R(J(s)) for each singular points y(s) ∈ C, (18)
i.e., dim N (J(s)) = dim N (J(s)) − 1 and dim N (J(s)) ≥ 2. This case shows that the tangent vector at a singular point has a nonzero component at ˙β(s) and can be expected to appear in the solution curve C of (6). We denote B(s) as the matrix B given by (15) at the point y(s) ∈ C.
Our strategy for detecting the singularity s∗of C is summarized in Algo-rithm 2. Let s1< s2 be two consecutive continuation method parameters
and µ(s1) and µ(s2) be the smallest eigenvalues in modulus of B(s1) and
B(s2), respectively. It is clear that if µ(s1) > 0 and µ(s2) < 0, then there
refine the interval (s1, s2). In the secant method loop, we use the inverse
power method (in Step (ii-c)) to compute the smallest eigenvalue. After convergence, we use Algorithm 3 described in the appendix to compute the tangent vectors at the singularity in Step (iii).
Algorithm 2 Detection of Singularity of C.
(i) Given µ(si) the smallest eigenvalue in modulus of B(si), i = 1, 2,
where µ(s1) > 0, µ(s2) < 0.
(ii) Perform the Secant Method until convergence,
(a) Compute y1(s0) := y(s0) = y(s1) +
tµ(s1)
µ(s2) − µ(s1)
, where t = y(s1) − y(s2).
(b) Perform the Newton Correction (14) 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, Go to (b).
(c) Compute µ(s0) of B(s0) with y(s0) = y`∞(s0) by using the
inverse power method.
(d) If |µ(s0)| <Tol, then perform (iii).
(e) If µ(s0) > 0, s1← s0, else s2← s0, go to (ii).
(iii) Compute the desired tangent vectors with y(s∗) = y`∞(s0) by the
methods suggested in [16, p.88-99].
It is worth noting that detecting the singular point of the solution curve C is equivalent to detecting the singularity of symmetric matrix B in (16). It can be shown that, provided the condition (18) holds, then the following statements are equivalent: (i) rank(J(s)) ≤ M − 1, (ii) N (J(s)) 6= {0}, (iii) B(s) is singular and there is a nonzero vector χ ∈ RM such that
(χ>, 0, 0)> ∈ N (B(s)). Therefore, in Step (ii-c) of Algorithm 2, instead of checking singularity of J(s) in (17a), we only need to check the singularity of the square symmetric matrix B(s) in (16).
We have proposed a hyperplane-constrained continuation method that deals with the numerical challenges C1 to C4 listed on page 8. Specifically,
2. The efficiency of Newton’s correction is improved, thanks to the bet-ter conditioned Jacobian matrix (14).
3. The bifurcation points can be detected accurately from (17). 4. The enhancements listed in the first two points lead to the
continua-tion method that is capable of following the desired solucontinua-tion curves. In the next section, we focus on the 3-coupled cases in both theoretical and numerical aspects. The results not only characterize the solutions of DNLS equations analytically, but demonstrate the bifurcation diagrams and visualize the theoretical predictions in Sections 3.
3
Bifurcation analysis for 3-coupled discrete
nonlinear Schr¨
odinger equations
In this section, we study the 3-coupled DNLS equations theoretically by determining the corresponding primal stalk solution curve and conducting a bifurcation analysis. Note that there is no bifurcation for 1-component DNLS equations due to the uniqueness of the positive solution [?]. Further-more, for 2-coupled DNLS equations, bifurcation analysis has been studied by Kuo, Lin, and Shieh in [?].
Lin and Wei [17] have analyzed the NLS equations (1) and the corre-sponding ground state solutions. Denoting βij = ζijβ (see (5)) and letting
Σ = 1 |β12| |β13| |β12| 1 |β23| |β13| |β23| 1 ,
some of their results regarding the 3-coupled NLS equations are categorized as follows.
Case 1 (all interactions are repulsive). If ζ12< 0, ζ13< 0 and ζ23<
0, then the ground state solution does not exist.
Case 2 (all interactions are attractive). If ζ12 > 0, ζ13 > 0, ζ23 > 0
and Σ is positive definite, then the ground state solution exists. Case 3 (two repulsive and one attractive interactions). If ζ12< 0,
ζ13 < 0, ζ23 > 0 and Σ is positive definite, then the ground state
solution does not exist.
Case 4 (two attractive and one repulsive interactions). If ζ12> 0,
ζ13> 0, ζ23< 0, β 1 and the ground state solution exists, then it
Now we use the same categories of ζij’s and consider the solution curve
C of (2) by letting
m = 3 and λ1= λ2= λ3= µ1= µ2= µ3= 1.
In such a case, the 3-coupled DNLS equations in (6) become Au1− u1+ u 3 1 + βζ12u 2 2 u1+ βζ13u 2 3 u1= 0, (19a) Au2− u2+ u 3 2 + βζ12u 2 1 u2+ βζ23u 2 3 u2= 0, (19b) Au3− u3+ u 33 + βζ13u1 2u3+ βζ23u 22u3= 0. (19c)
As Cases 1 and 2 are straightforward, we focus on the following two par-ticular settings of Cases 3 and 4:
ζ12= ζ13= −1, ζ23= 1, (20a)
ζ12= ζ13= 1, ζ23= −1. (20b)
In (20a), the 3-coupled DNLS equations of (6) become Au1− u1+ u 3 1 − βu 2 2 u1− βu 2 3 u1= 0, (21a) Au2− u2+ u 3 2 − βu 2 1 u2+ βu 2 3 u2= 0, (21b) Au3− u3+ u 33 − βu 2 1 u3+ βu 22u3= 0, (21c)
where β > 0. It is clear that if we set β := −β, then (21) describes the 3-coupled DNLS equations for the case (20b). Therefore, to investigate these two cases, we only need to consider 3-coupled DNLS equations (21) for β ∈ R.
In the following two theorems, we first explicitly determine the solutions located on the primal stalk and then discuss how other solution curves bifurcate from the primal stalk.
Theorem 1 The primal stalk of the solution curve
C = {y(s) = (x>(s), β(s))>| G(y(s)) = 0 is given in (21) and s ∈ R}, (22) for −1
3 ≤ β < 1, has the forms
u1= s 1 + 3β 1 + β − 2β2 ! u∗ and u2= u3= s 1 + β 1 + β − 2β2 ! u∗, (23)
where x(s) = (u>1(s), u>2(s), u>3(s))> and u∗ is the positive solution of
Proof. By letting
u2= u3= κu1 with κ > 0, (25)
it follows that equations (21b) and (21c) are identical. Thus, equations in (21) can be reduced to
Au1− u1+ (1 − 2βκ2)u 13 = 0,
Au1− u1+ (κ2− β + βκ2)u 13 = 0.
(26)
The system of equations in (26) has a positive solution u1, if 1 − 2βκ2=
κ2− β + βκ2. This implies that
κ = s
1 + β
1 + 3β. (27)
Substituting κ in (27) into the first equation of (26) we have
Au1− u1+
1 + β − 2β2
1 + 3β u
3
1 = 0. (28)
It can be easily verified that if u∗ is a solution of Au − u + u 3 = 0, then
u1=
s
1 + 3β
1 + β − 2β2u∗≡ ηu∗ (29)
solves (28).
The applicable range of β is determined by the following facts. Since κ → q 1 2 as β → 1 − (by (27)), we have u2= u3= κu1→ ∞ (30)
by (25) and (29). On the other hand, since
u1→ 0 as β → − 1 3 − (31) by (29), we have u2= u3→ r 3 2u∗ by (21b) and (21c).
Theorem 2 The primal stalk described by (23) 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
∗ ]] and u∗is the
Proof. Since (21) has a positive solution curve u2(β) = u3(β) = κu1(β),
for 0 < β < 1, where κ is defined in (27), the Jacobian matrix of (21) with respect to u is of the form
Gu(y(β)) = B1 E1 E1 E1 B2 E2 E1 E2 B2 (32) where B1= A − I + [[3u 12 − 2βu 2 2 ]], (33a) B2= A − I + [[(3 + β)u 22 − βu 2 1 ]], (33b) E1= −2β[[u1◦ u2]], (33c) E2= 2β[[u 22]]. (33d)
From (25), (27) and (29), we have
u2= u3= s 1 + β 1 + 3βu1= s 1 + β 1 + β − 2β2u∗, (34)
where u∗ is the positive solution of (24). Substituting (29) and (34) into
(33), we get B1= A − I + 3 + 7β − 2β2 1 + β − 2β2 [[u 2 ∗ ]], (35a) B2= A − I + 3 + 3β − 2β2 1 + β − 2β2 [[u 2 ∗ ]], (35b) E1= −2β p(1 + β)(1 + 3β) 1 + β − 2β2 [[u 2 ∗ ]], (35c) E2= 2β + 2β2 1 + β − 2β2[[u 2 ∗ ]]. (35d) 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 . (36) Hence σ(Gu(y(β))) = σ B1 E1 2E1 B2+ E2 ∪ σ(B2− E2). (37)
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 1 (2β+1)κη2 ⊗ 2βκη2[[u 2 ∗ ]] = 1 0 0 1 ⊗ B1− 0 1 2 a ⊗ 2βκη2[[u 2 ∗ ]] (38)
where κ and η are given by (27) and (29), respectively, and
a = 1 (2β + 1)κη2 = 1 − β p(1 + 3β)(1 + β). (39) Since a+√a2+8 2 and a−√a2+8 2 are eigenvalues of 0 1 2 a , from (35), (37), and (38), it follows that
σ(Gu(y(β))) = Λ1(β) ∪ Λ2(β) ∪ Λ3(β), (40) where Λ1(β) = σ A − I +4β + 3 2β + 1[[u 2 ∗ ]] , (41a) Λ2(β) = σ A − I + 3 + 7β − 2β 2 1 + β − 2β2 − (a + p a2+ 8)βκη2 [[u 2 ∗ ]] , (41b) Λ3(β) = σ A − I + 3 + 7β − 2β 2 1 + β − 2β2 − (a − p a2+ 8)βκη2 [[u 2 ∗ ]] . (41c) Since (a +pa2+ 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 ,
it holds that
Λ2(β) = σ(A − I + 3[[u ∗2]]). (42)
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]]).
(43) As β → 1−, we have that Λ1(β) → σ(A − I + 7 3[[u 2 ∗ ]]) (44) and 3+7β−2β1+β−2β22 − (a − √
a2+ 8)βκη2 → ∞. So there exists a β∗ with 0 <
β∗< 1 such that
Λ3(β) ⊂ R+, for β > β∗. (45)
If the number of nonnegative eigenvalues of Λ3(0) = σ(A−I+3[[u ∗2]]) is
p, then from (42)-(45) we see that the primal stalk of the solution curve C of (22) undergoes at least N − p bifurcation points at finite values 0 < β∗q < 1, q = 1, . . . , N − p.
4
Numerical results
In this section, we study numerical results of positive bound state solutions for 3-coupled DNLS equations (19) with λ1= λ2= λ3= µ1= µ2= µ3= 1
by using the hyperplane-constrained continuation method developed in Sec-tion 2. The initial point on the primal stalk of the soluSec-tion curve C of (19) is computed by the fixed point iteration method described [?]. A squared domain [−5, 5] × [−5, 5] with the grid size h = 0.2 is used in the compu-tations. The results, including solution profiles, bifurcation diagrams, and corresponding energies, are summarized in the following two simulations.
4.1
Simulation 1
We consider the case in which ζ12 = ζ13 = −1 and ζ23 = 1, as described
in (20a). That is, we assume that one attractive and two repulsive inter-actions occur among the components. This setting corresponds to Case 3 in Section 3. The positive bound state solutions of (19) are computed and denoted by the set of solution curves
C+ =(x>, β)>| G(x, β) = 0 is given in (19) with (46)
0.237 0.253 0.3850.387 β 1 3N
A
B
C
D
F
E
Figure 4: Solution profiles and a bifurcation diagram of the so-lution curves C+ = (x>, β)>| G(x, β) = 0 is given in (19) with
ζ12= ζ13= −1, ζ23= 1, for β ∈ R+}.
In Figure 4, we plot the bifurcation diagram of the (conceptual) solu-tion curves for β ∈ [0, 1.2]. The intersecsolu-tion points of the solusolu-tion curves indicate the bifurcation points. The primal stalk solution curve, which is also described analytically in (23), is plotted by the solid red curve. The curve starts from β = 0 and bifurcates at β = 0.237 and 0.385. The so-lution curves bifurcated from the primal stalk while β = 0.237 are plotted as dashed blue curves. This solution curve bifurcates again at β = 0.253. On the other hand, the solution curves bifurcated from the primal stalk while β = 0.385 are plotted as dotted green curves. This solution curve bifurcates again at β = 0.387.
Furthermore, we characterize the solution curves by showing the corre-sponding nodal domains of three positive bound state solutions. The nodal domains are attached in triples near the solution curves. Since the form of the solutions remain similar unless bifurcation occur, we simply show one representative nodal domain triple for each of solution curves. Note that in each of the nodal domain triples, the left, middle and right figures are the density plots of u1, u2 and u3, respectively. In particular, triple
A corresponds to the primal stalk (solid red curve); triples B, C, and D are associated with the dashed blue solution curves; triples E and F are associated with the dotted green solution curves.
The energy defined in (4) is computed for all solutions. In Figure 5, we plot the energy curves corresponding to the solution curves for β ∈ [0, 1.2] by using the same curve styles to indicate the corresponding solution curves. For example, the energies corresponding to the primal stalk solutions are plotted as the solid red curve. Furthermore, the nodal domains of u1,
1 15 20 25 30 35 40
(
)
E
x
A
B
C
D
E
F
0.237 0.253 0.385 0.387 β 0Figure 5: Energy curves of C+.
u2 and u3 (in the form of squared sums) are plotted in an overlapping
format to show their relative positions. These overlapping nodal domains are labeled from A to F to indicate their corresponding triples presented in Figure 4.
4.2
Simulation 2
We consider the case in which ζ12= ζ13= 1 and ζ23= −1, as described in
(20b). That is, we assume that one repulsive and two attractive interactions occur among the components (i.e. Case 4 in Section 3). The positive bound state solutions are computed and denoted by the set of solution curves
C− =(x>, β)>| G(x, β) = 0 is given in (19) with (47)
ζ12= ζ13= 1 and ζ23= −1, for β ∈ R+} .
We plot the bifurcation diagram and the energy curve in Figures 6 and 7, respectively, for β ∈ [0, 0.5]. Similar to the results reported in Simulation 1, nodal domains of positive bound state solutions are attached near the solution curve and the energy curve. Note that there is no bifurcation in the solution curve C−, but a turning point is found at β = 0.333.
4.3
Remarks on the two simulations
We highlight the following observations from Simulations 1 and 2 that are consistent with the solution characters obtained from theoretical analysis.
β 3N
A
B
0.333Figure 6: Solution profiles and a bifurcation diagram of the so-lution curves C− = (x>, β)>| G(x, β) = 0 is given in (19) with ζ12= ζ13= 1 and ζ23= −1, for β ∈ R+}. 0 15.5 16 16.5 17 17.5
(
)
E
x
βA
B
0.333• In Figure 4 of Simulation 1, the β’s in the primal stalk approach, but never reach, β = 1. In contrast, Figure 6 of Simulation 2 shows that a turning point is observed on the primal stalk at β = 0.333, which corresponds to the case of β = −13 in Theorem 1.
• As shown in Figure 5, the primal stalk energy curve (plotted as solid red) keeps rising as β increases. This phenomenon is consistent with (30).
• The computed solution profiles of C− not only have the property as
shown in (31), Figure 6 further shows that the u1 turns to negative
side after passing the turning point.
• As Theorem 2 discusses the number of bifurcation points, we note that it is proved in [17, Lemma 1] that the number of nonnegative eigenvalues of
4φ − φ + 3ω2
∗φ = λφ,
φ ∈ H2(Rn), is n + 1, where ω∗ is the unique solution of
4φ − φ + φ3= 0, φ > 0 in Rn, ω(x) → 0 as |x| → ∞.
In the 3-coupled DNLS equations (21) with a squared domain (n = 2), we can verify numerically that the number of nonnegative eigenvalues of Λ1(0) = σ(A − I + 3[[u ∗2]]) is 3, where u∗ is the positive solution
of (24).
5
Conclusion
Aiming at the coupled DNLS equations that are discretized from the NLS equations on the whole domain, we have developed a new hyperplane-constrained continuation method. The method overcomes the numerical difficulties that prevent standard continuation methods from working or from being efficient due to the ε-solutions. We have also analyzed the primal stalk solution curve for the 3-coupled DNLS equations and we have demonstrated numerical results showing the versatility of the bound state solutions.
Acknowledgments
This work is partially supported by the National Science Council, the Na-tional Center for Theoretical Sciences in Taiwan, and the Taida Institute of Mathematical Sciences.
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. Handbook of Nemerical Analysis, Edited by P. G. , 5, 1996.
[3] W. Z. Bao. Ground states and dynamics of multi-component Bose-Einstein condensates. SIAM Multiscale Modeling and Simulation, 2(2):210–236, 2004.
[4] 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.
[5] M. Berry. Large scale singular value computations. Internat. J. Su-percomputer Appl., 6(1):13–49, 1992.
[6] T. F. Chan. Deflation techniques and block-elimination algorithm for solving bordered singular systems. SIAM J. Sci. Stat. Compu., 5:121– 134, 1984.
[7] 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.
[8] 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.
[9] K. W. Chow. Periodic solutions for a system of four coupled nonlinear Schr¨odinger equations. Phys. Lett. A, 285:319–326, 2001.
[10] W. Govaerts. Stable solvers and block elimination for bordered sys-tems. SIAM J. Matrix Anal. Appl., 12:469–483, 1991.
[11] W. Govaerts and J.D. Pryce. Mixed block elimination for linear sys-tems with wider borders. IMA J. Numer. Anal., 13:161–180, 1993.
[12] Willy J. F. Govaerts. Numerical Methods for Bifurcations of Dynam-ical Equilibria Computations. SIAM, Philadelphia, 2000.
[13] F. T. Hioe. Solitary waves for n coupled nonlinear Schr¨odinger equa-tions. Phys. Rev. Lett., 82:1152–1155, 1999.
[14] 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.
[15] T. Kanna and M. Lakshmanan. Exact soliton solutions, shape chang-ing collisions, and partially coherent solitons in coupled nonlinear Schr¨odinger equations. Phys. Rev. Lett., 86:5043–5046, 2001.
[16] H. B. Keller. Lectures on Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin, 1987.
[17] 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.
[18] M. Mitchell, Z. Chen, M. Shih, and M. Segev. Self-trapping of partially spatially incoherent light. Phys. Rev. Lett., 77:490–493, 1996.
Appendix
Algorithm 3 Tangent Vectors at Singularity. (I) For dimN (J(s∗)) = 1 :
(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 [5] (or another suitable package); (ii) Form φ1= ¯ φ 0 and φ2= ¯ φ0 1 ;
(iii) Solve the real vector roots {(ˆµk, ˆνk)}2k=1 of a11µ2+ 2a12µν +
a22ν2with 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 [5] (or another suitable package);
(ii) Form φk= φ¯ (k) 0 ! , k = 1, . . . , ` and φ`+1= ¯ φ0 1 ;
(iii) Form trial tangent vectors ˙yk(s∗) = φk, k = 1, . . . , ` and ˙y`+1(s∗) =
