國立交通大學 應用數學系 碩士論文 橢圓與光學晶格位能井之非線性薛丁格方程 的數值研究 Numerical Studies in Nonlinear Schrödinger Equations with Elliptic and Optical Lattices of Trap Potentials

Transcription

1 國立交通大學 應用數學系 碩士論文 橢圓與光學晶格位能井之非線性薛丁格方程 的數值研究 Numerical Studies in Nonlinear Schrödinger Equations with Elliptic and Optical Lattices of Trap Potentials 研究生 : 莊佩錚 指導教授 : 張書銘博士 中華民國九十八年七月

2 橢圓與光學晶格位能井之非線性薛丁格方程的數值研究 Numerical Studies in Nonlinear Schrödinger Equations with Elliptic and Optical Lattices of Trap Potentials 研究生 : 莊佩錚 指導教授 : 張書銘博士 Student: Pei-Cheng Chuang Advisor: Dr. Shu-Ming Chang 國立交通大學應用數學系碩士論文 A thesis Submitted to Department of Applied Mathematics college of Science National Chiao Tung University in Partial Fulfillment of the Requirements for the Degree of Master in Applied Mathematices July 009 Hsinchu, Taiwan, Republic of China 中華民國九十八年七月

3 橢圓與光學晶格位能井之非線性薛丁格方程的數值研究 學生 : 莊佩錚 指導教授 : 張書銘博士 國立交通大學應用數學系 ( 研究所 ) 碩士班 摘 要 本論文主要運用數值方法來探討非線性薛丁格方程的解之穩定性分析, 其中這邊考慮的非線性薛丁格方程式是含有位能井 (trap potential) 首先, 我們先介紹非線性薛丁格方程式以及兩種不同的位能井 接著, 簡述我們所運用的數值方法 最後, 我們呈現一些數值模擬結果 經由我們的數值模擬, 發現在一些解的形式與位能井的選擇是可以觀察到解的穩定性變化 i

4 Numerical Studies in Nonlinear Schrödinger Equations with Elliptic and Optical Lattices of Trap Potentials Student: Pei-Cheng Chuang Advisors: Dr. Shu-Ming Chang Department (Institute) of Applied Mathematics National Chiao Tung University Abstract In this thesis, we would like to find out the stability of the solution which linearized around the special solution of the nonlinear Schrödinger equation. Here, we consider nonlinear Schrödinger equations with the potential part, optical lattices form and elliptic form. In the first part of the thesis, we introduce roughly the nonlinear Schrödinger equation and two kinds of trap potentials. Then, we show the mathematical analysis of the numerical method that we used. At the end, we give some numerical results of our numerical experience. ii

5 誌 謝 本篇論文的完成, 最要感謝的是我的指導老師張書銘博士, 於這兩年來的用心指導, 讓我得以對數值計算有更進一步的了解 除了研究上的指導, 亦在未來人生的規劃與為人處事上給予我許多關懷指引 老師, 謝謝您! 於口試期間, 亦承蒙黃聰明教授與吳金典博士給予我的建議及對疏漏處的導正, 使本論文得以更加嚴謹完整, 在此學生致上由衷的感謝 另遠在英國的芳嫻姐姐及我的室友們 : 郁凌 廷芳 維彧, 謝謝你們在我完成論文的階段中, 給予我的幫助 鼓勵與歡笑, 有了你們的陪伴, 讓我的研究生活更加順利 最後, 謝謝家人對我的支持與關心, 有了您們的鼓勵與幫助, 讓我得以專心完成本篇論文 僅將此小小成果與所有關心 幫助過我的長輩 老師 同學們一同分享這份喜悅, 謝謝大家! iii

6 Contents 中文摘要英文摘要誌謝 Contents List of Tables List of Figures i ii iii iv v vi 1 Introduction Nonlinear Schrödinger Equations Trap Potential V Optical Lattices Elliptic Potential Goals Mathematical Analysis 6.1 Perturbation in NLS Stability Numerical Method Matrix Form of Laplace operator Square Domain Unit disk Domain Algorithm Numerical Simulation Setting ω case µ case Conclusions 30 Reference 31 iv

7 List of Tables 1 The region of ɛ v

8 List of Figures 1 Optical Lattices Elliptic potential and m(x) = sin(πx 1) Optical lattices potential and m(x) = cos(πx 1) Boundary φ in ω case for λ = Boundary φ in ω case for λ = Boundary φ in µ case for λ = Boundary φ in µ case for λ = Solution of φ for ω 1 = 1, ω = 1, ɛ = 0.001, λ = Solution of φ for ω 1 = 1, ω = 1, ɛ = 0.001, λ = φ moves for λ = 1, ɛ = Spectrum of ω case for λ = ɛ of ω case for λ = Value of λ in ω case Spectrum of µ case for λ = ɛ of µ case for λ = Spectrum of µ case for λ = ɛ of µ case for λ = Value of λ in µ csae vi

9 1 Introduction In this thesis, we consider the nonlinear Schrödinger equation (NLS) with focusing parameter ɛ, i t ψ = ɛ ψ + V (x)ψ m(x) ψ ψ, where ψ(t, x) : R R 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. Firstly, we introduce nonlinear Schrödinger equations and the two kinds of trap potentials first. 1.1 Nonlinear Schrödinger Equations The Schrödinger equation is a basic equation in quantum theory. The study of Schrödinger equations plays an important role in modern physics. In 196, an Austrian physicist Erwin Schrödinger constructed the Schrödinger equation for explaining the active of particles. The famous Schrödinger equation is the form: Hψ = i ψ, where H is a Hamiltonian t operator, ψ is the wavefunction, represents Planck s constant over π, i is the imaginary unit. The Schrödinger equation is an equation to describe the possible distribution of atoms. Solutions to Schrödinger s equation describe not only atomic and subatomic systems, electrons and atoms, but also macroscopic systems, possibly even the whole universe. The nonlinear Schrödinger equation (NLS) is a nonlinear version of Schrödinger s equation in theoretical physics. In mathematical point of view, the NLS equation is a Schrödinger equation with nonlinear term. The nonlinear Schrödinger equation is the partial differential equation i t ψ = 1 xψ + c ψ ψ for the complex field ψ. There are many kinds of nonlinear Schrödinger equations. One general form of such equations [1] 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 ψ + ψ = ± ψ ψ. The sign on the right hand side is focusing nonlinearity, and the sign + is defocusing. 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, this can be used to model a confining magnetic trap for Bose-Einstein condensation. 1

10 1. 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 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(figure reference: [6]). 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 particles [11] recently. 1.. Elliptic Potential Let d d matrix M be symmetric, positive define. Let v R d and c R be arbitrary. The elliptic potential [] defines as V (u) = 1 ut Mu + u T v + c, where u is a column vector in R d. Because M is symmetric and psoitive definite, M 1 exists. Therefore, we can also write V (u) = 1 M 1 u + u T v + c. In R, the elliptic potential can be rewritten as V (x, y) = ax + bxy + cy + dx + ey + f, a, b, c, d, e, f R with b 4ac < 0. The shape of this potential looks like a bowl.

11 1.3 Goals In this thesis, we consider nonlinear Schrödinger equations with focusing parameter ɛ, i t ψ = ɛ ψ + V (x)ψ m(x) ψ ψ, () where ψ(t, x) : R R 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 thesis, we consider two kinds of different trap potentials and domains as follows: Case 1. (ω case) V (x) = (ω 1 x 1 ) + (ω x ), where ω 1, ω 0 and ω 1, ω = 1,,..., 5, and m(x) = sin(πx 1) with unit disk domain Ω = {(x 1, x ) x 1 + x < 1}. Case. (µ case) V (x) = 1 + sin (µ 1 x 1 ) + sin (µ x ), where µ 1, µ = 1,,..., 5 and m(x) = cos(πx 1) with square domain Ω = {(x 1, x ) 1 < x 1 < 1, 1 < x < 1}. Here, we have two kinds of trap potentials and function m(x) in Fig. 3 and Fig. 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 ω for the elliptic trap potential. And then we control the periodic for the optical lattices trap potential by changing the two parameters, µ 1 and µ. 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) = e iλt φ(x). The aim of this thesis is to find out the spectra of linearized operator which arises when the equation () 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 stability and instability of single-spike bound state can be obtained in [9, 10]. Next, in Section 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. 3

12 Figure : Elliptic potential and m(x) = sin(πx 1). 4

13 Figure 3: Optical lattices potential and m(x) = cos(πx 1). 5

14 Mathematical Analysis.1 Perturbation in NLS Recall the nonlinear Schrödinger equation: i t ψ = ɛ ψ + V (x)ψ m(x) ψ ψ. (3) We consider the special solutions of NLS equation (3): ψ(t, x) = e iλt φ(x), where λ = 0, 1, and general case. φ(x) is a real-valued function independent of time. Substitution ψ(t, x) = e iλt φ(x) into (3) and the φ satisfy the nonlinear elliptic equation: that is or ɛ φ + V φ m φ φ = λφ, ( ɛ + V m φ )φ = λφ (4) ( ɛ + (V + λ) m φ )φ = 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) = e iλt (φ(x) + h(t, x)). (6) The perturbation h(t, x) R satisfies an equation: t h = Lh + (nonlinear terms), where Lh = i{( ɛ + (V + λ) 3mφ )h}. (7) Remark. To explain how to obtain (7): In equation (3), replace ψ(x) by e iλt (ψ(x) + h(t, x)) then the left hand side of (3): And the right hand side of (3): i t ψ = i t (e iλt (φ + h)) = i(iλe iλt (φ + h) + e iλt t h) = e iλt ( λ(φ + h) + i t h). ɛ ψ + V (x)ψ m(x) ψ ψ = e iλt ( ɛ φ ɛ h + V (x)ψ + hφ m φ + h (φ + h)) = e iλt { ɛ φ + V φ ɛ h + hφ m( φ φ + φ ( h + h) + φh h + φ h + φ(h + h)h + h h)} = e iλt ( ɛ φ + V φ m φ φ) +e iλt ( ɛ h + V h m φ (h + h) m φ h) + e iλt O(h ). Notice that φ + h = (φ + h)(φ + h) = φ + φ( h + h) + h h. Then e iλt ( λ(φ + h) + i t h) = e iλt ( ɛ φ + V φ m φ φ) +e iλt ( ɛ h + V h m φ (h + h) m φ h) + O(h ). 6

15 So i t h = ( ɛ φ + (V + λ)φ m φ φ) +( ɛ h + (V + λ)h m φ (h + h) m φ h) + O(h ) = ɛ h + (V + λ)h m φ (h + h) m φ h + O(h ). 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φ ) Re h + i( ɛ + (V + λ) mφ ) Im h mφ Re h} = ( ɛ + (V + λ) mφ ) Im h + i{ ( ɛ + (V + λ) 3mφ ) Re h}. And we rewrite L as a matrix acting on where L + and L are self-adjoint. [ L = [ Re h Im h ], 0 L L + 0 ], (8) L + = ɛ + (V + λ) 3mφ, L = ɛ + (V + λ) mφ. (9) 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? [ ] 0 L Lemma 1. Let L = and Lh = ɛ h + (V + λ)h m φ L + 0 (h + h) m φ h, where L + and L are defined by (9). If ρ is an eigenvalue of L and [u T, v T ] 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, [ ] [ u u L = ρ w w So [ Then and (1)+i (11), therefore 0 L L + 0 ] [ u w ] ]. (10) [ u = ρ w ( ɛ + (V + λ) 3mφ )u = ρw (11) ]. ( ɛ + (V + λ) mφ )w = ρu. (1) ρ(u + iw) = i( ɛ + (V + λ) 3mq )u + ( ɛ + (V + λ) mq )w = i{( ɛ + (V + λ) 3mq )u + i( ɛ + (V + λ) mq )w} = i{( ɛ (u + iw) + (V + λ)(u + iw) mq (u + iw) mq u}. Hence u + iw is an eigenvactor of L corresponding to ρ. 7

16 . Stability Our aim of this thesis 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 = ẋ = f(x); x Rn, (13) where f : R n R n and x = [x 1, x,..., x n ] 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) = x 0 V is defined and remains in U for all t 0. Definition. If V can be chosen above so that, in addition to the properties for stability, we have lim t 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: ẋ + ḣ = f(x + h) = f(x ) + Df(x )h + O( h ). 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 h 1 (t) = e at (u cos bt w sin bt) and h (t) = e at (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.. 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 ẋ = 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. 8

17 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 O t (x) = x(t) and the initial value x(0) = x 0. Then the set O(x) = {O t (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 O t (ˆ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 O t (ˆx) O t τ(ˆx) (ˆx) < ɛ as t, then O t (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) = e iλt φ(x), we will solving the time-inedpent term φ(x). By (5), φ = φ(x 1, x ) satisfies [ ɛ φ + φ x 1 x The natural boundary condition is ] + (V (x 1, x ) + λ)φ m(x 1, x )φ 3 = 0. (14) lim φ(x 1, x ) = 0. (15) x Consider the same equation (5) in the unit disk domain Ω = {(x 1, x ) : x 1 + x < 1}, applying the polar coordinate transformation, x 1 = r cos θ, x = r sin θ, where r = x 1 + x, θ = tan 1 (x /x 1 ). We can rewrite (14), φ = φ(r, θ) in the polar coordinate system satisfies [ 1 φ ɛ (r r r r ) + 1 ] φ + (V (r, θ) + λ)φ m(r, θ)φ 3 = 0, (16) r θ with 0 < r < 1, 0 θ < π. And the boundary condition is lim φ(r, θ) = 0. (17) r Here, we have two equations and boundary conditions for solving φ in each domain. 9

18 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 Square Domain The Laplace operator is a seconed order differential operator. function φ in R. Then φ = φ x + φ y. Here, we consider the Here we use Finite-Difference method for the boundary value problem [1]. First step is to partition x into 1 = x 0 < x 1 < x < < x N 1 < x N = 1, and partition y into 1 = y 0 < y 1 < y < < y N 1 < y N = 1, where x i = 1 + i x, y j = 1 + j y and h = x = y = /N. We use the Taylor series in the variable x, y about x i, y j to generate the centered-difference formula φ x = f i 1,j f i,j + f i+1,j ( x) + O(( x) ), φ y = f i 1,j f i,j + f i+1,j ( y) + O(( y) ) for each i = 1,,..., N 1 and j = 1,,..., N 1, where f i,j = φ(x i, y j ). In difference-equation form, this result in the finite-difference method, with local truncation error of order O(( x) + ( y) ). Therefore [ fi 1,j f i,j + f i+1,j φ + f ] i,j 1 f i,j + f i,j+1 ( x) ( y) = f i 1,j + f i+1,j 4f i,j + f i,j 1 + f i,j+1 h. We can rewrite it as a matrix representation: where 1 h  I 0 I  I  = I  I 0 I  f 1 f. f N f N 1 AF, (N 1) (N 1) f j = [f 1,j, f,j,..., f N 1,j ] T,j = 1,,..., N 1, and I is an (N 1) (N 1) identity matrix., 10

19 3.1. Unit disk Domain Then we consider the equation (5) over a -dimensional disk. We use a discretization scheme [8] for equation (16). The Laplace operator is φ(r, θ) = r φ + θ φ = 1 ( r φ ) + 1 φ r r r r θ. 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, π]. Firstly, we consider the operator r φ on r [0, 1], r φ = 1 ( r r r φ ). r To partition r into 0 = r 0 < r 1 < r 1 < < r N 1 < r N = 1, r j = j r and r = 1 N, and to partition θ into 0 = θ 0 < θ 1 < θ < < θ M 1 < θ M = π, θ j = j θ and θ = π M, and denote f i,j = φ(r i, θ j ). For each θ j, r φ [ 1 r i 1 r i (f i+ 1,j f i 1,j) r i 1(f i 1 ( r),j f ] i 3,j). That is 1 ( r) r 1 r 1 r 1 r 1 0 r 1 1 r 3 r 3 (r 1 + r ) r r r N 1 (r r N 3 r N 3 N + r N 1 ) r N 1 r N 3 0 r N 1 r N 1 1 r N 1 (r N 1 + r N ) f 1,j f 3,j. f N 3,j f N 1,j Âf, where  = 1 ( r) a 1 c 1 0 b 1 a c b N a N 1 c N 1 0 b N 1 a N and f j = f 1,j f 3,j. f N 1,j, and b i = r i r i+ 1 a i = 1 (r i 1 + r i ) = i 1 + i r i 1 i 1 = i + 1, c i = r i r i 1 =, i = 1,..., N = 1 1, i = 1,..., N 1. i 1 11

20 Now we consider the operator θ φ = 1 r φ θ on θ [0, π]. Then for each r i 1, i = 1,,..., N [ Let and θ φ B = diag C = diag 1 r i 1 ( ( r 1 r 1 f(r i 1, θ j+1 ) f(r i 1 ( θ), θ j ) + f(r i 1 ], θ j 1 ). ) ( θ), ( θ), ( θ),..., r ( θ) N 1 1 ( θ), r 3 r 3 1 ( θ), r 5 r 5 1 ( θ),..., ) 1, r ( θ) N 1 f j = [f 1,j, f 3,j,..., f N 1,j] for j = 1,,..., M. Then we can rewrite the Laplace operator as  + B C C C  + B C C  + B C C C  + B Then the matrix A represents the Laplace operator. 3. Algorithm f 1 f. f M 1 f M AF. 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 R M M, q, m R M, 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 Aq = A[q 1, q,..., q M ] T, 1

21 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 (x 1, x ) and m(x 1, x ), respectively. We solve it by using an iteration method [7]: [ɛa + diag(v + λ)] q new = m q 3 old, (18) where q new and q old 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 q old > 0, and q old = 1. Step Solve (18), then obtain q new. Step 3 Let α new = 1 q new, q new = α new q new. Step 4 if (converge) then Output the solution (α new ) 1 q new. Stop. else Let q old = q new. Go to Step. 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 )]q = λq. (19) It turns into an eigenvalue problem. In Algorithm, it shows the iteration step. In Step II. Now we discretize L of (8) into an eigenvalue problem. [ ] [ ] u u L = ρ, (0) w w where [ L = 0 ɛa + diag(v + λ) diag(m q ) ɛa diag(v + λ) + diag(3m q ) 0 q and λ are obtain from Step I. We use ARPACK in MATLAB version R007a 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. Consider the solution form: ψ(t, x) = e iλt φ(x), λ = 0, 1, and general case. ]. 13

22 Algorithm Iterative algorithm for solving φ(x) in general solution Step 1 Choose an initial guess of q old > 0, and q old = 1. Step Solve (ɛa + diag(v m q old )) q new = λ q new, where λ is the smallest eigenvalue of ɛa + diag(v m q old ). Then obtain q new. Step 3 Let ˆq = (ɛa + diag(v m q new)) q new, α new = λ ˆq, q new = α new q new. Step 4 if (converge) then Output the solution q new and λ. Stop. else Let q old = q new. Go to Step. end if 4.1 Setting In ω case, the discretized matrix of Laplace operator has size NT by NT with N = r/ r and T = π/ θ. Here we use zero boundary condition and r = 1, N = 3, T = 64. And in µ case with square domain, the discretized natrix of Laplace operator has size N by N with N = /h and here N = 64. And so the matrix size of the operator L are NT by NT and N by N respectively. In our numerical method, we use finite-difference method to solve NLS equation. The finite-difference method has truncation error O(h ), where h is the grid size. In our experiment, O(h ) c 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 In Fig. 4, we plot the mean value of the boundary of φ for ɛ from 0 to We can see the mean of boundary will be less than 10 5 as ɛ smaller than So in that case, we only try ɛ [10 5, 0.004] 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 thesis, 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 14

23 Table 1: The region of ɛ. case µ case ω λ = 0 [10 5, 0.016] [10 5, 0.004] λ = 1 [10 5, 0.05] [10 5, 0.005] general λ [10 5, 0.001] [10 5, 0.001] Figure 4: Boundary φ in ω case for λ = 0. 15

24 Figure 5: Boundary φ in ω case for λ = 1. 16

25 Figure 6: Boundary φ in µ case for λ = 0. 17

26 Figure 7: Boundary φ in µ case for λ = 1. 18

27 Figure 8: Solution of φ for ω 1 = 1, ω = 1, ɛ = 0.001, λ = 1. the solution is stable while ɛ < ɛ, unstable when ɛ > ɛ. Since the numerical result of eigenvalues are approximartion of the exact mathematical 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 Since the regions of ɛ are small, and the spectrum of lineraization operation changes so fast. We move ɛ in a small step for each case. δɛ = 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. ω 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 ɛ = 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, x 1 = 1/4. 1. λ = 0: Testing ɛ in the region that we mentioned before, and find that the eigenvalues of the linearized operation L are all pure imaginary.. λ = 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. 1, it shows the ɛ of the parameters ω 1, ω = 1,,..., 5 19

28 Figure 9: Solution of φ for ω 1 = 1, ω = 1, ɛ = 0.001, λ = In elliptic potential case, we can find ɛ in e it φ(x) this solution form in our numerical method, but we can not find the ɛ in other solution form. At first we guess ɛ in e i0t φ(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. 0

29 Figure 10: φ moves for λ = 1, ɛ =

30 Figure 11: Spectrum of ω case for λ = 1.

31 Figure 1: ɛ of ω case for λ = 1. 3

32 Figure 13: Value of λ in ω case. 4

33 Figure 14: Spectrum of µ case for λ = µ case 1. λ = 0, we can find ɛ in this solution case. For µ 1, µ = 1,,..., 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, µ.. λ = 1, we also find ɛ. Figure 16, 17 show the eigenvalues and values of ɛ for each µ 1, µ. 3. Consider the same µ 1 and µ, 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 µ 1 or µ for different periodic of potential. We can see that the value of ɛ increasing when µ 1, µ go larger. To comparison e it φ(x) and e i0t φ(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 ), 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 ɛ < , λ is smaller than 1. And the minimum of trap potential is 1, that means the trap potential is no longer positive value. 5

34 Figure 15: ɛ of µ case for λ = 0. 6

35 Figure 16: Spectrum of µ case for λ = 1. 7

36 Figure 17: ɛ of µ case for λ = 1. 8

37 Figure 18: Value of λ in µ csae. 9

38 5 Conclusions In this thesis, we consider the NLS equation with the special solution e iλ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 ɛ. For future work, we are interested in improving the method of finding the spectrum of the linearization operator. 30

39 References [1] R. Burden and J. D. Faires. Numerical Analysis. Brooks Cole, 7 edition, 000. [] N. Cesa-Bianchi and G. Lugosi. Prediction, Learning, and Games. Cambridge University Press, 006. [3] S. M. Chang, S. Gustafson, K. Nakanishi, and T. P. Tsai. Spectra of linearized operators for NLS solitary waves. SIAM Journal on Mathematical Analysis, 39(4): , 007. [4] M. W. Hirsch, S. Smale, and R. L. Devaney. Differential equations, dynamical systems, and an introduction to chaos. Academic Press, 004. [5] P. Holmes and E. T. Shea-Brown. Stability. Scholarpedia, 1(10):1838, 006. [6] lattice. Optical lattice- wikipedia, the free encyclopedia. [7] 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):61 79, 00. [8] M. C. Lai. A note on finite difference discretizations for poisson equation on a disk. Numerical Methods for Partial Differential Equations, 17:199 03, 001. [9] T. C. Lin and J. Wei. Orbital stability of bound states of semi-classical nonlinear schrödinger equations with critical nonlinearity. SIAM Journal on Mathematical Analysis, 40(1): , 008. [10] T. C. Lin, Juncheng Wei, and Wei Yao. Orbital stability of bound states of nonlinear schrödinger equations with linear and nonlinear lattices. [11] M. P. MacDonald, G. C. Spalding, and K. Dholakia. Microfluidic sorting in an optical lattice. Nature, 46:41 44, 003. [1] C. Sulem and P.-L. Sulem. The nonlinear Schrodinger equation: Self-Focusing and Wave Collapse. Springer,