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 so- lution 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 ex- tensively studied along with many applications to quantum and statistical mechanics involving the Fokker-Planck and Schrödinger 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 dis- tributed 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 dis- cussed in this thesis preserves the orthogonality of the basis functions and provides simple expressions for the construction of discretized matrix op- erators. 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 im- plicit 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 prob- lems. 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 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . x Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 The Fokker-Planck equation . . . . . . . . . . . . . . . . . . 3 1.3 The Schrödinger equation . . . . . . . . . . . . . . . . . . . . 6 1.4 Overview of current numerical methods . . . . . . . . . . . . 7 1.5 Interpolation and quadrature . . . . . . . . . . . . . . . . . . 9 1.5.1 Orthogonal polynomials . . . . . . . . . . . . . . . . . 10 1.5.2 Sinc functions . . . . . . . . . . . . . . . . . . . . . . 13 1.5.3 Transformation of grid points . . . . . . . . . . . . . 14 1.6 Pseudospectral method . . . . . . . . . . . . . . . . . . . . . 15 1.7 Thesis objective . . . . . . . . . . . . . . . . . . . . . . . . . 17 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2 Spectral convergence of quadrature discretization method 28 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2 Discretization of the Fokker-Planck and Schrödinger eigen- value problems . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3 Electron relaxation . . . . . . . . . . . . . . . . . . . . . . . 35 2.4 Vibrational states of I2 . . . . . . . . . . . . . . . . . . . . . 36 2.5 The quartic bistable system . . . . . . . . . . . . . . . . . . . 37 iv Table of Contents 2.6 Convergence of eigenvalues; discussion of results . . . . . . . 38 2.7 Interpretation of the numerical results . . . . . . . . . . . . . 52 2.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3 Pseudospectral methods of solution of the Schrödinger equa- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2 Spectral convergence . . . . . . . . . . . . . . . . . . . . . . 67 3.3 Pseudospectral method . . . . . . . . . . . . . . . . . . . . . 69 3.4 The vibrational states of the Morse oscillator . . . . . . . . . 71 3.5 The vibration states of Ar2 . . . . . . . . . . . . . . . . . . . 73 3.6 Woods-Saxon potential . . . . . . . . . . . . . . . . . . . . . 74 3.7 Results and discussions . . . . . . . . . . . . . . . . . . . . . 74 3.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4 An efficient mapped pseudospectral method for weakly bound states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.2 Solution of the Schrödinger 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ödinger equation with the sinc func- tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 Interpolatory moving mesh . . . . . . . . . . . . . . . 132 5.1.3 Eigenvalues of the Smoluchowski and Kramers equa- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 5.2 Summary of results . . . . . . . . . . . . . . . . . . . . . . . 145 5.3 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . 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ödinger Operator . . . . 161 D Moving mesh method for the time dependent Smoluchowski equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 vi List of Tables 2.1 Convergence of the eigenvalues for the quartic potential, ǫ = 0.1. 45 2.2 Convergence of the eigenvalues for the quartic potential, ǫ = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.3 Convergence of the eigenvalues for the quartic potential, ǫ = 0.001. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1 Convergence of 103E10 for the I2 Morse potential. . . . . . . . 75 4.1 The parameters of the Tang-Toennies potentials for He2, Ne2, and Ar2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.2 The parameters of the Aziz-Slaman LM2M2 potential for He2. 103 4.3 Convergence of the Morse potential for Cs2. . . . . . . . . . . 106 4.4 The vibrational energy levels for He2, Ne2 and Ar2. . . . . . . 110 4.5 Convergence of the energy for Aziz-Slaman He2 potential. . . 117 5.1 Eigenvalues of Smoluchowski operator evaluated using differ- ent mappings . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.2 Eigenvalues of Smoluchowski operator evaluated using differ- ent mappings . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 5.3 The number of grid points required to obtain the eigenvalues 146 vii List of Figures 2.1 The relative error of eigenvalues versus the number of quadra- ture points for the electron thermalization problem. . . . . . 40 2.2 The relative error of eigenvalues versus the number of quadra- ture points for the I2 Morse potential. . . . . . . . . . . . . . 42 2.3 Variation of the potential in the Schrödinger equation for the quartic bistable system. . . . . . . . . . . . . . . . . . . . . . 43 2.4 The relative error of eigenvalues versus the number of quadra- ture points for the bistable system with ǫ = 0.01. . . . . . . . 48 2.5 The relative error of eigenvalues versus the number of quadra- ture points for the bistable system with ǫ = 0.001. . . . . . . 49 2.6 Variation of the expansion coefficients of the eigenfunctions for the electron relaxation problem. . . . . . . . . . . . . . . . 54 3.1 The expansion coefficients of the eigenfunctions of the I2 system. 68 3.2 The relative error of eigenvalues versus the number of quadra- ture points for the I2 Morse potential. . . . . . . . . . . . . . 77 3.3 The relative error of eigenvalues versus the number of quadra- ture points for the Ar2 Cahill and Parsegian potential. . . . . 79 3.4 The relative error of eigenvalues versus the number of quadra- ture points for the Woods-Saxon potential. . . . . . . . . . . 81 4.1 The mapping of collocation points from Laguerre grid to a modified grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.2 The He2, Ne2, and Ar2 potentials. . . . . . . . . . . . . . . . 99 4.3 The zero energy limit of the l = 0 phase shift for He2, Ne2, and Ar2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.4 The eigenvalues and the ghost level versus the number of collocation points for Cs2 Morse potential. . . . . . . . . . . . 108 4.5 The expansion coefficients and the eigenfunction of the ghost level for Cs2 Morse potential. . . . . . . . . . . . . . . . . . . 108 4.6 The ground state eigenfunction for Ne2. . . . . . . . . . . . . 110 4.7 The highest state eigenfunction of Ne2. . . . . . . . . . . . . . 111 viii List of Figures 4.8 The relative error of eigenvalues versus the number of collo- cation points for Ne2. . . . . . . . . . . . . . . . . . . . . . . 112 4.9 Variation of the expansion coefficients of the highest state eigenfunction of Ne2. . . . . . . . . . . . . . . . . . . . . . . . 114 4.10 The relative error of eigenvalues versus the number of collo- cation points for Ar2. . . . . . . . . . . . . . . . . . . . . . . . 115 4.11 The relative error of the eigenvalue versus the number of col- location points for He2. . . . . . . . . . . . . . . . . . . . . . 118 4.12 The eigenfunction of the Aziz-Slaman He2 potential. . . . . . 119 4.13 The eigenvalue and the ghost level versus the number of col- location points for Aziz-Slaman He2 potential. . . . . . . . . . 121 5.1 Unmapped sinc interpolation . . . . . . . . . . . . . . . . . . 127 5.2 Sinc interpolation with hyperbolic sine mapping . . . . . . . . 129 5.3 The reference grid and the corresponding mapping function . 131 5.4 Sinc interpolation with implicitly defined mapping . . . . . . 133 5.5 Solution of the Smoluchowski equation by interpolatory mov- ing mesh method . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.6 The evolution of grid points over time for the solution of Smoluchowski equation . . . . . . . . . . . . . . . . . . . . . . 137 5.7 The error of the solution of Smoluchowski equation for inter- polatory moving mesh method and unmapping sinc method . 137 5.8 Eigenfunctions of Smoluchowski equation solved by unmapped sinc method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 5.9 The maps constructed by rough sketches of eigenfunctions . . 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 grad- uate 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ödinger 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 nu- merical 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 ∂t = 1 k ∂2u ∂x2 (1.1) 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 ∂t = 1 k ∇2u, where ∇2 = ∂2/∂x2 + ∂2/∂y2 in 2D or ∂2/∂x2 + ∂2/∂y2 + ∂/∂z2 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 ∂y = ∇ · (D(x, u)∇u), (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 ∂u ∂t = D ∂2u ∂x2 + f(u), 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 influ- enced by other components in the system. The resulting reaction-diffusion equation for each of the component will contain the reaction term, or inter- action term, f , depending on all components. A two-component, u1(x, t) and u2(x, t), reaction-diffusion system has the form ∂u1 ∂t = D1∇2u1 + f1(u1, u2) ∂u2 ∂t = D2∇2u2 + f2(u1, u2) 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 inter- action terms of the Lotka-Volterra system, f1(u1, u2) = a1u1(1−b1u1−c1u2) and f2(u1, u2) = a2u2(1− b2u2− c2u1), model the logistic growth with com- petition. 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 ∂t + u ∂u ∂x = ν ∂2u ∂x2 , (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, ν(∂2u/∂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 solu- tions to the Fokker-Planck equation and the Schrödinger equation. Since the Fokker-Planck equation, the Schrödinger 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 proba- bility 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änggi 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 d2x dt2 + ν dx dt = F (x) + Γ(t), (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 dt = v, (1.5a) dv dt = −νv + F (x) + Γ(t), (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 ∂t = −v∂P ∂x − F (x)∂P ∂v + ν ∂(vP ) ∂v + νǫ ∂2P ∂v2 . (1.6) 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 prod- ucts [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 dt = 1 ν [F (x) + Γ(t)]. The probability density function P (x, t) = ∫∞ −∞ P (x, v, t) dv satisfies the Smoluchowski equation [81] ∂P ∂t = 1 ν [ −∂(FP ) ∂x + ǫ ∂2P ∂x2 ] . (1.7) A general 1D Fokker-Planck equation is given by [81] ∂ ∂t P (x, t) = ∂ ∂x [A(x)P (x, t)] + ∂2 ∂x2 [B(x)P (x, t)]. (1.8) The functions A(x) and B(x) are the drift and diffusion coefficients, re- spectively. The Smoluchowski equation in Eq. (1.7) is an example of a 1D Fokker-Planck equation. Other uses of the Fokker-Planck equation in- clude 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 B(x) exp ( − ∫ x A(x′) B(x′) dx′ ) . (1.9) with a normalization factor Kx. By setting P (x, t) = √ P0(x)p(x, t), Eq. (1.8) will have the self-adjoint form ∂p ∂t = 1√ P0 ∂ ∂x [ BP0 ∂ ∂x ( p√ P0 )] = LFPp. (1.10) A classical approach to the solution of a linear parabolic partial differen- tial 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 eigenfunc- tions, ψm(x), of the Fokker-Planck operator L FP, the solution p(x, t) is given by p(x, t) = ∑∞ m=0 ane λmtψm(x). It can be verified that L FP is a self-adjoint negative semi-definite operator so that · · · < λ2 < λ1 < λ0 = 0. Further- more, 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) = √ P0(x) + ∞∑ m=1 ane −|λm|tψm(x) 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 prob- lem LFPψ = λψ is important. It is worth noting that a Fokker-Planck equation can always be trans- formed to another Fokker-Planck equation with a constant diffusion coef- ficient with a suitable variable transformation [81]. If we transform x to a new variable y = ∫ x [B(x′)]−1/2 dx′ and define the new solutions to be p̃(y, t) = [B(x(y))]1/4p(x(y), t), Eq. (1.10) becomes ∂p̃ ∂t = 1 [ √ BP0]1/2 ∂ ∂y [√ BP0 ∂ ∂y ( p̃ [ √ BP0]1/2 )] = L̃FPp̃. (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 Ã(y) = P̃ ′0(y)/P̃0(y), analogous to Eq. (1.9). 1.3 The Schrödinger equation The time independent Schrödinger equation Hψ = − ~ 2 2µ d2ψ dy2 + V (y)ψ = Eψ, (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ödinger operator H is usually referred to as the Hamiltonian. The multi-dimensional Schrödinger equation Hψ = −∇2ψ + V (x)ψ is used in the studies of atoms or molecules. One of the applications for the Schrödinger 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 ex- act energy levels for the Morse potential are known, it serves as both a simple estimation of true energy levels and a benchmark for numerical tech- niques [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 L̃FPψ̃m = 1√ P̃0 d dy [ P̃0 d dy ( ψ̃m√ P̃0 )] = λ̃mψ̃m, (1.13) where λ̃m is the mth eigenvalues of L̃ FP, and ψ̃m(y) is the corresponding eigenfunction. Eq. (1.13) can be directly simplified to the form L̃FPψ̃m = d2ψ̃m dy2 − 1 2 P̃ ′′0 P̃0 − 1 4 ( P̃ ′0 P̃0 )2 ψ̃m = λ̃mψ̃m, (1.14) where ′ denotes d/dy. Equation (1.14) can be regarded as a Schrödinger equation. Therefore, the Schrödinger 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 Hamiltonian is exactly H = −L̃FP by comparing Eqs. (1.14) and (1.12). If we let λm be the eigenvalues of L FP, then λm = λ̃m = −Em. 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 ap- proximate solutions. For the solutions of eigenvalue problems such as the Schrödinger 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 prob- lem 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 eigen- value 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ödinger equation is rep- resented 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 ex- pansions. The differential equation is then represented (approximated) by a matrix problem and the eigenvalues are approximated by the eigenval- ues 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 di- agonlization 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 objec- tive 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 non- classical orthogonal polynomials. The implementation of all these methods are equivalent, and can be generalized as pseudospectral methods with or- thogonal 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 dissoci- ation 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 Chap- ter 4, we will apply transformation methods to both the polynomial and the sinc methods for the evaluation of eigenvalues in loosely bound states. Meth- ods 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))}Ni=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)}Nn=1 in the form f(x) = N∑ n=1 cnφn(x). The expansion coefficients cn are solved by setting f(xi) = g(xi). The func- tion 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)}Ni=1, and δij is the Kro- necker 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))}Ni=1 is represented by f(x) = N∑ j=1 g(xj)Cj(x). 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) = u n−1, the polynomial interpolating functions, called the Lagrange interpolating func- tion, are given by Ij(u) = N∏ i=1 i6=j u− ui uj − ui 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 f(u) = N∑ j=1 g(uj)Ij(u) (1.15) is exact, and the integral ∫ b a w(u)g(u) du = N∑ j=1 wjg(uj) (1.16) is exact for some function w(u) defined in an interval [a, b] if wj = ∫ b 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 set {uj}Nj=1, and the corresponding set of weights {wj}Nj=1, which satisfies∫ b a w(u)u k du = ∑N j=1wju k j for k = 0 . . . 2N−1. With uj and wj, the quadra- ture 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 re- currence relation Q0(u) = 1, Q1(u) = u− α0, (1.17a) Qn+1(u) = (u− αn)Qn(u)− βnQn−1(u), (1.17b) where the coefficients αn and βn are given by αn = 1 γn ∫ b a uw(u)[Qn(u)] 2 du, βn = γn γn−1 , (1.17c) γn = ∫ b a w(u)[Qn(u)] 2 du, (1.17d) satisfy the orthogonality ∫ b 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 de- terminant form Qn+1(u) = det u− α0 − √ β1 0 · · · 0 −√β1 u− α1 − √ β2 . . . ... 0 −√β2 u− α2 . . . 0 ... . . . . . . . . . −√βn 0 · · · 0 −√βn u− αn (1.18) 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∫ b a w(u)Pm−1(u)Pn−1(u) du = N∑ j=1 wjPm−1(uj)Pn−1(uj) = δmn for 1 ≤ m,n ≤ N . If Tjn = √wjPn−1(uj), then TtT = I = TTt, where the superscript t denotes the matrix transpose. Hence, T is an orthogonal ma- trix. By defining the Jacobi matrix Jmn = ∫ b a w(u)uPm−1(u)Pn−1(u) du = 11 Chapter 1. Introduction ∑N k=1wjujPm−1(uj)Pn−1(uj) for 1 ≤ m,n ≤ N and the diagonal matrix Uij = ujδij , the matrix J will satisfy J = T tUT, 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 J = α0 √ β1 0 · · · 0√ β1 α1 √ β2 . . . ... 0 √ β2 α2 . . . 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 = √ wjPn−1(uj). Since P0(uj) = γ −1/2 0 , the quadrature weights are given by wj = γ0([vj ]1) 2. The orthogonality ofT gives the equality ∑N n=1 √ wiwjPn−1(ui)Pn−1(uj) = δij . With √ wiwj replaced by wj, the polynomials Ij(u) = wj N∑ n=1 Pn−1(u)Pn−1(uj) 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 CPj (u) = √ w(u) w(uj) Ij(u), CPj (u) will be orthogonal with respect to the unit weight function, i.e.∫ b a C P i (u)C P j (u) du = ηjδij , where ηj = wj/w(uj). The constants ηj are the weights of the quadrature rule ∫ b a g(u) du ≈ N∑ j=1 ηjg(uj). 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] sinc(u) = sinπu πu if u 6= 0, 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 if −π ≤ ω ≤ π, 0 otherwise, which implies that ∫ ∞ −∞ sinc(u) exp(iωu) du = χ[−π,π](ω). (1.19) The orthogonality of translated sinc functions is derived by Parseval’s the- orem ∫ ∞ −∞ 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∫ ∞ −∞ g(u) du ≈ N∑ j=1 hg(uj) 13 Chapter 1. Introduction since ∫∞ −∞ Sj(u) du = h which can be derived by Eq. (1.19) for ω = 0. 1.5.3 Transformation of grid points The polynomial interpolation CPj (u) in Section 1.5.1 and the sinc interpo- lation 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 co- ordinates. 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 CPj (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 ηjδij = ∫ b a Ci(u)Cj(u) du = ∫ ρ−1(b) ρ−1(a) Ci(ρ(x))Cj(ρ(x))ρ ′(x) dx The new interpolating functions Ĉj(x) = √ ρ′(x) ρ′(xj) Cj(ρ(x)) will satisfy cardinality Ĉj(xi) = δij and orthogonality ∫ ρ−1(b) ρ−1(a) Ĉi(x)Ĉj(x) dx = η̂jδij , where η̂j = ηj/ρ ′(xj). The new quadrature rule becomes∫ ρ−1(b) ρ−1(a) g(x) dx ≈ N∑ j=1 η̂jg(xj). After mapping, the polynomial basis function becomes Ĉj(x) = √ ρ′(x)w(ρ(x)) ρ′(xj)w(ρ(xj)) Ij(ρ(x)) 14 Chapter 1. Introduction with quadrature weights η̂j = wj/[ρ ′(xj)w(ρ(xj))], and the sinc basis func- tion becomes Ĉj(x) = √ ρ′(x) ρ′(xj) Sj(ρ(x)) 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) = ∞∑ n=1 cnφn(x). Spectral convergence is achieved if the magnitudes of the expansion coef- ficients 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 condi- tions, 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∑ n=1 cn ∫ R φm(x)Lφn(x) dx = ∫ R φm(x)g(x) dx. (1.21) 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 inte- grals. The quadrature points are regarded to as the grid points for which the values of the solution are defined. The spectral method with quadra- ture evaluation of integral is referred to as the pseudospectral method. The main difference between spectral and pseudospectral methods is that spec- tral method is based on the expansion of basis functions, while pseudospec- tral method is based on a grid of points. 15 Chapter 1. Introduction The basis functions used in pseudospectral methods are usually interpo- lating 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 ψ(x) = N∑ j=1 √ ηjψ(xj) ( Cj(x)√ ηj ) . When the integral from the Galerkin method,∫ R Ci(x)√ ηi Lψ(x) dx = ∫ R Ci(x)√ ηi g(x) dx, 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∑ j=1 Lij √ ηjψ (m)(xj) = λm √ ηiψ (m)(xi), where Lij = (ηiηj) −1/2 ∫ RCi(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 ψ(x, t) = N∑ j=1 aj(t) ( 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 Lijaj(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 sin- gle approach. Although these pseudospectral methods have different back- grounds and theories, they are all equivalent to a method involving the ex- pansion 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, conver- gence 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ödinger eigen- value 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 trans- formations such as the hyperbolic sine mapping. A new flexible type of adap- tive 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. Numer- ical 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 Func- tional 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 mechan- ical 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ödinger 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. Com- put. 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 represen- tation 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 rep- resentation basis obtained by simultaneous diagonalization. J. Chem. Phys., 121:726, 2004. [28] L. Debnath. Nonlinear Partial Differential Equations for Scientists and Engineers. Birkhäuser, 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 nu- merical 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 non- local 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ón. Solutions through nonclassical potential symmetries for a generalized inhomogeneous nonlinear dif- fusion equation. Math. Meth. Appl. Sci., 31:753, 2008. [39] W. Gautschi. Orthogonal Polynomials: Computation and Approxima- tion. 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 mul- tiplicative 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 Equa- tions: 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änggi, 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. Mar- cel 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. Haewsungchar- ern, B. K. Bala, and J. Müller. 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ödinger 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ödinger 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 for- ward 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 perturba- tions 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 representa- tions 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ödinger equation. Chin. Phys. Lett., 21:1681, 2004. [65] X.-S. Liu and P.-Z. Ding. Numerical method based on Magnus ex- pansion and a new shooting method for eigenvalues of Schrödinger equation. Int. J. Quant. Chem., 103:149, 2005. [66] J. Lund and K. L. Bowers. Sinc Methods for Quadrature and Differ- ential Equations. SIAM, Philadelphia, 1992. [67] P. A. Markowich. Applied Partial Differential Equations. Springer- Verlag, Berlin, 2007. [68] L. A. Martino, A. Osella, C. Dorso, and J. L. Lanata. Fisher equation for anisotropic diffusion: simulating South American human disper- sals. 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ödinger equation by symplectic and trigono- metrically 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ánchez-Castellanos, Amezcua-Eccius, O. Álvarez 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ć 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éro-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 vibra- tional excitations in a heteronuclear diatomic molecule. J. Chem. Sci, 119:433, 2007. 25 Bibliography [89] B. Shizgal and R. Blackmore. Eigenvalues of the Boltzmann colli- sion 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ödinger 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ürmer, H. Köstler, and U. Rüde. 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şeli 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ázquez. 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 H2O. J. Chem. Phys., 106:6885, 1997. [104] E. T. Whittaker. On the functions which are represented by the expan- sions 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ödinger and Fokker-Planck equations: Comparison with sinc method 2.1 Introduction The solution of the Fokker-Planck and Schrödinger equation with a dis- cretization of the wave function on a grid of points has been of ongo- ing 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 discretiza- tion method (QDM), which is a pseudospectral (collocation) method based on nonclassical polynomials. For the solution of the Schrödinger equa- tion, Light and co-workers introduced the discrete variable representation (DVR) [27,38], which is often based on classical polynomials including Gaus- sian 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 con- volution (DSC) approach for solutions of the Fokker-Planck [68, 69] and Schrödinger 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 prob- lems 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ödinger 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 anal- ogous to techniques in neutron transport theory [24] and radiative trans- port [14] introduced over half a century ago. The evaluation of the matrix elements of the linear Boltzmann collision operator in some basis set (typ- ically 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 integrodif- ferential equation, Shizgal and Blackmore [8, 57] extended the formalism to the discretization of derivative operators. This technique was later applied to the Schrödinger 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 statis- tical mechanics [48]. 29 Chapter 2. Spectral convergence of quadrature discretization method The QDM differs from the other approaches in that it is based on non- standard polynomials orthogonal with respect to nonclassical weight func- tions. It originates from an interest in the solution of kinetic theory prob- lems [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 prob- lems is to consider the expansion of the distribution function in Laguerre (or Sonine) polynomials. For this application, with the small electron to mod- erator 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 con- vergence 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 antic- ipated. 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 eigenval- ues of the Fokker-Planck operator. The actual implementation of the QDM was later based on a grid defined by the quadrature associated with a non- classical 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ödinger equation is motivated by the relationship with the Fokker-Planck equation. The Schrödinger equations that are derived from a Fokker-Planck equation belong to the class of quantum problems in supersymmetric quantum me- chanics for which the ground state is known [18,29]. This chapter considers the solution of 1D Fokker-Planck and the Schrödinger 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 in- creased 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 multi- dimensional problems [36, 38, 45, 47]. The QDM has the potential to pro- vide 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 clas- sic 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 conver- gence 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 DSC- sinc 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ödinger equa- tion. 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 con- vergence 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ödinger 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ödinger equation is characterized by a triple well potential. This and similar bimodal systems are of consid- erable 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ödinger 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ödinger eigenvalue problems. The specific applica- tions 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 result- ing matrix representatives of the operators considered, and the fundamental reason for the rapid spectral convergence of the eigenvalues and eigenfunc- tions. 32 Chapter 2. Spectral convergence of quadrature discretization method 2.2 Discretization of the Fokker-Planck and Schrödinger 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 dx −B(x)d 2φm dx2 = λmφm, (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 Fokker- Planck 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}Ni=1 are the quadrature points defined by FN (xi) = 0, then the discrete representation of the operator L is [8, 63] LQDMij = − N∑ k=1 B(xk)[Dki + h(xk)δki][Dkj + h(xk)δkj], where the QDM representation of the derivative operator is given by [63] Dij = √ wiwj N∑ n=1 F ′ n−1(xi)Fn−1(xj), and h(x) = w′(x) 2w(x) − P ′ 0(x) 2P0(x) , where P0(x) = exp[− ∫ [A(x′)/B(x′)]dx′ − lnB(x)] is the equilibrium prob- ability 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 L̂QDMij = − N∑ k=1 B(xk)DkiDkj. (2.2) 33 Chapter 2. Spectral convergence of quadrature discretization method The eigenvalue problem, Eq. (2.1), can be transformed to a Schrödinger equation [63] with the transformation y = ∫ x [B(x′)]−1/2dx′ and we get Hψm(y) = −d 2ψm(y) dy2 + V (y)ψm(y) = λmψm(y), (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 HQDMij = N∑ k=1 DkiDkj + [V (yi)− V0(yi)]δij , where V0(y) = 1 2 w′′(y) w(y) − 1 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 ĤQDMij = N∑ k=1 DkiDkj. (2.5) 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 6= 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, D (1) ij = S ′ j,h(xi) = { (−1)i−j (i−j)h if i 6= j, 0 if i = j, (2.6) and D (2) ij = S ′′ j,h(xi) = { −2(−1)i−j (i−j)2h2 if i 6= j, − π2 3h2 if i = j. (2.7) which are the same expressions as given by Amore [1]. The explicit expres- sions, 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, LSCMij , is LSCMij = 1 h ∫ ∞ −∞ Si,h(x)LSj,h(x) dx and thus the discrete representation is given by LSCMij = −A(xi)D(1)ij +B(xi)D(2)ij , (2.8) which is not symmetric. However, the discrete representation of the Hamil- tonian is given by HSCMij = −D(2)ij + 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 = √ mv2/2kBTb, x ∈ [0,∞), where v is the electron speed, Tb is the temperature of the back- ground 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ödinger equation for this model, V (y) = y6/64− y2 + 15/4y2, is bounded and there are an infinite number of bound states. The weight function w(y) = y5 exp(−y4/16) can be derived from the su- perpotential [18,29] associated with the ground state eigenfunction for this potential [63]. Thus the eigenvalues can also be calculated from the numer- ical 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 dis- tribution of quadrature points for these weight functions is nonuniform. The calculation of the eigenvalues for this problem was considered by Wei [68] 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ödinger 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ödinger 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 λm = ~ 2α2 2µ ( m+ 1 2 )( 2 √ 2µD ~α − ( m+ 1 2 )) for m ≥ 0. The QDM uses the square of the ground state eigenfunction ψ0(y) = exp [ − √ 2µD ~ ( y + exp(−αy) α ) + αy 2 ] (2.10) as the weight function, i.e. w(y) = [ψ0(y)] 2. Thus the discretized Hamilto- nian with the QDM is HQDMij = ~ 2 2µ N∑ k=1 DkiDkj + λ0δij , where λ0 = V − ~22µ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 compu- tational 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, A(x) = x3 − x and B(x) = ǫ with x ∈ (−∞,∞), is P0(x) = exp[−x44ǫ + x 2 2ǫ ] [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 proba- bility 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ödinger equation is Hψm(y) = −ǫd 2ψm(y) dy2 + V (y)ψm(y) = λmψm(y) (2.11) with V (y) = (y3 − y)2 4ǫ − 3y 2 − 1 2 (2.12) 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 eigen- functions 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- cillator states, λ̂m = m and those for m even are triply degenerate except for 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(−y2/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 pro- vide 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 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 in- terval [0,∞). For the even eigenstates, we use web(y) = wa(y), y ∈ [0,∞), and for the odd eigenstates, wob (y) = y 2wa(y), y ∈ [0,∞). We occasion- ally 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 with this potential. The eigenvalue λ2 is a singlet state and tends to λ̂1 = 1 as ǫ→ 0. The states λ3, λ4 and λ5 make up a triplet which are degenerate for ǫ → 0 and tend to λ̂2 = 2. The remainder of the eigenvalue spectrum 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ödinger equations for several model systems as described in the pre- vious 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 ac- curate 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 eigenval- ues are calculated from the numerical diagonalization of Eq. (2.2). For the Schrödinger 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 eigenval- ues for the electron relaxation problem, the vibrational states for I2 modeled with a Morse potential and the quartic bistable system. The exact eigenval- 38 Chapter 2. Spectral convergence of quadrature discretization method ues, λexactm , used to evaluate the error are calculated to over 20 significant 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)m = log ∣∣∣∣∣λ (N) m − λexactm λexactm ∣∣∣∣∣ , (2.13) which is approximately equal to the number of significant figures in λ (N) m . For ease of notation, we henceforth write λ (N) n and ǫ (N) n as λn and ǫn, re- spectively. Although we compute the eigenvalues from a diagonalization of the discrete representations, LQDMij , Eq. (2.2), and H QDM ij , Eq. (2.5), it is 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 re- laxation 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ödinger equation (squares) is similar and spectral [9], and machine accuracy is achieved with approximately 20 basis functions, although accurate values with less preci- sion are obtained with fewer basis functions. We discuss the spectral conver- gence 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 else- where [2, 25, 37, 61]. For the Fokker-Planck operator, the weight function used is w(x) = x2 exp(−x2), whereas for the Schrödinger equation we use w(y) = y5 exp(−y4/16). The matrix representative of the Fokker-Planck op- erator 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ödinger equation. The convergence of the eigenvalues with the SCM (di- agonalization of Eq. (2.9)) is clearly considerably slower. The main reason 39 Chapter 2. Spectral convergence of quadrature discretization method 0 20 40 60 N -16 -12 -8 -4 0 QDM (se) QDM (fpe) SCM (se) RV A ǫ 1 0 20 40 60 N -16 -12 -8 -4 0 QDM (se) QDM (fpe) SCM (se) RV B ǫ 6 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ödinger equation (denoted by filled circles, xmax = 6). The QDM results are for the Schrödinger 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 fig- ures 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 conver- gence 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 re- sults 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 conver- gence 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 consider- ably 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 practi- cal 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 − λexactn for m = 1, 10, and 30 are at machine accuracy with N = 6, 24, and 70, respec- tively. This convergence is faster than the result reported by Mazziotti [42]. Figure 2.3 shows the variation of the potential in the Schrödinger, Eq. (2.12), and the distribution of the grid points for the quartic bistable prob- lems. 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[−(y4/4 − y2/2)/ǫ]/ exp(1/4ǫ) are almost uniformly distributed but with some localization in the vicinity 41 Chapter 2. Spectral convergence of quadrature discretization method 0 10 20 30 40 50 N -16 -12 -8 -4 0 SCM [-1.4, 1.4] QDM A ǫ 1 20 40 60 80 100 N -16 -12 -8 -4 0 SCM [-1.4, 2.4] QDM C ǫ 3 0 0 20 40 60 N -16 -12 -8 -4 0 SCM [-1.4, 1.4] QDM B ǫ 1 0 40 80 120 160 N -16 -12 -8 -4 0 SCM [-1.4, 3.9] QDM D ǫ 5 0 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 -2 -1 0 1 2 y -6 -4 -2 0 2 4 SCM QDM (a) V(y) A -1 0 1 y -4 0 4 SCM QDM (a) V(y) B QDM (b) -1 0 1 y -40 -20 0 20 40 QDM (a) QDM (c) SCM V(y) C Figure 2.3: Variation of the potential, V (y) = (y3 − y)2/4ǫ − (3y2 − 1)/2, in the Schrödinger 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[−(y4/4 − y2/2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞,∞), (b) web(y) = wa(y), y ∈ [0,∞) and (c) wc(y) = wa(y) + exp(−y2/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 web(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 su- persymmetric 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 quadra- ture points in this region, a Gaussian function is added so that the weight function becomes wc(y) = wa(y) + exp(−y2/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 simi- lar, 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 anal- ogous 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 eigen- values 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 λ4 λ5 ǫ1 ǫ2 ǫ3 ǫ4 ǫ5 QDMa 12 3.365796 0.939670 1.705240 2.743490 4.004819 -2.5 -1.9 -1.8 -1.2 -1.1 18 3.354699 0.927646 1.680910 2.601742 3.747467 -4.3 -3.5 -3.4 -2.6 -2.4 24 3.354531 0.927376 1.680273 2.595936 3.734290 -6.2 -5.4 -5.2 -4.4 -4.1 26 3.354529 0.927373 1.680266 2.595849 3.734062 -6.9 -6.0 -5.9 -5.0 -4.7 28 0.927372 1.680264 2.595827 3.734003 -7.6 -6.7 -6.5 -5.6 -5.3 30 0.927372 2.595822 3.733989 -8.3 -7.4 -7.2 -6.2 -5.9 32 2.595821 3.733986 -9.0 -8.1 -7.9 -6.9 -6.6 QDMb 3 3.367516 1.007699 1.764537 3.237921 4.994304 -2.4 -1.1 -1.3 -0.6 -0.5 6 3.355525 0.927402 1.684458 2.631757 3.807057 -3.5 -2.9 -2.6 -1.9 -1.7 9 3.354543 0.927372 1.680328 2.596849 3.735433 -5.4 - 4.5 -4.4 -3.4 -3.4 12 3.354529 1.680264 2.595835 3.733999 -7.5 -6.4 -6.5 -5.3 -5.4 15 2.595820 3.733985 -9.8 -8.5 -8.7 -7.3 -7.6 SCMc 12 6.825052 0.957566 1.817707 2.974622 4.093426 0.0 -1.5 -1.1 -0.8 -1.0 18 3.427580 0.926484 1.679049 2.597213 3.750052 -1.7 -3.0 -3.1 -3.3 -2.4 24 3.355037 0.927361 1.680236 2.595783 3.734026 -3.8 -4.9 -4.8 -4.8 -5.0 28 3.354524 0.927372 1.680265 2.595824 3.733991 -5.9 -6.4 -6.1 -5.9 -5.8 30 3.354531 1.680264 2.595821 3.733986 -6.1 -7.6 -7.7 -7.6 -6.7 32 3.354529 2.595820 3.733985 -7.8 -8.1 -7.8 -7.5 -7.2 a wa(y) = exp[−(y4/4− y2/2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞,∞). b For the even eigenstates, web(y) = wa(y), y ∈ [0,∞). For the odd eigenstates, wob (y) = y2wa(y), y ∈ [0,∞). c 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 QDMa 12 5.0833(-8) 1.866176 1.865861 3.9 -3.1 -4.2 24 3.6651(-11) 1.865757 1.865753 0.7 -3.2 -6.9 36 7.0354(-12) 1.388230 1.865752 1.865758 2.664871 -0.8 -0.4 -3.2 -5.6 -0.4 48 6.1809(-12) 0.994289 1.865735 1.865754 1.956370 -2.4 -1.6 -3.2 -6.6 -1.3 60 6.15499(-12) 0.968472 1.865337 1.869329 -2.8 -3.2 -3.4 -7.9 -2.9 72 6.15466(-12) 0.967870 1.864560 1.866993 -4.5 -5.3 -5.0 -9.8 -5.0 84 6.15465(-12) 0.967865 1.864542 1.866975 -3.4 -7.6 -7.3 -12.1 -7.3 90 0.967864 -2.5 -8.9 -8.5 -13.0 -8.5 QDMb 12 6.4259(-12) 1.256087 1.865747 1.865757 2.113341 -1.4 -0.5 -3.2 -5.7 -0.9 15 6.1656(-12) 0.990778 1.865720 1.865754 1.913825 -3.0 -1.6 -3.2 -6.8 -1.6 18 6.1405(-12) 0.969092 1.865601 1.875631 -2.5 -2.9 -3.2 -7.5 -2.3 24 6.1424(-12) 0.967879 1.864549 1.866982 -2.6 -4.8 -5.4 -9.6 -5.4 27 6.1436(-12) 0.967865 1.864542 1.866975 -2.6 -6.5 -6.6 -11.2 -6.6 30 6.1427(-12) 0.967864 -2.6 -7.4 -7.7 -12.6 -7.7 SCMc 12 7.4085(-1) 1.076821 3.336192 3.321671 1.574892 11.1 -0.9 -0.1 -0.1 -0.8 24 3.3865(-3) 0.967915 1.931199 1.930972 1.865051 8.7 -4.3 -1.4 -1.5 -3.0 36 -4.8093(-5) 0.967864 1.864066 1.864927 1.866629 6.9 -9.6 -3.6 -3.4 -3.7 48 1.7533(-8) 1.864542 1.865754 1.866975 3.5 -13.7 -6.8 -6.5 -6.8 60 7.9956(-12) -0.5 -13.6 -10.8 -10.5 -10.8 66 6.1435(-12) -2.6 -13.6 -12.6 -12.2 -12.5 a wa(y) = exp[−(y4/4 − y2/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 web(y) = wa(y), y ∈ [0,∞). For the odd eigenstates, wob (y) = y2wa(y), y ∈ [0,∞). c xmax = 1.5. 46 Chapter 2. Spectral convergence of quadrature discretization method Table 2.3: Convergence of the eigenvalues for the quartic potential, ǫ = 0.001. N λ2 λ3 λ4 λ5 ǫ2 ǫ3 ǫ4 ǫ5 QDMa 6 0.9980526 2.0000470 2.0200067 2.0694590 -3.0 -2.2 -1.8 -1.4 12 0.9969809 1.9878205 1.9880010 1.9881554 -6.1 -4.5 -4.3 -3.9 18 0.9969817 1.9878873 1.9878903 1.9878937 -7.1 -5.9 -6.5 -5.7 24 1.9878896 1.9878896 1.9878893 -8.3 -8.2 -8.5 -6.8 30 1.9878896 -8.7 -9.8 -8.9 -7.7 SCMb 12 3.4140030 21.7611343 21.7517512 3.4979300 0.4 1.0 1.0 -0.1 24 1.2875835 8.5716393 8.5704015 1.6476678 -0.5 0.5 0.5 -0.8 36 1.0279928 4.1568953 4.1569000 1.7776398 -1.5 0.0 0.0 -1.0 48 0.9984717 2.5344712 2.5345542 1.9649680 -2.8 -0.6 -0.6 -1.9 60 0.9970079 2.0928734 2.0928574 1.9872886 -4.6 -1.3 -1.3 -3.5 72 0.9969819 1.9995592 1.9995565 1.9878842 -6.8 -2.2 -2.2 -5.6 84 0.9969817 1.9884160 1.9884156 1.9878896 -9.4 -3.6 -3.6 -8.0 96 1.9878838 1.9878838 -12.3 -5.5 -5.5 -10.9 108 1.9878884 1.9878884 -13.2 -6.2 -6.2 -13.8 120 1.9878896 1.9878896 -13.9 -7.8 -7.8 -13.5 a wc(y) = exp[−(y4/4− y2/2)/ǫ]/ exp(1/4ǫ) + exp(−y2/2ǫ), y ∈ (−∞,∞). b xmax = 1.2. 47 Chapter 2. Spectral convergence of quadrature discretization method 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (a) QDM (b) SCM A ǫ 2 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (a) QDM (b) SCM C ǫ 4 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (a) QDM (b) SCM B ǫ 3 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (a) QDM (b) SCM D ǫ 5 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[−(y4/4 − y2/2)/ǫ]/ exp(1/4ǫ), y ∈ (−∞,∞). For the QDM(b), web(y) = wa(y), y ∈ [0,∞) for the even eigenstates and wob (y) = y 2wa(y), y ∈ [0,∞) for the odd eigenstates. For the SCM, xmax = 1.5. 48 Chapter 2. Spectral convergence of quadrature discretization method 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (c) QDM (b) SCM A ǫ 2 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (c) QDM (b) SCM QDM (a) C ǫ 4 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (c) QDM (b) SCM QDM (a) B ǫ 3 0 40 80 120 N -16 -12 -8 -4 0 4 QDM (c) QDM (b) SCM D ǫ 5 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 wob (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 func- tions 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ödinger 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 wob (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) pro- vides 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 wob (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 eigen- values 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 con- verge 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 con- vergence of the solutions of the Fokker-Planck and Shrödinger 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, |λ(N)m − λexactm | ≈ 10−qN . (2.14) Spectral convergence has been demonstrated in the results for which log |λ(N)m − λexactm | varies linearly with N . The convergence is from above the exact eigenvalue as shown in the tables consistent with a Rayleigh-Ritz varia- tional 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 eval- uated 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 effi- cient 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 mathemat- ical properties. We focus our attention on the electron relaxation problem from both the Fokker-Planck and Schrödinger 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 tridiag- onal matrix. Appendix B provides the details of the calculation of the matrix represen- tative 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 Fokker- Planck 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 eigenfunc- tions 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 web(y) = wa(y), y ∈ [0,∞), and wob (y) = y 2wa(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 eigen- values λ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 dif- fusion 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 off- diagonal 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−1Pn−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 10 20 30 n -16 -12 -8 -4 0 lo g |d n (m ) | m = 1 A m = 3 m = 5 0 10 20 30 n -16 -12 -8 -4 0 lo g |d n (m ) | m = 1 B m = 3 m = 5 Figure 2.6: Variation of the expansion coefficients d (m) n versus n of the mth eigenfunction of the Fokker-Planck operator (A) and the Hamiltonian in the Schrödinger 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 λ1 = a1 + c21 λ1 − a2 − c 2 2 λ1 − a3 − c 2 3 λ1 − a4 − ··· . (2.15) The other eigenvalues do not have a similar continued fraction representa- tions except for λN . The Nth approximation to λ1 is obtained by truncating the continued fraction with the term in c2N+1. The numerical value of λ (N) 1 can be determined by starting with the initial estimate λ (1) 1 = a1 and cal- culating successive approximations by substituting for λ1 on the right hand side of Eq. (2.15). The rate of convergence of this iterative process to λ (N) 1 is analogous to the convergence of λ (N) 1 to λ1 versus N . With Eq. (2.15), we obtain the approximation λ (2) 1 − λ(3)1 ≈ c21c 2 2 (λ1 − a3)(λ1 − a2)2 , where on the right hand side of the continued fraction we have set λ (2) 1 ≈ λ (3) 1 = λ1. Similarly we can calculate λ (3) 1 − λ(4)1 ≈ c21c 2 2c 2 3 (λ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) 1 − λ(4)1 λ (2) 1 − λ(3)1 = c23 (λ1 − a4)(λ1 − a3) = 10 −q. (2.16) With Eq. (2.16) and the value of λexact1 = 4.68340, we find that q = 1.05 is 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→∞ c 2 N−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 ma- trices 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 variation of the nth elements in the mth eigenvectors, d (m) n , 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 fig- ure is a demonstration of spectral convergence. The expansion coefficients for the eigenvectors, d (m) n , 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 cn−1d (1) n−1 + and (1) n + cnd (1) n+1 = λ exact 1 d (1) n . If the coefficient d (1) n+1 is neglected, we have the approximation d (1) n d (1) n−1 = cn−1 λexact1 − an . The slope of the graph in Fig. 2.6 is thus log(cn−1/|λexact1 −an|), the asymp- 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 ma- trix elements of the full matrix representative of the Schrödinger equa- tion for the same electron relaxation problem with the weight function w(y) = y5 exp(−y4/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ödinger 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] rela- tive to the phase space extent of the classical Hamiltonian corresponding to the quantum Hamiltonian operator. These analyses have been consid- ered 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 es- sential 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 spec- tral 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ödinger equation can be extremely rapid. For the potential in the Schrödinger equa- tion 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 equilib- rium 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 relax- ation 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ández, and H. Jones. A new approx- imation 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 infe- rior 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ödinger 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 representa- tion 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 semi- classical 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 repre- sentation. 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 optimiz- ing grid representations: the mapped Fourier method. Phys. Rev. E, 53:1217, 1996. [24] R. D. M. Garcia. The application of nonclassical orthogonal polynomi- als 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änggi, 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ácek, and L. S. Cederbaum. Schwinger and anomaly-free Kohn variational principles and a generalized Lanczos al- gorithm 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 Carte- sian coordinates and a rotational symmetry adapted Lanczos method. J. Chem. Phys., 119:6609, 2003. [46] R. Peyret. Spectral Methods for Incompressible Viscous Flow. Springer- Verlag, New York, 2002. [47] B. Poirier and J. C. Light. Phase space optimization of quantum repre- sentations: direct-product basis sets. J. Chem. Phys., 111:4869, 1999. [48] H. Risken. The Fokker-Planck Equation: Methods of Solution and Ap- plication. 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 func- tions. 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 reac- tions. 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 Boltz- mann collision operator for systems of two components at different tem- peratures. Chem. Phys., 6:54, 1974. [60] B. Shizgal, J. M. Lindenfeld, and R. Reeves. Eigenvalues of the Boltz- mann collision operator for binary gases: Mass dependence. Chem. Phys., 56:249, 1981. [61] B. Shizgal and D. R. A. McMahon. Electric field dependence of tran- sient 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ödinger 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ó, A. Nagy, T. Furtenbacher, and A. Császár. On one- dimensional discrete variable representations with general basis func- tion. 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 Fokker- Planck 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 H2O. 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 compli- cated 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 pseu- dospectral 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ödinger equation 3.1 Introduction There has been an ongoing effort by numerous researchers to develop ac- curate and efficient algorithms for the calculation of the eigenvalues of the 1D Schrödinger 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 Nu- merov [20], Runge-Kutta [1,21] or symplectic methods [30,31]. Other meth- ods 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şeli 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 N2 gives N eigenvalues of which a subset corresponds to the discrete eigenvalues for the problem. The spectral convergence of the eigen- values 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/s10910- 007-9341-8 65 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation eigenvalues versus N . Spectral convergence has been demonstrated recently by Lo and Shizgal [28] for several 1D problems. The present chapter is a con- tinuation 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 conver- gence of the eigenvalues of the Fokker-Planck equation can be obtained with a basis set defined by the equilibrium probability density as the weight func- tion. The Fokker-Planck equation can be transformed to a Schrödinger equa- tion with a potential such that the Hamiltonian belongs to the class of su- persymmetric 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 conver- gence [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 pseudospec- tral 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 previ- ous 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ödinger 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ödinger equation [12, 19]. If V (x) is the potential in a 1D Schrödinger equation, then the matrix elements of V (x) can be determined with a Gaus- 66 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation sian quadrature, that is, Vmn = ∫ b a φm(x)V (x)φn(x) dx ≈ N∑ k=1 ηkφm(xk)V (xk)φn(xk), (3.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), ∫ R w(x)Pm(x)Pn(x) dx = δmn, where R denotes the domain for x. With wk = ηkw(xk), there is a unitary transformation, Tkn = √ wkPn−1(xk), between the representation of a func- tion 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 trans- formation 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 princi- ple 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ödinger equation Hψ = − ~ 2 2µ d2ψ dx2 + V ψ = Eψ 67 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation 0 10 20 30 40 n -12 -8 -4 0 lo g|c n (m ) | A ψ1, s = 1 ψ10, s = 1 0 10 20 30 40 n -12 -8 -4 0 lo g|c n (m ) | B ψ1, s = 0.125 ψ10, s = 0.14 Figure 3.1: Variation of the expansion coefficients c (m) n 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 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation is generally considered with the expansion of the eigenfunctions in a basis set, ψ(x) ≈ √ w(x) N∑ n=1 cnPn−1(x). (3.2) The matrix representation of the Hamiltonian in this basis set is Hmn = (~2/2µ)Kmn + Vmn, where Kmn = − ∫ R √ w(x)Pm−1(x) d2 dx2 [√ w(x)Pn−1(x) ] dx, (3.3a) Vmn = ∫ R w(x)Pm−1(x)V (x)Pn−1(x) dx. (3.3b) The expansion coefficients cn in Eq. (3.2) can be considered as linear vari- ational 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 repre- sentation (VBR) as proposed by Light et al [25]. The eigenvalue estimates converge monotonically to the exact eigenvalues from above. Spectral con- vergence 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 poly- nomials orthogonal with respect to a chosen weight function, w(x) can be generated with the algorithm as described elsewhere [17, 18]. These poly- nomials define the quadrature points xk and weights wk for the quadrature rule ∫ R w(x)f(x) dx ≈ N∑ k=1 wkf(xk). (3.4) 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 ψ(x) ≈ N∑ k=1 Ik(x) √ w(x) w(xk) ψ(xk), 69 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation where the interpolating polynomial, Ik(x), is given by Ik(x) = wk N∑ n=1 Pn−1(xk)Pn−1(x) 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 evalu- ated and one of the cross terms is integrated by parts, Kmn can be written as Kmn = ∫ R w(x)P ′m−1(x)P ′ n−1(x) dx − ∫ R w(x)Pm−1(x)Ṽ (x)Pn−1(x) dx, where Ṽ (x) = 1 2 w′′(x) w(x) − 1 4 [ w′(x) w(x) ]2 is a reference potential [39]. The Hamiltonian matrix is therefore Hmn = ~ 2 2µ ∫ R w(x)P ′m−1(x)P ′ n−1(x) dx + ∫ R w(x)Pm−1(x) [ V (x)− ~ 2 2µ Ṽ (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 HQDMij = ~ 2 2µ N∑ k=1 DkiDkj + [ V (xi)− ~ 2 2µ Ṽ (xi) ] δij , (3.6) where Dij = √ wiwj N∑ n=1 P ′n−1(xi)Pn−1(xj) 70 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation is the discrete representation of the first derivative operator. Szalay [43] has provided explicit formulas forDij = ∫ R w(x)Pi(x)P ′j(x) dx and ∑N k=1DkiDkj =∫ R w(x)P ′i(x)P ′j(x) dx where Pi(x) = w −1/2 i Ii(x) (Table 1 of Ref. 43). For the QDM, the reference potential, Ṽ (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µ)Ṽ (x), the discrete representation of the Hamiltonian reduces to HQDMij = (~ 2/2µ) ∑N k=1DkiDkj , and is calculated exactly with the quadrature. Thus, the QDM preserves the variational as- pects 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ödinger equations that are isospectral with equivalent Fokker-Planck equations [33,39]. The eigenfunc- tion of the Fokker-Planck equation with zero eigenvalue is the equilibrium distribution function which is known, Peq(x) = exp (− ∫ xW (x′) dx′). For these problems, the ground-state of the Schrödinger equation is ψ0(x) =√ Peq(x). The function W is the superpotential in supersymmetric quan- tum mechanics [13, 39]. If there is no convenient choice for w(x) such that V (x) = (~2/2µ)Ṽ (x), one can choose a w(x) so that (~2/2µ)Ṽ (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 HQDMij = ~ 2 2µ 1 s2 N∑ k=1 DkiDkj + [ V (sxi + b)− ~ 2 2µ 1 s2 Ṽ (xi) ] δij , (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ödinger 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 Em = ~ 2 2µ α2 ( m+ 1 2 )( 2 α √ 2µ ~2 De − ( m+ 1 2 )) −D (3.9) for 0 ≤ m ≤ 77. The ground state eigenfunction is ψMorse0 (x) = exp [ − √ 2µ ~2 De ( x− xe + exp(−α(x− xe)) α ) + α(x− xe) 2 ] . (3.10) The Morse potential belongs to the family of shape invariant potentials of supersymmetric quantum mechanics [13]. If the weight function for the poly- nomial expansion is chosen to be w(x) = [ψMorse0 (x)] 2, the effective potential in Eq. (3.6) is simply (~2/2µ)Ṽ (xi) = VMorse(xi) − E0. With this choice of weight function, the discrete representation of the Hamiltonian, Eq. (3.6), reduces to HQDMij = (~ 2/2µ) ∑N k=1DkiDkj + E0δij . With scaling s and translation b, the QDM representation of the Hamil- tonian 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, La- guerre, 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 classi- cal polynomial basis such as Hermite polynomials. The DVR representation of the discrete Hamiltonian is HHermij = ~ 2 2µ 1 s2 KHermij + VMorse(sxi + b)δij (3.11a) with the kinetic energy operator KHermij = − { δij [ −2(N − 1) 3 − 1 2 + x2i 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ödinger equation but are 72 Chapter 3. Pseudospectral methods of solution of the Schrödinger 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)− d x6 + ex−6 (3.12) where x ∈ [0,∞) is the radial coordinate in Å and V (x) in eV. The pa- rameters in Eq. (3.12) are a = 1720 eV, b = 2.6920 Å−1, c = 0.2631 Å−1, d = 37.943 eV·Å6, e = 177588 Å12, ~2 = 0.0041801588 eV·u·Å2, 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, α = 1.685967091 Å−1, and xe = 3.761961562 Å. The weight function is thus w(x) = [ψMorse0 (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µ)Ṽ (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 Hamilto- nian in the Laguerre basis is given by HLagij = ~ 2 2µ 1 s2 KLagij + V (sxi + b)δij , (3.13a) where the discrete symmetric matrix representative of the kinetic energy is given by [2] KLagij = { 9/4x2i + Sii if i = j, (−1)i−j [(3/2)(xixj)−1/2(x−1i + x−1j ) + Sij] if i 6= j, (3.13b) and Sij = √ xixj N∑ k=1,k 6=i,j 1 xk(xk − xi)(xk − xj) . (3.13c) 73 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation 3.6 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 1 + t − u0t 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 if 0 ≤ x < L, 0 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) = [ψ sq 0 (x)] 2, where ψsq0 (x) = { c1 sin(c2x) if 0 ≤ x < L, exp(c3x) if x ≥ L, (3.14) c1 = exp(c3L)/ sin(c2L), 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 − u0L) = − √−ǫ0, exceeding u0, or ǫ0 ≈ −49.75457960982555. The QDM representation is given by Eq. (3.7), where (~2/2µ)Ṽ (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 es- timate obtained with the QDM with no translation or scaling of the quadra- ture 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ödinger equation Table 3.1: Convergence of 103E10 for the I2 Morse potential. N QDMa QDMb Hermitec 11 6.397074 6.432646 9.203165 12 5.929980 5.749709 6.515301 13 5.731211 5.634366 5.285059 14 5.653056 5.625372 5.229673 15 5.629265 5.623171 5.177875 16 5.624138 5.623269 5.524050 17 5.623357 5.623259 5.464870 18 5.623268 5.623260 5.615021 19 5.623260 5.619430 20 5.623601 21 5.622067 22 5.623292 23 5.623252 24 5.623235 25 5.623258 26 5.623257 27 5.623260 28 5.623259 29 5.623260 a s = 1, b = 0. b s = 1.05, b = 0. c s = 0.14, b = 0. 75 Chapter 3. Pseudospectral methods of solution of the Schrödinger 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 min- imum versus N at N = 15. The scaling s and translation b changes the weight function to [ψMorse0 ((x − b)/s)]2 so that V (x) − (~2/2µ)Ṽ (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 Hamilto- nian 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 po- tential given by Eq. (3.8). The negative of the relative error, defined by ǫm(N) = log ∣∣∣∣∣E (N) m − Eexactm Eexactm ∣∣∣∣∣ , represents approximately the number of significant figures for E (N) m . The exact eigenvalues Eexactm are given by Eq. (3.9). The QDM results for E1 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 conver- gence 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 impor- tant 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- tion, and thus we obtain the very rapid convergence of E (N) 1 , as shown in Fig. 3.2A. 76 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation 0 10 20 30 N -16 -12 -8 -4 0 lo g|( E 1 (N ) -E 1e x ac t )/E 1e x ac t | QDM, s = 1 Hermite, s = 0.125, b = 0 A 0 20 40 60 80 100 N -16 -12 -8 -4 0 lo g|( E 3 0(N ) -E 30 ex ac t )/E 30 ex ac t | QDM, s = 1 Hermite, s = 0.165, b = 0.1 QDM, s = 1.15 C 0 10 20 30 40 50 N -16 -12 -8 -4 0 lo g|( E 1 0(N ) -E 10 ex ac t )/E 10 ex ac t | QDM, s = 1 Hermite, s = 0.14, b = 0 QDM, s = 1.05 B 0 40 80 120 160 N -16 -12 -8 -4 0 lo g|( E 5 0(N ) -E 50 ex ac t )/E 50 ex ac t | QDM, s = 1 Hermite, s = 0.2, b = 0.2 QDM, s = 1.3 D Figure 3.2: Variation of log |(E(N)m − Eexactm )/Eexactm | versus the number 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ödinger equation In Fig. 3.1A we show the variation of the coefficients c (m) n for the mth eigenfunction versus n for m = 1 (triangles) and m = 10 (squares). The rapid exponential decrease of c (1) n versus n is an illustration of spectral con- vergence. 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. The variation of the coefficients c (m) n versus n for m = 1 and 10 in the Hermite expansion ψm(y) ≈ exp ( −1 2 [ y − b s ]2) N∑ n=1 c(m)n Q (s,b) n−1(y), (3.15) where y = sx + b is the scaled and translated coordinate, is shown in Fig. 3.1B with triangles and squares, respectively. In Eq. (3.15), Q (s,b) n (y) = νnHn((y− b)/s) are the scaled and translated Hermite polynomials normal- ized by the constants νn. Spectral convergence is also confirmed, but the rate of decrease of c (m) n is slower than with the QDM. For m = 10, the spectral convergence refers to the exponential decrease in the expansion co- efficients c (10) n 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, Eexactm , used for calculating the errors are evaluated in multiple precision 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 scal- ings 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 La- guerre 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 compu- tational 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 eigen- value is nearly optimized. 78 Chapter 3. Pseudospectral methods of solution of the Schrödinger equation 0 20 40 60 80 100 N -16 -12 -8 -4 0 lo g|( E 0 (N ) -E 0e x ac t )/E 0e x ac t | QDM, s = 1 Laguerre, s = 0.02, b = 0 A Laguerre, s = 0.03, b = 2.5 0 40 80 120 N -16 -12 -8 -4 0 lo g|( E 5 (N ) -E 5e x ac t )/E 5e x ac t | QDM, s = 1.3 Laguerre, s = 0.025, b = 0 C Laguerre, s = 0.05, b = 2.5 0 20 40 60 80 100 N -16 -12 -8 -4 0 lo g|( E 3 (N ) -E 3e x ac t )/E 3e x ac t | QDM, s = 1.05 Laguerre, s = 0.02, b = 0 B Laguerre, s = 0.035, b = 2.5 0 40 80 120 160 N -16 -12 -8 -4 0 lo g|( E 7 (N ) -E 7e x ac t )/E 7e x ac t | QDM, s = 2.0 Laguerre, s = 0.055, b = 0 D Laguerre, s = 0.1, b = 2.5 Figure 3.3: Variation of log |(E(N)m − Eexactm )/Eexactm | versus the number of 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ödinger 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 sub- traction of (~2/2µ)Ṽ (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 develop- ing 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 equa- tion [5,36]. The numerical implementation of the QDM and the DVR is with a pseudospectral approach that was developed by researchers in fluid me- chanics [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ödinger equation 0 20 40 60 N -16 -12 -8 -4 0 lo g|( E 0 (N ) -E 0e x ac t )/E 0e x ac t | A s = 1 s = 1.2 0 40 80 120 N -16 -12 -8 -4 0 lo g|( E 9 (N ) -E 9e x ac t )/E 9e x ac t | C s = 1.25 s = 1.5 0 20 40 60 80 N -16 -12 -8 -4 0 lo g|( E 4 (N ) -E 4e x ac t )/E 4e x ac t | B s = 1.1 s = 1.3 0 40 80 120 160 200 N -12 -8 -4 0 lo g|( E 1 3(N ) -E 13 ex ac t )/E 13 ex ac t | D s = 1.6 s = 1.9 Figure 3.4: Variation of log |(E(N)m − Eexactm )/Eexactm | versus the number of quadrature points N for the Woods-Saxon potential by QDM. (A): m = 0; (B): m = 4; (C): m = 9; (D): m = 13. 81 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 nu- merical solution of the Schrödinger 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ödinger 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ödinger 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 invari- ance, and exactly solvable potentials. Am. J. Phys, 56:163, 1988. [14] J. Echave and D. C. Clary. Potential optimized discrete variable repre- sentation. Chem. Phys. Lett., 190:225, 1992. [15] S. T. Epstein. The Variational Method in Quantum Chemistry. Aca- demic Press, New York, 1974. [16] F. M. Fernández, Q. Ma, and R. H. Tipping. Eigenvalues of the Schrödinger equation via the Riccati-Padé method. Phys. Rev. A, 40:6149, 1989. [17] R. D. M. Garcia. The application of nonclassical orthogonal polynomi- als in particle transport theory. Prog. Nucl. Energy, 35:249, 1999. [18] W. Gautschi. Orthogonal Polynomials: Computation and Approxima- tion. 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ödinger 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öm methods for the numerical solution of the Schrödinger 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 repre- sentations. 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ödinger equation by a high order perturbation method based on a linear reference potential. Comput. Phys. Com- mun., 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 representa- tions 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 dis- cretization method in the solution of the Schrödinger 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ödinger equation. J. Math. Chem., 40:257, 2006. [31] T. Monovasilis and T. E. Simos. Symplectic methods for the numerical integration of the Schrödinger equation. Comput. Mater. Sci., 38:526, 2007. [32] R. Peyret. Spectral Methods for Incompressible Viscous Flow. Springer- Verlag, New York, 2002. [33] H. Risken. The Fokker-Planck Equation: Methods of Solution and Ap- plication. 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ödinger 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ödinger equation. Appl. Math. Comput., 106:245, 1999. [42] T. E. Simos. A four-step exponentially fitted method for the numerical solution of the Schrödinger 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şeli 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şeli and M. Bahar Erseçen. 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ödinger 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ödinger equation. Comput. Phys. Commun., 160:23, 2004. [48] H. Wei. Ghost levels and near variational forms of the discrete variable representation: application to H2O. 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ödinger equation with arbitrary potential. Comput. Phys. Com- mun., 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 tem- peratures [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ödinger equation. Shizgal and co-workers [7,20,26–29] de- veloped 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 dis- crete variable representation (DVR) which is often employed with quadra- ture 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 poten- tial matrix elements need not be evaluated. Rather, the potential is eval- uated at the quadrature points but the accuracy of the method depends on whether the collocation provides accurate estimates of the matrix ele- ments [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 poten- tial 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 numer- ical 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 poten- tial 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 spec- trum. A mapped Fourier method has been extensively used for the determina- tion 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 Hamilto- nian, referred to as “ghost” levels, can occur [14,34,36]. In this chapter, we propose analogous mapping procedures for polyno- mial 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 ap- ply 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 Solution of the Schrödinger equation on a grid of quadrature points 4.2.1 The coefficient space representation The determination of the eigenvalues, Em, of the Schrödinger equation Hψm(x) = − ~ 2 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 ψm(x) ≈ ψ(N)m (x) = N∑ n=1 ĉ(m)n φ̂n(x). (4.2) The choice of basis set {φ̂n(x)}∞n=1 is derived from the orthogonality of a polynomial set {Pn(u)}∞n=0 given by ∫ b 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 ∫ b a w(u)Pm(u)Pn(u) du = ∫ ρ−1(b) ρ−1(a) ρ′(x)w(ρ(x))Pm(ρ(x))Pn(ρ(x)) dx = δ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) = √ ŵ(x)Pn−1(ρ(x)) with the weight function ŵ(x) = ρ ′(x)w(ρ(x)). A quadrature rule for some arbitrary function g(x) with this weight function is ∫ ρ−1(b) ρ−1(a) ŵ(x)g(x) dx = ∫ b a w(u)g(ρ−1(u)) du ≈ N∑ k=1 wkg(ρ −1(uk)) = N∑ k=1 wkg(xk). (4.4) The quadrature points {uk} and weights {wk} are associated with the poly- nomial PN (u) [12], and the quadrature points for x are xk = ρ −1(uk). Equa- tion (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 ŵ(x) in the form ∫ ρ−1(b) ρ−1(a) g(x) dx ≈ N∑ k=1 η̂kg(xk), (4.5) where the new quadrature weights are η̂k = wk/ŵ(xk). The matrix representation of the Hamiltonian in the new basis is ĤCmn = (~2/2µ)K̂Cmn + V̂ C mn, where K̂Cmn = − ∫ ρ−1(b) ρ−1(a) φ̂m(x)φ̂ ′′ n(x) dx, (4.6a) V̂ Cmn = ∫ ρ−1(b) ρ−1(a) φ̂m(x)V (x)φ̂n(x) dx, (4.6b) 90 Chapter 4. An efficient mapped pseudospectral method for weakly bound states and C denotes the coefficient space representation. Symmetrization of K̂Cmn is done with an integration by parts, i.e. K̂Cmn = ∫ ρ−1(b) ρ−1(a) φ̂′m(x)φ̂ ′ n(x) dx (4.7) with the assumption that the boundary term vanishes. With the explicit definition of φ̂n(x), K̂ C mn becomes K̂Cmn = ∫ ρ−1(b) ρ−1(a) [ρ′(x)]3w(ρ(x))P ′m−1(ρ(x))P ′ n−1(ρ(x)) dx − ∫ ρ−1(b) ρ−1(a) ρ′(x)w(ρ(x))Pm−1(ρ(x))Ṽ (x)Pn−1(ρ(x))) dx, (4.8a) where P ′(ρ(x)) = (dP (u)/du)|u=ρ(x), and Ṽ (x) = [ 1 2 w′′(ρ(x)) w(ρ(x)) − 1 4 ( w′(ρ(x)) w(ρ(x)) )2] [ρ′(x)]2 + 1 2 ρ′′′(x) ρ′(x) − 1 4 ( ρ′′(x) ρ′(x) )2 + w′(ρ(x))ρ′′(x) w(ρ(x)) . (4.8b) 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, K̂Cmn = N∑ k=1 wk[ρ ′(xk)] 2P ′m−1(uk)P ′ n−1(uk)− N∑ k=1 wkPm−1(uk)Ṽ (xk)Pn−1(uk). The potential energy matrix V̂ Cmn is evaluated similarly. Thus, the Hamilto- nian is given by ĤCmn = ~ 2 2µ N∑ k=1 wk[ρ ′(xk)] 2P ′m−1(uk)P ′ n−1(uk) + N∑ k=1 wkPm−1(uk) [ V (xk)− ~ 2 2µ Ṽ (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 The coefficients ĉ (m) n in Eq. (4.2) are given by ĉ(m)n = ∫ ρ−1(b) ρ−1(a) φ̂n(x)ψm(x) dx ≈ N∑ k=1 η̂kφ̂n(xk)ψm(xk), where the quadrature rule Eq. (4.5) has been used. Hence, Eq. (4.2) becomes ψm(x) ≈ ψ(N)m (x) = N∑ k=1 ψm(xk)Ĉk(x), (4.10) where the interpolating function Ĉk(x) is given by Ĉk(x) = η̂k N∑ n=1 φ̂n(xk)φ̂n(x) and satisfies the orthogonality ∫ ρ−1(b) ρ−1(a) Ĉk(x)Ĉl(x) dx = √ η̂kη̂lδkl and the cardinality Ĉ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 Ĉ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- tion ψ (N) m (x) is only an approximation to ψm(x) at values of x 6= xk. Also, the set Ĉk(x) depends on N . The matrix elements of the Hamiltonian, Eq. (4.1), in the physical space representation, denoted by P, is defined by ĤPij = (~ 2/2µ)K̂Pij + V̂ P ij , where K̂Pij = − 1√ η̂iη̂j ∫ ρ−1(b) ρ−1(a) Ĉi(x)Ĉ ′′ j (x) dx, (4.11a) V̂ Pij = 1√ η̂iη̂j ∫ ρ−1(b) ρ−1(a) Ĉi(x)V (x)Ĉj(x) dx. (4.11b) The kinetic energy matrix K̂Pij is symmetrized and discretized in a way similar to what is shown in Eqs. (4.7) and (4.8). The discretized Hamiltonian is given by ĤPij = ~ 2 2µ N∑ k=1 [ρ′(xk)] 2D (1) P [ki]D (1) P [kj] + [ V (xi)− ~ 2 2µ Ṽ (xi) ] δij . (4.12) 92 Chapter 4. An efficient mapped pseudospectral method for weakly bound states In Eq. (4.12), Ṽ is given by Eq. (4.8b), and the representation of the first derivative operator is [20] D (1) P [ij] = √ wiwj N∑ n=1 P ′n−1(ui)Pn−1(uj). The orthogonal transformation T̂ defined by T̂kn = √ η̂kφ̂n(xk) = √ wkPn−1(uk) relates the coefficient space and the physical space representations of the Hamiltonians in Eqs. (4.9) and (4.12). Thus, ĤPij = ∑N m=1 ∑N n=1 T̂imĤ C mnT̂jn, and ĉ(m)n = N∑ k=1 T̂kn √ η̂kψm(xk), (4.13a) √ η̂kψm(xk) = N∑ n=1 T̂knĉ (m) n . (4.13b) 4.3 The solution of the Schrödinger 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ödinger 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 sincu = sinπu πu if u 6= 0, 1 if u = 0. The set of sinc interpolating functions, or the Lagrange-sinc functions [3], are given by Sk(u) = sinc ( u− uk h ) (4.14) 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 uni- form 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 inter- polation. The mapped Lagrange-sinc functions are defined by Ĉk(x) = √ ρ′(x) ρ′(xk) Sk(ρ(x)) (4.15) and the orthogonality ∫ ρ−1(∞) ρ−1(−∞) Ĉi(x)Ĉ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 ≈ ∞∑ k=−∞ hg(xk). 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 ĤSij = ~ 2 2µ N∑ k=1 [ρ′(xk)] 2D (1) S [ki]D (1) S [kj] + [ V (xi)− ~ 2 2µ Ṽ S(xi) ] δij , (4.16) where S denotes the sinc representation and D (1) S [ij] = S ′ j(ui) = (−1)i−j (i− j)h if i 6= j, 0 if i = j, with Ṽ S(x) = 1 2 ρ′′′(x) ρ′(x) − 1 4 ( ρ′′(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) = x− b2 s2 , (4.17a) ρ(x) = s1 ln ( x− b2 s2 ) , (4.17b) ρ(x) = s1 arcsinh ( x− b2 s2 ) + b1. (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 s2 = (xmax − xmin) exp(uN/s1)− exp(u1/s1) , (4.18a) b2 = xmin − s2 exp(u1/s1). (4.18b) In addition, a parameter ŝ is introduced with s1 = (uN − u1)ŝ. (4.18c) The distribution of the grid points {xk}N−1k=2 is controlled indirectly by ŝ, where a small ŝ corresponds to a high concentration of grid points near 95 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 10 20 30 u 0 4 8 12 x = ρ- 1 (u ) A u1 uN xN = xmax x1 = xmin 0 10 20 30 u 0 4 8 12 x = ρ- 1 (u ) B u1 uN xN = xmax x1 = xmin 0 10 20 30 u 0 4 8 12 x = ρ- 1 (u ) C u1 uN xN = xmax x1 = xmin Figure 4.1: The mapping of collocation points using Eq. (4.17c) from La- guerre grid uk to a modified grid xk. (A) ŝ = 0.2, b2 = 5, x ∈ [1, 10]; (B) same as A except ŝ = 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 ŝ 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 2 ln ( c exp(u1/s1)− exp(uN/s1) c exp(−u1/s1)− exp(−uN/s1) ) , (4.19a) where c = xmax − b2 xmin − b2 , (4.19b) and s2 = xmax − xmin sinh ( uN−b1 s1 ) − sinh ( u1−b1 s1 ) . (4.19c) Equation (4.18c) is also used and ŝ controls the level of concentration of grid points near x = b2. A low value of ŝ corresponds to a high concentration, while a high value of ŝ gives a mapping close to a linear map. Figures 4.1A and 4.1B show the different distributions {xk} versus {uk} for different val- ues of ŝ. 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 re- duction 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 x6 − f8(bx)C8 x8 − f10(bx)C10 x10 , (4.21a) where the damping functions, fn, are given by f2n(bx) = 1− exp(−bx) 2n∑ 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 Å, 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). The discretized Hamiltonian, ĤLagij , for the Schrödinger equation, Eq. (4.1), is given by Eq. (4.12). The eigenfunction ψm(x) can be found by ψm(xk) = √ ρ′(xk)u 2 k exp(−uk) wk [Ψm]k, where uk and wk are the Laguerre quadrature points and weights, and Ψm is the eigenvector of ĤP, Eq. (4.12), associated with the eigenvalue Em. The 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. He2 Ne2 Ar2 A 41.96 199.5 748.3 b 2.523 2.458 2.031 C6 1.461 6.383 64.30 C8 14.11 90.34 1623 C10 183.6 1536 49060 2 3 4 5 6 x [Angstrom] -10 0 10 V( x) [m eV ] He2 Ne2 Ar2 A 2.6 2.8 3 3.2 x [Angstrom] -1 -0.8 -0.6 -0.4 -0.2 0 V( x) [m eV ] TT AS B TTY 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 -10 -8 -6 -4 -2 log(E) 0 2 4 6 8 δ 0 /pi TT Ar2 TT Ne2 TTY/AS He2 TT He2 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 − s2u1, 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 func- tion, given by w(u) = [ψMorse0 (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 Å, and α = 2.15179 Å−1. For Ar2, De = 0.0123670 eV, xe = 3.75669 Å, and α = 1.79072 Å−1. For this weight function, 1 2 w′′(u) w(u) − 1 4 ( w′(u) w(u) )2 = 2µ ~2 De[1− exp(−α(u− xe))]2 − α √ 2µ ~2 De + α2 4 . With the modified grid xk = ρ −1(uk), the Hamiltonian Ĥ QDM ij is given 101 Chapter 4. An efficient mapped pseudospectral method for weakly bound states by Eq. (4.12), where Ṽ QDM(x) = [ 2µ ~2 De[1− exp(−α(ρ(x) − xe))]2 − α √ 2µ ~2 De + α2 4 ] [ρ′(x)]2 + 1 2 ρ′′′(x) ρ′(x) − 1 4 ( ρ′′(x) ρ′(x) )2 + [ −2 √ 2µ ~2 De[1− exp(−α(ρ(x) − xe))] + α ] ρ′′(x). The eigenfunctions are ψm(xk) = η̂ −1/2 k [Ψm]k, where η̂k is given below Eq. (4.5), and [Ψm]k is the kth element of the eigenvector of Ĥ QDM associated 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 − s2u1 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], VTTY(x) = Dx 7/2β−1 exp(−2βx)− 12∑ n=3 f2n(bx) C2n x2n , (4.24) 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) 3C2n−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 VAS(x) = ǫV ∗ ( x rm ) , (4.25a) 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 189635.353 D 1.4088 α 10.70203539 ǫ/k 10.97 β -1.90740649 rm 2.9695 c6 1.34687065 Aa 0.0026000000 c8 0.41308398 r1 1.0035359490 c10 0.17060159 r2 1.4547903690 103 Chapter 4. An efficient mapped pseudospectral method for weakly bound states where V ∗(r) = A exp(−αr + βr2)− [ c6 r6 + c8 r8 + c10 r10 ] F (r) + Va(r), (4.25b) F (r) = exp [ − ( D r − 1 )2] if r < D, 1 if r ≥ D, (4.25c) Va(r) = Aa [ sin ( 2π r − r1 r2 − r1 − π 2 ) + 1 ] if r1 ≤ r ≤ r2, 0 if r < r1 or r > r2. (4.25d) with the constants as given in Table 4.2. The units of VAS and x are eV and Å, 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 Å−1. The other parameters used to define the weight function, Eq. (4.23), are De = 9.4659017(−4) eV and xe = 2.9720746 Å for VTTY, and De = 9.4532246(−4) eV and xe = 2.9695000 Å 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ödinger equa- tion 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 origi- nally 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 ex- actly 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 preci- sion. 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 Em = −α 2 ~ 2 2µ (√ 2µDe α~ − 1 2 −m )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 re- sults obtained by Willner et al [36] in Table 4.3. The QDM(arcsinh) provides results about two to three order of magni- tudes 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 re- ported by Willner et al have converged to only about two significant figures. 105 Chapter 4. An efficient mapped pseudospectral method for weakly bound states Table 4.3: Convergence of the Morse potential for Cs2. E exact m and ∆Em are in cm−1. m Eexactm ∆Em WDMSa QDM(arcsinh)b QDM(arcsinh)c N = 276 165 -3.26122 1.04(-4) -1.40(-6) -1.34(-6) 166 -2.10770 1.30(-4) -6.89(-7) -1.14(-6) 167 -1.20493 1.58(-4) -4.21(-7) -9.13(-7) 168 -0.552941 1.76(-4) -2.54(-7) -6.47(-7) 169 -0.151714 1.73(-4) -1.25(-7) -3.49(-7) 170 -1.25384(-3) 1.65(-4) 2.04(-4) 1.87(-4) N = 400 170 -1.25384(-3) -1.11(-5) -6.22(-13) -6.84(-13) a Fourier-sine results by Willner, Dulieu, and Masnou-Seeuws [36]. b ŝ = 0.2, b2 = 10, x ∈ [6, 80] for N = 276; ŝ = 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 ŝ = 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 fic- titious eigenvalues, called “ghost levels”, which do not belong to the spec- trum 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 for- mula 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- tation. It can be done by constructing ĤCmn, given by Eq. (4.9), only for 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 ′) m (x) = √ ρ′(x)w(ρ(x)) N ′∑ n=1 ĉ(m)n Pn−1(ρ(x)), where ĉ (m) n is the nth component of the corresponding eigenvector. Note that the unitary transformation between the coefficient space representa- tion 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 de- pends on this error should be sensitive to this change. Surprisingly, the 107 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 260 270 280 290 N -10 -8 -6 -4 -2 0 E A 260 270 280 290 N -10 -8 -6 -4 -2 0 E B 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. 0 100 200 n -16 -12 -8 -4 0 A lo g |ĉ n | 0 20 40 60 80 x -1 0 1 ψ( x) B 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ödinger 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)m = log ∣∣∣∣∣E (N) m − Eexactm Eexactm ∣∣∣∣∣ , is approximately the number of significant figures of the eigenvalue estimates E (N) m . The exact eigenvalues, Eexactm , are calculated with the mapped sinc 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. TTY He2 a AS He2 b TT Ne2 c TT Ar2 c E0 -1.08286(-4) -1.07315(-4) -2.07329 -10.4703 E1 -0.362465 -7.21707 E2 -1.56865(-3) -4.65637 E3 -2.74616 E4 -1.42543 E5 -0.608341 E6 -0.183899 E7 -0.0246251 a Tang-Toennies-Yiu, Eq. (4.24). b Aziz-Slaman LM2M2, Eq. (4.25). c Tang-Toennies, Eq. (4.21). 0 2 4 6 8 10 x ψ 0 (x) A QDM Lag 0 2 4 6 8 10 x ψ 0 (x) B sinc sinc(arcsinh) 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 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 5 10 15 20 25 x ψ 2 (x) A Lag Lag(ln) 0 5 10 15 20 25 x ψ 2 (x) B QDM QDM(ln) 0 5 10 15 20 25 x ψ 2 (x) C sinc(arcsinh) 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): ŝ = 0.25, x ∈ [2, 150]; (B) QDM: s2 = 1.8; QDM(ln): ŝ = 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 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 10 20 30 40 N -16 -12 -8 -4 0 ε 0 ( N ) sinc Lag QDM sinc(arcsinh) A 0 20 40 60 80 100 N -12 -8 -4 0 ε 2 (N ) sinc sinc(arcsinh) Lag Lag(ln) QDM QDM(ln) B Figure 4.8: Variation of ǫ (N) m for m = 0 and 2 versus the number of collo- cation points N for Ne2. Lag: xmin = 2 for all cases, (A) s2 = 0.06, (B), s2 = 0.3; Lag(ln): (B) ŝ = 0.25, x ∈ [2, 150]; QDM: xmin = u1 for all cases, (A) s2 = 1, (B) s2 = 1.8; QDM(ln): (B) ŝ = 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 Chapter 4. An efficient mapped pseudospectral method for weakly bound states improvement of the mapped sinc results as illustrated in Fig. 4.6B. The su- periority 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 be- yond 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 im- proves 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 La- guerre 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 im- proved, 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. Figure 4.9 shows the variation of the coefficients ĉ (m) 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- guerre, the variation of ĉ (m) 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 0 20 40 60 80 n -8 -4 0 A QDM QDM(ln) lo g ∣ ∣ ∣ĉ(2 ) n ∣ ∣ ∣ 0 20 40 60 80 n -8 -4 0 B Lag Lag(ln) lo g ∣ ∣ ∣ĉ(2 ) n ∣ ∣ ∣ Figure 4.9: Variation of the expansion coefficients ĉ (2) 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 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 10 20 30 40 N -16 -12 -8 -4 0 ε 0 ( N ) sinc sinc(arcsinh) Lag QDM A 0 20 40 60 N -16 -12 -8 -4 0 ε 3 (N ) sinc sinc(arcsinh) Lag QDM B 0 20 40 60 80 N -16 -12 -8 -4 0 ε 5 ( N ) sinc sinc(arcsinh) Lag Lag(ln) QDM QDM(ln) C 0 20 40 60 80 100 N -12 -8 -4 0 ε 7 (N ) sinc sinc(arcsinh) Lag Lag(ln) QDM QDM(ln) D Figure 4.10: Variation of ǫ (N) m for m = 0, 3, 5, and 7 versus the number of collocation pointsN 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) ŝ = 0.7, x ∈ [2.7, 13], (D) ŝ = 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) ŝ = 1.1, x ∈ [u1, 8.5], (C) ŝ = 0.6, x ∈ [u1, 13], (D) ŝ = 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 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 conver- gence 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 pro- vides 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 rmD, rmr1 and rmr2. Thus, with the mapped sinc grid we can only evaluate Eexact0 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, respec- tively. 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 Eexact0 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 poten- tial. N Laga Lag(arcsinh)b QDMc QDM(arcsinh)d sinc(arcsinh)e 10 -0.4077 0.0428 3.1934 -0.3581 -28.1437 20 0.5817 -0.1044 0.2667 -0.1056 -0.3888 30 0.0694 -0.1073 -0.0601 -0.1073 -0.1162 40 -0.0345 -0.0934 -0.1076 50 -0.0723 -0.1055 -0.1073 100 -0.1064 -0.1078 150 -0.1074 -0.1068 200 -0.1073 -0.1070 250 -0.1073 -0.1071 300 -0.1073 -0.1074 350 -0.1074 -0.1071 400 -0.1073 -0.1074 a s2 = 0.7, xmin = 1.5. b ŝ = 0.11, b2 = 4.3, x ∈ [1.5, 800]. c s2 = 0.5, xmin = u1. d ŝ = 0.12, b2 = 4.3, x ∈ [1.5, 800] e s2 = 0.8, b2 = 3, x ∈ [1.5, 800] 117 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 40 80 120 160 N -8 -6 -4 -2 0 2 ε 0 (N ) QDM(arcsinh) Lag(arcsinh) sinc(arcsinh) QDM Lag A 0 20 40 60 80 100 N -12 -8 -4 0 ε 0 (N ) QDM(arcsinh) Lag(arcsinh) sinc(arcsinh) QDM Lag B Figure 4.11: Variation of ǫ (N) 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): ŝ = 0.14, b2 = 3.1, x ∈ [1.1, 1800]; QDM: s2 = 0.6, xmin = u1; QDM(arcsinh): ŝ = 0.14, b2 = 3.5, x ∈ [1.1, 1800]; sinc(arcsinh): s2 = 0.6, b2 = 1.8, x ∈ [1.1, 1800]. 118 Chapter 4. An efficient mapped pseudospectral method for weakly bound states 0 200 400 600 800 x ψ 0 (x) A Lag Lag(arcsinh) 0 200 400 600 800 x ψ 0 (x) B QDM QDM(arcsinh) 0 200 400 600 800 x ψ 0 (x) C sinc(arcsinh) 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): ŝ = 0.11, b2 = 4.3, x ∈ [1.5, 800]; (B) QDM: s2 = 0.5, xmin = u1; QDM(arcsinh): ŝ = 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 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 ap- pear. 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 0 20 40 60 80 100 N 2 4 6 8 - lo g(- E) E0 ghost level A 0 20 40 60 80 100 N 2 4 6 8 - lo g(- E) E0 ghost level B 0 20 40 60 80 100 N 4 6 8 - lo g(- E) E0 ghost level C 0 20 40 60 80 100 N 4 6 8 - lo g(- E) E0 ghost level D 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 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 representa- tion 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 optimiz- ing grid representations: the mapped Fourier method. Phys. Rev. E, 53:1217, 1996. 122 Bibliography [12] W. Gautschi. Orthogonal Polynomials: Computation and Approxima- tion. Oxford, New York, 2004. [13] R. E. Grisenti, W. Schöllkopf, J. P. Toennies, G. C. Hegerfeldt, T. Köhler, 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éonard, M. Walhout, A. P. Mosk, T. Müller, M. Leduc, and C. Cohen-Tannoudji. Giant helium dimers produced by photoassoci- ation 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 weakly- bound 4Hex 20NeyH (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 dis- cretization method in the solution of the Schrödinger 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 Differen- tial 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öllkopf 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 func- tions. 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ödinger 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 H2O. J. Chem. Phys., 106:6885, 1997. [35] E. T. Whittaker. On the functions which are represented by the expan- sions 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 time- dependent 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 prob- 125 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, Ĉ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)Ĉj(x) is referred to as the interpolant of g(x). The Gaussian function G(x;µ, σ) = 1 σ √ 2π exp ( −1 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 in- terpolated 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 -2 -1 0 1 2 x 0 2 4 6 g(x ) A -2 -1 0 1 2 x -0.1 -0.05 0 0.05 0.1 g(x ) - f(x ) B Figure 5.1: Unmapped sinc interpolation of Eq. (5.2) using 160 grid points. (A) Grid point distribution; (B) Error of the interpolation. 127 Chapter 5. Conclusions and future directions The equidistribution principle is a method for placing the grid points, {xj}Nj=1, so that the quantity ∫ xj+1 xj M(x) dx is equally distributed [1]. The 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 ŝ 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 redistri- bution 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 way that the quantity ∫ xj+1 xj M(x;α) dx is distributed as evenly as possible. I therefore take the values of ŝ and b2 which minimize the function E(ŝ, b2) = N−1∑ j=1 (∫ xj+1 xj M(x;α) dx )2 where the physical grid is defined by xj = ρ −1(uj). The integrals and the minimization are evaluated by Riemann sum and steepest descent, respec- tively. For the interpolation of g(x), the above procedure yields the values ŝ = 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 -2 -1 0 1 2 x 0 2 4 6 g(x ) A -2 -1 0 1 2 x -0.1 -0.05 0 0.05 0.1 g(x ) - f(x ) B 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 interpo- lation. 129 Chapter 5. Conclusions and future directions Implicitly defined mapping An appropriately distributed grid can be constructed directly by the equidis- tribution principle without a predefined mapping function such as the hyper- bolic 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 pseudospec- tral methods. Since the mapping is constructed implicitly, it is referred to as an implicit mapping [1]. First, a reference grid {x̂j}Nj=1 is constructed by satisfying ∫ x̂j+1 x̂j M(x;α) dx = 1 N − 1 ∫ xmax xmin M(x;α) dx (5.4) 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) is required. By defining Û(x) as a piecewise linear function joining the points (x̂j , uj), the map ρ(x) is defined by fitting Û(x) by least squares. I propose ρ(x) to be of the form ρ(x) = x+ Nfit∑ n=1 cn sin ( nπ 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 E = ∫ xmax xmin { [ρ(x)− Û(x)]2 + β[ρ′′(x)]2 } 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 oscil- lations. 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 -2 -1 0 1 2 u -2 -1 0 1 2 x = ρ- 1 (u ) A -2 -1 0 1 2 u -2 -1 0 1 2 x = ρ- 1 (u ) B Figure 5.3: The reference grid x̂j (crosses) and the map ρ(x). The param- eters are N = 30, Nfit = 24, and α = 100. (A) β = 0; (B) β = 3× 10−5. 131 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 + β π 4n4 (xmax − xmin)3 )−1 × ∫ xmax xmin [Û (x)− x] sin ( nπ x− xmin xmax − xmin ) dx. Since Û(x) − x is a piecewise linear function, the integral on the right side can be evaluated analytically. The new physical grid, xj, is calculated nu- merically 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 inter- polation 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 time- dependent 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 be- haviours. 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 rein- terpolation 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 -2 -1 0 1 2 x 0 2 4 6 g(x ) A -2 -1 0 1 2 x -0.1 -0.05 0 0.05 0.1 g(x ) - f(x ) B 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 Chapter 5. Conclusions and future directions costs, a new grid is used only when it differs from the old grid by a signifi- cant amount. The self-adjoint form of the Smoluchowski equation given by Eq. (1.10) is ∂p ∂t = ǫ ν √ P0 ∂ ∂x [ P0 ∂ ∂x ( p√ P0 )] = LSp (5.6) when P (x, t) = √ P0(x)p(x, t). The solution p(x, t) is represented as an interpolation p(x, t) = N∑ j=1 aj(t) Cj(x) ηj , (5.7) where aj(t) = √ ηjp(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](x), x [n] j , and C [n] j (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. The grid points x [n] j 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 differ- ence 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 x [n] j , the monitor function can be considered as a piecewise function M [n+1](x;α) = [ 1 + α̂ ∣∣∣∣∣p(x [n] j+1, tn+1) + p(x [n] j , tn+1) 2 ∣∣∣∣∣ + α̂ ∣∣∣∣∣p(x [n] j+1, tn+1)− p(x[n]j , tn+1) x [n] j+1 − x[n]j ∣∣∣∣∣ ] 1 2 for x [n] j < x < x [n] j+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 re- 134 Chapter 5. Conclusions and future directions placed by the difference formula. Because the magnitude of p(x, t) changes over time, a new parameter α̂ = α/[maxj p(x [n] j , tn+1)] is introduced to keep the monitor function consistent at each time step. The reference grid x̂ [n+1] j is generated by the equidistribution principle ∫ x̂[n+1]j+1 x̂ [n+1] j M [n+1](x;α) dx = 1 N − 1 ∫ xmax xmin M [n+1](x;α) dx (5.8) as in Eq. (5.4). If the difference between the reference grids x̂ [n+1] j deter- mined by Eq. (5.8) and the reference grid x̂ [n] j used at the previous time step is small, i.e. max j ∣∣∣x̂[n+1]j − x̂[n]j ∣∣∣ ≤ ∆x̂ for some tolerance level ∆x̂, no new grid is required. One will therefore set x [n+1] j = x [n] j , ρ[n+1](x) = ρ[n](x), a [n+1] j (tn+1) = a [n] j (tn+1), and copy the values of all other quantities with [n] to [n+ 1]. The values of x̂ [n+1] j found by Eq. (5.8) will be abandoned, and will be set to x̂ [n+1] j = x̂ [n] j . On the other hand, if maxj ∣∣∣x̂[n+1]j − x̂[n]j ∣∣∣ exceeds the tolerance level, x̂ [n+1] j is used to construct ρ[n+1](x) as described in Section 5.1.1, and a new physical grid x [n+1] j is generated. All other quantities associated with x [n+1] j are required to be updated. The solution p(x [n+1] j , tn+1) is interpolated by p(x [n+1] j , tn+1) = N∑ i=1 a [n] i (tn+1)√ η̂ [n] i C [n] i (x [n+1] j ) and the updated expansion coefficient set to a [n+1] j (tn+1) = √ η̂ [n+1] j p(x [n+1] j , 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 -2 -1 0 1 2 x -1 0 1 2 3 4 P( x, 0) A -2 -1 0 1 2 x -1 0 1 2 3 4 P( x, 1. 6) C -2 -1 0 1 2 x -1 0 1 2 3 4 P( x, 0. 1) B -2 -1 0 1 2 x -1 0 1 2 3 4 P( x, 25 . 6) D 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 Chapter 5. Conclusions and future directions -2 -1 0 1 2 x 0.1 1 10 100 1000 t Figure 5.6: The evolution of grid points over time for the solution of Smolu- chowski equation. -2 -1 0 1 2 x -8 -4 0 4 8 P r ef (x, 1. 6) - P( x, 1. 6) A ×10-4 moving mesh unmapped -2 -1 0 1 2 x -8 -4 0 4 8 P r ef (x, 25 . 6) - P( x, 25 . 6) B ×10-4 moving mesh unmapped Figure 5.7: The error of the solution of Smoluchowski equation for in- terpolatory moving mesh method (N = 40) and unmapping sinc method (N = 141). (A) t = 1.6; (B) t = 25.6. 137 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 rep- resents 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 com- pared 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 solu- tion, 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 cal- culating the eigenvalues of highly diffuse states of the Schrödinger equation. In this section, I will demonstrate the use of implicit mapping for the evalu- ation 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 L S 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 Smolu- chowski matrix operator. In this example, N = 20 uniform grid points in the interval [−1.2, 1.2] are used. Among the eigenfunctions found by L̃Sij , 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 M(x;α) = [ 1 + α ∣∣∣∣∣ ψ̃(uj+1) + ψ̃(uj)2 ∣∣∣∣∣+ α ∣∣∣∣∣ ψ̃(uj+1)− ψ̃(uj)uj+1 − uj ∣∣∣∣∣ ] 1 2 for uj < x < uj+1 and j = 1, . . . , N − 1. Compared with Eq. (5.3), the sec- ond 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 compar- ison. 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 eigen- values 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 -1 0 1 x -0.4 0 0.4 0.8 1.2 A p̃ A (x ) -1 0 1 x -0.4 0 0.4 0.8 1.2 B p̃ B (x ) -1 0 1 x -0.4 0 0.4 0.8 1.2 C p̃ C (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 -1 0 1 u -1 0 1 A x = ρ − 1 A (u ) -1 0 1 u -1 0 1 B x = ρ − 1 B (u ) -1 0 1 u -1 0 1 C x = ρ − 1 C (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 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10 ρA(x) 12 -1.006404 -2.47368 -3.41230 -7.30289 -8.52820 14 -0.997546 -2.10652 -3.20082 -5.19971 -6.26385 16 -0.997455 -2.00468 -2.99653 -4.31843 -5.33439 18 -0.996985 -1.98875 -2.97411 -4.01143 -4.99822 20 -0.996982 -1.98789 -2.97266 -3.95399 -4.92765 22 -2.97265 -3.95118 -4.92340 24 -4.92337 26 -4.92340 ρB(x) 20 -2.52936 -2.54618 -4.34540 -4.59095 24 -1.99579 -2.03057 -4.65541 -4.74760 28 -2.00692 -2.03016 -3.98030 -3.99223 32 -1.98212 -1.99033 -3.99228 -4.00816 36 -1.98927 -1.98852 -3.95093 -3.95251 40 -1.98806 -1.98805 -3.95217 -3.95238 44 -1.98789 -1.98789 -3.95124 -3.95139 48 -3.95118 -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 λ6 λ7 λ8 λ9 λ10 ρC(x) 24 -1.021761 -2.40046 -2.36432 -2.86485 -3.70659 -4.41502 -4.67242 -8.34631 -9.07370 28 -1.001783 -2.07855 -2.12363 -2.32266 -3.24054 -4.72207 -4.69602 -6.04777 -7.12278 32 -0.997608 -2.00221 -1.99835 -2.08958 -3.03982 -4.03080 -4.10496 -4.85710 -5.69441 36 -0.997030 -1.99828 -2.00395 -2.01674 -2.98359 -3.98455 -3.97918 -4.27327 -5.14256 40 -0.996984 -1.98830 -1.98911 -1.99162 -2.97353 -3.97021 -3.97626 -4.03525 -4.96543 44 -0.996982 -1.98797 -1.98796 -1.98829 -2.97266 -3.95165 -3.95320 -3.96388 -4.92591 48 -1.98817 -1.98799 -1.98790 -2.97265 -3.95134 -3.95132 -3.95221 -4.92318 52 -1.98807 -1.98792 -1.98789 -3.95143 -3.95137 -3.95121 -4.92337 56 -1.98788 -1.98789 -3.95125 -3.95124 -3.95118 -4.92339 60 -1.98789 -3.95119 -3.95119 -4.92340 64 -3.95118 -3.95118 unmapped 40 -1.009628 -2.60022 -2.63149 -2.40029 -3.39238 -7.66165 -7.58504 -6.01067 -6.96642 48 -0.998489 -2.59873 -2.57044 -2.08817 -3.07145 -4.53141 -4.52771 -4.73743 -5.69731 56 -0.997096 -2.03357 -2.02456 -2.00389 -2.98787 -4.56529 -4.59437 -4.16818 -5.13911 64 -0.996987 -2.02917 -2.03092 -1.98931 -2.97401 -4.21517 -4.18406 -3.98992 -4.96098 72 -0.996982 -2.00925 -2.00198 -1.98796 -2.97272 -3.96703 -3.96886 -3.95542 -4.92706 80 -1.98808 -1.98807 -1.98789 -2.97265 -3.97783 -3.97444 -3.95138 -4.92359 88 -1.98949 -1.98834 -3.95461 -3.95325 -3.95118 -4.92340 96 -1.98741 -1.98792 -3.95144 -3.95140 104 -1.98789 -1.98789 -3.95139 -3.95125 112 -3.95118 -3.95118 120 -3.95117 128 -3.95118 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 [ −1 ǫ ( y2 2 + U(x) )] where Kxy is the normalization constant. Since Peq(x, y) can be written as Peq(x, y) = K exp(−y2/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(−y2/2ǫ) for y. By setting P (x, y) = √ P0(x)p(x, y), Eq. (1.6) becomes ∂p ∂t = − P ′ 0 2P0 yp− y ∂p ∂x + ν ∂ ∂y (yp) + U ′ ∂p ∂y + νǫ ∂2p ∂y2 = LKp (5.9) We consider the expansion of p(x, y) as p(x, y) = Nx∑ i=1 Ny∑ j=1 aij Ci(x)√ ηi CHermj (y)√ ηHermj (5.10) where aij = √ ηiηHermj p(xi, yj), η Herm j = wj/w(yj), C Herm j (y) are the scaled Hermite interpolating functions CHermj (y) = √ w(y) w(yj) N∏ k=1 k 6=j y − yk yj − yk , 144 Chapter 5. Conclusions and future directions w(y) = exp(−y2/2ǫ), yj and wj are the quadrature points and weights asso- ciated 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 + ( α ψ̃(xi+1)− ψ̃(xi) 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 y- directions, 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 pseu- dospectral methods for the solutions to the Fokker-Planck equation, Eq. (2.1), and the Schrödinger 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 cowork- ers [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. eigenvalue ρA(x) ρB(x) ρC(x) unmapped λ2 −0.501144 ± 1.320455i (40, 10) (56, 10) (104, 10) λ3 −0.501144 ± 1.320455i (40, 12) (56, 10) (104, 12) λ4 −0.617431 (32, 20) (76, 20) (112, 20) λ5 −0.998152 (40, 12) (60, 16) (112, 12) λ6 −0.998152 (48, 14) (60, 16) (112, 12) λ7 −1.00554 ± 2.63851i (40, 14) (52, 14) (104, 12) λ8 −1.00554 ± 2.63851i (40, 14) (52, 14) (104, 12) λ9 −1.23291 (36, 20) (84, 20) (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 corre- sponding 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 chap- ter 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 expo- nential 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ödinger 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ödinger 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ödinger equation a detailed com- parison 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 ground-state eigenfunction, ψ0(x), of the Schrödinger equation. In this case, ψ0(x) is represented exactly with the QDM, and the corresponding eigen- value is exact. This representation provides a very rapid convergence rate for the lower eigenstates. However, the convergence rate of the higher eigen- values 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 eigen- value 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 approx- imates ψ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 appro- priate 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 con- structing the weight function. Despite the fact that this eigenfunction was piecewise defined with discontinuous second derivative, the resulting poly- nomial 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 ap- plications focused on the vibrational states of diatoms with various poten- tials [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 Hamilto- nian 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 or- der 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 con- struct a mapping based on the behaviour of the eigenfunctions. For example in the Schrödinger 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 im- plicit 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. Conclusions and future directions 5.3 Future directions The interpolation with implicit mapping based on the equidistribution prin- ciple 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 di- minished, 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 re- gion 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 so- lution. 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 g(x, y) = Nx∑ i=1 Ny∑ j=1 g(xi, yj)C (x) i (x)C (y) j (y) where C (x) i and C (y) j 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 equidis- tribution principle for 2D adaptive grids is discussed in Ref. 6. In order to apply the adaptive grid to pseudospectral methods, a 2D mapping func- tion 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 Equa- tions. 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 dimen- sions. 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 Differen- tial Equations. SIAM, Philadelphia, 1992. [9] A. Madzvamuse and P. K. Maini. Velocity-induced numerical solu- tions 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 Ap- plication. Springer-Verlag, New York, 1984. [11] B. D. Shizgal and H. Chen. The quadrature discretization method (QDM) in the solution of the Schrödinger 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 Fokker- Planck 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 H2O. 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) = √ γnFn(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−1Qn−2(x) (A.1) for n ≥ 2, where Q0(x) = 1, Q1(x) = x− α0, αn = γ−1n ∫ 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 w(x) d dx [ w(x)B(x) dφm(x) dx ] . The matrix elements of the Fokker-Planck operator, Lmn, are thus Lmn = 1√ γm−1γn−1 ∫ ∞ 0 Qm−1(x) d dx [ w(x)B(x)Q′n−1(x) ] 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 w(x)Qm−1(x)Q ′ n−1(x) dx − 2√ γm−1γn−1 ∫ ∞ 0 x2w(x)Qm−1(x)Q ′ n−1(x) dx + 1√ γm−1γn−1 ∫ ∞ 0 xw(x)Qm−1(x)Q ′′ n−1(x) dx. (A.4) 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 x2Q′n−1 in the Qn basis, thus we write x2Q′n−1(x) = bnQn(x) + bn−1Qn−1(x) + · · · + b0Q0(x). (A.5) Hence, this integral vanishes when m > n+ 1. For m = n+ 1, −2√ γnγn−1 ∫ ∞ 0 x2w(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 x n in x2Q′n−1(x) is n− 1. So bn = n− 1. For m = n, the integral is given by −2 γn−1 ∫ ∞ 0 x2w(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 re- currence relation in Eq. (A.1). We first express Qn(x) as the polynomial Qn(x) = ∑n k=0An,kx k, where the first three coefficients are An,n = 1, An,n−1 = − ∑n−1 k=0 αk and An,n−2 = ∑n−1 k=0 ∑k−1 l=0 αkαl − ∑n−1 m=1 βm, so that from the expansion in Eq. (A.5), x2Q′n−1(x) = (n− 1)xn + (n− 2)An−1,n−2xn−1 + · · ·+An−1,1x2, bnQn(x) = bnx n + bnAn,n−1x n−1 + · · · + bnAn,0, bn−1Qn−1(x) = bn−1x n−1 + · · ·+ bn−1An−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 − bnAn,n−1 = −(n− 2) n−2∑ k=0 αk + (n− 1) n−1∑ k=0 αk = (n− 1)αn−1 + n−2∑ k=0 αk. 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 Lmn = −2(n− 1)αn−1 − 2 ∑n−2 k=0 αk when m = n > 1, −2(n− 1)√βn when m = n+ 1, −2(m− 1)√βm when m = n− 1, 0 otherwise. (A.8) The matrix operator L is given explicitly by L = 0 0 0 0 0 0 · · · 0 −4.9761 −1.3032 0 0 0 · · · 0 −1.3032 −11.258 −3.1214 0 0 · · · 0 0 −3.1214 −18.689 −5.3220 0 · · · 0 0 0 −5.3220 −27.137 −7.8425 . . . 0 0 0 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 Lmn = − 1√ γ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+1Qn+1(x) + bn−1Qn−1(x) + bn−3Qn−3(x) + · · ·+ { b1Q1(x) for even n, b0Q0(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√ γn+1γn−1 ∫ ∞ −∞ w(x)[Qn+1(x)] 2 dx = −bn+1√ γn+1γn−1 γn+1 = −bn+1 √ βn+1βn. 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−3xn−4 + · · · ] = (n− 1)xn+1 + [(n− 3)An−1,n−3 − (n− 1)]xn−1 + · · · bn+1Qn+1(x) = bn+1x n+1 + bn+1An+1,n−1x n−1 + · · · bn−1Qn−1(x) = bn−1x n−1 + · · · Equating the coefficients of xn−1 gives bn−1 = [(n− 3)An−1,n−3 − (n− 1)]− bnAn+1,n−1 = −(n− 3) n−2∑ k=1 βk − (n− 1) + (n− 1) n∑ k=1 βk = (n− 1)(βn + βn−1 − 1) + 2 n−2∑ k=1 βk. (B.2) 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 Lmn = −(β2 + β1 − 1) when m = n = 2, −(n− 1)(βn + βn−1 − 1)− 2 ∑n−2 k=1 βk when m = n > 2, −(n− 1)√βn+1βn when m = n+ 2, −(m− 1)√βm+1βm when m = n− 2, 0 otherwise, (B.3) 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 L = 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 . . . 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 Leven = 0 0 0 0 0 · · · 0 −1.9577 −0.4024 0 0 · · · 0 −0.4024 −3.9122 −0.9861 0 · · · 0 0 −0.9861 −5.8627 −1.7090 . . . 0 0 0 −1.7090 −7.8079 . . . ... ... ... . . . . . . . . . , (B.5) Lodd = −0.0101 −0.1407 0 0 0 · · · −0.1407 −2.0518 −0.5967 0 0 · · · 0 −0.5967 −4.0964 −1.2175 0 · · · 0 0 −1.2175 −6.1445 −1.9673 . . . 0 0 0 −1.9673 −8.1968 . . . ... ... ... . . . . . . . . . . (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ödinger Operator If the weight function w(y) satisfies V (y) = V0(y) as defined in Eq. (2.4), the self-adjoint Schrödinger operator in Eq. (2.3) is given by Hψm(y) = − 1√ w(y) d dy [ w(y) d dy ( ψ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ödinger 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ödinger 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,0Q0(x) +Cn,1Q1(x) + · · · + Cn,n−2Qn−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)∑ k=0 Cm,kCn,k (C.2) 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ödinger Operator Eq. (C.1) and is given by H = 0 0 0 0 0 0 · · · 0 4.7255 0.4163 0.3626 −0.1302 0.0750 · · · 0 0.4163 10.246 0.6053 1.1485 −0.4904 · · · 0 0.3626 0.6053 16.794 0.3615 2.4556 · · · 0 −0.1302 1.1485 0.3615 24.491 −0.4420 · · · 0 0.0750 −0.4904 2.4556 −0.4420 33.408 · · · ... ... ... ... ... ... . . . . (C.3) 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 func- tions orthogonal with weight function web(y) = P0(y), y ∈ [0,∞), for ǫ = 0.01 is given by Heven = 0 0 0 0 0 0 · · · 0 1.9204 0.3054 0.0137 −5.8(−8) 2.8(−7) · · · 0 0.3054 3.7174 0.7557 0.0412 −5.2(−6) · · · 0 0.0137 0.7557 5.3693 1.3264 0.0866 · · · 0 −5.8(−8) 0.0412 1.3264 6.8363 2.0169 · · · 0 2.8(−7) −5.2(−6) 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 wob (y) = y 2P0(y), y ∈ [0,∞), 162 Appendix C. Matrix Representation of the Schrödinger Operator for ǫ = 0.01 is given by Hodd = 0.0101 0.1405 0.0072 0 0 0 · · · 0.1405 1.9960 0.4987 0.0263 −5.1(−5) 1.0(−5) · · · 0.0072 0.4987 3.8705 0.9835 0.0588 −2.0(−4) · · · 0 0.0263 0.9835 5.6186 1.5740 0.1078 · · · 0 −5.1(−5) 0.0588 1.5740 7.2185 2.2626 · · · 0 1.0(−5) −2.0(−4) .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− λ0I)× 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 interpo- lated with the implicit mapping generated based on 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](x), x [n] j , and C [n] j (x), respectively. The discretization of Eq. (5.6) gives the Smoluchowski operator matrix LS[n][ij] = − ǫ ν N∑ k=1 [ D (1) ki ρ ′ [n](x [n] k ) + ĝ [n](x [n] k )δki ] × [ D (1) kj ρ ′ [n](x [n] k ) + ĝ [n](x [n] k )δkj ] , (D.1) whereD (1) ij is given by Eq. (2.6) and ĝ[n](x) = 1 2 [ ρ′′[n](x)/ρ ′ [n](x) + P ′ 0(x)/P0(x) ] . With the quadrature weights η [n] j , the coefficients aj(t) from Eq. (5.7) asso- ciated with x [n] j are denoted by a [n] j (t) = √ η [n] j p(x [n] j , t). Therefore, by using the expansion p(x, t) = N∑ j=1 a [n] j (t) C [n] j (x)√ η [n] j , the Smoluchowski equation in Eq. (5.6) becomes da [n] j /dt = ∑N j=1 L S [n][ij]a [n] j . Based on the backward difference method, the coefficients a [n] j (t) are evalu- ated by solving a [n] i (tn) = N∑ j=1 [ δij −∆tnLS[n][ij] ] a [n] j (tn+1) (D.2) 164 Appendix D. Moving mesh method for the time dependent Smoluchowski equation for a [n] j (tn+1), where ∆tn = tn − tn−1. If a new grid is used at tn+1, i.e. x [n+1] j 6= x[n]j , a[n+1]j (tn+1) has to be found. The solution of Eq. (5.6) and (1.7) is therefore given by p(x, tn) ≈ N∑ j=1 a [n] j (tn)√ η [n] j C [n] j (x). and P (x, tn) = √ 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ácek, 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. [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
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066507/manifest