A numerical study of the spectrum of the nonlinear Schrodinger equation

Hele tekst

(1)A Numerical Study of the Spectrum of the Nonlinear Schrödinger Equation by Carel Petrus Olivier. Presented in fulfillment of the requirement for the degree of. Master of Science (Applied Mathematics) At Department of Applied Mathematics, Faculty of Natural Sciences, Stellenbosch University Supervisors: Prof B.M. Herbst Department of Applied Mathematics University of Stellenbosch Prof I.V. Barashenkov Department of Applied Mathematics University of Cape Town December 2008.

(2) i. Declaration By submitting this thesis electronically, I declare that the entirety of the work contained therein is my own, original work, that I am the owner of the copyright thereof (unless to the extent explicitly otherwise stated) and that I have not previously in its entirety or in part submitted it for obtaining any qualication.. 26 November 2008. c 2008 Stellenbosch University Copyright . All rights reserved.

(3) Abstract The NLS is a universal equation of the class of nonlinear integrable systems. The aim of this thesis is to study the NLS numerically. More specically, an algorithm is developed to calculate its nonlinear spectrum. The nonlinear spectrum is then used as a diagnostic for numerical studies of the NLS. The spectrum consists of a discrete part, further subdivided into the main part, the auxiliary part, and the continuous spectrum. Two algorithms are developed for calculating the main spectrum. One is based on Floquet theory, rst implemented by Overman. [12]. The other is a direct. calculation of the eigenvalues by Herbst and Weideman. [16].. These algorithms. are combined through the marching squares algorithm to calculate the continuous spectrum. All ideas are illustrated by numerical examples.. ii.

(4) Opsomming Die NLS is 'n universele vergelyking, een van die nie-lineˆ ere integreerbare stelsels. Die doel van hierdie tesis is om die NLS numeries te bestudeer.. In die besonder. word 'n algoritme ontwikkel om die nie-lineˆ ere spektrum numeries te bereken. Die spektrum word dan gebruik as 'n instrument vir die numeriese studie van die NLS. Die spektrum bestaan uit 'n diskrete spektrum, en is verder verdeel in 'n hoofspektrum, 'n hulpspektrum, en 'n kontinue spektrum. Twee algoritmes word ontwikkel om die hoofspektrum te bereken.. Die een is gebasseer op Floquet teorie, en was. eerste ge¨ mplimenteer deur Overman. [12]. Die ander is 'n direkte berekening van. die eiewaardes deur Herbst en Weideman [16]. Hierdie algoritmes word gekombineer deur die marching squares algoritme om die kontinue spektrum te bereken. Alle idees word deur numeriese voorbeelde ge¨ llustreer.. iii.

(5) Acknowledgements To my father for his support, my mother for her encouragement and prayers, and my brother Francois for all the `error correcting'.. iv.

(6) Contents Abstract. ii. Opsomming. iii. Acknowledgements. iv. Contents. v. List of Abbreviations. vii. 1 Introduction. 1. 1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1. 1.2. Historical Background. . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. 1.3. Brief Overview. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3. 2 Analytical Properties of the NLS. 5. 2.1. Conservation Laws. . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2.2. The Hamiltonian Structure of the NLS. 2.3. Integrability of the NLS. 2.4. The Spectrum of the NLS. 2.5. 5. . . . . . . . . . . . . . . . . .. 6. . . . . . . . . . . . . . . . . . . . . . . . . .. 8. . . . . . . . . . . . . . . . . . . . . . . . .. 9. Floquet Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 10. 3 The Split-Step Method for the NLS. 13. 3.1. The Fourier Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 13. 3.2. The Semi Discrete NLS. . . . . . . . . . . . . . . . . . . . . . . . . .. 16. 3.3. Hamiltonian Properties of the Semi Discrete System . . . . . . . . . .. 16. 3.4. The First Order Split-Step Method. . . . . . . . . . . . . . . . . . . .. 17. 3.5. The Second Order Split-Step Method . . . . . . . . . . . . . . . . . .. 20. v.

(7) vi. CONTENTS. 4 Numerical Calculation of the Spectrum. 21. 4.1. The Discrete Spectrum Method (DSM) . . . . . . . . . . . . . . . . .. 21. 4.2. Continuous Spectrum Method (CSM) . . . . . . . . . . . . . . . . . .. 25. 4.3. Marching Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 25. 5 The Periodic NLS. 30. 5.1. Plane Wave Solution. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 5.2. Symmetric Perturbation Analysis. 30. . . . . . . . . . . . . . . . . . . . .. 32. 5.2.1. Linear Analysis . . . . . . . . . . . . . . . . . . . . . . . . . .. 32. 5.2.2. Perturbation Analysis. . . . . . . . . . . . . . . . . . . . . . .. 34. 5.3. Numerical Calculation of the Spectrum . . . . . . . . . . . . . . . . .. 36. 5.4. Application of Algorithm 1 . . . . . . . . . . . . . . . . . . . . . . . .. 39. 5.4.1. Homoclinic Structure of the NLS. . . . . . . . . . . . . . . . .. 39. 5.4.2. Asymmetric Perturbations . . . . . . . . . . . . . . . . . . . .. 40. Evolution of the Discrete Spectrum . . . . . . . . . . . . . . . . . . .. 44. 5.5.1. Symmetric and Asymmetric Perturbations. . . . . . . . . . . .. 44. 5.5.2. Homoclinic Chaos . . . . . . . . . . . . . . . . . . . . . . . . .. 45. 5.5. 6 The NLS on the Innite Line 6.1. The Single Soliton. 6.2. Proposed Algorithm. 6.3. Soliton Collisions. 47. . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 48. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 48. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 51. 6.4. The Stationary Bound State Soliton . . . . . . . . . . . . . . . . . . .. 53. 6.5. Soliton Content of Non-Soliton Potentials . . . . . . . . . . . . . . . .. 55. 7 Conclusion. 59. Bibliography. 60. A Proof of Theorems of Floquet Theory. A1. B Symplectic Integrators. B1. C Proof of a Trigonometric Series Identity. C1. D Proof of Theorem 5.1. D1.

(8) List of Abbreviations NLS. Nonlinear Schr¨ odinger Equation. KdV. Korteweg-de Vries Equation. ODE. Ordinary Dierential Equation. IST. Inverse Scattering Transform. DSM. Discrete Spectrum Method. CSM. Continuous Spectrum Method. MSA. Marching Squares Algorithm. FFT. Fast Fourier Transform. IFFT. Inverse Fast Fourier Transform. vii.

(9) Chapter 1 Introduction 1.1. Motivation. In recent times, the importance of nonlinearity has become more and more apparent. It arises in the modelling of a wide variety of subjects like physics, chemistry, biology, economics and transportation. For nonlinear systems, the principal of superposition does not apply. This complicates mathematical analysis. As a result, computers play an important role in the analysis of these systems. Due to the computer revolution of the last 60 years, some major breakthroughs were made in the understanding of nonlinear systems. During this time, computers were used to verify analytical results and to conduct numerical experiments. These experiments led to many important discoveries including chaos and solitons. Solitons are localized waves that propagate at a constant velocity while retaining their shape. Moreover, the shape and speed of solitons are not altered by collisions with other solitons. Solitons are associated with integrable mathematical models. These models include the Korteweg-de Vries equation (KdV), the Nonlinear Schr¨ odinger equation (NLS), and the sine-Gordon equation. In this thesis, we study the NLS for two dierent sets of boundary conditions. The periodic NLS is given by. iψt + ψxx + 2i |ψ|2 ψ = 0, ψ (x + L, t) = ψ (x, t) , ψ (x, 0) = f (x) . The NLS on the innite line is given by. 1. (1.1).

(10) CHAPTER 1.. 2. INTRODUCTION. iψt + ψxx + 2 |ψ|2 ψ = 0, ψ (x, t) → 0 as x → ±∞ ψ (x, 0) = f (x) .. (1.2) fast enough,. The NLS is used to describe a wide variety of physical phenomena. We highlight two of them, namely deep water waves and bre optics. For deep water waves, the NLS is used to describe the envelope of the wave.. For. example, the notorious rogue wave is described by the NLS [1]. It explains how it is possible for these waves to travel long distances without dispersing into the ocean. A very important application of soliton theory is in the eld of bre optics. Here solitons can be used to transmit data at an extremely high rate.. Because of the. increasing demand for information, and the increased availability of information, there is much interest in these developments.. 1.2. Historical Background. Soliton theory has a very interesting history. It all started in 1834 when John Scott Russell discovered the solitary wave. In a famous quote of soliton theory, he describes this discovery:. I was observing the motion of a boat which was rapidly drawn along a narrow channel which it had put in motion; it accumulated round the prow of the vessel in a state of violent agitation, then suddenly leaving it behind, rolled forward with great velocity, assuming the form of a large solitary elevation, a rounded, smooth and well-dened heap of water, which continued its course along the channel apparently without change of form or diminution of speed. I followed it on horseback, and overtook it still rolling on at a rate of some eight or nine miles an hour, preserving its original gure some thirty feet long and a foot to a foot and a half in height. Its height gradually diminished, and after a chase of one or two miles I lost it in the windings of the channel.. [2]. Note that this type of wave is called a solitary wave. He faced a lot of scepticism from the scientic community.. Some of the great scientists of the time, like Airy. and Stokes, raised their doubts about the existence of such a wave. But during the 1870's it was nally accepted. In 1895 Korteweg and De Vries derived an important equation called the Kortewegde Vries equation (KdV). The KdV was derived to model the motion of waves in shallow water. This equation was pivotal in the development of soliton theory [3]..

(11) CHAPTER 1.. 3. INTRODUCTION. The next breakthrough was achieved by Fermi, Pasta and Ulam (FPU) in 1955 during their investigation of heat conductivity of a solid [4]. They used a computer to study a one-dimensional lattice consisting of masses coupled by springs.. They expected. that, after suciently long time, the energy will spread equally across all the modes. To their surprise, the energy ows back and forth between some low-order modes, and eventually returns to the lowest mode. After this, the process repeats itself. This result inspired Kruskal and Zabusky to approach the FPU problem from the continuum perspective. In 1965 they showed that it is related to the KdV [5]. By imposing periodic boundary conditions, they investigated a sinusoidal wave numerically. The initial wave broke up into smaller waves with dierent amplitudes, with the larger amplitude waves traveling at higher velocities. As a result, after suciently long time, these waves overtook the smaller ones (as a result of the periodicity). After the collision the larger wave would reappear in front of the smaller wave. More importantly, the waves maintained their amplitude and velocity after collision. They called these special waves solitons, because of their particle-like collision properties. In 1967, Gardner, Greene, Kruskal and Miura made a very important contribution by solving the KdV analytically [6]. They did this by developing a method called the inverse scattering transform (IST). Zakharov and Shabat extended the IST to the NLS in 1972 [7]. This was done by constructing a Lax pair for the NLS. The corresponding eigenvalue problem, named after them, was generalised by Ablowitz, Kaup, Newell and Segur (AKNS) in 1973 [8, 9]. As a result, the IST was extended to various equations, including the sineGordon equation. During this time considerable advancements were made in numerical analysis. An integrable nite dierence scheme was developed for the NLS in 1975 by Ablowitz and Ladik [10]. In 1986 Herbst and Weideman analysed the split-step method for the NLS. [11]. This method is a very fast pseudo-spectral scheme which preserves the. Hamiltonian structure of the NLS. In the same year Overman developed a computer code to calculate the spectrum of the sine-Gordon equation [12]. It was also used to study the homoclinic structure of the periodic NLS [13, 14, 15]. In this contribution a numerical scheme is developed to calculate the continuous spectrum for the periodic NLS. Instead of using a grid on the complex plane, the discrete spectrum is calculated through an eigenvalue problem, derived by Herbst and Weideman [16]. This is combined with the marching squares algorithm (MSA) to construct the contour of the continuous spectrum. By doing this the number of operations is reduced.. 1.3. Brief Overview. In Chapter 2 we focus on the analytical properties of the NLS. We look at some of its conservation laws, and discuss its Hamiltonian structure and integrability. Finally,.

(12) CHAPTER 1.. 4. INTRODUCTION. we look at the Zakharov-Shabat eigenvalue problem, which is used to dene the spectrum. Chapter 3 is devoted to the development of the pseudo-spectral split-step method. This method is used throughout the thesis to solve the NLS numerically. In order to calculate the spectrum numerically, we develop three numerical methods in Chapter 4. In the following chapters we combine these methods to calculate the spectrum. The rst method is the Discrete Spectrum Method (DSM). The second is the Continuous Spectrum Method (CSM). The third method is the marching squares algorithm (MSA), which is used to combine the DSM and the CSM. Chapter 5 is devoted to the periodic NLS, given by (1.1). We study its homoclinic structure analytically by doing a perturbation analysis on the plane wave solution. This is followed by developing Algorithm 1 which calculates the spectrum based on the methods of Chapter 4. Finally we apply the split-step method and Algorithm 1 to study the homoclinic structure of the periodic NLS numerically. In Chapter 6 we look at solitons for the NLS on the innite line, given by (1.2). Algorithm 2 is developed to calculate its spectrum in this regime.. We implement. it on the single soliton, bound state solitons, and on soliton collisions. We end the chapter by investigating the spectrum of non-soliton potentials, and the possibility of calculating their `soliton content' numerically..

(13) Chapter 2 Analytical Properties of the NLS In order to study the NLS, we need a thorough understanding of its mathematical properties. This gives a good foundation enabling us to perform a numerical study. It also provides a framework to evaluate the accuracy of the numerical schemes. This chapter focuses on three important properties of the NLS with periodic boundary conditions (1.1). The rst property considered is the conservation laws of the NLS. Two conservation laws are derived.. The second property we discuss is the. Hamiltonian structure of the NLS. This implies an underlying symplectic structure for the phase space. The third property is the integrability of the periodic NLS. The integrability of the NLS empowers us to establish the concept of the spectrum of the NLS. The chapter is concluded with a discussion of the spectrum, as well as the development of Floquet theory. Floquet theory is useful for calculating the spectrum of the NLS.. 2.1. Conservation Laws. It is well known that the NLS has innitely many conservation laws [9]. This property is signicant, since it is a necessary condition for the integrability of the NLS. In this section we derive two important conservation laws of the NLS. The rst conservation 2 law shows that the area under the graph of the modulus squared, |ψ| , is preserved for all time. To nd this conservation law, we look for a function ∗ Multiplying (1.1) with ψ we get. I. iψ ∗ ψt + ψ ∗ ψxx + 2 |ψ|4 = 0. If we subtract (2.1) from its complex conjugate, we get. 5. which is constant over time.. (2.1).

(14) CHAPTER 2.. ANALYTICAL PROPERTIES OF THE NLS. ∗ |ψ|2t = iψ ∗ ψxx − iψψxx .. 6. (2.2). x over one period, and using integration by parts. By integrating (2.2) with respect to on the right hand side, we get. dI = 0, dt. (2.3). where. ˆ. L. |ψ|2 dx.. I=. (2.4). 0 Note that the boundary terms, which result from the integration by parts, is zero as a result of the periodicity.. Therefore. I. is a conserved quantity, and the rst. conservation law is given by (2.3) and (2.4). It follows that the area under the graph of the modulus squared is preserved for all time. Our next constant of motion is useful for nding the Hamiltonian of the NLS. Mul∗ tiplying (1.1) with ψt gives. i |ψ|2t + ψt∗ ψxx + 2 |ψ|2 ψψt∗ = 0.. (2.5). If we add (2.5) to its complex conjugate we get. ∗ ψt ψxx + ψt∗ ψxx + |ψ|4. t. = 0.. (2.6). To obtain the conservation law, we need to integrate (2.6) over one period with respect to. x.. If we use integration by parts on the rst two terms we obtain. ˆ 0 dJ which shows that dt. = 0,. L. − |ψx |2 t + |ψ|4 t dx = 0,. where. ˆ. L. |ψx |2 − |ψ|4 dx.. J=. (2.7). 0. 2.2. The Hamiltonian Structure of the NLS. The Hamiltonian structure of dynamical systems have been researched extensively since it rst appeared in the 19th century.. These systems accurately portray a. wide variety of physical systems characterised by energy conservation. An important.

(15) CHAPTER 2.. 7. ANALYTICAL PROPERTIES OF THE NLS. property of Hamiltonian systems is their preservation of a symplectic structure in the phase space. This property is utilised in Chapter 3 to solve the NLS numerically. We start by dening the Hamiltonian for a general innite dimensional system. Let. p (x, t). and. q (x, t). be dened on. 0 < x < L.. Let. pt = f1 (p, px , . . . , q, qx , . . . ) , qt = f2 (p, px , . . . , q, qx , . . . ) . Suppose there exists a function. h (p, px , . . . , q, qx , . . . ) ˆ. (2.8). such that. L. h (p, px , . . . , q, qx , . . . ) dx. H (p, q) =. (2.9). 0 satises the equations. ∂q δH = , ∂t δp. (2.10). ∂p δH = − . ∂t δq. (2.11). Then the system (2.8) is called a Hamiltonian system in canonical form. Here. δH is δp. the variational derivative dened by. H (p + η, q) − H (p, q) lim = →0 where. η. ˆ 0. L. δH η dx, δp. is an arbitrary function. It follows that (see e.g. [17]). n ∞ ∂ ∂h δH (p, q) X = − , δp ∂x ∂pn n=0 where. pn =. (2.12). ∂np δH . The same result applies to . ∂xn δq. In order to show that the NLS is a Hamiltonian system, we need to nd a function H that satises (2.10). Since |ψ|2 = ψψ ∗ , let p ≡ ψ ∗ and q ≡ ψ . If we substitute this into (1.1), and take its complex conjugate, we obtain the following system:. qt = iqxx + 2iq 2 p . pt = −ipxx − 2ip2 q. (2.13).

(16) CHAPTER 2.. 8. ANALYTICAL PROPERTIES OF THE NLS. Since the Hamiltonian is time invariant, the conserved quantities can be tested. It is easily veried that the function NLS. However, dierentiation of. I. J. δJ δp. where of. H. qt. in (2.4) is not the Hamiltonian function for the. in (2.7) yields. ∂h ∂ ∂h + ∂p ∂x ∂px ∂ qx = −2q 2 p − ∂x = iqt , =. is given by (2.13). If we let. H = −iJ ,. (2.14). it follows from (2.14) that this choice. satises (2.10). A similar calculation reveals that. H. also satises (2.11). We. conclude that the NLS (2.13) is a Hamiltonian system in canonical form, with the Hamiltonian given by. ˆ. L. px qx − p2 q 2 dx,. H (p, q) = −i. (2.15). 0 where. p = q∗.. 2.3. Integrability of the NLS. Integrability is a key concept of dierential equations. The integrability of systems places limits on the behaviour of their solutions. It also gives us deeper insight into the characteristics of the equation, and provides us with important tests to measure the accuracy of numerical solutions. Moreover, integrable systems have the property that they can be solved explicitly. Especially for innite dimensional nonlinear systems, these solutions provide us with vital information about the structure of the system. To show that the NLS is integrable, one has to nd a Lax pair corresponding to the NLS. To do this, we need to dene the Lax pair. Consider a general system given by (2.8). depending on a complex coecient also assumed that. ζ. ζ,. Q be two matrix functions, functions p, px , · · · and q, qx , · · · . It is. Let. and the. P. and. is time independent. Let. v=. v1 (x, t) v2 (x, t). .. If (2.8) arises as the compatibility condition of the eigenvalue problems. vx = P v, and. (2.16).

(17) CHAPTER 2.. 9. ANALYTICAL PROPERTIES OF THE NLS. vt = Qv, then the matrices. P. and. Q. (2.17). are called the Lax pair of (2.8). Zakharov and Shabat. [7] showed that the Lax pair of the NLS, written as the system (2.13), is given by. P =. −iζ q −p iζ. ,. (2.18). and. Q= where. q=ψ. and. p = ψ∗.. i(pq − 2ζ 2 ) iqx + 2ζq ipx − 2ζp −i(pq − 2ζ 2 ). ,. (2.19). Note that the NLS solution appears as a potential in the. Lax pair. Therefore solutions of the NLS are also referred to as potentials.. 2.4. The Spectrum of the NLS. The IST plays a very important role in the understanding of the dynamics of the NLS. By using the Lax pair, the NLS can be solved through the IST. To apply the method, we transform the potential evolves linearly over time.. ψ. onto its scattering data. The scattering data. The evolved scattering data can then be transformed. inversely, which yields the evolved potential. For a detailed discussion, see [9]. In this thesis we concentrate on the numerical calculation of the spectrum.. The. spectrum is time invariant, and consists of the continuous and discrete spectrum. The discrete spectrum (main spectrum) is a subset of the continuous spectrum which forms part of the scattering data for the periodic NLS. Therefore a lot of information about the solution of a given potential is contained in the discrete spectrum. The continuous spectrum provides information on the behaviour of the auxiliary spectrum (I) and (II). In [18] it is shown that the auxiliary spectra (I) and (II) form part of the scattering data for the periodic NLS. Section 2.3 shows that the NLS arises as the compatibility condition of the two linear operators. We write (2.16) and (2.17) as. L(x,t) v = 0.. (2.20). Here. L. (x). =. ∂ ∂x. + iζ ψ∗. −ψ ∂ − iζ ∂x. ,. (2.21).

(18) CHAPTER 2.. 10. ANALYTICAL PROPERTIES OF THE NLS. and. L. (t). =. ∂ ∂t. − i(|ψ|2 − 2ζ 2 ) −iψx∗ + 2ζψ ∗. In (2.21) and (2.22) the complex number. ζ. ∂ ∂t. −iψx − 2ζψ + i(|ψ|2 − 2ζ 2 ). .. (2.22). is a parameter.. The spectrum of the NLS is dened in terms of the spectrum of the spatial operator, L(x) . The discrete spectrum and the continuous spectrum are dened as. σc L(x) σd. . =. . ζ ∈ C | L(x) v = 0, |v|. bounded. ∀x ,. (2.23). . L(x) = ζ ∈ C | L(x) v = 0, v (x + 2L) = v (L) ∀ x .. (2.24). In order to proceed we need to divert and develop some Floquet theory.. 2.5. Floquet Theory. Floquet theory investigates rst-order systems of the form. v0 = P (x)v, with the property that the matrix means that. P (x + L) = P (x). for. (2.25). P (x) is periodic with minimal period L. This −∞ < x < ∞. Floquet theory is useful for. determining stability properties. A famous application is for determining the stability regions for Mathieu's equation. We state the main theorems of Floquet theory [19]. The proofs are given in Appendix A.. Theorem 2.1: Given the system with periodic coecient matrix P as in. Then the system has at least one non-trivial solution v = χ(x) such that χ(x + L) = µχ(x), where µ is constant. For any fundamental (solution) matrix. Φ(x). of the above system, all solutions can. be written as a linear combination of the columns of in Appendix A. The values of of the matrix. E,. µ,. .. (2.25). Φ,. as is stated in Proposition 2. corresponding to Theorem 2.1, are the eigenvalues. where. Φ(x + L) = Φ(x)E.. Theorem 2.2: For the system with periodic coecient matrix P as in. (2.26). , let E be given by (2.26). Suppose E has n distinct eigenvalues µ1 , µ2 · · · µn . Then all solutions are bounded whenever |µj | = 1 for all j = 1, 2, · · · , n. Our last theorem relates the matrix. P. in (2.25) with the eigenvalues:. (2.25).

(19) CHAPTER 2.. 11. ANALYTICAL PROPERTIES OF THE NLS. Theorem 2.3: If P is an n × n matrix, then ˆ. L. µ1 µ2 · · · µn = exp. . tr [P (s)] ds ,. 0. where a repeated eigenvalue µ is counted according to its multiplicity. In (2.23) and (2.24) the spectrum is dened in terms of the equation. L(x) v = 0.. This. can be written as the system. . v10 = −iζv1 + ψv2 . v20 = −ψ ∗ v1 + iζv2. (2.27). If we rewrite (2.27) as in (2.25), we get. P = Since. ψ. is. L. −iζ ψ −ψ ∗ iζ. .. (2.28). periodic, we see that Floquet theory applies to (2.27). We dene the. fundamental matrix. M (x; ψ, ζ). as the unique solution matrix satisfying. M (0, ψ, ζ) = I. The eigenvalues. µ1. and. µ2. are the eigenvalues of. E. (2.29) where. M (L; ψ, ζ) = M (0; ψ, ζ)E = E, i.e., the eigenvalues of. M (L; ψ, ζ).. (2.30). Let. 4 (ψ, ζ) = tr (M (L; ψ, ζ)) . Then. µ1. and. µ2. (2.31). can be calculated as the solution of the quadratic equation. λ2 − 4λ + det(M ) = 0. Since tr [P. (s)] = 0. for all. s,. (2.32). it follows from Theorem 2.3 that. µ1 µ2 = 1. From the solution of (2.32) it can be shown that det (M ). (2.33). = 1.. The roots of (2.32). solve the quadratic equation. (λ − µ1 ) (λ − µ2 ) = 0 =⇒ λ2 − (µ1 + µ2 ) λ + 1 = 0.. (2.34).

(20) CHAPTER 2.. 12. ANALYTICAL PROPERTIES OF THE NLS. From (2.32) and (2.34) it follows that. 4 = µ1 + µ2 .. (2.35). p 1 4 ± 42 − 4 . 2. (2.36). The solution of (2.32) is given by. µ1,2 =. ζ to be a spectral point, v has to be bounded. Since µ1 µ2 = 1, we see from Theorem 2.2 that, in order to obtain a bounded solution, we require |µ1 | = |µ2 | = 1. This is only true if 4 is real and −2 ≤ 4 ≤ 2, as is shown From (2.23) we know that for. below:. Complex 4 (ζ) : follows that. v be a bounded eigenvector. From (2.32) and Theorem 2.2 it |µ1 | = |µ2 | = 1. From (2.33) we see that µ2 = µ11 . Let µ1 = c + id. Let. It follows from (2.35) that. 4 = µ1 + µ2 = c + id + Note that. c2 + d 2 = 1. since. |µ1 | = 1.. ∆ = c + id + which shows that. |4 (ζ)| > 2 :. ∆. 1 . c + id. Therefore. 1 c − id = c + id + 2 = 2c, c + id c + d2. is real whenever. v. is bounded.. From (2.36) we see that the square root term is real.. that one of the Therefore. ζ. |µ|. values is larger than 1, which means that. v. This means. is unbounded.. is not part of the spectrum.. |4 (ζ)| < 2 :. From (2.36) we see that the square root term is a non-zero pure imag∗ inary number. This implies that µ2 = µ1 . Since µ1 µ2 = 1, it follows that |µ1 | = |µ∗1 | = |µ2 |. Applying Theorem 2.2, we see that the solutions for will be bounded. Therefore ζ is part of the spectrum.. v. |4 (ζ)| = 2 : From (2.36) µ1 = µ2 , Theorem 2.1. we see that the square root term disappears. shows that periodic. solutions exist. This implies that. ζ. (4 = 2). or anti periodic. Since. (4 = −2). is a discrete spectral point.. In conclusion, we can rewrite (2.23) and (2.24) in terms of the Floquet discriminant, given by. σc L(x) σd L. . = {ζ ∈ C | ∆. is real and. (x). = {ζ ∈ C | ∆. is real and. − 2 ≤ ∆ ≤ 2} , ∆ = ±2} .. (2.37) (2.38).

(21) Chapter 3 The Split-Step Method for the NLS In the study of nonlinear dynamics, numerical methods play a very important role. Numerical solutions provide insight into problems that cannot be solved analytically. When we choose a numerical scheme, it is important to choose a method that is fast and accurate. Sometimes it is possible to choose a method that preserves some of the properties of the original problem. Based on this criteria, the split-step method is an excellent choice for solving the periodic NLS. The split-step method is a pseudo-spectral symplectic numerical scheme. The pseudospectral discretisation produces a Hamiltonian system of ordinary dierential equations (ODEs).. This allows us to use a symplectic integrator, which preserves the. symplectic structure of the solution. We start the chapter by developing discrete Fourier theory. This is used to discretise the spatial dimension of the NLS, which yields a Hamiltonian system of ODEs in Section 3.2 and Section 3.3. From this we derive the split-step method in Section 3.4 and Section 3.5.. 3.1. The Fourier Series. Consider a suciently smooth periodic function it as a Fourier series. ∞ X. f (x) =. f (x) with period L.. fˆn eiµn x ,. We can express. (3.1). n=−∞ where. µn =. 2πn , L. (3.2). and. 1 fˆn = L. ˆ. L. f (x) e−iµn x dx. 0. 13. (3.3).

(22) CHAPTER 3.. 14. THE SPLIT-STEP METHOD FOR THE NLS. We can formally dierentiate (3.1) to get. ∞ X. 0. f (x) = i. µn fˆn eiµn x ,. (3.4). n=−∞ and. f 00 (x) = −. ∞ X. µ2n fˆn eiµn x .. (3.5). n=−∞ The symbol. F. denotes the mapping of the periodic function. f. onto its Fourier coef-. cients, i.e.. F (f ) = fˆ,. (3.6). where. n ˆ f := fˆn. for. o. −∞<n<∞ .. (3.7). The inverse maps the frequency domain back onto the function space, i.e.. F −1 fˆ = f.. (3.8). Using this notation, we can write (3.4) and (3.5) as. f 0 (x) = iF −1 [µF {f (x)}] , f 00 (x) = −F −1 µ2 F {f (x)} ,. (3.9) (3.10). where. n µF {f (x)} = µn fˆn and. n µ F {f (x)} = µ2n fˆn 2. for. for. o −∞<n<∞ ,. (3.11). o. −∞<n<∞ .. Discrete Case To apply the Fourier series numerically, we need a discrete analogue of these results. To do this, we dene a grid on [0, L]. If we choose an even number of grid points L with constant length h = , our grid is given by N. [x0 , x1 , · · · , xN −1 ] ,. N. (3.12).

(23) CHAPTER 3.. where If. xj = jh. fj = f (xj ),. THE SPLIT-STEP METHOD FOR THE NLS. for. 15. j = 0, · · · , N − 1.. then. fj = fj+N. follows from the periodic boundary conditions. The. discrete Fourier transform is given by. N/2−1 2πijn fˆn e N ,. X. fj =. (3.13). n=−N/2 where. N −1 2πijn 1 X ˆ fn = fj e− N . N j=0. (3.14). Dening the vectors. h iT ˆ ˆ ˆ f = f−N/2 , · · · , fN/2−1. and. f = [f0 , · · · , fN −1 ]T ,. we can write (3.14) as. ˆ f = Ff ,. (3.15). where. Fnj. =. 1 − 2πijn e N . N. (3.16). The inverse is given by (3.13), i.e.. f, f = F−1ˆ. (3.17). where. −1 Fjn. =e. 2πijn N. .. (3.18). The fast Fourier transform (FFT) is a fast algorithm to calculate the discrete Fourier coecients. For the split-step method, the FFT is used to compute (3.15). Likewise the inverse fast Fourier transform (IFFT) is used to compute (3.17)..

(24) CHAPTER 3.. 3.2. 16. THE SPLIT-STEP METHOD FOR THE NLS. The Semi Discrete NLS. In this section we discretise the periodic NLS (1.1) with respect to Fourier theory.. Since. ψ. is space periodic for all. Fourier series. ∞ X. ψ (x, t) =. t,. x. by using the. we use (3.13) to write. ψ. as a. ψˆn (t) eiµn x ,. n=−∞ where. 1 ψˆn (t) = L. ˆ. L. ψ (x, t) e−iµn x dx. 0. From (3.5) we get that, for a xed time,. ψxx = −. ∞ X. µ2n ψˆn eiµn x. n=−∞ 2 −1. = −F In order to discretise with respect to. x,. µ F [ψ] ,. (3.19). we choose a grid. [x0 , x1 , . . . , xN −1 ] , where. (3.20). xj = jh, and let Ψj (t) ≈ ψ (xj , t) for j = 0, · · · , N − 1.. If we discretise (3.19),. we get. 2 ˙ j = −iF−1 Ψ µn An + 2i |Ψj |2 Ψj , j for. j = 0, · · · , N −1.. −1 Here Fj is the. j 'th row of F−1 (i.e.. the vector. (3.21). h. −1 ,··· F j,− N 2. An = Fn Ψ (t) , Ψ (t) = [ψ0 (t) , · · · , ψN −1 (t)],. and Fn is the. n'th. 2. (3.22). row of the matrix F.. In conclusion, we discretised (1.1) to obtain a system of ODEs. (0) (3.21), with initial condition Ψ = [f0 , · · · , fN −1 ]T , where fj = f. 3.3. , F−1 j, N −1. This is given by. (xj ) .. Hamiltonian Properties of the Semi Discrete System. As shown in Chapter 2, the NLS is Hamiltonian. It is interesting to note that the Hamiltonian structure is preserved by using a spectral discretisation. To show that (3.21) has a Hamiltonian structure, we need to nd a function singular antisymmetric matrix. J. H (z). such that (3.21) can be written as. and a non-. i. ),.

(25) CHAPTER 3.. THE SPLIT-STEP METHOD FOR THE NLS. z˙ = J∇H (z) , where. z = [Ψ, Ψ∗ ]T .. 17. (3.23). The Hamiltonian of (3.21) is given by (see [20]). N. H(z) = −iN. −1 2 X. µ2n |An |2 + i. N −1 X. |Ψj |4 .. (3.24). j = 0, . . . , N − 1.. (3.25). j=0. n=− N 2 This implies that. ∂H ψ˙ j = ∂ψj∗. and. ∂H ψ˙ j∗ = − ∂ψj. for. This can be veried by dierentiating (3.24), which yields. . N. −1 2 X. N −1 X. !. 2πikn ∂H 1 ∂  + iΨ2j ∂ Ψ∗j 2 = −iN  µ2n An ∗ Ψ∗k e N ∗ ∗ ∂Ψj ∂Ψj N k=0 ∂Ψj n=− N 2   N −1 2 X 2πijn µ2n An e N  + 2i |Ψj |2 Ψj = −i  n=− N 2. . = −iFj−1 µ2n An + 2i |Ψj |2 Ψj ˙ j. = Ψ A similar calculation shows that. ˙ ∗ = − ∂H , Ψ j ∂Ψj. which proves that (3.24) is the Hamil-. tonian of (3.21).. 3.4. The First Order Split-Step Method. The symplectic integrator is a numerical method for solving initial value problems of Hamiltonian systems. It is designed to preserve the symplectic structure of the associated Hamiltonian system. The system (3.23) can be written in terms of the canonical Poisson brackets, given by. n X ∂F ∂G ∂F ∂G {F, G} = − , ∂q ∂p ∂p ∂q i i i i i=1.

(26) CHAPTER 3.. where. F. and. 18. THE SPLIT-STEP METHOD FOR THE NLS. G. are any functions of. z.. It follows that. z˙ = DH z,. (3.26). DH = {·, H} .. (3.27). where. From (3.24) we see that. H = H1 + H2 ,. where N. −1 2 X. H1 = −iN. µ2n |An |2 ,. (3.28). n=− N 2 and. H2 = i. N −1 X. |Ψj |4 .. (3.29). j=0 It is shown in [20] that. DH1. N −1 2 X ∂ 2 2 ∗ ∂ , =i µ n An ∗ − µ n A n ∂An ∂An N. n=−. (3.30). 2. and. DH2 = 2i. N −1 X. 2. |Ψj |. Ψ∗j. j=0. ∂ ∂ 2 − |Ψj | Ψj . ∂Ψ∗j ∂Ψj. (3.31). The formal solution of (3.26) is given by. z (t + τ ) = eτ (DH ) z (t) = eτ (DH1 +DH2 ) z (t) .. (3.32). To solve (3.32) numerically, we split the right-hand side of the equation which gives. z (t + τ ) = eτ DH1 eτ DH2 z (t) + O τ 2 . . The. O (τ 2 ). term results from the fact that. DH1. and. DH2. (3.33). do not commute, and is. derived in Appendix B. To solve (3.33), we need to solve the following two equations:. z˜ (t + τ ) = eτ DH1 z (t) , z (t + τ ) = eτ DH2 z˜ (t + τ ) .. (3.34) (3.35).

(27) CHAPTER 3.. 19. THE SPLIT-STEP METHOD FOR THE NLS. We start with (3.34) by writing it in dierential form as. z˙ = DH1 z. It follows from (3.30) that the. An. and. A∗n. terms decouple. Therefore. A˙ n = −iµ2n An , which is solved analytically i.e.. 2. An (t + τ ) = An (t) e−iµn τ .. (3.36). From (3.36) we obtain the solution of (3.34), given by. ˜ j (t + τ ) = F−1 Ψ j. ( X. ) −iµ2n τ. (Fnj Ψj (t)) e. .. (3.37). j To solve (3.35), we also rewrite it in dierential form as. z˙ = DH2 z. From (3.31) it follows that the. Ψj. and the. Ψ∗j. terms decouple as well.. Therefore. (3.35) reduces to. ˙ j = −2i |Ψj |2 Ψj . Ψ. (3.38). From (3.38) it follows that. d |Ψj |2 = 0. dt Therefore we can solve (3.38) analytically which gives. 2. Ψj (t + τ ) = Ψj (t) e−2i|Ψj | τ .. (3.39). This scheme is called the rst order split-step method. We can summarise it in the following steps.. Step 1: A (t) = FΨ (t) Step 2: A (t + τ ) = A (t) exp {−iµ2 τ }.

(28) CHAPTER 3.. 20. THE SPLIT-STEP METHOD FOR THE NLS. ˜ (t + τ ) = F−1 A (t + τ ) Step 3: Ψ

(29)

(30) 2

(31) ˜ ˜ Step 4: Ψ (t + τ ) = Ψ (t + τ ) exp −2i

(32) Ψ (t + τ )

(33)

(34) τ. 3.5. The Second Order Split-Step Method. To derive the second order split-step method, we split. eDH. as. eDH = e(c1 +c2 )DH1 +(d1 +d2 )DH2 , where. c1 + c2 = 1 = d 1 + d 2 .. (3.40). The second order split-step scheme then becomes. τ τ Ψ (t + τ ) = e 2 DH1 eτ DH2 e 2 DH1 Ψ (t) + O τ 3 . This and higher order symplectic schemes are derived in the Appendix B.. (3.41).

(35) Chapter 4 Numerical Calculation of the Spectrum This chapter is devoted to solving the spectrum numerically. To do so, we develop three numerical methods, namely the discrete spectrum method (DSM), the continuous spectrum method (CSM), and the Marching Squares Algorithm (MSA). Both the DSM and the CSM are based on the Floquet theory described in Chapter 2. In Section 4.1 we develop the DSM. The DSM calculates the discrete spectrum of a potential by solving an eigenvalue problem. Since the discrete spectrum forms part of the scattering data, the results of the DSM provides important information about the solution of the potential. In addition, the discrete spectrum forms the boundaries of the continuous spectrum.. However, the DSM cannot be used to calculate the. continuous spectrum. In Section 4.2 we develop the CSM. For a potential. ψ,. the CSM calculates the value. of the Floquet discriminant ∆ (ζ) for any ζ ∈ C. One can use this to nd values 2 of ζ for which 0 ≤ |∆| ≤ 4, which is the continuous spectrum as dened in (2.23). Therefore the CSM can be used to calculate the continuous spectrum. Section 4.3 is dedicated to the MSA, which is combined with the CSM and the DSM to calculate the continuous spectrum in Chapters 5 and 6.. 4.1. The Discrete Spectrum Method (DSM). In Chapter 2 we saw that the denition of the discrete spectrum for a given potential,. ψ,. is given by. . σd (ψ) = ζ ∈ C | L(x) v = 0, v (x + 2L) = v (L) ∀ x , where. 21. (4.1).

(36) CHAPTER 4.. 22. NUMERICAL CALCULATION OF THE SPECTRUM. L. (x). =. d dx. −ψ ∂ − iζ ∂x. + iζ ψ∗. .. (4.2). This amounts to looking for eigenvalues of (4.2) corresponding to periodic and antiperiodic eigenvectors. We can rewrite (4.1) in the form of an eigenvalue problem,. . Dx −ψ −ψ ∗ −Dx. . v1 v2. . = −iζ. v1 v2. with periodic or antiperiodic boundary conditions. Here dv ferentiation operator Dx is given by Dx : v → . dx. ,. v = [v1 , v2 ]T ,. x = [x0 , x1 , · · · , xN −1 ]T. We discretise (4.3) by choosing a grid. (4.3). and the dif-. on the interval. [0, L]. as dened by (3.12). Let. . D −Ψ −Ψ∗ −D. . v1 v2. . v1 v2. = −iζ. vi = [vi,0 , vi,1 , · · · , vi,N −1 ]T , vi,j = vi (xj ) the diagonal matrix where Ψjj = ψ (xj ).. Here is. . for. i = 1, 2. .. (4.4). and. j = 0, · · · , N − 1. Ψ. In order to reduce (4.4) to a discrete version of (4.3), we need a matrix. ∂vi (xj ), where (Dvi )j is the j 'th row of the vector Dvi . (Dvi )j ≈ ∂x dene D such that the periodic boundary conditions are imposed. Once we have determined the matrix problem for matrices.. D,. D. such that. We also want to. (4.4) reduces to a standard eigenvalue. There are many numerical methods available to solve this. type of problem. MATLAB's. eig. function was used to solve for the eigenvalues of. (4.4) numerically throughout the thesis. The choice of the matrix. D. plays a very important role in the accuracy of our. results. We show two methods to construct. D.. The rst one uses Finite Dierences,. the second one uses Fourier methods, which is essentially a pseudo-spectral method.. Central Dierences y be a periodic function dened on [0, L]. yj = y (xj ) . From the central dierence formula. Let. yj0 ≈. Let. y = [y0 , y1 , · · · , yN −1 ]T ,. we have. yj+1 − yj−1 . 2h. Because of the periodic boundary conditions, we have result, if we use. D. in (4.4), the solutions of. as a system of equations we get. where. v1. and. v2. yj = yj+N. for all. j.. As a. are periodic. By writing this.

(37) CHAPTER 4.. 23. NUMERICAL CALCULATION OF THE SPECTRUM. . 0 1 −1 0 0 −1.     1   0 D= 2h   ...    0 1. 0 1 0 ... 0 . . .. .. ... 0 0. ··· ··· 0. ··· ··· ···. .. ... .. ... .. ... .. ... .. 0 0 1 .. ... .. ... .. .. ... .. ···. ···. −1 0 0 −1. −1 0 0. .     .  .  . . 0    1  0. Although this matrix is easy to construct, its results are not suciently accurate.. Spectral Method Let us derive a matrix. D. based on the Fourier transform. In Chapter 3 we saw. yj. can be approximated with the discrete Fourier series. N. yj =. −1 2 X. yˆn e. 2πijn N. ,. (4.5). n=− N 2 and the derivative. yj0. with. N. −1 2 X 2πijn 2πin 0 yˆn e N , yj = L N n=−. (4.6). 2. where. yˆn = L = 2π . In order to retain yˆ− N = 0. (4.6) then becomes Let. 2. N −1 2πijn 1 X yj e− N . N j=0. the symmetry of the Fourier coecients, we set.

(38) CHAPTER 4.. 24. NUMERICAL CALCULATION OF THE SPECTRUM. N. yj0 =. −1 2 X. inˆ yn e. 2πijn N. n=− N +1 2 N. 1 = N. =. By using the identity. yj0 As a result,. yj0. 1 N. −1 2 X. N −1 X. in. n=− N +1 2 N −1 X. . ! − 2πikn N. e. yk e. 2πijn N. k=0. . N. −1 2 X. ine. . 2πi(j−k)n N.  yk .. n=− N +1 2. k=0. eix − e−ix = 2isinx. we get.  N N −1 −1 2 2 X X 2 2π (j − k) n   = − nsin yk . N k=0 n=1 N yk. becomes a linear combination of the. (4.7). terms. In Appendix C it is. shown that. N −1 2 X 2πn (j − k) − nsin = N n=1 N. 1 2. (−1)j−k cot (j−k)π N 0. if if. m 6= 0 . m=0. (4.8). If we apply this to (4.7), we get that. yj0 =. N −1 X. Djk yk ,. (4.9). k=0 where. Djk = The. Djk. 1 2. (−1)j−k cot (j−k)π N 0. values form the dierentiation matrix. D,. if if for. j 6= k . j=k j, k = 0, · · · , N − 1.. We can extend this to a function of arbitrary period. If y (x) is dened on we can make the linear transformation x ˜ = 2π x. Since x˜ ∈ [0, 2π], we get L. Therefore. (4.10). x ∈ [0, L],. x dy dy d˜ 2π dy = = . dx d˜ x dx L d˜ x. (4.11). dy 2π ≈ Dy. dx L. (4.12).

(39) CHAPTER 4.. 4.2. 25. NUMERICAL CALCULATION OF THE SPECTRUM. Continuous Spectrum Method (CSM) ψ,. in terms of the. − 2 ≤ ∆ ≤ 2} .. (4.13). In (2.37), we dened the continuous spectrum of a potential, Floquet discriminant. σc (ψ) = {ζ ∈ C | ∆ (ζ). is real and. In this section we calculate the Floquet discriminant,. ∆,. as a function of. ζ.. Recall. from Chapter 2 that. 4 (ψ, ζ) = tr (ML ) , where. ML = M (L; ψ, ζ) ,. and. M. . where. (4.14). is the fundamental matrix solution for the system. v10 = −iζv1 + ψv2 , v20 = −ψ ∗ v1 + iζv2. (4.15). M (0; ψ, ζ) = I .. Therefore the rst column of. ML. is the solution of (4.15) at. x = L,. where the initial. condition is given by. . If the function. ψ. v1 (0) v2 (0). . =. 1 0. .. (4.16). is known, we can integrate (4.15) numerically from. yields the rst column of. ML .. To calculate the second column of. 0. ML ,. to. L,. which. we integrate. (4.15) with initial condition. . v1 (0) v2 (0). . =. 0 1. .. Using the initial conditions (4.16) and (4.17), the discriminant of. (4.17). ζ. is given by (4.14).. There are many numerical methods available to calculate the initial value ODEs described above. MATLAB's. ode45 function was used to calculate (4.15) throughout. the thesis.. 4.3. Marching Squares. The MSA is an algorithm designed to calculate a contour on a 2D domain. It was originally developed for Image Processing to calculate contours numerically.. The. idea behind the algorithm is to construct rectangles (or squares), and to determine where the contour intersects their sides..

(40) CHAPTER 4.. 26. NUMERICAL CALCULATION OF THE SPECTRUM. We explain the method through the following example: plane, dened by the points. (xi , yj ),. where. Given a grid on the. i, j = −3, · · · , 3.. Let. xj = yj = 0.3j .. xy In. addition, suppose we are given the function value at each point, where the function is dened by. f (x, y) = x2 + y 2 .. (4.18). Suppose we want to use this information to nd the contour given by. f = 1.. (4.19). The rst step of the MSA is to distinguish between points on dierent sides of the contour. We do this in Figure 4.1 by marking the grid points satisfying. f >1. with. solid dots, and points on the other side with open dots. It is important to note that Figure 4.1 shows the domain of. f.. 1.5. 1. y. 0.5. 0. −0.5. −1. −1.5 −1.5. −1. −0.5. 0. 0.5. 1. 1.5. x Figure 4.1:. Points on Dierent Sides of the Contour. To initiate the algorithm, two neighboring points on dierent sides of the contour are required. Let us choose the points. (−1.2, 0). and. (−0.9, 0). (at 9 o'clock on the. circle in Figure 4.1). We use these points to form our initial square, as is shown in Figure 4.2..

(41) CHAPTER 4.. NUMERICAL CALCULATION OF THE SPECTRUM. 27. 0.4 0.35 0.3 0.25. y. 0.2 0.15 0.1 0.05 0 −0.05 −0.1 −1.3. −1.2. −1.1. −1. −0.9. −0.8. x Figure 4.2:. Initial rectangular grid. From Figure 4.2 we see that the contour passes through both the top and bottom sides of the square. One can determine this by looking at the pattern of the four dots.. Since the top left and right corners are on dierent sides of the contour, it. follows that the contour intersects the line between the points. The same holds for the corners on the bottom of the square. We use this information to construct a new square in such a way that it is intersected by the contour. We can do this either in an upward or downward direction. For this example, we follow the contour upwards. The result is that the contour is tracked in a clockwise direction. Since it is known that the contour intersects the side on the top of the square, we can use the two top corners to construct a new square. This process is repeated until the entire contour is tracked. seven steps.. Figure 4.3 shows the rst.

(42) CHAPTER 4.. NUMERICAL CALCULATION OF THE SPECTRUM. 28. 1.2. 1. y. 0.8. 0.6. 0.4. 0.2. 0. −1.2. −1. −0.8. −0.6. −0.4. −0.2. 0. x Figure 4.3:. Seven iterations on the circular contour. For the clockwise direction we use Figure 4.4 to choose the location of the next square. This choice ensures that the next square is intersected by the contour. For example in Figure 4.3 we see that the direction of the `marching' is either right or upwards. We can see that all the squares in Figure 4.3 correspond to those in the rst and last columns of Figure 4.4. For any rectangular set of points, the direction is determined by comparing it to Figure 4.4.. In some applications it may be necessary to follow the contour in both directions. The easiest way to change the direction is to change the sign of the equation. We illustrate this through the example above.. The contour. f =1. is equivalent to the contour given by. g = −1, where. g = −f .. (4.20). But when we apply the MSA, the dierence between (4.19) and. (4.20) is the following: If. f > 1,. g < −1. f < 1. The. it follows that. dots become open dots. The same holds for. Therefore all the coloured result is that the direction. changes from clockwise to anti-clockwise (see Figure 4.4)..

(43) CHAPTER 4.. NUMERICAL CALCULATION OF THE SPECTRUM. 2. 2. 2. 2. 0. 0. 0. 0. −2 −2. 0. 2. −2 −2. 0. 2. −2 −2. 0. 2. −2 −2. 2. 2. 2. 2. 0. 0. 0. 0. −2 −2. 0. 2. −2 −2. 0. 2. −2 −2. 0. 2. −2 −2. 2. 2. 2. 2. 0. 0. 0. 0. −2 −2. 0. 2. −2 −2. 0. Figure 4.4:. 2. −2 −2. 0. 2. −2 −2. MSA direction guide. 0. 2. 0. 2. 0. 2. 29.

(44) Chapter 5 The Periodic NLS Dynamical systems with homoclinic orbits present a numerical challenge. Even for the two dimensional case, numerical solutions near the homoclinic orbit can become chaotic. This is the result of the sensitivity associated with the homoclinic orbit. It was shown that the periodic NLS (1.1) has a homoclinic structure [13, 14, 15]. This is the perfect setting to test the numerical methods derived in Chapters 3 and 4. The plane wave solution forms the basis of the chapter. Section 5.1 focuses on the plane wave solution and its spectrum. This is followed by Section 5.2 which involves a perturbation analysis performed on the plane wave solution. These perturbations are symmetric. The results of the perturbation analysis point us in the direction of the homoclinic structure associated with the periodic NLS. Section 5.3 is dedicated to the development of Algorithm 1, which calculates the continuous spectrum of the periodic NLS. This is done by combining the DSM, the CSM and the MSA of Chapter 4. In Section 5.4 we study the implications of the perturbation analysis of Section 5.2 numerically. This is done by using the split-step method of Chapter 3 to solve the NLS numerically.. We also implement Algorithm 1 to investigate the homoclinic. structure. The evolution of the spectrum is examined in Section 5.5.. This is the result of. numerical errors, since the analytical spectrum is time invariant.. The homoclinic. chaos of the NLS is demonstrated by this evolution.. 5.1. Plane Wave Solution. The plane wave solution is given by. 2 ψ˜ (x, t) = ae2ia t ,. 30. (5.1).

(45) CHAPTER 5.. 31. THE PERIODIC NLS. a is

(46) any positive real number. We can see that the solution is space independent

(47)

(48) ˜

(49) that

(50) ψ (x, t)

(51) = a for all x and t.. where and. For further analysis, we need to calculate the discrete spectrum of. ψ˜.. It turns out. that it can be solved explicitly. This is done by calculating the Floquet discriminant,. ∆ (ζ). In Section 4.2, the discriminant is dened as. 4 (ψ, ζ) = tr (ML (ψ)) . In (5.2),. Mx (ψ). (5.2). is the unique fundamental solution matrix of (4.15) at. x,. satisfying. M0 (ψ) = I, where. ψ. (5.3). d ∆= dt 0. For the plane wave solution,. is the potential. In [15] it is shown that. can be calculated at. t=. 0. Therefore the spectrum Mx (a) is the fundamental. solution of. Lv = 0,. (5.4). where. L=. 1 0 0 −1. . ∂ −a ∂x. . 0 1 1 0. . + iζ. 1 0 0 −1. .. (5.5). Note that (5.5) is a linear operator. The general solution of (5.4) is given by. . where. v1 (x) v2 (x). k 2 = a2 + ζ 2 .. . = c1. 1 i (ζ + k) a. . ikx. e. +. 1 i (ζ − k) a. . e−ikx ,. (5.6). From (5.6), the initial condition (5.3) yields. . iζ  coskx − k sinkx Mx (a; ζ) =  a − sinkx k.  a sinkx  k . iζ coskx + sinkx k. (5.7). Therefore, by (5.2), the Floquet discriminant is given by. p 2 2 a +ζ L . ∆ (ζ) = 2cos Since the continuous spectrum is given by. 0 ≤ ∆2 ≤ 4,. uous spectrum for the plane wave solution is given by. (5.8). (5.8) shows that the contin-.

(52) CHAPTER 5.. 32. THE PERIODIC NLS. σc = {ζ |ζ ∈ R ∪ [−ia, ia]} . The discrete spectrum is given by spectrum is the set of. ζ 's. ∆ = ±2,. (5.9). and from (5.8) follows that the discrete. satisfying cos. p. a2+ ζ 2 L = ±1.. (5.10). This set is given by.

(53) πn 2

(54) σd = ζ

(55)

(56) ζ 2 = − a2 L. for. n = 0, 1, · · ·. .. (5.11). For our analysis we are particularly interested in imaginary double points. A double point belongs to the set. σdouble.

(57)

(58) d∆ =0 . = ζ

(59)

(60) ∆ (ζ) = 2, dζ. (5.12). From (5.8) it follows that the imaginary double points are given by.

(61) πn 2

(62) − a2 ζ

(63)

(64) ζ 2 = L. 5.2. where. Lπ n< , n∈N a. .. (5.13). Symmetric Perturbation Analysis. 5.2.1 Linear Analysis We investigate perturbations of the form. ψ (x, t) = ψ˜ (x, t) (1 + (x, t)) , where. ψ˜. is given by (5.1) and. || 1.. (5.14). If we substitute (5.14) into (1.1) and retain. linear terms only, we get that. t = ixx + 2ia2 ( + ∗ ) . We perturb only the. n'th. (5.15). mode by choosing. (x, t) = ˆ−n (t) e−iµn x + ˆn (t) eiµn x ,. (5.16). 2πn . By substituting (5.16) into (5.15), and using the linear independence L of the exponential functions, we get that for. µn =.

(65) CHAPTER 5.. 33. THE PERIODIC NLS. . (ˆn )t ˆ∗−n t. . =. 2ia2 − iµ2n 2ia2 −2ia2 −2ia2 + iµ2n. . ˆn ˆ∗−n. .. (5.17). The eigenvalues of this system are given by. λ± = ±µn. p. 4a2 − µ2n .. (5.18). The linear system (5.17) is unstable whenever the real part of the eigenvalue is positive.. From (5.18) it follows that the perturbation (5.14) is linearly unstable. whenever. µ2n < 4a2 .. (5.19). We see that this inequality is the same as the condition for imaginary double points in (5.13). It follows that (5.14) is unstable whenever there exists imaginary double points. From (5.18) it also follows that, for unstable perturbations, the plane wave solution contains a stable and an unstable manifold.. To nd these manifolds, we need to. analyse the eigenvectors. The eigenvectors of (5.17), corresponding to. v± =. λ± ,. are given by. µ2n − iλ± µ2n + iλ±. .. We can solve the linear system (5.17), which yields. . ˆn (t) ˆ−n (t). where. r±. . and. λ t 2 λ t 2 iφ+ iφ− + − (µ − iλ ) r e e + (µ − iλ ) r e + + − − n n λ t 2 e λ t , = 2 −iφ+ −iφ− + (µn − iλ+ ) r+ e e + (µn − iλ− ) r− e e − φ± are arbitrary constants. If we choose one of r± = 0 , eλ∓ . In this case we have. (5.20). then both. ˆ±n. depend only on. . ˆn (t) ˆ−n (t). . λ t 2 iφ± ± (µ − iλ ) r e ± ± n e λ t . = 2 −iφ± (µn − iλ± ) r± e e ±. If we substitute this into (5.16), we get. (x, t) = 2r± µ2n − iλ± cos (µn x + φ± ) eλ± t .. (5.21). By substituting (5.21) into (5.14) it follows that the initial conditions corresponding to the stable and unstable manifolds are given by. ψ (x, 0) = a + 0 µ2n − iλ± cos (µn x + φ± ) .. (5.22).

(66) CHAPTER 5.. 34. THE PERIODIC NLS. 5.2.2 Perturbation Analysis To further investigate (5.22), we turn to the discrete spectrum and its corresponding eigenvectors.. From (5.19) we know that the instability is associated with the. imaginary eigenvalues. Consider a perturbed initial condition of the form. ψ (x, 0) = a + exp (iφ) cos (µn x) .. (5.23). Since we want to investigate the splitting of a double point, the appropriate singular perturbation expansion is given by. 1. ζ = ζ0 + 2 ζ1 + ζ2 + · · · ,. (5.24). and. 1. v = v(0) + 2 v(1) + v(2) + · · · .. (5.25). If we substitute (5.24) and (5.25) into (5.4), and gather terms of various orders, we get. Lv(0) = 0 ( Lv(1) =. ( Lv(2) =. where. L. at. O 0. . (5.26). (0). −iζ1 v1 (0) −iζ1 v2 (0). at. 1 O 2. (0). (5.27). (1). −iζ2 v1 + q1 v2 − iζ1 v1 (0) (0) (1) −iζ2 v2 + q1∗ v1 − iζ1 v2. at. O 1 ,. (5.28). is given by (5.5), and. q1 = exp (iφ) cos (µn x) .. (5.29). Equation (5.26) is the same eigenvalue problem as the plane wave equation's. From (5.13) we know that. ζ02 = kn2 − a2 , kn = πn . In addition, we need to solve for the eigenfunctions. L (0) solution of v is given by (5.6), which can be written as where. (5.30) The general.

(67) CHAPTER 5.. 35. THE PERIODIC NLS. v(0) = A+ φ+ + A− φ− .. (5.31). In (5.31),. ±. φ =e and. A±. ±ikx. . . 1 i (ζ ± k) a. ,. (5.32). are arbitrary numbers.. We need to derive a solvability condition in order to solve (5.27) and (5.28). This is done in the following theorem:. Theorem 5.1: The solvability condition for Lv = F, where L is given by F1 F2. F=. (5.30). and. , is given by. ˆ. L. ± F 1 φ± 2 + F2 φ1 dx = 0.. (5.33). 0. The proof is shown in Appendix D. Since (5.27) depends only on. e±ikx , its solvability. condition becomes. . 2ζ0 L ζ1 a. 0 It follows from (5.34) that. ζ1 = 0.. 0 2ζ0 L ζ1 a. . A+ A−. Therefore all. . ζ1. = 0.. (5.34). terms in (5.28) are zero. As a. result, the solvability condition of (5.28) becomes. ". 2. 2ζ0 n) eiφ + 12 e−iφ ζ − 12 (ζ+k a2 a 2 2 2ζ0 1 (ζ−kn ) iφ ζ − 2 a2 e + 12 e−iφ a 2. For non-trivial solutions of. A± ,. #. A+ A−. = 0.. (5.35). we require that the determinant of the matrix in. (5.35) is zero. This yields. 16ζ02 2 2 2 2 ζ = 2 cos (2φ) − k + ζ . 0 2 n a2 a2 In (5.36),. ζ0. (5.36). refers to real or imaginary double points. Consider a real double point ζ0 . This implies ζ02 > 0. Since ζ02 = kn2 − a2 , it follows that a2 < kn2 . In this case, (5.36) becomes.

(68) CHAPTER 5.. 36. THE PERIODIC NLS. 1 8ζ02 1 = 8ζ02. ζ22 =. 2 a cos2φ − kn2 + ζ02. 2 a cos2φ − 2kn2 − a2 " 2 # a2 kn = cos2φ + 1 − 2 . 2 8ζ0 a Since. a2 < kn2 ,. it follows that. ζ22 < 0. for any choice of. φ.. This shows that the real. double points only split in the imaginary direction, which implies that there is no homoclinic boundary associated with the real double points, as is shown in [14]. In contrast, consider an imaginary double point 2 2 2 zero, it follows that (Imζ0 ) + kn = a . If we let tanθ. =. Imζ0. kn. ζ0 .. Since. ζ02 = kn2 − a2. is less than. ,. (5.37). then (5.36) becomes. and. From (5.38) we see that. 16ζ02 2 2 2 ζ = 2 cos (2φ) − 2 cos θ − sin θ 2 a2 = 2cos (2φ) − 2cos (2θ) = −4 sin (φ + θ) sin (φ − θ) a2 ζ22 = − 2 sin (φ + θ) sin (φ − θ) . 4ζ0 ζ2. can be real or purely imaginary, depending on. associated with a homoclinic boundary. If. φ = θ, then ζ2 = 0.. (5.38). φ.. This is. This is a perturbation. in the direction of the unstable manifold. To investigate these results numerically, we develop an algorithm based on the methods of Chapter 4.. 5.3. Numerical Calculation of the Spectrum. The spectrum is dened in terms of the function is dened in terms of the domain of ∆, namely T (x) as the set σ L = A1 A2 , where. ζ.. ∆ : C ⇒ C,. where the solution set. We can write the spectrum (4.13). A1 = {ζ ∈ C | Im [4 (ζ)] = 0} ,.

(69) CHAPTER 5.. 37. THE PERIODIC NLS. and. A2 = {ζ ∈ C | −2 ≤ Re [4 (ζ)] ≤ 2} . To visualise this, consider Figure 5.1. The set. A1. yields a contour in the complex. plane. This is represented by the straight line across the gure. The set. A2. yields a. surface. This is represented by the shaded area. The spectrum is the intersection of the two sets, marked with the bold line. To calculate the spectrum, we implement a MSA on the set. A1 .. We then calculate. Re [4 (ζ)], which shows us whether the point on the contour belongs to the spectrum. In terms of Figure 5.1, we use the real part to determine whether the point on lies in the shaded area. A2 .. 3.5. This is summarised in Algorithm 1.. Re(∆) = 2. Spectrum. 3 2.5. Im ζ. 2 1.5 1. Re(∆) = − 2. 0.5 0. Imag(∆) = 0 −0.5 −0.5. 0. 0.5. Figure 5.1:. 1. 1.5. Re ζ. 2. 2.5. 3. Continuous spectrum. Algorithm 1 Continuous Spectrum Use the DSM to calculate imaginary discrete spectral points for j=1:n Find. A1. near. ζj. using Step 1. Determine appropriate direction using Step 2. |∆| ≤ 2 or |∆|new < |∆|prev Use Step 3 to determine new point on. while end end. A1. ζ1 , · · · , ζn. A1.

(70) CHAPTER 5.. 38. THE PERIODIC NLS. Let us take a closer look at the Steps 1-3 of Algorithm 1:. Step 1: Find points on dierent sides of the contour Im(∆) = 0 In order to start the MSA, we need two points on dierent sides of the contour In theory, the points calculated by the DSM are on. A1 .. A1 .. However, due to numerical. errors this is often not the case. Since these dierences are relatively small, we can use an iterative procedure to nd the contour. A1 .. For example, consider Figure 5.2. Our initial point is denoted by the open dot. We choose two points, one above and one to the right of this point, marked with `x'. The distances of these points are determined by the grid size of the MSA. To nd out if these points are on dierent sides of the contour, we calculate each point's corresponding. ∆ value.. This is done by using the CSM. Figure 5.2 corresponds. to a case where all three points lie on the same side of. A1 , which means that the sign ∆ values to t straight lines. of Im (∆) is the same for all three points. We use these. horizontally and vertically. This is represented by the dotted lines, which is used to estimate Im (∆). = 0,. marked by the black dot. The process is continued until two. points on dierent sides of. A1. are found.. 0.6 0.55 0.5 0.45. Im ζ. 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.4. 0.45. Figure 5.2:. 0.5. 0.55. 0.6. 0.65. Re ζ. 0.7. 0.75. 0.8. 0.85. Fitting a straight line on a discrete point.

(71) CHAPTER 5.. 39. THE PERIODIC NLS. Step 2: Determine the Direction Before we start the MSA iteration, we need to determine whether the orientation is positive or negative.. A positive orientation refers to a clockwise direction.. example, in Figure 5.1, if we generate any point on the contour toward the shaded area. A2. A1 ,. For. the direction. needs to be followed.. To do this, we use the two points calculated in Step 2 to form an initial rectangle. We use this to march a single step in each direction. In both cases we calculate the value of Let. ∆. ∆+. for the point on. A1 .. This is done by using linear interpolation.. be the discriminant of the point with positive orientation, and. negative one.. |∆+ | < |∆− |,. If both. |∆± | < 2,. we need to follow both directions.. ∆−. the. Otherwise if. we follow a positive orientation and vice versa.. Step 3: Apply the MSA At this stage we found a rectangle around the contour direction to follow.. A1 ,. and established which. In this step we apply the MSA to nd the new points in the. domain. For each new point, we calculate the discriminant, which determines if the point is a part of the spectrum.. We also use it to decide whether to continue or. suspend the marching.. 5.4. Application of Algorithm 1. In this section we apply Algorithm 1 to study the spectrum numerically. The solutions are shown for two types of initial conditions. The rst type is given by (5.23), i.e. symmetric perturbations. These perturbations are used to demonstrate the homoclinic structure of the NLS. The second type of initial conditions uses asymmetric perturbations that was studied in. [21].. 5.4.1 Homoclinic Structure of the NLS The results of Section 5.2 show the existence of a homoclinic structure in the periodic NLS. Indeed, analytical expressions of these homoclinic orbits are derived in [22] by using Hirota's method. A useful 2D analogue of the homoclinic structure is given by Dung's equation of the form. x0 = y y 0 = x (1 − x2 ) .. (5.39).

(72) CHAPTER 5.. 40. THE PERIODIC NLS. The phase space of (5.39) is shown in Figure 5.3. an intuitive guide.. Note that we use it only as. For a detailed discussion of the relationship between Dung's. equation and the periodic NLS, see. [23].. From Figure 5.3 it follows that:. •. The xed point at the origin has a stable and an unstable manifold, i.e. the xed point is a saddle point.. •. The stable and unstable manifolds are connected to form a homoclinic orbit. •. A small perturbation around the origin is either inside, outside or on the homoclinic manifold. This is determined by the direction of the perturbation.. For the NLS, the plane wave solution plays the role of a xed point. The plane wave solution has a stable and an unstable manifold according to (5.18).. Therefore it. acts as a saddle point. In the perturbation analysis we derived the initial condition (5.23). This initial condition has a parameter. φ,. which determines the `direction' of. the perturbation. To illustrate the homoclinic structure of the NLS, consider (5.23). If the period of √ π in (5.37). In Figure 5.4, we show (1.1) is given by L = 2 2π , it follows that θ = 4 −3 the results for = 10 . Here the discrete spectrum is indicated by `x'. We see that the results for the spectrum match our analytical results. Figure 5.4 (a) shows the evolution inside the homoclinic orbit.. The period of the. solution in Figure 5.4 (b) is approximately twice as big as the period of Figure 5.4 (a). This corresponds to a solution outside the homoclinic orbit. Figure 5.4 (c) shows a solution near the homoclinic orbit. Note that it is impossible to follow a homoclinic orbit numerically. We discuss this in more detail in Section 5.5.. 5.4.2 Asymmetric Perturbations In [21], a perturbation analysis is performed on initial conditions of the form. ψ (x, 0) = ψ˜ + eiθ1 cos (µn x) + rn eiθ2 sin (µn x) .. (5.40). It is shown that the imaginary double points split as follows:. ζ1 =.  a2 [sin (θ1 + φ) sin (φ − θ1 )  4ζ02         + rn2 sin (θ2 + φ) sin (φ − θ2 )    + irn sin (θ2 − θ1 ) sin (2φ)]       0. , if. k=. µn 2. if. k 6=. µn 2. (5.41).

(73) CHAPTER 5.. where tanφ. 41. THE PERIODIC NLS. = Im (ζ0 ) /kn ,. and. ζ = ζ0 + ζ1 + · · · .. Figure 5.5 shows the results obtained from our numerical experiments. The numerical solution was obtained by using the second order split-step scheme (3.41). The spectrum was calculated by Algorithm 1. Once again the results match the analytical results. Note that the spectrum in Figure 5.4 and 5.5 were calculated at. t = 0.. This. shows that we can gain important information about the solution of a perturbed initial condition by looking at its spectrum.. 1. y. 0.5. 0. −0.5. −1. −2. −1.5. −1. −0.5. 0. 0.5. 1. 1.5. 2. x Figure 5.3:. Phase space for the two dimensional homoclinic system. (5.39).

(74) CHAPTER 5.. 42. THE PERIODIC NLS. 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35 −2. −1.5. −1. −0.5. 0. Re ζ. 0.5. 1. 1.5. 2 −3. x 10. (a) 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35 −2. −1.5. −1. −0.5. 0. Re ζ. 0.5. 1. 1.5. 2 −3. x 10. (b) 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35. −3. −2. −1. 0. Re ζ. 1. 2. 3 −3. x 10. (c) Evolution and spectrum for symmetric perturbations of the form (5.23) . (a) π is inside the homoclinic orbit (φ = 0) , (b) is outside the homoclinic orbit φ = 2 , and (c) π is near the homoclinic orbit φ = 4 . Figure 5.4:.

(75) CHAPTER 5.. 43. THE PERIODIC NLS. 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35 −2. −1.5. −1. −0.5. 0. Re ζ. 0.5. 1. 1.5. 2 −3. x 10. (a). 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35 −2. −1.5. −1. −0.5. 0. Re ζ. 0.5. 1. 1.5. 2 −3. x 10. (b). 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35 −2. −1.5. −1. −0.5. 0. Re ζ. 0.5. 1. 1.5. 2 −3. x 10. (c). 0.355. Im ζ. 0.354. 0.353. 0.352. 0.351. 0.35 −2. −1.5. −1. −0.5. 0. Re ζ. 0.5. 1. 1.5. 2 −3. x 10. (d). Numerical solutions for initial condition (5.40). In (a), θ1 = 0 and θ2 = − π6 . In (b), θ1 = 0 and θ2 = π6 . In (c), θ1 = π6 and θ2 = π2 . In (d), θ1 = 5π6 and θ2 = π2 . Figure 5.5:.

(76) CHAPTER 5.. 5.5. 44. THE PERIODIC NLS. Evolution of the Discrete Spectrum. One of the properties of the discrete spectrum is that it is time invariant. This means that the evolution of the potential does not change the discrete spectrum. For the numerical calculation of the spectrum, the situation is dierent. In this section we consider the inuence of time on the numerical results of the spectrum.. 5.5.1 Symmetric and Asymmetric Perturbations Figure 5.6 (a) shows the evolution of the spectrum for a perturbation inside the homoclinic orbit, corresponding to Figure 5.4 (a). remains constant until. t ≈ 11.. It is clear that the spectrum. This is followed by an increase in the imaginary. direction, the same direction as the split of the double point at. t = 0.. Notice how. the spectrum settles back to its former position. In Figure 5.6 (b) we see the evolution of the eigenvalues for a perturbation outside the homoclinic orbit, corresponding to Figure 5.4 (b). the real direction. Once again it remains there until increase in the imaginary part of. |ζ − ζ0 |,. Here, our initial split is in. t ≈ 11.. For this case we see an. while the real part of. |ζ − ζ0 | goes to zero.. Once the hump fades, the spectrum returns to its former position. This is similar to the perturbation inside the homoclinic orbit.. −4. x 10 15. 0. Real. 10. Real. Imaginary Imaginary. −1. 10. −2. 10. | ζ−ζ0 |. | ζ−ζ0 |. 10. −3. 10. −4. 10 5. −5. 10. −6. 10 0 0. 5. 10. 15. 20. 25. Time (a). 30. 35. 40. 45. 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. 50. Time (b). Evolution of discrete spectrum |ζ − ζ0 |. (a) corresponds to a potential inside the homoclinic orbit. (b) corresponds to a potential outside the homoclinic spectrum. Figure 5.6:.

No results found