行政院國家科學委員會專題研究計畫 成果報告
格若士-比塔烏斯基方程式的行進波解之穩定性研究
研究成果報告(精簡版)
計 畫 類 別 ： 個別型 計 畫 編 號 ： NSC 98-2115-M-009-007- 執 行 期 間 ： 98 年 08 月 01 日至 99 年 07 月 31 日 執 行 單 位 ： 國立交通大學應用數學系（所） 計 畫 主 持 人 ： 張書銘 計畫參與人員： 碩士班研究生-兼任助理人員：林育賢 碩士班研究生-兼任助理人員：陳勛暉 碩士班研究生-兼任助理人員：段俊旭 碩士班研究生-兼任助理人員：羅健峰 大專生-兼任助理人員：賴又菱 報 告 附 件 ： 出席國際會議研究心得報告及發表論文 處 理 方 式 ： 本計畫可公開查詢中 華 民 國 99 年 08 月 24 日
1
Introduction
In this project, we consider the nonlinear Schr¨odinger equation (NLS) with focusing pa-rameter ,
i∂tψ = −∆ψ + V (x)ψ − m(x)|ψ|2ψ,
where ψ(t, x) : R × R2 → C is the wavefunction, V (x) is the trap potential, m(x) is the
mass and s-wave scattering length term, and > 0. We would like to consider two kinds of trap potentials and domains in this problem.
First, we introduce nonlinear Schr¨odinger equations and the two kinds of trap poten-tials first.
1.1
Nonlinear Schr¨
odinger Equations
The Schr¨odinger equation is a basic equation in quantum theory. The study of Schr¨odinger equations plays an important role in modern physics. In 1926, an Austrian physicist Erwin Schr¨odinger constructed the Schr¨odinger equation for explaining the active of particles. The famous Schr¨odinger equation is the form: Hψ = i~∂ψ∂t, where H is a Hamiltonian operator, ψ is the wavefunction, ~ represents Planck’s constant over 2π, i is the imagi-nary unit. The Schr¨odinger equation is an equation to describe the possible distribution of atoms. Solutions to Schr¨odinger’s equation describe not only atomic and subatomic systems, electrons and atoms, but also macroscopic systems, possibly even the whole universe.
The nonlinear Schr¨odinger equation (NLS) is a nonlinear version of Schr¨odinger’s equation in theoretical physics. In mathematical point of view, the NLS equation is a Schr¨odinger equation with nonlinear term. The nonlinear Schr¨odinger equation is the partial differential equation
i∂tψ = − 1 2∂ 2 xψ + c|ψ| 2ψ
for the complex field ψ. There are many kinds of nonlinear Schr¨odinger equations. One general form of such equations [12] would be
i∂tψ + ∆ψ = f (ψ, ¯ψ).
These equations (particularly the cubic NLS equation) arise as model equations from several areas of physics: nonlinear optics, quantum condensates.
Here, we review some properties of NLS equations. Let’s consider the equation i∂tψ + ∆ψ = ±|ψ|2ψ.
The sign “−” on the right hand side is focusing nonlinearity, and the sign “+” is defo-cusing. When we add the potential part to the NLS equation, then it becomes
i∂tψ + ∆ψ = −m|ψ|p−1ψ + V ψ, (1)
where V is real and time independent. The equation (1) is sometimes referred to as a Gross-Pitaevskii equation. We can look at the dynamics of the Bose-Einstein condensate from the time-dependent Gross-Pitaevskii equation.
The behavior of ψ is determined by the potential V of the NLS equation. There is a special case of potential V = ±|x|2, this can be used to model a confining magnetic trap
1.2
Trap Potential V
In the previous subsection, we have introduced some property of NLS equation. Now, we introduce the trap potential. Trap potential is to trap atoms on the minima potential. The following are two forms of trap potentials in this study, which are optical lattices and elliptic form.
1.2.1 Optical Lattices
An optical lattice is essentially an artificial crystal of light. A periodic intensity pattern that is formed by the interference of two or more laser beams. The shape of optical lattices looks like an egg carton in Fig. 1. It is called an optical lattice, since the periodic arrangement of trapping sites resembles a crystalline lattice. In Fig. 1, atoms are cooled and congregate in the minimal potential. The well depth and the periodicity are two important parts to affect the potential shape.
Figure 1: Optical Lattices.
Besides trapping cold atoms, optical lattices are also use in sorting microscopic parti-cles [11] recently.
1.2.2 Elliptic Potential
Let d × d matrix M be symmetric, positive define. Let v ∈ Rd and c ∈ R be arbitrary.
The elliptic potential [2] defines as V (u) = 1
2u
T
M u + uTv + c,
where u is a column vector in Rd. Because M is symmetric and positive definite, M12
exists. Therefore, we can also write V (u) = 12||M12u||2 + uTv + c. In R2, the elliptic
potential can be rewritten as
V (x, y) = ax2+ bxy + cy2+ dx + ey + f, a, b, c, d, e, f ∈ R with b2− 4ac < 0. The shape of this potential looks like a bowl.
1.3
Goals
In this project, we consider nonlinear Schr¨odinger equations with focusing parameter , i∂tψ = −∆ψ + V (x)ψ − m(x)|ψ|2ψ, (2)
where ψ(t, x) : R × R2 → C is the wavefunction, V (x) is the trap potential, m(x) is
the mass and s-wave scattering length term, and > 0. Such equation occurs in many physics, including nonlinear optics, quantum physics, and water waves.
In this project, we consider two kinds of different trap potentials and domains as follows:
Case 1. (ω case)
V (x) = (ω1x1)2+ (ω2x2)2,
where ω1, ω2 ≥ 0 and ω1, ω2 = 1, 2, ..., 5, and m(x) = 1 + 12 sin(2πx1) with unit disk
domain Ω = {(x1, x2)|x21+ x22 < 1}.
Case 2. (µ case)
V (x) = 1 + sin2(µ1x1) + sin2(µ2x2),
where µ1, µ2 = 1, 2, ..., 5 and m(x) = 1 + 12cos(2πx1) with square domain Ω =
{(x1, x2)| − 1 < x1 < 1, −1 < x2 < 1}.
Here, we have two kinds of trap potentials and function m(x) in Fig. 3 and Fig. 2, respectively. We focus on different situation in ω case and µ case. We concentrate on the well depth of the potential by varying the parameters ω1 and ω2 for the elliptic trap
potential. And then we control the periodic for the optical lattices trap potential by changing the two parameters, µ1 and µ2. For keeping the shapes of trap potentials, we
use different types of domains for ω case and µ case.
The NLS equation has a special solution: ψ(t, x) = eiλtφ(x). The aim of this project
is to find out the spectra of linearized operator which arises when the equation (2) is linearized around the special solution ψ. The goal is to study the stability of the solution by the spectrum of its linear operator.
There are many literature on NLS equation. The necessary condition for orbital sta-bility and instasta-bility of single-spike bound state can be obtained in [9, 10].
Next, in Section 2 we will present the mathematical analysis about the NLS equation. And then in Section 3 we will discuss the numerical method. At last, some numerical results will be shown in Section 4.
2
Mathematical Analysis
2.1
Perturbation in NLS
Recall the nonlinear Schr¨odinger equation:
i∂tψ = −∆ψ + V (x)ψ − m(x)|ψ|2ψ. (3)
We consider the special solutions of NLS equation (3): ψ(t, x) = eiλtφ(x), where λ =
0, 1, and general case. φ(x) is a real-valued function independent of time. Substitution ψ(t, x) = eiλtφ(x) into (3) and the φ satisfy the nonlinear elliptic equation:
that is
(−∆ + V − m|φ|2)φ = −λφ (4) or
(−∆ + (V + λ) − m|φ|2)φ = 0. (5) Equation (4) and (5) will be used in different solution form and algorithm.
To study the stability of the special solution form, we consider solutions of the NLS equation of the form
ψ(t, x) = eiλt(φ(x) + h(t, x)). (6) The perturbation h(t, x) ∈ R satisfies an equation:
∂th = Lh + (nonlinear terms),
where
Lh = −i{(−∆ + (V + λ) − 3mφ2)h}. (7)
Remark. To explain how to obtain (7):
in equation (3), replace ψ(x) by eiλt(ψ(x) + h(t, x)) then the left hand side of (3):
i∂tψ = i∂t(eiλt(φ + h))
= i(iλeiλt(φ + h) + eiλt∂th)
= eiλt(−λ(φ + h) + i∂th).
And the right hand side of (3):
−∆ψ + V (x)ψ − m(x)|ψ|2ψ = eiλt(−∆φ − ∆h + V (x)ψ + hφ − m|φ + h|2(φ + h))
= eiλt{−∆φ + V φ − ∆h + hφ
−m(|φ|2φ + φ2(¯h + h) + φh¯h + φ2h + φ(h + ¯h)h + h2¯h)}
= eiλt(−∆φ + V φ − m|φ|2φ)
+eiλt(−∆h + V h − m|φ|2(h + ¯h) − m|φ|2h) + eiλtO(h2). Notice that |φ + h|2 = (φ + h)(φ + h) = |φ|2+ φ(¯h + h) + h¯h. Then
eiλt(−λ(φ + h) + i∂th) = eiλt(−∆φ + V φ − m|φ|2φ)
+eiλt(−∆h + V h − m|φ|2(h + ¯h) − m|φ|2h) + O(h2). So
i∂th = (−∆φ + (V + λ)φ − m|φ|2φ)
+(−∆h + (V + λ)h − m|φ|2(h + ¯h) − m|φ|2h) + O(h2) = −∆h + (V + λ)h − m|φ|2(h + ¯h) − m|φ|2h + O(h2). If we consider h(t, x) as a complex perturbation, ie. h(t, x) ∈ C, then Lh = L(Re h + i Im h)
= −i{(−∆ + (V + λ) − mφ2) Re h + i(−∆ + (V + λ) − mφ2) Im h − 2mφ2 Re h} = (−∆ + (V + λ) − mφ2) Im h + i{−(−∆ + (V + λ) − 3mφ2) Re h}.
And we rewrite L as a matrix acting on Re h Im h , L = 0 L− −L+ 0 , (8) where L+ = −∆ + (V + λ) − 3mφ2, L−= −∆ + (V + λ) − mφ2. (9)
L+ and L− are self-adjoint.
In the following Lemma, we will explain how we could rewrite the linerization term as (8). Are eigenvalues and eigenvectors of (8) and the original form the same or not? Lemma 1. Let L = 0 L− −L+ 0 and Lh = −∆h + (V + λ)h − m|φ|2(h + ¯h) − m|φ|2h, where L
+ and L− are defined by (9). If ρ is an eigenvalue of L and [uT, vT]T
be the eigenvector, then ρ is also an eigenvalue of L and u + iv is an eigenvector of L corresponding to ρ.
Proof. Since ρ is an eigenvalue of L, that is, L u w = ρ u w . (10) So 0 L− −L+ 0 u w = ρ u w . Then −(−∆ + (V + λ) − 3mφ2)u = ρw (11) and (−∆ + (V + λ) − mφ2)w = ρu. (12) (12) + i×(11), therefore
ρ(u + iw) = −i(−∆ + (V + λ) − 3mq2)u + (−∆ + (V + λ) − mq2)w = −i{(−∆ + (V + λ) − 3mq2)u + i(−∆ + (V + λ) − mq2)w} = −i{(−∆(u + iw) + (V + λ)(u + iw) − mq2(u + iw) − 2mq2u}. Hence u + iw is an eigenvactor of L corresponding to ρ.
2.2
Stability
Our aim of this project is to study the stability of these solutions. We show the stability by solving eigenvalue problem. So in this subsection, we recall the definition of some stability, and the relation between eigenvalues and stability [4, 5].
At first, we give some notations and definitions. Consider an ODE dx
dt = ˙x = f (x); x ∈ R
n, (13)
where f : Rn → Rn and x = [x
1, x2, ..., xn]T. If f (x∗) = 0 for all t, the point x∗ is called
an equilibrium point.
Definition 1. x∗ is a Lyapunov stable equilibrium if for every neighborhood U of x∗ there is a neighborhood V ⊆ U of x∗ such that every solution x(t) with x(0) = x0 ∈ V is defined
and remains in U for all t ≥ 0.
Definition 2. If V can be chosen above so that, in addition to the properties for stability, we have limt→∞x(t) = x∗ then we say that x∗ is asymptotically stable.
An equilibrium is called neutrally stable if it is Lyapunov stable but not asymptotically stable.
In studying the stability of x∗, we consider x∗ plus a small perturbation h(t), ie, x(t) = x∗+ h(t), where |h(t)| << 1. Subtitute x(t) into (13) and expand f (x) by Taylor series: ˙x∗+ ˙h = f (x∗+ h) = f (x∗) + Df (x∗)h + O(|h2|). The notation Df (x∗) is the n × n
Jacobian matrix of partial derivative of a vector-valued function f .
The eigenvalues and eigenvectors of the matrix Df (x∗) determine the general solution. In studying stability we want to know whether the solution grows, stays constant, or decay to 0 as t → ∞. It can be answered by evaluating the eigenvalues.
If λ is a real eigenvalue with eigenvector v, then there is a solution to the linearization equation of the form: h(t) = cveλt. If λ = a ± ib is a complex conjugate pair with eigenvectors v = u ± iw (where u, w are real), then h1(t) = eat(u cos bt − w sin bt) and
h2(t) = eat(u cos bt + w sin bt) are two linearly independent solutions. In both cases, the
real part of λ almost determines stability. Any solution of the linearized equation can be written as a linear combination of terms of these forms. We can obtain the following conclusions:
1. If all eigenvalues of Df (x∗) have negative real parts, then |h(t)| → 0 as t → ∞ for all solutions.
2. If there exists one eigenvalue of Df (x∗) has a positive real part, then there is a solution h(t) with |h(t)| → +∞ as t → ∞.
3. If some pair of complex-conjugate eigenvalues have zero real parts with distinct imaginary parts, then the corresponding solutions for |h(t)| oscillate as t → ∞ and neither decay nor grow as t → ∞.
Moreover, if x∗ is an equilibrium of ˙x = f (x) and all the eigenvalues of the matrix Df (x∗) have strictly negative real parts, then x∗is asymptotically stable. If at least one eigenvalue has strictly positive real part, then x∗ is unstable.
We had discussed the solution of the ODE is an equilibrium. The statement of stability may be extend to non-constant orbits of ODE.
Here, let Ot(x) = x(t) and the initial value x(0) = x0. Then the set O(x) = {Ot(x) :
0 ≤ t} is called orbit.
Definition 3. Let two orbits O(x) and O(ˆx). If there is a reparameterization of time ˆt(t) such that |Ot(ˆx) − Oˆt(t)(ˆx)| < for all t ≥ 0,then we say two orbits O(x) and O(ˆx) are
-close.
Definition 4. An orbit O(x) is orbitally stable if for any > 0, there is a neighborhood V of x so that, for all ˆx ∈ V , O(x) and O(ˆx) are -close. If additionally V may be chosen so that, for all ˆx ∈ V , there exists a constant τ (ˆx) so that |Ot(ˆx) − Ot−τ (ˆx)(ˆx)| < as
t → ∞, then Ot(x) is asymptotically stable.
The linearization skill in general orbit is similar to the the pervious linearization method. Also discuss the eigenvalues of the linearization operator.
3
Numerical Method
In this section, we discretized the Laplace operator by finite-difference method first, and we will show it on two kinds of domain respectively. Then we present the numerical method for solving φ(x) from (5) first. And finally, we will describe a numerical method to compute the spectra of L from (7).
Now, we recall our main question: to study the stability of the special solution form. For our solution form ψ(t, x) = eiλtφ(x), we will solving the time-inedpent term φ(x).
By (5), φ = φ(x1, x2) satisfies − ∂ 2 ∂x2 1 φ + ∂ 2 ∂x2 2 φ + (V (x1, x2) + λ)φ − m(x1, x2)φ3 = 0. (14)
The natural boundary condition is lim
|x|→∞φ(x1, x2) = 0. (15)
Consider the same equation (5) in the unit disk domain Ω = {(x1, x2) : x21+ x22 < 1},
applying the polar coordinate transformation,
x1 = r cos θ, x2 = r sin θ,
where r = px2
1+ x22, θ = tan −1(x
2/x1). We can rewrite (14), φ = φ(r, θ) in the polar
coordinate system satisfies − 1 r ∂ ∂r(r ∂φ ∂r) + 1 r2 ∂2φ ∂θ2 + (V (r, θ) + λ)φ − m(r, θ)φ3 = 0, (16) with 0 < r < 1, 0 ≤ θ < 2π. And the boundary condition is
lim
r→∞φ(r, θ) = 0. (17)
3.1
Matrix Form of Laplace Operator
Now, we are going to describe the matrix form of the Laplace operator. In this study, we have two kinds of domains, unit disk domain and square domain. We will show them respectively.
3.1.1 Square Domain
The Laplace operator is a seconed order differential operator. Here, we consider the function φ in R2. Then ∆φ = ∂ 2φ ∂x2 + ∂2φ ∂y2.
Here we use Finite-Difference method for the boundary value problem [1].
First step is to partition x into −1 = x0 < x1 < x2 < · · · < xN −1 < xN = 1, and
partition y into −1 = y0 < y1 < y2 < · · · < yN −1 < yN = 1, where xi = −1 + i∆x,
yj = −1 + j∆y and h = ∆x = ∆y = 2/N . We use the Taylor series in the variable x, y
about xi, yj to generate the centered-difference formula
∂2φ ∂x2 =
fi−1,j − 2fi,j + fi+1,j
(∆x)2 + O((∆x) 2),
∂2φ
∂y2 =
fi−1,j − 2fi,j + fi+1,j
(∆y)2 + O((∆y) 2
)
for each i = 1, 2, ..., N − 1 and j = 1, 2, ..., N − 1, where fi,j = φ(xi, yj).
In difference-equation form, this result in the finite-difference method, with local trun-cation error of order O((∆x)2+ (∆y)2). Therefore
−∆φ ' − fi−1,j − 2fi,j + fi+1,j (∆x)2 +
fi,j−1− 2fi,j+ fi,j+1
(∆y)2
= −fi−1,j + fi+1,j− 4fi,j+ fi,j−1+ fi,j+1
h2 .
We can rewrite it as a matrix representation:
1 h2 ˆ A −I 0 −I Aˆ −I . .. ... ... −I Aˆ −I 0 −I Aˆ f1 f2 .. . fN−2 fN− 1 ≡ AF, where ˆ A = 4 −1 0 −1 4 −1 . .. ... ... −1 4 −1 0 −1 4 (N −1)×(N −1) , fj = [f1,j, f2,j, ..., fN −1,j]T,j = 1, 2, ..., N − 1, and I is an (N − 1)2 × (N − 1)2 identity matrix.
3.1.2 Unit Disk Domain
Then we consider the equation (5) over a 2-dimensional disk. We use a discretization scheme [8] for equation (16). The Laplace operator is
∆φ(r, θ) = ∆rφ + ∆θφ = 1 r ∂ ∂r r∂φ ∂r + 1 r2 ∂2φ ∂θ2.
We discretize the Laplace operator in two parts. The first step is to discrstize the operator ∆rφ on r ∈ [0, 1], and then to discretize the operator ∆θφ on θ ∈ [0, 2π].
Firstly, we consider the operator −∆rφ on r ∈ [0, 1],
−∆rφ = − 1 r∂r r∂φ ∂r . To partition r into 0 = r0 < r1 2 < r1 < · · · < rN − 1 2 < rN = 1, rj = j∆r and ∆r = 1 N, and
to partition θ into 0 = θ0 < θ1 < θ2 < · · · < θM −1 < θM = 2π, θj = j∆θ and ∆θ = 2πM,
and denote fi,j = φ(ri, θj). For each θj,
−∆rφ '
" − 1
ri−1 2
ri(fi+12,j− fi−12,j) − ri−1(fi−12,j − fi−32,j)
(∆r)2 # . That is 1 (∆r)2 r1 r1 2 −r1 r1 2 0 −r1 r3 2 1 r3 2 (r1 + r2) −rr23 2 . .. . .. . .. −rN −2 rN − 3 2 1 rN − 3 2 (rN −2+ rN −1) −rrN −1 N − 32 0 −rN −1 rN − 1 2 1 rN − 1 2 (rN −1+ rN) f1 2,j f3 2,j .. . fN −3 2,j fN −1 2,j ≡ bAf , where b A = 1 (∆r)2 a1 c1 0 b1 a2 c2 . .. ... . .. bN −2 aN −1 cN −1 0 bN −1 aN and fj = f1 2,j f3 2,j .. . fN −1 2,j , ai = 1 ri−1 2 (ri−1+ ri) = i − 1 + i i −12 = 2, i = 1, . . . , N and bi = − ri ri+1 2 = −1 + 1 2i + 1, ci = − ri ri−1 2 = −1 − 1 2i − 1, i = 1, . . . , N − 1.
Now we consider the operator −∆θφ = − 1 r2 ∂2φ ∂θ2
on θ ∈ [0, 2π]. Then for each ri−1
2, i = 1, 2, ..., N −∆θφ ' " − 1 r2 i−12 f (ri−1 2, θj+1) − 2f (ri− 1 2, θj) + f (ri− 1 2, θj−1) (∆θ)2 # . Let B = diag 2 r2 1 2 (∆θ)2, 2 r2 3 2 (∆θ)2, 2 r2 5 2 (∆θ)2, . . . , 2 r2 N −12(∆θ) 2 ! and C = diag −1 r2 1 2 (∆θ)2, −1 r2 3 2 (∆θ)2, −1 r2 5 2 (∆θ)2, . . . , −1 r2 N −1 2 (∆θ)2 ! , fj = [f1 2,j, f 3 2,j, . . . , fN − 1 2,j] > for j = 1, 2, ..., M. Then we can rewrite the Laplace operator as
b A + B C C C A + Bb C . .. . .. . .. C A + Bb C C C A + Bb f1 f2 .. . fM−1 fM ≡ AF.
Then the matrix A represents the Laplace operator.
3.2
Algorithm
Now, we describe a numerical method to compute the spectrum of the linear operator L defined by (7) for > 0. There are two steps in our numerical method: First, we will solve the nonlinear problem (14) for φ. Second, compute the spectrum of the discretized linearized operator around φ.
Step I. Compute φ and it is energy minimizing among all solution of (14)–(17).
Step II. Compute the spectra of the discretized linearized operator L. Since h is no longer real, so the linearized operator is a little different from L in the equation (7). Here, we denote some notation of the following. For A ∈ RM ×M, q, m ∈ RM, m ∗ q
denote the Hadamard product of m and q, and q p means the p-time Hadamard product
of q. And diag(q) represents the diagonal matrix of q. In Step I. First, we represent −∆q as
where q is an approximation of φ. As for the form of matrix A, it was shown in previous two subsections. Then the discretization of the equation (14) (and (16)) become
Aq + (V + λ) ∗ q − m ∗ q 3
= 0,
where V and m are approximation of the function V (x1, x2) and m(x1, x2), respectively.
We solve it by using an iteration method [7]:
[A + diag(V + λ)]˜qnew = m ∗ q 3
old, (18)
where ˜qnew and qold are unknown and known vector. The iteration step is shown in
Algorithm 1 .
Algorithm 1 Iterative algorithm for solving φ(x)
Step 1 Choose an initial guess of qold> 0, and ||qold||2 = 1.
Step 2 Solve (18), then obtain ˜qnew.
Step 3 Let αnew = ||˜qnew1 ||
2, qnew = αnewq˜new.
Step 4
if (converge) then
Output the solution (αnew)
1
2qnew. Stop.
else
Let qold= qnew.
Go to Step 2. end if
When λ in a general case, λ is no longer a constant. λ varies from , so the numerical method has a little difference form. From (4), the discretization form becomes
[A + diag(V − m ∗ q 2
)]q = −λq. (19)
It turns into an eigenvalue problem. In Algorithm 2, it shows the iteration step. In Step II. Now we discretize L of (8) into an eigenvalue problem.
L u w = ρ u w , (20) where L = 0 A + diag(V + λ) − diag(m ∗ q2)
−A − diag(V + λ) + diag(3m ∗ q2) 0
. q and λ are obtain from Step I. We use ARPACK in MATLAB version R2007a to solve the linear algebraic eigenvalue problem and obtain eigenvalues ρ of L near origin.
4
Numerical Simulation
For each potential case, we summary the numerical results for three solution forms. Con-sider the solution form: ψ(t, x) = eiλtφ(x), λ = 0, 1, and general case.
Algorithm 2 Iterative algorithm for solving φ(x) in general solution Step 1 Choose an initial guess of qold> 0, and ||qold||2 = 1.
Step 2 Solve
(A + diag(V − m ∗ q 2
old))˜qnew = −λ˜qnew,
where −λ is the smallest eigenvalue of A + diag(V − m ∗ q 2
old). Then obtain ˜qnew.
Step 3 Let αnew = |−λ|
||˜qnew||2, qnew = αnewq˜new.
Step 4
if (converge) then
Output the solution qnew and λ. Stop.
else
Let qold= qnew.
Go to Step 2. end if
4.1
Setting
In ω case, the discretized matrix of Laplace operator has size N T by N T with N = r/∆r and T = 2π/∆θ. Here we use zero boundary condition and r = 1, N = 32, T = 64. And in µ case with square domain, the discretized natrix of Laplace operator has size N2 by
N2 with N = 2/h and here N = 64. And so the matrix size of the operator L are 2N T by 2N T and 2N2 by 2N2, respectively.
In our numerical method, we use finite-difference method to solve NLS equation. The finite-difference method has truncation error O(h2), where h is the grid size. In our experiment, O(h2) ≈ c · 10−3.
The region of that we consider about is larger than 10−5, that is because the value is smaller than 10−5 we treat it as zero. We should control the region of to satisfy the boundary condition, which means the region of would not be large. In the region that we considered, the boundary mean values of all cases are less than 10−5. In Fig. 4, we plot the mean value of the boundary of φ for from 0 to 0.003. We can see the mean of boundary will be less than 10−5 as smaller than 0.0024. So in that case, we only try ∈ [10−5, 0.0024] in our numerical experiment.
Besides the boundary condition, we also take care of the shape of φ, here we narrow down the focus on the solitary solution for this part only. As goes larger, sometimes the shape would change. The following statement like ∗ are all in the region of .
In the following table(Table 1), we write down the region of for each case. In each potential case, for our converient, we choose the intersection of region for different parameters.
In this project, we focus on the parameter changing, and to study when the solution of the NLS would be unstable or stable. The notation ∗ denote as the bifurcation of stable and unstable.
We find that a pair of purely imaginary eigenvalues will collide at the origin and split into a pair of real eigenvalues for < ∗. That is, when less than ∗, the spectrum of the linearized operator are all pure imaginary. As larger than ∗, the spectrum of the linearized operator has at least one pair of eigenvalue with nonzero real part. That means
Table 1: The region of . case µ case ω λ = 0 [10−5, 0.016] [10−5, 0.0024] λ = 1 [10−5, 0.025] [10−5, 0.005] general λ [10−5, 0.001] [10−5, 0.001]
Figure 8: Solution of φ for ω1 = 1, ω2 = 1, = 0.001, λ = 1.
the solution is stable while < ∗, unstable when > ∗.
Since the numerical results of eigenvalues are approximartion of the exact mathemat-ical values. There exist some error of the result. We determined that the eigenvalue has nonzero real part when the absolute value of the real part of eigenvalus is larger than 10−3.
Since the regions of are small, and the spectrum of lineraization operation changes so fast. We move in a small step for each case. δ = 10−5. To prevent losing some pairs of pure complex eigenvalue from turning to real eigenvalues, we also move the target of the algorithm for finding eigenvalues.
This method is not an efficient method but an easy way to catch the changing of spectra.
4.2
ω case
In ω case, the maxima of φ is not at center of the disk; instead, it is on the right of the center. As go larger and the potential change, the maxima of φ seem to moving to the center of domain. In figure 10, we fix = 0.0002 and changing the well depth of potential V (x). We find that the location of maxima φ changes. It moves to the center of the domain.
We know trap potential V (x) and m(x) would effect the solution φ. The NLS equation in our model is an focusing case, that is the atoms would stay near the maxima of m(x), ie, x1 = 1/4.
1. λ = 0: Testing in the region that we mentioned before, and find that the eigen-values of the linearized operation L are all pure imaginary.
2. λ = 1: We find there is an ∗ in this solution case. The ∗ represent the bifurcation of pure imaginary to has a pair of eigenvalue with nonzero real part. In Fig. 11 shows the spectrum of L . And in Fig. 12, it shows the ∗ of the parameters ω1, ω2 = 1, 2, ..., 5
Figure 9: Solution of φ for ω1 = 1, ω2 = 1, = 0.001, λ = 0.
3. In elliptic potential case, we can find ∗in eitφ(X) this solution form in our numerical
method, but we can not find the ∗ in other solution form. At first we guess ∗ in ei0tφ(x) form of solution can be found in smaller . But we do not find that in the region of as we mentioned before.
4. λ is general case: In this solution case, we can not find the ∗. When is quite small, the value of λ is already negative. That means the trap potential is not positive at all, V (x) < 0 on the center of the disk. In Figure 13 it shows the value of λ as changes.
Figure 14: Spectrum of µ case for λ = 0.
4.3
µ case
1. λ = 0, we can find ∗ in this solution case. For µ1, µ2 = 1, 2, ..., 5, we find that the
eigenvalues are pure imaginary when < ∗, and turn into a pair real value. Figure 14 plot the eigenvalues around original. And in Figure 15 show the ∗ for every µ1, µ2.
2. λ = 1, we also find ∗. Figure 16, 17 show the eigenvalues and values of ∗ for each µ1, µ2.
3. Consider the same µ1 and µ2, these two solution forms can be seen as changing
the well depth of potential (from (5)). The ∗ are different. So we know that the potential well depth may change the stability of solution.
In optical lattices potential case, changing µ1or µ2for different periodic of potential.
We can see that the value of ∗ increasing when µ1, µ2 go larger. To comparison
eitφ(x) and ei0tφ(x) these two forms of solutions, we can see that the second kind
of ∗ is less than the other ones.
4. λ is general case, there is no ∗ in this solution case. When is quit small( less than 1.2 × 10−4), the value of λ is positive(In Figure 18). This result is consistent with the condition of Lin in [9](λ is positive for small enough ). But λ become negative when gets larger, and the value of λ gets smaller as gets larger. When < 3.3 × 10−4, λ is smaller than −1. And the minimum of trap potential is 1, that means the trap potential is no longer positive value.
5
Conclusions
In this project, we consider the NLS equation with the special solution eiλtφ(x). We
mainly show the spectrum of the linerization operator L. By numerical computation, we can conclude some results: when < ∗ the eigenvalues are pure imaginary, and the eigenvalues turn to real as > ∗. That is, there is a bifurcation of stable and unstable, when λ = 0, 1 in optical lattices potential case and λ = 1 in elliptic potential case.
Here, we use Matlab code to find the spectrum of the linearization operator. It is a convenient way to solve the eigenvalue problem. Since the method of finding eigenvalue in Matlab is to search some eigenvalues which are near the target. And the distribution of eigenvalues changes fast when near ∗.
References
[1] R. Burden and J. D. Faires. Numerical Analysis. Brooks Cole, 7 edition, 2000. [2] N. Cesa-Bianchi and G. Lugosi. Prediction, Learning, and Games. Cambridge
Uni-versity Press, 2006.
[3] S. M. Chang, S. Gustafson, K. Nakanishi, and T. P. Tsai. Spectra of linearized oper-ators for NLS solitary waves. SIAM Journal on Mathematical Analysis, 39(4):1070– 1111, 2007.
[4] M. W. Hirsch, S. Smale, and R. L. Devaney. Differential equations, dynamical sys-tems, and an introduction to chaos. Academic Press, 2004.
[5] P. Holmes and E. T. Shea-Brown. Stability. Scholarpedia, 1(10):1838, 2006.
[6] T. M. Hwang and W. Wang. Analyzing and visualizing a discretized semilinear elliptic problem with neumann boundary conditions. Numerical Methods for Partial Differential Equations, 18(3):261–279, 2002.
[7] M. C. Lai. A note on finite difference discretizations for poisson equation on a disk. Numerical Methods for Partial Differential Equations, 17:199–203, 2001.
[8] T. C. Lin and J. Wei. Orbital stability of bound states of semi-classical nonlinear schr¨odinger equations with critical nonlinearity. SIAM Journal on Mathematical Analysis, 40(1):365–381, 2008.
[9] T. C. Lin, Juncheng Wei, and Wei Yao. Orbital stability of bound states of nonlin-ear schr¨odinger equations with linear and nonlinear lattices. Journal of Differential Equations, page to appear, 2010.
[10] M. P. MacDonald, G. C. Spalding, and K. Dholakia. Microfluidic sorting in an optical lattice. Nature, 426:421–424, 2003.
[11] C. Sulem and P.-L. Sulem. The nonlinear Schrodinger equation: Self-Focusing and Wave Collapse. Springer, 1999.
出席國際學術會議心得報告
報告人姓名 張書銘 服務機構及職稱 國立交通大學 助理教授 會議時間 地點 2010.06.27 – 2010.07.01 德國、柏林 本會核定 補助文號 NSC 98-2115-M-009-007 會議 名稱 (中文) 第八屆特徵值問題精確解國際會議(英文) 8th International Workshop on Accurate Solution of Eigenvalue Problems (IWASEP 8) 報告內容包括下列各項： 一、參加會議經過 本人於 2010 年 2 月開始規劃要參加 2010 年 6 月底的 IWASEP 8-第八屆特徵值 問題精確解國際會議，此國際會議是每兩年固定時間主辦的一個特徵值問題計算的 重要聚會，此會議雖然僅為期四天，但是這四天的學術演講都是單場的，沒有平行 的演講場次。因此，每個參與者都是同領域的學者，演講者更是在此方面的研究的 翹楚者。 這次的參與是本人第一次以參與者身份出席國際會議，雖然沒有在外語演講的 經驗上有所累積，但在與外國人的參與者之互動上，本人有更多的機會與其對談， 並且進一步認識了荷蘭的學者 Michiel Hochstenbach 教授(Dept. of Math. And Computer Science, TU Eindhoven, Netherlands)。期望今年或是明年，能夠邀請 Prof. Hochstenbach 來台灣訪問。 二、與會心得 首先要感謝行政院國家科學委員會給予本人出席國際會議的差旅補助，讓本人 有機會去見識國際會議，參與過程本人更增加不少國際觀與吸收到專業的研究知 識。 當中本人與 Prof. Hochstenbach 有相當多次的交談，其研究領域相當廣泛 (Numerical linear algebra, Eigenvalue problems, Model reduction, Ill-posed problems, inverse problems, regularization, Image analysis, Matrix functions, Linear systems, Scientific Computing, Nonlinear systems, Boundary element method, Numerical crack propagation, Industrial
mathematics via LIME)。期望今年或是明年，能夠邀請 Prof. Hochstenbach 來台 灣訪問。
三、考察參觀活動 無。
四、建議 無。 五、攜回資料名稱及內容 大會演講論文全文一冊及會議演講摘要手冊（內容包括此次會議的所有演講場 次與講題），收集完整的大會演講內容，相當豐富。 六、其他 無。
98 年度專題研究計畫研究成果彙整表
計畫主持人：張書銘 計畫編號： 98-2115-M-009-007-計畫名稱：格若士-比塔烏斯基方程式的行進波解之穩定性研究 量化 成果項目 實際已達成 數（被接受 或已發表） 預期總達成 數(含實際已 達成數) 本計畫實 際貢獻百 分比 單位 備 註 （ 質 化 說 明：如 數 個 計 畫 共 同 成 果、成 果 列 為 該 期 刊 之 封 面 故 事 ... 等） 期刊論文 0 0 0% 研究報告/技術報告 0 0 0% 研討會論文 0 0 0% 篇 論文著作 專書 0 0 0% 申請中件數 0 0 0% 專利 已獲得件數 0 0 0% 件 件數 0 0 0% 件 技術移轉 權利金 0 0 0% 千元 碩士生 3 0 100% 研 究 生 兼 任 助 理 參 與 本 計 畫 , 在 程 式 設 計 能 力 上 均有成長. 博士生 0 0 0% 博士後研究員 0 0 0% 國內 參與計畫人力 （本國籍） 專任助理 0 0 0% 人次 期刊論文 0 0 0% 研究報告/技術報告 0 0 0% 研討會論文 0 0 0% 篇 論文著作 專書 0 0 0% 章/本 申請中件數 0 0 0% 專利 已獲得件數 0 0 0% 件 件數 0 0 0% 件 技術移轉 權利金 0 0 0% 千元 碩士生 0 0 0% 博士生 0 0 0% 博士後研究員 0 0 0% 國外 參與計畫人力 （外國籍） 專任助理 0 0 0% 人次其他成果
