THE QUADRATURE DISCRETIZATION METHOD AND ITS APPLICATIONS By Heli Chen B. Sc. (Mechanics), Peking University, 1987 M. Sc. (Mechanical Engineering), University of Calgary, 1993 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of Mathematics) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1998 © Heli Chen, 1998 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The University of British Columbia Vancouver, Canada Date Oct, /I , rfU DE-6 (2/88) Abstract A discretization method referred to as the quadrature discretization method is in troduced and studied in this thesis. The quadrature discretization method is a spectral method based on a grid of points that coincide with the points of a quadrature. The quadrature is based on a set of nonclassical polynomials orthogonal with respect to some weight function over some interval. The method is flexible with respect to the choice of the weight function and quadrature points so that the optimum accuracy and conver gence in the solution of differential equations and/or partial differential equations can be obtained. The properties of the quadrature discretization method are studied and compared with classical spectral methods as well as finite difference methods. Several analytic model problems in one and three dimension are studied. The quadrature dis cretization method competes well with classical spectral methods and is far more superior than the finite difference method. In some cases it provides significant improvement in the accuracy and convergence of the solution of the problem over other methods. The main objective of the quadrature discretization method is to determine the weight function that defines the polynomial basis set and hence the grid points that provide opti mum convergence in a given application. The quadrature discretization method is applied to a large class of time dependent Fokker-Planck equations. Several weight functions are used and the results are compared with several other methods. The weight functions that have often provided rapid convergence of the eigenvalues and eigenfunctions of the Fokker-Planck operator are the steady solutions at infinite time. The quadrature discretization method is also employed in the solution of Schrodinger equations. The weight functions that are used are related to the ground state wave ii functions if known, or some approximate form. The eigenvalues and eigenfunctions of four different potential functions discussed extensively in the literature are calculated and compared with the published values. The eigen-problem of a two-dimensional Schrodinger equation with the Henon-Heiles potential is also calculated with the quadrature dis cretization method. The rate of convergence of the eigenvalues and eigenfunctions of the Schrodinger equations is very rapid with this approach. in Table of Contents Abstract ii Table of Contents v List of Tables vi List of Figures ix Acknowledgement xi 1 Introduction 1 2 The Quadrature Discretization Method (QDM) 14 2.1 Classical polynomials and the pseudospectral method 15 2.2 The quadrature discretization method 28 2.2.1 Construction of derivative matrices2.2.2 Generation of orthogonal polynomials, quadrature weights and points 33 2.3 Time integration 42 2.4 Numerical tests and results 45 2.4.1 One-dimensional Poisson equation 42.4.2 An analytic example 46 2.4.3 Three-dimensional Poisson problems . 60 3 Applications of the QDM to the Solution of Fokker-Planck Equation 66 3.1 Introduction 6iv 3.2 The QDM matrix representation of the Fokker-Planck operator 69 3.3 Eigenfunction expansion 72 3.4 Temporal discretization 3 3.5 Applications and results 4 3.6 Summary 102 4 Applications of the QDM to the Solution of Schrodinger equation 103 4.1 Introduction 104.2 The equivalence of the Fokker-Planck eigenvalue problem and the Schrodinger equation 106 4.3 The QDM matrix representation of the Schrodinger equations 108 4.4 Applications and results 110 4.5 Summary 145 Summary and future work 142 Bibliography 148 v List of Tables 2.1 Errors for the numerical approximation of u" = sin(7rx) 18 2.2 Maximum error in an and 8n calculated by the Stieltjes's procedure ... 36 2.3 Maximum error in the numerical solution of an and (3n for the Hermite polynomials 38 2.4 Maximum error for the numerical approximation of points and weights . 40 2.5 Standard error for the numerical approximation of u" = sin(7rx) 45 2.6 Numerical error in the eigenvalues for the Hermite and QDM methods . . 48 2.7 Numerical error in the eigenvalues for the Legendre method 49 2.8 Domain generated from the QDM and Hermite weight function 52 2.9 Maximum time step in the FE time integration 55 2.10 Comparison of the numerical error of the eigenvalues for the Hermite and Legendre methods 59 2.11 Examples of Poisson problems 61 2.12 Standard Error and CPU time for the numerical approximations for ex ample A 62 2.13 Comparison of the standard Errors and CPU times for the three numerical approximations 3 3.1 Convergence of the eigenvalues for the quartic potential 78 3.2 Comparison of the rate of convergence of the eigenvalues for the quartic potential 80 3.3 Comparison of the eigenvalues for large e for the quartic potential 81 vi 3.4 Variation of Ai versus e for the quartic potential 82 3.5 Convergence of the eigenvalues for the optical bistability model. ..... 88 3.6 Comparison of the rate of the convergence of the eigenvalues for the optical bistability model 90 3.7 Convergence of the eigenvalues for the climate model 92 4.1 Convergence of the eigenvalues with V(y) = y2 + J^I, w\{y) exp(—aV2)- 114 4.2 Convergence of the eigenvalues with V(y) = y2 + 2, w\{y) = exp( —ay2). 115 4.3 Convergence of the eigenvalues with V(y) = y2+ 1+\VQ0 2, u>i(y) = exp(—ay2). 116 4.4 Convergence of the eigenvalues with V(y) = y2 + T+^2> 9 = 1) wi(v) = exP(-y2(l + T^Y/2) , • 117 4.5 Convergence of the eigenvalues with V(y) = y2 + , g = 10, w\(y) = exP(-y2(l + TT^)1/2) 118 4.6 Convergence of the eigenvalues with V(y) = y2 + 1^GY2, 9 = 100, vh{y) = exp(-y2(l + T^-)'/2) ' . .' 119 4.7 Convergence of the eigenvalues with V(y) = y2 + T+^pi 9 = 1) wi{u) — exV(-aiy2)/(l+gy2)a* 120 4.8 Convergence of the eigenvalues with V(y) = y2 + j^~2, 9 = 10, w\(y) = exV(-aW2)/(l+gy2)^ 121 4.9 The convergence of eigenvalues with V(y) = y2 + ; 9 = 100, wi(y) = exp(-aiy2)/(l + ^2)^ 122 4.10 (ai, a2) used for Tables 4.7-4.9 123 4.11 Comparison of the results of Ai with V(y) = y2 + 124 4.12 Convergence of the eigenvalues for SE with V(y) = y6 — 3y2 130 4.13 Convergence of the eigenvalues of even parity with V(y) = \y2 + 2y4 + |y6.132 4.14 Convergence of the eigenvalues with V(y) = y2 + ey4 134 vii 4.15 Convergence of the eigenvalues with V(y) = y2 + ey4 calculated by fitting weight function to ground state eigenfunction 135 4.16 Convergence of the eigenvalues, An^, for the Henon-Heiles potential with Hermite points 137 4.17 Convergence of the eigenvalues, An>;, for the Henon-Heiles potential with new points 139 4.18 Eigenvalues for the Henon-Heiles potential 141 vm List of Figures 2.1 Maximum error E°° in variation with N for the solution of u" = sin(7rx). 19 2.2 Eigenvalues of the 1st derivative matrix with Chebyshev discretization . . 21 2.3 Eigenvalues of the 1st derivative matrix with Legendre discretization. . . 22 2.4 Largest eigenvalues of the 1st derivative matrix with Chebyshev and Leg endre discretization 23 2.5 Largest eigenvalues of the 2nd derivative matrix with Chebyshev and Leg endre discretization 5 2.6 Largest eigenvalues of the 1st and 2nd derivative matrices calculated by the KoslofF modified method 27 2.7 Largest eigenvalues of the 1st and 2nd derivative matrices calculated by the Hermite method 9 2.8 Errors in an and Bn calculated by the Stieltjes's procedure 37 2.9 Minimum spacing of points versus number of quadrature points N. . . . 41 2.10 Numerical approximation of the eigenvalues for the analytic FPE problem. 51 2.11 Largest numerical eigenvalues for the analytic FPE example 54 2.12 Numerical solution of the time dependent FPE with N = 6 56 2.13 Numerical solution of the time dependent FPE with TV = 15 57 2.14 Exact solution for example B at z — 0 64 2.15 Standard error for the solution of example B 65 3.1 The first ten eigenfunctions for the quartic potential model 76 ix 3.2 The variation of the lowest nonzero eigenvalue, Ai versus e for the quartic potential 83 3.3 The potential U{x) in the equilibrium distribution for the optical bistabil-ity problem. 6 3.4 The first ten eigenfunctions of the Fokker-Planck operator for the optical bistability problem 87 3.5 The potential and the equilibrium distribution for the climate model ... 93 3.6 The first ten eigenfunctions of the Fokker-Planck operator for the climate model 95 3.7 The variation of the lowest nonzero eigenvalue, Ai for the climate model. 96 3.8 Time dependence of the probability density function for the climate model 98 3.9 The variation of the global average temperature 99 3.10 Eigenfunctions of the Fokker-Planck operator for electron relaxation in Xe. 100 3.11 Variation of the maximum eigenvalue versus the grid size. 101 4.1 The nonpolynomial oscillator (NPO) potential 112 4.2 Variation of the error in Ai for the NPO potential 126 4.3 Ground state eigenfunction for the NPO potential with g = \ — 100, N = 25 127 4.4 Ground state eigenfunction for the NPO potential 128 4.5 Contour plots of the eigenfunctions of the Schrodinger equation for the Henon-FIeiles potential 13x Acknowledgement I wish to express my sincere thanks to my research supervisor, Dr. B. Shizgal. It has been a pleasure to have worked with him, and his support, direction and assistance will always be greatly appreciated. This thesis is dedicated to my families and Michael. xi Chapter 1 Introduction Differential equations are widely used in almost all areas of science, engineering for modeling and forecasting. However, most of them can hardly be solved in closed form. Therefore, accurate and efficient numerical algorithms are required to provide solutions to many different types of differential equations. The major numerical techniques widely used include finite difference method (FD), finite element method (FEM), and spectral methods. The main objective of these numerical methods is to approximate derivatives by algebraic expressions involving the solution evaluated on a grid. Ordinary differential equations (ODE) are thus reduced to algebraic equations and partial differential equations (PDE) (PDE) are reduced to systems of, ODEs. In this section, we shall confine our discussion to linear differential equations. A differential equation can be written as where u(x) is the solution of the equation, and L is the linear differential operator. Discretizing the variable £ on a grid with N points, Eq. (1.0.1) can be approximated by algebraic equations with LN, the matrix representative of the differential operator L, and UJV, the approxi mate solution at the grid points. With different numerical techniques, the choice of the 1 Lu = -f{x) Chapter 1. Introduction 2 gridding and construction of the derivative and differential matrix may vary. Similarly, a time-dependent partial differential equation can be written in general as du(x,t)/dt = Lu + f(x), (1.0.3) where t is the temporal variable, x is the spatial variable, f(x) is a source function, u(x, t) is the solution and L is the differential operator which contains all the spatial derivatives of u. The steady solution of Eq. (1.0.3) is Eq. (1.0.1). Eq. (1.0.3) must be coupled with an initial condition u(x,0) and suitable boundary conditions. For simplicity, we only consider one spatial dimension in this chapter. Similar statements can be applied to the higher dimensions. By discretizing the spatial variable x, the partial differential equation, Eq. (1.0.3) can be approximated by a system of first order ODEs, -^- = LN + fN. (1.0.4) We shall often confine the discussion of time discretizations to Eq. (1.0.3) with f(x) = 0, i.e., we consider -^ = LNuN. (1.0.5) There are several ways to solve the ODE system Eq. (1.0.5). The solution of this time dependent PDE can be evaluated by an eigenfunction expansion. The formal solution (with f(x) = 0) is then of the form u(x,t) = ^cneA"Vn(x), (1.0.6) n=0 where An and <f>n(x) are the eigenvalues and eigenfunctions of the eigenvalue problem L(f)n(x) - \n<f>n(x), (1.0.7) the coefficients cn are determined from the initial condition, u(x,0), and are given by cn = J u(x, 0)(j)n(x)dx. (1.0.8) Chapter 1. Introduction 3 The convergence and accuracy of the eigenvalues and eigenfunctions of the matrix Lyv determine the convergence and accuracy of the solution. The convergence as well as the rate of the convergence of the eigenvalues and eigenfunctions depend on the choice of the numerical methods and spatial gridding used to construct the differential discretization matrix representation in the spatial variable. The second method is to discretize the ODE system Eq. (1.0.5) in time, and thus integrate the ODE system numerically. This approach is applicable for both linear and nonlinear problems. In most applications, the temporal discretization uses conventional finite difference methods, including implicit, explicit and semi-implicit or implicit-explicit schemes. Among the factors which influence the choice of a time-discretization are the accuracy, stability, storage requirements and work demands of the methods. The explicit scheme is easy to implement and relatively less costly, but has a strong stability restriction which is determined by the eigenvalue spectrum of the differential operator L. The popular explicit schemes are the Adams-Bashforth methods which includes the simple first order forward Euler (FE) method. The implicit scheme provides the best stability but involves a matrix inverse and is expensive, especially for nonlinear problems. It is the best choice for stiff (large condition number) linear problems which encounter instability while using an explicit scheme. A related set of implicit methods are the Adams-Moulton methods which include the first order backward Euler (BE) and the second order Crank-Nicolson (CN) method. Implicit-explicit schemes become popular for the time integration of nonlinear problems. The idea of these type of schemes is to an apply implicit scheme for the linear term and an explicit scheme for the nonlinear term of the problem. Another popular class of temporal discretization methods are the Runge-Kutta Methods. Some standard references in these temporal discretization methods are the books by Gear (1971) [1], Lambert (1973) [2], Shampine and Gordon (1976) [3], Dekker (1984) [4], Butcher (1987) [5], and Hairer (1989) [6]. Chapter 1. Introduction 4 The basic idea above is to discretize the differential operator L in physical space, i.e. approximate the differential operator directly by discretization. Most finite difference and collocation methods use this idea. Another popular approach is to approximate the solution in transform space. It involves expansion of the solution by a finite sum. For example, the solution u(x) of Eq. (1.0.1) can be expanded in a set of orthogonal basis functions ^fc(x) (k = 1,2,...,^), that is k=N u(x) = J2 CkMx)- (1-0-9) k=i The transform coefficients are constant and can be solved by an algebraic equation Ac =b (1.0.10) where c is the coefficient vector of c/;, and A and b are the transform matrix and vector determined by minimizing the residual function R(x; Ck) = Lu + / or its related function-als. The methods used in this approach include finite element methods, ta,u and Galerkin spectral methods. Finite difference methods for solving differential equations consist of replacing each of the derivatives in the differential operator by an appropriate difference-quotient approx imation (such as du(x)/dx [u(x + h) — u(x — h)]/2h, where h is a small grid spacing). The difference quotient is generally chosen such that a certain order of truncation error is maintained. The method approximates derivatives of a function by local arguments and is designed to be exact for polynomials of low order. The matrix representative of the derivative operator is generally banded. The FD algorithm is easy to program and implement, but relatively low in accuracy. A higher order of accuracy for the finite difference method is possible and has strong restriction to the stability of the solution. The FEM is frequently used to solve partial differential equations that occur in en gineering applications, especially, in the areas of solid mechanics, elasticity, and compu tational fluid dynamics [7-17]. The idea of the finite element method is to divide the Chapter 1. Introduction 5 interval in x into a number of sub-interval and approximate the solution by a sum of basis functions which are piecewise polynomials of low fixed degree. The differential ma trix derived is sparse because only a handful of basis functions are non zero in a given sub-interval. Since the approximate solutions are piecewise polynomials, the pieces or the sub-intervals can be easily chosen to fit the geometry of the problem. The FEM is therefore very useful in solving problems with complex geometry in multi-dimensional problems. The FEM is closely related to the finite difference method. The discretization of the differential equation are difference equations, but derived by a different approach. In some cases (especially in one-dimensional cases), they can be equivalent. The disvan-tage of the FEM is the relatively low accuracy because each basis function is a polynomial of low degree. Spectral methods involve the expansion of the solution by a finite sum N u(x) « uN{x) = ak4>k{x) (1.0.11) k=0 in a set of orthogonal basis functions (trial functions), (j>k{x). When this series is substi tuted into Eq. (1.0.1), we can define a residual function R(x] a0, ai, ...,a^) = LUN + /. (1.0.12) The spectral method is constructed using test functions to minimize the residual i.e. the coefficients in the series expansion are chosen so that the residual is orthogonal to the space spanned by the basis of test functions, tpk{x): J R(x] ao, ai,ajs[)il>k{x)dx = 0. (1.0.13) This is analogous to techniques used in the FEM. The choice of the test functions dis tinguishes between the spectral Galerkin, tau and collocation (also referred to as pseu-dospectral) methods [18,19]. In the Galerkin approach, the test functions are the same as Chapter 1. Introduction 6 the trial functions and individually satisfy the boundary condition. Spectral tau meth ods are similar to Galerkin methods, but none of the test functions need to satisfy the boundary conditions. Hence a supplementary set of equations is required to apply the boundary condition. In the collocation method, the test functions are shifted Dirac delta functions centered at so-called collocation points. For many problems, especially nonlin ear ones, the collocation ( or pseudospectral) method is the easiest to implement, and the most efficient [19-23]. Spectral methods use global basis functions which are polynomials (or trigonometric polynomials) of high degree and generate solutions with high accuracy. Pseudospectral method can be viewed as the limit of FD methods when the order of accuracy tends to infinity [24]. Unlike low order FD method and FEM which result in large, sparse matrices; the differential matrix representative by the spectral method is usually full. There are numerous numerical schemes available in the numerical computation. Most of them belong to the three methods discussed above, or schemes which are a combination of several techniques. The geometry of the domain and required level of accuracy are the key factors in selecting among these approaches. The FEM is particularly well suited to problems in very complex geometries, whereas spectral methods can offer superior accu racies and cost efficiencies, mainly in simple geometries. Although sometimes a mapping of a complicated physical domain onto a simple computational domain is possible, the mapping must be smooth to preserve the high-order accuracy and exponential conver gence rates associated with the spectral method. Finite difference methods perform well over a broad range of accuracy requirements and (moderately complex) domains. For years, many numerical analysts have been working on the spectral domain decomposi tion or spectral element methods [25-32]. The techniques can apply to problems with complicated domain and provide superior accuracy and convergence over the FEM. The advantage of spectral methods is that they are more accurate than the FD method Chapter 1. Introduction 7 and the FEM so that far fewer grid points are required. A principal disadvantage is that in solving time dependent problems in finite boundaries, they are often subject to tight stability restrictions when an explicit temporal discretization methods is used. For a specific problem, a particular numerical scheme may be selected to meet stabil ity and convergence requirements, as well as the physical and geometric aspects of the problem. In the context of solving time-dependent PDEs, spectral methods have notable ad vantages over other numerical methods. First, the error between the numerical solution and exact solution decays exponentially versus number of grid points N. This high accuracy allows one to treat problems with very much fewer grid points which would otherwise require an enormous number of mesh points by other numerical methods. For multi-dimensional problems, the use of relatively coarse grids that suffice for most accu racy requirements offer superior efficiency in time and memory. Besides this potentially high accuracy, spectral methods are also powerful for many cases in which solutions and variable coefficients are nonsmooth or even discontinuous [19,22,33,34]. Spectral methods and pseudospectral methods based on the expansion of the solu tion in orthogonal polynomials have been successfully applied in many areas such as simulations of turbulence flow [18,19,35], computational fluid dynamics [18-20,29,36-41], weather prediction [19,42,43], kinetic theory [44-47], and quantum mechanics [48-51], etc. It has been demonstrated numerically that the condition number, i.e. the ratio of the largest to the smallest eigenvalue of the spectral approximation matrix is large. It is known that finite difference and finite element methods have an 0(N2) condition num ber. For spectral methods, the largest eigenvalue of the first derivative operator grows as 0(N2), and of the second derivative operator grows as 0(N4), as N increases, leading to Chapter 1. Introduction 8 poorly conditioned matrices [18,19,29,36,39,52-58]. This effect upon the numerical com putations of the eigenvalues of the derivative operator is caused by the round off error depending on machine precision of the computer. One cause of the roundoff error in us ing the spectral differentiation matrix is attributed to inaccurate computation of certain matrix elements, particularly certain large elements in the upper and lower right corner of the matrix [52,53,56]. Methods for reducing the error of Chebyshev pseudospectral approximation of derivative matrix are discussed by several researchers [52,53,55,56,59]. Several researchers [29,60] also point out that the poor condition of the derivative ma trices is related to the minimum spacing of the grid points. The large condition number of the spectral derivative matrix is not necessarily a problem in the accuracy of the so lution of differential equations. However, the large condition number of the derivative matrix brings stringent stability restrictions to solve time dependent PDEs by temporal discretization, since in actual computations, the stability limit depends on the precision of the numerical solution of eigenvalues. On the other hand, the differential matrices by a spectral method are usually dense. Direct inversion of these matrices is usually expen sive for large matrices. Iterative schemes are sometimes a practical necessity. Therefore, reduction of the condition number of spectral differential matrices is necessary for stabil ity of iterative methods and time integration. There have been numerous studies on the preconditioning of spectral matrices. The idea of preconditioning with second-order finite difference methods is described by Orszag [25]. Some other preconditioners based on the FD [61,62], the FEM [31,63,64], and the spectral methods [65,66] have been proposed and successfully tested in the literature. Besides preconditioning techniques, algorithms for direct improvement of the condition number of derivative matrices are also discussed by Heinrichs [57], and Kosloff and Tal-Ezer [60]. Most spectral and pseudospectral methods currently used are restricted to traditional or "classical" basis functions such as trigonometric functions and polynomial solutions Chapter 1. Introduction 9 of Sturm-Liouville eigenvalue problems. The details of these polynomials are described in most texts on spectral methods [18,19,29,36,38,39]. The choice of basis functions is dictated by the weight function, w(x) and the domain of the problem considered. For example, for [—1,1] and w(x) = 1, w(x) = (1 — x)a(l + x)P and w(x) = ^/=p", one has the "classical" Legendre, Jacobi, and Chebyshev polynomials, respectively. Laguerre and Hermite polynomials are used for w(x) = e~x on [0, co] and for w(x) = e~*2 on [-co, oo], respectively. When the basis set is chosen, the associated grid points in a pseudospectral method are fixed, and adaptive gridding techniques are difficult but possible [40,67]. In many cases, the choice of grid points is essential to the convergence of the eigenvalues and eigenfunctions of the differentiation matrix. For some physical problems, spectral solution with classical polynomials may result in poor convergence or even inaccurate results. Therefore more flexibility in choice of grid points is desired. Many problems of physics and engineering lead naturally to the solution of a partial differential equations in an unbounded domain. The condition at infinity given by a certain asymptotic behavior of the solution can be obtained by asymptotic analysis at infinity. However, for numerical approximation, it is not easy to interpret the behavior at infinity. Most numerical methods such as finite difference, finite element methods and spectral methods with trigonometric functions and Jacobi-type polynomials as basis functions can only apply to the problem with finite domain. One of the most widely used techniques is to restrict the computation to a finite domain by truncating the domain and imposing "artificial boundary" conditions. Another treatment consists of mapping the unbounded domain into a bounded one and using the standard discretization techniques mentioned above. The more attractive technique is to use basis functions defined on an unbounded domain such as Hermite polynomials (or Gaussian-type polynomials) on (—co,oo) and Laguerre polynomials on (0,oo). For those polynomials, the collocation points spread all over the infinite domain with increasing N, and there is no restriction Chapter 1. Introduction 10 on the size of the domain is required and thus no artificial boundary cutoff is needed. As to the approximation by spectral methods, some results and a comparison with these techniques can be found in the papers by several workers [29,67-69]. Although use of Hermite and Laguerre polynomials bring some promising results for solving unbounded boundary problems, the direct spectral approach may not produce good approximation. Scaling and shifting factors may be necessary in practical applications of the classical spectral methods. On the other hand, many physics problem have some special charac teristics and need special treatment when solving these models. For example, the solution of the optical bistability model for laser physics is usually characterized by two stable states. Use of the Hermite spectral method may result in slow convergence and/or poor accuracy. Thus, it is necessary to find some nonclassical basis functions that can relate to the physics of the problem very well and generate good convergence and accuracy of the solution. Recently, Shizgal and coworkers started using nonclassical polynomials on an infinite domain to solve some physics related problems [43,45-47,51,70-72]. They received very encouraging results. Another reason for using those physics related basis functions is that they have very attractive properties from the numerical point of view. There is evidence that the condition number of the first and second Hermite differen tiation matrices are 0(\/rN) and O(N), respectively [29,68,73]. Numerical results from our work also showed that the condition number for a second differentiation matrix can be improved to be between O(N) and 0(N2) depending on the physical problem and polynomial basis set chosen. This is better than second order finite difference which is 0(N2), and far better than 0(N4), for example, by Chebyshev spectral method. For these reasons, it is necessary to have a spectral method which is not restricted to the "classical" orthogonal basis functions and associated points. The purpose of this thesis is to study a fast, and accurate spectral method, namely, the quadrature discretization method (QDM), and apply the method to solve several Chapter 1. Introduction 11 differential equations, notably, the Fokker-Planck equation (FPE) and the Schrodinger equations (SE). nonclassical polynomials as basis functions in a spectral approximation. The idea of the QDM is initially introduced by Shizgal in solving eigenvalues of Lorentz Fokker-Planck equation [74]. A set of nonclassical polynomials called "speed" polynomials based on weight function x2 exp(—x2) was used in this work. The method was later developed fur ther by Shizgal and Blackmoreby using speed and some other weight functions [71,72,75]. Their generation of polynomials was based on a recursion subject to round-off errors and hence their approach was limited to a small set of weight functions. In this thesis, by introducing the Stieltjie's procedure, it is easy to construct accurately a set of nonclassi cal polynomials defined by any weight function. For a specific problem, a specific weight function may be chosen to optimize the accuracy and convergence of the solution. The details of the QDM is described in Chapter 2. Another main interest of this thesis is to solve differential equations in an unbounded domain. Two types of second order differential equations in an unbounded domain are discussed and solved numerically by the QDM. The first type of equation is the Fokker-Planck equation (FPE) of the form: with homogeneous boundary conditions at infinity. P(x,t) is related to some probability density function of a system. A(x) and B(x) are referred to as the drift and diffusion coefficients and depend on the particular application considered. The FPE was introduced by Fokker and Planck to describe the Brownian motion of particles [76,77]. For the past several decades, there has been an ongoing interest by numerous researchers in the description of nonequilibrium phenomena modeled with a In this thesis we introduce the quadrature discretization method (QDM), which allows Chapter 1. Introduction 12 FPE. This interest continues unabated to the present date [78]- [91]. The basis for many of these models is Brownian diffusion in a potential characterized by Gaussian white noises. This leads to a time dependent linear FPE with drift and diffusion coefficients which can be nonlinear functions of the independent variable of interest. The theoretical basis for this approach has been provided in several standard references [92-94]. The traditional method for the solution of the FPE is a spectral method which usu ally involves the expansion of the probability density function (PDF) in a suitable basis set, and the reduction of the differential equation to a set of algebraic equations for the expansion coefficients. An alternate approach involves the discretization of the PDF on a grid of points. This discrete approach in the solution of differential and/or inte gral equations has been used by researchers in other fields, notably neutron transport [95], radiative transfer [96], and computational fluid dynamics [18,36]. Fourier series or Chebyshev polynomials are almost exclusively chosen as basis functions in the applica tion of the pseudospectral approach. Other popular discretization schemes are based on the finite-difference technique such as those proposed by Chang and Cooper [97], Larsen et al [98] and Epperlein [99] primarily for the solution of a nonlinear FPE that arises in plasma physics. Park and Petrosian [100] have recently provided a detailed comparison of several different methods to the solution of FPEs applied to astrophysical problems. Another type of equation to be studied in this thesis is the Schrodinger equation (SE) that arises in quantum mechanics and molecular physics. The equation has the form with V(y) as the potential. There have been numerous studies on the solution of the SE with different methods for several potential function V(y) [101]- [133]. The methods in clude variational and spectral techniques with Gaussian-type basis functions [29,101-111], perturbative approach [112-117], Pade approximation [118-122], FD methods [123-125], Chapter 1. Introduction 13 and some other methods [126-131]. A very similar approach to the QDM for the solution of the SE referred to as the Discrete Variable Representation (DVR) was developed by Light and co-workers [132], and has been used by several other researchers [133]. The QDM and DVR differ in the use of polynomial basis functions. The QDM uses nonclas-sical polynomials and the DRV uses classical polynomial, majorly the Hermite (or scaled Hermite) polynomials. Hermite functions are the exact eigenfunctions of the quantum mechanical harmonic oscillator which has one well in the potential. This close connec tion with the physics makes them a nature choice of basis functions for the quantum mechanical problems. However for many problems with more than one well in potential Hermite polynomials are barely a good approximation of the solutions, as we will show later in this thesis. Therefore, new polynomial sets suitable to the physical nature of the solution are needed to improve the accuracy and convergence. The purpose of this thesis is to study the quadrature discretization method and apply it to solve second order partial differential equations with unbounded boundaries. It is shown in this thesis that the QDM is easy to implement. The method provides high accuracy and rapid convergence, and is very efficient in solving high dimensional PDEs. Unlike most spectral methods and FD methods, it is possible for the QDM discretization to preserve the symmetry for the self-joint joint differential operator. Furthermore, the flexibility of the QDM in choosing weight function provides one the opportunity to achieve the best solution with the least computation work. In Chapter 2, the basic ideas of the QDM are introduced. The QDM is tested and compared with other numerical methods. The applications of the QDM to the FPEs and SEs as well as the implementations and the results are presented in Chapter 3 and Chapter 4, respectively. Finally, a summary of the completed research and suggestions for the future work are given in Chapter 5. Chapter 2 The Quadrature Discretization Method (QDM) The QDM, originally introduced by Shizgal et al as discrete ordinate method [71, 72], is a numerical method for the efficient solution of integral and differential equations. The QDM can be used to solve a large class of differential and integral equations [71,72] and has also been applied to fluid dynamics problems [40]. The QDM is a spectral method based on the discretization of the solution on a grid of points (collocation points) that coincide with the points of a quadrature. The quadrature is based on a set of nonclassical polynomials orthogonal with respect to a specified weight function. In this chapter, the fundamentals of the QDM are introduced and discussed. It differs from traditional spectral methods in that it is based on nonclassical polynomials. In the next section, several aspects of classical polynomials and spectral/pseudospectral methods are discussed. The numerical solution of differential equations is based on the approximate representation of the first derivative operator by a finite derivative matrix of dimension equal to the number of grid points. Higher order derivatives are represented by powers of the derivative matrix. The QDM is based on the use of nonclassical polynomials which are generated by a method developed by Gautschi [134]. Several model problems are studied with the QDM, and the results are compared with those for the classical pseudospectral method and the finite difference method. In the comparison of exact and approximate solutions, we introduce the errors E\ (average error), E2 (standard error) and E°° (maximum error) defined by 14 Chapter 2. The Quadrature Discretization Method (QDM) 15 N (2.0.1) 1 N (2.0.2) and E°° = max{\u1 — u\ ; = 1,2, ...,#}, (2.0.3) where uf is the approximate solution and u\ is the exact solution. 2.1 Classical polynomials and the pseudospectral method Classical polynomials referred to in this thesis are the polynomials used by most re searchers in the application of spectral methods. They include Chebyshev, Legendre, Jacobi, Laguerre and Hermite polynomials as well as their related polynomials. The details of these polynomials are discussed in most standard references on orthogonal polynomials and spectral methods [18,39,36,135]. Here we only give a brief summary of these polynomials and some of their properties. Classical polynomials used in spectral methods are related to the polynomial solutions of eigenvalue problems of Sturm-Liouville equation in some interval /. An are the eigenvalues and un{x) are the corresponding eigenfunctions. The coefficients a(x), h{x) and w(x) are real-valued functions, and w(x) is referred to as the weight function. With particular choices for these coefficients, several classical polynomial sets are obtained. For Jacobi polynomials, a(x) = (1 - x)a+1(l + xY+\ b{x) = 0, w{x) = (l-x)a(l + x)^ (2.1.2) Chapter 2. The Quadrature Discretization Method (QDM) 16 for Va; £ (—1,1), and a > — 1, ft > —I. Legendre and Chebyshev polynomials are Jacobi polynomials with a = (5 = 0 and a = ft = — \ respectively. For Laguerre polynomials, we have the coefficients a(x) = xa+1e~x i b(x) = 0, w(x) = xae~\ (2.1.3) for \/x G (0, oo), a > — 1, whereas for Hermite polynomials, these are given by, a(x) = w(x) = e~x\ b(x) = 0, (2.1.4) for 6 ( — 00,00). The polynomials satisfy a three recurrence relation [39] Un+1(x) = (pnX + On)un(x) + TnUn--L, (2.1.5) where u0(x) and u\(x) are the first two polynomials in a given set. The recurrence relations for Legendre, Chebyshev and Hermite polynomials are as follows: Legendre 2n + 1 n Ln+1(x) = ———xLn(x) —Ln_i(x), (2.1.6) Tl -\- 1 77. —I— _L with Lo(x) = 1, Li(x) = x; Chebyshev Tn+1{x) = 2xTn{x)-Tn^(x), (2.1.7) with To(x) = 1, T\(x) = x; Hermite Hn+1(x) = 2xHn(x)-2nHn-1, (2.1.8) with H0(x) = 1, H^x) = 2x; If the zeros of the Nth order polynomial are {xi,i = 1,2,.., N}, a Gauss quadrature formula, for an arbitrary smooth function f(x), gives N / f(x)w(x)dx m J2f(xj)w3, (2.1.9) JI 3=1 Chapter 2. The Quadrature Discretization Method (QDM) 17 where Wj = JIlj(x)w(x)dx, are the weights, and Xj are the quadrature points of the integration formula. If f(x) is a polynomial of degree 27V — 1 or less, Eq. (2.1.9) is exact [39,136]. The function lj(x) is the interpolating polynomial uniquely defined by the condition ro, iii^j lj(xi) = { (2.1.10) ll, iii = j The quadrature points and weights for the classical polynomials and related Gauss-type quadrature formula can be found in standard references in spectral methods [18,19]. The classical pseudospectral method is based on the discretization of a smooth func tion f(x) on a grid of these quadrature points {xi, i = 0,1, 2,N}. Derivatives of the function are approximated by the analytic derivatives of the interpolating polynomial lj(x). The first derivative can be represented in a matrix form as N (DNf)(Xi) = £(Av),?/(*-.;), (i = 0,1,2,AO, . (2.1.11) 3=0 where DN is the derivative matrix with the entries (Dw)ij = l'j(xi). Higher order deriva tive matrices can be obtained by multiplying DN by itself a suitable number of times. The spectral derivative matrix is usually full. The explicit expression of the derivative matrices based on the quadrature points is described in the references by Canuto [18], Peyret [137,138], and Funaro [39,68]. The advantage of the spectral method is its spectral accuracy, that is, errors for numerical approximation decrease exponentially with increasing number of grid points. Consider a second order differential equation d2u . , N . . — = sm(7rx) (2.1.12) with homogeneous Dirichlet boundary condition u( — 1) = u(l) = 0. The problem has an analytic solution u(x) = — ir2sin(7rai). We solve the problem with a Legendre spectral Chapter 2. The Quadrature Discretization Method (QDM) 18 Legendre FD N Ex E2 E°° Ex E2 E°° 4 4.09(-02) 5.01(-02) 6.13(-02) 1.56(-01) 1.91(-01) 2.34(-01) 8 2.40(-04) 2.83(-04) 4.46(-04) 3.66(-02) 4.01(-02) 5.30(-02) 16 1.50(-11) 1.72(-11) 2.91(-11) 8.68(-03) 9.46(-03) 1.30(-02) 32 2.73(-14) 3.23(-14) 6.03(-14) 2.11(-03) 2.31(-03) 3.22(-03) 64 1.21(-14) 1.52(-14) 4.36(-14) 5.19(-04) 5.73(-04) 8.04(-04) 128 1.29(-04) 1.43(-04) 2.01(-04) 256 5.01(-05) 512 1.26(-05) Table 2.1: Errors for the numerical approximation of u" = sin(7rx) with a Legendre spectral method and a finite difference method. method and a 2nd order central difference method. Table 2.1 gives the errors of the numerical solutions of Eq. (2.1.12) with both methods. It can be seen from the results that for the Legendre method, the error decays very fast. The maximum error decays to O(10~14) within 32 points. Increasing TV beyond this point would not reduce the errors significantly, and the errors remain close to machine accuracy. For the FD method, the error decays at a much slower rate. With 128 grid points, the error is O(10-4) and with 512 grid points, it is O(10-5). This indicates that to obtain 4 decimal places accuracy for the solution of Eq. (2.1.12), the FD method requires more than 100 grid points, while for the Legendre spectral method, approximately 8 points are needed. As a matter of fact, the error for the FD method decays at the rate of 0(N~2). For the Legendre spectral method the errors decay at an exponential rate as the number of grid points N increases. To show this point, we plot the logarithm of the errors versus logio(N) in Fig. 2.1. The exponential accuracy of the spectral method makes it very attractive for solving multi-dimensional problems since much coarser grids are needed in comparison with the finite difference method. Chapter 2. The Quadrature Discretization Method (QDM) 19 Chapter 2. The Quadrature Discretization Method (QDM) 20 As mentioned in Chapter 1, spectral derivative matrices based on classical polyno mials have some drawbacks. The main disadvantage of the spectral discretization is the large value of the largest eigenvalue of the differentiation matrix, which can result in an instability when iteration and/or explicit time integration methods are applied. The largest eigenvalue of the derivative matrices with Chebyshev and Legendre collocation methods grows very fast as number of grid points N increases. Consider the eigenvalue problem of the advection operator ^j^ = X»Mx) on (-1,1) (2.1.13) subject to the boundary condition <M1) = 0 (2.1.14) The eigenvalues are computed with Chebyshev and Legendre collocation methods. The boundary condition is imposed by removing the last row and column of the derivative matrix DM- The resulting (A^ —1) X (N — 1) is diagonalized numerically. Figs. 2.2 and 2.3 show the eigenvalues of the first derivative matrix with N—8, 16, 32, 64. Although most of the eigenvalues converge to a curve in the left half-plane between —iN and iN, a few spurious eigenvalues diverge at a much faster rate. The modulus of the largest eigenvalue |AJV| m relation with the number of grid points is plotted in Fig. 2.4. The largest eigenvalue grows like 0(N2) as N —> oo. For the Le gendre method, |A/v| ~ 0.0S0N2, and for the Chebyshev method, \XN\ w 0.089N2. For the second derivative matrix, we consider the eigenvalue problem of the diffusion operator ^M = Xn(f)n(x) on (-1,1) (2.1.15) with the boundary conditions of Dirichlet type: </>n(l) = </>„( —1) = 0. The problem has the analytical solution, XN = — (nir/2)2 and (j)n(x) = siny/X^x, n = 0,1,2,.... Chapter 2. The Quadrature Discretization Method (QDM) 21 s> oh 16 32 ro 128 64 h Oh -64 -128 I 1 1 1 1 o 1 1 N=32--o -- o -1 0 1,1, 1 512 256 -128 -64 0 Real 64 128 -256 h -512 -512 -256 256 512 Figure 2.2: Eigenvalues of the 1st derivative matrix with Chebyshev discretization with N=8, 16, 32, 64. Chapter 2. The Quadrature Discretization Method (QDM) 22 Figure 2.3: Eigenvalues of the 1st derivative matrix with Legendre discretization with N=8, 16, 32, 64. Chapter 2. The Quadrature Discretization Method (QDM) 23 4.00 0.00 1 ' 1 ' 1 ' 1 ' 1 ' 1 0.40 0.80 1.20 1.60 2.00 2.40 log 10 (N) Figure 2.4: Eigenvalues of 1st derivative matrix with Chebyshev and Legendre discretiza tion. The solid line is for Chebyshev. The dashed line is for Legendre. Chapter 2. The Quadrature Discretization Method (QDM) 24 The largest eigenvalues calculated with Chebyshev and Legendre collocation methods are given in Fig. 2.5. The boundary conditions are imposed by removing the first and last rows and columns of the derivative matrix D^. The resulting (TV — 2) x (N — 2) is diagonalized numerically. The largest eigenvalue grows like 0(N4) as N —» co, and |An| = 0.026Ar4 and |An| = 0.047iV4 for the Legendre and Chebyshev method, respectively. Several researchers showed that the effect upon the numerical computation of the eigenvalues of the derivative matrix is caused by round off error [52,53,55,56,59]. This property of the derivative matrix does not play a role in the solution of a differential equation. As seen in Fig. 2.1, the error of the solution of the differential equation Eq. (2.1.12) decreases at a very fast rate with increasing. However, for the time-dependent differential equation, the large value of the largest eigenvalue causes an instability in the time integration when explicit time-stepping techniques are applied. For the first order explicit Euler method, for example, the time step At to ensure the stability of the time integration may be 0(N~4) for a 2nd order time-dependent differential equation. To reduce the order of growth rate of the largest eigenvalue, preconditioning of the derivative matrix is usually used. The topic of the preconditioning is out of the scope of this thesis. Methods for preconditioning can be found in many references [66,139-141]. Another way to improve the condition of the spectral derivative matrix is to use mapping techniques. Several researchers suggested that [29,60] the growth rate of the largest eigenvalue is related to the minimum spacing Asmm between the adjacent grid points. The largest eigenvalue (in modulus) of the first derivative matrix is of the order of 1/Axmin, and for the second derivative matrix it is of the order of 1/A2xmin. Most classical points in finite domains, for example, Chebyshev and Legendre points, have minimum spacing of 0(1/N2) such that the 1st and 2nd derivative matrices grow like 0(N2) and 0(N4), respectively. Kosloff and Tal-Ezer proposed a modified Chebyshev method to improve this condition [60]. They mapped the Chebyshev points in [-1,1] Chapter 2. The Quadrature Discretization Method (QDM) 25 Figure 2.5: Eigenvalues of the 2nd derivative matrix with Chebyshev and Legendre dis cretization. The solid line is for Chebyshev. The dashed line is for Legendre. Chapter 2. The Quadrature Discretization Method (QDM) 26 to a new set of grid points in [-1,1] which has minimum spacing of 0(1/N) such that the largest eigenvalues of 1st and 2nd derivative matrices improve to 0(N) and 0(N2), respectively. Their algorithm is similar to those of the Fourier method but it provides accurate solutions for nonperiodic problems. The detailed description of the method can be found in their paper [60]. Fig. 2.6 give the largest eigenvalues related to the modified Chebyshev approximation versus the number of grid points N. Classical Chebyshev and Legendre polynomials are defined in a finite domain (-1,1). For a problem in an arbitrary finite domain (a,b), one can either map the domain into (-1,1) by using change of variable, or use scaled Chebyshev and Legendre polynomials. For problems with unbounded domain, one needs to truncate artificially the infinite domain into a finite domain in order to apply the Chebyshev or Legendre spectral method. The accuracy and convergence of the solution are directly related to the domain truncation. For a certain number of grid points, the larger the domain, the coarser the grid, and the slower the convergence. If the domain is too small, the error caused by domain truncation may exceed the accuracy requirement. The Hermite points for the infinite domain naturally spread over the whole domain with increasing ./V and no artificial boundary cutoff is involved. As shown in several literatures [29,68,73], Hermite points have another advantage that the minimum spacing is 0(l/y/~N) such that the largest eigenvalues of the first and second derivative matrices are 0(y/~N) and O(N) as N —> oo. This property of Hermite polynomials and their close connection to the physics usually make them the natural candidate for solving problems in an infinite domain. However, it may not be always the best choice as we will see later in the thesis. The largest eigenvalues of the derivative matrix for the Hermite points versus number of grid points N is plotted on a logarithm scale in Figure 2.7. Numerical fitting shows that the largest eigenvalue satisfies AJV ~ Q.2SN0^ for the first derivative matrix and XN ~ 0.3SN1'27 for the second derivative matrix. These orders are slightly Chapter 2. The Quadrature Discretization Method (QDM) 27 Figure 2.6: Largest eigenvalues of the 1st and 2nd derivative matrices calculated by the Kosloff modified Chebyshev method. Chapter 2. The Quadrature Discretization Method (QDM) 28 larger than the asymptotic results. 2.2 The quadrature discretization method 2.2.1 Construction of derivative matrices Consider a set of polynomials, Fn(x), orthogonal with respect to a weight function w(x), that is, J w(x)Fn(x)Fm(x)dx = 8nm. (2.2.1) There are as many quadrature rules as there are sets of orthogonal polynomials. The familiar quadratures are those based on the classical polynomials discussed in the last section. Theoretically, a set of orthogonal polynomials can be constructed for any interval and any weight function. The major development of the QDM is to generate, with the method introduced by Gautschi [134], new sets of points and weights from nonclassical polynomials chosen to maximize the rate of convergence of the solution of a specific differential equation. The QDM is based on the discrete representation of a function j\x) by its values at the set of N quadrature points Xi, that is, f(xi). We expand f(x) in the {Fn} orthogonal basis set, N f{x) = J2anFn(x), (2.2.2) n=l where an are the expansion coefficients of f(x) and are given by, an = jw(x)f(x)Fn(x)dx. (2.2.3) The expansion coefficients can be evaluated numerically with the quadrature rule in Eq.(2.1.9), and N an = Y, W3f(x,)Fn(xj). (2.2.4) Chapter 2. The Quadrature Discretization Method (QDM) 29 0.00 1 1 1 1 1 ' 1 0.40 0.80 1.20 1.60 2.00 log 10 (N) Figure 2.7: Largest eigenvalues of the 1st and 2nd derivative matrices calculated by the Hermite method. Chapter 2. The Quadrature Discretization Method (QDM) 30 Eq. (2.2.4) gives the transform from the physical space {/(x,)} to the transform space {an}. If we use Eq.(2.2.4) to evaluate f(x) at the quadrature points X{ we have another transform that does the reverse transformation which is from the transform space {an} to the physical space {/(a^)}, N f{xi) = '£anFn{xi). (2.2.5) 71=1 Substituting Eq. (2.2.4) into Eq. (2.2.2), we obtain f(x) = J2J:wlf(xl)Fn(x1)Fn(x). (2.2.6) 1=1i=i The above equation indicates that f(x) can be calculated at any point from its values f(xj), that is, N /W = E'iW/W. (2-2-7) where the interpolation polynomial lj(x) is given by N 71=1 and satisfies (0, ifi^j Ijixi) = <| (2.2.9) 1, ifi = j If Eq. (2.2.7) is differentiated, one obtains DF{XI) -EA./fe), (2.2.10) where A,-j are the entries of the derivative matrix A and are defined by = /;(.'.•,) = E '^./^(x,)/.;:!,-;). (2.2.11; 71 = 1 It is easy to notice that the transformations between the physical space and transform space shown in Eq. (2.2.4) and Eq. (2.2.5) are not symmetric. To symmetrize the transformations, we consider an alternative function / defined by /(•'•) = V™j(x). (2.2.12) Chapter 2. The Quadrature Discretization Method (QDM) 31 If the expansion coefficients of f(x) in the basis set Fn(x) are an, then there is a unitary (or orthogonal) transformation between the physical and transform spaces, that is, N f(X*) = E V^Fn{ (2.2.13) n=l N ^^E^-'N/W. (2.2.14) i=i where the matrix elements of the symmetric transformation T are T;n = s/w^n^Xi). If in Eq. (2.2.13), Xi is replaced with x, and Eq. (2.2.14) is used for an, one obtains a A^th-order interpolation, given by /» = E4(*)/(*;), (2-2-15) 3 = 1 where the interpolation polynomial Ij(x) is given by AT Ij{x) = y/w^J]^2Fn{xj)Fn(x). (2.2.16) 71=1 The algorithm for differentiation in the discrete basis is given by differentiating Eq. (2.2.15) and using Eq. (2.2.16) one has that df{ Xi N dx . . .7 = 1 - E (2-2.17) where the derivative operator, Dij = y/W{W3T'j(xi) is given by N Di3 = ^JT3 E Fn(x3)F^{xt). (2.2.18) n=l We refer to the derivative matrix D with the elements Dij as the QDM modified derivative matrix. We have the following relation between the derivative matrices A,-7 and Dij &n = vS'MAr (2.2.19) The matrix element of the derivative operator d/dx in the transform space (also called Galerkin 's matrix) is dnm = j w(x)Fn{x)F'm{x)dx. (2.2.20) Chapter 2. The Quadrature Discretization Method (QDM) 32 If we evaluate it with the quadrature rule in Eq.(2.1.9), we have that, N dnm = E WkFn(xk)F'n(xk). (2.2.21) For some specific physical models, using the transform matrix have an attractive advan tage that the resulting differential matrix in the transform space may be banded [29,142]. If we transform the derivative matrix d back to the physical space with transform matrix T, we obtain the QDM derivative matrix D = TdTT, that is. N N N Dv = EE Bn(Xl)^/u7z[J2 wkFn(xk)F^(xk)}Fm(x3)^IJ]. (2.2.22) n=l m=l k=l The sum over n yields 6ik and the sum over k can subsequently be performed and one obtains Dij given by Eq. (2.2.18) The second derivative matrix A^2' or D^2' can be obtained simply by multiplying it by itself. It is easy to show that the matrix elements for A^2' and D^2) satisfies Ag> = y/wj/wM?- (2-2.23) The application to differential equations is based on the algorithm for numerical differentiation, defined by, df N lf-U*, = EAvf(x3), (2-2.24) i=i or df N V^[^],=,, = E A,/(^)v%, (2.2,25) where the matrix A„ and are defined in equations (2.2.11) and (2.2.18) respectively. Generally speaking, the matrix representation of a differential equation can be written in a simple way by replacing derivatives with Aij or Dij and function with their values at set of points Xi. The advantage of introducing the new derivative matrix D will be demonstrated later in this chapter and in the applications as shown in chapter 3 and chapter 4. Chapter 2. The Quadrature Discretization Method (QDM) 33 The QDM method is analogous to the collocation method by using the classical sets of polynomials (e.g. Jacobi polynomial, Hermite polynomial, etc.) and their corresponding weight functions to generate weights and points to construct the differential matrix. The advantage of the QDM over the traditional pseudospectral method is that it allows the use of arbitrary orthogonal polynomial basis sets or equivalently arbitrary weight functions. This give us a much wider choice of the collocation points and weights to apply to a specific problem. The main idea is to choose a weight function related to the physical nature of the problem to be solved such that the solution converges rapidly. 2.2.2 Generation of orthogonal polynomials, quadrature weights and points There have been several discussions on the generation of a set of polynomials orthogonal to each other with respect to some weight function w(x). This is a numerically unstable problem [143,75] and a practical approach, referred to as the discretized Stieltjes's pro cedure, has been provided recently by Gautschi [134]. Basis set generation; The Stieltjes's procedure: The generation of a set of polynomials orthogonal with respect to some weight function w(x) has been discussed in several texts [135] and papers [71,72,75]. For a nonclassical weight function, the usual Schmidt procedure [144] which involves orthogonalization of a given member of the set to all the functions of lower order and then normalized, is highly unstable due to roundoff errors and is not practical. The Schmidt procedure is analogous to methods based on the moments of the weight function. The best a/pproach is one based on the three term recurrence relation of the polynomials {Qn{x)} with Qo = 1, given by Qn+i{x) = (x - an)Qn(x) - BnQn-^x). (2.2.26) Chapter 2. The Quadrature Discretization Method (QDM) 34 It can be shown that the two recurrence relations Eq.(2.1.5) and Eq.(2.2.26) are equivalent and satisfy cr, n Pv 1 (2.2.27) Pn = —, (2.2.28) PnPn-l anc Qn = — • (2.2.29) Pn-\pn-1---P\ For the classical polynomials, the coefficients an and Pn are known and exact. For example, an = 0, Pi = 1/3 and Pn = n2/(4n2 — 1) (n > 1) for the Legendre case; an = 0, Pi = 1/2 and Pn = 1/4 (n > 1) for the Chebyshev case; and an = 0, pi = 1/2 and Pn = n/2 (n > 1) for the Hermite case. The polynomials generated by Eq. (2.2.26) are not normalized. The corresponding normalized polynomials Pn(x) = <5n(^)/v/7n satisfy the recurrence relation xPn(x) = ^/Pn+iPn+1(x) + anPn(x) + ^JPnPn_,(x), (2.2.30) where the normalization factors are given by 7n = / w(x)Q2n(x)dx. (2.2.31) For the QDM, we are interested in the non-classical sets of orthogonal polynomials and an, Pn are usually unknown and calculated numerically. The recurrence coefficients in Eq. (2.2.26) are given by _ J w(x)xQn(x)dx an~ f w(x)Ql(x)dx > KLLil) and f w(x)Ql(x)dx / w(x)Q2n_1(X)dy Chapter 2. The Quadrature Discretization Method (QDM) 35 One of the practical approaches of numerical computation of an and Bn in Eqs. (2.2.32) and (2.2.33) is the Stieltjes's procedure introduced by Gautschi [134]. The method in volves the accurate calculation of the integrals in Eqs. (2.2.32) and (2.2.33) by subdivid ing the domain of interest into many subdomains and evaluating the contribution from each subdomain with a high order quadrature. Gautschi discusses the use of Gauss-Legendre and Feijer quadrature rules. The calculation begins with Q-i = 0, Q0 = 1 and ao = A) = 0. Equation (2.2.26) can then be used to generate Qi and then ot\ and B\ with Eqs. (2.2.32) and (2.2.33). In this way all the recurrence coefficients are calculated. The quadrature weights and points are calculated as discussed elsewhere [135,75]. The dis cretized Stieltjes's procedure was used to calculate the recurrence coefficients for several nonclassical weight functions reported in previous papers. It is found to work remarkably well as long as care is taken to ensure the accuracy of the integrals by a judicious choice for the number of subdivisions and the order of the quadrature in each interval. Since there is no analytic solution of an and 3n for most non-classical recurrence relations, the accuracy of the numerical calculation of an and Bn will be essential to the later calculation of quadrature points and weights as well as the derivative matrices. To investigate the accuracy of the Stieltjes's procedure we calculate the an and 8n for up to 100 quadrature points for the classical polynomials with the Stieltjes's procedure and compare the results with the known analytic solution. The maximum errors E°° in an and Bn computed by Stieltjes's procedure for several classical polynomials with TV quadrature points are shown in Table 2.2. The error curves for the an and Bn with N = 50 is given in Fig. 2.8. As shown in the table, the maximum errors in Bn for both Legendre and Hermite polynomials and in an for all the polynomials are very small and at most O(10~13). The maximum error in Bn for the Chebyshev polynomials are 7.76 x 10~5. It is unchanged with changing because the maximum error eventually occurs at the first coefficient B^. The reduced accuracy of Chapter 2. The Quadrature Discretization Method (QDM) 36 N Chebyshev 5 9 528058591081317e- 17 7 763988590681060e- 05 10 5 386806816078257e- 16 7 763988590681060e- 05 20 9 710341644410283e- 16 7 763988590681060e- 05 50 1 805638645361598e- 15 7 763988590681060e- 05 100 2 490836730411662e- 15 7 763988590681060e- 05 Legendre 5 7 070397677170780e- 16 3 053113317719181e- 15 10 7 070397677170780e- 16 3 053113317719181e- 15 50 4 679789671731484e- 15 3 053113317719181e- 15 100 6 021950962055382e- 15 3 164135620181696e- 15 Hermite 5 4 996208974320000e- 15 9 103828801930001e- 15 10 4 996208974320000e- 15 4 840572387370000e- 14 50 2 609915956900000e- 14 3 375077994860000e- 13 70 5 530643324450000e- 14 4 760636329590000e- 13 100 5 530643324450000e- 14. 7 034373084020000e- 13 Table 2.2: Maximum error in ctn and f3n calculated by the Stieltjes's procedure. the recurrence coefficients /3n's for the Chebyshev polynomial is due to the less accurate intergrations. Recall the weight function for Chebyshev polynomial is w(x) = m interval [-1,1] and it is divergent as x tends to the end points. When one integrates a polynomial with respect to this weight function numerically, a large round-off error is involved and results in poor accuracy of the integration. For all the computations, the maximum errors are insensitive to the changes in N. Theoretically, for polynomials defined in an infinite domain, the integrals for a's and /3's are evaluated in the entire domain. When numerical integration is applied, one has to truncate the domain into a finite one. If the truncated domain for the integration is Chapter 2. The Quadrature Discretization Method (QDM) 37 .20.00 1 1 ' ' ' 1 1 1 ' 1 0.00 10.00 20.00 30.00 40.00 50.00 n Figure 2.8: Errors in an and Bn calculated by the Stieltjes's procedure with N=50. The dash dotted curve is for Chebyshev polynomials, the dashed curve is for the Legendre polynomials and the solid curve is for the Hermite polynomials. Chapter 2. The Quadrature Discretization Method (QDM) 38 N E°° 01 N [-20, 20] 5 4 996208974320000e-15 9.103828801930001e-15 10 4 996208974320000e-15 4.840572387370000e-14 50 2 609915956900000e- 14 3.375077994860000e-13 70 5 530643324450000e- 14 4.760636329590000e-13 90 5 530643324450000e- 14 7.034373084020000e-13 100 5 530643324450000e- 14 7.034373084020000e-13 [-5, 5] 5 2 346651055977530e- 15 7.246219748680005e-06 10 1 263572315381132e- 14 2.919377614958041e-02 50 3 598600987376839e- 14 1.874888749856461e+01 100 8 320267776295266e- 14 4.374979589455547e+01 [-50, 50] 5 6 658312592764631e- 16 1.154631945610163e-14 10 1 382306109940420e- 15 3.907985046680551e-14 50 1 895834758617558e- 14 1.492139745096210e-13 70 3 565939958333915e- 14 2.060573933704291e-13 >91 overflow Table 2.3: Maximum error in the numerical solution of an and Bn for the Hermite poly nomials on three intervals [-20,20], [-5, 5] and [-50, 50]. Chapter 2. The Quadrature Discretization Method (QDM) 39 too small, the truncation error for the integration is large and results in low accuracy for the calculation of an and fin. If the domain is too large, the integrals can be very large and cause number overflow for the computer. Careful choice of domain cutoff for the numerical integration is therefore needed for the accuracy and stability of the results for the numerical calculation of an and /3n. The computer overflow in the integration of low degree polynomials is rare so that it is not necessary a worry if number of quadrature points required is not large. The error estimation of an and ftn for the Hermite polyno mials calculated by Stieltjes's procedure with several choices of integration intervals are given in Table 2.3. Calculation of quadrature points and weights If the an and f3n for the recurrence relation are known, the quadrature points which are the roots of the Nth polynomial constructed from the recurrence can be calculated by diagonalizing the tri-diagonal matrix [136] / J = 0 y/p\ a2 0 0 0 0 0 0 \ V (2.2.34) 0 0 0 ... ayy-i /SJV-I 0 0 0 ... PN-I CJV J The points are the eigenvalues of the matrix of 3. The corresponding weights are given in the terms of the first component of the ith normalized eigenvector of the matrix J [136,134]. In Table 2.4 we give the points and weights of several previous cases studied in the above section. In order to study the effect of the error of the calculated an and /3n, we Chapter 2. The Quadrature Discretization Method (QDM) 40 N poo poo Chebyshev 5 2.954894016793386e-05 20 7.757441346845083e-06 50 3.121430352637233e-06 100 1.569879214957659e-06 7.806354429684692e-05 2.321553559997414e-05 9.635848417718917e-06 4.898200317991069e-06 Legendre 5 2.220446049250313e-50 1.415534356397075& 100 3.885780586188049e-15 3.774758283725532e-15 15 5.384581669432010e-15 15 7.507883204027621e-15 Hermite [-20, 20] 5 4.440892098500626e-50 5.773159728050815e-70 7.460698725481052e-100 9.592326932761353e-15 1.032507412901395e-14 14 1.271205363195804e-14 14 1.068589661201713e-14 14 1.276756478318930e-14 [-5, 5] 5 1.495381356519943e-07 5.903613908841976e-08 50 4.192398352706889e+00 5.962550119618734e-02 100 8.408353932924946e+00 8.056316459982088e-02 [-50, 50] 5 2.664535259100376e-15 2.109423746787797e-15 50 5.417888360170763e-14 5.828670879282073e-15 70 4.618527782440652e-14 5.190292640122607e-15 Table 2.4: Maximum error for the numerical approximation of points and weights. Chapter 2. The Quadrature Discretization Method (QDM) 41 0.80 1.20 1.60 2.00 log10(N) Figure 2.9: Minimum spacing of points versus number of quadrature points N. The dash dotted curve is for Chebyshev points; the dashed curve is for Legendre points and the solid curve is for Hermite points. Chapter 2. The Quadrature Discretization Method (QDM) 42 compare the results calculated directly from the exact an and 8n, as well as the exact points and weights where it is applicable. • The results from Table 2.4 show no sign of the extra error involved for the calculated points and weights except the error from the errors of an and 8. If we can produce the recurrence coefficient accurately enough, the accuracy of the points and weights will not be affected. The minimum spacing between adjacent quadrature points is plotted in Figure 2.9, The minimum spacing satisfies Aminx ~ 9.87/N2, Aminx « 12.2/A^2, and Aminx & 2.22/v^/V for the Chebyshev, Legendre, and Hermite points, respectively. Again, the results show that the minimum spacing is decreased at rate of 0(1/N2) for the Chebyshev and Legendre points and 0(1 /s/N) for the Hermite points as number of quadrature points N increases. 2.3 Time integration In most applications to the PDE, the spatial discretization uses the QDM but the tem poral discretization uses conventional finite difference. Consider a time dependent PDE |^ = /(u,i) <>0, (2.3.1) where the general (linear or nonlinear) operator / contains the spatial part of the PDE. After discretizing the spatial operator by the QDM, the semi-discretized version of Eq. (2.3.1) becomes where U(t) is the set of grid values of u(x,t), F is the discretization representation of the operator /. In this section, we confine our discuss to the linearized version of (2.3.2), Chapter 2. The Quadrature Discretization Method (QDM) 43 i.e., we consider ft = LU' <2-3-3» where L is the matrix resulting from the linearization of (2.3.2). If F(U) is linear, i.e., L = F, and L is a diagonalizable matrix, we can solve Eq. (2.3.3) by eigenfunction expansion. Eigenfunction expansion The eigenvalue problem of Eq. (2.3.3) is defined by L(f>n(x) = -en<f)n(x), (2.3.4) where en are the eigenvalues and 4>n(x) are the eigenfunctions. The linear time-dependent equation, admits a solution of the form oo u(x,t) = 5>ne-e"Vn(x), (2.3.5) 71=0 where the coefficients cn are determined from the initial condition, u(x, 0), and are given by cn = J u(x,0)<f>n(x)dx. (2.3.6) In actual calculations, the solution u(x,t) given by Eq. (2.3.5) must be truncated at n = JV, and evaluated at quadrature points. It is essential that the lower eigenvalues and eigenfunctions converge accurately and rapidly in order to achieve a fast accurate approximation of the solution. Since no numerical approximation involved in temporal variable t, the accuracy of the solution of Eq. (2.3.3) is of the same order as the numerical method applied to approx imate the differential operator L. In the case that the eigenvalues are easy to calculate accurately and converge very fast, the eigenfunction expansion technique is a super choice Chapter 2. The Quadrature Discretization Method (QDM) 44 for the solution of the time dependent problem. The convergence could be difficult at small times depending on the choice of the initial condition. For most applications to solve time dependent problem, temporal discretization is a more popular approach. Temporal discretization Most temporal discretization uses a finite difference method. The eigenvalue spectrum of the discretization matrix in space determines the stability of the time-discretization. The region of the absolute stability of a numerical method is defined for the problem f = w <2 3 7> to be the set of all XAt such that \\U\ \ is bounded as t —• oo. A method is called A-stable if the region of the absolute stability includes the region Re{\At} < 0. The most simple and popular schemes widely used in the applications are the first order forward Euler (FE) and backward Euler (BE) method, as well as the second or der Crank-Nicolson (CN) method. The discretization formula for Eq. (2.3.3) of these methods are as following: for the forward Euler method, jjn+l = Tjn + AtLn ^ g For the backward Euler method, Un+1 = Un + AtLn+\ (2.3.9) and for the Crank-Nicolson method, Tjn+l =zrjn + }-&t[Ln+1 + Ln}. (2.3.10) The FE method is an explicit method and the absolute stability region is |1 + AAi| < 1. If the eigenvalue A is real, then —2 < XAt < 0. Both FE and CN are implicit methods and are A-stable, i.e. they are absolutely stable in the entire left-half plane. Chapter 2. The Quadrature Discretization Method (QDM) 45 N QDM1 QDM2 QDM3 QDM4 4 7.08(-03) 3.11(-03) 9.26(-02) 9.62(-02) 6 1.59(-02) 4.47(-03) 9.26(-03) 1.18(-02) 8 5.25(-04) 1.99(-04) 4.01(-04) 6.10(-04) 10 2.66(-05) 5.43(-06) 1.04(-05) 1.85(-05) 12 4.43(-07) 9.86(-08) 1.84(-07) 3.75(-07) 14 8.77(-09) 1.29(-09) 2.36(-09) 5.43(-09) 16 8.59(-ll) 1.28(-11) 2.31(-11) 5.90(-ll) 20 1.20(-14) 4.44(-15) 3.48(-15) 6.78(-15) 25 4.95(-14) 7.35(-15) 1.05(-14) 1.23(-14) Table 2.5: Standard error for the numerical approximation of u" = sin(7ra;) for the four QDM methods. The detailed discussion of these methods and other high order and more complicated numerical schemes are out of the scope of this thesis and can be found in standard references [1-6]. 2.4 Numerical tests and results In this section, we give three analytic examples to study the properties of the QDM. 2.4.1 One-dimensional Poisson equation The first example we discuss is the one-dimensional Poisson problem in section 2.1 defined by Eq. (2.1.12). We compared the convergence of the solution of the Poisson problem between the classical Legendre method and a second order FD method. Since the in tention of the QDM is to use nonclassical base functions, in this section we applied four nonclassical polynomial sets based on four weight functions. They are wx(x) = sin2(7roj), w>2(x) — ex2, wz(x) = e~x2, and 104(2) = 1 — x2. We refer to the QDM based on these weight functions as QDM1, QDM2, QDM3 and QDM4, respectively. The standard errors of the numerical solution for the four methods are given in Table Chapter 2. The Quadrature Discretization Method (QDM) 46 2.5. All the methods work remarkly well and the rate of convergence is exponential. With 20 quadrature points, the numerical approximations of the solution is at least 14 digits accurate. The rate of convergence for the QDM competes well with that for the classical Legendre method and the method is much faster than the FD method (see Table 2.1). We also study the largest eigenvalue of the derivative operators based on the nonclassical QDM, the result is similar to the classical methods, that is, the largest eigenvalue of the first derivative matrix grows like 0(N2) and of the second derivative matrix grows like 2.4.2 An analytic example In order to demonstrate the applicability of the QDM and some of its properties, we consider a time-dependent Fokker-Planck equation defined by The physical interpretation of P(x,t) is the probability density function. It is defined in x £ (—00,00), and satisfies P(±oo,£)=0. With a ^-function initial condition at xo, Eq. (2.4.1) has an analytic solution 0{N4). where (2.4.2) where z = e 1/4. The equilibrium solution as t —>• 00 is given by (2.4,4) If we set P(x,t) = Po(x)g(x,t), then the equation for g(x,t) is given by, (2.4.5) Lg(x,t), Chapter 2. The Quadrature Discretization Method (QDM) 47 It is easily shown that the operator L is self-adjoint with respect to the equilibrium so lution PQ{X). 1. Eigenvalue problem. First, we consider the eigenvalue problem of Eq. (2.4.5): -L<f>n(x) = [-Aix)-^) + B(X)^}MX) = -KMX), (2-4.6) where An and <f)n are the eigenvalues and eigenfunctions, respectively. The exact eigen values of Eq. (2.4.6) are An = n/4, 7i = 0,l,2,.... (2.4.7) Three set of polynomials are applied in the QDM to solve the eigenvalue problem Eq. (2.4.6). They are scaled Legendre polynomials with weight function w\(x) = 1 in a finite interval, scaled Hermite polynomials with weight function w(x) = e~4x'2, and nonclassical polynomials with weight function w(x) = Po(x). We refer to the methods based on these polynomials as Legendre, Hermite, and QDM, respectively. The discretization matrix LN of the differential operator in Eq. (2.4.4) is in the form -LN = -AA + BA2 (2.4.8) where, A and B are the diagonal coefficient matrices of A(x) and B(x), A is the derivative matrix with its entries defined in Eq. (2.2.11). If we use the equilibrium solution Po(x) as weight function, the Galerkin matrix LN of the operator L can be obtained in a simple form LN = BDDT, ' (2.4.9) where D is the modified QDM derivative matrix with its elements defined in Eq. (2.2.18). The detailed derivation of the matrix is given in Chapter 3. We refer this approach of Chapter 2. The Quadrature Discretization Method (QDM) 48 N Ai A5 Aio Al5 A20 Exact 0.25 1.25 2.5 3.75 5.0 Error Hermite 4 .218(-01) 6 .123(-02) .465(+01) 10 .592(-06) .524(+00) .186(+02) 15 .205(-12) .682(-02) .278(+01) .319(+02) 20 .104(-13) .889(-04) .126(+00) .616E+01) .458(4-02) 25 .905(-14) .146(-07) .504(-01) .654(+00) .103(4-02) 30 .133(-13) .184(-08) .192(-01) .919(+00) .171(4-01) 40 .508(-13) .389(-12) .681(-04) .574(+00) .168(4-01) 50 .143(-13) .913(-13) .777(-07) .285(+00) .141(4-01) QDM 4 .205(-02) 6 .467(-04) .204(+00) 10 .530(-08) .242(-02) .576(+01) 15 .722(-15) .308(-05) .282(+00) .124(+02) 20 .999(-15) .212(-07) .104(-01) .140(4-01) .213(4-02) 25 .000(4-00) .206(-ll) .506(-04) .230(4-00) .358(4-02) 30 .666(-15) .422(-14) .453(-07) .152(-01) .100(4-01) 40 .555(-15) .666(-14) .151(-13) .767(-06) .227(-01) 50 .355(-14) .000(+00) .444(-14) .439(-12) .577(-05) MQDM 4 .154(-03) 6 .257(-05) .388(+00) 10 .347(-07) .216(-01) .306(+01) 15 .158E-08) .180(-06) .408(+00) .808(4-01) 20 .140(-11) .999(-09) .447(-02) .151(4-01) .153(4-02) 25 .934(-13) .107(-11) .193(-03) .402(-02) .369(4-01) 30 .433(-14) .102(-13) .163(-07) .315(-01) .559(4-00) 40 .278(-15) .666(-15) .431(-13) .312(-05) .759(-02) 50 .472(-14) .688(-14) .160(-13) .276(-ll) .240(-05) Table 2.6: Numerical error in the eigenvalues of the analytic example for the Hermite and QDM methods. Chapter 2. The Quadrature Discretization Method (QDM) 49 N Ai A5 Aio Al5 A20 Exact 0.25 1.25 2.5 3.75 5.0 Error Legendrea 4 .295(-02) 6 .102(-03) .500(+00) 10 .143(-06) .230(+00) .845(4-01) 15 .105(-11) .192(-01) .259(4-01) .207(4-02) 20 .156(-13) .140(-02) .978(4-00) .823(4-01) .386(4-02) 25 .354(-10) .125(-04) .395(4-00) .453(4-01) .174(4-02) 30 .491(-08) .242(-06) .854(-01) .268(4-01) .113(4-02) 40 .462(-04) .691(-03) .128(-02) .672(4-00) .466(4-02) 50 .414(+00) .341(+00) .583(4-00) .523(4-00) .137(4-01) Legendreb 4 .250(-01) 6 .262(-02) .761(+00) 10 .468(-06) .292(+00) .186(4-01) 15 .200(-12) .568(-01) .928(4-00) .606(4-01) 20 .247(-14) .566(-05) .424(4-00) .128(4-01) .127(4-02) 25 .164(-14) .143(-07) .656(-01) .796(4-00) .118(4-02) 30 .239(-14) .770(-ll) .199(-03) .845(-02) .511(4-01) 40 .711(-14) .134(-12) .128(-09) .109(-03) .440(-00) 50 .799(-14) .222(-12) .695(-12) .102(-05) .728(-00) Legendre0 4 .105(+00) 6 .313(-01) .901(+00) 10 .297(-04) .739(+00) .255(4-00) 15 .137(-09) .494(4-00) .168(4-01) .133(4-01) 20 .927(-14) .314(4-00) .139(4-01) .257(4-01) .424(4-01) 25 .255(-14) .465(-01) .109(4-01) .229(4-01) .343(4-01) 30 .555(-14) .942(-05) .916(4-00) .194(4-01) .317(4-01) 40 .389(-15) .443(-ll) .388(4-00) .135(4-01) .250(4-01) 50 .178(-14) .254(-10) .216(-05) .841(4-00) .183(4-01) Legendre points in a (-2,2); 6 (-3,3); c (-4,4). Table 2.7: Numerical error in the eigenvalues of the analytic example for the Legendre method. Chapter 2. The Quadrature Discretization Method (QDM) 50 discretization as modified QDM (MQDM). The main advantage of using the MQDM representation is that the differential matrix is simple to construct. The matrix repre sentative is symmetric and the eigenvalues are all real. The matrix representations in Eq. (2.4.8) and Eq. (2.4.9) are diagonalized to give the approximate eigenvalues and eigenfunctions. The convergence of the eigenvalue A1; A5, Aio, A15, and A2o is shown in Table 2.6 and 2.7 for each method. The exact eigenvalues are given in the top of the tables. Table 2.6 compares the errors of these eigenvalues calculated by the Hermite, QDM and MQDM. As we can see from the table, QDM and MQDM give relatively the same rate of convergence for most eigenvalues, and the error of eigenvalues decays very fast. As A^ increases, the convergence of the first eigenvalue Ai for the MQDM is slightly slower than that for the QDM. To obtain 14 decimal places accuracy of Ai it requires A^ = 30 points for the MQDM in comparison with N = 15 points for the QDM. The Hermite method calculates small eigenvalues (e.g. A1; A5) very well but the results for the larger eigenvalues (e.g. Ai5, A2o) are not as good. The eigenvalue A20 is only accurate to 0(1) for the Hermite method with A^ = 50 in comparison with O(10~5) for the QDM and MQDM. For the QDM, the points spread over the interval (—00,00). As N increases, the domain is getting larger. Since we use the equilibrium function Po(x) as the weight function, the resulting domain reflects the natural property of the solution. At the boundary points, the rate of decay of the the solution is the same as the exact solution. We may refer to this domain as natural domain of the problem. With the same number of quadrature points on a larger domain, some points are beyond the natural domain and can not contribute to improve the accuracy. Therefore the rate of convergence is expected to be slower. If the domain is much smaller than the natural domain, the boundary condition may not be accurately applied, which results in low accuracy in the solution. Therefore, choosing the weight function is very crucial for the convergence of Chapter 2. The Quadrature Discretization Method (QDM) 51 Hermite, N=8 OOOOOOOOO Hermite, N=32 0.5 -0.5 -1 0 • o o o o ooaasao 0 -100 -50 Legendre, N=8 Legendre, N=32 0.5 CO c ca E -0.5 o o o oo o o o o -2 -1 Real Figure 2.10: Numerical approximation of the eigenvalues for the analytic FPE problem. Chapter 2. The Quadrature Discretization Method (QDM) 52 N QDM Hemite 4 (-1.834, 1.834) (-1.010, 1.010) 6 (-2.077, 2.077) (-1.326, 1.326) 10 (-2.375, 2.375) (-1.834, 1.834) 15 (-2.605, 2.605) (-2.344, 2.344) 20 (-2.764, 2.764) (-2.775, 2.775) 25 (-2.885, 2.885) (-3.155, 3.155) 30 (-2.983, 2.983) (-3.498, 3.498) 40 (-3.136, 3.136) (-4.106, 4.106) 50 (-3.252, 3.252) (-4.632, 4.632) Table 2.8: Domain generated from the QDM and Hermite weight function for the analytic FPE example. the solution. We will discuss more about the choice of the weight function in Chapter 3 and Chapter 4. Table 2.8 gives the domains for the QDM and Hermite method versus number of grid. To apply the Legendre points, one has to truncate the infinite domain ( —oo, oo) into a finite domain (—c, c). Table 2.7 gives the numerical error of these eigenvalues calculated by Legendre method in the intervals a: (-2,2), b: (-3,3), and c: (-4,4), respectively. We refer to the above three Legendre methods as Legendre", Legendre6 and Legendrec method, respectively. As we can see from the table, the Legendre6 method on interval (-3,3) is the best among the three Legendre methods. As N increase, the error of the eigenvalues calculated by the Legendre" method decays initially, and then increase. This is because the truncated domain is too small, and with N increasing, the truncation error caused by the boundary cut off accumulates and is getting larger and larger. The Legendre0 method is in general slower than Legendre6 method, since the domain (-4,4) is wider than it requires and some points near the boundary are not contributing. In comparison with the QDM and MQDM, the convergence of the Legendre6 method is Chapter 2. The Quadrature Discretization Method (QDM) 53 slower, especially for the large eigenvalues, but it is better than the Hermite method. For the Hermite and. Legendre methods, some of the eigenvalues of the differential matrices are imaginary as shown in Fig. 2.10 with N — 8 and N = 32. For the QDM and MQDM, all the eigenvalues are real. This property is obviously important with regard to CPU time and memory, especially for higher dimensional problems. Fig. 2.11 compares the largest eigenvalues of the differential matrix versus number of quadrature points for each method. The largest eigenvalue grows like C^TV1-76) for the QDM, 0(A^1-85) for the MQDM, C^/V1-2) for the Hermite, and 0(N2) for the Legendre method as ./V increases (N < 100). It worth to point out that, although the largest eigenvalue of the second derivative matrix for the Legendre method grows at much faster pace ( 0(N4)), this property is not shown in this application. 2. Time dependent solution We intend to solve the Eq. (2.4.1) with a ^-function at x0 as initial condition. The application of the eigenvalue expansion stated in section 2.3 in the solution the the time-dependent problem with x0 chosen as one of the quadrature points is straight forward. In this section our interest is to apply the temporal propagation techniques such as forward Euler or/and backward Euler method to solve the time-dependent the FPE problem. We set x0 = 1. Since it is very hard to discretize the ^-function numerically based on the quadrature points, in order to compare the numerical result with the exact solution, we use the exact solution at time t=l as the initial distribution. Both forward and backward Euler are used to solve the problem. Since the real part of the eigenvalues are all negative for the QDM, Hermite and Legendre differential matrices, backward Euler method is always stable, the choice of the time step is depended on the accuracy requirement of the solution. For the forward Euler method, the time step At has to satisfy the condition At < 2/\max, where Amax is the largest eigenvalue of the differential Figure 2.11: Largest numerical eigenvalues for the analytic FPE example. The results with the QDM, MQDM, Hermite and Legendre are denoted by solid lined with asterisks, crosses, solid circles and diamonds, respectively. Chapter 2. The Quadrature Discretization Method (QDM) 55 N QDM Hermite Legendre 10 0.24 0.094 0.45 15 0.12 0.056 0.20 20 0.075 0.039 0.11 25 0.051 0.030 0.071 30 0.037 0.024 0.048 35 0.028 0.020 0.035 40 0.022 0.017 0.026 50 0.015 0.013 0.016 60 0.011 0.010 0.011 70 0.0085 0.0093 0.0084 80 0.0067 0.0079 0.0063 90 0.0054 0.0070 0.0050 Table 2.9: Maximum time step in the FE time integration. matrix. Table 2.9 gives the maximum time step for the forward Euler method for the QDM, Hermite and Legendre differential matrix. The numerical solution of Eq. (2.4.1) solved by the forward Euler method with Ai = 0.01 is given in Fig. 2.12 and Fig. 2.13 with N = 6 and TV = 15 quadrature points respectively. In the figures, the solid line is the exact solution. As we can see from Fig. 2.12, with only six quadrature points, the QDM approximates the exact solution very well. Among the three discretization methods in spatial dimension, the QDM approxi mation is the best. The Legendre method is better than the Hermite method, but still a little bit off the exact solution. Hermite method approximate poorly to the exact solu tion for large times. With increasing number of quadrature points to = 15 (Fig. 13), all the numerical approximations can not distinguish from the exact solution in the figure. 3. An equivalent problem If we make a simple change of variable y = \ sinh x and let P(y, t) — \ cosh xP(x, t) Chapter 2. The Quadrature Discretization Method (QDM) 56 Figure 2.12: Numerical solution of the time dependent FPE with N = 6. The solid lines are the exact solution. The plus signs are the results with the QDM, the asterisks are the result for the Legendre method and the circles are the results with the Hermite method. Chapter 2. The Quadrature Discretization Method (QDM) 57 Figure 2.13: Numerical solution of the time dependent FPE with TV = 15. The solid lines are the exact solution. The plus signs are the results with the QDM, the asterisks are the result for the Legendre method and the circles are the results with the Hermite method. Chapter 2. The Quadrature Discretization Method (QDM) 58 in Eq. (2.4.1) and (2.4.2), one finds the equation, dP(y,t) = l,dP{y,t) ldP\y,t) dt 4l dy 2 dy2 (2.4.10) The equilibrium solution of this equation is dP(y,t) = ldP(y,t) ldP2(y,t) dt 4 dy 8 dy2 The eigenvalues of its eigen-problem are rc/4, (n = 0,1,N) and the eigenfunctions are the normalized Hermite polynomials Hn(y). We solve the eigen-problem numerically by the QDM and MQDM with w(y) = Po(y) and Legendre methods in the interval (-6,6). Since the eigenfunctions are polynomials, we expect accurate results in numerical ap proximations. Table 2.10 gives the convergence of the eigenvalues for the three methods. For the QDM and MQDM, the solutions are exact. While for the Legendre method, the eigenvalues are exact for small N and the error is getting larger with increasing N after a certain level of N. The poor approximations are a result of the domain cutoff for the Legendre method. Furthermore we should point out that, the Galerkin matrix in transform space (see Eq. (3.2.12), Chapter 3 for detail) is diagonal with the Hermite polynomials as basis set. The diagonal elements are the eigenvalues of the eigen-problem of Eq. (2.4.12). The eigenfunction matrix is the transform matrix defined in Eq. (2.2.15). In comparison with the current eigen-problem, the convergence of the eigenvalues of the eigen-problem Eq. (2.4.6) is much slower because the eigenfunctions of the problem are not polynomials. Since the QDM and MQDM are exact for the problem, it is easy to know that the largest eigenvalue is O(N). For the Legendre method we applied, the largest eigenvalue Chapter 2. The Quadrature Discretization Method (QDM) N domain Ai A5 Aio Al5 A20 Exact 0.25 1.25 2.5 3.75 5.0 Error Hermite 4 (-2.020, 2.020) .555(-16) 10 (-3.668, 3.668) .222(-15) .222(-15) .111(-13) 20 (-5.550, 5.550) .111(-14) .266(-14) .311(-14) .755(-14) .109(-12) 30 (-6.996, 6.996) .178(-14) .155(-14) .666(-14) .933(-14) .400(-13) 40 (-8.213, 8.213) .888(-15) .799(-14) .933(-14) .933(-14) .320(-13) 50 (-9.284, 9.284) .355(-14) .666(-15) .266(-14) .355(-14) .258(-13) 70 (-11.14, 11.14) .160(-13) .622(-14) .311(-14) .444(-14) .977(-14) 90 (-12.74, 12.74) .738(-14) .533(-14) .191(-13) .933(-14) .888(-15) Legendre 4 (-6, 6) .111(-15) 10 .222(-15) .205(-12) .129(-11) 20 .833(-15) .566(-12) .183(-09) .141(-09) .768(-10) 30 .350(-14) .169(-12) .133(-09) .317(-08) .549(-07) 40 .173(-13) .187(-12) .223(-09) .647(-08) .342(-06) 50 .533(-14) .571(-13) .127(-10) .331(-07) .128(-05) 70 .279(-13) .400(-09) .250(4-00) .496(4-00) .994(4-00) 90 .192(-13) .240(+01) .449(4-01) .100(4-01) .979(4-00) Table 2.10: Comparison of the numerical error of the eigenvalues of Eq. (2.4.12) for Hermite and the Legendre methods. Chapter 2. The Quadrature Discretization Method (QDM) 60 is marginally O(N). One should notice that the largest eigenvalue approximated by the Legendre method for this problem is much smaller than the largest eigenvalue of the second derivatives (which is 0(N4)). It occurs to us that the condition of a differential matrix is related to both differential operator and boundary conditions. It doesn't have to be of the same order as of the highest derivative in the differential operator. For the time integration by the explicit forward Euler method, the time step restriction should be OiN'1). As we recall the maximum time step is around 0(N2) for the time integration for the solution of Eq. (4.2.1) and (4.2.2). For the above examples, it is interesting to point out that although the condition number (or the largest eigenvalue) of the first and second derivative matrices depends on the minimum spacing of the grid as we discussed previously, it doesn't seem true for the differential operator. For the Legendre points, the minimum spacing is 0(iV2), and the condition number of the second derivative matrix is 0(N4). However, the condition numbers of the differential matrix representatives for the above examples are 0(N2) and O(N), respectively. 2.4.3 Three-dimensional Poisson problems In one dimensional examples, we already show that the QDM provides very high accuracy and fast convergence in the solution of differential equations. In comparison with the FD method, the QDM usually requires much coarser grid of points for the solution to converge to the same accuracy. This indicates a great advantage of the QDM in solving high dimensional problems. Our interest in this section is to solve the three dimensional Poisson equation of the form, (2.4.13) Chapter 2. The Quadrature Discretization Method (QDM) 61 f(x,y,z) Boundary Condition u(x,y,z) A 3g(^+J/+2) exact solution e(x+y+z) B -207T2^(x, y, z)* + sin(47ra;) sin(27ry)a(z)* u\Bc = 0 * <j>(x, y, z) = sin(47ra;) sm(2%y)[l - e-5^-1)], a(z) = 10(1 - lO^e-5^2-1)]. Table 2.11: Examples of Poisson problems subject to Dirichlet boundary conditions. The specific examples in Table 2.11 are con sidered. Numerical solution of example A is calculated by a second order central differencing FD method and the QDM on a IBM RS6000 computer. Two set of polynomials are applied in the QDM. One is the classical Chebyshev polynomials. Since the main in tention of the QDM is to use nonclassical polynomial basis set, we also applied a set of nonclassical polynomials associated with an arbitrary weight function w(x) = e6x2. For convenience, we use the same set of polynomials and the same number of grid points in each dimension. As we will see later in Chapter 4, sometimes it is necessary to use differ ent grids in different dimension in order to obtain superior convergence of the solution. The standard errors (E2) of the solution and the CPU time of the calculation.for each numerical method are given in Table 2.12. In the table, N is the number of grid points. For the QDM with nonclassical polynomials, the first CPU time is for the calculation of the points and the QDM derivative matrix associated with the weight function w(x) = e6x2. The second CPU time is for solving the discretized algebraic equation. As seen from the table, the error of the solution calculated with the FD method decays at the rate Chapter 2. The Quadrature Discretization Method (QDM) 62 FD Chebyshev QDM •N E2 CPU N E2 CPU N E2 CPU1 CPU2 4 .754(-2) 0.01s 4 .486(-3) 0.01s 4 .113(-2) 1.72s 0.00s 10 .143(-2) 0.03s 6 .263(-5) 0.03s 6 .882(-5) 1.76s 0.00s 20 .366(-3) 0.39s 8 .731(-8) 0.05s 8 .254(-7) 1.77s 0.01s 30 .163(-3) 1.99s 10 .136(-10) 0.10s 10 .478(-10) 1.81s 0.03s 40 .920(-4) 6.75s 12 .122(-12) 0.21s 12 .145(-12) 1.83s 0.08s 50 .589(-4) 16.5s 16 .522(-13) 0.50s 16 .631(-13) 1.86s 0.15s 70 .300(-4) 61.1s 90 .182(-4) 182.s Table 2.12: Standard Error and CPU time for the numerical approximations of example A of the Poisson problems. of 0(1/N2) as N increases, whereas the error of the solution calculated with the QDM decays at a much faster rate. With N = 10, or 10 x 10 x 10 grid, E2 = 0.136 x 10~10 for the Chebyshev approximation, and E2 = 0.478 x 10_1° for the QDM. For the FD method, E2 = 0.589 x 10~4 with N = 50, or 50 x 50 x 50 grid. To achieve accuracy of O(10-12), only a 12 x 12 x 12 grid is needed, and it costs less than a second for both QDM and Chebyshev method. However, for the FD method, with 90 x 90 x 90 grid, it takes about 3 minutes and one can only obtain accuracy of O(10~4). Comparing the CPU time, the QDM is much faster than the FD method if the same accuracy of the solution is required. Although the FD method performs relatively poorly in the calculation of Example A in comparison with the QDM, the result calculated by it is still acceptable. In some cases, the FD method may fail while the QDM can do very well. To show this point, we studied example B of the Poisson problems described in Table 2.11. Fig. 2.14 plots the exact solution of this problem at z = 0, i.e., u(x, y, 0) = (1 — e5) sin(47ra;) sin(27rj/). As seen in the figure, the solution of the problem is highly oscillated in x and y dimensions. The Poisson problem is solved by the three methods used in Example A. Fig. 2.15 shows the Chapter 2. The Quadrature Discretization Method (QDM) 63 TV FD Chebyshev QDM E2 CPU E2 CPU E2 CPU 10 .218(02) 0.04s .154(02) 0.06s .272(02) 1.77s 16 .732(01) 0.19s .121(00) 0.18s .375(00) 1.89s 20 .452(01) 0.44s .144(-02) 0.44s .492(-02) 2.33s 26 .261(01) 1.19s .449(-06) 1.18s .186(-05) 3.22s 30 .194(01) 2.20s .947(-09) 2.07s .460(-08) 4.19s 34 .150(01) 3.67s .137(-11) 3.49s .719(-11) 5.59s 40 .108(01) 7.15s .267(-12) 6.73s .849(-12) 8.94s 50 .688(00) 17.4s 70 .349(00) 36.5s 90 .211(00) 118.s Table 2.13: Comparison of the standard Error E2 and CPU times for the three numerical approximations of example B of the Poisson problems. standard error of the numerical solution versus logwN for the three numerical methods. The comparison of errors is given in Table 2.13. The results indicate that the numerical solution calculated by the FD method con verges very slowly. Even with a 90 x 90 x 90 grid, the FD solution barely reaches 1 decimal place accuracy. While for the Chebyshev method and the QDM, the numerical solution can achieve 11 decimal place accuracy with 40 x 40 x 40 grid. For both examples, it is evident that the QDM based on nonclassical polynomials competes well with the traditional Chebyshev method in the accuracy and convergence of the solution of PDEs. Chapter 2. The Quadrature Discretization Method (QDM) Figure 2.14: Exact solution for example B at Chapter 2. The Quadrature Discretization Method (QDM) 65 4.00 -16.00 1 1 1 1 1 1 0.80 1.20 1.60 2.00 log10(N) Figure 2.15: Standard error of the solution for example B. Chapter 3 Applications of the QDM to the Solution of Fokker-Planck Equation 3.1 Introduction The Fokker-Planck equation was introduced by Fokker and Planck to describe the Brow-nian motion of particles [76,77]. For the past several decades, there has been an on going interest by numerous researchers in the description of nonequilibrium phenomena modeled with a Fokker-Planck equation (FPE). This interest continues unabated to the present date [78] - [91]. The basis for many of these models is the Brownian diffusion in a potential characterized by Gaussian white noise. This leads to a time dependent linear FPE with drift and diffusion coefficients which can be nonlinear functions of the independent variable of interest. The theoretical basis for this approach has been pro vided in several standard references [92] - [94]. A general Fokker-Planck equation in one dimension is of the form dP(x,t) dA(x)P(x,t) d2B(x)P(x,t) dt ~ dx + d*x • [6AA) where P(x,t) is related to probability density function of x and A(x) and B(x) > 0 are referred to as the drift and diffusion coefficients, and depend on the particular application considered. The detail of the derivation of the Fokker-Planck equation can be found in the book by Risken [93]. The main objective in many of studies in Fokker Planck equation is the approach of some initial probability density function (PDF) to equilibrium, which in some ap plications may be characterized by two or more stable states. Examples of systems 66 Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 67 studied with a FPE include, model systems [78]- [80], [145]- [146], electron relaxation in gases [147], reactive systems [83,148,149], polymer dynamics [84], optical bistability [70,150,151], nucleation [152], dielectric relaxation [153], climate models [154], biological applications [155], astrophysical problems [86,100,156], economics [89], ionospheric ap plications [87,91], plasma physics [88], nuclear dynamics [90], and numerous other appli cations. These Fokker-Planck equations are strictly linear although the coefficients may be nonlinear functions of the independent variable. The main objective of the present chapter is to provide an efficient and rapid numerical method of solution of this large class of Fokker-Planck equations. The traditional method for the solution of the FPE usually involves the expansion of the probability density function in a suitable basis set, and the reduction of the differen tial equation to a set of algebraic equations for the expansion coefficients. This is referred to as spectral method of solution as discussed at length by Funaro [39] and Canuto et al [18,36]. An alternate approach involves the discretization of the PDF on a grid of points. This discrete approach in the solution of differential and/or integral equations has been used by researchers in other fields, notably neutron transport [95], radiative transfer [96], and computational fluid dynamics [18,36]. If the grid corresponds to the roots of the Nth order polynomial of some basis set, the solution procedure is referred to as a pseudospec tral approach as described by Fornberg [19]. Fourier series or Chebyshev polynomials are almost exclusively chosen as basis functions in the application of the pseudospectral ap proach. Other popular discretization schemes are based on the finite-difference technique such as those proposed by Chang and Cooper [97], Larsen et al [98] and Epperlein [99] primarily for the solution of a nonlinear FPE that arises in plasma physics. Park and Petrosian [100] have recently provided a detailed comparison of several different methods of solution of Fokker-Planck equations applied to astrophysical problems. The purpose of the present work is to apply the QDM to the solution of several Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 68 Fokker-Planck equations. The FPE describes the relaxation of some non-equilibrium distribution, P(x,t), to a "steady distribution", Po(x), at infinite time which can be written in the form, This steady distribution, characterized by the "potential" U(x), is the eigenfunction of the Fokker-Planck operator with zero eigenvalue, that is the "ground" state. In the QDM, the steady distribution given by Eq. (3.1.2) is employed as the weight function for the polynomials used to generate the quadrature grid. In the present work, we consider these weight functions together with alternate choices in order to maximize the rate of convergence of the solutions. We use Gautsch's Stielties procedure [134] to generate the polynomial sets orthogonal with arbitrary weight function. This permits the construction of polynomial basis sets and the corresponding quadrature grids for any appropriate weight function that accelerates the convergence of the solutions. Although the present work is restricted to Markovian processes and a linear FPE, it is anticipated that the present methodology will be easily adapted to nonlinear and multidimensional problems. The distinction between the FPE considered here with drift and diffusion coefficients which are nonlinear functions of the independent variable x and "truly" nonlinear equations was made clearer by Drozdov and Morillo [157]. For the non-Markovian situations, several workers have derived approximate FPEs in one and two dimensions [158]. We anticipate that the QDM will provide an accurate and efficient approach to these problems. (3.1.2) with (3.1.3) Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 69 3.2 The QDM matrix representation of the Fokker-Planck operator Equation (3.1.1) can be rewritten as dP(x,t) d dt dx A(x)P(x,t) + dB{x)P{x,t) dx (3.2.1) and has a stationary solution at infinite time given by Eq. (3.1.2). If we set P(x,t) = Po{x)g(x,t), then the equation for g(x,t) is given by, dg(x,t) dt -AW^^l + Bixf9^ dx = -Lg(x,t), dx2 (3.2.2) It is easily shown that the operator L is self-adjoint with respect to the equilibrium solution Po(x). Direct application of the QDM derivative matrix to the "spatial" operator L gives the matrix "representative" N N N Lgi = -A(xi) Aij9j + B(xi) J2 J2 A3kg{xk)-j=l j=l k=l (3.2.3) which is not symmetric and its eigenvalues may be imaginary. It is therefore important to construct a symmetric representative of the operator L. Let's consider the matrix representative of operator L in the transform space, Lr J w(x)Fn(x)LFm(x)dx, (3.2.4) where Fn(x) is the set of polynomials orthonormal with respect to the weight function w(x), that is J w(x)Fn(x)Fm(x)dx = Snm. (3.2.5) For convenience, it will be useful to introduce a second set of orthogonal functions, defined by Q n — \ to a; Po(x Fn{x) (3.2.6) Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 70 which are orthogonal with respect to equilibrium solution PQ{X), that is, J P0(x)Qn(x)Qm(x)dx = cw (3.2.7) The matrix elements of the operator L in the basis set {Qn(x)} are given by Lnm = J P0(x)Qn(x)LQm(x)dx, (3.2.8) = J P0(x)Qn(x)(A(x)Q'm(x)-B(x)Q'm(x))dx. Integrating by parts on Eq. (3.2.8), Lnm = -B(x)P0{x)Qn{x)Q'm(x) (3.2.9) + J A(x)P0(x)Qn(x)Q'm(x) + (B(x)P0(x)Q'n(x)Qff_(x)dx, = -B(x)P0(x)Qn{x)Q'm + J (A(x)P0(x) + (B(x)P0(x))')Qn(x)Q'm{x) + B(x)P0(x)Qn(x)'Q'm(x)dx, and using the facts that Po(x) = 0 at the boundaries and A(x)P0(x) + dB(x)P°(x) = Q for all x on the defined interval, we obtain that Lnm = J P0(x)B(x)Q'n(x)Q'm(x)dx. (3.2.10) In terms of the polynomial basis {Fn(x)} and weight function w(x), upon which the quadrature rule is based, we have, substituting Eq. (3.2.6) into Eq. (3.2.10), that Lnm = J B(x)w{x)[-^ + q(x)]Fn(x)[~ + q(x)]Fm(x)dx, (3.2.11) where w'(x) P^{x) = ^ - ^-T. (3.2.12) HK ' 2w(x) 2P0{x) V ; If the weight function w{x) could be chosen equal to PQ{X) then q(x) = 0. Applying the quadrature rule Eq. (2.1.9), we have that N Lnm ~ E B{xk)wk[F'n{xk) + q(xk)Fn(xk)][Fm(xk) + q(xk)Fm(xk)], (3.2.13) k=i Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 71 where xk, k = 1,2, ...,N are the quadrature points which are the zeros of FN(X). The matrix representative Lnm is transformed back to the physical space with transformation matrix T defined in Chapter 2, that is, T,-n = y/wlFn[xi) and we obtain N N Lij — ^ ' / ~] TinLnmTmj, (3.2.1.4) n=l m=l Substituting Eq (3.2.13) and transformation matrix T into Eq.(3.2.14), that is, N N N Lv = EE \/wlFn(xi)Y^B(xk)wk[Fll(xk)+q{xk)Fn(xk)][F^ 71=1 771 = 1 k = \ and using the fact that and (3.2.15) wkFn(xk)Fn(x,) = 8%k (3.2.16) N 71=1 N J2wkFn{xk)Fm(xk) = 6NM, (3.2.17) k=i one finds that N Lij ~ zZ B(xk)[Dkl + q(xK)6ki][Dkj + q(xk)Skj], (3.2.18) fc=i where are the elements of derivative matrix defined in Eq. (2.2.18). It is obvious that this matrix representation is symmetric. For the case when the weight function w(x) = Po(x), or q(x) = 0, the matrix elements of operator L, Lij, is simplified as TV Lij = E B(xk)DkiDkj (3.2.19) k=\ The great advantage of the QDM is that the matrix representation of the Fokker-Planck operator is easily constructed and evaluated for arbitrary coefficients A(x) and B(x). Since any set of orthogonal polynomials could be employed, it provides more opportunity to choose the weight function or basis set to optimize the matrix for rapid convergence and high accuracy of the solutions. It is expected that the convergence would be rapid for w(x) = PQ(X) although this is not always the case. The use of the modified Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 72 derivative operator D^ rather than derivative operator A,-,- provides a symmetric matrix representation of FP operator. It is essential for the stability in the time integration of the ODE system ^^. = -YtLii9{xi,t). (3.2.20) .7 = 1 3.3 Eigenfunction expansion The formal solution of the Fokker-Planck equation, Eq. (3.2.1), may be evaluated by eigenfunction expansion. If <j)n(x) are the eigenfunctions of operator L in Eq. (3.2.3), the eigenvalue problem of the FP equation is defined by L<f>n(x) = -en<j>n(x), (3.3.1) where en are the eigenvalues and <f>n(x) are the eigenfunctions. The linear time-dependent equation, Eq. (3.2.1), admits a solution of the form oo P(x,t) = P0(x) £ cne-^Vn(x), (3.3.2) where the coefficients cn are determined from the initial condition, P(x, 0), and are given by cn = J P(x, 0)(f)n(x)dx. (3.3.3) The eigenvalues in Eq. (3.3.1) and Eq. (3.3.2) are in units of reciprocal time. The ground state eigenfunction is the equilibrium distribution Po(x) with zero eigenvalue. The eigenvalues en are real and satisfy en > 0, for n > 0. In actual calculations, the solution P(x,t) given by Eq. (3.3.2) must be truncated at n = N. It is essential that the lower eigenvalues and eigenfunctions converge accurately and rapidly in order to achieve a fast accurate approximation of the solution. To obtain rapid convergence of the eigenvalues and eigenfunctions of the FP operator, the choice of the weight function or the basis functions plays an important role. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 73 Since no numerical approximation involved in temporal variable t, the accuracy of the solution (3.3.2) is as the same order as the numerical method applied to approximate the differential operator L. 3.4 Temporal discretization Numerical schemes for temporal discretization of the Fokker-Planck equations have been studied by several researchers [97-100]. Most of the methods used are finite difference methods. Chang and Cooper [97] used centered differences on the diffusion term and weighted differences on the drift term of the FP equation. Their method corresponds to a standard finite difference scheme for solving time-dependent FP equations and was used in many applications [154]. Larsen et al [98] then generalized the Chang-Cooper method to allow for more efficient solution of the nonlinear FP equation using larger time steps. More recently, Epperlein [99] developed a scheme which extended the standard Chang-Cooper scheme by providing not only number density conservation but also energy conservation and allowed much larger time steps for stable and accurate solutions. The schemes are implicit in the temporal discretization. In a review of numerical methods for solving FP equations of stochastic acceleration, Park et al [100] compared several finite difference schemes, including the Chang-Cooper method, the method by Larsen et al, and some semi-implicit schemes based on the above two schemes. They pointed out that the fully implicit Chang-Cooper method is the most robust finite difference method. Other methods suffer from instability and accuracy problems when dealing with some FP equations. Although there is some success of these methods, the finite difference methods generally encounter some difficulties that would be easy to solve or avoid in the QDM. First, the FP equation is defined over a infinite interval. To a,void this problem for the finite difference methods, one has to evaluate the equation over a finite interval Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 74 and impose boundary conditions on the newly defined boundaries. This implementation usually would cause some boundary effects in the numerical calculation. Second, the current finite difference methods available are at most 2nd order accurate in the spatial dimension, while the QDM is spectral accurate. Furthermore, since the differential matrix constructed in the finite difference schemes is usually not positive definite, an implicit or semi-implicit time discretization has to be performed to obtain a stable solution. With the QDM, the differential matrix of the FP operator is symmetric and the eigenvalues of the matrix are positive. This indicates that one can even use explicit time discretization to obtain a stable solution of the problem if the condition number of the matrix is not very large. For a implicit method, there is no restriction to the choice of time steps for a stable solution. Finally, the QDM differential matrix is very easy to construct and the method is very easy to implement. The finite difference schemes seem more complex to implement. At present, for simplicity, we only use the first order implicit and explicit Euler discretization in the time intergration. 3.5 Applications and results In this section, we apply the methodology of the previous sections to several different systems described with a FPE with bistable states. Our main objective is to determine the optimum basis sets with the QDM formalism for the evaluation of the eigenvalues and eigenfunctions of the Fokker-Planck operator. The time dependent solutions are then expressed in terms of an eigenfunction expansion. (1) The Quartic Potential: A very popular model system is the quartic potential, Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 75 defined by choosing, A(x) — x3 — x and B(x) = e with x € [—00, 00] . The equilibrium PDF, Po(x), exhibits a bimodal form indicative of two stable states that occur for x = ±1. This system has been considered by many authors in a study of the role of fluctuations in systems far from equilibrium and the subsequent evolution of such systems. Several workers have sought numerical and semianalytical solutions to this system. Eigenvalues for this system were reported by Dekker and van Kampen [146] with a finite difference approximation, by Risken [93] with a matrix continued fraction technique and Indira et al [145] with a finite difference scheme to solve the time dependent FPE. Blackmore and Shizgal [71] employed the QDM, with a quadrature based on the steady equilibrium PDF as the weight function. Since the equivalent Schrodinger equation (see section 4.2 for detail) has three wells at x = 0 as well as at x = ±1 (see Fig. 2 of reference [71]), in this work we use three different weight functions and compare the results for each. We choose (a) the equilibrium PDF wa(x) = Po(x) = exp[—(x4/4 — x2/2)/e], (b) a Gaussian weight function centered at the origin, Wb(x) = exp(—x2/2e), (c) a second narrower Gaussian weight function centered at the origin, wc(x) = exp(—4x2/e), and (d) the sum of the equilibrium PDF and a Gaussian weight function peaked at the origin, Wd(x) = PQ(X) + exp(—x2/2.37e). The motivation for these different choices becomes clear when the variation of the eigenfunctions are seen. The widths of the Gaussian functions were chosen by trial and error to optimize the rate of convergence. The eigenfunctions for e = 0.01 and e = 0.001 are plotted in Fig. 3.1. The dotted curves are for e = 0.01 and the solid curves are for e = 0.001. It is clear that the eigenfunctions are symmetric for n even and antisymmetric for n odd.- For e = 0.01, some of the eigenfunctions are localized either near the origin and/or near ±1, whereas others are spread over the domain. For e = 0.001, the eigenfunctions are generally concentrated either near x = ±1 and/or near x = 0. The eigenfunctions are concentrated in the Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 76 Figure 3.1: First ten eigenfunctions for the quartic potential model. The dotted curves are for e = 0.01 and the solid curves are for e = 0.001. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 77 regions of the minima in the potential of the Schrodinger equation [71]. For e = 0.01, the eigenfunctions are not as localized as for e = 0.001. The eigenfunctions are concentrated in the three regions corresponding to the extrema of the quartic potential, i.e., x = 0 and x = ±1, respectively. Therefore, for small e this suggests choosing a weight function Wd(x) with peaks at three extrema of the potential. This results in a set of quadrature points concentrated in these three regions where the eigenfunctions are localized. It is anticipated that this choice will provide a rapid convergence of the eigenvalues and eigenfunctions. In the previous work by Blackmore and Shizgal [71], only the weight function wa(x) was used. The convergence of the lower order eigenvalues with several different weight functions for e = 1.0, 0.01 and 0.001 is shown in Table 3.1. The last entry in each column is the eigenvalue converged to the figures shown. For e = 1.0, the eigenfunctions are not very localized and with the equilibrium PDF, wa(x), as weight function the rate of convergence of the eigenvalues is quite rapid . For the next set of results with e = 0.01, a Gaussian weight function centered at the origin is used and the convergence of the eigenvalues is slower than for e = 1.0. The convergence to 5 significant figures is achieved with 60 quadrature points. The first eigenvalue becomes extremely small with decreasing e and the convergence is very.slow and this behavior is discussed further with regard to the other applications. As e —> 0, the eigenvalues A3, A4 and A5 are triply degenerate [71] and tend towards 2. Similarly, the eigenvalues A7, A8 and A9 tend towards 4. Hence, with decreasing e the calculation of the splitting of these nearly degenerate eigenvalues becomes more difficult, as seen for A3 —> A5. Results for e = 0.01 with the equilibrium distribution as weight function give the same final converged values for the eigenvalues but the rate of convergence is different. The difference in the results for these weight functions is demonstrated further for e = 0.001. The results for e = 0.001 in Table 3.1 are reported for three different weight functions, Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 78 N Ai A2 A3 A4 A5 A6 A7 A8 A9 € = 1.0a 10 0.79216 3.5537 6 8623 11.065 15.910 23.238 30.185 46.074 15 0.79209 3.5489 6 8441 10.829 15.445 20.607 26.644 33.051 42.023 20 0.79209 3.5489 6 8436 10.827 15.415 20.543 26.161 32.291 38.846 30 6 8435 15.415 20.540 26.153 32.219 38.709 40 26.153 32.219 38.709 e = 0.016 10 0.96787 1.8659 2.6791 3.4024 4.0276 15 0.96786 1.8658 2.6776 3.3652 3.8922 20 2.6768 3.3328 3.3441 3.7107 30 1.7728 2.6765 2.9702 3.3321 3.5456 40 1.7206(-3) 1 8642 1.8708 2.6772 3.3050 3.3689 3.4526 50 6.9871(-6) 1 8645 1.8670 2.6772 3.3045 3.3655 3.4515 60 6.0771(-8)* 1.8670 3.3044 3.3655 3.4515 e = 0.001° 4 1.9960 2.0060 6 1.9880 1.9882 3.9999 4,0268 8 1.9879 1.9879 3.9529 3.9547 10 1.9879 3.9512 3.9513 15 3.9512 3.9512 e = 0.0016 4 0.99698 1 9880 6 1 9879 2.9727 3.9524 8 2.9726 3.9513 10 3.9512 e = 0.001d 6 1.0096 2 0220 2.0383 8 0.82578 1 9897 2.0024 3.6293 4.0967 10 1.0016 1 9912 1.9933 4.1423 4.1748 17.0192 18.3206 15 0.99953 1 9879 1.9879 1.9900 3.1069 3.9524 3.9546 4.3375 20 0.99698 1.9885 2.9723 3.9511 3.9512 4.0025 30 1.9879 2.9727 3.9512 3.9541 40 2.9726 3.9512 Table 3.1: Convergence of the eigenvalues for the quartic potential. aw(x) = Po[x); bw(x) = eT^IM; dw(x) = P0{x) + e-* K™^). PQ(X) = -*t) I * \s the equilibrium solution of the Fokker-Planck equation. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 79 the equilibrium distribution (wa), a Gaussian centered on the origin (wb), and the sum of the equilibrium and a Gaussian at the origin (wc). With the equilibrium distribution, only the eigenvalues A4, A5, A8 and A9, are calculated (see Fig. 3.1: solid curves), whereas for the Gaussian weight function centered at the origin, only the eigenvalues A2, A3, A6, A7 and A10, are calculated. In both cases, the convergence is very rapid. This demonstrates that the choice of weight function which determines the distribution of grid points is extremely important. The eigenfunctions <f>4, <TS5, (f>$ and <T59 in Fig. 3.1 have their greatest variation near x = ±1 and the equilibrium distribution provides the required distribution of grid points for the rapid convergence obtained. By contrast, the eigenfunctions </>2, </>3, (f>6 and 4>7 have their greatest variation near x = 0 and the Gaussian distribution centered at the origin provides the required distribution of grid points for the rapid convergence of these eigenfunctions. The eigenfunction pairs cf>3 and CT55, and </37 and <f>9, are antisymmetric (see Fig. 3.1) and the corresponding eigenvalues are (nearly) degenerate (see Table 3.1). The final entry in Table 3.1 is for e = 0.001 and a weight function which is the sum of the equilibrium distribution and a Gaussian weight function centered at the origin. In this case, all the eigenvalues are calculated and the convergence is moderately rapid, more so than with only the equilibrium distribution as weight function. A summary of the rate of convergence of the eigenvalues is presented in Table 3.2. Here, A^, Ni, Nc and Nd are the number of quadrature points associated with the four weight functions shown at the bottom of the table. The entries in the table are the approximate numbers of points required to calculate the eigenvalues shown to 5 significant figures. As can be seen in the results, the rate of convergence of the eigenvalues is extremely rapid for e = 1.0 with both weight function wa(x) and wc(x), for e = 0.01 with weight function iW(,(aj), and for e = 0.001 with weight function wj,(x). For e = 1.0, the convergence is twice as slow with the weight function Wb(x) because the Gaussian function employed is too narrow and there is an insufficient number of points in the Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 80 Ai A2 A3 A4 A5 Aio Al5 A20 A25 e = 1.0 A 0.79209 3.5489 6.8435 10.827 15.415 45.599 85.415 133.03 187.38 Na 15 15 30 20 20 40 40 50 60 Nb 50 70 50 70 70 90 110 140 160 Nc 20 20 20 15 20 30 50 40 50 € = 0.01 A 6.16(-12) 0.96786 1.8645 1.8658 1.8670 3.9435 5.9608 8.7931 12.269 Na 50 80 80 40 70 90 100 110 100 Nb 80 15 50 20 50 60 60 70 60 e = 0.001 A 1.2016(-109)(°) 0.99698 1.9879 1.9879 1.9879 4.9234 7.8013 9.6867 11.545 Na ** ** 8 8 ** 20 ** 30 Nb 4 6 ** ** 10 ** 20 ** Nd 20 15 15 30 40 40 50 60 (a)The value is estimated by Kramers' approximation. **More than 200 quadrature points are required for convergence. Table 3.2: Comparison of the rate of convergence of the eigenvalues for the quartic potential. Na,'Nb, Nc, and Nj, are the number of quadrature points for weight function Wa{X) = P0(x), Wb(x) = e-*2/(2£), Wc(x) = e-2/(0.25e) and Wd(xj = pQ(xj + e-x2/(2.37e)) respectively. interval where the eigenfunctions are large. The first eigenvalue decreases extremely rapidly with decreasing e and the value for e = 0.001 in Table 3.2 is determined with the Kramers' approximation; see Eqs. (3.5.2) and (3.5.3). Generally speaking, for large e, the eigenvalues and eigenfunctions can.converge very fast with a weight function equal to either the equilibrium distribution or a Gaussian cen tered near the origin (as long as a proper width is selected). For e = 0.01, the eigenvalues converge faster with the equilibrium function P0(x), i.e. wa(x), as weight function than with a Gaussian weight function Wb(x) shown in the Table 3.2. For e = 0.001, the rate of convergence improves dramatically for certain eigenvalues with wa(x) and for others with Wb(x). But both weight functions fail to provide accurate results for other eigenvalues. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 81 e Ai A2 A3 QDM 50 9.041867 30.45742 56.83964 100 13.04778 43.49627 81.02501 200 18.71484 61.94008 115.2331 300 23.06393 76.09384 141.4834 400 26.73058 88.02646 163.6139 Wu and Kapral 50 9.04 30.5 56.8 200 13.1 43.5 81.0 200 18.7 61.9 115 300 23.1 76.1 142 400 26.7 88.0 164 Table 3.3: Comparison of the eigenvalues for large e for the quartic potential. w{x) = P0{x) = exp[-(£ - f )/e], N = 20. By adding these two weight functions together all the eigenvalues can be calculated with a moderately rapid convergence. The first 25 eigenvalues can be calculated in this way to 5 significant figures with 60 points. For large e, the eigenvalues and eigenfunctions can converge very fast with either the equilibrium solution or a Gaussian centered in the middle (as long as a proper width is selected) as weight function. Wu and Kapral [149] employed this Fokker-Planck equation to study the barrier crossing dynamics in terms of the spectral properties of the Fokker-Planck operator. They studied the chemical rate problem in terms of the eigenvalue spectrum of the Fokker-Planck operator and a second projected operator defined in their paper. They were primarily interested in the very low barrier limit, that is for e >> 1. They expanded the eigenfunctions in a set of scaled Hermite polynomials and diagonalized the matrix representative of the Fokker-Planck operator in this basis set. They mention that the choice of scaling parameter is particularly important. Wu and Kapral reported the first Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 82 e yQDM ^Kramers ^Kramers _ yQDMyyQUM 0.01 6.1596(-12) 6.2518(-12) 0.0150 0.02 1.6232(-6) 1.6776(-6) 0.0335 0.04 8.0874(-4) 8.6901(-4) 0.0752 0.05 2.7761(-3) 3.0331(-3) 0.0926 0.06 6.3146(-3) 6.9792(-3) 0.1052 0.08 1.7774(-2) 1.9779(-2) 0.1128 0.1 3.3545(-2) 3.6951(-2) 0.1015 0.15 8.2136(-2) 8.5024(-2) 0.0352 0.2 1.3478(-1) 1.2897(-1) -0.0431 0.25 1.8717(-1) 1.6560(-1) -0.1152 0.3 2.3802(-l) 1.9564(-1) -0.1781 0.4 3.3407(-l) 2.4095(-l) -0.2787 0.45 3.7934(-l) 2.5828(-l) -0.3191 0.5 4.2294(-l) 2.7303(-l) -0.3544 Table 3.4: Variation of Ai versus e for the quartic potential: Comparison with the Kramers' estimate. three eigenvalues for several large values of e. These are shown in Table 3.3 in comparison with the QDM results with the equilibrium distribution as weight function and N = 20. For these large values of e, the QDM is extremely efficient and the results in Table 3.3 are converged to the number of significant figures shown. The more difficult numerical problem is for smaller values of e. The Wu and Kapral paper is but one of many models of chemical reactions viewed as the rate of passage of reactants in one well over the barrier to the second well which represents the product states. Usually, the initial distribution is taken to be the equilib rium distribution restricted to one well that is in the region —oo < x < 0. The solution of the Fokker-Planck equation is then given by Eq. (3.3.2) and the rate of change of the density of reactants as the integral of the PDF over half the domain. The rate becomes Figure 3.2: The variation of the lowest nonzero eigenvalue, Xi versus e for the quartic potential. The solid line is the converged result with the QDM. The dashed line is the Kramers' approximation, Eqs. (3.5.2) and (3.5.3). Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 84 purely exponential after an initial transient if the smallest non-zero eigenvalue is well separated from the higher eigenvalues, i.e. Aa << A2 < A3 The rate coefficient is then approximated by 1/Ai. A well known approximation for the reaction dynamics in a bistable potential, that is, for Ai is the result obtained by Kramers over fifty years ago [159]. This asymptotic result derived elsewhere [93,159] is given by, Aj « ^y/U"(0)\U''(l)\e-^-u^e. (3.5.2) With the explicit expression for U(x) for the quartic potential, we get that, Aa « -e'^4t. (3.5.3) 7T A comparison of the converged QDM results for Ai with the approximate result, Eq. (3.5.3), is shown in Table 3.4 and Figure 3.2. The graph demonstrates the very rapid decrease in AT with decreasing e consistent with the form given by Eq. (3.5.3). Departures from Kramers' result can be either positive or negative and become significant for e > 0.2. (2) Optical Bistability. Optical bistability refers to the phenomenon which can arise when light is transmitted in an optical cavity containing a medium. For particular conditions, the transmitted light intensity, It, is a nonlinear function of the incident light intensity, I;. The relationship between the incident and transmitted light intensities depends critically on the parameter C = ctL/2T where a is the linear absorption coefficient, L is the length of the cavity and T is the mirror transmissivity coefficient. When C exceeds a critical value, the It versus li is then discontinuous and this system is characterized by two stable states. Further detailed discussions of this phenomenon can be found elsewhere [150,151]. The behavior of this physical system can, within certain constraints, be modeled with Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 85 a one-dimensional FPE with coefficients A{x) = y — x 2C (3.5.4) 1 4- x 2 ' and 2qx2 (3.5.5) (1 + z2)2' The dimensionless input and output amplitudes are denoted by y and x £ [0, oo], re spectively. The parameter q = C/2NS where Ns is the saturation photon number. The stationary distribution PQ(X) can be bimodal for particular values of the parameters, in particular q and y. The variation of the potential in the equilibrium distribution, Eq. (3.1.2), is shown in Fig. 3.3 for q = 0.4, and different y values for which there exists either a single state potential or a bistable potential. Previous work by Blackmore et al [70] employed a weight function which approximated the bimodal PQ(X) as a sum of two Gaussians. In this work, we choose wa(x) = PQ(X) and a second weight function made up of two equilibrium distributions, that is, w\,(x) — ^o(^)+0.33Po(0.33a;+6.6). The quadrature weights and points are calculated as discussed in Sections 2.2 and 2.3. The first ten eigenfunctions for C = 8, q = 0.4, and y = 8 and y — 9 are shown in Fig. 3.4. The dotted curves are for y = 8.0 and the solid curves are for y = 9.0. The eigenfunctions for y = 9.0 in Fig. 3.4 are localized in two regions; one is near x = 0 and the other near x = 6.6. On the other hand, the equilibrium distribution (the n = 0 eigenfunction in Fig. 3.4) has only one peak near x = 6.6. If the equilibrium distribution were used as the weight function, the quadrature points generated would be concentrated in this region. The eigenfunctions for y = 9 are more localized in the two regions than for y = 8. A comparison of the numerical convergence of the eigenvalues, An, with weight func tions wa(x) and Wb(x) is shown in Table 3.5 for C = 8, q = 0.4 and several y values. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 86 Figure 3.3: The potential U(x) in the equilibrium distribution for the optical bistability problem, C = 8 and q = 0.4. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 87 Figure 3.4: The first ten eigenfunctions of the Fokker-Planck operator for the optical bistability problem with C = 8, q = 0.4. The dotted curves are for y = 8.0 and the solid curves are for y = 9.0. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 88 N Ai A2 A3 A4 A5 A6 A7 A8 A9 y = 7.5a 10 0.7779 3 4958 5 4983 6.7301 9 2306 13.2532 19.6935 29.7897 46.2350 20 0.5809 0 9297 1 3063 1.7289 2 1979 2.7500 3.7418 5.4299 5.7352 30 0 9291 1 3033 1.6991 2 1135 2.5445 2.9903 3.4510 3.9434 40 2.5440 2.9886 3.4455 3.9132 y = 8.0a 10 1.609(-5) 0.3773 0 7078 1.0355 1 4857 2 1358 3 7191 7 3766 16.4200 20 5.494(-7) 0.3771 0 6692 0.7342 0 9843 1 2403 1 5440 1 9160 2.3760 30 6.653(-7) 0 5255 0.6980 0 8822 1 0626 1 3013 1 5612 1.8527 40 4.642(-7) 0 5210 0.6976 0 8729 1 0468 1 2742 1 5206 1.7880 45 4.641(-7) 0 5209 0.6976 0 8728 1 0467 1 2738 1 5195 1.7860 50 1 2738 1 5195 1.7858 y = 9.0" 2 0.6606 4 0.6605 1.3178 1.9735 6 1.9718 2.6224 3.2768 8 2.6223 3.2693 3.9130 10 3.9126 y = 9.0b 6 0.1721 0.6613 1.3510 40.5944 8 0.04931 0.6606 1.3219 2.0198 5.0508 10 0.1482 0.6606 1.3172 1.3946 2.0090 3.3048 21.3499 20 0.2064 0.6605 1.3178 1.4933 1.9718 2.6223 2.9576 30 0.2063 1.3178 1.4849 1.9718 2.6223 2.9039 35 1.4849 2.8388 40 2.8383 2698 3 9180 3 2693 3 9126 3 2693 3 9126 Table 3.5: Convergence of the eigenvalues for the optical bistability model with q = 0.4. aw(x) = PQ(X). bw(x) = Po(x) + 0.33Po(0.33x + 6.6). Po(x) is the normalized equilibrium solution of the Fokker-Planck equation. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 89 As can be seen from the results, the rate of convergence of the eigenvalues is extremely rapid for the smaller y values with both weight functions. For y = 9, the convergence of some of the eigenvalues is extremely rapid with wa(x), whereas other eigenvalues are not calculated. With the weight function.Wb(x), we are able to calculate the first nine eigenvalues to 5 significant figures with 40 quadrature points. A comparison of the rate of convergence of several weight functions used for con structing the QDM matrix in the present and in the previous work is shown in Table 3.6. For the weight functions we used in our present work, wa(x) and Wb(x), the rate of convergence of the eigenvalues is more than twice as fast as that for weight functions wc(x) = x2exp(—x2) and Wd(x) = x2exp(—x2) + x2exp[—(x — 6.6)2/0.1)] used in the previous work by Blackmore et al [70]. For q = 0.4 and y = 7.5 and 8, the convergence is the same with both wa(x) and Wb(x). The so-called "speed" quadrature points with wc(x) gives poorer convergence. For y = 9, wa(x) is superior for some of the eigenvalues and worse for others. For example, for y = 9 and q = 0.4, the eigenvalues AT, A4 and A7 converge with 90 to 110 points, while eigenvalues A2, A3, A5, A6, A8 and A9 converge with less than 10 points. This is because that the eigenfunctions for y = 9 are concentrated in two regions; one is near x = 0, the other is near x = 6.6. However, the peak of the weight function wa(x) = Po(x) at x = 6.6 dominates and, with small N, most of the quadrature points are distributed near this region. Since there are not enough points or none in the x = 0 region, the eigenvalues corresponding to the eigenfunctions in this region are not calculated. With increasing N, there are more points near x = 0, and the "missing" eigenvalues are recovered. With Wb as weight function, the grid points are better distributed and all the eigenvalues converge quickly. The results for q = 1 are also shown in Table 3.6 and the results are similar to those obtained for q =0.4. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 90 Ai A2 A3 A4 A5 A6 A7 As A9 y = 7.5 q = 0.4 A 0.5809 0.9291 1.3033 1.6991 2.1135 2.5440 2.9886 3.4455 3.9132 Na 20 30 30 30 30 40 40 40 40 Nb 20 30 30 30 30 40 40 40 40 Nc 30 40 40 50 60 70 80 80 80 2/= 8.0 q = 0.4 A 4.642(-7) 0.3771 0.5209 0.6976 0.8728 1.0467 1.2738 1.5195 1.7858 Na 40 20 50 40 50 50 50 50 50 Nb 40 20 50 40 50 . 50 50 50 50 Nc 60 60 40 70 70 70 80 80 90 y =9.0 q = 0.4 A 0.2063 0.6605 1.3178 1.4849 1.9718 2.6223 2.8388 3.2693 3.9126 Na 90 4 4 100 6 6 110 8 10 Nb ' 30 20 20 30 20 20 40 30 30 Nc 30 100* 100 30 100* 100* 30 100* 100* Nd 30 60 70 40 90* 90* 50 90* 90* y = 7.5 q = 1.0 A 0.3037 0.6883 1.1146 1.5723 2.0519 2.5446 3.0443 3.5542' 4.0872 Na 20 20 20 20 30 30 30 30 30 Nb 20 20 20 20 30 30 30 30 30 Nc 30 30 40 40 40 40 50 50 60 y = 8.0 q = 1.0 A 1.514(-3) 0.3343 0.4922 0.7460 1.0443 1.3761 1.7367 2.1225 2.5324 Na 20 30 30 30 30 30 30 30 30 Nb 20 30 30 30 30 30 30 30 30 Nc 40 40 40 50 50 60 60 60 70 y = 9.0 q = 1.0 A 0.3166 0.6592 1.3099 1.8157 1.9516 2.5835 3.2045 3.2588 3.8136 Na 50 4 6 60 6 8 10 70 10 Nb 40 20 20 30 30 40 50 40 50 Nc 30 80 100 30 100 100 100 100 100 * Converge to essentially 3 figures. Table 3.6: Comparison of the rate of the convergence of the eigenvalues for the optical bistability model. awa(x) = P0(x); bwb(x) = Po(^) + 0.33Po(0.33x + 6.6); cwc(x) = x2e~x\ dwd(x) = x2e-* + x2e-(x-6.6)Vo.i_ Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 91 (3) Climate Models. Nicolis and Nicolis [154] presented a stochastic model for climatic transitions. The independent variable x in the Fokker-Planck equation represents the globally averaged surface temperature, T. The time dependent PDF of the temperature is then given by Eq. (3.1.5) with B(T) = q2/2 and A{T) = Q[l-a(T)]-toT\ (3.5.6) with l-a(r) = 7i,r<T1, (3.5.7) l-a(T) = 7o + /?r,T1 < x < T2; l-a(T) = 72,T>T2. The coefficient A(T) is the difference between the solar influx Q[l — a(T)] (Q being related to solar constant, taken to be 340Wm-2 with a the albedo) and the infrared cooling rate, ecrT4, (e being the emissivity and o the Stefan constant). In Eq. (3.5.7), 70, 71 and 72 are constants and ft is a temperature feedback coefficient. The values of the parameters used in this model are defined as Q = 340, e = 0.61, o =5.67xl0~8, T2 = 297.0, 7, = 61.2, 72 = 0.75, and for the continuity of A(x), 70 and T\ are defined by 70 — 72 — /?T2 and Ti - (71 - 7o)//3, respectively. For appropriate values of these parameters the system is bistable, and the correspond ing climate potential has two minima at T = T± and a maximum at T = T0 which lies between T+ and T_. It is important to mention that the defined A(T) is a nonsmooth piecewise continuous function. The equilibrium distribution can be written in the form: P0(T) = exp(--2U(T)), (3.5.8) Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 92 Ai A2 A3 A4 A5 Aio Al5 A20 A25 10 1.3990(-4) 0.57262 1.0695 1.1826 1.8395 20 3.6683(-5) 0.51358 0.65862 1.0430 1.1819 3.5430 7.6004 30 3.3949(-5) 0.49091 0.61516 0.94902 1.1703 2.9259 5.9582 10.019 15.110 40 3.3889(-5) 0.48983 0.61350 0.94177 1.1687 2.7497 5.2499 8.6076 12.757 50 3.3888(-5) 0.48981 0.61347 0.94160 1.1687 2.7335 5.0674 8.0389 11.684 60 0.48980 0.61347 0.94160 2.7331 5.0466 7.8278 11.045 70 5.0458 7.7945 10.865 80 7.7931 10.824 90 10.822 Table 3.7: Convergence of the eij ^envalues for the climate model with 8 = 0.007123 9 = 7.13. where U(T) is the climate potential defined by U(T) = — JT A(T')dT'. (3.5.9) Nicolis and Nicolis were particularly interested in the long time evolution of this system, namely, the time scale determined by the lowest non-zero eigenvalue of the Fokker-Planck operator. They employed a finite difference scheme proposed by Chang and Cooper [97] and integrated the FPE numerically to determine this long time depen dence. Figure 3.5 displays the graphs for steady state function, PQ(T), and the potential, U(T), for several values of 8 with q = 7.13 yr_1K2. With an increase in 8, the steady state moves from one peak at T+ to two peaks at T+ and T_, and then one peak at T_. In this way, there is a transition of the global average temperature from a warm to a cold state. The form of the potential changes from one well near T+ to two wells near T+ and T_, and then again to one well near T_. For convenience, we first rescale Eq. (3.5.6) with a linear transformation given by, T = 38.5x + 243.5 to position the two minima around x = ±1 and the maximum around x = 0. The results are studied versus the parameters 8 and q which determine the nature of the potential. The QDM is applied to this system with the equilibrium function PQ(X) as weight Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 93 Figure 3.5: (A) The potential U{T) and (B) the equilibrium distribution P0(T) for the climate model, q = 7.13 yr-1K2. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 94 function. In Fig. 3.6, we show the first ten eigenfunctions for q = 7.13yr_1 K2 and 0 = 0.007123, corresponding to the symmetric potential and symmetric stationary distri bution function in Fig. 3.5. The convergence of the eigenvalues corresponding to these eigenfunctions is shown in Table 3.7, and 25 eigenvalues can be calculated to 5 significant figures with 90 quadrature points. The behavior of this system is similar to the other bistable systems. The calculations of the eigenvalues converge very quickly for small barriers, that is, for q large. The convergence is slower as q decreases as was the case for the quartic potential for small e. The calculations also converge quickly for values of 0 which yield a potential with a single well at either cold or hot temperatures. As with the other bistable systems, the long time behavior of the system is determined by the first eigenvalue which can be very well separated from the higher eigenvalues. The variation of log10(Ai) versus 0 for several q values is shown in Fig. 3.7. The decrease in Ai as 0 increases is because of the increasing barrier between the two states until 0 = 0.007123 where Aj reaches a minimum. It then increases for increasing 0 as the barrier again decreases until there is only one well at high temperatures. The behavior shown in Fig. 3.7 is referred to as "critical slowing down" and is determined in the present work in a more direct approach than by Nicolis and Nicolis. The minimum value of Ai in Fig. 3.7 which occurs for 0 = 0.007123, decreases as q decreases and as the barrier between the two wells increases. The dashed curves in Fig. 3.7 are the estimates of Ai determine with Kramers' approximation, Eq. (3.5.2), modified to take into account the contribution from both potential minima. The approximation is surprisingly accurate for the range of q and 0 in Fig. 3.7. For all the bistable systems studied in this chapter, the time dependent probability density functions are given by Eq. (3.3.2). We show only one example in Fig. 8 of the probability density for the climate model (q = 7.13 and 0 =0.007123) with an initial Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 95 Figure 3.6: The first ten eigenfunctions of the Fokker-Planck operator for the climate model q = 7.13 yr^K2, and 0 = 0.007123. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 96 Figure 3.7: The variation of the lowest nonzero eigenvalue, Aa for the climate model. The solid line is the converged result with the QDM. The dashed line is the Kramers' approximation, Eq. (3.5.2). Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 97 delta-function distribution and reduced temperature i = 1.5 or T(0) = 301K. The dis tribution function initially broadens, with the peak moving to lower x and then there is a bifurcation with a second peak appearing at x = ±1. With increasing time there is a decrease in the peak at x = 1 and an increase in the peak at x = — 1. At equilibrium, the two peaks are equal in magnitude. The time dependent solution of the probability density function P(x, t) for the climate model for q = 7.13 yr_1K2, and 8 = 0.007123 with S function at x = 1.5 ( at high temperature T=401K) as initial condition is shown in Fig. 8. Implicit backward Eular scheme was used to calculate the time evolution. The One peak initial distribution is quickly transform into distribution with two peaks at x = ±1 (T=250K and T=282K) and stablized gradually. The time dependence of the global average temperature corresponding to the distri bution in Fig. 3.8 is shown by the solid curve in Fig. 3.9. The time scale is determined by Ai listed in Table 3.7. Because of the very large separation in the value of Ai relative to the higher eigenvalues listed in Table 3.7, the initial transient is not discernible in Fig. 3.9. There is a very rapid, almost instantaneous decrease in the temperature to about 280K and then the gradual slow decrease. The dashed curve is the single exponential decay with decay constant equal to Ai. Comparison between the QDM and other numerical methods of solution of Fokker-Planck equation are available in the papers by Chen et al [45] and Leung et al [47]. Some of the results are shown in Figure 3.10 and Figure 3.11. Detailed discussion can be found in the papers. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 98 t=0 2.00 Figure 3.8: Time dependence of the probability density function for the climate model for q = 7.13 yr_1K2, and 8 = 0.007123. The initial distribution is a delta function at T = 301 if corresponding to x = 1.5. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 99 Figure 3.9: The variation of the global average temperature corresponding to the distri bution in Fig. 7. The solid line is the result all terms in the eigenfunction expansion whereas the dashed line is the result with just two terms corresponding to A0 and Ai. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 100 Figure 3.10: Eigenfunctions of the Fokker-Planck operator for electron relaxation in Xe, E/n = 0.25 Td. The solid lines represent the exact results. The circle symbols are the results with the QDM and the plus signs are the results with a finite difference method. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 101 Figure 3.11: Variation of the maximum eigenvalue A^ versus the grid size, TV, for a finite difference approach and the QDM. Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation 102 3.6 Summary In the present chapter, we have demonstrated that accurate eigen-solutions of different Fokker-Planck equations can be determined with the QDM. The main theme of this work is with regard to the choice of the weight function that determines the distribution of quadrature points in the QDM. We have demonstrated that the convergence of the eigen functions and eigenvalues of the Fokker-Planck operator is rapid with the appropriate choice of weight function. For the bistable systems studied here, which include the quar tic potential, optical bistability and a climate model, the first non-zero eigenvalue can be extremely small and nearly degenerate with the zero eigenvalues that characterizes the equilibrium distribution. The small difference between these nearly identical eigenvalues can be determined with the QDM. The slow approach to equilibrium is then character ized by the smallest non-zero eigenvalue. It is anticipated that the gridding technique described in this work will find important applications to multidimensional problems as well as nonlinear problems which may include time dependent forcing terms and/or time dependent drift and diffusion coefficients in the Fokker-Planck equation. Chapter 4 Applications of the QDM to the Solution of Schrodinger equation 4.1 Introduction There have been increased interest in the solution of the quantum mechanical problems with a discretization of the wave function in the Schrodinger equation. The traditional methods for solving a Schrodinger equation usually involve the expansion of the wave function in a suitable basis set, and the reduction of the differential equation to a set of algebraic equations for the expansion coefficients. The discrete approach in the solu tion of other differential and/or integral equations has been used by researchers in other fields, notably neutron transport [95], radiative transport [96], and computational fluid dynamics [18,36]. In kinetic theory, a discretization method of solution of the linear Boltzmann equation was introduced by Shizgal [75]. For this problem, the velocity dis tribution function for atomic species for a model reactive system was evaluated at set of quadrature points based on a set of "speed" polynomials orthogonal with weight function x2exp(—x2) on the interval [0, oo]. There have been numerous papers on the solution of the elementary one-dimensional Schrodinger equation, with different methods and several choices for the potential function V(y) [101]- [131]. The details of the potentials in these references are discussed later with a rational given for those choices for benchmarking the QDM against other numerical methods. Some of 103 Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 104 the potentials studied include the nonpolynomial oscillator (NPO)potential of the form, V(y) = y2 + (4.1.2) i +gy This potential is chacterized by the parameters A and g. It is close to a harmonic oscilla tor except a deep and narrow anharmonic well near the origin for large A and g. In this case numerical calculation of the solution is extremely difficult. Mitra [101] employed Hermite polynomials as basis functions and reduced the Schrodinger equation to matrix form. Mitra obtained the eigenvalues and eigenfunctions by numerical diagonalization and reported numerical results for the first three eigenvalues. Kaushal [112] described a perturbative approach and compared with the previous numerical results. Bessis and Bessis [104] demonstrated that the matrix elements of the potential with Hermite basis functions can be done analytically and the numerical integrations by Mitra are unnec essary. Flessas [126] showed that for particular relationships between A and g there are some exact results for the eigenvalues of this potential. For example, if A = — Ag — 2g2, then E1 = l-2g and for A = -7g2-6g±g^25g2 - 12g + 4, then E2 = (9g + X)/g. These results are useful for benchmarking different numerical methods. Hautot [105] reconsid ered the calculation of the matrix elements of the Hamiltonian for this potential in the Hermite basis set. Lai and Lin [118] reported additional exact solutions not discovered by Flessas, and also introduced a Pade approximant analysis. Fack and Vanden Berghe [124] employed several different finite difference schemes to solve for the eigenvalues and eigenfunctions for this problem. They employed a fine grid of points and diagonalized matrices of dimensions 200 x 200. They compared their results with available numerical results of previous workers, as well as for models with known exact results. Varshni [128] and Witwit [121,110] extended the earlier work to a three-dimensional version of this potential. Scherrer et al [129] employed the continued fraction approach by Risken employed in the solution of the Fokker-Planck equation. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 105 We have also considered the potential given by, V{y) = y6-3y2 (4.1.3) considered by Sinha et al [131]. This potential belongs to the class of potentials that arise in supersymmetric quantum mechanics [114,115], and are the same class that results in the transformation of the Fokker-Planck equation to the Schrodinger equation [51,160]. These authors consider a comparison of the SWKB results [114,115] and an exact calcu lation of the eigenvalues from a direct integration of the Schrodinger equation. In this application, the ground state is known and can be used as the weight function.in the QDM. Kaluza [130] considered the anharmonic sextic oscillator defined by the potential, V(y)=1-y2 + 2y4+i-y6 (4.1.4) Kaluza employed an analytical Lanczos procedure to generate the tridiagonal matrix representative of the Hamiltonian for this potential. Since the algorithm is analogous to a Schmidt orthogonalization procedure, it suffers from considerable roundoff error. This problem was alleviated by using multiple precision arithmetic. Braun et al [125] employed a spectral method based on Chebyshev polynomials to study the same potential and was able to reproduce the numerical results of Kaluza and extend the precision of many of the higher eigenvalues up to 18 significant figures with 512 grid points. A fourth potential that we consider in this study is of the form, V(y) = y2 + ey4 (4.1.5) which has been studied by several workers. Banerjee et al [102] and Banerjee [103] em ployed a non-perturbative method with the product of scaled Gaussian and a polynomial as weight function to calculate the eigenvalues for this potential for various values of Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 106 e. Fernandez et al [106] and Arteca et al [108] applied a variational method to obtain the eigenvalues and compared with Banerjee's results. Fernandez and Castro [120] ob tained the eigenvalues of this potential by solving the corresponding Riccati equation with Pade approximants. Recently, Fernandez and Tipping [122] improved the solution of the Riccati equation for this potential with a separation of the eigenfunctions into odd and even parity. Fack and Vanden Berghe employed a finite difference method to solve this problem. Witwit [111,117] extended the work to two and three dimensional problems. We also include an application of the QDM to a two-dimensional Schrodinger equation with the Henon-Heiles potential that has been used by several researchers as a benchmark of different methods [161-168]. This model has also been used to study the chaotic behaviour [169-173]. In this first instance, we choose to apply the QDM with Hermite points and weights. We also employ a grid determined by other nonclassical weight function. The result are in very good agreement with those of other researchers. 4.2 The equivalence of the Fokker-Planck eigenvalue problem and the Schrodinger equation The eigenvalue problems of Fokker-Planck equation are equivalent to Schrodinger equa tions: Recall the eigenvalue problem of FPE presented in Eq. (3.3.1). If the independent variable, x, is transformed to a new variable, y, (4.2.1) and we define ipn(y) by My) = \/Po{x{y))Mx(y)), (4-2-2) Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 107 with Po(y) = yB(x(y))Po(x(y))i then the Fokker-Planck eigenvalue equation, Eq. (3.3.1), is transformed into a Schrodinger equation -^j^ + v(y)My) = ^n(y), (4-2.3) #V>u(y) = tnipn(y)-The potential function in the Schrodinger equation is derived from the drift and diffusion coefficients in the Fokker-Planck equation, that is VM = \w\y) - (4.2.4) where The potential functions obtained in this way belong to the class of potentials that occur in supersymmetric quantum mechanics [114,115]. If \f~P~o is differentiated twice, it is clear that it is the eigenfunction with a zero eigen value of the Schrodinger equation Eq. (4.2.3). The Schrodinger equation, Eq. (4.2.3), can in turn be transformed into a different Fokker-Planck equation. The equivalence of the Fokker-Planck equation with the Schrodinger equation has been discussed by several au thors [160]. It can be shown that the time-dependent Fokker-Planck equation equivalent to this stationary Schrodinger equation is given by dP(y,t) _ d dt dy with stationary solution given by W>{y)P(y,t)+d-rm (4.2.6) P0(y)=exP(- j W(y)dy). (4.2.7) The Fokker-Planck eigenvalue problem is given by u/( ,d<j)n d2cf)n Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 108 where the drift coefficient, A(y) = W(y), and the diffusion coefficient, B(y) = 1. The Fokker-Planck operator, L, is self-adjoint with the equilibrium function Po{y)-4.3 The QDM matrix representation of the Schrodinger equations The matrix representation of the Schrodinger equation, Eq. (4.2.3), has similar form to that of the Fokker-Planck equation. Consider the matrix elements of the Hamiltonian for a basis set {Sn(y)} orthogonal with unit weight function, that is, Hnm = -J Sn(y)S'm(y)dy + J Sn(y)V(y)Sm(y)dy. (4.3.1) With an integration by parts in the first integral, we have that Hnm = J S'n(y)S'm(y)dy + Vnm, (4.3.2) where Vnm = J Sn(y)V(y)Sm(y)dy. We now introduce a second polynomial set, {Fn} orthogonal with weight function w(y), that is, Sn(y) = y/^y)Fn{y), (4,3.3) where w(y) = exp(— / W(y')dy'). Eq. (4.3.2) can be rewritten as /in1 in' ™(y)[K + 7TF^Fn + ^ Fn\dy + Vnm. (4.3.4) If one of the cross terms in the integrand above is integrated by parts one gets that Hnm = J wF'nF'mdy + [Km - Vnm], (4.3.5) where V(y) = \w\y) -\w\y). (4.3.6) If the matrix representative Hnm is transformed back to the discrete representation in physical space with the transformation T, that is, N N Hij — } ] } ] TinHnmTmj, (4.3.7) n=0 m=0 Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 109 one finds that N Hi, = E DkiDkj + [V{yi) - V(yt))Si3. (4.3.8) If the potential of interest can be factorized in accordance with Eq. (4.2.4) then a possible choice of weight function would be given by the equilibrium distribution function or the square of the ground state wave function, Eq. (4.2.7). For this choice the term [V{yi) - V{Di)]Sij is zero since V(y) = V(y). The method can be applied to higher dimensional problems with product space of one-dimensional bases. For a two-dimensional Schrodinger equation, d2 d2 1 ~Q^~]py~ + V(X>y) ^nm(x,y) = Xnm^nm{x,y) (4.3.9) the eigenfunctions are represented by a two-dimensional grid constructed from the prod uct space of orthogonal polynomials in x and y. The matrix representative of the two-dimensional Hamiltonian for bases set {Xn(x), Ym(y)} orthogonal with unit weight func tion is given by 82 d2 = - JJ Xn,{x)Ym,{y)(— + — )Xn(x)Ym(y)dxdy + 11 Xnl(x)Ym,(y)V(x,y)Xn(x)Ym(y)dxdy (4.3.10) With an integration by parts, we obtain that -hfn'mi nm Sm'mJ X'n,(x)X'n(x)dx +Sn>n J Ym,(y)Ym(y)dy+ Vnlmi,nm, (4.3.11) where Vn'm\nm - If Xnl(x)Ym'{y)V(x,y)Xn(x)Ym{y)dxdy. As for the one-dimensional case, consider polynomial sets {Gn(x)}, {Em(y)} orthog onal with weight function u(x),v(y) respectively, that is, Xn(x) = ^u(x)Gn(x) (4.3.12) Ym(y) = yJv~ty~)Em(y) (4.3.13) Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 110 Eq. (4.3.11) can be rewritten as Hn'm',nm = ^m'm J u(x)G'n,(x)G'n(x)dx + $n'n J v(x)E'm,(y)E'm(lj)dlJ +[Vn ], (4.3.14) where J(\u\x)^-U'{x))u(x)G'nl{x)G'n{x)dx +*n'n /{\v2{y) - \v\y))v{y)E'm,{y)E'm{y)dy, (4.3.15) and U(x),V(y) satisfy U'(x) = -}nu(x), (4.3.16) V'(y) = -\nv(y). (4.3.17We obtain the QDM representation of the Hamiltonian by transforming matrix i/rn'n,m'm to the discrete representation Hijtki, that is, Hij,ki = hi E Dk'iDk'3 + Sij Dk'kDkn + [V(xi, yk) - V{xt, yk)]8ij8ki, (4.3.1.8) fc'=i fc'=i where V(x, y) = \u\x) - \u\x) + \v\y) - l-V\y). (4.3.19) and Nx, and are the numbers of quadrature points chosen in actual applications. 4.4 Applications and results The QDM has been applied to several one-dimensional Schrodinger equations by Shizgal and Chen [51]. In this work, the QDM was applied to several one-dimensional and a two dimensional Schrodinger equations. The eigenvalues and eigenfunctions were calculated Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 111 with different choice of weight functions and were compared with the results by other authors. The main purpose of this study is to consider the solution of the Schrodinger equation with the QDM, and to study the rate of the convergence of the eigenvalues versus the number of grid points (equivalently basis functions) for different weight functions. The basis functions, Fn(x), are orthonormal with respect to the weight function, w(x). Our interest is to try to suggest the weight function that provides optimal convergence of the eigenvalues. We consider four different one-dimensional potentials in the Schrodinger equation that have received considerable attention in the literature over the past decade. If the convergence for one-dimensional problems can be optimized, there would be a con siderable savings in computer time when applied to two and three dimensional problems. This has been demonstrated by Shizgal and Chen [51] in the application of the QDM to the two-dimensional Henon-Heles potential. The first potential that we have chosen and which has been studied extensively [101]-[129] is the NPO model, Eq. (4.1.2), shown in Fig. 4.1 as the solid curves. The dashed curves are the harmonic potentials, V(y) = y2 + X/g, for A = g = 100 and A = g = 10 (upper curve) and for A = 10 and g = 100 (lower curve); the potential departs from harmonic in the vicinity of the origin. The deep narrow anharmonic well near the origin gets deeper and narrower with increasing A and g. Many of the previous calculations have emphasized the calculation of the ground state eigenvalue for large g. For situations where the potential is close to harmonic, it would appear useful to use the scaled Hermite polynomials as basis functions based on the weight function, w\(y) = exp(—ay2), where a is a scaling parameter. For this NPO potential, we have carried out an extensive analysis of the behavior versus the two parameters g and A and for different weight functions. The results are summarized in Tables 4.1-4.10. In Tables 4.1-4.3 with g = 1, 10 and 100, we use the Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 112 Figure 4.1: The nonpolynomial oscillator (NPO) potential, V(y) = y2 + 1 ^gy2. A and g equal to (a) 10, 10 (b) 100, 100 and (c) 10, 100. The dash lines are the corre sponding harmonic potential V(y) — y2 + X/g. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 113 weight function for scaled Hermite polynomials and vary the scaling parameter a for each of the first 5 eigenvalues so as to get the value of a that yields the most rapid conver gence. The QDM is implemented as discussed in the previous papers [51] by constructing the orthogonal polynomials for the chosen weight function with the algorithm described by Gautschi [134]. The quadrature points are then determined [135] and the eigenvalues calculated from the numerical diagonalization of the QDM representative of the Hamil tonian, Eq. (4.3.8). The results are shown for A = 1, 10 and 100 in each table. In Table 4.1, we reproduce exactly (to 9 significant figures) the harmonic oscillator eigenvalues for A = 0. With increasing g, it is seen that the eigenvalues are getting increasingly equally spaced consistent with an harmonic potential. The underlined portion of each eigenvalue indicates the convergence to that number of significant figures. For g = 1, 10 and 100 in Tables 4.1-4.3, we get convergence of the eigenvalues with 25-45, 60-120, 170-180 quadrature points, respectively. The convergence is clearly much slower for the large values of g. The results in Table 4.3 for the largest eigenvalues are converged to no less that 3 significant figures. The slow convergence for large g is due to the narrow anharmonic form of the potential near the origin; see Fig. 4.1. For the results in Table 4.1, the values of a were chosen arbitrarily. The interest in this study is to develop techniques to optimize the convergence by selecting a weight function related in some way to the potential. Mitra [101] chose a = \A + A and Bessis and Bessis [104] suggested a = yl + A/(l + 0.5g). In Tables 4.4-4.6, we show the results analogous to those in Tables 4.1-4.3 using the value of a suggested by Bessis and Bessis. It is clear that the convergence in Tables 4.1-4.3 is faster than the convergence In Tables 4.4-4.6. We have extended the previous efforts by employing a weight function chosen empir ically but taking into account the form of the potential. Our previous experience [51] Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 114 N Ai A2 7.41750590 A = 100 26.70397902 26.70596477 26.70596566 26.70596563 26.70596563 A3 g=_0 5.00000000 g=l 5.62484751 5.58979112 5.58977905 5.58977894 5.58977893 5.58977893 g= 1 10.73364911 10.70106074 10.70102615 10.70102563 10.70102557 10.70102558 10.70102558 g = 1 41.44872496 41.44110330 41.44109963 41.44109976 41.44109975 41.44109975 A4 7.00000000 7.83801690 7.64836479 7.64820406 7.64820127 7.64820124 7.64820124 13.61840145 13.38898345 13.38834923 13.38832431 13.38832353 13.38832349 13.38832349 53.83672948 53.83909110 53.83909383 53.83909326 53.83909327 53.83909326 53.83909326 A5 9.00000000 10.32756366 9.68550819 9.68407574 9.68404264 9.68404202 9.68404202 16.71237123 15.82571826 15.81924074 15.81888806 15.81887214 15.81887152 15.81887149 15.81887149 64.45752724 64.18782502 64.18745791 64.18744198 64.18744105 64.18744100 64.18744100 A = 0 5 1.00000000 3.00000000 A = 1 10 1.23249101 3.51099389 20 1.23235080 3.50738872 25 1.23235072 3.50738837 30 1.23235072 3.50738835 35 3.507388340 A = 10 10 2.78256744 7.41859167 20 2.78233128 7.41750446 25 2.78233044 7.41750609 30 2.78233052 7.41750588 35 2.78233052 7.41750590 40 45 50 10 9.35966852 20 9.35941813 25 9.35941803 30 9.35941803 35 40 45 Table 4.1: Convergence of the eigenvalues with V(y) = y2 + ^ n • w\(y) = exp(—a!/2)5 where a is chosen for the fatest convergence. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 115 N Ai A2 A3 A4 As A = 1 g= 10 10 1.11770702 3.54168906 6.69900727 11.03978681 16.69475171 30 1.05932983 3.08883073 5.09038673 7.13612019 9.26980877 50 1.05929698 3.08809133 5.08285715 7.09048160 9.08892124 60 1.05929690 3.08809085 5.08284796 7.09037430 9.08805809 70 1.05929689 3.08809085 5.08284769 7.09037053 9.08801960 80 1.05929688 5.08284768 7.09037041 9.08801815 90 1.05929688 5.08284768 7.09037041 9.08801810 100 9.08801810 A = 10 g= 10 10 1.65877686 4.53929108 8.04585051 13.12551490 19.66357857 30 1.58013523 3.88195452 5.85711306 8.03082593 10.30455803 50 1.58002278 3.87904292 5.83286153 7,90413992 9.88876928 70 1.58002235 3.87903684 5.83276776 7.90315755 9.88233330 80 1.58002233 3.87903683 5.83276755 7.90315433 9.88230079 90 1.58002233 3.87903683 5.83276753 7.90315417 9.88229884 100 5.83276753 7.90315416 9.88229873 .110 7.90315416 9.88229873 A = 100 g= 10 10 5.82541635 12.16555870 15.97213490 22.13479362 29.83816388 50 5.79394465 11.57221790 13.62913696 15.99309324 17.99876164 70 5.79394241 11.57219684 13.62877371 15.98848089 17.97250413 90 5.79394231 11.57219677 13.62877143 15.98843454 17.97208972 100 5.79394230 11.57219678 13.62877142 15.98843423 17.97208598 110 5.79394230 11.57219678 13.62877142 15.98843421 17.97208565 ' 120 15.98843421 17.97208562 130 17.97208562 Table 4.2: Convergence of the eigenvalues with V(y) where cv is chosen for the fatest convergence. y2 + vh(y) = exp(-m/2 Chapter 4. Applications of the QDM to the Solution of Schrodinger equation N Ai A2 A3 A4 A5 A = 1 g= 100 10 1.74034331 6.61289499 14.45434406 25.73416848 40.57425497 30 1.04840358 3.33049421 6.20396393 9.96961039 14.75340147 50 1.01083691 3.04192292 5.19309999 7.64468189 10.57282857 100 1.00841233 3.00987806 5.00989049 7.01481472 9.03654508 150 1.00841061 3.00983181 5.00927636 7.00985617 9.00959190 170 1.00841060 3.00983177 5.00927557 7.00984578 9.00949511 180 1.00841060 3.00983177 5.00927553 7.00984517 9.00948856 A = 10 g = 100 10 2.12557689 8.03895659 17.66231586 31.55531547 49.82378424 50 1.09321568 3.19606870 5.54663516 8.42645412 11.98070336 100 1.08408954 3.09891916 5.09892856 7.13621223 9.24759634 150 1.08406343 3.09831922 5.09279892 7.09883559 9.09763231 170 1.08406335 3.09831722 5.09276616 7.09850083 9.09530236 180 1.08406335 3.09831706 5.09276332 7.09846755 9.09503285 A = 100 g = 100 10 2.92175390 9.34160581 19.45189904 34.28680807 53.64778578 50 1.84742726 4.11049745 6.47955464 9.57139193 13.28590654 100 1.83638157 3.98422018 5.93857806 8.04492347 10.17201242 150 1.83633621 3.98310435 5.92841712 7.98535022 9.95499695 170 1.83633594 3.98309903 5.92834037 7.98458485 9.95023642 180 1.83633590 3.98309857 5.92833282 7.98449794 9.94960676 Table 4.3: Convergence of the eigenvalues with V(y) = y2 + 1+^2. w\(y) = exp where a is chosen for the fatest convergence. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 1 N X1 A2 A = 1 3.50666367 3.50738781 3.50738831 3.50738835 3.50738835 A = 10 7.41837822 7.41750412 7.41750587 7.41750590 7.41750590 A = 100 26.70572641 26.70596964 26.70596558 26.70596563 26.70596563 A3 g= 1 5.59128149 5.58978006 5.58977901 5.58977894 5.58977893 5.58977893 g=l 10.73118613 10.70105591 10.70102563 10.70102558 10.70102558 g= 1 41.44628014 41.44114465 41.44110119 41.44109978 41.44109975 41.44109975 A4 7.64562064 7.64819920 7.64820110 7.64820123 7.64820124 7.64820124 13.60393371 13.38888589 13.38832411 13.38832349 13 38832349 53 91385775 53 84147035 53 83917850 53 83909702 53 83909346 53 83909328 53 83909327 53 83909326 53.83909326 A5 9.69097632 9.68404558 9.68404226 9.68404204 9.68404202 9.68404202 16.67001724 15.82484591 15.81888465 15.81887151 15.81887149 15 81887149 64 86511926 64 23043100 64 19022920 64 18763807 64 18745616 64 18744228 64 18744111 64 18744101 64 18744100 64.18744100 10 1.23272180 30 1.23235100 40 1.23235074 50 1.23235073 60 1.23235072 70 1.23235072 10 2.78258502 20 2.78233137 30 2.78233053 40 2.78233052 45 2.78233052 50 10 9.35945915 15 9.35941761 20 9.35941803 25 9.35941803 30 35 40 45 50 55 Table 4.4: Convergence of the eigenvalues with V(y) = y2 + 1 ^ 2 • Weight functi u>i{y) = exp(-yyi + T^)-Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 118 N Ai A2 As A4 As A = 1 g= 10 10 1.06515662 3.08692234 5.08700592 7.08838586 9.09234744 30 1.06003407 3.08794408 5.08337022 7.09012096 9.08856179 50 1.05946708 3.08805698 5.08296831 7.09031284 9.08814362 100 1.05930820 3.08808859 5.08285570 7.09036658 9.08802645 150 1.05929829 3.08809057 5.08284868 7.09036993 9.08801914 170 1.05929756 3.08809071 5.08284816 7.09037018 9.08801860 180 1.05929736 3.08809075 5.08284802 7.09037025 9.08801845 A = 10 g= 10 10 1.61407526 3.87252286 5.85545056 7.90700483 9.98909626 30 1.58268033 3.87852707 5.83441119 7.90231956 9.88401555 50 1.58046212 3.87895272 5.83303946 7.90301643 9.88258285 100 1.58003812 3.87903382 5.83277730 7.90314922 9.88230893 150 1.58002355 3.87903660 5.83276829 7.90315378 9.88229952 170 1.58002282 3.87903674 5.83276784 7.90315400 9.88229905 180 1.58002265 3.87903677 5.83276773 7.90315406 9.88229894 A = 100 g = io 10 5.89164179 11.65464995 14.22630311 17.92322840 22.47689280 30 5.79569188 11.57183960 13.62953798 15.99205694 17.99649871 50 5.79404301 11.57217532 13.62879829 15.98841237 17.97213192 100 5.79394280 11.57219667 13.62877155 15.98843409 17.97208577 130 5.79394234 11.57219677 13.62877143 15.98843420 17.97208563 140 5.79394232 11.57219677 13.62877142 15.98843420 17.97208562 150 5.79394231 11.57219677 13.62877142 15.98843421 17.97208562 160 5.79394230 11.57219678 15.98843421 170 5.79394230 11.57219678 Table 4.5: Convergence of the eigenvalues with V(y) = y2 + 1 _^gy2 • Weight function Wl(y) = exp(-y2yiT^J^). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 119 N Ai A3 10 1.00943021 30 1.00902768 50 1.00883028 100 1.00860665 150 1.00851699 170 1.00849611 180 1.00848760 A = 1 3.00981139 3.00981943 3.00982338 3.00982785 3.00982964 3.00983006 3.00983023 g = 100 5.00980562 5.00959630 5.00949367 5.00937742 5.00933081 5.00931996 5.00931554 7.00981397 7.00982620 7.00983220 7.00983900 7.00984172 7.00984236 7.00984261 9.00989938 9.00973614 9.00965609 9.00956543 9.00952908 9.00952062 9.00951717 10 1.09402833 30 1.08994543 50 1.08798805 100 1.08583151 150 1.08499665 170 1.08480631 180 1.08472929 A = 10 3.09811836 3.09819977 3.09823881 3.09828179 3.09829842 3.09830221 3.09830374 g = 100 5.09791649 5.09580087 5.09478843 5.09367433 5.09324340 5.09314518 5.09310544 7.09814755 7.09827100 7.09833030 7.09839559 7.09842085 7.09842660 7.09842893 9.09889054 9.09723424 9.09644473 9.09557603 9.09524005 9.09516348 9.09513250 10 1.92323022 30 1.87989757 50 1.86200758 100 1.84549528 150 1.84038615 170 1.83936741 180 1.83897462 A = 100 3.98167955 3.98225971 3.98260632 3.98292370 3.98302127 3.98304068 3.98304816 100 5.97436235 5.94937345 5.94066820 5.93271091 5.93026369 5.92977655 5.92958881 8.01099094 7.98317686 7.98370029 7.98417970 7.98432709 7.98435641 7.98436771 10.12548823 9.96554325 9.95876339 9.95257021 9.95066627 9.95028732 9.95014127 Table 4.6: Convergence of the eigenvalues with V(y) = y2 -f 1^y2 • Weight function ™1(y) = exp(-yyi + TTA-^-). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 120 N 5 10 15 20 25 30 00 (c) 5 10 15 20 25 30 1.39754248 1.23347542 1.23235072 1.23235072 1.23235072 1.23235353 2.78138892 2.78233156 2.78233052 2.78233052 A = 1 6.00256900 3.50334632 3.50738845 3.50738835 3.50738835 3.50738835 3.50739706 A = 10 8.72184392 7.41816173 7.41750593 7.41750590 7.41750590 g = 1 11.66691170 5.91911961 5.58977876 5.58977894 5.58977893 5.58977893 5.58977892 5.58983355 g = 1 14.67163572 10.81174060 10.70102881 10.70102558 27.12138709 7.73245910 7.64821025 7.64820124 7.64820124 7.64820121 7.64906899 29.92249451 13.48916964 13.38872711 13.38832349 10.70102558 13.38832349 52.46889173 13.42453838 9.68403519 9.68404205 9.68404202 9.68404202 9.68404195 55.15913549 18.62437460 15.82253275 15.81887215 15.81887149 15.81887149 2.78233054 7.41767206 10.70448059 13.39000325 5 10 15 20 25 30 (0 9.35957820 9.35941835 9.35941803 9.35941803 A = 100 26.76092127 26.70595968 26.70596563 26.70596563 g = 1 41.56662303 41.44117242 41.44109980 41.44109975 41.44109975 60.49260798 53.84491548 53.83909597 53.83909326 53.83909326 91.26562732 64.45670875 64.18766807 64.18744157 64.18744100 64.18744100 9.35941803 26.70596563 41.44109975 53.83909296 Table 4.7: Convergence of the eigenvalues with V(y) = y2 + 1^ 2 • Weight function W\(y) = exp(—Oiiy2)/(l + gy2)"2, where ot\ and a2 are given in Table 10. Results from (b) Fack and Vanden Berghe (1985) [124]; (c) Lai and Lin (1982) [118]. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 121 N Ai A2 A3 A4 A5 A = 1 g= 10 10 1.06078000 3.10750137 5.34002953 7.84826441 12.39506550 12 1.05934121 3.08925905 5.10521689 7.20306392 9.76442476 15 1.05929700 3.08810693 5.08309896 7.09455928 9.11338782 20 1.05929690 3.08809085 5.08284789 7.09037308 9.08806709 25 1.05929688 3.08809085 5.08284767 7.09037041 9.08801812 30 1.05929688 5.08284768 7.09037041 9.08801810 35 5.08284768 9.08801810 A = 10 g= 10 10 1.67530513 4.31996615 9.75230993 10.43741586 29.54884661 20 1.58002638 3.87916818 5.83491979 7.91869454 9.97591517 25 1.58002232 3.87903807 5.83278723 7.90346720 9.88457979 30 1.58002233 3.87903684 5.83276775 7.90315728 9.88233966 35 1.58002233 3.87903683 5.83276753 7.90315420 9.88229915 40 3.87903683 5.83276753 7.90315416 9.88229873 45 7.90315416 9.88229873 (b) 1.58002233 3.87903683 5.83276752 7.90315413 9.88229866 A = 100 g= 10 10 7.48981433 8.03655640 44.09323078 49.71269815 149.70360371 20 5.79394731 11.57682425 13.70481854 16.27467472 19.37491150 25 5.79394193 11.57225713 13.62958358 16.00289112 18.03283044 30 5.79394232 11.57219704 13.62878139 15.98861076 17.97411182 35 5.79394230 11.57219678 13.62877147 15.98843650 17.97211015 40 5.79394230 11.57219678 13.62877142 15.98843422 17.97208593 45 13.62877142 15.98843421 17.97208562 50 15.98843421 17.97208562 Table 4.8: Convergence of the eigenvalues with V(y) = y 2 + ^2 • Weight function wi(y) = exp(—aiy2)/(l + gy 2) 012, where ct\ and a2 are given in Table 10. (b) results from Fack and Vanden Berghe (1985) [124]. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 122 N Ai A2 A3 A4 A5 A = 1 g= 100 10 1.04129101 3.29082677 6.43712676 10.95519923 55.86316463 20 1.00841165 3.00986701 5.00988241 7.01509885 9.04372955 25 1.00841060 3.00983203 5.00928092 7.00992522 9.01020204 30 1.00841060 3.00983177 5.00927556 7.00984573 9.00949652 35 3.00983177 5.00927551 7.00984496 9.00948602 40 5.00927551 7.00984495 9.00948596 45 7.00984495 9.00948596 A = 10 g = 100 10 1.45426410 5.01654653 29.97472996 72.98922478 243.24844347 20 1.08416996 3.10035001 5.11504606 5.81559287 .7.21912972 30 1.08406336 3.09831782 5.09277943 7.09864757 9.09661478 40 1.08406334 3.09831700 5.09276191 7.09844919 9.09486638 45 1.08406334 3.09831700 5.09276189 7.09844907 9.09486470 50 5.09276190 7.09844907 9.09486466 55 5.09276190 9.09486466 A = 100 g = 100 10 .25273328 58.55206092 526.35484897 565.24148666 1785.60704400 20 1.89936536 4.43121455 8.05177942 12.63918879 27.70284628 30 1.83635795 3.98374475 5.93607533 8.03086862 10.15714058 40 1.83633587 3.98309869 5.92833557 7.98453363 9.95000108 50 1.83633584 3.98309834 5.92832858 7.98444358 9.94916197 60 1.83633583 3.98309834 5.92832857 7.98444352 9.94916096 65 1.83633583 5.92832857 7.98444352 9.94916096 0>) 1.83633444 3.98309836 5.92832790 7.90315413 9.88229866 Table 4.9: Convergence of the eigenvalues with V(y) = y2 + 1+YGY2 • Weight function wi(y) = exP(~~aiy2)/(l + 9y2)a21 where a\ and a2 are given in Table 10. (b) Results from Fack and Vanden Berghe (1985) [124]. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 123 g/A 1 10 100 1 (1,10) (1.2,10) (3,12) 10 (1.4,6) (2,8) (2,14) 10 (1.4,6) (2,8) (2,14) 100 (2,6) (2.4,8) (2.5,16) Table 4.10: (ai, a2) used for Tables 4.7-4.9. has suggested that a useful choice of weight function would be derived from the "super-potential" associated with the potential. This would require the solution of the Riccati equation [120] which is as difficult if not more so than the solution of the Schrodinger equation. However, we have also shown that this choice of weight function is not always the best choice. Nevertheless, we have used an empirical weight function of the form, w1(y)=exp{-a1y2)/(l+gy2)a2. (4.4.1) The results obtained with this weight function are shown in Tables 4.7-4.9. In Table 4.10, we list the values of ct\ and a2 in the weight function. For all pairs of A and g, we obtain convergence of the eigenvalues to 9-10 significant figures with no more than 60 quadrature points. It is useful to compare the convergence of A5 for A = 100 and g =100 in Table 4.9 and Table 4.3. In Table 4.3, A5 is converged to 9.950 with 180 quadrature points whereas it is converged to 9.94916096 with 50 points in Table 4.9. This demonstrates the usefulness of the QDM and the use of arbitrary weight functions to accelerate the convergence. This could mean a great decrease in computational times for two and three dimensional problems. In Table 4.11, we compare the present results for Ai with the results reported in the literature by other workers. The methods used by others have been summarized in the introduction to the thesis. The weight function used is of the form given by Eq. (4.4.1) Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 124 10 100 500 = 1 QDM 1.232350723406 2 782330515932 9.359418026324 21.65874769959 a 1.23235072 2 78233052 9.35941803 b 1.23235072 c 1.24213 d 1.23235353 2 78233054 9.35941803 e 1.23237205 2 782330 9.35941803 21.6587477 f 1.23235 2 78233 9.3594 g= 10 QDM 1.059296880862 1 580022327392 5.793942300193 16.73274738223 a 1.05929688 1 58002233 5.79394230 b 1 58002233 e 1.05929700 1 5800249 5.793947 16.73919 f 1.05929 1 58002 5.794 g = 100 QDM 1.008410597947 1 084063335494 1.836335833449 5.083683913501 a 1.00841060 1 08406334 1.83633583 b 1.83633444 c 1 08411 1.8411 e 1.0084106 1 0840543 1.8363850 5.0840857 f 1.00841 1 08406 1.8364 g = 500 QDM 1.001849154630 1 084063335494 1.18486023962 1.92317625551 a c 1.18451 1.92255 e 1.0018491 1 0184910 1.1848632 1.9232260 Table 4.11: Comparison of the results of Ai with V(y) = y2 + 1+^2 • a. Scherrer, Risken and Leiber (1988) [129], b. Fack and Vanden Berghe (1985) [124], c. Chaudhuri and Mukherjee (1983) [107], d. Lai and Lin (1982) [118], e. Bessis and Bessis (1980) [104], f. Mitra (1978) [101]. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 125 with values of c*i and a2 which are chosen empirically for different values of A and g. The QDM results shown in this table are converged to the significant figures shown; either 12 or 13. These converged results are considered to be accurate based on our experience in the calculation for analytic models. The most difficult parameter region is for g = 500 and as can be seen from the results in the table, we ha,ve achieved remarkable convergence with g = 500. The only other work to compare with are the results by Bessis and Bessis [104] and by Chaudhuri and Mukherjee [107]. The QDM results are far superior to the previous results. Figure 4.2 shows the variation of the error in Ai for the NPO model (g = A = 10) versus the number of quadrature points, n, for four different weight functions. Aexaci is defined as the eigenvalue converged to 14 significant figures calculated with the QDM. The fourth weight function (d) gives the most rapid convergence. The significant improvement with weight function (d) over the scaled Gaussian weight function can be explained with the variation of the eigenfunction shown in Figures 4.3 and 4.4. The convergence of the eigenvalue depends on the accurate determination of the eigenfunction near the origin. Figures 4.3 and 4.4 show the details of the ground state eigenfunction near the origin. The solid curves are determined with the new weight function, Eq. (4.4.1), and N = 140. This is considered to be very close to the exact result. The other results are obtained with N = 25. Figure 4.3 is for g = A = 100 for three different weight functions; Hermite polynomials (*), scaled Hermite polynomials (+) and the new weight function (o). Figure 4.3B shows the eigenfunction on a small scale near the origin. From the figure, we see that the points generated from the scaled Gaussian weight function can not describe the rapid variation of the eigenfunction near the origin. How ever, the new weight function, with a denser grid of quadrature points near the origin where the potential (Fig. 4.1) and the eigenfunction vary rapidly, is better. It is clear that the results with the new weight function, Eq. (4.4.1), gives the best convergence. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 126 0.00 -16.00 1 ' 1 1 1 ' 1 1 1 0 40 80 120 160 n Figure 4.2: Variation of the error in Ai, AAi = |Ai — Ajraci| for the NPO potential versus the number of grid points n for different weight functions, (a) w(y) = exp(—y2), (b) w(y) = exp(-j/yi + (c) w(y) = exp(-5.8j/2), (d) w(y) = ef^; A = g = 10 Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 127 Figure 4.3: Ground state eigenfunction for the NPO potential with g = X — 100, N = 25, with different weight functions. (*) w(y) — exp(—y2), (+) w(y) — exp( —17?y2) and (o) w(y) = e(i+g~y2)* • The solid curve is for the last weight function with N=140; (A) full scale; (B) small scale near the origin. Figure 4.4: Ground state eigenfunction for the NPO potential. A and g are equal to (a) 10, 10 (b) 100, 100 (c) 10, 100 and ( ) 0, 0. (A) full scale; (B) small scale near the origin; (C) small scale at large positions from the origin. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 129 Figure 4.4 shows the behavior near the origin for three different pairs of values of g and A. Figure 4.4B shows the small scale behavior near the origin whereas Fig. 4.4C shows the small scale behavior at large positions from the origin. The dashed curve is the result with the Gaussian weight function. If the potential belongs to the class of potentials in supersymmetric quantum mechan ics [114,115], then the ground state eigenfunction is known with the eigenvalue equal to zero. This is the case for the potential given by Eq. (4.1.3), considered by Sinha et al [131]. The weight function that corresponds to the superpotential is of the form, w2(y) = exp(-2/4/4) (4.4.2) The basis set was determined following the prescription by Gautschi [134] and the quadrature points as described in the earlier papers. For this choice of weight func tion, V(y) = V(y) and the representative of the Hamiltonian in the QDM representation is from Eq. (4.3,8) given by Hij = Y2k DkiDkj- We have studied the convergence of the eigenvalues for this potential with three different weight functions one of which corre sponds to the Hermite polynomials, W2(y) defined earlier with a = 5, and another given by, w2(y) = exp(-y4/4 - 5y2) (4.4.3) The results with the three weight functions are shown in Table 4.12. The overall conver gence is very similar with all three weight functions, although the second weight function Eq. (4.4.2) appears to give marginally faster convergence, in particular for the first eigen value. Our results are consistent with the results of Sinha et al [131] to the precision that they report in their paper. The third potential (Eq. (4.1.4)) chosen was studied by Braun et al [125] and Kaluza [130] and also belongs to the class of potentials in supersymmetric quantum mechanics. Kaluza chose basis functions such that the matrix representative of the Hamiltonian is Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 130 N Ai A3 A5 Aio w2(y) = exp(-5y2) 5 3.20578381 19.34421619 10 1.92166391 11.48428663 24.44773756 15 1.93541230 11.67877474 25.22960745 72 68872477 20 1.93548442 11.68098869 25.25435384 71 64137641 25 1.93548209 11.68097117 25.25461676 71 57368183 30 1.93548210 11.68097087 25.25460450 71 57923539 35 11.68097089 25.25460490 71 57902800 40 25.25460488 71 57903698 45 71 57903668 50 71 57903669 w2(y) = exp(-2/4/4) 5 1.95003306 13.51720225 10 1.93549705 11.68815652 25.58769695 15 1.93548226 11.68108903 25.26571988 75 81549114 20 1.93548210 11.68097109 25.25463882 72 04071624 25 11.68097089 25.25460546 71 58445530 30 25.25460488 71 57920993 35 71 57903737 40 71 57903670 45 71 57903669 w2(y) = exp(-y4/4 - by2) 5 4.54466778 23.73100470 10 2.23089971 13.39054786 31.59035642 15 1.94701006 11.78606709 25.89623300 83 65936104 20 1.93570651 11.68371077 25.28051990 73 07062002 25 1.93548392 11.68099725 25.25493023 71 64923337 30 1.93548212 11.68097108 25.25460771 71 57992969 35 1.93548210 11.68097089 25.25460489 71 57904422 40 11.68097089 25.25460488 71 57903671 45 71 57903669 Table 4.12: Convergence of the eigenvalues for SE with V(y) = y6 — 3y2. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 131 tridiagonal. The generation of the basis set is essentially a Gram-Schmidt orthogonal-ization which is subject to considerable round-off errors [75,143]. Kaluza avoids these numerical difficulties by using symbolic algebraic techniques in mathematics. For ar bitrary weight functions, this analytic approach is not feasible, whereas the Gautschi algorithm is generally convergent. Braun et al [125] employed a spectral method of solu tion based on Chebyshev polynomials on a finite interval where the cutoff at y = 8 is an additional parameter. They used up to 512 grid points and report the first 48 eigenvalues up to 18 significant figures. We have chosen the weight function, w3(y) = exp(-2y2 - y4/2) (4.4.4) and determined the matrix representative of the Hamiltonian in the "polynomial basis" representation, Eq. (4.3.5) with Vnm — Vnm = 2. The matrix elements of the Hamiltonian are determined with the quadrature define by the weight function, Eq. (4.4.4). Because of the symmetry of the potential, the eigenfunctions are of either even or odd parity. The matrix Hnm of dimension N X N can be decomposed into two matrices for the odd and even eigenfunctions each of dimension |x y. Since the matrix Hnm is pentadiagonal, the submatrices of even and odd parity are tridiagonal. The convergence of the eigenvalues from the numerical diagonalization of these tridiagonal matrices is rapid. The first 48 eigenvalues converge to 13 significant figures with less than 100 grid points and agree with the results by Braun et al. The fourth potential studied is given by Eq. (4.1.5). This potential is not in the class of potentials in supersymmetric quantum mechanics. We have in the first instance used scaled Hermite polynomials and the associated quadrature points to determine the eigenvalues with Eq. (4.3.8). The convergence of the lower order eigenvalues is shown in Table 4.14 for three values of e. The scaling is very important in order that the grid points are distributed over the region of y where the eigenfunctions are concentrated. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 132 N Ai A7 1 3 5 8 10 15 20 25 30 35 1.00000000000000000 1.00000000000000000 15.832389169799 15.124216267224 15.118931530866 15.118929992544 15.118929986242 15.118929986242 40.6232236023546 36.367167641896 36.343021051640 36.342716214160 36.342716212413 36.342716212413 66.261603950851 62.648395926012 62.356049424923 62.356028944861 62.356028944603 62.356028944604 62.356028944604 N ^20 A 30 MO A 48 30 40 50 55 60 65 70 75 80 85 90 95 100 310.4920471524 309.4993497820 309.4993484837 309.4993484837 848.8060217068 588.5806628599 566.4282265701 566.4026817440 566.4026355012 566.4026354734 566.4026354734 1346.579274312 947.4614543288 893.9968790569 872.0907745529 868.2562193165 868.1457422322 868.1452015357 868.1452006773 868.1452006767 868.1452006767 1597. 1364. 1248. 1183. 1149. 1138. 1137. 1137. 1137. 1137. 1137. 421054106 247596709 445773964 544197185 943901457 668487703 541785229 522672203 522588690 522588541 522588541 Table 4.13: Convergence of the eigenvalues of even parity with V(y) = \y2 + 2y4 + \y&. w3(y) = exp(-2y2-y4/2). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 133 With the notion that the optimal weight function should be the square of the ground state eigenfunction, we have fitted, to polynomials, the ground state eigenfunction deter mined previously with Hermite quadrature points. This is an alternative to solving the Riccati equation for the superpotential [106]. The fit is reasonably accurate but V(y) is not exactly equal to V(y). In Table 4.15, we show the convergence of the eigenvalues with this alternate weight function. The results with this weight function show a moderate improvement in the rate of convergence. We find for example that with the new weight function for e = 100, Ai is converged to 9 significant figures with 15 points whereas 20 points are required with scaled Hermite polynomials. Similarly, A3 is converged to 8 significant figures with 20 points whereas 25 are required with scaled Hermite polynomi als. The choice of weight function is clearly important for the rapid convergence of the eigenvalues. The QDM is also applied to a two-dimensional Schrodinger equation with Henon-Heiles potential, given by 1pnl(x,lj) = \nl^nl(x,y), (4.4,5) where V(x,y) = |(x2 + y2) + Xx(y2 — \x2) and A = A/0.0125, consistent with the choice of previous workers. Since the parameter A is not too large and the Hamiltonian is close to the two-dimensional harmonic oscillator, we choose, in the first instance, the discretization based on Hermite polynomials in both the x and the y directions. The QDM representation of the Hamiltonian given by Eq. (4.3.18) is easily constructed, with the potential function evaluated at the quadrature points in each dimension. The numerical diagonalization of the QDM matrix yields the eigenvalues and eigenfunctions. The convergence of the two basis sets in each variable is shown in Table 4.16. Nx and Ny are the numbers of quadrature points in x and y dimension, respectively. In the tables, the eigenvalues are labeled with the principle quantum number n and angular 1 d2 2d2 I ill 2d2y + V(x,y) Chapter 4. Applications of the QDM to the Solution of Schrodinger equation N Ai A3 As Aio € = 10(a) 12 2.44917485 16.63595545 35.88068953 94.30085478 15 2.44917408 16.63591955 35.88506209 95.81165911 20 2.44917407 16.63592150 35.88517148 96.15949348 25 2.44917407 16.63592149 35.88517122 96.15623411 30 16.63592149 35.88517122 96.15626312 35 96.15626298 40 96.15626298 e = 100<b) 10 4.99945382 34.87447875 75.72914876 253.32604009 12 4.99941563 34.87402295 75.88739267 201.40793502 15 4.99941758 34.87398862 75.87689375 205.27637088 20 4.99941755 34.87398427 75.87700463 204.79428957 25 4.99941755 34.87398426 75.87700403 204.79476335 30 34.87398426 75.87700403 204.79477459 35 204.79477451 40 204.79477451 e = lOOOOW 10 22.86146298 160.68335404 350.84170426 1022.19210882 12 22.86161867 160.68601691 350.38352262 924.84691394 15 22.86160889 160.68588347 350.43503532 944.02953926 20 22.86160887 160.68591272 350.43589703 947.71986787 25 22.86160887 160.68591261 350.43589621 947.68562278 30 160.68591261 350.43589622 947.68596392 35 350.43589622 947.68596166 40 947.68596167 45 947.68596167 Table 4.14: Convergence of the eigenvalues with V(y) = y2 + ey4. ^w4(y) = exp( ^w4{y) = exp(-l(h/2); ^wA{y) = exp(-6(h/2). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 135 N Ai no 10 12 15 20 25 30 35 2.44917318 2.44917406 2.44917407 2.44917407 e = lQ(a> 16.63603391 16.63593038 16.63592170 16.63592149 16.63592149 35.86694240 107.31938413 35.88380588 35.88516632 35.88517122 35.88517122 98.84824260 96.71828902 96.16096863 96.15625913 96.15626298 96.15626298 12 15 20 25 30 35 40 4.99941762 4.99941755 4.99941755 6 = lQQfc) 34.87397375 34.87398436 34.87398426 34.87398426 75.87733275 75.87701004 75.87700401 75.87700403 75.87700403 210.04422203 205.20313119 204.79819433 204.79477654 204.79477452 204.79477451 204.79477451 10 12 15 20 25 30 35 40 22.86160088 22.86160897 22.86160887 22.86160887 e = 10000<c) 160.68728162 160.68596913 160.68591446 160.68591261 160.68591261 350.26068143 350.42650209 350.43583997 350.43589612 350.43589622 350.43589622 1055.94778633 973.50076873 952.12677503 947.72238259 947.68593841 947.68596166 947.68596167 947.68596167 Table 4.15: Convergence of the eigenvalues with V(y) = y2 + ty4 calculated by fit ting weight function to ground state eigenfunction. ^w4(y) = exp( — (y4 + 5y2)/2); ^w4(y) = exp(-(2y4 + 6y2)); ^w4(y) = exp(-(50j/4 + 25y2)). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 136 momentum quantum number I, as discussed by Noid and Marcus [174]. The convergence 28 basis functions were required in each dimension. We cannot comment on the rate of convergence of other methods, since these are not reported by the earlier authors. The contour plots of several eigenfunctions are shown in Fig. 4.5 . The Czv symmetry of the eigenfunction is clearly evident and essentially converged with 32 new basis func tions. It is surprising that this small number of basis functions can resolve the details of the eigenfunctions as shown in Fig. 4.5. Since the basic idea of the QDM is to suggest nonclassical basis functions for which the rate of convergence is increased, we have also used the grid points based on the polynomials orthogonal with respect to the weight function, The grid points in y remain defined by the Hermite polynomials. We find some improvement in the rate of convergence, as shown in Table 4.17 for a selected number of eigenvalues. Since the potential is similar to a two-dimensional harmonic oscillator by virtue of the choice of the value A, it is not surprising that the Hermite basis set works as well as it does. The rates of convergences of the eigenvalue in the x and y dimensions are similar for both weight functions we applied. In Table 4.18, we compare, with the results of Feit et al [164], where the eigenvalues occur in almost degenerate pairs. Their results were calculated with 128 x 128 Fourier basis functions. Our results in this table for Hermite quadrature points converge to the significant figures shown with Nx = Ny = 50. With the new quadrature defined by the weight function given by Eq. (4.4.6), this convergence is achieved with Nx = Ny = 32. This demonstrated a significant saving of computational time in comparison with other researchers. As can be seen from of the lower-order eigenvalues is very rapid. Although the higher eigenvalues converge slower, the convergence is evidently still rapid. Even for the 80th eigenvalue , Ai2>2, only (4.4.6) Chapter 4. Applications of the QDM to the Solution of Schrodinger equation Ny/NX 4 6 8 10 12 16 20 24 28 A2,o (4th) 4 2.9692 2.9605 2.9601 2.9601 6 2.9632 2.9566 2.9563 2.9563 8 2.9631 2.9566 2.9563 2.9562 A5jl (16th) 4 7.0182 5.9903 5.9126 5.9092 5.9089 6 5.9505 5.8914 5.8677 5.8322 5.8256 8 5.8784 5.8311 5.8279 5.8242 5.8192 10 5.8716 5.8204 5.8179 5.8176 5.8176 12 5.8713 5.8198 5.8174 5.8171 5.8170 16 5.8712 5.8198 5.8174 5.8170 5.8170 AM (44th) 4 13.2657 12.6973 12.2964 6 11.0336 10.1869 10.1688 10.1686 10.1686 8 11.1240 9.6935 9.1708 9.1287 9.1284 9.1284 10 10.0448 9.0536 9.0347 9.0264 9.0241 9.0240 12 13.1807 9.9003 9.0499 9.0334 9.0247 9.0223 9.0222 16 12.4574 9.8905 9.0497 9.0330 9.0243 9.0218 9.0217 A9,_9 (54th) 4 15.4681 14.5918 14.5446 6 12.2319 12.1574 11.6740 11.5832 11.5749 8 11.2850 10.7096 10.6113 10.4985 10.4785 10.4774 10 12.4479 10.6346 10.1248 10.0684 10.0643 10.0640 10.0639 12 11.8715 10.5047 10.1240 10.0441 10.0379 10.0375 10.0375 16 15.5873 11.4840 10.3447 10.1237 10.0425 10.0361 10.0356 10.0355 20 14.9696 11.4788 10.3234 10.1237 10.0425 10.0360 10.0355 10.0354 Aio.io (65th) 4 17.8647 17.2548 17.0406 6 14.1919 13.5202 13.2560 13.1758 13.1734 8 12.3695 12.1856 11.8637 11.8527 11.8502 11.8489 10 12.3214 11.5307 11.2580 11.1959 11.1959 11.1699 11.1609 12 14.4481 11.9189 11.1145 11.0690 11.0547 11.0532 11.0531 11.0531 16 13.4192 11.6239 11.1072 11.0642 11.0515 11.0506 11.0502 11.0501 20 18.1812 13.2653 11.6110 11.1055 11.0637 11.0512 11.0499 11.0498 11.0497 24 17.7635 13.2631 11.6095 11.1051 11.0636 11.0512 11.0498 11.0497 11.0497 Ai2,2 (80th) 4 26.0197 21.9893 21.0597 6 16.2730 16.0380 15.8520 15.5594 8 17.2961 14.4322 13.7976 13.4736 13.4400 13.4335 10 17.7483 13.5335 13.0533 12.8019 12.6324 12.6245 12.6207 12 14.3616 12.8665 12.5890 12.3091 12.2442 12.2176 12.2143 16 16.1730 13.4822 12.4343 12.3128 12.1503 12.1477 12.1465 12.1463 20 25.3965 15.9680 13.3034 12.3997 12.2177 12.0693 12.0666 12.0655 12.0651 24 21.9153 15.6562 13.2932 12.3940 12.2139 12.0689 12.0660 12.0652 12.0651 28 21.3602 15.5597 13.2926 12.3932 12.2138 12.0689 12.0660 12.0651 12.0650 Table 4.16: Convergence of the eigenvalues, Ani/, for the Henon-Heiles potential wi Hermite points". "Weight function, u(x) = exp(—x2) and v(y) = exp(—y2). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 138 Figure 4.5: Contour plots of the eigenfunctions of the Schrodinger equation for the Henon-Heiles potential with n and / equal to (a) 2,0 (b) 6,0 (c) 9,-9 and (d) 10,6. new quadrature points were used with Nx = Ny — 32. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 139 Ny/NX 4 6 8 10 12 16 20 24 A2,o (4th) 4 2.9616 2.9601 6 2.9575 2.9563 8 2.9575 2.9562 A6,i (16th) 4 7.0082 5.9439 5.9089 5.9089 6 5.9531 5.8738 5.8283 5.8248 8 5.8413 4.9878 4.9865 4.9865 10 5.8764 5.8185 5.8176 5.8172 12 5.8761 5.8179 5.8171 5.8170 16 5.8761 5.8179 5.8170 5.8170 A8,8 (44th) 4 12.7920 12.2620 6 11.1456 10.1687 10.1686 10.1686 8 11.2597 9.7190 9.1286 9.1284 9.1284 10 10.0603 9.0446 9.0250 9.0241 9.0240 12 13.2.729 10.0034 9.0436 9.0232 9.0222 9.0222 16 12.5442 9.9984 9.0433 9.0227 9.0218 9.0217 A9 _9 (54th) 4 14.7114 6 12.2173 11.7764 11.5794 8 11.3010 10.6143 10.5439 10.4775 10 12.5191 10.6767 10.0691 10.0644 10.0639 12 11.9125 10.5473 10.0420 10.0379 10.0375 16 15.6786 11.5413 10.3703 10.0403 10.0361 10.0355 20 15.0393 11.5353 10.3499 10.0403 10.0360 10.0354 Aio.10 (65th) 4 17 1066 6 14.1440 13.2565 13 1733 8 12.3709 11.9700 11.8518 11 8489 10 12.3221 11.5196 11.1959 11.1959 11 1608 12 14.4451 11.9199 11.0714 11.0553 11.0531 11 0531 16 13.4728 11.7027 11.0702 11.0517 11.0504 11 0501 20 18.2794 13.3615 11.6918 11.0697 11.0514 11.0498 11 0497 A12,2 (80th) 4 23 8554 20.5669 6 16.0422 15 6786 15.2330 8 17.2781 14.3895 13.4600 13 4334 13.4323 10 17.9498 13.6380 12.9474 12.6293 12 6209 12.6180 12 14.4303 12.8445 12.3391 12.2266 12 2142 12.2141 16 16.2757 13.5128 12.4913 12.1806 12.1469 12 1463 12.1462 20 25.6830 15.9901 13.3842 12.4590 12.1705 12.0660 12 0651 12.0651 24 21.9864 15.7582 13.3790 12.4532 12.1701 12.0654 12 0651 12.0650 Table 4.17: Convergence of eigenvalues, An)/, for the Henon-Heiles potential with new points6. ^Weight function, u(x) = exp(—x2 4- §Ax3) and v(y) = exp(—y2). Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 140 the results in Table 4.18, we have obtained excellent agreement with the previous results. For n < 11 our results are smaller than those of Feit et al, whereas for n > 11 our results are slight larger. 4.5 Summary In this chapter, we have provided an extensive study of the use of the QDM in the solution of the Schrodinger equation for several one-dimensional and two-dimensional potential functions considered recently by several other researchers. The main theme of this work is to determine the optimal set of basis functions, equivalently the weight function, that provides rapid convergence of the eigenvalues versus the number of basis functions or grid points. Although this work is more restricted to one-dimensional problems, the extension to two and three dimensions is straightforward [51] as shown for the Henon-Heiles problem given in this chapter. The eigenvalues can be determined by the numerical diagonalization of the representative of the Hamiltonian in either the polynomial or the discrete basis. The work in this chapter generally employed the discretized version of the Hamiltonian at a set of points that correspond to the quadrature points associated with the chosen weight function. The distribution of grid points is determined by the weight function, which controls the convergence of the eigenvalues and eigenfunctions. We have demonstrated in this chapter the flexibility of the QDM in that arbitrary weight functions can be employed to improve the rate of convergence. In some cases, the improvement is remarkable such as for the nonpolynomial oscillator. Chapter 4. Applications of the QDM to the Solution of Schrodinger equation 141 n 1 Feit et al. QDMa QDM6 3 3 3.9825 3.982 417 ncc -3 3.9859 3.985 761 nc 5. 3 5.8672 5.867 015 nc -3 5.8816 5.881 446 nc 6 6 6.9991 6.998 932 nc -6 6.9996 6.999 387 nc 7 3 7.6979 7.697 721 nc -3 7.7371 7.736 885 nc 8 6 8.8116 8.811 327 nc -6 8.8154 8.815 188 nc 9 9 10.0356 10.035 413 nc -9 10.0359 10.035 592 nc 10 6 10.5727 10.572 480 nc -6 10.5907 10.590 470 nc 11 3 11.1603 11.160 258 11.160 259 -3 11.3253 11.325 231 nc 11 9 11.7497 11.749 519 nc -9 11.7525 11.752 297 nc 12 6 12.3335 12.333 785 12.333 786 -6 12.2771 12.277 192 nc 12 12 12.7474 12.748 445 12.748 520 -12 13.0310 13.032 062 13.032 065 13 3 13.0868 13.086 873 nc -3 13.0800 13.081 196 13.081 199 Table 4.18: Eigenvalues for the Henon-Heiles potential. aHermite qudrature points based on u(x) = exp(—x2) and v(y) = exp(—y2), with Nx = Ny = 50. 6New quadrature points based on u[x) = exp(—x2 + §A£3) and v(y) = exp(—y2), with Nx = Ny = 32. cnc indicates no change from the QDM° results. Chapter 5 Summary and future work In previous chapters, we have presented the quadrature discretization method (QDM) and have applied it to some model problems as well as a large class of Fokker-Planck equations and Schrodinger equations. It has been shown that the QDM is easy to imple ment. The method provides high accuracy and rapid convergence, and is very efficient in solving high dimensional PDEs. Unlike most spectral methods and FD methods, it is possible for the QDM discretization to preserve the symmetry for the self-joint joint dif ferential operator. Furthermore, by introducing Stieltjes's procedure, one can generate a set of nonclassical polynomials for any weight function, which is a significant improvment from previous work by Shizgal and Blackmore [71,72,74]. The flexibility of the QDM in choosing weight function provides one the opportunity to achieve the best solution with the least computing work. The QDM is a discretization procedure based on a grid of points that coincide with the quadrature points defined by a weight function over some interval. The major advan tage of the QDM over the classical spectral method is that the QDM allows one to use nonclassical polynomial basis set such that more choices a,re provided in the solution of differential equations. In the case that the basis functions are classical, the QDM is equiv alent to the classical pseudospectral method. By introducing the Stieltjes's procedure, it is easy for one to construct a set of nonclassical polynomials defined by an arbitrary weight function accurate^. For a specific problem, a specific weight function close to the solution, which is usually nonclassical, may be chosen to optimize the accuracy and 142 Chapter 5. Summary and future work 143 convergence of the solution. Therefore, the use of nonclassical basis sets may provide not only more choices but also superior accuracy and convergence than the classical spectral method for the solution of a given problem. The modified QDM (MQDM) derivative matrix is also introduced in the thesis. It is the Galerkin matrix of the derivative operator based on the nonclassical polynomial basis set. The MQDM has similar properties as the QDM. Furthermore, the MQDM matrix representative of a second order self-joint differential operator is symmetric. In the application to the model problems, such as the one-dimensional and three-dimensional Poisson equations, the QDM works remarkable well and competes favorably with the classical spectral method in both accuracy and the rate of convergence of the solution. In comparison with the finite difference method, the numerical solution of the QDM is more accurate and converges at a much faster rate. In the solution of the three dimensional Poisson problem, the QDM requires much smaller number of grid points and less CPU time than the FD method to converge to a given accuracy. In the solution of the model problem, an analytic time-dependent Fokker-Planck equations, the QDM with the equilibrium solution as weight function converges faster than the classical Hermite and Legendre methods. All the eigenvalues of the Fokker-Planck discretization matrix calculated by the QDM (and MQDM) are real. By contrast, spurious imaginary eigenvalues occur for the classical Hermite and Legendre methods. The condition number of the FP operator is less than or equal to 0(N2) for the QDM and classical methods. The main objective of the QDM is to determine the weight function that defines the polynomial basis sets and hence the grid points which provide optimum convergence in a given application. In the applications to the Fokker-Planck and Schrodinger equations, the matrix representative of the differential operator is very easy to construct and the Chapter 5. Summary and future work 144 matrix is symmetric. All the eigenvalues, except the zero eigenvalue, are real and nega tive. Spurious imaginary eigenvalues do not occur. The convergence of the eigenvalues and eigenfunctions is extremely rapid with a judicious choice of the weight function. The QDM is applied to three one-dimensional Fokker-Planck equations referred to as the optical bistable potential, the quartic potential and the climate potential, respectively. Comparisons have been made to investigate the influence of the choice of weight function or basis set on the rate of convergence of the eigenvalues with several weight functions. The result indicates that the choice of weight function (or the polynomial basis set) is important with respect to the rate of convergence of the eigenvalues and eigenfunctions. In general, to achieve fast convergence of the eigenvalues and eigenfunctions of the Fokker-Planck operator, the equilibrium solution can be chosen as weight function if the potential barriers are small. If the potential is symmetric, the weight function chosen as a Gaussian function centered on some symmetry axis with an appropriate width can be used. In the case that the potential barriers are very large, the weight function should be chosen such that the quadrature points that are generated are concentrated near the minima and maxima of the potential, or in the wells of the equivalent Schrodinger potential. The time dependent solution of the Fokker-Planck equation has been calculated for the climate model by both explicit and implicit time discretization. For implicit time discretization, since the eigenvalues of the differential matrix are negative, the solution is always stable and the selections of the time step is only controlled by the accuracy requirement. For explicit time discretization, the time step is decided by both stability and accuracy requirements. For a stable solution, one need to consider the condition number of the differential matrix. The larger the condition number, the smaller the time step can be chosen. We find that for all three Fokker-Planck equations we have studied, the condition number of the differential matrix constructed with equilibrium solution as weight function, grows as 0(NP), with 1 < p < 2, as N —> co. Chapter 5. Summary and future work 145 In the Fokker-Planck equation for the climate model, the drift coefficient function is piecewise linear but not smooth. By applying the QDM, we obtain very good convergence of the eigenvalues and eigenfunctions. This indicates that the QDM can work well for nonsmooth functions. The QDM is also applied to the one-dimensional Schrodinger equations with several potentials. The QDM in the solution of the Schrodinger equation is closely related to the methods for the solution of the Fokker-Planck equation. The matrix representative of the Hamiltonian is symmetric and all the eigenvalues are real. Different choices of weight functions in the QDM solution of the Schrodinger equations are compared and studied. The ground state wave function (if it is known) or some approximate form is chosen as weight function in the QDM to solve the problem. The convergence is very rapid in comparison with the classical Hermite method. For the Schrodinger equation with the nonpolynomial oscillator potential, the Hermite method converges very slowly in some cases, while the QDM with appropriate weight function converges extremely rapidly. In the application to the two dimensional Schrodinger equation with the Henon-Heiles potential, very rapid convergence of the eigenvalues and eigenfunctions is obtained with selected weight functions. The basis sets are constructed with different weight functions in each dimension, so as to optimaize the convergence of the eigenvalues and the resolution of the eigenfunctions. The results are in very good agreement with previous results by others. In general, the QDM is very efficient to the solution of the differential equations discussed in this thesis. The QDM also includes most classical spectral methods when one choose to use classical weight function (or basis set). Ability to use nonclassical basis set in the QDM gives one more flexibility in the solution of differential equations. With the QDM, it is possible to obtain the superior convergence and accuracy of the solution of some problems, for which any classical spectral method fails to deliver. The QDM can Chapter 5. Summary and future work 146 be easily extended to the solution of high dimensional differential equations. The present work has concentrated on the basic development of the QDM and its application to some model problems as well as Fokker-Planck and Schrodinger equations. For a complete understanding, improvement and broader applications of the QDM, some future work will be considered. (i) Adaptive selection of the weight function to optimize the solution. In general we may consider a function close to the solution as a weight function. In the case that no clue about the solution is given, one may choose an arbitrary weight function to generate a solution and then use this solution as a new weight function to solve the solution again. (ii) The nonclassical polynomials are generated by the Stieltjes's procedure. The accuracy of polynomials and hence the points and weights as well as the derivative matrix depends on the accuracy of numerical calculation of the recurrence coefficients by the Stieltjes's procedure. In some cases, the Stieltjes's procedure may not be accurate. Therefore other possibilities may be considered and studied to improve the accuracy of calculation of the nonclassical polynomials. (iii) The current QDM is based on the nonclassical polynomial basis set. To extend the class of the nonclassical basis set, a possibility of other nonpolynomial basis sets, for example the rational polynomial basis set, etc. may be considered. (iv) Both Fokker-Planck and Schrodinger equations are parabolic type in the infinite interval. They have similar forms and properties as advection-diffusion equations used extensively for the fluid dynamic problems. The techniques used in solving Fokker-Planck and Schrodinger equations may apply to this class of problems. Therefore studies of the QDM in the solution of the advection-diffusion equation is our main task in the near future. (v) All the problems we discussed in this thesis are linear differential equations. Some studies in spectral method in the solution of nonlinear problem (for example, the Burger's Chapter 5. Summary and future work 147 equation) has been carried out. Our interest in the future is to study how to apply the QDM to the solution of the nonlinear differential equations. (vi) When compared with the FD method, the spectral method is superior in accuracy and convergence of rate. However finite difference is the most popular numerical schemes used by scientists and engineers. One of the reasons is that it is easy to use and standard software routines are usually available and ready to use. In the future, we may consider to provide a general guideline for the choice of the weight function for a given class of problems and to develop standard software routines in order that one can apply the QDM easily. Bibliography [1] C.W. Gear, Numerical Initial value problems in ordinary differential equations. (Prentice-Hall, New Jersy, 1971). [2] J.D. Lambert, Computational methods in ordinary differential equations. (John Wiley k Sons, New York, 1973). [3] L.F. Shampine and M.K. Gordon, Computer solution of ordinary differential, equa tions: the initial value problem. (Freeman, San Francisco, 1976). [4] K. Dekker, Stability of Runge-Kutta methods for stiff nonlinear differential equa tions. (Elsevier Science, New York, 1984). [5] J.C. Butcher, The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods. (John Wiley & Sons, New York, 1987). [6] E. Hairer, The numerical solution of differential-algebraic systems by Runge-Kutta methods. (Springer-Verlag, 1989). [7] D. Gottlieb, M. Gunzburger, and E. Turkel, SIAM J. Numer. Anal. 19(4), 671 (1982). [8] R.E. Ewing, Compuf. Meth. in Appl. Mech. Engrg. 47, 74 (1984). [9] R.E. Ewing and R.F. Heinemann, Comput. Meth. in Appl. Mech. Engrg. 47, 161 (1984). [10] V. Girailt and 0. Raviart, Finite element methods for Navior-Stokes equations. (Springer, Berlin, 1896) [11] T. Peterson, SIAM J. numer. Anal. 28(1), 133 (1991). [12] R. Rannacher, Aamm Z. angew. Math. Mech. 73(9), 203 (1993). [13] P.P. Chavez, A two-dimensional finite element analysis of the stationary semicon ductor device equations, [microform]. Thesis. University of British Columbia. (1997) [14] M.A. Crisfield, Non-linear finite element analysis of solids and structures. (Wiley, New York, 1997) 148 Bibliography 149 [15] J. Bonet, Nonlinear continuum mechanics for finite element analysis. (Cambridge, New York, 1997). [16] V. Thomee, Galerkin finite element methods for parabolic problems. (Springer, Berlin, 1997) [17] L.T. Tenek, Finite element analysis for composite structures. Solid mechanics and its applications; v. 59. (Kluwer Academic Publishers, Boston, 1998). [18] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral methods in fluid dynamics. (Spring-Verlag, Berlin, 1987). [19] B. Fornberg, A practical guide to pseudospectral methods. (Cambridge University Press, Cambridge, 1996). [20] M.Y. Hussaini, M.D. Salas and T.A. Zang, Advances in computational transonics. (Pineridge Press, Swansea, 1983). [21] E. Tadmor, SIAM J. Numer. Anal. 26(1), 30 (1989). [22] J.M. Augenbaum, Appl. Numer. Math. 5, 459 (1989). [23] A. Noullez and M. Vergassola, J. Sci. Comput. 9(3), 259-281 (1994). [24] B. Fornberg, SIAM J. Numer. Anal. 27(4), 904 (1990). [25] A. Orszag, J. Comput. Phys. 37, 70 (1980). [26] K.Z. Korczak and A.T. Pattera, J. Comput. Phys. 62, 361 (1980). [27] M.G. Macaraeg and CL. Streett, Appl. Numer. Math. 2, 95 (1986). [28] D. Funaro, Numer. Math. 52, 329 (1988). [29] J.P. Boyd, Chebychev and Fourier Spectral Methods. Lecture Notes in engineering. (Springer, Berlin, 1989). [30] W. Heinrichs, Appl. Math. 37(6), 401 (1992). [31] M.O. Deville and E. H. Mund, J. Comp. Phys. 95, 359 (1991). [32] E.T. Olsen, and J. Douglas, Numer. Math. 69, 333 (1995). [33] T.D. Taylor, R.B. Myers and J.H. Albert, Comput. Fluids 9(4), 469 (1981). [34] C. L. Streett, AIAA-83-1949-CP: Proceedings of the sixth AIAA computational Fluid Dynamics conference, Danvers, Massachusatts, USA (1983). Bibliography 150 [35] A. Orszag and L.C. Kell, J. Fluid. Mech. 96, 159 (1980). [36] D. Gottlieb and S. Orszag, Numerical analysis of spectral methods: theory and applications. (SIAM, Philadephia, 1977). [37] A. Zebib, J. Comput. Phys. 53, 443 (1984). [38] R.G. Voigt, D. Gottlieb, and M.Y. Hussaini, Spectral method for partial differential equations. (SIAM, Philadelphia, 1984). [39] N. Funaro, Polynomial approximation of differential equations. Lecture notes in physics. (Springer-Verlag, New York, 1992). [40] G. Mansell, W. Merryfield, B. Shizgal and U. Weinert, Comput. Meth. Appl. Mech. Engrg. 104, 295 (1993); H. H. Yang, B. R. Seymour and B. D. Shizgal, Comp. Fluids 23, 829 (1994); H. H. Yang and B. Shizgal, Comput. Meth. Appl. Mech. Engrg. 118, 47 (1994). [41] I.M. Kuria and P.E. Raad, Comput. Meth. Appl. Mech. Engrg. 120, 163 (1995). [42] B. Machenhauer, Spectral method in Numerical Methods on Atmospheric Models, 1 (seminar proceeding), European Center for Medium Range Weather Forecast, Reading, UK (1991). [43] B.D. Shizgal and H. Chen, J. Chem. Phys. 107(19), 8051 (1997). [44] R.E. Robson, Aust. J. Phys. 46, 465 (1993). [45] H. Chen, K. Leung and B. D. Shizgal, Proceeding of the 20th international Sym posium on Rarefied Gas Dynamics, 273 (1996). [46] K. Leung, Electron and Ion relaxation in atomic and molecular moderators. PhD. thesis. (University of British Columbia, Vancouver, 1997). [47] K. Leung, B.D. Shizgal and H. Chen, J. Math. Chem. In press. [48] R. Kosloff and H. Tal-Ezer, J. Chem. Phys. 81, 3967 (1984). [49] M. Berman, R. Kosloff and H. Tal-Ezer, J. Phys. A 25, 1283 (1992). [50] W. Zhu, Y. Huang and D.J. Kouri, Chem. Phys. left. 217(1), 73 (1994). [51] B.D. Shizgal and H. Chen, J. Chem. Phys. 104(11), 4137 (1996). [52] K.S. Breuer and R.M. Everson, J. Comput. Physics 99, 56 (1992). [53] L.N. Trefethen and M.R. Trummer, SIAM J. Numer. Anal. 24(5), 1008 (1987). Bibliography 151 [54] J.A.C. Weideman and L.N. Trefethen, SIAM J. Numer. Anal. 25(6), 1279 (1988). [55] W.S. Don and A. Solomonoff, SIAM J. Sci. Comput. 16(6), 1253 (1995). [56] A. Bayliss, A. Class and B.J. Matkowsky, J. Compt. Phys. 116, 380 (1994). [57] W. Heinrichs, Math. Comput. 53, 103 (1989). [58] W. Heinrichs, Numer. Math. 56, 25 (1989). [59] E. Rothman, Proc. 2nd symposium on high-performance computing (Durand and Dabaghi, eds), Montpellier, France, 7 (1991). [60] D. Kosloff and H. Tal-Ezer, J. Comput Phys. 104, 457 (1993). [61] T.A. Zang, Y.S. Wong and M.Y. Hussaini, J. Compt. Phys. 54, 489 (1984). [62] D. Funaro and W. Heinrichs, Numer. Math. 58, 399 (1990). [63] M.O. Deville and E.H. Mund, SIAM J. Stat. Comput. 13(2), 596 (1992). [64] C. Canuto and A. Quarteroni, J. Comput. Phys. 60, 315 (1985). [65] T. Phillips, T. Zang and M.Y. Hussaini, IAM J. Numer. Analysis 6, 273 (1986). [66] D. Funaro, SIAM J. Numer. Anal. 24(5), 1024 (1987). [67] T. Tang, SIAM J. Sci. Comput. 14(3), 594 (1993). [68] D. Funaro and 0. Kavian, Math. Comput. 57(196), 597 (1991). [69] V. Iranzo and A. Falques, Comput. Meth. in Appl. Mech. Engrg. 98, 105 (1992). [70] R. Blackmore, U. Weinert and B. Shizgal, Trans. Theory and Stat. Phys. 15, 181 (1986). [71] R. Blackmore and B.D. Shizgal, Phys. Rev. A3l(3), 1855 (1985). [72] B. Shizgal and R. Blackmore, J. Comput. Phys. 55, 313 (1984). [73] J.A.C. Weideman, Numer. Math. 61, 409 (1992). [74] B. Shizgal, J. Chem. Phys. 70(4), 1948 (1978). [75] B. Shizgal, J. Computat. Phys. 41, 309 (1981) [76] A.D. Fokker, Ann. Physik 43, 810 (1914). Bibliography 152 [77] M. Planck. Sitzber. Akad. Wiss. 324 (1917). [78] A. N. Malakhov and A. L. Pankratov, Physica A 229, 109 (1996). [79] T. Blum and A. J. McKane, J. Phys. A: Math. Gen. 29, 1859 (1996). [80] H. Risken and T. Leiber, Phys. Rev. A 40, 1582 (1989). [81] J. B. Avalos and I. Pagonabarraga, Phys. Rev. E 52, 5881 (1995). [82] Y. Jia and J. Li, Phys. Rev. E 53, 5764 (1996). [83] B. Nowakowski, Phys. Rev. A 53, 2964 (1996). [84] S. Stepanow, Phys. Rev. E 54, 2209 (1996). A. Perico, R. Pratolongo, K. F. Freed and A. Szabo, J. Chem. Phys. 101, 2554 (1994); R. Pratolongo, A. Perico, K. F. Freed and A. Szabo, J. Chem. Phys. 102, 4683 (1995); G. J. Moro, J. Chem. Phys. 103, 7514 (1995) [85] A. P. Blokhin and M. F. Gelin, Mol. Phys. 87, 455 (1966); Physica A 229, 501 (1996). [86] T. Theuns, Month. Not. R. Astron. Soc. 279, 827 (1996). [87] O. Lies-Svendsen and M. H. Rees, J. Geophys. Res. 101, 2415 (1996). [88] H. G. Jhang and C. S. Chang, Phys. Plasmas 2, 3917 (1995). [89] P. Wilmott, S. Howison and J. Dewynne The Mathematics of Financial Derivatives (Cambridge University Press, Cambridge, 1995). [90] Y.Abe, S. Ayik, P.-G. Reinhard and E. Suraud, Phys. Reports 275, 49 (1996). [91] M. Schulz, Physica Scripta 37, 632 (1988). [92] N. G. van Kampen, Stochastic Processes in Chemistry and Physics (North Holland, Amsterdam, 1981). [93] H. Risken, The Fokker-Planck Equation, 2nd ed. (Springer, Berlin, 1984). [94] C. W. Gardiner, Stochastic Methods, (Springer, Berlin, 1983) [95] J. J. Duderstadt and W. R. Martin, Transport Theory (Wiley, New York, 1979); K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, Reading, Mass., 1967); B. Davison, Neutron Transport (Oxford University Press, Oxford, 1957). Bibliography 153 [96] V. Kourganov, Basic Methods in Transfer Problems (Oxford University Press, Ox ford, 1963). [97 [98 [99 [100 [101 [102 [103 [104 [105 [106 [107 [108 [109 [110 [111 [112 [113 [114; [115 [116 [117 J. S. Chang and G. Cooper, J. Computat. Phys. 6, 1 (1970). E. W. Larsen, C. D. Levermore, G. C. Pomraning and J. G. Sanderson, J. Compu tat. Phys. 61, 359 (1985). E. M. Epperlein, J. Computat. Phys. 112, 291 (1994); Laser and Particle Beams 12, 257 (1994). B. T. Park and V. Petrosian, Astrophys. J. Suppl. Series 103, 255 (1996); Astro-phys. J. 446, 699 (1995). A. K. Mitra, J. Math. Phys. 19, 2018 (1978). K. Banerjee, S.P. Bhatnagar, V.Choudhry and S.S.Kanwal, Proc. R. Soc. Lond. A 360, 575 (1978). K. Banerjee, Proc. R. Soc. Lond. A 364 (1978) 265 N. Bessis and G. Bessis, J. Math. Phys. 21, 2760, (1980). A. Hautot, J. Computat. Phys. 39, 72 (1981). F. M. Fernandez, G. A. Arteca and E. A. Castro, Physica A 122, 37 (1983). R. N. Chaudhuri and B. Mukherjee, J. Phys. A: Math. Gen. 16, 4031 (1983). G. A. Arteca, F. M. Fernandez and E. A. Castro, J. Math. Phys. 25, 932, 1984. M. R. M. Witwit, Pramana, J. Phys. 43, 279 (1994). M. R. M. Witwit, J. Computat. Phys. 129, 220 (1996). M. R. M. Witwit, J. Math. Chem. 19, 75 (1996). R. S. Kaushal, J. Phys. A: Math. Gen. 12, 253 (1979). N. Bessis, G. Bessis and G. Hadinger, J. Phys. A: Math. Gen. 16, 497 (1983). 'R. Dutt, A. Khare and U. P. Sukhatme, Am. J. Phys. 56, 163 (1988). A. Comtet, A. D. Bandrauk and D. K. Campbell, Phys. Letters B 150, 159 (1985). R. Dutt, A. Mukherjee and Y. P. Varshni, Phys. Rev. A 52, 1750 (1995). M. R. M. Witwit, J. Math. Chem. 20, 273 (1996). Bibliography 154 [118] C. S. Lai and H. E. Lin, J. Phys. A: Math. Gen. 15 495 (1982). [119] C. R. Handy, J. Phys. A: Math. Gen. 18, 3593 (1985). [120] F. M. Fernandez and E. A. Castro, J.Phys. A: Math. Gen. 20, 5541 (1987) [121] M. R. M. Witwit, J. Phys. A: Math. Gen. 24, 5291 (1991). [122] F. M. Fernandez and R. H. Tipping, Can. J. Phys. 74, 697 (1996). [123] S. Galicia and J.P. Killingbeck, Phys. Lett. A 71, 17 (1979). [124] V. Fack and G. Vanden Berghe, J. Phys. A: Math. Gen. 18 3355 (1985). [125] M. Braun, S. A. Sofianos, D. G. Papageorgiou and I. E. Lagaris, J. Computat. Phys. 126 315 (1996). [126] G. P. Flessas, Phys. Let. A 83 (1981) 121; ibid. 100, 383 (1984). [127] D. Singh and Y. P. Varshni, Phys. Rev. A 28, 2606 (1983). [128] Y. P. Varshni, Phys. Rev. 36 3009 (1987). [129] H. Scherrer, H. Risken and T. Leiber, Phys. Rev. A 38, 3949 (1988). [130] M. Kaluza, Comp. Phys. Comm. 79, 425 (1994). [131] A. Sinha, R. Roychoudhury and Y. P. Varshni, Can. J. Phys. 74, 39 (1996). [132] J.C. Light, I.P. Hamilton and J.V. Lill, J. Chem. Phys. 82(3) 1400 (1985); R.M. Whitnell and J.C. Light, ibid. 89(6), 3674 (1989); S.E. Choi and J.C. Light, ibid. 92(4), 2129 (1990); ibid. 97(10), 7031-7054 (1992); Z. Back and J.C. Light, Annu. Rev. Phys. Chem. 40, 469 (1989) [133] H. Wei and T.J. Carrington, Chem. Phys. 101, 1343 (1994); K. Sakimoto and K. Onda, ibid. 100, 1171 (1994); S. M. Auerbach and W.H. Miller, ibid. 100, 1103 (1994); J.K.G. Watson, Can. J. Phys. 72, 238 (1994); J. Szalay, J. Chem. Phys. 99, 1978 (1993); J. Tennyson and J.R. Henderson, ibid. 91, 3815 (1989). [134] W.J. Gautschi, Computat. Appl. Math. 12, 61 (1985). [135] P. Nevai, editor Orthogonal Polynomials: Theory and Practice (Kluwer Academic, Norwell, MA, 1990); G.Freud, Orthogonal Polynomials (Pergamon, New York, 1971); T. S. Chihara, An Introduction of Orthogonal Polynomials, (Gordon and Breach, New York, 1978). Bibliography 155 [136] P. Davis and P. Rabinowitz, Method of Numerical Integration. Academic, New York (1975). [137] R. Peyret and T.D. Taylor, Computational methods for fluid flow (Springer-Verlag, New York, 1983). [138] R. Peyret, Handbook of computational fluid mechanics. (Academic Press, London, 1996). [139] D.J. Evans, Preconditioning methods : analysis and applications. (Gordon and Breach Science, New York, 1983). [140] I. Gohberg, M. Hanke and I. Koltracht, SIAM J. Numer. Anal. 31(2), 429 (1994). [141] B. Koren and B. van Leer, Adv. Comput. Math 4, 127 (1995). [142] H. Chen and B.D. Shizgal, J. Math. Chem. In press. (1998) [143] A. S. Clarke and B. Shizgal, J. Computat. Phys. 104, 140 (1993). [144] J.D. Jackson, Mathematics for Quantum Mechanics. 93-94, (Academic Press, New York 1975). [145] C. Indira, M. C. Valsakumar, K. P. N. Murthy and G. Ananthakrisna, J. Stat. Phys. 33, 181 (1983). [146] K. Voigtlaender and H. Risken, Chem. Phys. Lett. 105, 506 (1984); H. F. Ouyang, Z. Q. Huang, and E. J. Ding, Phys. Rev. E 50, 2491 (1994); A. Debosscher, Phys. Rev. A 42, 4485 (1990); A. Debosscher, Phys. Rev. A 44, 908 (1991); G. Hu and Q. Zheng, Phys. Lett. A 110, 68 (1985); G. Hu, Phys. Lett. A 110, 253 (1985); H. Dekker and N. G. Van Kampen, Phys. Lett A 73, 374 (1979). [147] B. Shizgal and D. R. A. McMahon, Phys. Rev. A 32, 3669 (1985). [148] X. G. Wu and R. Kapral, J. Chem. Phys. 91, 5528 (1989) [149] I. L'Heureux, Phys. Rev. E 51, 2787 (1995); V. I. Melnikov, Phys. Rep. 209, 1 (1991); P. Hanggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990); P. Talkner, Chem. Phys. 180, 199 (1994). [150] R. Bonifacio and L. A. Lugiato, in Disspative Systems in Quantum Optics, edited by R. Bonifacio (Springer Verlag, Berlin, 1982); C. M. Bowden, M. Ciftan and H. R. Robl, editors Optical Bistability (Plenum Press, New York, 1981). Bibliography 156 [151] J. C. Englund, W. C. Schieve, W. Zurek and R. F. Gragg, in Optical Bistability edited by C. M. Bowden, M. Ciftan and H. R. Robl (Plenum Press, New York, 1981); A. R. Bulsara, W. C. Cheive and R. F. Gragg, Phys. Lett. A 68, 294 (1978); P. Hanggi, A. R. Bulsara and R, Janda, Phys. Rev. A 22, 671 (1980); A. Schenzle and H. Brand, Phys. Rev. A 20, 1628 (1979). [152] V. A. Shneidman, Phys. Rev. A 44, 8441 (1991); D. T. Wu, J. Chem. Phys. 97, 1922, 2644 (1992); I. Edrei and M. Gitterman, J. Chem. Phys. 85 190, (1986); G. Shi, H. Seinfeld, and K. Okuyama, Phys. Rev. A 41, 2101 (1990); B. Shizgal and J. C. Barrett, J. Chem. Phys. 91, 6505 (1989); L. Demeio and B. Shizgal, J. Chem. Phys. 98, 5713 (1993). [153] W. T. Coffey and D. S. F. Crothers, Phys. Rev. A 54, 4768 (1996); W. T. Coffey, D. S. F. Crothers, Yu. P. Kalmykov, E. S. Massawe a,nd J. T. Waldron, Phys. Rev. E 49, 1869 (1994). [154] C. Nicolis and G. Nicolis, Tellus 33, 225 (1981); C. Nicolis, Solar Physics 70, 473 (1981). [155] D. R. Chialvo and A. V. Apkaran, J. Stat. Phys. 70, 375 (1993); T. Heskes, J. Phys. A: Math. Gen. 27, 5145 (1994). [156] K. R. Yawn, B. N. Miller and A. Maier, Phys. Rev. E 52, 3390 (1995); L. Spitzer, Jr., Dynamical Evolution of Globular Clusters, (Princeton University Press, Prince ton, N.J., 1987); J. Binney and S. Tremaine, Galactic Dynamics, (Princeton Uni versity Press, Princeton, N.J., 1987) [157] A. N. Drozdov and M. Morillo, Phys. Rev. E 54, 3304 (1996). [158] P. Hanggi and P. Jung, Adv. Chem. Phys. 89, 239 (1995); J. Kuczka, P. Hanggi and A. Gadomski, Phys. Rev. E 51, 2933 (1995); A. J. Bray, A. J. McKane and T. J. Newman, Phys. Rev. A 41, 657 (1990); M. M. Dygas, B. J. Matkowsky and Z. Schuss, J. Chem. Phys. 83, 597 (1985); J. Masoliver, B. J. West and K.Lindenberg, Phys. Rev. A 35, 3086 (1987); R. F. Fox, Phys. Rev. A 37, 911 (1988); S. A. Adelman, J. Chem. Phys. 64, 124 (1976); [159] H. A. Kramers, Physica 7, 284 (1940). [160] M. J. Englefield, Physica A 167, 877 (1990); H. Okamoto, J. Phys. A: Math. Gen. 23, 5535 (1990), T: Miyazawa, Phys. Rev. A 39, 1447 (1989); M. J. Englefield, J. Stat. Phys. 52, 369 (1988); H. Risken and K. Voigtlander, Z. Phys. B 54 253, (1984). [161] D.W.Noid and R.A.Marcus, J. Chem. Phys. 67(2), 559 (1977). Bibliography 157 [162] D.W.Noid, M. L. Koszykowski, M. Tabor, and R.A.Marcus, J. Chem. Phys. 72(2), 6169 (1980). [163] M.J. Davis and E. J. Heller, J. chem. Phys. 71(8), 3383 (1979). [164] M.D.Feit, J.A.Fleck, Jr., and A.Steiger, J. comput. Phys. 47, 412 (1982). [165] R. Kosloff and H. Val-Ezer. Chem. Phys. Letter 223, 127 (1986). [166] G. C. Groenboom and H. M. Buck, J. Chem. Phys. 92(7), 4374 (1990). [167] R.Q. Chen and H. Guo, J. chem. Phys. 105(4), 1311 (1996). [168] O.A. Sharafeddin, Chem. Phys. Lett. 247, 470 (1995). [169] B.I. Henry and J. Grindlay, Can. J. Phys. 75(8), 517 (1997). [170] P.K. Chattaraj96 and S. Sengupta, Indian J. Pure Appl. Phys. 34(8) 518 (1996). [171] S. Sengupta and P.K. Chattaraj96, Phys. Lett. A 215, 119 (1996). [172] W.M. Vieira and P.S. Letelier, Phys. Rev. Lett. 76, 1409 (1996). [173] P. Hatchison and R.E. Wyatt, Phys. Rev. A 23, 1567 (1981). [174] Noid, D.W. and Marcus, R.A. J Chem. Phys. 62(6), 2119 (1975).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The quadrature discretization method and its applications
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The quadrature discretization method and its applications Chen, Heli 1998
pdf
Page Metadata
Item Metadata
Title | The quadrature discretization method and its applications |
Creator |
Chen, Heli |
Date | 1998 |
Date Issued | 2009-06-19T23:52:56Z |
Description | A discretization method referred to as the quadrature discretization method is introduced and studied in this thesis. The quadrature discretization method is a spectral method based on a grid of points that coincide with the points of a quadrature. The quadrature is based on a set of nonclassical polynomials orthogonal with respect to some weight function over some interval. The method is flexible with respect to the choice of the weight function and quadrature points so that the optimum accuracy and convergence in the solution of differential equations and/or partial differential equations can be obtained. The properties of the quadrature discretization method are studied and compared with classical spectral methods as well as finite difference methods. Several analytic model problems in one and three dimension are studied. The quadrature discretization method competes well with classical spectral methods and is far more superior than the finite difference method. In some cases it provides significant improvement in the accuracy and convergence of the solution of the problem over other methods. The main objective of the quadrature discretization method is to determine the weight function that defines the polynomial basis set and hence the grid points that provide optimum convergence in a given application. The quadrature discretization method is applied to a large class of time dependent Fokker-Planck equations. Several weight functions are used and the results are compared with several other methods. The weight functions that have often provided rapid convergence of the eigenvalues and eigenfunctions of the Fokker-Planck operator are the steady solutions at infinite time. The quadrature discretization method is also employed in the solution of Schrodinger equations. The weight functions that are used are related to the ground state wave functions if known, or some approximate form. The eigenvalues and eigenfunctions of four different potential functions discussed extensively in the literature are calculated and compared with the published values. The eigen-problem of a two-dimensional Schrodinger equation with the Henon-Heiles potential is also calculated with the quadrature discretization method. The rate of convergence of the eigenvalues and eigenfunctions of the Schrodinger equations is very rapid with this approach. |
Extent | 6700172 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-06-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080038 |
URI | http://hdl.handle.net/2429/9493 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- ubc_1998-345408.pdf [ 6.39MB ]
- [if-you-see-this-DO-NOT-CLICK]
- [if-you-see-this-DO-NOT-CLICK]
- Metadata
- JSON: 1.0080038.json
- JSON-LD: 1.0080038+ld.json
- RDF/XML (Pretty): 1.0080038.xml
- RDF/JSON: 1.0080038+rdf.json
- Turtle: 1.0080038+rdf-turtle.txt
- N-Triples: 1.0080038+rdf-ntriples.txt
- Original Record: 1.0080038 +original-record.json
- Full Text
- 1.0080038.txt
- Citation
- 1.0080038.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 48 | 13 |
United States | 43 | 12 |
Morocco | 7 | 1 |
United Kingdom | 7 | 0 |
Canada | 6 | 3 |
Vietnam | 5 | 0 |
India | 4 | 0 |
Iraq | 4 | 0 |
Iran | 2 | 0 |
Japan | 2 | 0 |
France | 1 | 0 |
Russia | 1 | 0 |
Algeria | 1 | 0 |
City | Views | Downloads |
---|---|---|
Washington | 30 | 0 |
Unknown | 29 | 12 |
Hangzhou | 20 | 0 |
Henderson | 6 | 0 |
Guangzhou | 5 | 0 |
Hanoi | 5 | 0 |
Changsha | 5 | 0 |
Shenzhen | 4 | 12 |
Ashburn | 3 | 0 |
Beijing | 3 | 1 |
Guiyang | 2 | 0 |
Ürümqi | 2 | 0 |
Kaul | 2 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080038/manifest