Pseudospectral Methods in Quantum and Statistical Mechanics by Joseph Quin Wai Lo B.Sc., Simon Fraser University, 2002 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Mathematics) The University of British Columbia (Vancouver, Canada) July, 2008 c Joseph Quin Wai Lo 2008 Abstract The pseudospectral method is a family of numerical methods for the solution of differential equations based on the expansion of basis functions defined on a set of grid points. In this thesis, the relationship between the distribution of grid points and the accuracy and convergence of the solution is emphasized. The polynomial and sinc pseudospectral methods are extensively studied along with many applications to quantum and statistical mechanics involving the Fokker-Planck and Schr¨ odinger equations. The grid points used in the polynomial methods coincide with the points of quadrature, which are defined by a set of polynomials orthogonal with respect to a weight function. The choice of the weight function plays an important role in the convergence of the solution. It is observed that rapid convergence is usually achieved when the weight function is chosen to be the square of the ground-state eigenfunction of the problem. The sinc method usually provides a slow convergence as the grid points are uniformly distributed regardless of the behaviour of the solution. For both polynomial and sinc methods, the convergence rate can be improved by redistributing the grid points to more appropriate positions through a transformation of coordinates. The transformation method discussed in this thesis preserves the orthogonality of the basis functions and provides simple expressions for the construction of discretized matrix operators. The convergence rate can be improved by several times in the evaluation of loosely bound eigenstates with an exponential or hyperbolic sine transformation. The transformation can be defined explicitly or implicitly. An explicit transformation is based on a predefined mapping function, while an implicit transformation is constructed by an appropriate set of grid points determined by the behaviour of the solution. The methodologies of these transformations are discussed with some applications to 1D and 2D problems. The implicit transformation is also used as a moving mesh method for ii Abstract the time-dependent Smoluchowski equation when a function with localized behaviour is used as the initial condition. iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Abstract List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . x Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . 1.1 Preliminaries . . . . . . . . . . . . . . . 1.2 The Fokker-Planck equation . . . . . . 1.3 The Schr¨ odinger equation . . . . . . . . 1.4 Overview of current numerical methods 1.5 Interpolation and quadrature . . . . . . 1.5.1 Orthogonal polynomials . . . . . 1.5.2 Sinc functions . . . . . . . . . . 1.5.3 Transformation of grid points . 1.6 Pseudospectral method . . . . . . . . . 1.7 Thesis objective . . . . . . . . . . . . . . . . . . . . . . . . 1 1 3 6 7 9 10 13 14 15 17 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2 Spectral convergence of quadrature discretization method 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Discretization of the Fokker-Planck and Schr¨ odinger eigenvalue problems . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Electron relaxation . . . . . . . . . . . . . . . . . . . . . . . 2.4 Vibrational states of I2 . . . . . . . . . . . . . . . . . . . . . 2.5 The quartic bistable system . . . . . . . . . . . . . . . . . . . 28 28 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 35 36 37 iv Table of Contents 2.6 2.7 2.8 Convergence of eigenvalues; discussion of results . . . . . . . Interpretation of the numerical results . . . . . . . . . . . . . Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 52 57 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3 Pseudospectral methods of solution of the Schr¨ odinger equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Spectral convergence . . . . . . . . . . . . . . . . . . . . . . 3.3 Pseudospectral method . . . . . . . . . . . . . . . . . . . . . 3.4 The vibrational states of the Morse oscillator . . . . . . . . . 3.5 The vibration states of Ar2 . . . . . . . . . . . . . . . . . . . 3.6 Woods-Saxon potential . . . . . . . . . . . . . . . . . . . . . 3.7 Results and discussions . . . . . . . . . . . . . . . . . . . . . 3.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 65 67 69 71 73 74 74 80 82 4 An efficient mapped pseudospectral method for weakly bound states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.2 Solution of the Schr¨ odinger equation on a grid of quadrature points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.2.1 The coefficient space representation . . . . . . . . . . 89 4.2.2 The physical space representation . . . . . . . . . . . 92 4.3 The solution of the Schr¨ odinger equation with the sinc functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.4 Explicit transformations of grid points . . . . . . . . . . . . . 95 4.5 Weakly bound vibrational states of He2 , Ne2 , and Ar2 . . . . 97 4.5.1 Tang-Toennies potential . . . . . . . . . . . . . . . . 98 4.5.2 Tang-Toennies-Yiu and Aziz-Slaman He2 potentials . 102 4.5.3 Justifications for the computational domains . . . . . 104 4.6 A benchmark calculation for Cs2 ; elimination of ghost levels 105 4.7 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 109 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5 Conclusions and future directions . . . . . . . . . . . . . . . 125 5.1 Preliminary investigations in mappings and other directions 125 5.1.1 Interpolation with mapping . . . . . . . . . . . . . . . 126 v Table of Contents 5.1.2 5.1.3 5.2 5.3 Interpolatory moving mesh . . . Eigenvalues of the Smoluchowski tion . . . . . . . . . . . . . . . . Summary of results . . . . . . . . . . . Future directions . . . . . . . . . . . . . . . . . . . . . . . . . and Kramers equa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 138 145 150 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Appendices A Lorentz Fokker-Planck Operator . . . . . . . . . . . . . . . . 154 B Fokker-Planck Operator for the Bistable System . . . . . . 157 C Matrix Representation of the Schr¨ odinger Operator . . . . 161 D Moving mesh method for the time dependent Smoluchowski equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 vi List of Tables 2.1 2.2 2.3 Convergence of the eigenvalues for the quartic potential, ǫ = 0.1. 45 Convergence of the eigenvalues for the quartic potential, ǫ = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Convergence of the eigenvalues for the quartic potential, ǫ = 0.001. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1 Convergence of 103 E10 for the I2 Morse potential. . . . . . . . 4.1 The parameters of the Tang-Toennies potentials for He2 , Ne2 , and Ar2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 The parameters of the Aziz-Slaman LM2M2 potential for He2 . 103 Convergence of the Morse potential for Cs2 . . . . . . . . . . . 106 The vibrational energy levels for He2 , Ne2 and Ar2 . . . . . . . 110 Convergence of the energy for Aziz-Slaman He2 potential. . . 117 4.2 4.3 4.4 4.5 5.1 5.2 5.3 75 Eigenvalues of Smoluchowski operator evaluated using different mappings . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 Eigenvalues of Smoluchowski operator evaluated using different mappings . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 The number of grid points required to obtain the eigenvalues 146 vii List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 3.1 3.2 3.3 3.4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 The relative error of eigenvalues versus the number of quadrature points for the electron thermalization problem. . . . . . The relative error of eigenvalues versus the number of quadrature points for the I2 Morse potential. . . . . . . . . . . . . . Variation of the potential in the Schr¨ odinger equation for the quartic bistable system. . . . . . . . . . . . . . . . . . . . . . The relative error of eigenvalues versus the number of quadrature points for the bistable system with ǫ = 0.01. . . . . . . . The relative error of eigenvalues versus the number of quadrature points for the bistable system with ǫ = 0.001. . . . . . . Variation of the expansion coefficients of the eigenfunctions for the electron relaxation problem. . . . . . . . . . . . . . . . The expansion coefficients of the eigenfunctions of the I2 system. The relative error of eigenvalues versus the number of quadrature points for the I2 Morse potential. . . . . . . . . . . . . . The relative error of eigenvalues versus the number of quadrature points for the Ar2 Cahill and Parsegian potential. . . . . The relative error of eigenvalues versus the number of quadrature points for the Woods-Saxon potential. . . . . . . . . . . The mapping of collocation points from Laguerre grid to a modified grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . The He2 , Ne2 , and Ar2 potentials. . . . . . . . . . . . . . . . The zero energy limit of the l = 0 phase shift for He2 , Ne2 , and Ar2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The eigenvalues and the ghost level versus the number of collocation points for Cs2 Morse potential. . . . . . . . . . . . The expansion coefficients and the eigenfunction of the ghost level for Cs2 Morse potential. . . . . . . . . . . . . . . . . . . The ground state eigenfunction for Ne2 . . . . . . . . . . . . . The highest state eigenfunction of Ne2 . . . . . . . . . . . . . . 40 42 43 48 49 54 68 77 79 81 96 99 100 108 108 110 111 viii List of Figures 4.8 4.9 4.10 4.11 4.12 4.13 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 The relative error of eigenvalues versus the number of collocation points for Ne2 . . . . . . . . . . . . . . . . . . . . . . . Variation of the expansion coefficients of the highest state eigenfunction of Ne2 . . . . . . . . . . . . . . . . . . . . . . . . The relative error of eigenvalues versus the number of collocation points for Ar2 . . . . . . . . . . . . . . . . . . . . . . . . The relative error of the eigenvalue versus the number of collocation points for He2 . . . . . . . . . . . . . . . . . . . . . . The eigenfunction of the Aziz-Slaman He2 potential. . . . . . The eigenvalue and the ghost level versus the number of collocation points for Aziz-Slaman He2 potential. . . . . . . . . . Unmapped sinc interpolation . . . . . . . . . . . . . . . . . . Sinc interpolation with hyperbolic sine mapping . . . . . . . . The reference grid and the corresponding mapping function . Sinc interpolation with implicitly defined mapping . . . . . . Solution of the Smoluchowski equation by interpolatory moving mesh method . . . . . . . . . . . . . . . . . . . . . . . . . The evolution of grid points over time for the solution of Smoluchowski equation . . . . . . . . . . . . . . . . . . . . . . The error of the solution of Smoluchowski equation for interpolatory moving mesh method and unmapping sinc method . Eigenfunctions of Smoluchowski equation solved by unmapped sinc method . . . . . . . . . . . . . . . . . . . . . . . . . . . . The maps constructed by rough sketches of eigenfunctions . . 112 114 115 118 119 121 127 129 131 133 136 137 137 140 141 ix Acknowledgements I would like to express my gratitude to Professor Bernie Shizgal. His expert guidance and extreme patience have been invaluable throughout my graduate study and research. I appreciate my parents for taking care of me all the time. I am also grateful to my beloved girlfriend, Joanne, for her endless support, encouragement, and patience. x Statement of Co-Authorship The studies of QDM and sinc methods in Chapter 2 were suggested by Professor Shizgal. The transformation methods in Chapters 3 – 5 are my own work. Applications to Fokker-Planck and Schr¨ odinger equations were decided together by Professor Shizgal and me. All numerical simulations and results in this thesis were done by me. I wrote the first draft of each chapter, and subsequent drafts were done together with Professor Shizgal. xi Chapter 1 Introduction 1.1 Preliminaries The prediction and analysis of physics phenomena often require the numerical solution of differential equations subject to well defined boundary conditions [28, 94]. The development of numerical methods for the solution of a large variety of problems is of practical importance. Before specifying the systems of interest in this thesis, it is useful to provide a brief summary. One of the simplest partial differential equations is the heat equation [94]. By letting x be the position on a rod, t be time, u(x, t) be the temperature at x, and the constant k be the diffusivity, the one-dimensional (1D) heat equation ∂u 1 ∂2u (1.1) = ∂t k ∂x2 models the temperature of the rod at any time t if an initial condition u(x, 0) and appropriate boundary conditions are provided. Different types of boundary conditions represent different situations. For example, if the ends of the rod are connected to heat baths of constant temperatures, the temperatures at the end points are assumed to be constant for all t. If the ends of the rod are insulated, there will be no heat flow at the ends of the rod, and the boundary condition becomes ∂u/∂x = 0 at the endpoints. For the diffusion of heat on a uniform plate or solid, the temperature u(x, t) at the position x = [x, y] in two dimensions (2D) or [x, y, z] in three dimensions (3D) is governed by the heat equation ∂u 1 = ∇2 u, ∂t k where ∇2 = ∂ 2 /∂x2 + ∂ 2 /∂y 2 in 2D or ∂ 2 /∂x2 + ∂ 2 /∂y 2 + ∂/∂z 2 in 3D is the Laplacian. The heat equation can be generalized as a diffusion equation which models the diffusion of a substance in a media when u is considered as the concentration of the substance at x [82,100,108]. If the rate of diffusion 1 Chapter 1. Introduction is not a constant but depends on the position or the concentration, the diffusion equation is given by ∂u = ∇ · (D(x, u)∇u), ∂y (1.2) where ∇ = [∂/∂x, ∂/∂y] in 2D or [∂/∂x, ∂/∂y, ∂/∂z] in 3D is the gradient operator, and D(x, u) is the diffusion coefficient. The study of fluid motion in porous media [30, 38, 51, 85], tumor development [40, 49, 75], and image noise reduction, [54, 67, 95, 107] are some of the major applications to Eq. (1.2). A slight variation to the heat equation, Eq. (1.1), is the reaction-diffusion equation ∂2u ∂u = D 2 + f (u), ∂t ∂x where D is the diffusion coefficient, and f (u) is the reaction term [44,74,75]. As a simple example for the reaction-diffusion equation, Fisher’s equation [68, 77, 101], where f (u) = au(1 − bu), models the growth of the population of a species with a net reproduction rate of a and a carrying capacity of 1/b. The reaction term f (u) = au(1 − bu) represents the logistic growth of the population. Fisher’s equation is widely used as a demonstration of traveling waves occurring in combustion [25, 36, 86]. In many multi-component systems, the growth of one component is influenced by other components in the system. The resulting reaction-diffusion equation for each of the component will contain the reaction term, or interaction term, f , depending on all components. A two-component, u1 (x, t) and u2 (x, t), reaction-diffusion system has the form ∂u1 = D1 ∇2 u1 + f1 (u1 , u2 ) ∂t ∂u2 = D2 ∇2 u2 + f2 (u1 , u2 ) ∂t where x is any-dimensional position. One of the examples directly related to Fisher’s equation is the Lotka-Volterra system with diffusion, which is used for predicting the population of two competing species [32,60,75]. The interaction terms of the Lotka-Volterra system, f1 (u1 , u2 ) = a1 u1 (1−b1 u1 −c1 u2 ) and f2 (u1 , u2 ) = a2 u2 (1 − b2 u2 − c2 u1 ), model the logistic growth with competition. The parameters ai , 1/bi , and ci for i = 1, 2 represent the net birth 2 Chapter 1. Introduction rates, the carrying capacities, and the abilities for competing with the other species, respectively. Another example related to the diffusion equation is Burgers’ equation ∂u ∂2u ∂u +u = ν 2, ∂t ∂x ∂x (1.3) which is used for the simulation of 1D fluid flow with flow velocity u(x, t) [18, 28, 34, 101, 106]. The first term, ∂u/∂t, represents the change of the flow with respect to t. For a steady flow, ∂u/∂t = 0. The second term, u(∂u/∂x), is the advection term, which represents the change of u along the flow. The right hand side, ν(∂ 2 u/∂x2 ) is the viscosity term, with ν as the viscosity coefficient. In case of ν = 0, Eq. (1.3) is used as a simple model for the development of shock waves. Burgers’ equation can be converted to a heat equation ∂φ/∂t = ν∂ 2 φ/∂x2 by the Cole-Hopf transformation u = −2ν(∂φ/∂x)/φ [28]. In this thesis, we will focus on the numerical methods for the solutions to the Fokker-Planck equation and the Schr¨ odinger equation. Since the Fokker-Planck equation, the Schr¨ odinger equation, and the differential equations discussed above are all parabolic partial differential equations, the discussions of the numerical methods in this thesis will be useful for the understanding of physical phenomena described by the above differential equations. 1.2 The Fokker-Planck equation The motion of particles influenced by an external fluctuating force is referred to as Brownian motion. The position and the velocity of these particles can be described by a probability distribution. The Fokker-Planck equation was developed to study the evolution of Brownian motion in terms of a probability density function [81, 99]. Kramers explained the kinetics of chemical reactions using Brownian motion in a potential as a model, and derived the Kramers equation which is a 2D Fokker-Planck equation [45,58,81]. H¨ anggi et al [47] and Pollak et al [79] have published reviews on the history, the analysis, and the development of Kramers’ reaction theory. 3 Chapter 1. Introduction Consider the 1D Brownian motion of a particle in a background fluid. The equation of motion of the particle is given by d2 x dx = F (x) + Γ(t), +ν 2 dt dt (1.4) where ν is the damping constant, F (x) = −dU (x)/dx is the force per unit mass due to an external potential U (x), and Γ(t) is a randomly fluctuating force per unit mass. Equation (1.4) is called the Langevin equation [81]. The stochastic force Γ(t) is assumed to be Gaussian distributed with zero mean since the average motion described by Eq. (1.4) can be shown to be independent of Γ(t). It is also assumed that the collisions of different fluid particles are independent. In other words, the correlation of the fluctuating force at different times t and t′ is zero, i.e. Γ(t)Γ(t′ ) = 2νǫδ(t − t′ ), where δt is the unit impulse function. The parameter ǫ controls the strength of the fluctuation. The Langevin equation, Eq. (1.4), can be solved statistically by simulating an ensemble of particles when it is written in the form dx = v, dt dv = −νv + F (x) + Γ(t), dt (1.5a) (1.5b) and Γ(t) is generated by a random number generator. Since the motion is influenced by a stochastic force, the motion itself is also stochastic. The probability density function describing the position and the velocity of the ensemble of particles, P (x, v, t), satisfies Kramers equation ∂P ∂P ∂P ∂(vP ) ∂2P = −v − F (x) +ν + νǫ 2 . (1.6) ∂t ∂x ∂v ∂v ∂v Originally, Kramers developed Eq. (1.6) for reaction kinetics, and the position x is regarded as the reaction coordinate [58]. The potential U (x) is often modeled with a bistable function, for example U (x) = x4 /4 − x2 /2, with one well corresponding to the reactants and another to the products [6, 12, 13, 23, 33, 41, 43, 91]. An example is given by Blackmore and Shizgal for the trans-gauche isomerization of n-butane [14]. For large damping, ν → ∞, the distribution in v is close to Maxwellian, and the second derivative in Eq. (1.4) can be neglected. The equation of 4 Chapter 1. Introduction motion in this limit becomes dx 1 = [F (x) + Γ(t)]. dt ν ∞ −∞ P (x, v, t) dv The probability density function P (x, t) = Smoluchowski equation [81] ∂P 1 ∂(F P ) ∂2P = − +ǫ 2 . ∂t ν ∂x ∂x satisfies the (1.7) A general 1D Fokker-Planck equation is given by [81] ∂ ∂2 ∂ P (x, t) = [A(x)P (x, t)] + 2 [B(x)P (x, t)]. ∂t ∂x ∂x (1.8) The functions A(x) and B(x) are the drift and diffusion coefficients, respectively. The Smoluchowski equation in Eq. (1.7) is an example of a 1D Fokker-Planck equation. Other uses of the Fokker-Planck equation include the studies of laser and plasma physics [31, 69, 78, 81], and climate systems [35, 76]. The equilibrium density function P0 (x) of Eq. (1.8) is P0 (x) = Kx exp − B(x) x A(x′ ) ′ dx . B(x′ ) with a normalization factor Kx . By setting P (x, t) = (1.8) will have the self-adjoint form 1 ∂ ∂p ∂ =√ BP0 ∂t ∂x P0 ∂x p √ P0 = LFP p. (1.9) P0 (x)p(x, t), Eq. (1.10) A classical approach to the solution of a linear parabolic partial differential equation such as the Fokker-Planck equation in Eq. (1.10) is to apply a separation of variables. In terms of the eigenvalues, λm , and the eigenfunctions, ψm (x), of the Fokker-Planck operator LFP , the solution p(x, t) is given λm t ψ (x). It can be verified that LFP is a self-adjoint by p(x, t) = ∞ m m=0 an e negative semi-definite operator so that · · · < λ2 < λ1 < λ0 = 0. Furthermore, the ground-state eigenfunction is given by ψ0 (x) = P0 (x) and all eigenfunctions are orthogonal. Therefore, the solution p can be written in the form ∞ p(x, t) = an e−|λm |t ψm (x) P0 (x) + m=1 5 Chapter 1. Introduction and converges to P0 (x) as t → ∞. The eigenvalues then provide the rate of the evolution of the solution, and hence the study of the eigenvalue problem LFP ψ = λψ is important. It is worth noting that a Fokker-Planck equation can always be transformed to another Fokker-Planck equation with a constant diffusion coefficient with a suitable variable transformation [81]. If we transform x to x a new variable y = [B(x′ )]−1/2 dx′ and define the new solutions to be 1/4 p˜(y, t) = [B(x(y))] p(x(y), t), Eq. (1.10) becomes 1 ∂ √ ∂ ∂ p˜ = √ BP0 1/2 ∂t ∂y ∂y [ BP0 ] p˜ √ [ BP0 ]1/2 ˜ FP p˜. =L (1.11) By comparing Eq. (1.11) with Eq. (1.10), the new diffusion coefficient is ˜ B(y) = 1, and the new equilibrium solution is P˜0 (y) = B(x(y))P0 (x(y)). ˜ It follows that the new drift coefficient is A(y) = P˜0′ (y)/P˜0 (y), analogous to Eq. (1.9). 1.3 The Schr¨ odinger equation The time independent Schr¨ odinger equation 2 Hψ = − d2 ψ + V (y)ψ = Eψ, 2µ dy 2 (1.12) where is Planck’s constant, µ is the mass of a particle, V (y) is the potential function, and E is the energy of the state of a system, serves as the basis of quantum mechanics [26, 80, 98]. The solution ψ(y), called the wavefunction, describes the state of the system. The Schr¨ odinger operator H is usually referred to as the Hamiltonian. The multi-dimensional Schr¨ odinger equation 2 Hψ = −∇ ψ + V (x)ψ is used in the studies of atoms or molecules. One of the applications for the Schr¨ odinger equation, Eq. (1.12), is to calculate the vibrational energy levels of diatomic molecules [3]. The Morse potential V (y) = De [1 − exp(−αy)]2 is a simple model for approximating the molecular potential energy when the nuclei of the atoms are displaced away from their equilibrium position by length y [17, 73, 88]. Since the exact energy levels for the Morse potential are known, it serves as both a simple estimation of true energy levels and a benchmark for numerical techniques [7, 10, 11, 15, 20, 50, 62, 65, 70, 72, 84, 88, 96, 102, 103, 105]. In Chapters 6 Chapter 1. Introduction 3 and 4, we will use the Morse potential to construct the basis functions for the calculation of vibrational energy levels when a real potential is used. The problem can be extended to the studies of vibrations in triatomic or polyatomic molecules [22, 27, 59, 83]. The time independent Fokker-Planck equation, Eq. (1.11), is given by ˜ FP ψ˜m = L 1 d ˜ d P0 dy P˜0 dy ψ˜m P˜0 ˜ m ψ˜m , =λ (1.13) ˜ m is the mth eigenvalues of L ˜ FP , and ψ˜m (y) is the corresponding where λ eigenfunction. Eq. (1.13) can be directly simplified to the form 2 2ψ ˜m ˜ ′′ 1 P˜ ′ P d 1 0 ˜ m ψ˜m , ˜ FP ψ˜m = ψ˜m = λ L − 0 − (1.14) dy 2 2 P˜0 4 P˜0 where ′ denotes d/dy. Equation (1.14) can be regarded as a Schr¨ odinger equation. Therefore, the Schr¨ odinger equation is also used as an alternative method for solving the eigenvalues of Fokker-Planck equations in Chapter 2. If = 1, µ = 1/2, and V (y) = P˜0′′ (y)/2P˜0 (y) − (1/4)[P˜0′ (y)/P˜0 (y)]2 , the ˜ FP by comparing Eqs. (1.14) and (1.12). If Hamiltonian is exactly H = −L ˜ m = −Em . we let λm be the eigenvalues of LFP , then λm = λ 1.4 Overview of current numerical methods Since the exact solutions to most differential equations cannot be represented by simple closed forms, numerical methods are required for obtaining approximate solutions. For the solutions of eigenvalue problems such as the Schr¨ odinger equation, Eq. (1.12), there are two major types of numerical methods: the shooting method, and the basis function expansion. For the shooting method, the solution of a boundary value problem is found by evaluating two initial value problems [2,16,48], which can be solved by many integration schemes such as the Euler, Runge-Kutta or Numerov methods [46]. However, since the eigenvalues are unknown, the shooting method finds the eigenvalues by refining from initial guesses one at a time. One way to obtain the initial guesses is to set up a matrix eigenvalue problem by applying the finite difference method to the boundary value problem 7 Chapter 1. Introduction and solve for the eigenvalues. One of the eigenvalues will then be put back to the boundary value problem. The solution is then integrated from both end points using the Runge-Kutta, Numerov, or multistep methods to a matching point. The difference from the two parts of the solution at the matching points will be used as a parameter for the calculation of the eigenvalue correction. This procedure is repeated until the correction is within the tolerance level [21, 50, 52, 55, 56]. The shooting method is more widely used for quantum mechanics problems when the Schr¨ odinger equation is represented by the Hamiltonian canonical form [53, 64, 72, 84]. This method is referred to as the symplectic shooting method. One of the advantages for the shooting method is that the accuracy of eigenvalues can be improved with a relatively low increase in computational cost if explicit integration schemes can be used, which makes it a good method for refinements of eigenvalues. However, the shooting method can find only one eigenvalue at a time. Also, there may be a problem for numerical stability when an explicit integration scheme is used. Another method for eigenvalue problems is based on basis function expansions. The differential equation is then represented (approximated) by a matrix problem and the eigenvalues are approximated by the eigenvalues of the matrix. Fourier, sinc, and polynomial bases are some of the popular bases used in eigenvalue problems [1, 8, 24, 61, 70, 71, 92, 97, 102]. These methods are in general referred to as spectral methods. If the basis functions are defined on a grid of points, then these methods are referred to as collocation-spectral methods, or pseudospectral methods [19]. With pseudospectral methods, eigenvalues can be found by a single matrix diagonlization without the need of correction schemes. Also, the errors of the eigenvalues decrease exponentially if an appropriately chosen basis is used. However, the computational cost increases dramatically as the number of grid points or the number of basis functions increases. The main objective for developing an efficient pseudospectral method is to find a set of basis function so that accurate eigenvalues can be obtained with as few grid points as possible. A few pseudospectral schemes have been developed independently. Light and coworkers introduced the discrete variable representation (DVR) [62,63], and Baye and coworkers introduced the Lagrange mesh method [9], both for the evaluation of quantum mechanical problems. These methods are usually based on classical polynomials, but can be applied to sinc, Fourier, and other bases as well. Shizgal et al developed the quadrature discretization method 8 Chapter 1. Introduction (QDM) [14, 89, 90] for the solution of Boltzmann equation, and applied the same method to Sturm-Liouville and Fokker-Planck equations using nonclassical orthogonal polynomials. The implementation of all these methods are equivalent, and can be generalized as pseudospectral methods with orthogonal interpolating basis function. For the calculation of vibrational energy eigenstates of diatomic molecules, polynomial methods with appropriately chosen weight functions are usually the most efficient for eigenvalues well below the dissociation energy, as shown in Chapter 4 and in Ref. 92. However, for eigenvalues close to the dissociation energy, the corresponding eigenfunctions are very loosely bound, and polynomial methods become less efficient. The efficiency of existing basis functions can sometimes be improved with transformation of coordinates. By doing so, the grid points can be relocated to the positions where the solution or the eigenfunctions can be better described. The transformation method has been widely used in sinc bases [57,66,93,105] but seldom mentioned in polynomial methods. In Chapter 4, we will apply transformation methods to both the polynomial and the sinc methods for the evaluation of eigenvalues in loosely bound states. Methods for determining suitable transformations will be developed in Chapter 5. In the following sections, we will discuss the basics of the pseudospectral method used in this thesis. 1.5 Interpolation and quadrature Interpolation is a method for estimating unknown values which lie within the range of a discrete set of given data points [4, 5, 29]. Precisely, for an unknown function g(x) with given data points {(xi , g(xi ))}N i=1 , xi < xj for i < j, one can estimate the value g(x) for x ∈ (x1 , xN ) by interpolation. A function, f (x), which interpolates g(x) using N given points is usually written as an expansion of N basis functions {φn (x)}N n=1 in the form N f (x) = cn φn (x). n=1 The expansion coefficients cn are solved by setting f (xi ) = g(xi ). The function f (x) is called the interpolant of g(x). 9 Chapter 1. Introduction It is usually convenient to introduce a new set of basis functions Cj (x) where each interpolates the data points {(xi , δij )}N i=1 , and δij is the Kronecker delta, and has the value 1 if i = j and 0 otherwise. The functions Cj (x) are the interpolating functions, or the cardinal functions, and satisfy Cj (xi ) = δij , referred to as the cardinality condition. The interpolation of g(x) with given {(xi , g(xi ))}N i=1 is represented by N f (x) = g(xj )Cj (x). j=1 We focus on two types of interpolations based on orthogonal polynomials and sinc functions. 1.5.1 Orthogonal polynomials Polynomials are often used for interpolation. If we set φn (u) = un−1 , the polynomial interpolating functions, called the Lagrange interpolating function, are given by N u − ui Ij (u) = uj − ui i=1 i=j which are polynomials of degree N − 1, denoted by PN −1 [5]. It is thus clear that if g(u) ∈ Pk for k ≤ N − 1, the interpolation N f (u) = g(uj )Ij (u) (1.15) j=1 is exact, and the integral N b w(u)g(u) du = a wj g(uj ) (1.16) j=1 b is exact for some function w(u) defined in an interval [a, b] if wj = a w(u)Ij (u) du. This method of integral evaluation is a quadrature rule, and wj are the quadrature weights. 10 Chapter 1. Introduction The Lagrange interpolation given by Eq. (1.15) can be defined for any choice of discrete points ui . Therefore it is possible to consider a particular N set {uj }N j=1 , and the corresponding set of weights {wj }j=1 , which satisfies b k a w(u)u k du = N j=1 wj uj for k = 0 . . . 2N −1. With uj and wj , the quadrature rule in Eq. (1.16) gives an exact integral value for a larger family of g(u) ∈ Pk , k ≤ 2N − 1. The actual values of uj and wj will be given later. For w(u) ≥ 0, the polynomials Qn (u) generated by the three-term recurrence relation Q0 (u) = 1, Q1 (u) = u − α0 , (1.17a) Qn+1 (u) = (u − αn )Qn (u) − βn Qn−1 (u), (1.17b) where the coefficients αn and βn are given by αn = 1 γn b uw(u)[Qn (u)]2 du, βn = a b γn = w(u)[Qn (u)]2 du, γn , γn−1 (1.17c) (1.17d) a b satisfy the orthogonality a w(u)Qm (u)Qn (u) du = γn δmn [37, 39]. With Eqs. (1.17a) and (1.17b), the polynomial Qn (u) can be represented in determinant form √ u − α0 − β1 0 ··· 0 √ √ .. .. − β1 u − α1 − β2 . . √ .. Qn+1 (u) = det (1.18) . 0 − β2 u − α2 0 √ .. .. .. .. . . . . − βn √ 0 ··· 0 − βn u − αn Since we assume that the quadrature rule in Eq. (1.16) is exact for g(u) ∈ P2N −1 with the appropriate choice of uj and wj , the orthogonality for the √ normalized polynomials Pn (u) = Qn (u)/ γn becomes N b w(u)Pm−1 (u)Pn−1 (u) du = a wj Pm−1 (uj )Pn−1 (uj ) = δmn j=1 √ for 1 ≤ m, n ≤ N . If Tjn = wj Pn−1 (uj ), then Tt T = I = TTt , where the superscript t denotes the matrix transpose. Hence, T is an orthogonal mab trix. By defining the Jacobi matrix Jmn = a w(u)uPm−1 (u)Pn−1 (u) du = 11 Chapter 1. Introduction N k=1 wj uj Pm−1 (uj )Pn−1 (uj ) for 1 ≤ m, n ≤ N and the diagonal matrix Uij = uj δij , the matrix J will satisfy J = Tt UT, and it is clear that uj are the eigenvalues of J. The corresponding eigenvectors vj are the column vectors of Tt . It is easy to show that √ α0 β1 0 √ √ β1 α1 β2 √ J= 0 β2 α2 .. .. .. . . . 0 ··· 0 ··· .. . .. . .. . βN −1 0 .. . 0 βN −1 αN −1 . Hence, with Eq. (1.18), uj are the roots of both QN (u) and PN (u). The nth √ component of the jth eigenvector of J has the form [vj ]n = wj Pn−1 (uj ). −1/2 Since P0 (uj ) = γ0 , the quadrature weights are given by wj = γ0 ([vj ]1 )2 . The orthogonality of T gives the equality N n=1 √ δij . With wi wj replaced by wj , the polynomials √ wi wj Pn−1 (ui )Pn−1 (uj ) = N Ij (u) = wj Pn−1 (u)Pn−1 (uj ) n=1 of degree N −1 satisfy both the cardinality Ij (ui ) = δij and the orthogonality b a w(u)Ii (u)Ij (u) du = wj δij . If we further define the interpolating functions as w(u) CjP (u) = Ij (u), w(uj ) CjP (u) will be orthogonal with respect to the unit weight function, i.e. b P P a Ci (u)Cj (u) du = ηj δij , where ηj = wj /w(uj ). The constants ηj are the weights of the quadrature rule N b a g(u) du ≈ ηj g(uj ). j=1 12 Chapter 1. Introduction 1.5.2 Sinc functions Another popular type of interpolation based on the sinc function has been well studied by Whittaker [104] and Shannon [87]. The sinc function is defined by [66, 93] sin πu if u = 0, sinc(u) = πu 1 if u = 0. The integral form of the sinc function with translation j, sinc(u) = 1 2π π exp(−iωu) dω, −π indicates that the Fourier transform of sinc(u) is the characteristic function χ[−π,π](ω) = 1 0 if −π ≤ ω ≤ π, otherwise, which implies that ∞ −∞ sinc(u) exp(iωu) du = χ[−π,π] (ω). (1.19) The orthogonality of translated sinc functions is derived by Parseval’s theorem ∞ −∞ sinc(u) sinc(u − j) du = 1 2π π exp(−ijω) dω = δ0j (1.20) −π if j is an integer. With these formulas, we can define the sinc interpolating functions Sj (u) = sinc u − uj h along the uniform grid of N points uj = umin +(j −1)h when the width of the grid is h = (umax − umin )/(N − 1). The interval [umin , umax ] is assumed to be the interval of interpolation. These Sj (u) satisfy the cardinality Sj (ui ) = δij ∞ by the definition of sinc function, and the orthogonality −∞ Si (u)Sj (u) du = hδij by Eq. (1.20). By using the interpolation g(u) ≈ N j=1 g(uj )Sj (u), the integral can be approximated by the quadrature rule N ∞ −∞ g(u) du ≈ hg(uj ) j=1 13 Chapter 1. Introduction since 1.5.3 ∞ −∞ Sj (u) du = h which can be derived by Eq. (1.19) for ω = 0. Transformation of grid points The polynomial interpolation CjP (u) in Section 1.5.1 and the sinc interpolation Sj (u) in Section 1.5.2 are based on predefined grids uj . However, this interpolation will be inefficient if the grid points are at the places away from the characteristics of the function to be interpolated. As a result, it is necessary to investigate the possibilities of grid redistribution. The redistribution of grid points can be done by a transformation of coordinates. With an introduction of a map u = ρ(x), a new set of grid points xj = ρ−1 (uj ) is obtained [66,93]. The function g(x) can then be interpolated over this new grid, and the efficiency will depend on the choice of the map. To ensure ρ−1 (u) exists, the map ρ(x) must satisfy ρ′ (x) > 0 for all x in the domain of interpolation. Suppose the original interpolating functions, for example CjP (u) or Sj (u), are represented by a generic notation Cj (u). Then with a simple substitution u = ρ(x), the orthogonality of Cj (u) becomes ρ−1 (b) b ηj δij = Ci (u)Cj (u) du = a Ci (ρ(x))Cj (ρ(x))ρ′ (x) dx ρ−1 (a) The new interpolating functions Cˆj (x) = ρ′ (x) Cj (ρ(x)) ρ′ (xj ) −1 ρ (b) will satisfy cardinality Cˆj (xi ) = δij and orthogonality ρ−1 (a) Cˆi (x)Cˆj (x) dx = ηˆj δij , where ηˆj = ηj /ρ′ (xj ). The new quadrature rule becomes N ρ−1 (b) ρ−1 (a) g(x) dx ≈ ηˆj g(xj ). j=1 After mapping, the polynomial basis function becomes Cˆj (x) = ρ′ (x)w(ρ(x)) Ij (ρ(x)) ρ′ (xj )w(ρ(xj )) 14 Chapter 1. Introduction with quadrature weights ηˆj = wj /[ρ′ (xj )w(ρ(xj ))], and the sinc basis function becomes ρ′ (x) Cˆj (x) = Sj (ρ(x)) ρ′ (xj ) with quadrature weights ηˆj = h/ρ′ (xj ). 1.6 Pseudospectral method The spectral method for the solution of a differential equation Lψ(x) = g(x) defined on the domain R with homogeneous boundary conditions involves an expansion of basis functions {φn (x)}∞ n=1 in the form [19] ∞ ψ(x) = cn φn (x). n=1 Spectral convergence is achieved if the magnitudes of the expansion coefficients cn decrease exponentially. The truncated expansion N n=1 cn φn (x) will give a good approximation to ψ(x) if cn are small enough for all n > N . If φn (x) are normalized and satisfy the same homogeneous boundary conditions, the coefficients cn can then be solved by the Galerkin method, which requires that for all m ≤ N , R φm (x)Lψ(x) dx = R φm (x)g(x) dx, or with the expansion, N cn n=1 φm (x)Lφn (x) dx = R φm (x)g(x) dx. (1.21) R Equation (1.21) is a simple matrix equation Lc = g, where Lmn = R φm (x)Lφn (x) dx is a matrix element of the differential operator L, and gm = R φm (x)g(x) dx. For simplicity, quadrature rules are used for the evaluation of the integrals. The quadrature points are regarded to as the grid points for which the values of the solution are defined. The spectral method with quadrature evaluation of integral is referred to as the pseudospectral method. The main difference between spectral and pseudospectral methods is that spectral method is based on the expansion of basis functions, while pseudospectral method is based on a grid of points. 15 Chapter 1. Introduction The basis functions used in pseudospectral methods are usually interpolating functions Cj (x), mapped or unmapped, which satisfy the cardinality √ Cj (xi ) = δij at the quadrature points xi . With normalized basis Cj (x)/ ηj , the solution of the differential equation Lψ(x) = g(x) is expanded in the form N Cj (x) √ ψ(x) = ηj ψ(xj ) √ . ηj j=1 When the integral from the Galerkin method, R Ci (x) √ Lψ(x) dx = ηi R Ci (x) √ g(x) dx, ηi is evaluated by quadrature, the resulting equation becomes Lψ(x)|x=xi = g(xi ). Hence, the pseudospectral method is often considered as a collocation method. Eigenvalue problems Lψ (m) (x) = λψ (m) (x) are solved in a similar way by pseudospectral method. After evaluation of matrix elements by quadrature, the eigenvalue problem becomes N √ √ Lij ηj ψ (m) (xj ) = λm ηi ψ (m) (xi ), j=1 where Lij = (ηi ηj )−1/2 R Ci (x)LCj (x) dx. The eigenvalues λm are found by √ diagonalizing L. The corresponding eigenvectors are ηj ψ (m) (xj ). For time dependent problem ∂ψ(x, t)/∂t = Lψ(x, t), the expansion N ψ(x, t) = aj (t) j=1 Cj (x) √ ηj . √ is used, where aj (t) = ηj ψ(xj , t). With the pseudospectral method, the differential equation become a′j (t) = N j=1 Lij aj (t). Time integration can be done with other methods such as the backward Euler or Runge-Kutta methods. 16 Chapter 1. Introduction 1.7 Thesis objective One of the main objectives of this thesis is to present the equivalence of all the different pseudospectral methods discussed in Section 1.4 from a single approach. Although these pseudospectral methods have different backgrounds and theories, they are all equivalent to a method involving the expansion of orthogonal interpolating functions discussed in Sections 1.5 and 1.6. This method will serve as the central principle in all developments and calculations in the rest of the thesis. Another objective is to study the performance of the pseudospectral method for different problems using different basis functions. Throughout the thesis, I will demonstrate the exponential convergence of eigenvalues using various pseudospectral methods. Exponential, or spectral, convergence refers to the convergence of the eigenvalue calculated by N quadrature points, λ(N ) , to the exact eigenvalue, λexact , at an exponential rate, that is λ(N ) − λexact ≈ 10−qN for some positive constant q [37, 42]. In Chapter 2, I will provide detailed comparisons between polynomial and sinc methods using applications to Fokker-Planck and Schr¨ odinger eigenvalue problems including electron thermalization in a background gas, bistable models for chemical reactions, and vibrational states of I2 approximated by the Morse potential. The exact ground-states eigenfunctions of all these problems are known, and the weight functions of the orthogonal polynomials can therefore be chosen to be the square of the ground-state eigenfunctions. Wei claimed that the sinc method yields a faster rate of convergence to the exact eigenvalues than these polynomial methods. I will show in Chapter 2 that my results are opposite to Wei’s conclusion. For vibrational states of noble gas diatoms studied in Chapters 3 or 4, the ground-state eigenfunctions are unknown. Because the potential functions of these systems are close to the Morse potential with appropriately fitted parameters, we will see that polynomials based on the Morse ground-state eigenfunction gives better convergence rates than other classical methods for the lower-states eigenvalues of noble gas systems. The last objective is to study the improvements of the performance of existing basis functions with a transformation of coordinates. This can be 17 Chapter 1. Introduction understood as the redistribution of interpolation grid points. In Chapter 3, I will show that for highly excited vibrational states of I2 where the eigenfunctions are loosely bound, a linear mapped (scaling and translation) polynomial basis yields a much faster convergence rate over the same set of polynomials without mapping. The idea of transformation is extended to nonlinear mappings in Chapter 4. With a hyperbolic sine mapping, the convergence rates for eigenvalues obtained in He2 system are more than ten times faster than without mapping. In Chapter 5, we will discuss the optimization of explicitly defined transformations such as the hyperbolic sine mapping. A new flexible type of adaptive transformation will also be introduced, which depends only on the shape of the function to be interpolated. I will apply this adaptive transformation to a time-dependent problem and develop an adaptive moving mesh scheme for the psuedospectral method. Furthermore, the possibility of applying this adaptive transformation to eigenvalue problems will be presented. 18 Bibliography [1] P. Amore. A variational sinc collocation method for strong coupling problems. J. Phys. A: Math. Gen., 39:L349, 2006. [2] U. M. Ascher, R. M. M. Mattheij, and D. Russell, Robert. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Prentice Hall, Englewood Cliffs, 1988. [3] P. W. Atkins. Molecular Quantum Mechanics, 3rd Ed. Oxford, New York, 1997. [4] K. Atkinson and W. Han. Theoretical Numerical Analysis: a Functional Analysis Framework. Springer-Verlag, New York, 2001. [5] K. E. Atkinson. An Introduction to Numerical Analysis. John Wiley & Sons, New York, 1978. [6] S. K. Banik, J. R. Chaudhuri, and D. S. Ray. The generalized Kramers theory for nonequilibrium open one-dimensional systems. J. Chem. Phys., 112:8330, 2000. [7] T. Barakat, K. Abodayeh, and O. M. Al-Dossary. Exact solutions for vibrational levels of the Morse potential via the asymptotic iteration method. Czech. J. Phys., 56:583, 2006. [8] D. Baye. Lagrange-mesh method for quantum-mechanical problems. Phys. Stat. Sol. (b), 243:1095, 2006. [9] D. Baye and P. H. Heenen. Generalised meshes for quantum mechanical problems. J. Phys. A: Math. Gen., 19:2041, 1986. [10] D. Baye and M. Vincke. Lagrange meshes from nonclassical orthogonal polynomials. Phys. Rev. E, 59:7195, 1999. [11] O. Bayrak and I. Boztosun. 19 Bibliography [12] A. M. Berezhkovskii, V. Zitserman, and A. Polimeno. Numerical test of Kramers reaction rate theory in two dimensions. J. Chem. Phys., 105:6342, 1996. [13] D. J. Bicout and A. Szabo. First passage times, correlation functions, and reaction rates. J. Chem. Phys., 106:10292, 1997. [14] R. Blackmore and B. Shizgal. A solution of Kramers equation for the isomerization of n-butane in CCl4 . J. Chem. Phys., 83:2934, 1985. [15] M. Braun, S. A. Sofianos, D. G. Papageorgiou, and I. E. Lagaris. An efficient Chebyshev-Lanczos method for obtaining eigensolutions of the Schr¨ odinger equation on a grid. J. Comput. Phys., 126:315, 1996. [16] R. L. Burden. Numerical Analysis. Brooks/Coles, Pacific Grove, 2001. [17] C. E. Burkhardt. Vibration-rotation coupling in a Morse oscillator. Am. J. Phys, 75:686, 2007. [18] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, New York, 1988. [19] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods: Fundamentals in Single Domains. Springer-Verlag, Berlin, 2006. [20] H. Chen and B. D. Shizgal. A spectral solution of the Sturm-Liouville equation: comparison of classical and nonclassical basis set. J. Comput. Appl. Math., 136:17, 2001. [21] X. Chen. The shooting method for solving eigenvalue problems. J. Math. Anal. Appl., 203:435, 1996. [22] O. Christiansen. Vibrational structure theory: new vibrational wave function methods for calculation of anharmonic vibrational energies and vibrational contributions to molecular properties. Phys. Chem. Chem. Phys., 9:2942, 2007. [23] W. T. Coffey, Y. P. Kalmykov, and S. V. Titov. Solution of the master equation for Wigner’s quasiprobability distribution in phase space for the Brownian motion of a particle in a double well potential. J. Chem. Phys., 127:074502, 2007. 20 Bibliography [24] D. T. Colbert and W. H. Miller. A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method. J. Chem. Phys., 96:1982, 1992. [25] J. Coville and L. Dupaigne. Propagation speed of travelling fronts in non local reation-diffusion equations. Nonlin. Anal., 219:219, 2006. [26] A. Das. Lectures on Quantum Mechanics. Hindustan, New Delhi, 2003. [27] R. Dawes and T. Carrington. A multidimensional discrete variable representation basis obtained by simultaneous diagonalization. J. Chem. Phys., 121:726, 2004. [28] L. Debnath. Nonlinear Partial Differential Equations for Scientists and Engineers. Birkh¨ auser, Boston, 1997. [29] P. Deuflhard and A. Hohmann. Numerical Analysis: A First Course in Scientific Computation. Walter de Gruyter, Berlin, 1995. [30] V. M. Devigne, I. S. Pop, C. J. van Duijn, and T. Clopeau. A numerical scheme for the pore-scale simulation of crystal dissolution and precipation in porous media. SIAM J. Numer. Anal., 46:895, 2008. [31] Y.-L. Dong, B. Zhao, and J. Zheng. Numerical investigation of nonlocal electron transport in laser produced plasmas. Chin. Phys., 16:3742, 2007. [32] S. R. Dunbar. Travelling wave solutions of diffusive Lotka-Volterra equations. J. Math. Bio., 17:11, 1983. [33] J. Dunkel, W. Ebeling, L. Schimansky-Geier, and P. Hangii. Kramers problem in evolutionary strategies. Phys. Rev. E, 67:061118, 2003. [34] H. M. El-Hawary and E. O. Abdel-Rahman. Numerical solution of te generalized Burger’s equation via spectral/spline methods. Appl. Math. Comput., 170:267, 2005. [35] G. Feng, X. Gao, W. Dong, and J. Li. Time-dependent solutions of the Fokker-Planck equation of maximally reduced air-sea coupling climate model. Chaos, Solitons & Fractals, 37:487, 2008. [36] L. Ferragut, M. I. Asensio, and S. Monedero. A numerical method for solving convection-reaction-diffusion multivalued equations in fire spread modelling. Adv. Eng. Soft., 38:366, 2007. 21 Bibliography [37] G. Freud. Orthogonal Polynomials. Pergamon, Oxford, 1971. [38] M. L. Gandarias and M. S. Bruz´ on. Solutions through nonclassical potential symmetries for a generalized inhomogeneous nonlinear diffusion equation. Math. Meth. Appl. Sci., 31:753, 2008. [39] W. Gautschi. Orthogonal Polynomials: Computation and Approximation. Oxford, New York, 2004. [40] A. Gerisch and M. A. J. Chaplain. Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion. J. Theo. Bio., 250:684, 2008. [41] G. Goswami, P. Majee, P. K. Ghosh, and B. C. Bag. Colored multiplicative and additive non-Gaussian noise-driven dynamical system: mean first passage time. Physica A, 374:549, 2007. [42] D. Gottlieb and S. A. Orszag. Numerical Analysis of Spectral Methods: Theory and Applications. SIAM, Philadelphia, 1977. [43] D. Goulding, D. Curtin, T. Piwonski, J. Houlihan, J. P. Gleeson, and G. Huyet. Kramers’ law for a bistable system with time-delayed noise. Phys. Rev. E, 76:031128, 2007. [44] P. Grindrod. Theory and Application of Reaction-Diffustion Equations: Patterns and Waves. Clarendon, New York, 1996. [45] D. T. Haar. Master of Modern Physics: The Scientific Contributions of H. A. Kramers. Princeton, Princeton, 1998. [46] E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems. Springer-Verlag, Berlin, 1993. [47] P. H¨ anggi, P. Talkner, and M. Borkovec. Reaction-rate theory: fifty years after Kramers. Rev. Mod. Phys., 62:251, 1990. [48] J. D. Hoffman. Numerical Methods for Engineers and Scientists. Marcel Dekker, New York, 2001. [49] A. Iomin. Toy model of fractional transport of cancer cells due to self-entrapping. Phys. Rev. E, 73:061918, 2006. [50] H. Ishikawa. Numerical methods for the eigenvalue determination of second-order ordinary differential equations. J. Comput. Appl. Math., 208:404, 2007. 22 Bibliography [51] S. Janjai, N. Lamlert, P. Intawee, B. Mahayothee, M. Haewsungcharern, B. K. Bala, and J. M¨ uller. Finite element simulation of drying of mango. Biosys. Eng., 99:523, 2008. [52] H. W. Jones and J. L. Jain. Eigenvalues to arbitrary precision for one-dimensional Schr¨ odinger equations by the shooting method using integer arithmetic. Int. J. Quant. Chem., 80:842, 2000. [53] Z. Kalogiratou, T. Monovasilis, and T. E. Simos. Computation of the eigenvalues of one-dimensional Schr¨ odinger equation by symplectic methods. Int. J. Quant. Chem., 106:795, 2006. [54] J. H. Karl. Introduction to Digital Signal Processing. Academic, San Diego, 1989. [55] J. P. Killingbeck, N. A. Gordon, and M. R. M. Witwit. Stable forward shooting for eigenvalues and expectation values. Phys. Lett. A, 206:279, 1995. [56] J. P. Killingbeck and G. Jolicard. The eighth order Numerov method. Phys. Lett. A, 261:40, 1999. [57] V. Kokoouline, O. Dulieu, R. Kosloff, and F. Masnou-Seeuws. Mapped Fourier methods for long-range molecules: application to perturbations in the Rb2 (0+ u ) photoassociation spectrum. J. Chem. Phys., 110:9865, 1999. [58] H. A. Kramers. Brownian motion in a field of force and the diffusion model of chemical reactions. Physica VII, 4:284, 1940. [59] H. S. Lee and J. C. Light. Molecular vibrations: iterative solution with energy selected bases. J. Chem. Phys., 118:3458, 2003. [60] W.-T. Li and S.-L. Wu. Traveling waves in a diffusive predator-prey model with holling type-III functional response. Chaos, Solitons & Fractals, 37:476, 2008. [61] J. C. Light and T. Carrington. Discrete variable representations and their utilization. Adv. Chem. Phys., 114:263, 2000. [62] J. C. Light, I. P. Hamilton, and J. V. Lill. Generalized discrete variable approximation in quantum mechanics. J. Chem. Phys., 82:1400, 1985. 23 Bibliography [63] J. V. Lill, G. A. Parker, and J. C. Light. Discrete variable representations and sudden models in quantum scattering theory. Chem. Phys. Lett., 89:483, 1982. [64] X.-S. Liu, Y.-H. Chi, and P.-Z. Ding. Symplectic schemes and the shooting method for eigenvalues of the Schr¨ odinger equation. Chin. Phys. Lett., 21:1681, 2004. [65] X.-S. Liu and P.-Z. Ding. Numerical method based on Magnus expansion and a new shooting method for eigenvalues of Schr¨ odinger equation. Int. J. Quant. Chem., 103:149, 2005. [66] J. Lund and K. L. Bowers. Sinc Methods for Quadrature and Differential Equations. SIAM, Philadelphia, 1992. [67] P. A. Markowich. Applied Partial Differential Equations. SpringerVerlag, Berlin, 2007. [68] L. A. Martino, A. Osella, C. Dorso, and J. L. Lanata. Fisher equation for anisotropic diffusion: simulating South American human dispersals. Phys. Rev. E, 76:031923, 2007. [69] C. Masoller. Noise-induced resonance in delayed feedback systems. Phys. Rev. Lett., 88:034102, 2002. [70] D. A. Mazziotti. Spectral difference methods for solving the differential equations of chemical physics. J. Chem. Phys., 117:2455, 2002. [71] R. Meyer. Trigonometric interpolation method for one-dimensional quantum-mechanical problems. J. Chem. Phys., 52:2053, 1970. [72] T. Monovasilis, Z. Kalogiratou, and T. E. Simos. Computation of the eigenvalues of the Schr¨ odinger equation by symplectic and trigonometrically fitted symplectic partitioned Runge-Kutta methods. Phys. Lett. A, 372:569, 2008. [73] P. M. Morse. Diatomic molecules according to the wave mechanics II: vibrational levels. Phys. Rev., 34:57, 1929. [74] J. D. Murray. Mathematical Biology I: an Introduction, 3rd Ed. Springer-Verlag, New York, 2003. [75] J. D. Murray. Mathematical Biology II: Spatial Models and Biomedical Applications, 3rd Ed. Springer-Verlag, New York, 2003. 24 Bibliography [76] C. Nicolis. Long-term climatic transitions and stochastic resonance. J. Stat. Phys., 70:3, 1993. [77] D. Olmos and B. D. Shizgal. A pseudospectral method of solution of Fisher’s equation. J. Comput. Appl. Math., 219:219, 2006. [78] A. G. Peeters and D. Strintzi. The Fokker-Planck equation, and its application in plasma physics. Ann. Phys. (Berlin), 17:142, 2008. [79] E. Pollak and P. Talkner. Reaction rate theory: what it was, where is it today, and where is it going? Chaos, 15:026116, 2005. [80] A. I. M. Rae. Quantum Mechanics, 5th Ed. Taylor & Francis, New York, 2008. [81] H. Risken. The Fokker-Planck Equation: Methods of Solution and Application. Springer-Verlag, New York, 1984. [82] P. Rosenau. Fast and superfast diffusion processes. Phys. Rev. Lett., 74:1056, 1995. ´ [83] C. A. S´ anchez-Castellanos, Amezcua-Eccius, O. Alvarez Bajo, and R. Lemus. A local-normal description of vibrational excitations of pyramidal molecules in terms of Morse oscillators. J. Mol. Struc. (Theochem), 247:140, 2008. [84] J. M. Sanz-Serna and A. Portillo. Classical numerical integrators for wave-packet dynamics. J. Chem. Phys., 104:2349, 1996. [85] S. Savovi´c and A. Djordjevich. Numerical solution of the diffusion equation describing the flow of radon through concrete. Appl. Radiat. Isot., 66:552, 2008. [86] O. S´ero-Guillaume, S. Ramezani, J. Margerit, and D. Calogine. On large scale forest fires propagation models. Int. J. Therm. Sci., 47:680, 2008. [87] C. E. Shannon. Communication in the presence of noise. Proc. IRE, 37:10, 1949. [88] S. Sharma, P. Sharma, and H. Singh. Quantum control of vibrational excitations in a heteronuclear diatomic molecule. J. Chem. Sci, 119:433, 2007. 25 Bibliography [89] B. Shizgal and R. Blackmore. Eigenvalues of the Boltzmann collision operator for binary gases: relaxation of anisotropic distributions. Chem. Phys., 77:417, 1983. [90] B. Shizgal and R. Blackmore. A discrete ordinate method of solution of linear boundary value and eigenvalue problems. J. Comput. Phys., 55:313, 1984. [91] B. Shizgal and R. Blackmore. Discrete ordinate method of solution of a Fokker-Planck equation with a bistable potential. Chem. Phys. Lett., 109:242, 1984. [92] B. D. Shizgal and H. Chen. The quadrature discretization method (QDM) in the solution of the Schr¨ odinger equation with nonclassical basis functions. J. Chem. Phys., 104:4137, 1996. [93] F. Stenger. Numerical Methods Based on Sinc and Analytic Functions. Springer-Verlag, New York, 1993. [94] W. A. Strauss. Partial Differential Equations: An Introduction. John Wiley & Sons, New York, 1992. [95] M. St¨ urmer, H. K¨ ostler, and U. R¨ ude. A fast full multigrid solver for applications in image processing. Numer. Linear Algebra Appl., 15:187, 2008. [96] V. Szalay. The generalized discrete variable representation. An optimal design. J. Chem. Phys., 105:6940, 1996. [97] H. Ta¸seli and H. Alıcı. The scaled Hermite-Weber basis in the spectral and pseudospectral pictures. J. Math. Chem., 38:367, 2005. [98] W. Thirring. Quantum Mathematical Physics: Atoms, Molecules and Large Systems, 2nd Ed. Springer-Verlag, Berlin, 2002. [99] N. G. van Kampen. Stochastic Processes in Physics and Chemistry, 3rd Ed. Elsevier, Amsterdam, 2007. [100] J. L. V´ azquez. The Porous Medium Equation: Mathematical Theory. Oxford, New York, 2007. [101] A.-M. Wazwaz. Analytic study on Burger, Fisher, Huxley equations and combined forms of these equations. Appl. Math. Comput., 195:754, 2008. 26 Bibliography [102] G. W. Wei. Solving quantum eigenvalue problems by discrete singular convolution. J. Phys. B: At. Mol. Opt. Phys., 33:343, 2000. [103] H. Wei. Ghost levels and near variational forms of the discrete variable representation: application to H2 O. J. Chem. Phys., 106:6885, 1997. [104] E. T. Whittaker. On the functions which are represented by the expansions on the interpolation theory. Proc. Roy. Soc. Edinburgh, 35:181, 1915. [105] K. Willner, O. Dulieu, and F. Masnou-Seeuws. Mapped grid methods for long-range molecules and cold collisions. J. Chem. Phys., 120:548, 2004. [106] W. L. Wood. An exact solution for Burger’s equation. Commun. Numer. Meth. Engng., 22:797, 2006. [107] J. Q. Xia, J. Y. Lo, K. Yang, C. E. Floyd, and J. M. Boone. Dedicated breast computed tomography: volume image denoising via a partial diffusion equation based technique. 35:1950, 2008. [108] R. S. Zola, M. K. Lenzi, L. R. Evangelista, E. K. Lenzi, L. S. Lucena, and L. R. de Silva. Exact solutions for a diffusion equation with a nonlinear external force. Phys. Lett. A, 372:2359, 2008. 27 Chapter 2 Spectral convergence of the quadrature discretization method in the solution of the Schr¨ odinger and Fokker-Planck equations: Comparison with sinc method 2.1 Introduction The solution of the Fokker-Planck and Schr¨ odinger equation with a discretization of the wave function on a grid of points has been of ongoing interest for over two decades [1, 4, 5, 13, 17, 21, 23, 27, 32, 33, 38, 47, 50, 53, 57, 62, 63, 66, 68–70]. Shizgal and co-worker [53, 57, 62, 63] developed a discrete ordinate method later referred to as the quadrature discretization method (QDM), which is a pseudospectral (collocation) method based on nonclassical polynomials. For the solution of the Schr¨ odinger equation, Light and co-workers introduced the discrete variable representation (DVR) [27,38], which is often based on classical polynomials including Gaussian and Fourier basis functions [1, 9, 64] and has been used exclusively by several researchers [23, 32, 50, 72]. The basic ideas of Fourier methods have been reviewed by Kosloff [33]. Wei introduced the discrete singular convolution (DSC) approach for solutions of the Fokker-Planck [68, 69] and Schr¨ odinger equations [70]. Although Wei’s DSC-sinc method is discussed A version of this chapter has been published. Joseph Lo and Bernie D. Shizgal. J. Chem. Phys., 125:194108, 2006. 28 Chapter 2. Spectral convergence of quadrature discretization method in terms of wavelets and signal analysis it is the same as the sinc collocation method (SCM) presented recently by Amore [1]. A Lagrange mesh method was introduced by Baye et al [5]. Other approaches for quantum problems include the work by Balint-Kurti and Pulay [4] and the Fourier based method of Colbert and Miller [17]. These discrete methods employed by different researchers have similar features but also many subtle differences. Mazziotti recently introduced a new class of spectral difference method that yield sparse matrices and spectral convergence [42]. The most important feature of any numerical methodology is the rate of convergence of the solution versus the number of grid points or the number of basis functions. The main objective of the present chapter is to provide a detailed study of the rate of convergence of the QDM [53, 57, 62, 63] for the Fokker-Planck and Schr¨ odinger equations. We compare the results with the sinc method by Wei [68]. Boyd [10] has recently compared the efficiency of DSC-sinc [68] and related Fourier based methods with finite difference algorithm. In this chapter, we provide a demonstration of the spectral convergence of the eigenvalues obtained with the QDM for several model systems and provide some insight as to the origin of the rapid convergence. Spectral convergence refers to the exponential convergence of solutions of differential equations versus the number of basis functions or the number of grid points. The mathematical background of spectral methods is discussed at length elsewhere [9, 12, 26, 30, 46]. The use of a Gaussian quadrature was developed by Shizgal [54–56,60] to efficiently discretize the integral operator in the Boltzmann equation analogous to techniques in neutron transport theory [24] and radiative transport [14] introduced over half a century ago. The evaluation of the matrix elements of the linear Boltzmann collision operator in some basis set (typically Laguerre or Hermite polynomials) leads to a cumbersome algebraic problem and the final expressions involve considerable round-off errors in their numerical evaluation [31, 40, 59, 67]. The quadrature approach is far more efficient and accurate. As the Boltzmann equation is an integrodifferential equation, Shizgal and Blackmore [8, 57] extended the formalism to the discretization of derivative operators. This technique was later applied to the Schr¨ odinger equation [62] and the Fokker-Planck equation [63] and referred to as the QDM. Fokker-Planck equations arise as a limit of the Boltzmann equation [3, 59] and in many problems in nonequilibrium statistical mechanics [48]. 29 Chapter 2. Spectral convergence of quadrature discretization method The QDM differs from the other approaches in that it is based on nonstandard polynomials orthogonal with respect to nonclassical weight functions. It originates from an interest in the solution of kinetic theory problems [54–57,60] and was later applied to quantum problems. The important distinction between these classes of problems is that in kinetic theory the “ground state” eigenvalue λ0 = 0, and the corresponding eigenfunction is the equilibrium Maxwellian distribution. The QDM was developed for the solution of a Fokker-Planck equation for the relaxation of electrons in inert gas moderators [52]. The traditional approach for such kinetic theory problems is to consider the expansion of the distribution function in Laguerre (or Sonine) polynomials. For this application, with the small electron to moderator mass ratio, the expansion in the Laguerre polynomials in the electron energy converges extremely slowly. Shizgal [52] showed that if the “speed” polynomials orthogonal with respect to weight function w(x) = x2 exp(−x2 ), where x is the reduced particle speed, are chosen as basis functions, the convergence of the eigenvalues of the Fokker-Planck operator is extremely rapid. The reason provided in Ref. 52 was that the reduced particle speed was the more appropriate independent variable than the reduced energy as discussed earlier [49, 51]. Moreover, the matrix representation of the Fokker-Planck operator is tridiagonal in this new set, and thus rapid convergence is anticipated. The essential philosophy of the QDM stems from the observation that the choice of the equilibrium distribution function as a weight function in the definition of basis set leads to a rapid convergence of the eigenvalues of the Fokker-Planck operator. The actual implementation of the QDM was later based on a grid defined by the quadrature associated with a nonclassical weight function. The method is thus classified as a pseudospectral (collocation) method [8, 63]. The eigenfunction for the ground state of the Fokker-Planck operator is the equilibrium distribution, which is known. The excited state eigenvalues are then determined by using the square of the ground state eigenfunction as the weight function. Orthogonal polynomials and associated quadrature points are then constructed from this weight function. The choice of the weight function in the determination of the eigenfunctions of the Schr¨ odinger equation is motivated by the relationship with the Fokker-Planck equation. The Schr¨ odinger equations that are derived from a Fokker-Planck equation belong to the class of quantum problems in supersymmetric quantum mechanics for which the ground state is known [18, 29]. This chapter considers the solution of 1D Fokker-Planck and the Schr¨ odinger 30 Chapter 2. Spectral convergence of quadrature discretization method eigenvalue problems. The study of 1D problems has been and continues to be of interest [1, 4, 19, 27, 36, 38, 41, 45, 47, 68, 73, 74] with the objective that any acceleration of the convergence for 1D problems would yield increased efficiencies for higher-dimensional problems. Many applications to multi-dimensional systems employ the tensor product of 1D basis sets and the use of symmetries to reduce the dimensionality of the matrices to be diagonalized. There have been several discussions of the importance of the choice of 1D basis functions for direct product bases in applications to multidimensional problems [36, 38, 45, 47]. The QDM has the potential to provide efficient basis functions (or quadrature points) in each dimension of a multi-dimensional problems so as to significantly reduce the size of the final matrices involved. The QDM has previously been applied to the classic 2D Henon-Heiles problem [62] and performs as well as or better than other methods [11, 47, 68]. It was also applied to the solutions of Kramers equation [7] and to other 2D or 3D problems in kinetic theory and fluid dynamics [16, 58, 76, 77]. We limit the objectives of this chapter to the study of spectral convergence for nonclassical polynomials for several 1D problems. We thus focus our attention on the 1D problems previously studied by Wei [68] so as to study the spectral convergence of the QDM and to compare with this DSCsinc method. This approach is exactly the same as the one described in the recent paper by Amore [1] referred to as the SCM and we prefer this nomenclature. The first problem considered arises from the relaxation of electrons in atomic moderators assuming a hard sphere cross section for the electron-atom momentum transfer cross section [61]. This Fokker-Planck eigenvalue problem can be transformed to an equivalent Schr¨ odinger equation. We solve both eigenvalue problems with nonclassical weight functions as done previously [7, 61–63] but here present results showing the spectral convergence of the QDM in comparison with the convergence rates obtained with the SCM. It is important to note that the rate of convergence of the eigenvalues is generally more rapid for the lower states and slower for the higher states [62, 63]. Wei [68, 69] did not show a comparison of the rates of convergence but suggested that the SCM is superior to the QDM (Table 1 of Ref. 68 and the discussion after Eq. (42)). Table 1 in this paper by Wei [68] suggests that 70 quadrature points are required to achieve the convergence shown for all the eigenvalues. The value of λ1 = 4.86340 in this table is taken from Ref. 62, where it is clearly shown that only six speed polynomial basis functions and not 70 are required to get convergence to six significant figures. Wei compares the QDM value in his Table 1 and shows 31 Chapter 2. Spectral convergence of quadrature discretization method λ1 = 4.683395 with N = 60 for the DSC-Sinc method. The tabulation of these results in this way suggests that the SCM is marginally superior to the QDM, which it is not. The purpose of this chapter is to illustrate the rapid convergence of the QDM and the much slower convergence of the SCM. Since the SCM is designed for problems with x ∈ (−∞, ∞), it is anticipated that the results for the electron relaxation problem will not provide spectral convergence even though the boundary condition at x = 0 is satisfied. The Morse potential has often been used as a benchmark of numerical methods in the solution of the Schr¨odinger equation [4, 6, 11, 15, 39, 65, 70–73]. We also study in this chapter the rates of convergence of the QDM and SCM method as applied to the vibrational states of I2 with a Morse potential. Another application is concerned with the Fokker-Planck operator for the so-called quartic potential [8, 63] for which the equilibrium distribution is bimodal. The corresponding Schr¨ odinger equation is characterized by a triple well potential. This and similar bimodal systems are of considerable interest in modeling dynamical evolutionary problems and chemical reactions [20, 22, 28, 34, 35]. The emphasis in this chapter, not discussed previously, is on the dependence of the rate of convergence on the choice of basis functions and thus on the distribution of grid points. We show that the discussions by Wei [68–70] are incomplete and that the convergence of the QDM, with a nonuniform grid distribution, is considerably faster than the SCM, with a uniform grid distribution, for solutions of the Fokker-Planck equation on [0, ∞). The SCM can provide convergence rates comparable with the QDM for solutions of the Schr¨ odinger equation on (−∞, ∞) in some instances. A complete comparison is provided in this chapter. In Section 2.2, we briefly review the discretization procedures for the Fokker-Planck and Schr¨ odinger eigenvalue problems. The specific applications to the electron relaxation problem, the vibrational states in the Morse potential for I2 , and the quartic bistable potential are discussed in Sections 2.3 – 2.5. We present the numerical results in Section 2.6. In Section 2.7 we provide a discussion concerning the interpretation of the QDM in terms of the phase space coverage of the basis sets used, the structure of the resulting matrix representatives of the operators considered, and the fundamental reason for the rapid spectral convergence of the eigenvalues and eigenfunctions. 32 Chapter 2. Spectral convergence of quadrature discretization method 2.2 Discretization of the Fokker-Planck and Schr¨ odinger eigenvalue problems As the basic methodology of the calculations reported in this chapter have been presented in previous papers [7, 8, 61–63], only a very brief summary of the final results are provided here. Thus, the first problem considers the determination of the eigenvalues, λm , and eigenfunctions, φm (x), of the 1D Fokker-Planck eigenvalue problem of the form Lφm = A(x) dφm d2 φm − B(x) = λm φm , dx dx2 (2.1) where A(x) and B(x) are the drift and diffusion coefficients and take on different explicit forms in each application. The time dependent FokkerPlanck equation for the probability density, P (x, t), can be solved in terms of these eigenvalues and eigenfunctions. The QDM is based on a grid of quadrature points associated with a set of polynomials, {Fn (x)}, orthogonal with respect to weight function w(x). If {xi }N i=1 are the quadrature points defined by FN (xi ) = 0, then the discrete representation of the operator L is [8, 63] N =− LQDM ij B(xk )[Dki + h(xk )δki ][Dkj + h(xk )δkj ], k=1 where the QDM representation of the derivative operator is given by [63] Dij = √ N ′ wi wj Fn−1 (xi )Fn−1 (xj ), n=1 and h(x) = w′ (x) P ′ (x) − 0 , 2w(x) 2P0 (x) where P0 (x) = exp[− [A(x′ )/B(x′ )]dx′ − ln B(x)] is the equilibrium probability density [62]. The operator L is self-adjoint with the scalar product defined with the weight function w(x) = P0 (x) for which h(x) = 0. For this choice of weight function, the eigenvalues are given by the numerical diagonalization of the symmetric matrix N ˆ QDM = − L ij B(xk )Dki Dkj . (2.2) k=1 33 Chapter 2. Spectral convergence of quadrature discretization method The eigenvalue problem, Eq. (2.1), can be transformed to a Schr¨ odinger x ′ −1/2 ′ equation [63] with the transformation y = [B(x )] dx and we get Hψm (y) = − d2 ψm (y) + V (y)ψm (y) = λm ψm (y), dy 2 (2.3) where V (y) is defined in terms of A(x) and B(x) and given explicitly by Eqs. (10) – (12) of Ref. 63. As shown in Refs. 62 and 63, the discrete representation of the Hamiltonian H in Eq. (2.3) is N QDM Hij = k=1 where V0 (y) = Dki Dkj + [V (yi ) − V0 (yi )]δij , 1 w′′ (y) 1 − 2 w(y) 4 w′ (y) w(y) 2 . (2.4) If the weight function w(y) is derived from the “superpotential” [18, 29] associated with the ground state eigenfunction [63] so that V0 (y) = V (y), the discrete representation of the Hamiltonian reduces to N ˆ QDM = H ij Dki Dkj . (2.5) k=1 with no explicit reference to the potential. The SCM [1, 9, 17, 50, 68–70] is based on a uniform grid with spacing h, xj = jh, and basis functions chosen to be the sinc interpolation function Sj,h (x) = sinc(π(x − jh)/h) = sin[π(x − jh)/h][π(x − jh)/h] for x = xj and unity for x = xj . For the problems considered here on the infinite and semi-infinite domains, the interval is cutoff at appropriate values xmin and xmax . The first and second matrix derivative operators on this uniform grid are given by, (−1)i−j ′ if i = j, (1) (i−j)h Dij = Sj,h(xi ) = (2.6) 0 if i = j, and i−j (2) ′′ Dij = Sj,h(xi ) = − 2(−1) (i−j)2 h2 π2 − 3h 2 if i = j, if i = j. (2.7) which are the same expressions as given by Amore [1]. The explicit expressions, Eqs. (2.6) and (2.7), were not provided by Wei in his papers [68–70]. 34 Chapter 2. Spectral convergence of quadrature discretization method The sinc functions are a set of basis functions orthogonal with unit weight function. The matrix representation of the Fokker-Planck operator, LSCM , ij is 1 ∞ = Si,h (x)LSj,h (x) dx LSCM ij h −∞ and thus the discrete representation is given by (1) (2) LSCM = −A(xi )Dij + B(xi )Dij , ij (2.8) which is not symmetric. However, the discrete representation of the Hamiltonian is given by (2) SCM Hij = −Dij + V (xi )δij (2.9) and is symmetric. Even though the SCM is simple to use, it may not be an efficient numerical scheme as the infinite or semi-infinite domains have to be made finite with N grid points in [xmin , xmax ] that are uniform and not related to the potential. 2.3 Electron relaxation One particular application arises for the relaxation of electrons in gases [61] for which the electron-atom cross section is taken to be a hard sphere and B(x) = x and A(x) = 2x2 − 3 where the reduced speed is x = mv 2 /2kB Tb , x ∈ [0, ∞), where v is the electron speed, Tb is the temperature of the background gas, and kB is the Boltzmann constant. The speed polynomials, orthogonal with weight function w(x) = x2 exp(−x2 ) for which h(x) = 0 provide a very rapid convergence of the eigenvalues [52] calculated from the numerical diagonalization of Eq. (2.2). The potential in the Schr¨ odinger equation for this model, V (y) = y 6 /64− 2 + 15/4y , is bounded and there are an infinite number of bound states. The weight function w(y) = y 5 exp(−y 4 /16) can be derived from the superpotential [18, 29] associated with the ground state eigenfunction for this potential [63]. Thus the eigenvalues can also be calculated from the numerical diagonalization of Eq. (2.5). In principle, the eigenvalues defined by Eqs. (2.1) and (2.3) are isospectral although the convergence rates of the eigenvalues determined from Eqs. (2.2) and (2.5) will be different. The distribution of quadrature points for these weight functions is nonuniform. The calculation of the eigenvalues for this problem was considered by Wei [68] y2 35 Chapter 2. Spectral convergence of quadrature discretization method with the SCM in comparison with the earlier studies. The work by Wei did not consider a complete study of the convergence properties of the SCM in comparison with the QDM which is presented in this chapter. 2.4 Vibrational states of I2 A model for the vibrational states of diatomic molecules that has been employed as a benchmark of numerical methods for the solution of the Schr¨ odinger equation [6, 11, 15, 39, 42, 65, 71–73] is the Morse potential V (y) = D(1 − exp(−αy))2 . For I2 , the potential parameters are D = 0.0224 a.u. and α = 0.9374 a.u. In the Schr¨ odinger equation,−( 2 /2µ)(d2 ψm /dx2 ) + V ψm = λm ψm , we take Planck’s constant = 1 a.u. and a reduced mass µ = 119406 a.u. This potential has 77 eigenvalues given exactly by √ 2 α2 1 2 2µD 1 λm = m+ − m+ 2µ 2 α 2 for m ≥ 0. The QDM uses the square of the ground state eigenfunction √ 2µD exp(−αy) αy ψ0 (y) = exp − y+ + (2.10) α 2 as the weight function, i.e. w(y) = [ψ0 (y)]2 . Thus the discretized Hamiltonian with the QDM is QDM Hij = 2 2µ N Dki Dkj + λ0 δij , k=1 2 where λ0 = V − 2µ V0 and V0 arises from the weight function given by Eq. (2.4). The SCM discretization is performed in the computational domain [xmin , xmax ] so the grid points are given by xj = xmin + (j − 1)h, where h = (xmax − xmin )/(N − 1). Wei [70] has reported the application of the SCM to this system and chose xmin = −1.4 and xmax = 1.4. This computational domain is appropriate for lower eigenvalues but larger xmax values are required for the calculation of the higher eigenvalues, namely λ30 and λ50 . In the previous papers by Wei, the eigenvalues only up to m = 20 [71] and m = 24 [70] were calculated. 36 Chapter 2. Spectral convergence of quadrature discretization method 2.5 The quartic bistable system The steady distribution for the time dependent Fokker-Planck equation ∂P (x, t)/∂t = LP (x, t) with the drift and diffusion coefficients defined by, 4 2 A(x) = x3 − x and B(x) = ǫ with x ∈ (−∞, ∞), is P0 (x) = exp[− x4ǫ + x2ǫ ] [8, 63]. In the limit ǫ → 0, the steady distribution is bimodal with narrow peaks at x = ±1. This quartic bistable system finds useful applications as a model for reactive systems [8,20,22,28, 29,34, 35, 48,68]. The eigenfunctions of the linear Fokker-Planck operator, Eq. (2.1), are highly localized and the convergence of the eigenvalues versus N can be slow especially for small ǫ. The eigenvalue spectrum was calculated with the QDM by Blackmore and Shizgal [8] with a weight function corresponding to the equilibrium probability density, wa (y) = P0 (y)/ exp(1/4ǫ). The matrix representation of the Fokker-Planck operator with the polynomials defined by the weight function wa (x) is pentadiagonal, Eq. (B.3), and from the structure of the matrix, Eq. (B.4), it is clear that there is no coupling between eigenfunctions of even and odd symmetries. The corresponding Schr¨ odinger equation is Hψm (y) = −ǫ d2 ψm (y) + V (y)ψm (y) = λm ψm (y) dy 2 (2.11) with (y 3 − y)2 3y 2 − 1 − (2.12) 4ǫ 2 and there is a factor of ǫ in front of the second derivative operator in Eq. (2.11). The potential is characterized by three wells at y = 0 and √ at y = ± 6 + 3 1 + 18ǫ/3 → ±1 as ǫ → 0. Owing to the symmetry of the potential, the eigenvalues occur as a series of singlets and triplets with eigenfunctions of even and odd symmetries (see Eqs. (B.5) and (B.6)) [8,34,35,48]. In the limit ǫ → 0, the eigenvalue spectrum coincides with the harmonic osˆ m = m and those for m even are triply degenerate except for cillator states, λ m = 0 which is doubly degenerate. For small values of ǫ, the calculation of the splitting of these degenerate levels is a difficult computational problem. Shizgal and Chen [63] employed nonclassical weight functions to improve the rate of convergence. These authors showed that the addition of a Gaussian weight function centered at y = 0, that is wc (y) = wa (y) + exp(−y 2 /2ǫ), is required to resolve the eigenvalues localized in the middle well. Wei [68] applied the SCM to this system and presented an incomplete comparison with these previous works. The main purpose of the present study is to provide a detailed comparison of the QDM and SCM with particular attention paid to the rate of convergence. In the present work, we introduce a novel V (y) = 37 Chapter 2. Spectral convergence of quadrature discretization method technique for calculating the eigenvalues corresponding to eigenfunctions of even or odd symmetry by using weight functions over the semi-infinite interval [0, ∞). For the even eigenstates, we use wbe (y) = wa (y), y ∈ [0, ∞), and for the odd eigenstates, wbo (y) = y 2 wa (y), y ∈ [0, ∞). We occasionally refer collectively to both weight functions as wb (y). We focus on the five lowest eigenvalues, namely λ1 – λ5 . The first nonzero eigenvalue λ1 , which becomes very small as ǫ decreases, is related to the rate coefficient for the passage of particles from one side of the quartic potential to the other side [8, 20, 22, 28, 34, 35, 48]. This eigenvalue, nearly degenerate with the eigenvalue λ0 = 0, represents the rate coefficient for the reaction modeled ˆ1 = 1 with this potential. The eigenvalue λ2 is a singlet state and tends to λ as ǫ → 0. The states λ3 , λ4 and λ5 make up a triplet which are degenerate ˆ 2 = 2. The remainder of the eigenvalue spectrum for ǫ → 0 and tend to λ then alternates between singlet and triplet states. The rate of convergence versus the size of the basis set or equivalently the number of grid points is not uniform for all the eigenvalues. The convergence is sensitive to the nature of the eigenfunctions and in particular their symmetry. 2.6 Convergence of eigenvalues; discussion of results We have applied the QDM and the SCM to a solution of the Fokker-Planck and Schr¨ odinger equations for several model systems as described in the previous sections. We present in this section a detailed comparison of the two discretization methods. As discussed in Section 2.2, the quadrature points with the QDM based on nonclassical weight functions are non-equidistant. The choice of basis functions and thus the distribution of the quadrature points play an important role on the rate of convergence of the solution. One requires the appropriate weight function in order to construct an accurate and rapidly convergent numerical scheme. For the Fokker-Planck equation, a good choice of weight function is the equilibrium solution, that is, w(x) = P0 (x). Since h(x) = 0 for this weight function, the eigenvalues are calculated from the numerical diagonalization of Eq. (2.2). For the Schr¨ odinger equation, the weight function w(y) is chosen such that V0 (y) = V (y) and the matrix representative of the Hamiltonian given by Eq. (2.5) is diagonalized. We study the rate of convergence of the eigenvalues for the electron relaxation problem, the vibrational states for I2 modeled with a Morse potential and the quartic bistable system. The exact eigenval38 Chapter 2. Spectral convergence of quadrature discretization method ues, λexact , used to evaluate the error are calculated to over 20 significant m figures in extended precision with up to 150 quadrature points for the QDM. As a measure of the accuracy of the eigenvalues for a particular N we report (N ) ) ǫ(N m λm − λexact m = log , λexact m (2.13) (N ) which is approximately equal to the number of significant figures in λm . (N ) (N ) For ease of notation, we henceforth write λn and ǫn as λn and ǫn , respectively. Although we compute the eigenvalues from a diagonalization of QDM the discrete representations, LQDM , Eq. (2.2), and Hij , Eq. (2.5), it is ij often instructive to transform to the polynomial basis representation. If the matrix representatives in the polynomial basis sets are diagonally dominant and sparse, then rapid convergence is anticipated. The eigenvalue spectrum of the electron relaxation problem is completely discrete and there is an infinite number of bound states. Figure 2.1 shows a comparison of the rates of convergence of λ1 and λ6 for the electron relaxation problem with the QDM and the SCM. We choose these eigenvalues as representative of the lower order eigenvalues. The convergence with the QDM for the Fokker-Planck equation (triangles) and Schr¨ odinger equation (squares) is similar and spectral [9], and machine accuracy is achieved with approximately 20 basis functions, although accurate values with less precision are obtained with fewer basis functions. We discuss the spectral convergence of λ1 in Section 2.7. Also shown in Fig. 2.1 (the curve labeled RV) are the results with the representation of Fokker-Planck operator in Laguerre basis set in the speed variable as derived by Risken and Voigtlaender [49]. This representation also depends on a scaling variable, which we choose as α = 0.5. A discussion of the use of such scaling factors was presented elsewhere [2, 25, 37, 61]. For the Fokker-Planck operator, the weight function odinger equation we use used is w(x) = x2 exp(−x2 ), whereas for the Schr¨ w(y) = y 5 exp(−y 4 /16). The matrix representative of the Fokker-Planck operator with the nonclassical polynomials defined with w(x) = x2 exp(−x2 ) is a diagonally dominant tridiagonal matrix [52] as shown in Eq. (A.8) and (A.9) in Appendix A. The discretization of the Fokker-Planck eigenvalue problem with the SCM yields a non-symmetric matrix, Eq. (2.8), which possesses many complex eigenvalues. Hence the rate of convergence is very slow and we restrict our comparisons with the SCM to those based on the Schr¨ odinger equation. The convergence of the eigenvalues with the SCM (diagonalization of Eq. (2.9)) is clearly considerably slower. The main reason 39 Chapter 2. Spectral convergence of quadrature discretization method A 0 B 0 SCM (se) SCM (se) -4 ǫ1 ǫ6 -4 -8 RV -8 RV QDM (fpe) -12 -16 QDM (fpe) QDM (se) 0 20 40 N -12 60 -16 QDM (se) 0 20 40 60 N Figure 2.1: Variation of ǫm versus the number of quadrature points N for the electron thermalization problem for (A) m = 1 and (B) m = 5. The SCM results are for the Schr¨ odinger equation (denoted by filled circles, xmax = 6). The QDM results are for the Schr¨ odinger equation (se, denoted with squares) and the Fokker-Planck equation (fpe, denoted with triangles). The curve labeled RV is the result with the Laguerre polynomials in speed [49]. 40 Chapter 2. Spectral convergence of quadrature discretization method for the slow convergence of the SCM is that the sinc functions are orthogonal with respect to the unit weight function on (−∞, ∞), and for the electron relaxation problem the domain is [0, ∞). With the QDM, we need N = 7 to get λ1 = 4.683395 converged to seven significant figures and N = 17 for λ6 = 40.05238. Table 1 in the paper by Wei [68] suggesting that N = 70 is required for the convergence of these eigenvalues is incorrect. We confirm Wei’s results [68] that to calculate these eigenvalues to seven significant figures requires N = 60 for the SCM. The results in Fig. 2.1 show the rapid decrease in the error versus N for the QDM and the much slower convergence with the SCM. The superiority of the QDM over the SCM for this application is thus clear. The plus symbols in Fig. 2.1 correspond to the results obtained with the Laguerre basis set in speed and appropriately scaled. The convergence of the QDM and the SCM for the vibrational states of I2 modeled with the Morse potential is shown in Fig. 2.2. The convergence for the QDM with the square of the ground state eigenfunction as the weight function w(y), Eq. (2.10), is extremely rapid for λ1 . For λ10 , the QDM requires only 24 points to reach machine precision. It is considerably faster than the convergence for SCM with the computational domain [−1.4, 1.4] used by Wei [70]. The convergence for λ30 is somewhat faster with QDM than with the SCM and it is about the same for λ50 . Although it is useful to show the spectral convergence to machine accuracy in Fig. 2.2, it is the convergence to moderate accuracy for smaller N that is of practical interest. To achieve six significant figures for these vibrational states, λ1 = 8.52997 × 10−4 , λ10 = 5.62326 × 10−3 , λ30 = 1.40897 × 10−2 , and λ50 = 1.96125 × 10−2 a.u. require N = 4, 19, 58 and 119 with the QDM. The corresponding values for the SCM are N = 33, 48, 84 and 133. The QDM is clearly superior to the SCM for the lower and moderately excited states, but the results become comparable for the highest state studied. In the previous papers by Wei, the eigenvalues only up to m = 20 [71] and m = 24 [70] were calculated. We also find that the differences λm − λexact for n m = 1, 10, and 30 are at machine accuracy with N = 6, 24, and 70, respectively. This convergence is faster than the result reported by Mazziotti [42]. Figure 2.3 shows the variation of the potential in the Schr¨ odinger, Eq. (2.12), and the distribution of the grid points for the quartic bistable problems. For ǫ = 0.1 in Fig. 2.3A, there are two moderately deep wells in the vicinity of y = ±1 and a shallow well at y = 0. The quadrature points for the supersymmetric weight function wa (y) = exp[−(y 4 /4 − y 2 /2)/ǫ]/ exp(1/4ǫ) are almost uniformly distributed but with some localization in the vicinity 41 Chapter 2. Spectral convergence of quadrature discretization method A 0 B 0 SCM [-1.4, 1.4] -4 -4 -8 ǫ10 ǫ1 QDM -8 SCM [-1.4, 1.4] QDM -12 -16 -12 0 10 20 30 40 -16 50 0 20 N 40 60 N C 0 D 0 SCM [-1.4, 2.4] -4 ǫ50 ǫ30 -4 -8 QDM -8 QDM SCM [-1.4, 3.9] -12 -16 20 -12 40 60 80 N 100 -16 40 80 120 160 N Figure 2.2: Variation of ǫm for (A) m = 1, (B) m = 10, (C) m = 30, and (D) m = 50 versus the number of quadrature points N for the I2 Morse potential. The computational domain [xmin , xmax ] for the SCM are also shown. 42 Chapter 2. Spectral convergence of quadrature discretization method 4 A V(y) 2 0 QDM (a) -2 SCM -4 -6 -2 -1 0 1 2 y B V(y) C V(y) 40 4 20 0 0 QDM (a) QDM (a) QDM (c) QDM (b) -4 -20 SCM -1 SCM 0 y 1 -40 -1 0 1 y Figure 2.3: Variation of the potential, V (y) = (y 3 − y)2 /4ǫ − (3y 2 − 1)/2, in the Schr¨ odinger equation for the quartic bistable system: (A) ǫ = 0.1, (B) ǫ = 0.01, (C) ǫ = 0.001. The distributions of grid points are shown for the SCM and the QDM for the different weight functions: (a) wa (y) = exp[−(y 4 /4 − y 2 /2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞, ∞), (b) wbe (y) = wa (y), y ∈ [0, ∞) and (c) wc (y) = wa (y) + exp(−y 2 /2ǫ), y ∈ (−∞, ∞). 43 Chapter 2. Spectral convergence of quadrature discretization method of y = ±1. The grid points for the SCM are uniform as shown. By contrast, for ǫ = 0.01 in Fig. 2.3B, there are three deep wells in the vicinity of y = ±1 and at y = 0. The quadrature points with wa (y) are concentrated near the outermost potential wells and less dense near the potential about y = 0. For wbe (y) = wa (y), y ∈ [0, ∞) the quadrature points are denser near the origin and in the potential well near y = 1. The SCM points are distributed uniformly. In Fig. 2.3C with ǫ = 0.001, the quadrature points with the supersymmetric weight function, wa (y), are concentrated only in the regions near the left and right wells of the potential. There are no quadrature points in the middle well in the vicinity of y = 0. To redistribute some quadrature points in this region, a Gaussian function is added so that the weight function becomes wc (y) = wa (y) + exp(−y 2 /2ǫ) [15], with the result that the quadrature points are redistributed from the outermost potential wells to the well about the origin. We show the rate of convergence of the eigenvalues for the quartic bistable system in Tables 2.1 – 2.3 and in Figs. 2.4 and 2.5 for ǫ = 0.1, 0.01, and 0.001. The computational domain is [−xmax , xmax ]. A portion of the results in Tables 2.1 and 2.2 for QDMa appeared in [8], but here we closely follow the convergence of the triplet, λ3 , λ4 and λ5 , and emphasize the rate of convergence in terms of the error ǫm . The convergence rates of the eigenvalues calculated with the QDM and the SCM for ǫ = 0.1 are similar, as shown in Table 2.1. For this value of ǫ, λ1 is not too small and λ3 , λ4 and λ5 are well separated. As shown in Appendix B, the calculation of the eigenvalues associated with the eigenfunctions of even and odd symmetries can be uncoupled and calculated separately. So the values of N for wa (y) in Table 2.1 can be considered effectively as N/2 and the convergence is analogous to that shown for wb (y). It is surprising that λ2 converges as quickly as shown in the table since with wb (y) the vanishing of the eigenfunction at the origin is not consistent with the actual boundary condition. In Table 2.1, we show that the value of λ1 is converged to 7 significant figures for the QDM with N = 26 for wa (y) and N = 12 for wb (y) whereas we find that N = 32 is required for the SCM to get comparable accuracy. The manner in which Wei [68] compares the convergence (see Table II of his paper) is inappropriate. He lists the converged values of a number of selected eigenvalues and lists one value of N for all. He incorrectly suggests that N = 60 is required to get the specified precision for λ1 with the QDM. For ǫ = 0.01, the calculation of the eigenvalues becomes more difficult owing to the smallness of λ1 and the near degeneracy of the eigenvalues 44 Chapter 2. Spectral convergence of quadrature discretization method Table 2.1: Convergence of the eigenvalues for the quartic potential, ǫ = 0.1. The errors ǫm are given by Eq. (2.13). N λ1 × 102 λ2 λ3 12 18 24 26 28 30 32 3.365796 3.354699 3.354531 3.354529 0.939670 0.927646 0.927376 0.927373 0.927372 0.927372 3 6 9 12 15 12 18 24 28 30 32 a b c 3.367516 3.355525 3.354543 3.354529 6.825052 3.427580 3.355037 3.354524 3.354531 3.354529 1.007699 0.927402 0.927372 0.957566 0.926484 0.927361 0.927372 λ4 ǫ1 ǫ2 ǫ3 ǫ4 ǫ5 1.705240 1.680910 1.680273 1.680266 1.680264 λ5 QDMa 2.743490 4.004819 2.601742 3.747467 2.595936 3.734290 2.595849 3.734062 2.595827 3.734003 2.595822 3.733989 2.595821 3.733986 -2.5 -4.3 -6.2 -6.9 -7.6 -8.3 -9.0 -1.9 -3.5 -5.4 -6.0 -6.7 -7.4 -8.1 -1.8 -3.4 -5.2 -5.9 -6.5 -7.2 -7.9 -1.2 -2.6 -4.4 -5.0 -5.6 -6.2 -6.9 -1.1 -2.4 -4.1 -4.7 -5.3 -5.9 -6.6 1.764537 1.684458 1.680328 1.680264 QDMb 3.237921 4.994304 2.631757 3.807057 2.596849 3.735433 2.595835 3.733999 2.595820 3.733985 -2.4 -3.5 -5.4 -7.5 -9.8 -1.1 -2.9 - 4.5 -6.4 -8.5 -1.3 -2.6 -4.4 -6.5 -8.7 -0.6 -1.9 -3.4 -5.3 -7.3 -0.5 -1.7 -3.4 -5.4 -7.6 1.817707 1.679049 1.680236 1.680265 1.680264 SCMc 2.974622 4.093426 2.597213 3.750052 2.595783 3.734026 2.595824 3.733991 2.595821 3.733986 2.595820 3.733985 0.0 -1.7 -3.8 -5.9 -6.1 -7.8 -1.5 -3.0 -4.9 -6.4 -7.6 -8.1 -1.1 -3.1 -4.8 -6.1 -7.7 -7.8 -0.8 -3.3 -4.8 -5.9 -7.6 -7.5 -1.0 -2.4 -5.0 -5.8 -6.7 -7.2 wa (y) = exp[−(y 4 /4 − y 2 /2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞, ∞). For the even eigenstates, wbe (y) = wa (y), y ∈ [0, ∞). For the odd eigenstates, wbo (y) = y 2 wa (y), y ∈ [0, ∞). xmax = 2.2. 45 Chapter 2. Spectral convergence of quadrature discretization method Table 2.2: Convergence of the eigenvalues for the quartic potential, ǫ = 0.01. N λ1 λ2 λ3 λ4 λ5 ǫ1 ǫ2 ǫ3 ǫ4 ǫ5 3.9 0.7 -0.8 -2.4 -2.8 -4.5 -3.4 -2.5 -0.4 -1.6 -3.2 -5.3 -7.6 -8.9 -3.1 -3.2 -3.2 -3.2 -3.4 -5.0 -7.3 -8.5 -4.2 -6.9 -5.6 -6.6 -7.9 -9.8 -12.1 -13.0 -0.4 -1.3 -2.9 -5.0 -7.3 -8.5 -1.4 -3.0 -2.5 -2.6 -2.6 -2.6 -0.5 -1.6 -2.9 -4.8 -6.5 -7.4 -3.2 -3.2 -3.2 -5.4 -6.6 -7.7 -5.7 -6.8 -7.5 -9.6 -11.2 -12.6 -0.9 -1.6 -2.3 -5.4 -6.6 -7.7 11.1 8.7 6.9 3.5 -0.5 -2.6 -0.9 -4.3 -9.6 -13.7 -13.6 -13.6 -0.1 -1.4 -3.6 -6.8 -10.8 -12.6 -0.1 -1.5 -3.4 -6.5 -10.5 -12.2 -0.8 -3.0 -3.7 -6.8 -10.8 -12.5 QDMa 12 24 36 48 60 72 84 90 12 15 18 24 27 30 12 24 36 48 60 66 5.0833(-8) 3.6651(-11) 7.0354(-12) 6.1809(-12) 6.15499(-12) 6.15466(-12) 6.15465(-12) 6.4259(-12) 6.1656(-12) 6.1405(-12) 6.1424(-12) 6.1436(-12) 6.1427(-12) 7.4085(-1) 3.3865(-3) -4.8093(-5) 1.7533(-8) 7.9956(-12) 6.1435(-12) 1.388230 0.994289 0.968472 0.967870 0.967865 0.967864 1.256087 0.990778 0.969092 0.967879 0.967865 0.967864 1.076821 0.967915 0.967864 1.866176 1.865757 1.865752 1.865735 1.865337 1.864560 1.864542 1.865861 1.865753 1.865758 1.865754 2.664871 1.956370 1.869329 1.866993 1.866975 1.865747 1.865720 1.865601 1.864549 1.864542 QDMb 1.865757 2.113341 1.865754 1.913825 1.875631 1.866982 1.866975 3.336192 1.931199 1.864066 1.864542 SCMc 3.321671 1.574892 1.930972 1.865051 1.864927 1.866629 1.865754 1.866975 a wa (y) = exp[−(y 4 /4 − y 2 /2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞, ∞); λ1 is calculated by diagonalizing the matrix Lodd of dimension N/2 × N/2 defined by Eqs. (B.3) and (B.6). b For the even eigenstates we (y) = w (y), y ∈ [0, ∞). For the odd eigenstates, wo (y) = y 2 w (y), y ∈ a a b b [0, ∞). c x max = 1.5. 46 Chapter 2. Spectral convergence of quadrature discretization method Table 2.3: 0.001. Convergence of the eigenvalues for the quartic potential, ǫ = N λ2 λ3 6 12 18 24 30 0.9980526 0.9969809 0.9969817 2.0000470 1.9878205 1.9878873 1.9878896 12 24 36 48 60 72 84 96 108 120 a b 3.4140030 1.2875835 1.0279928 0.9984717 0.9970079 0.9969819 0.9969817 21.7611343 8.5716393 4.1568953 2.5344712 2.0928734 1.9995592 1.9884160 1.9878838 1.9878884 1.9878896 λ4 λ5 QDMa 2.0200067 2.0694590 1.9880010 1.9881554 1.9878903 1.9878937 1.9878896 1.9878893 1.9878896 ǫ2 ǫ3 ǫ4 ǫ5 -3.0 -6.1 -7.1 -8.3 -8.7 -2.2 -4.5 -5.9 -8.2 -9.8 -1.8 -4.3 -6.5 -8.5 -8.9 -1.4 -3.9 -5.7 -6.8 -7.7 SCMb 21.7517512 3.4979300 8.5704015 1.6476678 4.1569000 1.7776398 2.5345542 1.9649680 2.0928574 1.9872886 1.9995565 1.9878842 1.9884156 1.9878896 1.9878838 1.9878884 1.9878896 0.4 -0.5 -1.5 -2.8 -4.6 -6.8 -9.4 -12.3 -13.2 -13.9 1.0 0.5 0.0 -0.6 -1.3 -2.2 -3.6 -5.5 -6.2 -7.8 1.0 0.5 0.0 -0.6 -1.3 -2.2 -3.6 -5.5 -6.2 -7.8 -0.1 -0.8 -1.0 -1.9 -3.5 -5.6 -8.0 -10.9 -13.8 -13.5 wc (y) = exp[−(y 4 /4 − y 2 /2)/ǫ]/ exp(1/4ǫ) + exp(−y 2 /2ǫ), y ∈ (−∞, ∞). xmax = 1.2. 47 Chapter 2. Spectral convergence of quadrature discretization method 4 4 A B QDM (a) 0 SCM 0 QDM (a) QDM (b) ǫ3 -4 ǫ2 -4 -8 -8 SCM -12 -16 0 -12 40 80 -16 120 QDM (b) 0 40 N 80 120 N 4 4 C D 0 0 QDM (a) SCM -4 ǫ5 ǫ4 -4 QDM (a) -8 -8 SCM QDM (b) -12 -16 -12 0 40 80 N 120 -16 QDM (b) 0 40 80 120 N Figure 2.4: Variation of ǫm for (A) m = 2, (B) m = 3, (C) m = 4, and (D) m = 5 versus the number of quadrature points N for the bistable system with ǫ = 0.01. For QDM(a), wa (y) = exp[−(y 4 /4 − y 2 /2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞, ∞). For the QDM(b), wbe (y) = wa (y), y ∈ [0, ∞) for the even eigenstates and wbo (y) = y 2 wa (y), y ∈ [0, ∞) for the odd eigenstates. For the SCM, xmax = 1.5. 48 Chapter 2. Spectral convergence of quadrature discretization method 4 4 A B SCM QDM (b) 0 0 QDM (a) -4 -4 QDM (c) ǫ3 ǫ2 SCM -8 QDM (b) -8 QDM (c) -12 -16 -12 0 40 80 -16 120 0 40 N 80 120 N 4 4 C D SCM 0 QDM (b) 0 QDM (a) -4 -4 -8 -8 QDM (b) -12 -16 SCM ǫ5 ǫ4 QDM (c) QDM (c) -12 0 40 80 N 120 -16 0 40 80 120 N Figure 2.5: Variation of ǫm for (A) m = 2, (B) m = 3, (C) m = 4, and (D) m = 5 versus the number of quadrature points N for the bistable system with ǫ = 0.001. For QDM weight functions, see the caption for Fig. 2.4. For the SCM, xmax = 1.2. 49 Chapter 2. Spectral convergence of quadrature discretization method λ3 , λ4 and λ5 . In the limit ǫ → 0, λ2 → 1 and λ3 , λ4 and λ5 → 2. The variation of some of the lower order eigenfunctions for this quartic bistable system were shown in Fig. 1 in Ref. 63. Owing to the symmetry of the problem [8, 34, 35], the eigenfunctions are either concentrated in the two outside wells near y = ±1, the central well at y = 0 or in all three. The eigenfunction associated with λ2 is localized in the middle well whereas the eigenfunction associated with λ4 is concentrated in the two outer wells and those of λ3 and λ5 are concentrated in all three wells. The results for the convergence of these eigenvalues are summarized in Table 2.2 and in Fig. 2.4. The eigenvalues converge very differently and no one value of N suffices for all the eigenvalues. As mentioned previously, the values of N in the table for wa (y) are effectively N/2 because of the separability of the eigenvalue problem into odd and even eigenfunctions. In Table 2.2, the most noticeable features of the results is that the convergence of λ1 is far more rapid with the QDM than for the SCM. The λ1 value is close to machine accuracy and reliable values cannot be calculated in double precision for ǫ 0.008. The values reported for QDMa were determined by diagonalizing the matrix of dimension N/2 × N/2 in Eq. (B.3) in Appendix B for the odd eigenstates. While the lower order estimates for this eigenvalue with the QDM are close to the converged value, the lower order estimates with the SCM are very poor. As can be seen from the results in the table, we can calculate λ1 with wbo (y) to three significant figures with N = 18 whereas the estimates with the SCM are unconverged until N = 66. The value converged to 8 significant figures is λ1 = 6.1546530 × 10−12 with 45 nonclassical basis functions based on wa (y). Wei [68] reports λ1 = 1.278 ×10−10 in Table 2.3 of his paper with N = 52 which is 2 orders of magnitude too large. In the text, Wei incorrectly compares a value of λ1 = 1.28 ×10−8 with λ1 = 6.077 ×10−8 reported by Shizgal and Chen [63] (their Table I) which was clearly stated as being unconverged. The footnote to that table mentions that a value of λ1 = 6.16 ×10−12 was obtained with N = 80. The value of λ1 = 6.154650×10−12 in Table III in Wei’s paper [68] is taken from Table III of Blackmore snd Shizgal [8] and is for N = 70 rather than N = 100 as stated by Wei. The results with the SCM (for the Schr¨ odinger equation) depend on the choice of xmax . The values of xmax in the tables and figure captions are chosen as a compromise between a fast convergence and the attainability of a converged eigenvalue close to machine precision. Figure 2.4 summarizes the spectral convergence of these lower order eigenvalues for ǫ = 0.01. In each case the half-range weight function wb (y) (open squares) outperforms the other choices, except for the unattainability 50 Chapter 2. Spectral convergence of quadrature discretization method of machine accuracy for λ2 . The initial convergence is comparable with the SCM but, owing to the imposition of the incorrect boundary condition at the origin, machine accuracy is not achieved. The convergence with wb (y) is very rapid and faster than the SCM for the other eigenvalues. For λ3 , the initial estimates with wa (y) are good but the convergence is very flat until about N = 60. The SCM performs better than the QDM with wa (y) for this eigenvalue. However, the weight function wbo (y) provides rapid convergence for this case because there is a node for the corresponding eigenfunction at the origin. This is not the case for the eigenfunction associated with λ2 . The SCM does not perform as well although it achieves machine accuracy for smaller N than the QDM with wa (y). For λ4 , the QDM with wa (y) provides initially good estimates but then achieves machine accuracy slowly. The QDM with wa (y) doe not provide a rapid convergence for λ5 and SCM performs better, but the QDM with wbo (y) provides the best convergence. The numerical results shown in Table 2.3 complement the results in Fig. 2.4. Table 2.3 and Fig. 2.5 show a comparison of the convergence of the lower order eigenvalues with ǫ = 0.001, for the QDM, with the weight function wc (y) and wb (y) in comparison with the results for the SCM. Wei [68] did not report results for this value of ǫ. The quadrature points generated by the supersymmetric weight function, wa (y), are localized in the left and right wells. Since there are no quadrature points near the middle well, the eigenvalues associated with the eigenfunctions ψ2 (y) and ψ5 (y) which have major characteristics in the vicinity of y = 0 [63] are not found. There is a somewhat similar behaviour for the half-range polynomials defined by wb (y) for which eigenvalues λ2 and λ5 are calculated but only with N greater than about 80, and the ensuing convergence is slow. To calculate these eigenvalues efficiently, one requires a different weight function, wc (y), with an additional Gaussian centered at y = 0. The convergence with this weight function shown in Table 2.3 and in Fig. 2.5 is very rapid. With wc (y) the number of quadrature points in the vicinity of y = 0 increases significantly at the expense of some points in the outer wells. Thus the eigenvalues converge faster than with wa (y) or with wb (y). The convergence of λ3 and λ4 is extremely rapid for the QDM and all three weight functions. The SCM is based on a uniform grid and does not resolve well the details of the localized eigenfunctions and the convergence rate is not comparable to the rate with the QDM. 51 Chapter 2. Spectral convergence of quadrature discretization method 2.7 Interpretation of the numerical results The main objective of the present chapter is to compare the rates of convergence of the solutions of the Fokker-Planck and Shr¨ odinger equations obtained with the QDM and the SCM. These results have been compared and discussed in the previous section. In this section, we are interested in the analysis of the spectral convergence of the QDM as shown in Figs. 2.1, 2.2, 2.4, and 2.5. Spectral convergence as mentioned in the introduction refers to the exponential decrease of the errors in the eigenvalues versus the size of the basis set or the number of quadrature points, N , that is, ) exact |λ(N | ≈ 10−qN . m − λm (2.14) (N ) Spectral convergence has been demonstrated in the results for which log |λm − λexact | varies linearly with N . The convergence is from above the exact m eigenvalue as shown in the tables consistent with a Rayleigh-Ritz variational approach [30]. We consider here an analysis based on the polynomial basis set representation of the operators. It is important to keep in mind the definitions of the matrix representatives of the Fokker-Planck operator and the Hamiltonian as given by Eq. (A.3) and (C.1) in Appendices A and C, respectively. For the Fokker-Planck operator for the electron relaxation problem and the bistable system, the matrix representatives can be evaluated analytically as shown in Eqs. (A.8) and (B.3) in Appendices A and B. There is a unitary transformation between the basis set representation, Eq. (A.8), and the pseudospectral representation, Eq. (2.2). Although the pseudospectral representation (equivalently collocation or DVR) is an efficient numerical methodology, an analysis of the numerical results from the basis set representation of the operators is, in fact, the more fundamental and insightful approach. In this section, we provide some insights as to why the QDM provides in many instances a very rapid convergence. We cannot hope to provide in this chapter a detailed analysis for all the examples studied. The numerical analysis of the results is very difficult for the nonclassical polynomial basis functions used here as there is very little known concerning their mathematical properties. We focus our attention on the electron relaxation problem from both the Fokker-Planck and Schr¨ odinger equations. In Appendix A, we show, for the electron relaxation problem, the details of the derivation of the tridiagonal matrix representative of the Fokker-Planck operator, Eq. (A.8), in the speed polynomial basis set based on w(x) = x2 exp(−x2 ). The numerical values of the lower order matrix elements are also shown in Eq. 52 Chapter 2. Spectral convergence of quadrature discretization method (A.9) and it is clear that the matrix is diagonally dominant. Much of the analysis of the spectral convergence presented here is based on this tridiagonal matrix. Appendix B provides the details of the calculation of the matrix representative of the Fokker-Planck operator for the bistable system with the bimode weight function wa (x) = P0 (x)/ exp(1/4ǫ) and P0 (x) = exp[−(x4 /4ǫ) + (x2 /2ǫ)]. We evaluate analytically the matrix representative of this FokkerPlanck operator in terms of the recurrence coefficients in the three-term recurrence relations for these polynomials. From the structure of the matrix shown in Eq. (B.3) in Appendix B the N × N matrix, Eq. (B.4), can be split into two N/2 × N/2 matrices, Eqs. (B.5) and (B.6), that are diagonally dominant. Diagonalization of these matrices gives separately the eigenfunctions of odd and even symmetries and the corresponding eigenvalues. One might expect that the convergence will be rapid since these matrices are diagonally dominant with diagonal elements close to the converged values of the eigenvalues. While this is the case for some of the eigenvalues, the convergence of the degenate eigenfunctions for small ǫ of the order of 0.01 is particularly slow with the weight function wa (y) = P0 (y). For this reason, two sets of half-range polynomials based on wbe (y) = wa (y), y ∈ [0, ∞), and wbo (y) = y 2 wa (y), y ∈ [0, ∞), were employed. These basis sets yield the diagonally dominant matrices in Eqs. (C.4) and (C.5) in Appendix C. The first basis set is for the even eigenfunctions and the second basis set is for the odd eigenfunctions. The rapid convergence of the eigenvalues for these basis sets is shown in Table 2.2. The convergence of the nearly degenerate eigenvalues λ3 and λ5 is faster with this basis set than the one defined with wa (y). All of the examples considered here lead to matrix representatives that are diagonally dominant, and we would expect a rapid convergence of the eigenfunctions and eigenvalues. For the Fokker-Planck equation and a diffusion coefficient, B(x), which is a polynomial, the matrix representative is banded. For the electron relaxation problem, B(x) = x and the matrix representative is tridiagonal, Eqs. (A.8) and (A.9). If the diagonal and offdiagonal matrix elements of this tridiagonal matrix are denoted by an and cn respectively, the eigenvalues for order N are determined as the roots of the characteristic polynomial, PN (λ) = 0, where the polynomials Pn (λ) can be determined by iteration of a three-term recurrence relation, that is, Pn (λ) = (λ − an )Pn−1 − c2n−1 Pn−2 with P0 = 1 and P1 = λ − a1 . We focus on λ1 and find that it can be 53 Chapter 2. Spectral convergence of quadrature discretization method 0 0 A B m=5 m=5 -4 -8 m=1 m=3 -12 -16 log |dn(m)| log |dn(m)| -4 0 10 -8 m=1 m=3 -12 20 30 -16 0 10 n 20 30 n (m) Figure 2.6: Variation of the expansion coefficients dn versus n of the mth eigenfunction of the Fokker-Planck operator (A) and the Hamiltonian in the Schr¨ odinger equation (B) for the electron relaxation problem; m = 1, 3, and 5. 54 Chapter 2. Spectral convergence of quadrature discretization method represented as a continued fraction [75] of the form c21 λ1 = a1 + λ1 − a2 − . c22 λ1 − a3 − (2.15) c23 λ1 − a4 − ··· The other eigenvalues do not have a similar continued fraction representations except for λN . The N th approximation to λ1 is obtained by truncating (N ) the continued fraction with the term in c2N +1 . The numerical value of λ1 (1) can be determined by starting with the initial estimate λ1 = a1 and calculating successive approximations by substituting for λ1 on the right hand (N ) side of Eq. (2.15). The rate of convergence of this iterative process to λ1 (N ) is analogous to the convergence of λ1 to λ1 versus N . With Eq. (2.15), we obtain the approximation (2) (3) λ1 − λ1 ≈ c21 c22 , (λ1 − a3 )(λ1 − a2 )2 (2) where on the right hand side of the continued fraction we have set λ1 (3) λ1 = λ1 . Similarly we can calculate (3) (4) λ1 − λ1 ≈ ≈ c21 c22 c23 . (λ1 − a4 )(λ1 − a3 )2 (λ1 − a2 )2 With the spectral convergence given by Eq. (2.14), we can estimate q with the ratio (3) (4) λ1 − λ1 c23 = = 10−q . (2.16) (2) (3) (λ1 − a4 )(λ1 − a3 ) λ1 − λ1 With Eq. (2.16) and the value of λexact = 4.68340, we find that q = 1.05 is 1 in very good agreement with the slope of the linear portion of the graph in Fig. 2.1, which is q = 0.990. One can show by induction that the rate of convergence depends on the quantity limN →∞ c2N −1 /|(λ1 − aN )(λ − aN −1 )|, which we estimate to be 0.088 and with this we get q = 1.06. We have demonstrated explicitly spectral convergence albeit for such tridiagonal matrices was discussed in the appendix of the paper by Meyer et al [43]. These authors prove convergence if c2N −1 /|(λ1 − aN )(λ − aN −1 )| ≤ 1.5. The rapid convergence of the QDM is related to the small values of the off-diagonal elements relative to the diagonal elements [44]. In Fig. 2.6A, we show the 55 Chapter 2. Spectral convergence of quadrature discretization method (m) variation of the nth elements in the mth eigenvectors, dn , which serve as the expansion coefficients of the mth eigenfunction, with m = 1, 3, and 5. The very rapid linear decrease of the expansion coefficients shown in the figure is a demonstration of spectral convergence. The expansion coefficients (m) for the eigenvectors, dn , can also be represented as continued fractions. For our purposes here, we consider a simple estimate of the slope of the graph in Fig. 2.6A for m = 1, and we have that the eigenvector expansion coefficients satisfy a three-term recurrence relation given by (1) (1) exact (1) cn−1 dn−1 + an d(1) dn . n + cn dn+1 = λ1 (1) If the coefficient dn+1 is neglected, we have the approximation (1) dn (1) dn−1 = cn−1 exact λ1 − an . The slope of the graph in Fig. 2.6 is thus log(cn−1 /|λexact − an |), the asymp1 totic value of which we estimate to be −0.53, in good agreement with the slope in Fig. 2.6 (for m = 1) which is −0.52. In Appendix C, we show the numerical values of the lower order matrix elements of the full matrix representative of the Schr¨ odinger equation for the same electron relaxation problem with the weight function w(y) = y 5 exp(−y 4 /16), Eqs. (C.2) and (C.3). Although the matrix is full, the off-diagonal elements are much smaller than the diagonal elements and the convergence of the eigenvalues is actually very similar to that for the tridiagonal matrix. The two matrices give the same eigenvalues in the limit N → ∞. A similar demonstration of spectral convergence is shown in Fig. 2.6B for the corresponding eigenfunctions of the Schr¨ odinger equation for the electron relaxation problem in Section 2.2. There have been numerous discussions in the literature regarding the “phase space” coverage of a chosen basis set [13, 23, 32, 33, 36, 45, 47] relative to the phase space extent of the classical Hamiltonian corresponding to the quantum Hamiltonian operator. These analyses have been considered for the most part for Fourier basis functions or harmonic oscillator eigenfunctions (Hermite polynomials). In this chapter, there is no defining Hamiltonian problem for the basis sets used and thus a direct phase space analysis cannot be considered. We can mention that for the basis sets used here the coverage is exact for the ground state and thus it should be fairly 56 Chapter 2. Spectral convergence of quadrature discretization method good for the excited states. That the representation of the operators is in some sense “optimal” is seen from the diagonally dominant matrices that result and the rapid spectral convergence that is obtained. One method for easily optimizing basis sets is to scale the independent variable [1, 2, 25, 61] or equivalently the quadrature points and weights. Scaling of the Laguerre polynomials in speed as discussed by Risken and Voigtlaender [49] is essential to get the spectral convergence shown in Fig. 2.1. A scaling of the Fourier-Hermite solutions of the Vlasov equation [25] is required to obtain convergent solutions. Some researchers have developed criterion for the best choice of scaling parameters [1, 2]. Further work concerning the choice of basis functions in this chapter and in what sense they are optimal is ongoing. 2.8 Summary We have demonstrated in this chapter for three model problems the spectral convergence of the QDM, a collocation (pseudospectral) methods based on nonclassical polynomial basis functions. With the appropriate choice of the weight function, the convergence of the eigenvalues of the Schr¨ odinger equation can be extremely rapid. For the potential in the Schr¨ odinger equation corresponding to the electron relaxation problem as well as for the Morse potential for I2 the QDM [62, 63] outperforms the SCM proposed by Wei [68]. For the quartic bistable problem, the nature of the eigenfunctions in the triple well potential parameterized by a diffusion coefficient ǫ in the corresponding Fokker-Planck equation varies considerably so that no one weight function in the QDM can provide a rapid convergence for all the eigenvalues. For ǫ > 0.1, the weight function corresponding to the equilibrium distribution, P0 (y), provides rapid convergence for a large number of eigenvalues. For intermediate values, ǫ ≈ 0.01, a half-range weight function, wb (y), y ∈ [0, ∞) provides a rapid convergence of the eigenvalues for which the corresponding eigenfunctions possess a node at the origin. For ǫ ≈ 0.001 a weight function with a Gaussian about the origin added to P0 (y) provides very rapid convergence. We have also provided an analysis of the spectral convergence of λ1 and the associated eigenfunction for the electron relaxation problem. 57 Bibliography [1] P. Amore. A variational sinc collocation method for strong coupling problems. J. Phys. A: Math. Gen., 39:L349, 2006. [2] P. Amore, A. Aranda, F. M. Ferm´ andez, and H. Jones. A new approximation method for time-dependent problems in quantum mechanics. Phys. Lett. A, 340:87, 2005. [3] K. Anderson and K. E. Shuler. On the relaxation of the hard-sphere Rayleigh and Lorentz gas. J. Chem. Phys., 40:633, 1964. [4] G. G. Balint-Kurti and P. Pulay. A new grid-based method for the direct computation of excited molecular vibrational states: test application to formaldehyde. Theochem: J. Mol. Struct., 341:1, 1995. [5] D. Baye, M. Hesse, and M. Vincke. The unexplained accuracy of the Langrange-mesh method. Phys. Rev. E, 65:026701, 2002. [6] D. Baye and M. Vincke. Lagrange meshes from nonclassical orthogonal polynomials. Phys. Rev. E, 59:7195, 1999. [7] R. Blackmore and B. Shizgal. A solution of Kramers equation for the isomerization of n-butane in CCl4 . J. Chem. Phys., 83:2934, 1985. [8] R. Blackmore and B. D. Shizgal. Discrete-ordinate method of solution of Fokker-Planck equations with nonlinear coefficients. Phys. Rev. A, 31:1855, 1985. [9] J. P. Boyd. Chebyshev and Fourier Spectral Methods. Dover, New York, 2000. [10] J. P. Boyd. A proof that the discrete singular convolution (DSC) / Lagrange-distributed approximating function (LDAF) method is inferior to high order finite differences. J. Comput. Phys., 214:538, 2006. [11] M. Braun, S. A. Sofianos, D. G. Papageorgiou, and I. E. Lagaris. An efficient Chebyshev-Lanczos method for obtaining eigensolutions of the Schr¨ odinger equation on a grid. J. Comput. Phys., 126:315, 1996. 58 Bibliography [12] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, New York, 1988. [13] M. C. Cargo and R. G. Littlejohn. Phase space deformation and basis set optimization. Phys. Rev. E, 65:026703, 2002. [14] S. Chandrasekhar. Radiative Transfer. Dover, New York, 1949. [15] H. Chen and B. D. Shizgal. A spectral solution of the Sturm-Liouville equation: comparison of classical and nonclassical basis set. J. Comput. Appl. Math., 136:17, 2001. [16] H. Chen, Y. Su, and B. D. Shizgal. A direct spectral collocation poisson solver in polar and cylindrical coordinates. J. Comput. Phys., 160:453, 2000. [17] D. T. Colbert and W. H. Miller. A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method. J. Chem. Phys., 96:1982, 1992. [18] A. Comtet, A. D. Bandrauk, and D. K. Campbell. Exactness of semiclassical bound state energies for supersymmetric quantum mechanics. Phys. Lett. B, 150:159, 1985. [19] O. Dulieu, R. Kosloff, F. Masnou-Seeuws, and G. Pichler. Quasibound states in long-range alkali dimers: grid method calculations. J. Chem. Phys., 107:10633, 1997. [20] J. Dunkel, W. Ebeling, L. Schimansky-Geier, and P. Hangii. Kramers problem in evolutionary strategies. Phys. Rev. E, 67:061118, 2003. [21] J. Echave and D. C. Clary. Potential optimized discrete variable representation. Chem. Phys. Lett., 190:225, 1992. [22] K. S. Fa and E. K. Lanzi. Power law diffusion coefficient and anomalous diffusion: analysis of solutions and first passage time. Phys. Rev. E, 67:061105, 2003. [23] E. Fattal, R. Baer, and R. Kosloff. Phase space approach for optimizing grid representations: the mapped Fourier method. Phys. Rev. E, 53:1217, 1996. [24] R. D. M. Garcia. The application of nonclassical orthogonal polynomials in particle transport theory. Prog. Nucl. Energy, 35:249, 1999. 59 Bibliography [25] L. Gibelli and B. D. Shizgal. Spectral convergence of the Hermite basis function solution of the Vlasov equation: the free-streaming term. J. Comput. Phys., 219:477, 2006. [26] D. Gottlieb and S. A. Orszag. Numerical Analysis of Spectral Methods: Theory and Applications. SIAM, Philadelphia, 1977. [27] I. P. Hamilton and J. C. Light. On distributed Gaussian bases for simple model multidimensional vibrational problems. J. Chem. Phys., 84:306, 1986. [28] P. H¨ anggi, P. Talkner, and M. Borkovec. Reaction-rate theory: fifty years after Kramers. Rev. Mod. Phys., 62:251, 1990. [29] Y. He, Z. Cao, and Q. Shen. Bound-state spectra for supersymmetric quantum mechanics. Phys. Lett. A, 326:315, 2004. [30] R. N. Hill. Rates of convergence and error estimation formulas for the Rayleigh-Ritz variational method. J. Chem. Phys., 83:1173, 1985. [31] K. D. Knierem, M. Waldman, and E. A. Mason. Moment theory of electron thermalization in gases. J. Chem. Phys., 77:943, 1982. [32] V. Kokoouline, O. Dulieu, R. Kosloff, and F. Masnou-Seeuws. Mapped Fourier methods for long-range molecules: application to perturbations in the Rb2 (0+ u ) photoassociation spectrum. J. Chem. Phys., 110:9865, 1999. [33] R. Kosloff. In J. Zang and R. E. Wyatt, editors, Dynamics of Molecular and Chemical Reactions. Decker, New York, 1996. [34] R. S. Larson and M. Kostin. Kramers’s theory of chemical kinetics: eigenvalue and eigenfunction analysis. J. Chem. Phys., 69:4821, 1978. [35] R. S. Larson and M. Kostin. Friction and velocity in Kramers’ theory of chemical kinetics. J. Chem. Phys., 72:1392, 1980. [36] H. S. Lee and J. C. Light. Molecular vibrations: iterative solution with energy selected bases. J. Chem. Phys., 118:3458, 2003. [37] K. Leung, B. D. Shizgal, and H. Chen. The quadrature discretization method (QDM) in comparison with other numerical method of solution of the Fokker-Planck equation for electron thermalization. J. Math. Chem., 24:291, 1998. 60 Bibliography [38] J. C. Light and T. Carrington. Discrete variable representations and their utilization. Adv. Chem. Phys., 114:263, 2000. [39] J. C. Light, I. P. Hamilton, and J. V. Lill. Generalized discrete variable approximation in quantum mechanics. J. Chem. Phys., 82:1400, 1985. [40] M. J. Lindenfeld and B. Shizgal. Matrix elements of the Boltzmann collision operator for gas mixtures. Chem. Phys., 41:81, 1979. [41] L. Lombardini and B. Poirier. Rovibrational spectroscopy calculations of neon dimer using a phase space truncated Weyl-Heisenberg wavelet basis. J. Chem. Phys., 124:144107, 2006. [42] D. A. Mazziotti. Spectral difference methods for solving the differential equations of chemical physics. J. Chem. Phys., 117:2455, 2002. [43] H. D. Meyer, J. Hor´acek, and L. S. Cederbaum. Schwinger and anomaly-free Kohn variational principles and a generalized Lanczos algorithm for nonsymmetric operators. Phys. Rev. A, 43:3587, 1991. [44] H. D. Meyer and S. Pal. A band-Lanczos method for computing matrix elements of a resolvent. J. Chem. Phys., 91:6195, 1989. [45] J. Montgomery and B. Poirier. Eigenspectra calculations using Cartesian coordinates and a rotational symmetry adapted Lanczos method. J. Chem. Phys., 119:6609, 2003. [46] R. Peyret. Spectral Methods for Incompressible Viscous Flow. SpringerVerlag, New York, 2002. [47] B. Poirier and J. C. Light. Phase space optimization of quantum representations: direct-product basis sets. J. Chem. Phys., 111:4869, 1999. [48] H. Risken. The Fokker-Planck Equation: Methods of Solution and Application. Springer-Verlag, New York, 1984. [49] H. Risken and K. Voigtlaender. Solutions of the Fokker-Planck equation describing the thermalization of neutrons in a heavy gas moderator. Z. Phys. B: Condens. Matter, 54:253, 1984. [50] C. Schwartz. High-accuracy approximation techniques for analytic functions. J. Math. Phys., 26:411, 1985. [51] C. S. Shapiro and N. Corngold. Approach to equilibrium of a neutron gas. Phys. Rev., 137:A1686, 1965. 61 Bibliography [52] B. Shizgal. Eigenvalues of the Lorentz Fokker-Planck equation. J. Chem. Phys., 70:1948, 1979. [53] B. Shizgal. A Gaussian quadrature procedure for use in the solution of the Boltzmann equation and related problems. J. Comput. Phys., 41:309, 1981. [54] B. Shizgal. Nonequilibrium time dependent theory of hot atom reactions. III. Comparison with EstrupVWolfgang theory. J. Chem. Phys., 74:1401, 1981. [55] B. Shizgal. Discrete versus continuum relaxation modes of a hard sphere gas. Can. J. Phys., 62:97, 1984. [56] B. Shizgal and R. Blackmore. Eigenvalues of the Boltzmann collision operator for binary gases: relaxation of anisotropic distributions. Chem. Phys., 77:417, 1983. [57] B. Shizgal and R. Blackmore. A discrete ordinate method of solution of linear boundary value and eigenvalue problems. J. Comput. Phys., 55:313, 1984. [58] B. Shizgal and R. Blackmore. A collisional kinetic theory of a plane parallel evaporating planetary atmosphere. Planet. Space Sci., 34:279, 1986. [59] B. Shizgal and J. M. Fitzpatrick. Matrix elements of the linear Boltzmann collision operator for systems of two components at different temperatures. Chem. Phys., 6:54, 1974. [60] B. Shizgal, J. M. Lindenfeld, and R. Reeves. Eigenvalues of the Boltzmann collision operator for binary gases: Mass dependence. Chem. Phys., 56:249, 1981. [61] B. Shizgal and D. R. A. McMahon. Electric field dependence of transient electron transport properties in rare-gas moderators. Phys. Rev. A, 32:3669, 1985. [62] B. D. Shizgal and H. Chen. The quadrature discretization method (QDM) in the solution of the Schr¨ odinger equation with nonclassical basis functions. J. Chem. Phys., 104:4137, 1996. [63] B. D. Shizgal and H. Chen. The quadrature discretization method in the solution of the Fokker-Planck equation with nonclassical basis function. J. Chem. Phys., 107:8051, 1997. 62 Bibliography [64] F. Stenger. Numerical Methods Based on Sinc and Analytic Functions. Springer-Verlag, New York, 1993. [65] V. Szalay. The generalized discrete variable representation. An optimal design. J. Chem. Phys., 105:6940, 1996. [66] V. Szalay, G. Czak´ o, A. Nagy, T. Furtenbacher, and A. Cs´ asz´ ar. On onedimensional discrete variable representations with general basis function. J. Chem. Phys., 119:10512, 2003. [67] L. A. Viehland. Velocity distribution function and transport coefficients of atomic ions in atomic gases by a Gram-Charlier approach. Chem. Phys., 179:71, 1994. [68] G. W. Wei. Discrete singular convolution for the solution of the FokkerPlanck equation. J. Chem. Phys., 110:8930, 1999. [69] G. W. Wei. A unified approach for the solution of the Fokker-Planck equation. J. Phys. A: Math. Gen., 33:4935, 2000. [70] G. W. Wei. Solving quantum eigenvalue problems by discrete singular convolution. J. Phys. B: At. Mol. Opt. Phys., 33:343, 2000. [71] G. W. Wei, D. S. Zhang, D. J. Kouri, and D. K. Hoffman. Lagrange distributed approximating functionals. Phys. Rev. Lett., 79:775, 1997. [72] H. Wei. Ghost levels and near variational forms of the discrete variable representation: application to H2 O. J. Chem. Phys., 106:6885, 1997. [73] H. Wei and T. Carrington. The discrete variable representation of a triatomic Hamiltonian in bond length-bond angle coordinates. J. Chem. Phys., 97:3029, 1992. [74] H. Wei and T. Carrington. Discrete variable representations of complicated kinetic energy operators. J. Chem. Phys., 101:1343, 1994. [75] H. A. Yamani and M. S. Abdelmonem. The analytic inversion of any finite symmetric tridiagonal matrix. J. Phys. A: Math. Gen., 30:2889, 1997. [76] H. H. Yang, B. R. Seymour, and B. D. Shizgal. A Chebyshev pseudospectral multi-domain method for steady flow past a cylinder up to Re = 150. Computers & Fluids, 23:829, 1994. 63 Bibliography [77] H. H. Yang and B. Shizgal. Chebyshev pseudospectral multi-domain technique for viscous flow calculation. Comp. Meth. Appl. Mech. Eng., 118:47, 1994. 64 Chapter 3 Pseudospectral methods of solution of the Schr¨ odinger equation 3.1 Introduction There has been an ongoing effort by numerous researchers to develop accurate and efficient algorithms for the calculation of the eigenvalues of the 1D Schr¨ odinger equation for several different potentials. Although these 1D calculations are not computationally intensive, their improvement will find useful application to multi-dimensional problems. Examples of this research endeavor are the papers by Simos and coworkers that are concerned with either higher order algebraic methods [42], and direct integrations with a Numerov [20], Runge-Kutta [1,21] or symplectic methods [30,31]. Other methods include a Ricatti-Pade approach [16], a Chebyshev-Lanczos method [6] and a Hill determinant method [50] to mention just a few. Other references to similar studies were provided in previous papers [9, 10]. Pseudospectral methods include the quadrature discretization method [9,39,40], the discrete variable representation [24], and the Lagrange mesh method [2–4]. Other pseudospectral methods based on classical bases have also been discussed by Ta¸seli and coworkers [44, 45]. The pseudospectral methods evaluate the eigenfunctions on grid of points which coincide with the quadrature points for the weight function chosen. The diagonalization of the discrete matrix representative of the Hamiltonian of dimension N 2 gives N eigenvalues of which a subset corresponds to the discrete eigenvalues for the problem. The spectral convergence of the eigenvalues refers to the exponential decrease of the error in the approximate A version of this chapter has been accepted for publication. Joseph Q. W. Lo and Bernie D. Shizgal. J. Math. Chem., available online at http://dx.doi.org/10.1007/s10910007-9341-8 65 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation eigenvalues versus N . Spectral convergence has been demonstrated recently by Lo and Shizgal [28] for several 1D problems. The present chapter is a continuation of these earlier studies but we here consider a shift and scaling of the quadrature points which we demonstrate can accelerate the convergence. Shizgal and coworkers [28, 35, 40] have demonstrated that rapid convergence of the eigenvalues of the Fokker-Planck equation can be obtained with a basis set defined by the equilibrium probability density as the weight function. The Fokker-Planck equation can be transformed to a Schr¨ odinger equation with a potential such that the Hamiltonian belongs to the class of supersymmetric quantum mechanics for which the ground state is known [39]. Chen and Shizgal [9] and Lo and Shizgal [28] obtained rapid convergence of the eigenvalues of the Hamiltonian operator when the square of the ground-state eigenfunction is used for the weight function. Occasionally modifications of these weight functions are still required to improve convergence [9, 28]. In most applications, such nonclassical basis sets were used. This method is referred to as the quadrature discretization method (QDM) since it is implemented as a pseudospectral (collocation) method. The QDM is generally based on quadrature points defined by nonclassical polynomials orthogonal with respect to a specific weight function. There has been a long history of the application of numerical pseudospectral methods [8, 32]. Shizgal [36] developed a pseudospectral method based on a Gaussian quadrature to accurately discretize the integral operator in the Boltzmann equation. By contrast, the calculation of the matrix elements of the Boltzmann collision operator in the basis set can lead to considerable round-off errors [27]. Shizgal and Blackmore [37] then applied this method to the solution of differential equations. The approach follows on previous works employing similar methods in neutron transport [17]. The QDM has been used to solve the time dependent Fokker-Planck equation [5, 40], the Poisson equation [10], the advection-diffusion equation [38], and the Schr¨ odinger equation [9, 28, 39]. An alternate method developed by Light and coworkers [24–26], which is referred to as the discrete variable representation (DVR), was originally based on the numerical evaluation of matrix elements of the potential in the Schr¨ odinger equation [12, 19]. If V (x) is the potential in a 1D Schr¨ odinger equation, then the matrix elements of V (x) can be determined with a Gaus- 66 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation sian quadrature, that is, N b Vmn = a φm (x)V (x)φn (x) dx ≈ ηk φm (xk )V (xk )φn (xk ), (3.1) k=1 where {φn } is a set of orthonormal basis functions, and {xk } and {ηk } are appropriate sets of quadrature points and weights, respectively [11]. The DVR is also a pseudospectral method in which the solution is evaluated at a set of grid points analogous to the QDM. The DVR was thus introduced by Light et al [25,26] and applied to several quantum problems [14,22,24,34,43]. We consider a set of polynomials Pn (x) orthogonal with respect to a weight function w(x), w(x)Pm (x)Pn (x) dx = δmn , R where R denotes the domain for x. With wk = ηk w(xk ), there is a unitary √ transformation, Tkn = wk Pn−1 (xk ), between the representation of a function in a basis set (that is the coefficients cn in ψ(x) = n cn φn (x) with φn (x) = w(x)Pn−1 (x)) and the discrete representation, ψ(xk ). The transformation of Vmn , Eq. (3.1), with T gives the diagonal representation of the potential in the discrete representation, that is, V (xk ). Further details are provided in [5, 39, 40] and Section 3.2. Thus the QDM was developed originally for kinetic theory problems whereas the DVR was introduced for the accurate evaluation of potential matrix elements. In Section 3.2, we briefly explain the variational principle and spectral convergence. Section 3.3 discusses the formulation of the pseudospectral method and the differences between QDM and DVR. Several applications are presented in Sections 3.4 to 3.6 with discussions of results in 3.7. 3.2 Spectral convergence The solution of the Schr¨ odinger equation 2 Hψ = − d2 ψ + V ψ = Eψ 2µ dx2 67 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation A B 0 0 ψ10, s = 0.14 -4 log|cn(m)| log|cn(m)| ψ10, s = 1 -8 -4 -8 ψ1, s = 1 ψ1, s = 0.125 -12 -12 0 10 20 30 40 0 10 n 20 30 n (m) Figure 3.1: Variation of the expansion coefficients cn versus n for ψm (x) of the I2 system for m = 1 and 10 for both QDM and Hermite. (A): QDM with b = 0; (B): Hermite with b = 0. 68 40 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation is generally considered with the expansion of the eigenfunctions in a basis set, N ψ(x) ≈ w(x) cn Pn−1 (x). (3.2) n=1 The matrix representation of the Hamiltonian in this basis set is Hmn = ( 2 /2µ)Kmn + Vmn , where Kmn = − w(x)Pm−1 (x) R Vmn = d2 dx2 w(x)Pn−1 (x) dx, w(x)Pm−1 (x)V (x)Pn−1 (x) dx. (3.3a) (3.3b) R The expansion coefficients cn in Eq. (3.2) can be considered as linear variational parameters. The extremum of R ψ(x)Hψ(x) dx/ R ψ(x)ψ(x) dx with respect to cn is equivalent to the diagonalization of the finite matrix Hmn [15]. Thus, we can refer to this basis set as the variational basis representation (VBR) as proposed by Light et al [25]. The eigenvalue estimates converge monotonically to the exact eigenvalues from above. Spectral convergence refers to the exponential decrease of the coefficients cn versus n as shown in Fig. 3.1 and in Fig. 6 of Ref. 28. 3.3 Pseudospectral method Pseudospectral methods were popularized by researchers interested in the numerical solution of problems in fluid dynamics [8, 32]. The set of polynomials orthogonal with respect to a chosen weight function, w(x) can be generated with the algorithm as described elsewhere [17, 18]. These polynomials define the quadrature points xk and weights wk for the quadrature rule N R w(x)f (x) dx ≈ wk f (xk ). (3.4) k=1 If the coefficients cn in Eq. (3.2) are evaluated with the quadrature rule in Eq. (3.4), the expansion of ψ(x) is written as N ψ(x) ≈ Ik (x) k=1 w(x) ψ(xk ), w(xk ) 69 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation where the interpolating polynomial, Ik (x), is given by N Ik (x) = wk Pn−1 (xk )Pn−1 (x) n=1 and Ik (xj ) = δjk satisfies the cardinality. The expansion of ψ(x) in terms of normalized interpolating functions as basis functions was discussed by Baye and coworkers [2, 4], and is referred to as a Lagrange mesh. The QDM considers the symmetric kinetic energy matrix elements that results from an integration by parts, that is, Kmn = R d dx w(x)Pm−1 (x) d dx w(x)Pn−1 (x) dx provided that the boundary term vanishes. After the derivatives are evaluated and one of the cross terms is integrated by parts, Kmn can be written as Kmn = R ′ ′ w(x)Pm−1 (x)Pn−1 (x) dx − where w(x)Pm−1 (x)V˜ (x)Pn−1 (x) dx, R 1 w′′ (x) 1 w′ (x) V˜ (x) = − 2 w(x) 4 w(x) 2 is a reference potential [39]. The Hamiltonian matrix is therefore 2 Hmn = 2µ R ′ ′ w(x)Pm−1 (x)Pn−1 (x) dx 2 + R w(x)Pm−1 (x) V (x) − 2µ V˜ (x) Pn−1 (x) dx. (3.5) As shown elsewhere [9,39], the discrete representation of the Hamiltonian is obtained with Hmn and the transformation T and we have that QDM Hij = N 2 2µ 2 k=1 Dki Dkj + V (xi ) − where Dij = √ 2µ V˜ (xi ) δij , (3.6) N ′ Pn−1 (xi )Pn−1 (xj ) wi wj n=1 70 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation is the discrete representation of the first derivative operator. Szalay [43] has provided explicit formulas for Dij = R w(x)Pi (x)Pj′ (x) dx and N k=1 Dki Dkj = ′ ′ R w(x)Pi (x)Pj (x) dx −1/2 where Pi (x) = wi Ii (x) (Table 1 of Ref. 43). For the QDM, the reference potential, V˜ (x), written in terms of the weight function forms an important aspect of the development. If w(x) can be chosen such that V (x) ≡ ( 2 /2µ)V˜ (x), the discrete representation of the QDM Hamiltonian reduces to Hij = ( 2 /2µ) N k=1 Dki Dkj , and is calculated exactly with the quadrature. Thus, the QDM preserves the variational aspects of the VBR. By contrast, applications with the DVR can lead to errors arising from the inexactness of the quadrature, Eq. (3.1), and “ghost” levels have been reported [48,49]. There is a class of Schr¨ odinger equations that are isospectral with equivalent Fokker-Planck equations [33,39]. The eigenfunction of the Fokker-Planck equation with zero eigenvalue is the equilibrium x distribution function which is known, Peq (x) = exp − W (x′ ) dx′ . For these problems, the ground-state of the Schr¨ odinger equation is ψ0 (x) = Peq (x). The function W is the superpotential in supersymmetric quantum mechanics [13, 39]. If there is no convenient choice for w(x) such that V (x) = ( 2 /2µ)V˜ (x), one can choose a w(x) so that ( 2 /2µ)V˜ (x) is close to V (x). It has been shown [3] that scaling and translating of x can often improve the convergence of the eigenvalues. The scaled and translated Hamiltonian is given by QDM Hij 2 1 = 2µ s2 N k=1 2 Dki Dkj + V (sxi + b) − 1 ˜ V (xi ) δij , 2µ s2 (3.7) where xi has been multipled by the scaling factor s and translated by b. 3.4 The vibrational states of the Morse oscillator The Morse potential for I2 has been well studied [4, 6, 28, 29]. The potential is defined by VMorse (x) = De [1 − exp(−α(x − xe ))]2 − D, (3.8) 71 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation where De = 0.0224, D = 0, α = 0.9374, xe = 0, 2 = 1, µ = 119406, and x ∈ (−∞, ∞), all in atomic units. The exact eigenvalues are given by 2 Em = 2µ α2 m + 1 2 2 α 2µ 2 De − m + 1 2 −D (3.9) for 0 ≤ m ≤ 77. The ground state eigenfunction is ψ0Morse (x) = exp − 2µ 2 De x − xe + exp(−α(x − xe )) α + α(x − xe ) . (3.10) 2 The Morse potential belongs to the family of shape invariant potentials of supersymmetric quantum mechanics [13]. If the weight function for the polynomial expansion is chosen to be w(x) = [ψ0Morse (x)]2 , the effective potential in Eq. (3.6) is simply ( 2 /2µ)V˜ (xi ) = VMorse(xi ) − E0 . With this choice of weight function, the discrete representation of the Hamiltonian, Eq. (3.6), QDM reduces to Hij = ( 2 /2µ) N k=1 Dki Dkj + E0 δij . With scaling s and translation b, the QDM representation of the Hamiltonian is given by Eq. (3.7), where V (sxi + b) = VMorse (sxi + b). This scaling and translation of the grid points was not used in Ref. 28. The QDM is almost always considered with nonclassical polynomials, whereas the DVR approach often uses classical bases (Fourier, Hermite, Laguerre, etc). It is this different choice of basis set that distinguishes the two methods. Thus a typical DVR approach to this problem is to choose a classical polynomial basis such as Hermite polynomials. The DVR representation of the discrete Hamiltonian is Herm Hij = 2 1 Herm K + VMorse (sxi + b)δij 2µ s2 ij (3.11a) with the kinetic energy operator Herm Kij = − δij − 2(N − 1) 1 x2i − + 3 2 3 +(1 − δij ) 1 2 − 2 (xi − xj )2 . (3.11b) reported by Szalay [43]. Thus it is clear that the QDM and the DVR are pseudospectral methods of solution of the Schr¨ odinger equation but are 72 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation based on different weight functions and quadrature points. A comparison of the application of both methods to this system is presented in Section 3.7. 3.5 The vibrational states of Ar2 The Cahill and Parsegian Ar2 potential [7] V (x) = a exp(−bx)(1 − cx) − x6 d + ex−6 (3.12) ˚ and V (x) in eV. The pawhere x ∈ [0, ∞) is the radial coordinate in A rameters in Eq. (3.12) are a = 1720 eV, b = 2.6920 ˚ A−1 , c = 0.2631 ˚ A−1 , d = 37.943 eV·˚ A6 , e = 177588 ˚ A12 , 2 = 0.0041801588 eV·u·˚ A2 , and µ = 20 u. This potential has eight bound states. The weight function is chosen by approximating this potential with a Morse potential in Eq. (3.8) such that the x-intercept and the minimum of the well coincides with Eq. (3.12). We thus get the parameters De = 0.01239309488 eV, D = De , A. The weight function is α = 1.685967091 ˚ A−1 , and xe = 3.761961562 ˚ Morse thus w(x) = [ψ0 (x)]2 defined on x ∈ [0, ∞). With this weight function, the QDM representation of the Hamiltonian is given by Eq. (3.7), where ( 2 /2µ)V˜ (xi ) = VMorse(xi ) + D − α De 2 /2µ + ( 2 /2µ)(α2 /4). For the DVR, the classical Laguerre polynomials defined on x ∈ [0, ∞) and orthogonal with respect to w(x) = x2 exp(−x) are used. The Hamiltonian in the Laguerre basis is given by Lag = Hij 2 1 Lag K + V (sxi + b)δij , 2µ s2 ij (3.13a) where the discrete symmetric matrix representative of the kinetic energy is given by [2] Lag Kij = 9/4x2i + Sii if i = j, −1 −1 i−j −1/2 (−1) [(3/2)(xi xj ) (xi + xj ) + Sij ] if i = j, and Sij = √ xi xj N k=1,k=i,j 1 . xk (xk − xi )(xk − xj ) (3.13b) (3.13c) 73 Chapter 3. 3.6 Pseudospectral methods of solution of the Schr¨ odinger equation Woods-Saxon potential As a third sample problem, we also consider the calculation of the bound states of the Woods-Saxon potential which is defined by V (x) = u0 u0 t , − 1 + t a0 (1 + t)2 where t = exp((x − xe )/a0 ), u0 = −50, xe = 7, a0 = 0.6, 2 = 1, µ = 0.5, and x ∈ [0, ∞). This potential has been considered by Simos [41] with an improved finite difference scheme, by Wang [46, 47] using an improved Numerov method, and by Zakrzewski [51] using a power series method. Since this potential has a shape close to the square well potential Vsq (x) = u0 0 if 0 ≤ x < L, if x ≥ L, when L = 6.2, we consider the square of the ground-state eigenfunction of Vsq (x) as the weight function, i.e. w(x) = [ψ0sq (x)]2 , where ψ0sq (x) = c1 sin(c2 x) if 0 ≤ x < L, exp(c3 x) if x ≥ L, (3.14) √ √ c1 = exp(c3 L)/ sin(c2 L), c2 = ǫ0 − u0 , and c3 = − −ǫ0 . The ground-state eigenvalue ǫ0 for Vsq (x) is the smallest root of the equation √ √ √ ǫ0 − u0 cot( ǫ0 − u0 L) = − −ǫ0 , exceeding u0 , or ǫ0 ≈ −49.75457960982555. The QDM representation is given by Eq. (3.7), where ( 2 /2µ)V˜ (xi ) = Vsq (xi ) − ǫ0 . 3.7 Results and discussions We have calculated the eigenvalues for the I2 potential, Eq. (3.8), with the diagonalization of Eq. (3.7). The convergence of E10 for three different discretizations is shown in Table 3.1. In the first column, the eigenvalue estimate obtained with the QDM with no translation or scaling of the quadrature points converges monotonically to within seven significant figures of the exact value. The underlined portions of the estimates in the table show the converged values to five significant figures. In the second column of Table 74 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation Table 3.1: Convergence of 103 E10 for the I2 Morse potential. N 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 a b c QDMa 6.397074 5.929980 5.731211 5.653056 5.629265 5.624138 5.623357 5.623268 5.623260 QDMb 6.432646 5.749709 5.634366 5.625372 5.623171 5.623269 5.623259 5.623260 Hermitec 9.203165 6.515301 5.285059 5.229673 5.177875 5.524050 5.464870 5.615021 5.619430 5.623601 5.622067 5.623292 5.623252 5.623235 5.623258 5.623257 5.623260 5.623259 5.623260 s = 1, b = 0. s = 1.05, b = 0. s = 0.14, b = 0. 75 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation 3.1, we show the improved convergence with scaling s = 1.05 although it is no longer monotonic from above. The eigenvalue estimate exhibits a minimum versus N at N = 15. The scaling s and translation b changes the weight function to [ψ0Morse ((x − b)/s)]2 so that V (x) − ( 2 /2µ)V˜ (x) is no longer a constant, and the second integral in Eq. (3.5) cannot be evaluated exactly by quadrature. Hence, the QDM representative of the Hamiltonian in Eq. (3.7) is not equivalent to the basis set representation and the variational principle no longer holds. The eigenvalues calculated with the diagonalization of the discrete Hamiltonian in the scaled Hermite basis set, Eq. (3.11a), is shown in the third column of Table 3.1. The eigenvalues converge nonmonotonically as well and slower than the QDM results. Figure 3.2 shows the convergence of E1 , E10 , E30 , and E50 for the potential given by Eq. (3.8). The negative of the relative error, defined by (N ) ǫm (N ) = log exact Em − Em , exact Em (N ) represents approximately the number of significant figures for Em . The exact are given by Eq. (3.9). The QDM results for E exact eigenvalues Em 1 without scaling shown in Fig. 3.2A converge much faster than the results with DVR with scaled Hermite polynomials. In Figs. 3.2B-3.2D, the convergence of the eigenvalues E10 , E30 , and E50 evaluated with the QDM is also much faster than the DVR results with Hermite polynomials. In order to capture the behaviour of the eigenfunctions for the higher states which are more loosely bound, the scale factor introduced for QDM serves to expand the computational domain. The scaling of the QDM is particularly important for E50 . The translation of the QDM grid, b = (1 − s)x1 , is chosen to keep the lowest grid point unchanged while the length of the computational domain is scaled by s, which is optimized by trial and error. The results with Hermite polynomials are also optimized with the values of s and b by trial and error. The fast rate of convergence for QDM is anticipated because the mth eigenfunction ψm (x) is expanded in polynomials orthogonal with respect to the weight function w(x) = [ψ0 (x)]2 . The first excited state, ψ1 (x), is the simplest one to represent in the basis set generated with this weight func(N ) tion, and thus we obtain the very rapid convergence of E1 , as shown in Fig. 3.2A. 76 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation 0 0 B log|(E10(N)-E10exact)/E10exact| log|(E1(N)-E1exact)/E1exact| A -4 QDM, s = 1 Hermite, s = 0.125, b = 0 -8 -12 -16 0 10 20 QDM, s = 1 Hermite, s = 0.14, b = 0 -8 -12 QDM, s = 1.05 -16 30 N 0 -4 0 10 20 30 40 0 -4 QDM, s = 1 -8 -12 QDM, s = 1.15 -16 0 20 40 60 80 D log|(E50(N)-E50exact)/E50exact| log|(E30(N)-E30exact)/E30exact| C Hermite, s = 0.165, b = 0.1 100 N 50 N Hermite, s = 0.2, b = 0.2 -4 -8 QDM, s = 1 -12 QDM, s = 1.3 -16 0 40 80 120 160 N (N ) exact )/E exact | versus the number Figure 3.2: Variation of log |(Em − Em m of quadrature points N for the I2 Morse potential. QDM: Eq. (3.7) with b = (1 − s)x1 ; Hermite: Eq. (3.11). (A): m = 1; (B): m = 10; (C): m = 30; (D): m = 50. 77 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation (m) In Fig. 3.1A we show the variation of the coefficients cn for the mth eigenfunction versus n for m = 1 (triangles) and m = 10 (squares). The (1) rapid exponential decrease of cn versus n is an illustration of spectral convergence. The QDM yields diagonally dominant matrix representations in the polynomial basis set as discussed previously by Lo and Shizgal [28] and thus gives the rapid spectral convergence shown in Fig. 3.2. (m) The variation of the coefficients cn Hermite expansion ψm (y) ≈ exp − 1 y−b 2 s versus n for m = 1 and 10 in the 2 N (s,b) c(m) n Qn−1 (y), (3.15) n=1 where y = sx + b is the scaled and translated coordinate, is shown in Fig. (s,b) 3.1B with triangles and squares, respectively. In Eq. (3.15), Qn (y) = νn Hn ((y − b)/s) are the scaled and translated Hermite polynomials normalized by the constants νn . Spectral convergence is also confirmed, but the (m) rate of decrease of cn is slower than with the QDM. For m = 10, the spectral convergence refers to the exponential decrease in the expansion co(10) efficients cn versus n for n 10 beyond the maxima shown in Fig. 3.1. Figure 3.3 shows the convergence of the eigenvalues E0 , E3 , E5 , and E7 for Ar2 . The potential supports 8 bound states. The exact eigenvalues, exact , used for calculating the errors are evaluated in multiple precision Em with a higher order method and are assumed to be correct to at least 20 significant digits. In Fig. 3.3, we show that the QDM discretization gives a very rapid convergence rate for the four eigenvalues when appropriate scalings are used. The results for the scaled and translated Laguerre expansion are also shown in Fig. 3.3. Since the eigenfunctions are concentrated in the region of the potential well, in order to improve the convergence of the Laguerre expansion a translation of the grid points is essential. This changes the weight function to w(y) = ((y − b)/s)2 exp(−(y − b)/s) and the computational domain to y ∈ [b, ∞). This translation redistributes the Laguerre quadrature points concentrated near x = 0 to the region of the potential well. The value of b for each state must be chosen appropriately such that the domain truncation error incurred by neglecting the contribution from [0, b) is less than machine precision and the convergence rate of the eigenvalue is nearly optimized. 78 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation 0 0 Laguerre, s = 0.02, b = 0 -4 QDM, s = 1 -8 Laguerre, s = 0.03, b = 2.5 -12 -16 0 20 40 60 80 B log|(E3(N)-E3exact)/E3exact| log|(E0(N)-E0exact)/E0exact| A N 0 QDM, s = 1.05 -8 Laguerre, s = 0.035, b = 2.5 -12 -16 100 Laguerre, s = 0.02, b = 0 -4 0 20 40 60 80 0 -4 Laguerre, s = 0.025, b = 0 -8 Laguerre, s = 0.05, b = 2.5 -12 -16 0 40 80 D log|(E7(N)-E7exact)/E7exact| log|(E5(N)-E5exact)/E5exact| C QDM, s = 1.3 120 N 100 N Laguerre, s = 0.055, b = 0 -4 -8 QDM, s = 2.0 -12 Laguerre, s = 0.1, b = 2.5 -16 0 40 80 120 160 N (N ) exact )/E exact | versus the number of Figure 3.3: Variation of log |(Em − Em m quadrature points N for the Ar2 Cahill and Parsegian potential. QDM: Eq. (3.7) with b = (1 − s)x1 ; Laguerre: Eq. (3.13). (A): m = 0; (B): m = 3; (C): m = 5; (D): m = 7. 79 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation For the Woods-Saxon potential, we use the ground state eigenfunction, Eq. (3.14), of a square well potential as the weight function which generates the basis set and quadrature points used in the QDM. The diagonalization of the matrix representative of the H given by Eq. (3.7) gives approximate eigenvalues for the Woods-Saxon potential. No translation of grid points is required (b = 0) as both potentials give the same boundary conditions for the eigenfunctions, ψ(0) = ψ(∞) = 0. To evaluate the accuracy of the approximate eigenvalues we compare these with the exact numerical values reported by Ledoux [23]. In Fig. 3.4, the values of s for the curves labeled with the open circles were chosen so that ǫm (N ) ≈ −6 and remains more or less constant with increasing N . We suspect that the minimum error of 10−6 which does not decrease with an increase in N arises from the subtraction of ( 2 /2µ)V˜ (x), which is not continuous. If larger values of s are chosen (filled circles), the accuracy of the eigenvalue estimates is improved to machine precision. 3.8 Summary The quadrature discretization method (QDM) was originally developed to solve the time dependent Boltzmann and Fokker-Planck equations for kinetic theory problems [5, 35–37, 39, 40]. This motivated the concept of developing nonclassical basis sets for different problems based on weight functions that were associated with the equilibrium solution which is the ground state eigenfunction. The spectral convergence of this large class of eigenvalue problems was found to be very rapid. The DVR was developed from the need to calculate potential energy matrix elements numerically and quadrature methods were adopted, analogous to applications of the Boltzmann equation [5,36]. The numerical implementation of the QDM and the DVR is with a pseudospectral approach that was developed by researchers in fluid mechanics [8, 32]. In this chapter, we have demonstrated that the convergence of the vibrational states of Ar2 and I2 is accelerated with the appropriate choice of the basis set modified by scaling and translation of the coordinate. The convergence with Hermite or Laguerre basis functions which are often the basis sets of choice for many applications of the DVR is generally much slower. 80 Chapter 3. Pseudospectral methods of solution of the Schr¨ odinger equation 0 0 B log|(E4(N)-E4exact)/E4exact| log|(E0(N)-E0exact)/E0exact| A -4 s=1 -8 s = 1.2 -12 -16 0 20 40 -4 -8 s = 1.3 -12 -16 60 s = 1.1 0 20 40 N 0 60 -4 s = 1.25 -8 -16 s = 1.5 0 40 80 D log|(E13(N)-E13exact)/E13exact| log|(E9(N)-E9exact)/E9exact| C -12 80 N 120 0 -4 -8 s = 1.6 s = 1.9 -12 0 N 40 80 120 160 N (N ) exact )/E exact | versus the number of Figure 3.4: Variation of log |(Em − Em m quadrature points N for the Woods-Saxon potential by QDM. (A): m = 0; (B): m = 4; (C): m = 9; (D): m = 13. 81 200 Bibliography [1] Z. A. Anastassi and T. E. Simos. A family of exponentially-fitted Runge-Kutta methods with exponential order up to three for the numerical solution of the Schr¨ odinger equation. J. Math. Chem., 41:79, 2007. [2] D. Baye and P. H. Heenen. Generalised meshes for quantum mechanical problems. J. Phys. A: Math. Gen., 19:2041, 1986. [3] D. Baye, M. Hesse, and M. Vincke. The unexplained accuracy of the Langrange-mesh method. Phys. Rev. E, 65:026701, 2002. [4] D. Baye and M. Vincke. Lagrange meshes from nonclassical orthogonal polynomials. Phys. Rev. E, 59:7195, 1999. [5] R. Blackmore and B. D. Shizgal. Discrete-ordinate method of solution of Fokker-Planck equations with nonlinear coefficients. Phys. Rev. A, 31:1855, 1985. [6] M. Braun, S. A. Sofianos, D. G. Papageorgiou, and I. E. Lagaris. An efficient Chebyshev-Lanczos method for obtaining eigensolutions of the Schr¨ odinger equation on a grid. J. Comput. Phys., 126:315, 1996. [7] K. Cahill and V. A. Parsegian. Rydberg-London potential for diatomic molecules and unbonded atom pairs. J. Chem. Phys., 121:10839, 2004. [8] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, New York, 1988. [9] H. Chen and B. D. Shizgal. The quadrature discretization method in the solution of the Schr¨ odinger equation. J. Math. Chem., 24:321, 1998. [10] H. Chen and B. D. Shizgal. A spectral solution of the Sturm-Liouville equation: comparison of classical and nonclassical basis set. J. Comput. Appl. Math., 136:17, 2001. 82 Bibliography [11] P. J. Davis and P. Rabinowitz. Methods of Numerical Integration, 2nd ed. Academic Press, New York, 1984. [12] A. S. Dickinson and P. R. Certain. Calculation of matrix elements for one-dimensional quantum mechanical problems. J. Chem. Phys., 49:4209, 1968. [13] R. Dutt, A. Khare, and U. Sukhatme. Supersymmetry, shape invariance, and exactly solvable potentials. Am. J. Phys, 56:163, 1988. [14] J. Echave and D. C. Clary. Potential optimized discrete variable representation. Chem. Phys. Lett., 190:225, 1992. [15] S. T. Epstein. The Variational Method in Quantum Chemistry. Academic Press, New York, 1974. [16] F. M. Fern´ andez, Q. Ma, and R. H. Tipping. Eigenvalues of the Schr¨ odinger equation via the Riccati-Pad´e method. Phys. Rev. A, 40:6149, 1989. [17] R. D. M. Garcia. The application of nonclassical orthogonal polynomials in particle transport theory. Prog. Nucl. Energy, 35:249, 1999. [18] W. Gautschi. Orthogonal Polynomials: Computation and Approximation. Oxford, New York, 2004. [19] D. O. Harris, G. G. Engerholm, and W. D. Gwinn. Calculation of matrix elements for one-dimensional quantum-mechanical problems and the application to anharmonic oscillators. J. Chem. Phys., 43:1515, 1965. [20] Z. Kalogiratou, T. Monovasilis, and T. E. Simos. Numerical solution of the two-dimensional time independent Schr¨ odinger equation with Numerov-type methods. J. Math. Chem., 37:271, 2005. [21] Z. Kalogiratou and T. E. Simos. Construction of trigonometrically and exponentially fitted Runge-Kutta-Nystr¨ om methods for the numerical solution of the Schr¨ odinger equation and related problems – a method of 8th algebraic order. J. Math. Chem., 31:211, 2002. [22] H. Karabulut and E. L. Sibert. Trigonometric discrete variable representations. J. Phys. B: At. Mol. Opt. Phys., 30:L513, 1997. 83 Bibliography [23] V. Ledoux, M. Rizea, L. Ixaru, G. Vanden Berghe, and M. Van Daele. Solution of the Schr¨ odinger equation by a high order perturbation method based on a linear reference potential. Comput. Phys. Commun., 175:424, 2006. [24] J. C. Light and T. Carrington. Discrete variable representations and their utilization. Adv. Chem. Phys., 114:263, 2000. [25] J. C. Light, I. P. Hamilton, and J. V. Lill. Generalized discrete variable approximation in quantum mechanics. J. Chem. Phys., 82:1400, 1985. [26] J. V. Lill, G. A. Parker, and J. C. Light. Discrete variable representations and sudden models in quantum scattering theory. Chem. Phys. Lett., 89:483, 1982. [27] M. J. Lindenfeld and B. Shizgal. Matrix elements of the Boltzmann collision operator for gas mixtures. Chem. Phys., 41:81, 1979. [28] J. Lo and B. D. Shizgal. Spectral convergence of the quadrature discretization method in the solution of the Schr¨ odinger and Fokker-Planck equations: comparison with sinc methods. J. Chem. Phys., 125:194108, 2006. [29] D. A. Mazziotti. Spectral difference methods for solving the differential equations of chemical physics. J. Chem. Phys., 117:2455, 2002. [30] T. Monovasilis, Z. Kalogiratou, T. Monovasilis, and T. E. Simos. Trigonometrically fitted and exponentially fitted symplectic methods for the numerical integration of the Schr¨ odinger equation. J. Math. Chem., 40:257, 2006. [31] T. Monovasilis and T. E. Simos. Symplectic methods for the numerical integration of the Schr¨ odinger equation. Comput. Mater. Sci., 38:526, 2007. [32] R. Peyret. Spectral Methods for Incompressible Viscous Flow. SpringerVerlag, New York, 2002. [33] H. Risken. The Fokker-Planck Equation: Methods of Solution and Application. Springer-Verlag, New York, 1984. [34] B. I. Schneider and N. Nygaard. Orthogonal functions, discrete variable representation, and generalized Gauss quadratures. 106:10773, 2002. 84 Bibliography [35] B. Shizgal. Eigenvalues of the Lorentz Fokker-Planck equation. J. Chem. Phys., 70:1948, 1979. [36] B. Shizgal. A Gaussian quadrature procedure for use in the solution of the Boltzmann equation and related problems. J. Comput. Phys., 41:309, 1981. [37] B. Shizgal and R. Blackmore. A discrete ordinate method of solution of linear boundary value and eigenvalue problems. J. Comput. Phys., 55:313, 1984. [38] B. D. Shizgal. Spectral methods based on nonclassical basis functions: the advection-diffusion equation. Computers & Fluids, 31:825, 2002. [39] B. D. Shizgal and H. Chen. The quadrature discretization method (QDM) in the solution of the Schr¨ odinger equation with nonclassical basis functions. J. Chem. Phys., 104:4137, 1996. [40] B. D. Shizgal and H. Chen. The quadrature discretization method in the solution of the Fokker-Planck equation with nonclassical basis function. J. Chem. Phys., 107:8051, 1997. [41] T. E. Simos. A new finite difference scheme with minimal phase-lag for the numerical solution of the Schr¨ odinger equation. Appl. Math. Comput., 106:245, 1999. [42] T. E. Simos. A four-step exponentially fitted method for the numerical solution of the Schr¨ odinger equation. J. Math. Chem., 40:305, 2006. [43] V. Szalay. Discrete variable representations of differential operators. J. Chem. Phys., 99:1978, 1993. [44] H. Ta¸seli and H. Alıcı. The Laguerre pseudospectral method for the reflection symmetric Hamiltonians on the real line. J. Math. Chem., 41:407, 2007. [45] H. Ta¸seli and M. Bahar Erse¸cen. The scaled Hermite-Weber basis still highly competitive. J. Math. Chem., 34:177, 2003. [46] Z. Wang. A new effective algorithm for the resonant state of a Schr¨ odinger equation. Comput. Phys. Commun., 167:1, 2005. [47] Z. Wang, Y. Ge, Y. Dai, and D. Zhao. A Mathematica program for the two-step twelfth-order method with multi-derivative for the numerical 85 Bibliography solution of a one-dimensional Schr¨ odinger equation. Comput. Phys. Commun., 160:23, 2004. [48] H. Wei. Ghost levels and near variational forms of the discrete variable representation: application to H2 O. J. Chem. Phys., 106:6885, 1997. [49] K. Willner, O. Dulieu, and F. Masnou-Seeuws. Mapped grid methods for long-range molecules and cold collisions. J. Chem. Phys., 120:548, 2004. [50] M. R. M. Witwit and J. P. Killingbeck. A Hill-determinant approach to symmetric double-well potentials in two dimensions. Can. J. Phys., 73:632, 1995. [51] A. J. Zakrzewski. Highly precise solutions of the one-dimensional Schr¨ odinger equation with arbitrary potential. Comput. Phys. Commun., 175:397, 2006. 86 Chapter 4 An efficient mapped pseudospectral method for weakly bound states; vibrational states of He2, Ne2, Ar2, and Cs2 4.1 Introduction In recent years, there has been considerable research devoted to the study of molecular systems loosely bound by weak van der Waals forces [16, 19]. The interesting aspects of these systems are their behavior at very low temperatures [8, 17], the nature of the long range interactions between atoms [2, 31, 32], and the development of computational methods to analyze such systems [10, 15, 21, 36]. The nature of helium dimers and other molecular clusters has attracted considerable attention [13,24]. These systems present a challenge for the numerical evaluations of their bound states as the wave functions for the loosely bound states are very diffuse. There are presently several different numerical methods available for the solution of the Schr¨ odinger equation. Shizgal and co-workers [7,20,26–29] developed a pseudospectral method based on nonclassical polynomials which is referred to as the quadrature discretization method (QDM). Light and coworkers [18] developed a similar collocation method referred to as the discrete variable representation (DVR) which is often employed with quadrature points based on classical polynomials. Baye and coworkers [3–5] used basis sets based on interpolation functions analogous to sinc basis functions. He refers to this method as the Lagrange mesh method. Other methods A version of this chapter has been submitted for publication. Joseph Q. W. Lo and Bernie D. Shizgal. J. Phys. B 87 Chapter 4. An efficient mapped pseudospectral method for weakly bound states such as the Fourier [9–11, 15, 25] and sinc methods [1, 33] based on uniform grids have also been used. The main motivation for the introduction of the DVR is that potential matrix elements need not be evaluated. Rather, the potential is evaluated at the quadrature points but the accuracy of the method depends on whether the collocation provides accurate estimates of the matrix elements [18]. A detailed discussion of these aspects of the DVR has been presented by Wei [34] who demonstrated that “ghost” levels that are not a part of the true spectrum can appear. The QDM in some instances such as for the Morse potential can be based on a quadrature for which the potential function does not appear explicitly and in such cases ghost levels do not occur. The efficiency of a numerical method is in terms of the convergence rate of the eigenstates with respect to the number of grid points or equivalently the number of basis functions. An efficient method is particularly important for multi-dimensional problems. Hence, the development of efficient numerical algorithms is important. The optimization of the weight function that defines the orthogonal polynomials and the quadrature points is an ongoing objective in the further development of the QDM [20]. The calculation of the vibrational states of diatoms with the weight function given by the square of the ground-state eigenfunction of the Morse potential has been shown to be very effective for the lower states because of the similarity of the Morse potential and the actual potentials, especially in the region of the well [7, 27, 29]. However, asymptotically the Morse potential decays exponentially whereas the actual potentials generally decay as an inverse power. As a result, even though the grid points can be scaled and translated, the basis functions and quadrature points derived from a Morse weight function do not provide a rapid convergence for the highest states as compared with the lower states. The reason for this is that the vibrational eigenfunctions close to the dissociation limit extend far into the classically forbidden region. The motivation of this chapter is to find a method which can produce rapidly converged eigenvalues for all bound states in the spectrum. A mapped Fourier method has been extensively used for the determination of these very loosely bound states [11, 14, 15, 36]. The basic underlying principle is to create an adaptive grid which takes into account the local de 88 Chapter 4. An efficient mapped pseudospectral method for weakly bound states Broglie wavelength for a given potential. The advantage of this approach is that the adaptive grid points are concentrated near the potential well, where the eigenfunctions undergo large variations, and the grid also extends to larger radial distances so as to capture the features of the higher-state eigenfunctions in the classically forbidden region. As a result, this method can provide accurate estimates for very loosely bound vibrational states. However, some spurious eigenvalues not in the spectrum of the Hamiltonian, referred to as “ghost” levels, can occur [14, 34, 36]. In this chapter, we propose analogous mapping procedures for polynomial methods which can improve the evaluation of the very loosely bound eigenstates and develop a technique for the elimination of ghost levels. In Section 4.2, we define the coefficient and physical space representation and introduce several maps for polynomial expansions. In Section 4.3, we apply a mapping procedure to the sinc collocation method. Section 4.4 is a summary of the maps that we use. The potentials for He2 , Ne2 , and Ar2 and the basis functions are defined in Section 4.5. A benchmark problem is reported in Section 4.6 and a possible treatment of ghost levels is suggested. In Section 4.7, we discuss the results for He2 , Ne2 , and Ar2 with the mapped polynomial and sinc methods. 4.2 4.2.1 Solution of the Schr¨ odinger equation on a grid of quadrature points The coefficient space representation The determination of the eigenvalues, Em , of the Schr¨ odinger equation 2 Hψm (x) = − 2µ ′′ ψm (x) + V (x)ψm (x) = Em ψm (x) (4.1) for x ∈ [0, ∞) often involves an expansion of the eigenfunctions in the form N (N ) ψm (x) ≈ ψm (x) = ˆ cˆ(m) n φn (x). (4.2) n=1 The choice of basis set {φˆn (x)}∞ n=1 is derived from the orthogonality of a b ∞ polynomial set {Pn (u)}n=0 given by a w(u)Pm (u)Pn (u) du = δmn with the 89 Chapter 4. An efficient mapped pseudospectral method for weakly bound states introduction of a map u = ρ(x) such that ρ−1 (b) b w(u)Pm (u)Pn (u) du = a ρ′ (x)w(ρ(x))Pm (ρ(x))Pn (ρ(x)) dx ρ−1 (a) = δmn . (4.3) In Eq. (4.3), the right hand side has been expressed in terms of the variable x, where x = ρ−1 (u) is the inverse map. Thus the basis functions are φˆn (x) = w(x)P ˆ ˆ = ρ′ (x)w(ρ(x)). n−1 (ρ(x)) with the weight function w(x) A quadrature rule for some arbitrary function g(x) with this weight function is ρ−1 (b) b w(u)g(ρ−1 (u)) du w(x)g(x) ˆ dx = ρ−1 (a) a N ≈ wk g(ρ−1 (uk )) k=1 N = wk g(xk ). (4.4) k=1 The quadrature points {uk } and weights {wk } are associated with the polynomial PN (u) [12], and the quadrature points for x are xk = ρ−1 (uk ). Equation (4.4) is exact if g(ρ−1 (u)) is a polynomial in u of degree at most 2N − 1. The quadrature rule in Eq. (4.4) can be rewritten with the absence of the weight function w(x) ˆ in the form N ρ−1 (b) g(x) dx ≈ ρ−1 (a) ηˆk g(xk ), (4.5) k=1 where the new quadrature weights are ηˆk = wk /w(x ˆ k ). ( C = ˆ mn The matrix representation of the Hamiltonian in the new basis is H C + Vˆmn , where 2 /2µ)K C ˆ mn C ˆ mn K =− C Vˆmn = ρ−1 (b) ρ−1 (a) ρ−1 (b) φˆm (x)φˆ′′n (x) dx, φˆm (x)V (x)φˆn (x) dx, (4.6a) (4.6b) ρ−1 (a) 90 Chapter 4. An efficient mapped pseudospectral method for weakly bound states ˆC and C denotes the coefficient space representation. Symmetrization of K mn is done with an integration by parts, i.e. ˆC = K mn ρ−1 (b) ρ−1 (a) φˆ′m (x)φˆ′n (x) dx (4.7) with the assumption that the boundary term vanishes. With the explicit ˆ C becomes definition of φˆn (x), K mn ρ−1 (b) ˆC = K mn ρ−1 (a) ′ ′ [ρ′ (x)]3 w(ρ(x))Pm−1 (ρ(x))Pn−1 (ρ(x)) dx ρ−1 (b) − ρ′ (x)w(ρ(x))Pm−1 (ρ(x))V˜ (x)Pn−1 (ρ(x))) dx, (4.8a) ρ−1 (a) where P ′ (ρ(x)) = (dP (u)/du)|u=ρ(x) , and 1 w′′ (ρ(x)) 1 − V˜ (x) = 2 w(ρ(x)) 4 + w′ (ρ(x)) w(ρ(x)) 1 ρ′′′ (x) 1 − 2 ρ′ (x) 4 2 [ρ′ (x)]2 2 ρ′′ (x) ρ′ (x) + w′ (ρ(x))ρ′′ (x) . (4.8b) w(ρ(x)) The notation w′ (ρ(x)) denotes (dw(u)/du)|u=ρ(x) and the second derivative is similarly defined. The integration in Eq. (4.8) is done with the quadrature rule, Eq. (4.4), that is, N N ˆC = K mn k=1 ′ ′ wk [ρ′ (xk )]2 Pm−1 (uk )Pn−1 (uk ) − wk Pm−1 (uk )V˜ (xk )Pn−1 (uk ). k=1 C is evaluated similarly. Thus, the HamiltoThe potential energy matrix Vˆmn nian is given by C ˆ mn H = 2 2µ N ′ ′ wk [ρ′ (xk )]2 Pm−1 (uk )Pn−1 (uk ) k=1 N + k=1 2 wk Pm−1 (uk ) V (xk ) − 2µ V˜ (xk ) Pn−1 (uk ). (4.9) 91 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 4.2.2 The physical space representation (m) The coefficients cˆn cˆ(m) n in Eq. (4.2) are given by ρ−1 (b) = ρ−1 (a) N φˆn (x)ψm (x) dx ≈ ηˆk φˆn (xk )ψm (xk ), k=1 where the quadrature rule Eq. (4.5) has been used. Hence, Eq. (4.2) becomes N ψm (x) ≈ (N ) ψm (x) ψm (xk )Cˆk (x), = (4.10) k=1 where the interpolating function Cˆk (x) is given by N φˆn (xk )φˆn (x) Cˆk (x) = ηˆk n=1 √ and satisfies the orthogonality ρ−1 (a) Cˆk (x)Cˆl (x) dx = ηˆk ηˆl δkl and the cardinality Cˆk (xl ) = δkl . The representation, Eq. (4.10), is referred to as the physical space representation, since the coefficients are precisely ψm (x) evaluated at the quadrature points xk . The functions Cˆk (x) are also called the Lagrange functions by Baye [4]. Equation (4.10) is an expression of the interpolation of ψm (x) with the exact values at xk . However, the eigenfunc(N ) tion ψm (x) is only an approximation to ψm (x) at values of x = xk . Also, the set Cˆk (x) depends on N . ρ−1 (b) The matrix elements of the Hamiltonian, Eq. (4.1), in the physical space ˆ P = ( 2 /2µ)K ˆ P + Vˆ P , where representation, denoted by P, is defined by H ij ij ij ˆP = − K ij VˆijP = 1 ηˆi ηˆj 1 ηˆi ηˆj ρ−1 (b) ρ−1 (a) ρ−1 (b) Cˆi (x)Cˆj′′ (x) dx, Cˆi (x)V (x)Cˆj (x) dx. (4.11a) (4.11b) ρ−1 (a) ˆ P is symmetrized and discretized in a way The kinetic energy matrix K ij similar to what is shown in Eqs. (4.7) and (4.8). The discretized Hamiltonian is given by P ˆ ij H = 2 2µ N k=1 (1) (1) [ρ′ (xk )]2 DP [ki]DP [kj] + V (xi ) − 2 2µ V˜ (xi ) δij . (4.12) 92 Chapter 4. An efficient mapped pseudospectral method for weakly bound states In Eq. (4.12), V˜ is given by Eq. (4.8b), and the representation of the first derivative operator is [20] (1) DP [ij] = √ N ′ Pn−1 (ui )Pn−1 (uj ). wi wj n=1 ˆ defined by The orthogonal transformation T Tˆkn = √ ηˆk φˆn (xk ) = wk Pn−1 (uk ) relates the coefficient space and the physical space representations of the N ˆ ˆC ˆ ˆP = N Hamiltonians in Eqs. (4.9) and (4.12). Thus, H ij m=1 n=1 Tim Hmn Tjn , and N Tˆkn cˆ(m) = n ηˆk ψm (xk ), (4.13a) k=1 N Tˆkn cˆ(m) n . ηˆk ψm (xk ) = (4.13b) n=1 4.3 The solution of the Schr¨ odinger equation with the sinc functions The sinc interpolation was originally introduced by Whittaker [35], and has been widely studied [6,22,30] and used to solve the Schr¨ odinger equation [3, 23,33,36]. Lo and Shizgal [20] have recently reported a detailed comparison of sinc methods with the QDM based on nonclassical polynomials. The sinc function is defined by sin πu if u = 0, sinc u = πu 1 if u = 0. The set of sinc interpolating functions, or the Lagrange-sinc functions [3], are given by u − uk Sk (u) = sinc (4.14) h 93 Chapter 4. An efficient mapped pseudospectral method for weakly bound states and have the property Sk (ul ) = δkl , where {uk }∞ k=−∞ form an infinite uniform grid given by uk = u0 + kh for some starting point u0 and width h. ∞ These functions satisfy the orthogonality −∞ Si (u)Sj (u) du = hδij . The interpolation of an arbitrary function g(u) in the sinc basis is given by g(u) ≈ ∞ k=−∞ g(uk )Sk (u), and is an equality if the Fourier transform of g vanishes outside [−π/h, π/h]. The quadrature rule associated with the ∞ sinc expansion is −∞ g(u) du ≈ ∞ k=−∞ hg(uk ). The popularity of the sinc interpolation is due to the simplicity of the uniformly distributed grid points. The mapping used in Section 4.2.1 can also be applied to the sinc interpolation. The mapped Lagrange-sinc functions are defined by ρ′ (x) Sk (ρ(x)) ρ′ (xk ) Cˆk (x) = (4.15) −1 ρ (∞) and the orthogonality ρ−1 (−∞) Cˆi (x)Cˆj (x) dx = ηˆi ηˆj δij holds for this new basis, where ηˆi = h/ρ′ (xi ) and ρ−1 (±∞) = limu→±∞ ρ−1 (u). The new quadrature rule becomes ∞ ρ−1 (∞) ρ−1 (−∞) ′ ρ (x)g(x) dx ≈ hg(xk ). k=−∞ For computational purposes, the sum is truncated to N terms. If the computational domain is [x1 , xN ] = [xmin , xmax ], the width of the grid uk is given by h = [ρ(xmax ) − ρ(xmin )]/(N − 1), and the grid points are uk = ρ(xmin ) + (k − 1)h and xk = ρ−1 (uk ) for k = 1, . . . , N , and the solution is in the form given by Eq. (4.10). The discretization of the Hamiltonian is similar to that in Section 4.2.2, and is given by S ˆ ij H = 2 2µ N k=1 (1) where S denotes the sinc representation and i−j (−1) (1) ′ DS [ij] = Sj (ui ) = (i − j)h 0 with 2 (1) [ρ′ (xk )]2 DS [ki]DS [kj] + V (xi ) − 1 ρ′′′ (x) 1 V˜ S (x) = − 2 ρ′ (x) 4 2µ V˜ S (xi ) δij , (4.16) if i = j, if i = j, ρ′′ (x) ρ′ (x) 2 . 94 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 4.4 Explicit transformations of grid points The map ρ introduced in Sections 4.2 and 4.3 has a range [a, b] and is an increasing function of x such that ρ′ (x) > 0 for all x. This monotonicity allows the existence of the inverse map ρ−1 (u) for all u ∈ [a, b]. The map can be as simple as scaling and translation, that is ρ(x) = (x − b)/s, or a complicated function for which ρ−1 cannot be solved explicitly, and the points xk can only be found by numerically solving the equations uk = ρ(xk ). In this chapter, we consider three different types of maps for polynomial bases, x − b2 , s2 x − b2 , ρ(x) = s1 ln s2 x − b2 ρ(x) = s1 arcsinh + b1 . s2 ρ(x) = (4.17a) (4.17b) (4.17c) Equation (4.17a) is a scaling and translation of the grid points {xk }, and we refer it as a linear map. Equation (4.17b) is an exponential map, which concentrates the grid points for small x, and separates the points for large x. Equation (4.17c) is a hyperbolic sine map, which concentrates the grid points to the region close to x = b2 . We have found it useful from the numerical experiments reported in Section 4.7 to further control the computational domain and the distribution of grid points. We do so by restricting the computational domain of x to [xmin , xmax ], where xmin and xmax are freely chosen. If Eq. (4.17b) is solved for x and written as x1 = xmin and xN = xmax , it is easy to show that the solution of the resulting two equations give s2 and b2 as (xmax − xmin ) , exp(uN /s1 ) − exp(u1 /s1 ) b2 = xmin − s2 exp(u1 /s1 ). s2 = (4.18a) (4.18b) In addition, a parameter sˆ is introduced with s1 = (uN − u1 )ˆ s. (4.18c) −1 The distribution of the grid points {xk }N ˆ, k=2 is controlled indirectly by s where a small sˆ corresponds to a high concentration of grid points near 95 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 12 A xN = xmax x = ρ-1(u) 8 4 x1 = xmin 0 0 10 u1 20 12 30 uN u 12 B C 8 8 x = ρ-1(u) xN = xmax x = ρ-1(u) xN = xmax 4 4 x1 = xmin x1 = xmin 0 0 10 u1 20 u 30 uN 0 0 10 u1 20 30 uN u Figure 4.1: The mapping of collocation points using Eq. (4.17c) from Laguerre grid uk to a modified grid xk . (A) sˆ = 0.2, b2 = 5, x ∈ [1, 10]; (B) same as A except sˆ = 0.5; (C) same as A except b2 = 2. 96 Chapter 4. An efficient mapped pseudospectral method for weakly bound states xmin , and a large sˆ corresponds to a mapping close to a linear map. Similarly, with Eq. (4.17c) the parameters s2 and b1 can be written as b1 = s1 ln 2 c exp(u1 /s1 ) − exp(uN /s1 ) c exp(−u1 /s1 ) − exp(−uN /s1 ) where c= and s2 = xmax − b2 , xmin − b2 xmax − xmin sinh uN −b1 s1 − sinh , (4.19a) (4.19b) u1 −b1 s1 . (4.19c) Equation (4.18c) is also used and sˆ controls the level of concentration of grid points near x = b2 . A low value of sˆ corresponds to a high concentration, while a high value of sˆ gives a mapping close to a linear map. Figures 4.1A and 4.1B show the different distributions {xk } versus {uk } for different values of sˆ. The parameter b2 remains as a free parameter and controls the region where points are concentrated, as illustrated in Figs. 4.1A and 4.1C. For the sinc expansion, the corresponding maps are ρ(x) = x, (4.20a) ρ(x) = ln(x − b2 ), (4.20b) ρ(x) = arcsinh x − b2 s2 , (4.20c) where the free parameters are s2 and b2 . Note that Eq. (4.20a) is referred to as an unmapped sinc expansion. The sinc discretization of the Hamiltonian is given by Eq. (4.16). The transformations in Eqs. (4.20a) and (4.20c) are used and are referred to as “sinc” and “sinc(arcsinh)”, respectively. 4.5 Weakly bound vibrational states of He2 , Ne2 , and Ar2 There have been many studies on numerical methods for vibrational states of diatoms [10, 15, 21, 36]. The lower bound states can be easily evaluated with precision. The higher states, however, are more weakly bound and 97 Chapter 4. An efficient mapped pseudospectral method for weakly bound states require much wider domains in order to describe these eigenfunctions. The main objective of this chapter is to show that with the appropriate mapping of the grid points accurate results can be obtained with a considerable reduction in the number of grid points. 4.5.1 Tang-Toennies potential The vibrational states of He2 , Ne2 , and Ar2 are used as an application for a detailed comparison of the discretization methods mentioned above. The potential model is the one given by Tang and Toennies [31]. VTT (x) = A exp(−bx) − f6 (bx) C6 C8 C10 − f8 (bx) 8 − f10 (bx) 10 , 6 x x x (4.21a) where the damping functions, fn , are given by 2n f2n (bx) = 1 − exp(−bx) k=0 (bx)k . k! (4.21b) Both VTT and x in this potential are in atomic units with the parameters in Table 4.1. For our purpose, the units of energy and distance are converted to eV and ˚ A, respectively. The potential energy curves for He2 , Ne2 , and Ar2 are illustrated in Fig. 4.2A. With application of Levinson’s theorem as discussed by Shizgal [27], zero, three, and eight negative eigenvalues can be found for He2 , Ne2 , and Ar2 , respectively, as shown in Fig. 4.3. A classical basis set; Laguerre polynomials For the semi-infinite domain, a natural choice of basis functions are the Laguerre polynomials orthogonal with weight function w(u) = u2 exp(−u). ˆ Lag , for the Schr¨ The discretized Hamiltonian, H odinger equation, Eq. (4.1), ij is given by Eq. (4.12). The eigenfunction ψm (x) can be found by ψm (xk ) = ρ′ (xk )u2k exp(−uk ) [Ψm ]k , wk where uk and wk are the Laguerre quadrature points and weights, and Ψm ˆ P , Eq. (4.12), associated with the eigenvalue Em . The is the eigenvector of H linear and exponential maps in Eqs. (4.17a) and (4.17b) are referred to as 98 Chapter 4. An efficient mapped pseudospectral method for weakly bound states Table 4.1: The parameters of the Tang-Toennies potentials for He2 , Ne2 , and Ar2 , Eq. (4.21) [31]. All parameters are in atomic units. A b C6 C8 C10 He2 41.96 2.523 1.461 14.11 183.6 Ne2 199.5 2.458 6.383 90.34 1536 Ar2 748.3 2.031 64.30 1623 49060 0 A B 10 -0.2 AS V(x) [meV] V(x) [meV] TTY 0 He2 Ne2 -0.4 TT -0.6 -0.8 -10 Ar2 2 3 4 x [Angstrom] 5 6 -1 2.6 2.8 3 x [Angstrom] 3.2 Figure 4.2: (A) The Tang-Toennies potential for He2 , Ne2 , and Ar2 ; (B) The comparison of three different He2 potentials. (solid) Tang-Toennies (TT); (dotted) Tang-Toennies-Yiu (TTY); (dashed) Aziz-Slaman (AS). 99 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 8 TT Ar2 δ0/π 6 4 2 TT Ne2 TTY/AS He2 TT He2 0 -10 -8 -6 -4 -2 log(E) Figure 4.3: The zero energy limit of the l = 0 phase shift for He2 , Ne2 , and Ar2 . 100 Chapter 4. An efficient mapped pseudospectral method for weakly bound states “Lag” and “Lag(ln)”, respectively. The first grid point x1 = xmin is chosen to be 2.7 for Ar2 and 2.0 for Ne2 for the Laguerre expansion. Therefore, for Lag, b2 = xmin − s2 u1 , and the parameters for Lag(ln) are given by Eq. (4.18). A nonclassical basis set; Morse basis functions Since the Morse potential VMorse (u) = De [1 − exp(−α(u − xe ))]2 − De (4.22) is a reasonably accurate approximation to the potential VTT , the Morse ground-state eigenfunction is used to construct a nonclassical weight function, given by w(u) = [ψ0Morse (u)]2 = exp −2 2µ 2 De u − xe + exp(−α(u − xe )) α + α(u − xe ) . (4.23) The parameters De , α, and xe are chosen by equating the root and the minimum points of VMorse with V . For Ne2 , De = 3.65111(−3) eV, xe = 3.08986 ˚ A, and α = 2.15179 ˚ A−1 . For Ar2 , De = 0.0123670 eV, xe = 3.75669 −1 ˚ ˚ A, and α = 1.79072 A . For this weight function, 1 w′′ (u) 1 − 2 w(u) 4 w′ (u) w(u) 2 = 2µ D [1 − exp(−α(u − xe ))]2 − α 2 e 2µ D + 2 e α2 . 4 ˆ QDM is given With the modified grid xk = ρ−1 (uk ), the Hamiltonian H ij 101 Chapter 4. An efficient mapped pseudospectral method for weakly bound states by Eq. (4.12), where 2µ V˜ QDM (x) = 2 De [1 − exp(−α(ρ(x) − xe ))]2 − α + 1 ρ′′′ (x) 1 − 2 ρ′ (x) 4 + −2 2µ 2 ρ′′ (x) ρ′ (x) 2µ D + 2 e α2 [ρ′ (x)]2 4 2 De [1 − exp(−α(ρ(x) − xe ))] + α ρ′′ (x). −1/2 The eigenfunctions are ψm (xk ) = ηˆk [Ψm ]k , where ηˆk is given below Eq. ˆ QDM associated (4.5), and [Ψm ]k is the kth element of the eigenvector of H with the eigenvalue Em . The maps in Eqs. (4.17a) and (4.17b) are used, and are referred to as “QDM” and “QDM(ln)”, respectively. Because the grid uk is defined by the Morse potential, u1 can be chosen to have the same value as x1 , i.e. b2 = u1 − s2 u1 for QDM. The parameters for QDM(ln) are given by Eq. (4.18) with xmin = u1 . 4.5.2 Tang-Toennies-Yiu and Aziz-Slaman He2 potentials Whereas the Tang-Toennies He2 potential defined in Eq. (4.21) does not support a bound state, the potential provided by Tang, Toennies, and Yiu [32], 12 C2n VTTY (x) = Dx7/2β−1 exp(−2βx) − f2n (bx) 2n , (4.24) x n=3 does support one bound state as verified in Fig. 4.3. In Eq. (4.24), the damping functions, f2n , are given by Eq. (4.21b), D = 7.449, β = 1.3443, C6 = 1.461, C8 = 14.11, C10 = 183.5, b = 2β − [(7/2β) − 1]/R, C2n = (C2n−2 /C2n−4 )3 C2n−6 for n = 6, . . . , 12, and VTTY , x and all parameters are in atomic units. The second He2 potential is the LM2M2 potential provided by Aziz and Slaman [2] which was used by Lombardi et al [21]. The Aziz-Slaman LM2M2 potential is given by x VAS (x) = ǫV ∗ , (4.25a) rm 102 Chapter 4. An efficient mapped pseudospectral method for weakly bound states Table 4.2: The parameters of the Aziz-Slaman LM2M2 potential for He2 , Eq. (4.25) [2]. All parameters are dimensionless, except for ǫ/k in K and rm in ˚ A. A α β c6 c8 c10 189635.353 10.70203539 -1.90740649 1.34687065 0.41308398 0.17060159 D ǫ/k rm Aa r1 r2 1.4088 10.97 2.9695 0.0026000000 1.0035359490 1.4547903690 103 Chapter 4. An efficient mapped pseudospectral method for weakly bound states where c6 c8 V ∗ (r) = A exp(−αr + βr 2 ) − 6 + 8 + r r 2 exp − D − 1 r F (r) = 1 A sin 2π r − r1 − π + 1 a r2 − r1 2 Va (r) = 0 c10 F (r) + Va (r), r 10 if r < D, (4.25b) (4.25c) if r ≥ D, if r1 ≤ r ≤ r2 , (4.25d) if r < r1 or r > r2 . with the constants as given in Table 4.2. The units of VAS and x are eV and ˚ A, respectively. Figure 4.2B shows that for He2 both VTTY and VAS are marginally wider than VTT . This small increase in width allows one negative eigenvalue for VTTY and VAS as confirmed in Fig. 4.3. The main difference between VTTY and VAS , however, is that VTTY is infinitely differentiable on (0, ∞), while the second derivative of VAS does not exist at x/rm = D, r1 , and r2 . The implementations of Laguerre and sinc bases are essentially the same as discussed before. The QDM methodology has to be modified as the Morse potential with matched values of De , α and xe does not have a bound state and hence no weight function can be constructed. We therefore adjust the value of α to 1.5 ˚ A−1 . The other parameters used to define the weight A for function, Eq. (4.23), are De = 9.4659017(−4) eV and xe = 2.9720746 ˚ ˚ VTTY , and De = 9.4532246(−4) eV and xe = 2.9695000 A for VAS . The sinh map is used for the polynomial methods, and the mapping parameters are given by Eq. (4.19). 4.5.3 Justifications for the computational domains The domains of all the potentials provided in Sections 4.5.1 and 4.5.2 are x ∈ [0, ∞), whereas the domains of the modified bases discussed above are given by x ∈ [ρ−1 (a), ∞), where a = 0 for Laguerre, and a = −∞ for QDM and sinc. The disagreement between the domain of the Schr¨ odinger equation and the domain of the basis functions is insignificant from a numerical aspect, since the eigenfunction decays very rapidly as x decreases from the classically forbidden region to zero. For the Laguerre expansion, an appro- 104 Chapter 4. An efficient mapped pseudospectral method for weakly bound states priate choice of xmin will assure that this error is below machine accuracy. For the QDM, the Morse weight function given by Eq. (4.23) is originally defined with u ∈ (−∞, ∞), but decays as exp(− exp(−(u − xe ))) as u decreases from xe . If the domain is truncated to u ∈ [0, ∞), the integration by parts done in Eq. (4.7) will yield a boundary term w(0) which is not exactly zero. However, with an appropriate choice of ρ(x), w(0) is well below machine accuracy. For the sinc expansion, the solution outside the computational domain [xmin , xmax ] is assumed to be zero. Convergence to the exact eigenvalues may be achieved only when xmin → 0 and xmax → ∞, but appropriate choices of xmin and xmax will give a domain truncation error less than machine precision. 4.6 A benchmark calculation for Cs2 ; elimination of ghost levels We first test the mapping procedure with the Morse potential in Eq. (4.22), for which exact eigenvalues are known. The Cs2 Morse potential, previously considered by Willner et al [36] is used for this benchmark. The parameters are De = 0.016627 a.u., xe = 8.77a0 , and α = 0.372031199 a.u.. The Planck’s constant is taken to be 1 a.u. and the reduced mass is µ = 121135.9042132189 a.u.. The exact eigenvalues are given by √ 2 2µDe 1 α2 2 Em = − − −m . 2µ α 2 This potential contains bound states for m = 0 to 170. The highest six bound states obtained with the QDM(arcsinh) are compared with the results obtained by Willner et al [36] in Table 4.3. The QDM(arcsinh) provides results about two to three order of magnitudes more accurate than the mapped Fourier-sine method used by Willner et al for m = 165 to m = 169 on the short interval up to xmax = 80. Both methods are based on 276 grid points. To obtain the highest state, a wider interval with xmax = 500 and N = 400 is used. The result with the QDM(arcsinh) has converged to machine accuracy, while the results reported by Willner et al have converged to only about two significant figures. 105 Chapter 4. An efficient mapped pseudospectral method for weakly bound states exact and ∆E Table 4.3: Convergence of the Morse potential for Cs2 . Em m −1 are in cm . m exact Em 165 166 167 168 169 170 -3.26122 -2.10770 -1.20493 -0.552941 -0.151714 -1.25384(-3) ∆Em WDMSa QDM(arcsinh)b QDM(arcsinh)c N = 276 1.04(-4) -1.40(-6) -1.34(-6) 1.30(-4) -6.89(-7) -1.14(-6) 1.58(-4) -4.21(-7) -9.13(-7) 1.76(-4) -2.54(-7) -6.47(-7) 1.73(-4) -1.25(-7) -3.49(-7) 1.65(-4) 2.04(-4) 1.87(-4) 170 -1.25384(-3) -1.11(-5) N = 400 -6.22(-13) -6.84(-13) a Fourier-sine results by Willner, Dulieu, and Masnou-Seeuws [36]. bs ˆ = 0.2, b2 = 10, x ∈ [6, 80] for N = 276; sˆ = 0.14, b2 = 10, x ∈ [6, 500] for N = 400. c parameters same as b , but N ′ = N + 1 for ghost level removal. 106 Chapter 4. An efficient mapped pseudospectral method for weakly bound states The use of Lag(arcsinh) does not provide eigenvalues for m ≥ 165 with 276 points, but the error of E170 is approximately −6.03(−11) with 400 points when sˆ = 0.165 and b = 10. It suggests that the mapping of polynomials, if done properly, can be more efficient than mapped Fourier methods for very loosely bound eigenvalues. The untransformed QDM does not provide E165 to E169 with 276 points, and E170 with 400 points. We have demonstrated that the mapping procedure can greatly improve the convergence of the eigenvalues. Unfortunately, it sometimes yield fictitious eigenvalues, called “ghost levels”, which do not belong to the spectrum of the Hamiltonian. The problem of ghost levels has been discussed by Wei [34], Kallush et al [14], and Willner et al [36]. For the Cs2 potential, one ghost level is found with QDM(arcsinh). The position of the ghost level versus N is shown in Fig. 4.4A. The ghost level can be easily identified by the fact that it is not converging to any value. It can also be identified by the eigenfunction coefficients, shown in Fig. 4.5A, which do not decay as n increases. The eigenfunction oscillates in the physical space as shown in Fig. 4.5B. The problem of ghost levels arises from the error of the quadrature formula when approximating the matrix elements shown in Eq. (4.6) or (4.11), as discussed by Wei [34]. No ghost levels are found when the integrals are evaluated with high accuracy. To improve the accuracy of the quadrature, a larger set of quadrature points is used in the coefficient space represenˆ C , given by Eq. (4.9), only for tation. It can be done by constructing H mn 1 ≤ m, n ≤ N ′ for some N ′ ≤ N . We use this N ′ × N ′ Hamiltonian matrix to evaluate the eigenvalues Em . The eigenfunctions in the physical space representation are N′ (N ′ ) ψm (x) = cˆ(m) n Pn−1 (ρ(x)), ρ′ (x)w(ρ(x)) n=1 (m) where cˆn is the nth component of the corresponding eigenvector. Note that the unitary transformation between the coefficient space representation and the physical space representation given by Eq. (4.13) no longer holds if N ′ < N . Because we used a higher quadrature rule to construct the Hamiltonian matrix, the discretization error is decreased, and the ghost levels which depends on this error should be sensitive to this change. Surprisingly, the 107 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 0 A B -4 -4 E -2 E -2 -6 -6 -8 -8 -10 260 270 280 290 -10 260 270 N 280 290 N Figure 4.4: The eigenvalues Em and the ghost level in logarithmic scale versus the number of collocation points N for Cs2 Morse potential. The parameters are the same as in Table 4.3. (A) N ′ = N ; (B) N ′ = N − 1. A 0 B 1 ψ(x) log |ˆ cn | -4 -8 0 -12 -1 -16 0 100 200 n 0 20 40 60 80 x Figure 4.5: The expansion coefficients and the eigenfunction of the ghost level for Cs2 with the Morse potential parameters in Table 4.3. 108 Chapter 4. An efficient mapped pseudospectral method for weakly bound states ghost level can be eliminated by merely setting N ′ = N − 1. Figure 4.4B shows the position of the ghost level when N ′ = N −1 on the smaller domain x ∈ [6, 80]. The ghost level moves into the continuum for N ≥ 266. Since only one less basis function is considered, there is no large impact on the convergence of eigenvalues, as shown in Table 4.3 under QDM(arcsinh)1 . In the larger domain x ∈ [6, 500], the ghost level moves to the continuum when N ≥ 376. With Lag(arcsinh), there is one ghost level for the calculation with x ∈ [6, 500] and N = 400. If we set N ′ = N − 1, the ghost level does not appear with N = 400. The method we used to eliminate the ghost level does not require any modification of the Schr¨ odinger equation or the basis functions. 4.7 Results and discussion As mentioned previously, the noble gas pairs He2 , Ne2 and Ar2 support 1, 3 and 8 bound states with potentials as noted in Table 4.4. The table lists the values of the eigenstates as calculated with the mapped sinc method and are accurate to the significant figures shown. We begin the discussion with the results for Ne2 with the Tang-Toennies potential, Eq. (4.21). This potential can support three bound states and the ground state, ψ0 (x), and the uppermost loosely bound state, ψ2 (x), are shown in Figs. 4.6 and 4.7 with the grid point distributions as symbols and N = 20. We also show the convergence rates for these eigenfunctions in Fig. 4.8. The logarithmic value of the relative error, (N ) ) ǫ(N m = log exact Em − Em , exact Em is approximately the number of significant figures of the eigenvalue estimates (N ) exact , are calculated with the mapped sinc Em . The exact eigenvalues, Em method in multiple precision, and are accurate to at least 20 digits with N ≈ 100. The ground state eigenfunction, ψ0 (x), in Fig. 4.6A is a single bell shaped nodeless function localized more or less in the interval x ∈ [2, 6]. The rapid numerical convergence of the eigenvalues would require a dense mesh of grid points in this region. The convergence of E0 for N 15 shown in Fig. 4.8A is most rapid with QDM followed by Laguerre and sinc with the sinh mapping and finally sinc. The mapped sinc is the most rapid for N 15. These different rates of convergence can be understood in terms of the number of grid points in the interval x ∈ [2, 6] which explains the 109 Chapter 4. An efficient mapped pseudospectral method for weakly bound states Table 4.4: The vibrational energy levels for He2 , Ne2 and Ar2 in meV. E0 E1 E2 E3 E4 E5 E6 E7 a b c TTY He2 a -1.08286(-4) AS He2 b -1.07315(-4) TT Ne2 c -2.07329 -0.362465 -1.56865(-3) TT Ar2 c -10.4703 -7.21707 -4.65637 -2.74616 -1.42543 -0.608341 -0.183899 -0.0246251 Tang-Toennies-Yiu, Eq. (4.24). Aziz-Slaman LM2M2, Eq. (4.25). Tang-Toennies, Eq. (4.21). A B sinc(arcsinh) ψ0(x) ψ0(x) Lag sinc QDM 0 2 4 6 x 8 10 0 2 4 6 8 x Figure 4.6: The ground state eigenfunction for Ne2 , ψ0 (x). For all cases, N = 20. (A) Lag: s2 = 0.06, xmin = 2; QDM: s2 = 1; (B) sinc: x ∈ [2, 8.5]; sinc(arcsinh): s2 = 0.8, b2 = 2.1, x ∈ [2, 8.5]; the collocation points for different methods are shown by the symbols. 110 10 Chapter 4. An efficient mapped pseudospectral method for weakly bound states A B Lag 5 10 15 20 25 QDM 0 5 10 x 15 20 x C sinc(arcsinh) ψ2(x) 0 QDM(ln) ψ2(x) ψ2(x) Lag(ln) 0 5 10 15 20 25 x Figure 4.7: The highest state eigenfunction of Ne2 , ψ2 (x). For all cases, N = 40. (A) Lag: s2 = 0.3, xmin = 2; Lag(ln): sˆ = 0.25, x ∈ [2, 150]; (B) QDM: s2 = 1.8; QDM(ln): sˆ = 0.29, x ∈ [u1 , 150]; (C) sinc(arcsinh): s2 = 0.8, b2 = 2.1, x ∈ [2, 150]; the collocation points for different methods are shown by the symbols. 111 25 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 A sinc B 0 QDM -4 sinc Lag QDM -8 Lag ε2(N) ε0(N) -4 Lag(ln) -8 sinc(arcsinh) -12 sinc(arcsinh) -12 -16 0 10 20 30 N 40 QDM(ln) 0 20 40 60 80 N (N ) Figure 4.8: Variation of ǫm for m = 0 and 2 versus the number of collocation points N for Ne2 . Lag: xmin = 2 for all cases, (A) s2 = 0.06, (B), s2 = 0.3; Lag(ln): (B) sˆ = 0.25, x ∈ [2, 150]; QDM: xmin = u1 for all cases, (A) s2 = 1, (B) s2 = 1.8; QDM(ln): (B) sˆ = 0.29, x ∈ [u1 , 150]; sinc: (A) x ∈ [2, 8.5], (B) x ∈ [2, 150]; sinc(arcsinh): s2 = 0.8, b2 = 2.1 for all cases, (A) x ∈ [2, 8.5], (B) x ∈ [2, 150]. 112 100 Chapter 4. An efficient mapped pseudospectral method for weakly bound states improvement of the mapped sinc results as illustrated in Fig. 4.6B. The superiority of the QDM results arises from the choice of weight function as the square of the ground state eigenfunction of the Morse potential that is a close approximation to the Tang-Toennies potential. The eigenfunction ψ2 (x), shown in Fig. 4.7, has two nodes near the origin at x ≈ 3.40 and x ≈ 4.71 and a long diffuse portion that extends far beyond x > 25. For this eigenstate, we expect that convergence will be rapid if there is a sufficient number of points densely distributed in the region x ∈ [2, 6] as for ψ0 (x) as well as a less dense distribution for much larger x. The convergence shown in Fig. 4.8B is the worst with sinc and also slow for both QDM and Laguerre basis functions without mapping. The rate of convergence is greatly improved with the appropriate mappings and QDM mapped exponentially is the most rapid. For both Laguerre and QDM shown in Figs. 4.7A and 4.7B, the exponential map, Eq. (4.17b), redistributes the grid points so that more points are in the region near the nodes where a large variation occurs, and fewer points on the decaying tail. This map improves the efficiency of the interpolation, and hence a faster convergence is expected. We did not obtain a bound state with unmapped sinc with 50 grid points, but the sinh map redistributes the grid points analogous to the mapped polynomials (Fig. 4.7C) and the convergence is greatly improved. The convergence of the highest eigenvalue, E2 , with the unmapped Laguerre and QDM is comparable as in Fig. 4.8B. No convergence is observed with the sinc method without mapping. We used the exponential map for the polynomial methods, and found that the convergence is drastically improved, where mapped QDM results yield the fastest rate of convergence. With the sinh map applied to the sinc method, the rate of the convergence is similar to the mapped Laguerre expansion. (m) Figure 4.9 shows the variation of the coefficients cˆn for ψ2 (x) versus n for the QDM (Fig. 4.9A) and the Laguerre expansion (Fig. 4.9B) with and without mapping. With the exponential mapping for both QDM and La(m) guerre, the variation of cˆn versus n exhibits a maximum and then decays rapidly, more so for QDM than Laguerre. These results are consistent with the convergence of E2 shown in Fig. 4.8B. For Ar2 , the convergence rate of the ground, 3rd, 5th, and the highest state eigenvalues are illustrated in Fig. 4.10. Without mapping, the rate of 113 Chapter 4. An efficient mapped pseudospectral method for weakly bound states A 0 B 0 (2) Lag -4 log cˆn log cˆn (2) QDM QDM(ln) -4 Lag(ln) -8 -8 0 20 40 60 80 0 20 n 40 60 n (2) Figure 4.9: Variation of the expansion coefficients cˆn versus n of the highest state eigenfunction ψ2 (x) of Ne2 . (A) QDM: s2 = 1.8, b2 = −1.5; QDM(ln): s1 = 5, s2 = 3, b2 = −2.5; (B) Lag: s2 = 0.3, b2 = 2; Lag(ln): s1 = 46, s2 = 2.75, b2 = −0.75. 114 80 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 0 A -4 B sinc -4 sinc(arcsinh) ε3(N) ε0(N) Lag sinc -8 Lag -8 QDM -12 sinc(arcsinh) -12 QDM -16 0 10 20 30 -16 40 0 20 40 N 60 N 0 C sinc -4 D 0 sinc sinc(arcsinh) Lag(ln) Lag(ln) -8 QDM QDM ε7(N) ε5(N) -4 Lag Lag -8 -12 QDM(ln) -12 QDM(ln) -16 0 20 40 60 80 N sinc(arcsinh) 0 20 40 60 80 N (N ) Figure 4.10: Variation of ǫm for m = 0, 3, 5, and 7 versus the number of collocation points N for Ar2 . Lag: xmin = 2.7 for all cases, (A) s2 = 0.03, (B) s2 = 0.04, (C) s2 = 0.05, (D) s2 = 0.1; Lag(ln): (C) sˆ = 0.7, x ∈ [2.7, 13], (D) sˆ = 0.3, x ∈ [2.7, 40]; QDM: xmin = u1 for all cases, (A) s2 = 1, (B) s2 = 1.05, (C) s2 = 1.25, (D) s2 = 1.85; QDM(ln): (B) sˆ = 1.1, x ∈ [u1 , 8.5], (C) sˆ = 0.6, x ∈ [u1 , 13], (D) sˆ = 0.3, x ∈ [u1 , 40]; sinc: (A) x ∈ [2.7, 6], (B) x ∈ [2.7, 8.5], (C) x ∈ [2.7, 13], (D) x ∈ [2.7, 40]; sinc(arcsinh): s2 = 0.7, b2 = 3.3 for all cases, (A) x ∈ [2.7, 6], (B) x ∈ [2.7, 8.5], (C) x ∈ [2.7, 13], (D) x ∈ [2.7, 40]. 115 100 Chapter 4. An efficient mapped pseudospectral method for weakly bound states convergence is more rapid with QDM except for E7 for which Laguerre and QDM results are similar. The mapped sinc method improves the convergence of all the eigenvalues and the higher ones more so than the lower ones. It is important to notice in Fig. 4.10D that a converged value of E7 cannot be obtained with unmapped sinc and N ≤ 100. The results obtained with the mapped QDM converges at the fastest rate. The Tang-Toennies-Yiu, Eq. (4.24), and the Aziz-Slaman, Eq. (4.25), potentials are used for He2 . Both potentials support one very loosely bound state. These potentials serve as examples where the convergence of the “ground-state” eigenvalue can be very difficult to obtain, and mapping provides the ability to efficiently and accurately calculate the eigenvalue. The second derivative of the Aziz-Slaman He2 potential has three discontinuities at x equal to rm D, rm r1 and rm r2 . Thus, with the mapped sinc grid we can only evaluate E0exact to six significant figures with N ≥ 400. Figures 4.11A and 4.11B show the rates of convergence of the single bound state eigenvalue for the Aziz-Slaman and Tang-Toennies-Yiu potentials, respectively. Table 4.5 compares the convergence of the single eigenvalue for QDM and Laguerre as well as the mapped QDM, Laguerre and sinc methods for the Aziz-Slaman potential. It is important to note that the convergence with QDM and Laguerre expansion is comparable and very slow. Lombardi et al [21] also obtained a very slow convergence with DVR such that they obtained eigenvalues accurate to three significant figures with N 150 and four significant figures with N 400. However, the sinh mapping for both methods greatly improves the convergence. Similarly, there is a considerable improvement for the mapped sinc method where without the mapping the bound state cannot be obtained, analogous to the highest eigenvalues for Ne2 and Ar2 . The Tang-Toennies-Yiu He2 potential is an entire function and hence does not exhibit the slow convergence found for the Aziz-Slaman potential. The exact eigenvalue E0exact can therefore be evaluated to more than 20 significant figures. The convergence rates of the eigenvalues are shown in Fig. 4.11B, where the mapping greatly improves the convergence for all methods. For the Tang-Toennies-Yiu potential, the mapped sinc method provides a faster convergence than the polynomial methods. The ground state eigenfunction for the Aziz-Slaman potential for He2 is shown in Fig. 4.12 and the x variation exhibits two different spatial scales. There is a very rapid increase near the origin followed by a slow decay. The challenge of any numerical method is to accurately capture this behavior by distributing the points appropriately. This is clearly demonstrated in Fig. 116 Chapter 4. An efficient mapped pseudospectral method for weakly bound states Table 4.5: Convergence of the energy (in µeV) for Aziz-Slaman He2 potential. N 10 20 30 40 50 100 150 200 250 300 350 400 Laga -0.4077 0.5817 0.0694 -0.0345 -0.0723 -0.1064 -0.1074 -0.1073 -0.1073 -0.1073 -0.1074 -0.1073 Lag(arcsinh)b 0.0428 -0.1044 -0.1073 QDMc 3.1934 0.2667 -0.0601 -0.0934 -0.1055 -0.1078 -0.1068 -0.1070 -0.1071 -0.1074 -0.1071 -0.1074 QDM(arcsinh)d sinc(arcsinh)e -0.3581 -28.1437 -0.1056 -0.3888 -0.1073 -0.1162 -0.1076 -0.1073 a s2 = 0.7, xmin = 1.5. sˆ = 0.11, b2 = 4.3, x ∈ [1.5, 800]. c s = 0.5, x 2 min = u1 . ds ˆ = 0.12, b2 = 4.3, x ∈ [1.5, 800] e s = 0.8, b = 3, x ∈ [1.5, 800] 2 2 b 117 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 2 A B 0 Lag 0 Lag QDM sinc(arcsinh) ε0(N) ε0(N) -2 -4 Lag(arcsinh) QDM -4 QDM(arcsinh) -8 QDM(arcsinh) -6 Lag(arcsinh) -8 0 40 80 120 sinc(arcsinh) 160 N -12 0 20 40 60 80 N (N ) Figure 4.11: Variation of ǫ0 versus the number of collocation points N for He2 . (A) Aziz-Slaman: (Table 4.5); (B) Tang-Toennies-Yiu: Lag: s2 = 1, xmin = 1.1; Lag(arcsinh): sˆ = 0.14, b2 = 3.1, x ∈ [1.1, 1800]; QDM: s2 = 0.6, xmin = u1 ; QDM(arcsinh): sˆ = 0.14, b2 = 3.5, x ∈ [1.1, 1800]; sinc(arcsinh): s2 = 0.6, b2 = 1.8, x ∈ [1.1, 1800]. 118 100 Chapter 4. An efficient mapped pseudospectral method for weakly bound states A B QDM(arcsinh) ψ0(x) ψ0(x) Lag(arcsinh) Lag 200 400 600 800 0 200 x 400 600 x C sinc(arcsinh) ψ0(x) 0 QDM 0 200 400 600 800 x Figure 4.12: The eigenfunction of the Aziz-Slaman He2 potential. For all cases, N = 50. (A) Lag: s2 = 1, xmin = 1.5; Lag(arcsinh): sˆ = 0.11, b2 = 4.3, x ∈ [1.5, 800]; (B) QDM: s2 = 0.5, xmin = u1 ; QDM(arcsinh): sˆ = 0.12, b2 = 4.3, x ∈ [1.5, 800]; (C) sinc(arcsinh): s2 = 0.8, b2 = 3, x ∈ [1.5, 800]; the collocation points for different methods are shown by the symbols. 119 800 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 4.12 for all three methods. With the mappings, there is an increase in the density of points both near the origin as well as for large x. The grid for the mapped sinc method provides more points in this asymptotic region than the mapped polynomial methods, but the convergence of the eigenvalues is somewhat slower than with the polynomial methods as shown in Table 4.5. As is well-known, the numerical calculation of the Hamiltonian matrix elements, Eq. (4.6), are inexact and as a consequence ghost levels can appear. No ghost levels were found for Ne2 and Ar2 with all methods as the results converge to the exact eigenvalues. However, there is one ghost level with both the mapped QDM and the mapped Laguerre method for both He2 potentials. The position of the ghost levels for Aziz-Slaman potential with mapped QDM and Laguerre are shown in Fig. 4.13A and 4.13C. If the technique for the elimination of the ghost levels discussed in Section 4.6 is used, the ghost level moves up from approximately −10−4 to −10−7 with N ′ = N − 1. The position is shown in Fig. 4.13B and 4.13D. The ghost level can be completely eliminated if N ′ = N − 2. For Tang-Toennies-Yiu potential, the ghost level can be eliminated by setting N ′ = N − 1. 120 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 8 8 A B E0 6 ghost level 4 2 E0 -log(-E) -log(-E) 6 ghost level 4 0 20 40 60 80 2 100 0 20 40 N 60 80 C D ghost level 8 -log(-E) -log(-E) 8 E0 6 100 N E0 6 ghost level 4 4 0 20 40 60 N 80 100 0 20 40 60 80 N Figure 4.13: The eigenvalue E0 and the ghost level in logarithmic scale versus the number of collocation points N for Aziz-Slaman He2 potential. The parameters are the same as in Table 4.5. Lag(arcsinh): (A) N ′ = N , (B) N ′ = N − 1; QDM(arcsinh): (C) N ′ = N , (D) N ′ = N − 1. 121 100 Bibliography [1] P. Amore. A variational sinc collocation method for strong coupling problems. J. Phys. A: Math. Gen., 39:L349, 2006. [2] R. A. Aziz and M. J. Slaman. An examination of ab initio results for the helium potential energy curve. J. Chem. Phys., 94:8047, 1991. [3] D. Baye. Lagrange-mesh method for quantum-mechanical problems. Phys. Stat. Sol. (b), 243:1095, 2006. [4] D. Baye and P. H. Heenen. Generalised meshes for quantum mechanical problems. J. Phys. A: Math. Gen., 19:2041, 1986. [5] D. Baye, M. Hesse, and M. Vincke. The unexplained accuracy of the Langrange-mesh method. Phys. Rev. E, 65:026701, 2002. [6] J. P. Boyd. Chebyshev and Fourier Spectral Methods. Dover, New York, 2000. [7] H. Chen and B. D. Shizgal. A spectral solution of the Sturm-Liouville equation: comparison of classical and nonclassical basis set. J. Comput. Appl. Math., 136:17, 2001. [8] C. Chin and P. S. Julienne. Radio-frequency transitions on weakly bound ultracold molecules. Phys. Rev. A, 71:012713, 2005. [9] D. T. Colbert and W. H. Miller. A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method. J. Chem. Phys., 96:1982, 1992. [10] O. Dulieu, R. Kosloff, F. Masnou-Seeuws, and G. Pichler. Quasibound states in long-range alkali dimers: grid method calculations. J. Chem. Phys., 107:10633, 1997. [11] E. Fattal, R. Baer, and R. Kosloff. Phase space approach for optimizing grid representations: the mapped Fourier method. Phys. Rev. E, 53:1217, 1996. 122 Bibliography [12] W. Gautschi. Orthogonal Polynomials: Computation and Approximation. Oxford, New York, 2004. [13] R. E. Grisenti, W. Sch¨ ollkopf, J. P. Toennies, G. C. Hegerfeldt, T. K¨ ohler, and M. Stoll. Determination of the bond length and binding energy of the helium dimer by diffraction from a transmission grating. Phys. Rev. Lett., 85:2284, 2000. [14] S. Kallush and R. Kosloff. Improved methods for mapped grids: applied to highly excited vibrational states of diatomic molecules. Chem. Phys. Lett., 433:221, 2006. [15] V. Kokoouline, O. Dulieu, R. Kosloff, and F. Masnou-Seeuws. Mapped Fourier methods for long-range molecules: application to perturbations in the Rb2 (0+ u ) photoassociation spectrum. J. Chem. Phys., 110:9865, 1999. [16] D. M. Leitner, J. D. Doll, and R. M. Whitnell. Quantum mechanics of small Ne, Ar, Kr, and Xe clusters. J. Chem. Phys., 94:6644, 1991. [17] J. L´eonard, M. Walhout, A. P. Mosk, T. M¨ uller, M. Leduc, and C. Cohen-Tannoudji. Giant helium dimers produced by photoassociation of ultracold metastable atoms. Phys. Rev. Lett., 91:073203, 2003. [18] J. C. Light and T. Carrington. Discrete variable representations and their utilization. Adv. Chem. Phys., 114:263, 2000. [19] Y. D. Liu and P.-N. Roy. Energy levels and wave functions of weaklybound 4 Hex 20 Ney H (x + y = 2) system using Pekeris coordinates and a symmetry-adapted Lanczos approach. J. Chem. Phys., 121:6282, 2004. [20] J. Lo and B. D. Shizgal. Spectral convergence of the quadrature discretization method in the solution of the Schr¨ odinger and Fokker-Planck equations: comparison with sinc methods. J. Chem. Phys., 125:194108, 2006. [21] M. Lombardi, P. Barletta, and A. Kievsky. Variational estimates using a discrete variable representation. Phys. Rev. A, 70:032503, 2004. [22] J. Lund and K. L. Bowers. Sinc Methods for Quadrature and Differential Equations. SIAM, Philadelphia, 1992. [23] D. A. Mazziotti. Spectral difference methods for solving the differential equations of chemical physics. J. Chem. Phys., 117:2455, 2002. 123 Bibliography [24] W. Sch¨ ollkopf and J. P. Toennies. The nondestructive detection of the helium dimer and trimer. J. Chem. Phys., 104:1155, 1996. [25] C. Schwartz. High-accuracy approximation techniques for analytic functions. J. Math. Phys., 26:411, 1985. [26] B. Shizgal. A Gaussian quadrature procedure for use in the solution of the Boltzmann equation and related problems. J. Comput. Phys., 41:309, 1981. [27] B. D. Shizgal. The quadrature discretization method (QDM) in the calculation of the rotational-vibrational transitions in rare gas dimers. J. Mol. Struc. (Theochem), 391:131, 1997. [28] B. D. Shizgal and H. Chen. The quadrature discretization method (QDM) in the solution of the Schr¨ odinger equation with nonclassical basis functions. J. Chem. Phys., 104:4137, 1996. [29] B. D. Shizgal and H. Chen. The quadrature discretization method in the solution of the Fokker-Planck equation with nonclassical basis function. J. Chem. Phys., 107:8051, 1997. [30] F. Stenger. Numerical Methods Based on Sinc and Analytic Functions. Springer-Verlag, New York, 1993. [31] K. T. Tang and J. P. Toennies. The van der Waals potentials between all the rare gas atoms from He to Rn. J. Chem. Phys., 118:4976, 2003. [32] K. T. Tang, J. P. Toennies, and C. L. Yiu. Accurate analytical He-He van der Waals potential based on perturbation theory. Phys. Rev. Lett., 74:1546, 1995. [33] G. W. Wei. Solving quantum eigenvalue problems by discrete singular convolution. J. Phys. B: At. Mol. Opt. Phys., 33:343, 2000. [34] H. Wei. Ghost levels and near variational forms of the discrete variable representation: application to H2 O. J. Chem. Phys., 106:6885, 1997. [35] E. T. Whittaker. On the functions which are represented by the expansions on the interpolation theory. Proc. Roy. Soc. Edinburgh, 35:181, 1915. [36] K. Willner, O. Dulieu, and F. Masnou-Seeuws. Mapped grid methods for long-range molecules and cold collisions. J. Chem. Phys., 120:548, 2004. 124 Chapter 5 Conclusions and future directions 5.1 Preliminary investigations in mappings and other directions In this section, I provide a continuation of the work of Chapter 4, which focuses on the mapping of grid points. In Chapter 4, explicit mappings such as the exponential or the hyperbolic sine mappings, Eq. (4.17), were introduced. There were, however, no discussions about the optimization of the parameters. In this section, I will suggest a method for determining appropriate parameters for these explicit mappings. It is expected that each type of explicit mappings is appropriate for only a small class of functions. For example, in the case of interpolation using hyperbolic sine mapping only functions with large variations concentrated in one narrow region can be beneficial from the redistribution of grid points. To add flexibility to the mapping, I hope to develop a technique which can construct an appropriate mapping based on the behaviour of the function. In this section, I will focus only on the sinc interpolation [8, 13]. The motivation of the investigations in mapping methods is to find a simple way to reduce the number of grid points required to solve timedependent or eigenvalue problems when classic basis functions can be used. For 2D problems, the size of the matrices get larger so that the methods developed in this thesis for 1D problems will have important applications to 2D problems. In section 5.1.3, I will provide some preliminary results for the eigenvalues of the Kramers equation to illustrate that mapping will be useful in 2D problems. I will at first start with interpolation of functions, and apply the same technique for interpolation to time-dependent problems and eigenvalue prob125 Chapter 5. Conclusions and future directions lems. 5.1.1 Interpolation with mapping Consider the interpolation of a function g(x) in the domain x ∈ [xmin , xmax ]. The sinc interpolating functions, Sj (u), defined on the uniform grid uj = xmin + (j + 1)h and h = (xmax − xmin )/(N − 1) are defined in Eq. (4.14). With an introduction of a map u = ρ(x), the new mapped sinc interpolating functions, Cˆj (x), are given by Eq. (4.15), where the mapped grid points are xj = ρ−1 (uj ). The new grid xj is referred to as the physical grid, and the original grid uj is referred to as the computational grid. Without loss of generality, it is convenient to require the map to satisfy ρ(xmin ) = xmin and ρ(xmax ) = xmax . The map must also satisfy ρ′ (x) > 0 for x ∈ [xmin , xmax ]. ˆ The expansion f (x) = N j=1 g(xj )Cj (x) is referred to as the interpolant of g(x). The Gaussian function 1 1 G(x; µ, σ) = √ exp − 2 σ 2π x−µ σ 2 . (5.1) is used to illustrate of the use of mappings. I consider the interpolation of the function g(x) = 0.75G(x; 0, 0.5) + 0.25G(x; 0.5, 0.02) (5.2) in the interval x ∈ [−2, 2]. This function is very similar to the Gaussian bell curve except there is a narrow spike at x = 0.5. This function is interpolated by 160 uniformly distributed grid points with no mapping, i.e. u = ρ(x) = x. Figure 5.1A shows the function g(x) with the grid points. The error g(x) − f (x) is shown in Fig. 5.1B, with a maximum magnitude of approximately 0.1. The function g(x) can be interpolated with much fewer grid points if the grid points are dense near the region x = 0.5 and sparse elsewhere. One way to redistribute the grid is to uniformly place the grid points along the curve so that the arc lengths of the curve between any two consecutive grid points are the same. This is an example of grid redistribution based on the equidistribution principle. 126 Chapter 5. Conclusions and future directions 6 A B 0.1 4 g(x) g(x) - f(x) 0.05 2 0 -0.05 0 -0.1 -2 -1 0 x 1 2 -2 -1 0 1 x Figure 5.1: Unmapped sinc interpolation of Eq. (5.2) using 160 grid points. (A) Grid point distribution; (B) Error of the interpolation. 127 2 Chapter 5. Conclusions and future directions The equidistribution principle is a method for placing the grid points, xj+1 {xj }N M (x) dx is equally distributed [1]. The j=1 , so that the quantity xj function M (x) is referred to as the monitor function. The monitor function can be any function which reflects some property of g(x). For arc length equidistribution, M (x) = 1 + [g′ (x)]2 . Explicitly defined mapping Since a large variation of g(x), as shown in Fig. 5.1A, occurs at a narrow region near x = 0.5, a high concentration of grid points is required in order to describe the behaviour of g(x) correctly. A hyperbolic sine mapping [8], Eqs. (4.17c), (4.18c), and (4.19), with u1 = xmin , uN = xmax and b2 ≈ 0.5 will redistribute a uniform grid towards x = 0.5. The optimization of the mapping parameters sˆ and b2 can be done by the equidistribution principle. I propose a new monitor function M (x; α) = 1 + α|g(x)| + α|g′ (x)| (5.3) which includes the contribution from the magnitude of both the value of g(x) and its derivative. The parameter α controls the level of grid redistribution between the regions where g(x) ≈ 0 and g′ (x) ≈ 0, and the regions where the magnitude or the derivative of g(x) is large. Based on the equidistribution principle, the grid should be set up in a x way that the quantity xjj+1 M (x; α) dx is distributed as evenly as possible. I therefore take the values of sˆ and b2 which minimize the function N −1 2 xj+1 E(ˆ s, b2 ) = M (x; α) dx j=1 xj where the physical grid is defined by xj = ρ−1 (uj ). The integrals and the minimization are evaluated by Riemann sum and steepest descent, respectively. For the interpolation of g(x), the above procedure yields the values sˆ = 0.115 and b2 = 0.498 when α = 100 and N = 18. Figures 5.2A and 5.2B show the interpolation of g(x) and its error. The maximum error with 18 points under the sinh transformation is comparable to the error with 160 equally spaced grid points shown in Fig. 5.1B. 128 Chapter 5. Conclusions and future directions 6 A 4 0.05 g(x) - f(x) g(x) B 0.1 2 0 -0.05 0 -0.1 -2 -1 0 x 1 2 -2 -1 0 1 x Figure 5.2: Sinc interpolation of Eq. (5.2) with hyperbolic sine mapping using 18 grid points. (A) Grid point distribution; (B) Error of the interpolation. 129 2 Chapter 5. Conclusions and future directions Implicitly defined mapping An appropriately distributed grid can be constructed directly by the equidistribution principle without a predefined mapping function such as the hyperbolic sine function used above. This idea has been widely used in the finite difference and finite element [4, 9], and I attempt to apply it to pseudospectral methods. Since the mapping is constructed implicitly, it is referred to as an implicit mapping [1]. First, a reference grid {ˆ xj }N j=1 is constructed by satisfying x ˆj+1 xmax 1 M (x; α) dx = M (x; α) dx (5.4) N − 1 xmin x ˆj for all k with x ˆ1 = xmin and x ˆN = xmax . To simplify the evaluation of the integrals, M (x; α) can be considered as a piecewise constant function approximating Eq. (5.3). The reference grid can be used directly for a finite difference method. For the pseudospectral method, however, a differentiable mapping function ˆ (x) as a piecewise linear function joining the ρ(x) is required. By defining U ˆ (x) by least squares. I points (ˆ xj , uj ), the map ρ(x) is defined by fitting U propose ρ(x) to be of the form Nfit ρ(x) = x + cn sin nπ n=1 x − xmin xmax − xmin (5.5) so that ρ(xmin ) = xmin and ρ(xmax ) = xmax . The sum on the right hand side of Eq. (5.5) is analogous to a Fourier-sine expansion of order Nfit , and one can take advantage of the orthogonality of sine basis when solving for the coefficients cn . The coefficients cn are found by minimizing the error function xmax E= xmin ˆ (x)]2 + β[ρ′′ (x)]2 [ρ(x) − U dx, where ρ(x) is given by Eq. (5.5) and β is a smoothing parameter. Smoothing may be necessary to prevent any oscillation in ρ(x) produced after fitting. Figure 5.3A shows that for β = 0, the fitting generates a lot of small oscillations. The oscillations disappear, in Fig. 5.3B, when a small smoothing β = 3 × 10−5 is used. In both cases, N = 30, Nfit = 24, and α = 100. 130 Chapter 5. Conclusions and future directions A 2 1 x = ρ-1(u) 1 x = ρ-1(u) B 2 0 0 -1 -1 -2 -2 -2 -1 0 u 1 2 -2 -1 0 1 u Figure 5.3: The reference grid x ˆj (crosses) and the map ρ(x). The parameters are N = 30, Nfit = 24, and α = 100. (A) β = 0; (B) β = 3 × 10−5 . 131 2 Chapter 5. Conclusions and future directions Because of the orthogonality of sine basis, the system of equations ∂E/∂cn = 0 is uncoupled. The coefficients cn therefore have an explicit form cn = 2 xmax − xmin + β π 4 n4 (xmax − xmin )3 xmax × xmin −1 ˆ (x) − x] sin nπ [U x − xmin xmax − xmin dx. ˆ (x) − x is a piecewise linear function, the integral on the right side Since U can be evaluated analytically. The new physical grid, xj , is calculated numerically by solving the equations ρ(xj ) = uj . The interpolation and its error are plotted in Fig. 5.4. The error with 30 points is similar to the errors for unmapped and explicitly mapped interpolation discussed above. 5.1.2 Interpolatory moving mesh The method of adaptive grid and implicit mapping discussed in Section 5.1.1 provides an efficient method for interpolation as the distribution of the grid is optimized based on the behaviour of the function being interpolated. For a time-dependent function, it is possible to construct a grid for each time step so the grid moves as the function evolves along time. The timedependent grid is referred to as a moving mesh. The moving mesh method for initial value problems based on finite difference or finite element methods has been widely studied in the context of fluid dynamics [4, 9]. In contrast, similar adaptive methods for pseudospectral methods are seldom mentioned. In this section, the implicit mapping of the grid discussed in Section 5.1.1 is applied to the solution of the time-dependent Smoluchowski equation in Eq. (1.7). The eigenvalues of the same Smoluchowski operator were studied in Chapter 2. For this time-dependent problem, a moving mesh method is used to illustrate the effectiveness of the interpolation of the solution when the initial condition and the steady-state solution have very different behaviours. The simplest type of moving mesh involves the construction of a new grid based on the behaviour of the solution at each time step [14], and is therefore referred to as the interpolatory moving mesh method. Since reinterpolation of the solution and reconstruction of the Smoluchowski matrix operator is necessary with any change of the grid, to reduce computational 132 Chapter 5. Conclusions and future directions 6 A 4 0.05 g(x) - f(x) g(x) B 0.1 2 0 -0.05 0 -0.1 -2 -1 0 x 1 2 -2 -1 0 1 x Figure 5.4: Sinc interpolation of Eq. (5.2) with implicitly defined mapping using 30 grid points. The parameters are Nfit = 24, α = 100, and β = 3 × 10−5 . (A) Grid point distribution; (B) Error of the interpolation. 133 2 Chapter 5. Conclusions and future directions costs, a new grid is used only when it differs from the old grid by a significant amount. The self-adjoint form of the Smoluchowski equation given by Eq. (1.10) is ∂p ǫ ∂ ∂ = √ P0 ∂t ∂x ν P0 ∂x p √ P0 (5.6) P0 (x)p(x, t). The solution p(x, t) is represented as an when P (x, t) = interpolation N p(x, t) = aj (t) j=1 where aj (t) = = LS p Cj (x) , ηj (5.7) √ ηj p(xj , t). The initial condition is interpolated with the implicit mapping generated from the monitor function in Eq. (5.3). The maps, the physical grids, and the interpolating functions used at subsequent time steps tn are denoted by [n] [n] ρ[n] (x), xj , and Cj (x), respectively. All other quantities required to be redefined are denoted with [n]. The details for the discretization and time integration of Eq. (5.6) are given in Appendix D. [n] The grid points xj at tn are used to determine the solution p(x, t), Eq. (D.2), at tn+1 . Before proceeding to the next step tn+2 , one needs to decide whether a new grid is necessary. The decision is made based on the difference between the reference grid used to generate the current map ρ[n] (x) and the reference grid determined by the behaviour of the solution at current time step tn+1 . Since the solution p at tn+1 given by Eq. (D.2) is defined on the grid the monitor function can be considered as a piecewise function [n] xj , [n] M [n+1] (x; α) = 1 + α ˆ [n] p(xj+1 , tn+1 ) + p(xj , tn+1 ) 2 [n] +α ˆ [n] [n] p(xj+1 , tn+1 ) − p(xj , tn+1 ) [n] 1 2 [n] xj+1 − xj [n] for xj < x < xj+1 and j = 1, . . . , N − 1. Comparing with Eq. (5.3), the second term is replaced by the midpoint values, and the third term is re134 Chapter 5. Conclusions and future directions placed by the difference formula. Because the magnitude of p(x, t) changes [n] over time, a new parameter α ˆ = α/[maxj p(xj , tn+1 )] is introduced to keep the monitor function consistent at each time step. [n+1] The reference grid x ˆj is generated by the equidistribution principle [n+1] x ˆj+1 [n+1] x ˆj M [n+1] (x; α) dx = xmax 1 N −1 M [n+1] (x; α) dx (5.8) xmin [n+1] as in Eq. (5.4). If the difference between the reference grids x ˆj deter- [n] x ˆj mined by Eq. (5.8) and the reference grid used at the previous time step is small, i.e. [n+1] [n] max x ˆj −x ˆj ≤ ∆ˆ x j for some tolerance level ∆ˆ x, no new grid is required. One will therefore set [n+1] [n] [n+1] [n] xj = xj , ρ[n+1] (x) = ρ[n] (x), aj (tn+1 ) = aj (tn+1 ), and copy the [n+1] values of all other quantities with [n] to [n + 1]. The values of x ˆj by Eq. (5.8) will be abandoned, and will be set to [n+1] On the other hand, if maxj x ˆj [n+1] x ˆj [n] −x ˆj [n+1] x ˆj = found [n] x ˆj . exceeds the tolerance level, is used to construct ρ[n+1] (x) as described in Section 5.1.1, and a new [n+1] physical grid xj [n+1] is generated. All other quantities associated with xj [n+1] are required to be updated. The solution p(xj [n+1] p(xj , tn+1 ) N = i=1 [n] ai (tn+1 ) [n] , tn+1 ) is interpolated by [n] [n+1] Ci (xj ) ηˆi [n+1] and the updated expansion coefficient set to aj (tn+1 ) = [n+1] ηˆj [n+1] p(xj , tn+1 ). The following example illustrates the use of moving mesh for solving Eq. (5.6) in the case ǫ = 0.1 and ν = 1. The initial condition P (x, 0) = G(x; 0.5, 0.02) is considered, where G is the Gaussian function defined in Eq. (5.1). The parameters used for the moving mesh are xmin = −2.2, xmax = 2.2, α = 20, β = 10−4 , Nfit = 24, and ∆ˆ x = 0.01(xmax − xmin ). 135 Chapter 5. Conclusions and future directions 4 4 B 3 3 2 2 P(x,0.1) P(x,0) A 1 0 -1 1 0 -2 -1 0 1 -1 2 -2 -1 x 0 1 2 x 4 4 C D 3 P(x,25.6) P(x,1.6) 3 2 1 0 -1 2 1 0 -2 -1 0 x 1 2 -1 -2 -1 0 1 x Figure 5.5: Solution of the Smoluchowski equation by interpolatory moving mesh method. The distribution of grid points is plotted under each graph. (A) t = 0; (B) t = 0.1, (C) t = 1.6, (D) t = 25.6. 136 2 Chapter 5. Conclusions and future directions 1000 100 t 10 1 0.1 -2 -1 0 1 2 x Figure 5.6: The evolution of grid points over time for the solution of Smoluchowski equation. ×10-4 ×10-4 A moving mesh 4 0 -4 unmapped -8 B 8 Pref (x,25.6) - P(x,25.6) Pref (x,1.6) - P(x,1.6) 8 moving mesh 4 0 -4 unmapped -8 -2 -1 0 x 1 2 -2 -1 0 1 x Figure 5.7: The error of the solution of Smoluchowski equation for interpolatory moving mesh method (N = 40) and unmapping sinc method (N = 141). (A) t = 1.6; (B) t = 25.6. 137 2 Chapter 5. Conclusions and future directions The solution with N = 40 and ∆tn = 0.01 is shown in Fig. 5.5 with the grid distribution shown underneath. The dashed curve in Fig. 5.5D represents the steady state solution P0 (x). Figure 5.6 shows the distribution of the moving grid points at any time step. The sawtooth behaviour of the curves suggests that the grid is not updated at every time step. In fact, there are only 64 grid updates among 10000 time steps from t = 0 to t = 100 with the chosen tolerance level ∆ˆ x. The accuracy of the moving mesh solution with 40 grid points is compared with the solution calculated by 141 unmapped grid points in Fig. 5.7 at t = 1.6 and 25.6. Since the exact solution is unknown, a reference solution, denoted by Pref (x, t), solved by using 701 unmapped points is used for measuring the accuracy of both solutions. The value of the moving mesh solution at the reference grid points are found by interpolation. On the other hand, the grid points of the unmapped solution coincide with the grid points of Pref (x, t), and no interpolation is necessary. The larger error of the moving mesh solution for x < 0 shown in Fig. 5.7A is due to the much farther separation of grid points in the region than those grid points in the unmapped solution. The overall accuracy of the moving mesh solution with 40 points, however, is better than the unmapped solution with 140 points because the initial condition P (x, 0) is more accurately interpolated with mapping. 5.1.3 Eigenvalues of the Smoluchowski and Kramers equation Smoluchowski equation It has been shown in Chapter 4 that mapping is an efficient method for calculating the eigenvalues of highly diffuse states of the Schr¨ odinger equation. In this section, I will demonstrate the use of implicit mapping for the evaluation of eigenvalues of the Smoluchowski, Eq. (1.7), and Kramers equations, Eq. (1.6). Transformation of grid points is effective for small ǫ because the structure of these eigenfunctions is localized. The value used in this section is ǫ = 0.001. The eigenvalues λm of the Smoluchowski operator LS in Eq. (5.6), LS ψm (x) = λm ψm (x), are evaluated by finding the eigenvalues of the discretized matrix LS given by Eq. (D.1) without [n]. For implicit mapping, a preliminary 138 Chapter 5. Conclusions and future directions estimation of the grid distribution is required as explained in Section 5.1.1. This estimation requires the knowledge of the regions which the structure of the eigenfunctions ψm (x) are localized. The information for constructing the grid can be provided by very rough sketches of the eigenfunctions. To get the rough sketches, the eigenvalue problem is solved in the first place by a low order unmapped sinc method. With u = ρ(x) = x, the uniform grid uk is used for constructing the Smoluchowski matrix operator. In this example, N = 20 uniform grid points in ˜S , the interval [−1.2, 1.2] are used. Among the eigenfunctions found by L ij two of them used for the grid estimation are shown in Figs. 5.8A and 5.8B. These sketches of eigenfunctions are denoted by ψ˜A (x) and ψ˜B (x), and are scaled to a maximum value of 1. The reference grid x ˆj is found by Eq. (5.4) with the piecewise constant monitor function ˜ j+1 ) + ψ(u ˜ j) ˜ j+1 ) − ψ(u ˜ j) ψ(u ψ(u M (x; α) = 1 + α +α 2 uj+1 − uj 1 2 for uj < x < uj+1 and j = 1, . . . , N − 1. Compared with Eq. (5.3), the second term is replaced by the midpoint value, and the third term is replaced with the difference formula. The parameters used to construct the maps ρ(x) are α = 100, β = 10−5 and Nfit = 12. Figures 5.9A and 5.9B shows the maps ρA (x) and ρB (x) generated by ψ˜A and ψ˜B , respectively. A third map ρC (x) generated by ψ˜C (x) = ψ˜A (x) + ψ˜B (x) and scaled to a maximum value of 1 is also shown in Figs. 5.8C and 5.9C. Tables 5.1 and 5.2 show the eigenvalues with six significant figures found by the three mappings. The unmapped results are also included for comparison. For ρA (x), most of the grid points are concentrated in the vicinity of x = 0. Since the eigenfunctions ψ3 , ψ4 , ψ7 , and ψ8 have their main structure near x = ±1, the corresponding eigenvalues are not found with ρA (x). On the other hand, the eigenfunctions ψ2 , ψ5 , ψ6 , ψ9 , and ψ10 have structure localized near x = 0, the same map can calculate the corresponding eigenvalues very efficiently. The same argument applies to ρB (x) which the grid points are localized in the vicinity of x = ±1. This map can produce the eigenvalues λ3 , λ4 , λ7 , and λ8 . It is worth noting that using this map will require approximately two times the number of grid points than the map ρA (x) in order to calculate the eigenvalues to six significant figures. The reason is that the grid points 139 Chapter 5. Conclusions and future directions B 0.8 0.8 p˜B (x) 1.2 0.4 0.4 0 0 -0.4 -0.4 -1 0 1 -1 0 x 1 x C 1.2 0.8 p˜C (x) p˜A (x) A 1.2 0.4 0 -0.4 -1 0 1 x Figure 5.8: Eigenfunctions of Smoluchowski equation solved by unmapped sinc method with 20 grid points. (A) and (B) two eigenfunctions; (C) sum of the two eigenfunctions. 140 Chapter 5. Conclusions and future directions A B 1 x = ρ−1 B (u) x = ρ−1 A (u) 1 0 -1 0 -1 -1 0 1 -1 0 u 1 u C x = ρ−1 C (u) 1 0 -1 -1 0 1 u Figure 5.9: The maps constructed by rough sketches of eigenfunctions. 141 Chapter 5. Conclusions and future directions Table 5.1: Eigenvalues of Smoluchowski operator for ǫ = 0.001 evaluated using different mappings. N λ2 12 14 16 18 20 22 24 26 -1.006404 -0.997546 -0.997455 -0.996985 -0.996982 λ3 λ4 λ5 -2.47368 -2.10652 -2.00468 -1.98875 -1.98789 λ6 ρA (x) -3.41230 -3.20082 -2.99653 -2.97411 -2.97266 -2.97265 λ7 λ8 λ9 λ10 -7.30289 -5.19971 -4.31843 -4.01143 -3.95399 -3.95118 -8.52820 -6.26385 -5.33439 -4.99822 -4.92765 -4.92340 -4.92337 -4.92340 ρB (x) 20 24 28 32 36 40 44 48 -2.52936 -1.99579 -2.00692 -1.98212 -1.98927 -1.98806 -1.98789 -2.54618 -2.03057 -2.03016 -1.99033 -1.98852 -1.98805 -1.98789 -4.34540 -4.65541 -3.98030 -3.99228 -3.95093 -3.95217 -3.95124 -3.95118 -4.59095 -4.74760 -3.99223 -4.00816 -3.95251 -3.95238 -3.95139 -3.95118 142 Chapter 5. Conclusions and future directions Table 5.2: Eigenvalues of Smoluchowski operator for ǫ = 0.001 evaluated using different mappings. N λ2 λ3 λ4 λ5 24 28 32 36 40 44 48 52 56 60 64 -1.021761 -1.001783 -0.997608 -0.997030 -0.996984 -0.996982 -2.40046 -2.07855 -2.00221 -1.99828 -1.98830 -1.98797 -1.98817 -1.98807 -1.98788 -1.98789 -2.36432 -2.12363 -1.99835 -2.00395 -1.98911 -1.98796 -1.98799 -1.98792 -1.98789 -2.86485 -2.32266 -2.08958 -2.01674 -1.99162 -1.98829 -1.98790 -1.98789 40 48 56 64 72 80 88 96 104 112 120 128 -1.009628 -0.998489 -0.997096 -0.996987 -0.996982 -2.60022 -2.59873 -2.03357 -2.02917 -2.00925 -1.98808 -1.98949 -1.98741 -1.98789 -2.63149 -2.57044 -2.02456 -2.03092 -2.00198 -1.98807 -1.98834 -1.98792 -1.98789 -2.40029 -2.08817 -2.00389 -1.98931 -1.98796 -1.98789 λ6 ρC (x) -3.70659 -3.24054 -3.03982 -2.98359 -2.97353 -2.97266 -2.97265 unmapped -3.39238 -3.07145 -2.98787 -2.97401 -2.97272 -2.97265 λ7 λ8 λ9 λ10 -4.41502 -4.72207 -4.03080 -3.98455 -3.97021 -3.95165 -3.95134 -3.95143 -3.95125 -3.95119 -3.95118 -4.67242 -4.69602 -4.10496 -3.97918 -3.97626 -3.95320 -3.95132 -3.95137 -3.95124 -3.95119 -3.95118 -8.34631 -6.04777 -4.85710 -4.27327 -4.03525 -3.96388 -3.95221 -3.95121 -3.95118 -9.07370 -7.12278 -5.69441 -5.14256 -4.96543 -4.92591 -4.92318 -4.92337 -4.92339 -4.92340 -7.66165 -4.53141 -4.56529 -4.21517 -3.96703 -3.97783 -3.95461 -3.95144 -3.95139 -3.95118 -3.95117 -3.95118 -7.58504 -4.52771 -4.59437 -4.18406 -3.96886 -3.97444 -3.95325 -3.95140 -3.95125 -3.95118 -6.01067 -4.73743 -4.16818 -3.98992 -3.95542 -3.95138 -3.95118 -6.96642 -5.69731 -5.13911 -4.96098 -4.92706 -4.92359 -4.92340 143 Chapter 5. Conclusions and future directions associated with ρB (x) are concentrated in two different regions, while grid points associated with ρA (x) are only concentrated in a single region. In order to have grid points concentrated in all three regions x = 0, ±1, the map ρC (x) is found by using ψˆC = ψˆA + ψˆB , shown in Fig. 5.8C, which has characteristics in all three regions. In this case, all eigenvalues can be found. It requires more grid points than the previous two cases because grid points are distributed in three regions. Nevertheless, much fewer grid points are needed for the map ρC (x) than without mapping, as shown in Table 5.2. Kramers equation It can be easily verified that the equilibrium solution of Kramers equation in Eq. (1.6) is Peq (x, y) = Kxy exp − y2 + U (x) 2 1 ǫ where Kxy is the normalization constant. Since Peq (x, y) can be written as Peq (x, y) = K exp(−y 2 /2ǫ)P0 (x) for some constant K, where P0 (x) = Kx exp(−U (x)/ǫ) for some constant Kx is the equilibrium solution of the Smoluchowski equation, the solution P (x, y) can be expanded in the same way as in Smoluchowski equation for x, and in scaled Hermite functions orthogonal with respect to the weight function w(y) = exp(−y 2 /2ǫ) for y. By setting P (x, y) = P0 (x)p(x, y), Eq. (1.6) becomes ∂p P′ ∂p ∂ ∂p ∂2p = − 0 yp − y + ν (yp) + U ′ + νǫ 2 = LK p ∂t 2P0 ∂x ∂y ∂y ∂y (5.9) We consider the expansion of p(x, y) as Nx Ny p(x, y) = i=1 j=1 Ci (x) CjHerm (y) aij √ ηi η Herm (5.10) j where aij = ηi ηjHerm p(xi , yj ), ηjHerm = wj /w(yj ), CjHerm (y) are the scaled Hermite interpolating functions CHerm (y) = j w(y) w(yj ) N k=1 k=j y − yk , yj − yk 144 Chapter 5. Conclusions and future directions w(y) = exp(−y 2 /2ǫ), yj and wj are the quadrature points and weights associated with the scaled Hermite polynomials orthogonal with respect to w(y). For the eigenvalue problem LK ψ(x, y) = λψ(x, y), I again use ψˆA (x), ψˆB (x) and ψˆC (x) from the Smoluchowski equation, Fig. 5.8, for constructing the maps, in the x-direction, ρA (x), ρB (x) and ρC (x), respectively, Fig. 5.9. In this example, the arclength monitor function is used, M (x; α) = 1 + ˜ i+1 ) − ψ(x ˜ i) ψ(x α xi+1 − xi 2 1/2 for xi < x < xi+1 and i = 1, . . . , N − 1. The parameters are α = 8, β = 10−4 and Nfit = 12. Table 5.3 shows the number of grid points, in both the x and the ydirections, required for the mapped and unmapped results to obtain the eigenvalues with all significant figures shown in the table. The eigenvalues are sorted by the real parts. In all mapped cases, significantly fewer grid points are required to obtain the eigenvalues at the same accuracy than without mapping. 5.2 Summary of results The objective of the thesis was to demonstrate a unified approach on pseudospectral methods for the solutions to the Fokker-Planck equation, Eq. (2.1), and the Schr¨ odinger equation, Eq. (2.3). Both the polynomial and the sinc methods have been extensively discussed in the previous chapters. For both polynomial and sinc collocation (SCM) methods, a mapping technique was introduced to improve the performance of pseudospectral methods. There are numerous pseudospectral methods widely used in quantum chemistry: the discrete variable representation (DVR) introduced by Light and coworkers [7], the Lagrange meshes introduced by Baye and coworkers [3], and the quadrature discretization method (QDM) introduced by Shizgal and coworkers [11]. Even though the origins and the philosophies of these methods are different, their implementations are basically the same. They are all based on orthogonal basis functions which can also serve as interpolating functions, as discussed in Chapters 2 – 4. The grid points for 145 Chapter 5. Conclusions and future directions Table 5.3: The number of grid points required to obtain the eigenvalues correct to all significant figures given in the first column. The first and second numbers in parentheses represents the number of grid points in x and y-directions, respectively. λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 eigenvalue −0.501144 ± 1.320455i −0.501144 ± 1.320455i −0.617431 −0.998152 −0.998152 −1.00554 ± 2.63851i −1.00554 ± 2.63851i −1.23291 ρA (x) ρB (x) (40, 10) (40, 12) (32, 20) (40, 12) (48, 14) (40, 14) (40, 14) (36, 20) ρC (x) (56, 10) (56, 10) (76, 20) (60, 16) (60, 16) (52, 14) (52, 14) (84, 20) unmapped (104, 10) (104, 12) (112, 20) (112, 12) (112, 12) (104, 12) (104, 12) (112, 20) 146 Chapter 5. Conclusions and future directions interpolation coincide with the points used in quadrature formulas. This special property of the basis functions gives a very simple representation of differential operators derived by the quadrature rule, yet the expansion of the solution converges exponentially to the exact solutions. The mapping introduced in Chapter 4 preserves the same orthogonality and cardinality, and hence the mapped differential operators are closely related to the corresponding unmapped differential operators but with redistributed grid points. In Chapter 2, detailed comparisons of the QDM (polynomial method) and SCM, both without mapping, were done. The motivation of this chapter came from the overly acclaimed results for SCM demonstrated by Wei [17, 18]. I have shown that unmapped SCM, with a uniform grid, does provide exponential convergence for problems with infinite domains, but is generally much less efficient than QDM for which the weight function, and hence the quadrature points, can be tailored for specific problems. No exponential convergence for SCM is observed with semi-infinite domains, such as the electron relaxation problem discussed in Section 2.3. These observations were not discussed by Wei. Fokker-Planck equations can be converted to a Schr¨ odinger equation with a change of variable as discussed in [10, 12] and in Section 2.2. The convergence of the eigenvalues for the Fokker-Planck equations and their corresponding Schr¨ odinger equations have also been discussed in Chapter 2 for QDM. The exponential convergence for QDM was illustrated in Section 2.7. In Chapter 3, I provided for the Schr¨ odinger equation a detailed comparison of the QDM with other polynomial pseudospectral methods such as the Hermite and Laguerre polynomials. For QDM, the polynomials were generated by the weight function which is chosen to be the square of the odinger equation. In this case, ground-state eigenfunction, ψ0 (x), of the Schr¨ ψ0 (x) is represented exactly with the QDM, and the corresponding eigenvalue is exact. This representation provides a very rapid convergence rate for the lower eigenstates. However, the convergence rate of the higher eigenvalues can be improved by introducing scaling and translations of the grid points. For example, we observed that for the vibrational states of I2 , it requires less than 90 grid points with scaling and translation for the eigenvalue M = 50 to converge to machine accuracy, compared to 140 with the original grid. This is shown in Fig. 3.2. 147 Chapter 5. Conclusions and future directions In cases for which the exact ψ0 (x) is unknown, a function that approximates ψ0 (x) can be used to construct the weight function. In Section 3.5, the exact ψ0 (x) for Ar2 was not given. However, the potential in Eq. (3.12) can be nicely fitted by the Morse potential given in Eq. (3.8) with appropriate parameters. Hence, the ground-state eigenfunction for the Morse potential, Eq. (3.10), is expected to be close to ψ0 (x), and is a good choice for the construction of the weight function. The convergence rates of the eigenvalues shown in Fig. 3.3 using this weight function were considerably faster than with other classical weight functions such as the Laguerre basis. The idea of using such a substitution for the weight function was further illustrated in the example of the Woods-Saxon potential. The ground state eigenfunction for the square well potential, Eq. (3.14) was used for constructing the weight function. Despite the fact that this eigenfunction was piecewise defined with discontinuous second derivative, the resulting polynomial expansion still gave fast convergence. Chapter 4 served as an extension of Chapter 3 by considering non-linear mappings instead of linear mappings, i.e. scaling and translation. The applications focused on the vibrational states of diatoms with various potentials [2,15,16]. I have demonstrated that the exponential and the hyperbolic sine mapping can dramatically improve the convergence for highly excited loosely bound states. The eigenvalues for the six highest bound states of Cs2 shown in Table 4.3 could be calculated efficiently with the hyperbolic sine mapping applied to QDM using the Morse ground-state eigenfunction given in Eq. (4.23). The mapping gives even more accurate results for He2 , for which there is only one very loosely bound eigenstate. The hyperbolic sine mapping applied to polynomials provided four significant figures as shown in Table 4.5 with only 30 points. On the other hand, more than 400 points are required without mapping. This improvement can be explained in Fig. 4.12 by the properly redistributed grid points that is achieved with the mapping. Another issue discussed in Chapter 4 was the problem of ghost levels. Ghost levels are spurious eigenvalues found by numerical methods which are not physical eigenvalues within the spectrum of the Hamiltonian. The ghost levels are results of the inexactness of the matrix elements of the Hamiltonian when integrals are evaluated by quadrature rules [19]. The ghost levels can be eliminated, or at least moved from the discrete spectrum into the continuum, if a more accurate quadrature rule is used. I have shown that, for polynomials, ghost levels could be removed by using a sightly higher order quadrature for evaluating the matrix elements. This can be easily done 148 Chapter 5. Conclusions and future directions by merely considering the eigenvalues of a smaller submatrix at the upper left corner of the original discretized matrix operator. For the applications I have considered in Chapter 4, the required submatrices were at most two rows and two columns smaller than the original matrix. Hence, the impact of the convergence rate was negligible. The mapping of grid points discussed in Chapter 4 is equivalent to a change of variables in a differential equation. Hence, mapping is a universal technique and can be applied to many types of problems. The mapping itself is not required to be either exponential or hyperbolic sine. One can construct a mapping based on the behaviour of the eigenfunctions. For example in the Schr¨ odinger equation, the formulas given in Eqs. (4.9), (4.12), and (4.16) can be used with any type of mapping. The choice of the mapping can be based on the behaviour of the potential function. Section 5.1 introduced methods for the construction of efficient mapping for the interpolation of a function based on the equidistribution principle. The resulting grid points are placed so that some property of the function, such as the arc length, is distributed over the grid as evenly as possible. Suitable parameters of explicitly defined mapping, for example Eq. (4.17b) or (4.17c) used in Chapter 4, can be determined by the equidistribution principle. The mapping can also be constructed implicitly from a reference grid generated directly by the equidistribution principle using a least squares fit. Mapping is particularly useful when the function has localized behaviours such as the function given in Eq. (5.2). I have shown in Figs. 5.1, 5.2, and 5.4 that interpolation over mapped grids requires significantly fewer grid points than over a uniform grid. The idea of implicit mapping can also be used in time-dependent problems for which the grid moves as the function evolves over time, as shown in Figs. 5.5 and 5.6. I have also applied the implicit mapping to eigenvalue problems such as Fokker-Planck and Kramers equations. The improvements due to the mapping for time-dependent and eigenvalue problems are illustrated in Fig. 5.7 and Tables 5.1 – 5.3. 149 Chapter 5. 5.3 Conclusions and future directions Future directions The interpolation with implicit mapping based on the equidistribution principle can be applied to a wide range of functions. The quality of interpolation depends on the choice of the monitor function and the construction of the mapping function. A detailed analysis on the choice of the monitor function for different situations can be one of our future tasks. The construction of the mapping function also requires a deeper study. Using Eq. (5.5) as the form of the mapping function will lead to unwanted oscillations, as shown in Fig. 5.3A. With a smoothing parameter β, the oscillations can be diminished, but the resulting mapping function will not be as close to the original reference grid as before. This problem can be seen in Fig. 5.3B, when the function does not follow the points (crosses) closely in the flat region u ∈ [0, 1.5]. Therefore, another way to construct the mapping function could be with splines or parametric representations. In the interpolatory moving mesh method, the x-direction of the solution p(x, t) is interpolated on a grid which changes with the bahaviour of the solution. It is therefore expected that the interpolation is much more efficient than with a non-evolving grid. However, the time-integration is evaluated by backward Euler method, and the error is generally expected to be poor. Other methods such as the Runge-Kutta method can be used to improve the accuracy of the solution, but numerical stability is an important issue for forward methods. In the future, it will be important to find a method, perhaps an adaptive method, for a more accurate time integration. The moving mesh for finite difference and finite elements are usually done without reinterpolation of the solution [4, 5]. The time-dependent problem is solved simultaneously with a moving mesh partial differential equation, which is used to define new grid points at each time step. As a result, the solution is automatically defined on new grids and reinterpolation is not needed. The possibility of a similar algorithm for pseudospectral methods can be studied in the future. As shown in Tables 5.1 and 5.2, the implicit mapping for sinc method can improve the convergence rate of eigenvalues over the unmapped sinc method. However by comparing the results with Table 2.3, the QDM still provides the best convergence rate for this particular problem. We would like to study the performance of the implicit mapping with different types of basis functions. It will also be useful to apply implicit mapping on other 150 Chapter 5. Conclusions and future directions types of problems, including the vibrational energy problem discussed in Chapter 4. The simplest form of interpolation in two dimensions is to define the interpolating functions as tensor products of 1D interpolating functions. A function g(x, y) can be interpolated in the form Nx Ny g(x, y) = (x) (y) g(xi , yj )Ci (x)Cj (y) i=1 j=1 (x) (y) where Ci and Cj are 1D interpolating functions. In this case, any kind of mapping will be done separately for each variable. Two dimensional mapping will be studied in the future. The equidistribution principle for 2D adaptive grids is discussed in Ref. 6. In order to apply the adaptive grid to pseudospectral methods, a 2D mapping function is required. The construction of such a mapping function may involve Fourier expansions, thin plate splines, or parametric representations. The study of 2D mappings on moving mesh methods and eigenvalues problem will start after the development of 2D mapped interpolation. 151 Bibliography [1] U. M. Ascher, R. M. M. Mattheij, and D. Russell, Robert. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Prentice Hall, Englewood Cliffs, 1988. [2] R. A. Aziz and M. J. Slaman. An examination of ab initio results for the helium potential energy curve. J. Chem. Phys., 94:8047, 1991. [3] D. Baye. Lagrange-mesh method for quantum-mechanical problems. Phys. Stat. Sol. (b), 243:1095, 2006. [4] G. Beckett, J. A. Mackenzie, A. Ramage, and D. M. Sloan. On the numerical solution of one-dimensional PDEs using adaptive methods based on equidistribution. J. Comput. Phys., 167:372, 2001. [5] W. Huang, Y. Ren, and R. Russell. Moving mesh partial differential equations (MMPDEs) based on the equidistribution principle. SIAM J. Numer. Anal., 31:709, 1994. [6] W. Huang and D. Sloan. A simple adaptive grid method in two dimensions. SIAM J. Sci. Comput., 15:776, 1994. [7] J. C. Light and T. Carrington. Discrete variable representations and their utilization. Adv. Chem. Phys., 114:263, 2000. [8] J. Lund and K. L. Bowers. Sinc Methods for Quadrature and Differential Equations. SIAM, Philadelphia, 1992. [9] A. Madzvamuse and P. K. Maini. Velocity-induced numerical solutions of reaction-diffusion systems on continuously growing domains. J. Comput. Phys., 225:100, 2007. [10] H. Risken. The Fokker-Planck Equation: Methods of Solution and Application. Springer-Verlag, New York, 1984. [11] B. D. Shizgal and H. Chen. The quadrature discretization method (QDM) in the solution of the Schr¨ odinger equation with nonclassical basis functions. J. Chem. Phys., 104:4137, 1996. 152 Bibliography [12] B. D. Shizgal and H. Chen. The quadrature discretization method in the solution of the Fokker-Planck equation with nonclassical basis function. J. Chem. Phys., 107:8051, 1997. [13] F. Stenger. Numerical Methods Based on Sinc and Analytic Functions. Springer-Verlag, New York, 1993. [14] H. Z. Tang. A moving mesh method for the Euler flow calculations using a directional monitor function. Commun. Comput. Phys., 1:656, 2006. [15] K. T. Tang and J. P. Toennies. The van der Waals potentials between all the rare gas atoms from He to Rn. J. Chem. Phys., 118:4976, 2003. [16] K. T. Tang, J. P. Toennies, and C. L. Yiu. Accurate analytical He-He van der Waals potential based on perturbation theory. Phys. Rev. Lett., 74:1546, 1995. [17] G. W. Wei. Discrete singular convolution for the solution of the FokkerPlanck equation. J. Chem. Phys., 110:8930, 1999. [18] G. W. Wei. Solving quantum eigenvalue problems by discrete singular convolution. J. Phys. B: At. Mol. Opt. Phys., 33:343, 2000. [19] H. Wei. Ghost levels and near variational forms of the discrete variable representation: application to H2 O. J. Chem. Phys., 106:6885, 1997. 153 Appendix A Matrix Representation of the Lorentz Fokker-Planck Operator We consider {Qn (x)}∞ n=0 to be the set of polynomial orthogonal with respect to a weight function w(x) = P0 (x), i.e. w(x)Qm (x)Qn (x) dx = γn δmn , √ where Qn (x) = γn Fn (x) and {Fn (x)}∞ n=0 is a normalized set as defined in Section 2.2. These polynomials satisfy the three-term recursion formula Qn (x) = (x − αn−1 )Qn−1 (x) − βn−1 Qn−2 (x) (A.1) for n ≥ 2, where Q0 (x) = 1, Q1 (x) = x − α0 , αn = γn−1 xw(x)[Qn (x)]2 dx and βn = γn /γn−1 . Qn (x) is an nth degree monic polynomial, for which the coefficient of xn is unity. The self-adjoint Fokker-Planck operator in Eq. (2.1) is given by Lφm (x) = 1 d dφm (x) w(x)B(x) . w(x) dx dx The matrix elements of the Fokker-Planck operator, Lmn , are thus Lmn = √ ∞ 1 γm−1 γn−1 Qm−1 (x) 0 d w(x)B(x)Q′n−1 (x) dx dx (A.2) with normalization factors (γm−1 γn−1 )−1/2 in front. If the integral in Eq. (A.2) is evaluated by parts, we get Lmn = − √ ∞ 1 γm−1 γn−1 0 w(x)B(x)Q′m−1 (x)Q′n−1 (x) dx, (A.3) where the boundary terms are zero. Since L is self-adjoint, the matrix L is symmetric as is clear from Eq. (A.3). For the specific case of the Lorentz Fokker-Planck operator, A(x) = 2x2 − 3 and B(x) = x in Eq. (2.1), we have that w(x) = x2 exp(−x2 ). 154 Appendix A. Lorentz Fokker-Planck Operator Taking the derivative of the term in square bracket in Eq. (A.2), and using xw′ (x) = 2(1 − x2 )w(x), we get that Lmn = √ ∞ 3 γm−1 γn−1 −√ 0 2 w(x)Qm−1 (x)Q′n−1 (x) dx ∞ x2 w(x)Qm−1 (x)Q′n−1 (x) dx γm−1 γn−1 0 ∞ 1 +√ xw(x)Qm−1 (x)Q′′n−1 (x) dx. (A.4) γm−1 γn−1 0 It is sufficient to consider only the matrix elements for m > n because L is symmetric. Since both Q′n−1 (x) and xQ′′n−1 (x) are polynomials of degree n − 2 < m − 1, the first and third integrals in Eq. (A.4) vanish. To evaluate the second integral, we expand x2 Q′n−1 in the Qn basis, thus we write x2 Q′n−1 (x) = bn Qn (x) + bn−1 Qn−1 (x) + · · · + b0 Q0 (x). (A.5) Hence, this integral vanishes when m > n + 1. For m = n + 1, √ −2 γn γn−1 ∞ 0 x2 w(x)Qn (x)Q′n−1 (x) dx −2bn =√ γn γn−1 ∞ 0 w(x)[Qn (x)]2 dx = −2bn βn . (A.6) Since Qn (x) is a monic polynomial of degree n, the coefficient of xn in x2 Q′n−1 (x) is n − 1. So bn = n − 1. For m = n, the integral is given by −2 γn−1 ∞ 0 x2 w(x)Qn−1 (x)Q′n−1 (x) dx = −2bn−1 γn−1 ∞ 0 w(x)[Qn−1 (x)]2 dx = −2bn−1 . (A.7) We wish to express bn−1 in terms of the αk and βk in the three-term recurrence relation in Eq. (A.1). We first express Qn (x) as the polynomial n k Qn (x) = k=0 An,k x , where the first three coefficients are An,n = 1, n−1 k−1 n−1 An,n−1 = − n−1 k=0 αk and An,n−2 = k=0 l=0 αk αl − m=1 βm , so that from the expansion in Eq. (A.5), x2 Q′n−1 (x) = (n − 1)xn + (n − 2)An−1,n−2 xn−1 + · · · + An−1,1 x2 , bn Qn (x) = bn xn + bn An,n−1 xn−1 + · · · + bn An,0 , bn−1 Qn−1 (x) = bn−1 xn−1 + · · · + bn−1 An−1,0 . 155 Appendix A. Lorentz Fokker-Planck Operator Equating the coefficients of xn−1 we get bn−1 = (n − 2)An−1,n−2 − bn An,n−1 n−2 = −(n − 2) k=0 n−1 αk + (n − 1) αk k=0 n−2 = (n − 1)αn−1 + αk . k=0 It is understood that the sum is not considered when m = n = 1. Hence the entries of the symmetric operator L for the Lorentz Fokker-Planck equation is given by −2(n − 1)αn−1 − 2 n−2 k=0 αk when m = n > 1, −2(n − 1)√β when m = n + 1, n Lmn = (A.8) √ −2(m − 1) β when m = n − 1, m 0 otherwise. The matrix operator L 0 0 0 −4.9761 0 −1.3032 0 0 L= 0 0 0 0 .. .. . . is given explicitly by 0 0 0 −1.3032 0 0 −11.258 −3.1214 0 −3.1214 −18.689 −5.3220 0 0 .. . 0 0 0 0 ··· ··· ··· ··· . −5.3220 −27.137 −7.8425 . . . 0 −7.8425 −36.505 . . .. .. .. .. . . . . . (A.9) It is clear that the off-diagonal elements of this tridiagonal matrix are smaller than the diagonal elements and rapid convergence of the eigenvalues is obtained [2]. The lower order eigenvalues, λ1 = 4.6834, λ2 = 10.113, and λ3 = 16.430, are close to the diagonal elements. 156 Appendix B Matrix Representation of the Fokker-Planck Operator for the Bistable System For the bistable system, A(x) = x3 − x and B(x) = ǫ in Eq. (2.1). It follows that the weight function is w(x) = exp(−(x4 /4 − x2 /2)/ǫ), and w′ (x) = −(x3 − x)w(x)/ǫ. Since w(x) is even and the integrals are evaluated over (−∞, ∞), the coefficients αn are zero. Thus the polynomials Qn (x) are even when n is an even number, and odd when n is an odd number. We use the same notation as in Appendix A, but these are different polynomials. The matrix elements in Eq. (A.2) are 1 Lmn = − √ γm−1 γn−1 ∞ −∞ (x3 − x)w(x)Qm−1 (x)Q′n−1 (x) dx +√ ∞ ǫ γm−1 γn−1 −∞ w(x)Qm−1 (x)Q′′n−1 (x) dx. As in Appendix A, it is necessary to consider only m ≥ n. The second integral always vanishes because Q′′n−1 (x) is a polynomial of degree n − 3 < m − 1. For the first integral, the expansion in the Qn basis, (x3 − x)Q′n−1 (x) = bn+1 Qn+1 (x) + bn−1 Qn−1 (x) + bn−3 Qn−3 (x) + ··· + b1 Q1 (x) for even n, b0 Q0 (x) for odd n (B.1) 157 Appendix B. Fokker-Planck Operator for the Bistable System is used. If m = n + 2, then −√ ∞ 1 γn+1 γn−1 −∞ (x3 − x)w(x)Qn+1 (x)Q′n−1 (x) dx ∞ −bn+1 =√ w(x)[Qn+1 (x)]2 dx γn+1 γn−1 −∞ −bn+1 =√ γn+1 = −bn+1 βn+1 βn . γn+1 γn−1 By the same argument as in Appendix A below Eq. (A.6), the coefficient of xn+1 in (x3 − x)Q′n−1 (x) is n − 1, and therefore bn+1 = n − 1. If m = n, then − 1 ∞ γn−1 −∞ (x3 − x)w(x)Qn−1 (x)Q′n−1 (x) dx = −bn−1 γn−1 ∞ −∞ w(x)[Qn−1 (x)]2 dx = −bn−1 . To express the coefficient bn−1 in terms of βk , the expansion of Eq. (B.1) is used. The representation of Qn is given below Eq. (A.7) in Appendix A. (x3 − x)Q′n−1 (x) = (x3 − x) (n − 1)xn−2 + (n − 3)An−1,n−3 xn−4 + · · · = (n − 1)xn+1 + [(n − 3)An−1,n−3 − (n − 1)]xn−1 + · · · bn+1 Qn+1 (x) = bn+1 xn+1 + bn+1 An+1,n−1 xn−1 + · · · bn−1 Qn−1 (x) = bn−1 xn−1 + · · · Equating the coefficients of xn−1 gives bn−1 = [(n − 3)An−1,n−3 − (n − 1)] − bn An+1,n−1 n−2 = −(n − 3) k=1 n βk − (n − 1) + (n − 1) βk k=1 n−2 = (n − 1)(βn + βn−1 − 1) + 2 βk . (B.2) k=1 158 Appendix B. Fokker-Planck Operator for the Bistable System The sum is ignored when n ≤ 2. Hence, the symmetric operator L is given by when m = n = 2, −(β2 + β1 − 1) n−2 −(n − 1)(βn + βn−1 − 1) − 2 k=1 βk when m = n > 2, Lmn = −(n − 1) βn+1 βn when m = n + 2, (B.3) −(m − 1) βm+1 βm when m = n − 2, 0 otherwise, where the βk recursion coefficients depend on ǫ as the weight function w(x) depends on ǫ. The matrix operator L for ǫ = 0.1 is given explicitly by 0 0 0 0 0 0 ··· 0 −0.0101 0 −0.1407 0 0 ··· 0 0 −1.9577 0 −0.4024 0 ··· .. 0 −0.1407 . 0 −2.0518 0 −0.5967 . L= .. 0 . 0 −0.4024 0 −3.9122 0 .. 0 . 0 0 −0.5967 0 −4.0964 .. .. .. .. .. .. .. . . . . . . . (B.4) The matrix L is pentadiagonal, with zero off-diagonal entries. This result was previously reported by Blackmore and Shizgal [1] but the derivation was not provided. As discussed in this previous paper, these eigenfunctions of odd and even parity do not couple. The matrix defined by Eq. (B.3) and 159 Appendix B. Fokker-Planck Operator for the Bistable System Eq. (B.4) can be split into two matrices which are explicitly 0 0 0 0 0 ··· 0 −1.9577 −0.4024 0 0 ··· 0 −0.4024 −3.9122 −0.9861 0 · · · .. Leven = 0 , . 0 −0.9861 −5.8627 −1.7090 . 0 0 0 −1.7090 −7.8079 . . .. .. .. .. .. .. . . . . . . −0.0101 −0.1407 0 0 0 ··· −0.1407 −2.0518 −0.5967 0 0 ··· 0 −0.5967 −4.0964 −1.2175 0 ··· . Lodd = 0 0 −1.2175 −6.1445 −1.9673 . . . 0 0 0 −1.9673 −8.1968 . . .. .. .. .. .. .. . . . . . . (B.5) . (B.6) The spectrum of the Fokker-Planck operator in terms of a set of singlet and nearly degenerate triplets is clear from the structure of these matrices. 160 Appendix C Matrix Representation of the Schr¨ odinger Operator If the weight function w(y) satisfies V (y) = V0 (y) as defined in Eq. (2.4), the self-adjoint Schr¨ odinger operator in Eq. (2.3) is given by Hψm (y) = − 1 d d w(y) dy dy w(y) ψn (y) w(y) , where the form is recognized as the self-adjoint form of a Fokker-Planck operator with a diffusion coefficient equal to unity. Shizgal and Chen [3] discussed the transformation of the Schr¨ odinger equation to a corresponding Fokker-Planck equation. When the basis { w(y)Qn (y)}∞ n=0 , where Qn are orthogonal polynomials satisfying Eq. (A.1), are used, the matrix elements of the Schr¨ odinger operator, Hmn are given by Hmn = √ ∞ 1 γm−1 γn−1 0 w(y)Q′m−1 (y)Q′n−1 (y) dy (C.1) after evaluated by parts with vanishing boundary terms. To evaluate this integral, we expand Q′n−1 in the Qn basis, which gives Q′n−1 (x) = Cn,0 Q0 (x) + Cn,1 Q1 (x) + · · · + Cn,n−2 Qn−2 (x) for n ≥ 2. With the orthogonality relation, the integral in Eq. (C.1) is given by Hmn = √ 1 γm−1 γn−1 min(m−2,n−2) Cm,k Cn,k (C.2) k=0 for m ≥ 2 and n ≥ 2. The matrix elements vanish for m = 1 or n = 1. In principle, Cm,k can be expressed in terms of αk and βk , but the expressions will be too complicated for practical use. The matrix operator H for the electron relaxation problem is evaluated with a quadrature integration of 161 Appendix C. Matrix Representation of the Schr¨ odinger Operator Eq. (C.1) and is given 0 0 0 4.7255 0 0.4163 H = 0 0.3626 0 −0.1302 0 0.0750 .. .. . . by 0 0.4163 10.246 0.6053 1.1485 −0.4904 .. . 0 0 0 ··· 0.3626 −0.1302 0.0750 · · · 0.6053 1.1485 −0.4904 · · · 16.794 0.3615 2.4556 · · · . (C.3) 0.3615 24.491 −0.4420 · · · 2.4556 −0.4420 33.408 · · · .. .. .. .. . . . . The matrix is diagonally dominant and rapid convergence of the eigenvalues analogous to the rate with the matrix Eq. (A.9) is obtained and shown in Fig. 2.1. The matrix representation of the Hamiltonian operator for the bistable system which yields the even eigenfunctions with the half-range basis functions orthogonal with weight function wbe (y) = P0 (y), y ∈ [0, ∞), for ǫ = 0.01 is given by Heven = 0 0 0 0 1.9204 0.3054 0 0.3054 3.7174 0 0.0137 0.7557 0 −5.8(−8) 0.0412 0 2.8(−7) −5.2(−6) .. .. .. . . . 0 0 0 ··· 0.0137 −5.8(−8) 2.8(−7) · · · 0.7557 0.0412 −5.2(−6) · · · 5.3693 1.3264 0.0866 ··· . 1.3264 6.8363 2.0169 ··· 0.0866 2.0169 8.0068 ··· .. .. .. .. . . . . (C.4) The corresponding matrix representative of the Hamiltonian operator for the bistable system which yields the odd eigenfunctions with the half-range basis functions orthogonal with weight function wbo (y) = y 2 P0 (y), y ∈ [0, ∞), 162 Appendix C. Matrix Representation of the Schr¨ odinger Operator for ǫ = 0.01 is given by Hodd 0.0101 0.1405 0.0072 0.1405 1.9960 0.4987 0.0072 0.4987 3.8705 0 0.0263 0.9835 = 0 −5.1(−5) 0.0588 0 1.0(−5) −2.0(−4) .. .. .. . . . 0 0 0 ··· 0.0263 −5.1(−5) 1.0(−5) · · · 0.9835 0.0588 −2.0(−4) · · · 5.6186 1.5740 0.1078 ··· . 1.5740 7.2185 2.2626 ··· .1078 2.2626 8.6342 ··· .. .. .. .. . . . . (C.5) The matrix elements are calculated numerically with the integration of Eq. (C.1). It is clear that both matrices are diagonally dominant and the rapid convergence of the eigenvalues shown in Table 2.2 is understood. The matrix representative of the Hamiltonian for vibrational states of I2 modeled with a Morse potential and the basis functions defined by the weight function in Eq. (C.1) is given by (H − λ0 I) × 106 0 0 0 0 0 0 ··· 0 568.67 −32.293 1.4997 −0.0604 0.0022 · · · 0 −32.293 1135.5 −78.973 4.2384 −0.1911 · · · = 0 1.4997 −78.973 1700.5 −136.56 8.2008 · · · 0 −0.0604 4.2384 −136.56 2263.7 −203.25 · · · 0 0.0022 −0.1911 8.2008 −203.25 2825.0 · · · .. .. .. .. .. .. .. . . . . . . . in atomic units. As in the other model systems, this matrix representative is diagonally dominant and the rapid convergence is anticipated. 163 Appendix D Moving mesh method for the time dependent Smoluchowski equation The initial condition of the Smoluchowski equation, Eq. (5.6), is interpolated with the implicit mapping generated based on the monitor function in Eq. (5.3). The maps, the physical grids, and the interpolating functions [n] [n] used at subsequent time steps tn are denoted by ρ[n] (x), xj , and Cj (x), respectively. The discretization of Eq. (5.6) gives the Smoluchowski operator matrix LS [n] [ij] = − ǫ ν N (1) [n] [n] Dki ρ′[n] (xk ) + gˆ[n] (xk )δki k=1 (1) [n] [n] × Dkj ρ′[n] (xk ) + gˆ[n] (xk )δkj , (D.1) (1) where Dij is given by Eq. (2.6) and gˆ[n] (x) = 1 2 ρ′′[n] (x)/ρ′[n] (x) + P0′ (x)/P0 (x) . [n] With the quadrature weights ηj , the coefficients aj (t) from Eq. (5.7) asso[n] [n] ciated with xj are denoted by aj (t) = the expansion N p(x, t) = [n] [n] [n] aj (t) j=1 [n] ηj p(xj , t). Therefore, by using Cj (x) [n] ηj , [n] the Smoluchowski equation in Eq. (5.6) becomes daj /dt = Based on the backward difference method, the coefficients ated by solving [n] N ai (tn ) = j=1 [n] δij − ∆tn LS[n] [ij] aj (tn+1 ) [n] N S j=1 L[n] [ij]aj . [n] aj (t) are evalu- (D.2) 164 Appendix D. Moving mesh method for the time dependent Smoluchowski equation [n] for aj (tn+1 ), where ∆tn = tn − tn−1 . If a new grid is used at tn+1 , i.e. [n+1] [n] [n+1] xj = xj , aj (tn+1 ) has to be found. The solution of Eq. (5.6) and (1.7) is therefore given by N p(x, tn ) ≈ and P (x, tn ) = j=1 [n] aj (tn ) [n] ηj [n] Cj (x). P0 (x)p(x, tn ), respectively. 165 Bibliography [1] R. Blackmore and B. Shizgal. A solution of Kramers equation for the isomerization of n-butane in CCl4 . J. Chem. Phys., 83:2934, 1985. [2] H. D. Meyer, J. Hor´ acek, and L. S. Cederbaum. Schwinger and anomalyfree Kohn variational principles and a generalized Lanczos algorithm for nonsymmetric operators. Phys. Rev. A, 43:3587, 1991. [3] B. D. Shizgal and H. Chen. The quadrature discretization method in the solution of the Fokker-Planck equation with nonclassical basis function. J. Chem. Phys., 107:8051, 1997. 166
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Pseudospectral methods in quantum and statistical mechanics
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Pseudospectral methods in quantum and statistical mechanics Lo, Joseph Quin Wai 2008
pdf
Page Metadata
Item Metadata
Title | Pseudospectral methods in quantum and statistical mechanics |
Creator |
Lo, Joseph Quin Wai |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | The pseudospectral method is a family of numerical methods for the solution of differential equations based on the expansion of basis functions defined on a set of grid points. In this thesis, the relationship between the distribution of grid points and the accuracy and convergence of the solution is emphasized. The polynomial and sinc pseudospectral methods are extensively studied along with many applications to quantum and statistical mechanics involving the Fokker-Planck and Schroedinger equations. The grid points used in the polynomial methods coincide with the points of quadrature, which are defined by a set of polynomials orthogonal with respect to a weight function. The choice of the weight function plays an important role in the convergence of the solution. It is observed that rapid convergence is usually achieved when the weight function is chosen to be the square of the ground-state eigenfunction of the problem. The sinc method usually provides a slow convergence as the grid points are uniformly distributed regardless of the behaviour of the solution. For both polynomial and sinc methods, the convergence rate can be improved by redistributing the grid points to more appropriate positions through a transformation of coordinates. The transformation method discussed in this thesis preserves the orthogonality of the basis functions and provides simple expressions for the construction of discretized matrix operators. The convergence rate can be improved by several times in the evaluation of loosely bound eigenstates with an exponential or hyperbolic sine transformation. The transformation can be defined explicitly or implicitly. An explicit transformation is based on a predefined mapping function, while an implicit transformation is constructed by an appropriate set of grid points determined by the behaviour of the solution. The methodologies of these transformations are discussed with some applications to 1D and 2D problems. The implicit transformation is also used as a moving mesh method for the time-dependent Smoluchowski equation when a function with localized behaviour is used as the initial condition. |
Extent | 1658867 bytes |
Subject |
Pseudospectral methods Fokker-Planck equation Schroedinger equation Interpolation |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-08-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066507 |
URI | http://hdl.handle.net/2429/1298 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_fall_lo_quin_wai.pdf [ 1.58MB ]
- Metadata
- JSON: 24-1.0066507.json
- JSON-LD: 24-1.0066507-ld.json
- RDF/XML (Pretty): 24-1.0066507-rdf.xml
- RDF/JSON: 24-1.0066507-rdf.json
- Turtle: 24-1.0066507-turtle.txt
- N-Triples: 24-1.0066507-rdf-ntriples.txt
- Original Record: 24-1.0066507-source.json
- Full Text
- 24-1.0066507-fulltext.txt
- Citation
- 24-1.0066507.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066507/manifest