Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The quadrature discretization method and its applications Chen, Heli 1998

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1998-345408.pdf [ 6.39MB ]
Metadata
JSON: 831-1.0080038.json
JSON-LD: 831-1.0080038-ld.json
RDF/XML (Pretty): 831-1.0080038-rdf.xml
RDF/JSON: 831-1.0080038-rdf.json
Turtle: 831-1.0080038-turtle.txt
N-Triples: 831-1.0080038-rdf-ntriples.txt
Original Record: 831-1.0080038-source.json
Full Text
831-1.0080038-fulltext.txt
Citation
831-1.0080038.ris

Full Text

T H E Q U A D R A T U R E D I S C R E T I Z A T I O N M E T H O D A N D ITS APPLICATIONS By Heli Chen B. Sc. (Mechanics), Peking University, 1987 M . Sc. (Mechanical Engineering), University of Calgary, 1993  A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF PHILOSOPHY  in T H E FACULTY OF GRADUATE STUDIES  (Department of Mathematics) We accept this thesis as conforming to the required standard  T H E 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 department  or  by  his  or  scholarly purposes may be granted by the head of her  representatives.  It  is  understood  that  copying  my 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  DE-6 (2/88)  Oct,  /I  ,  rfU  Abstract  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 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 discretization 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  2  The Quadrature Discretization Method (QDM)  14  2.1  Classical p o l y n o m i a l s and the pseudospectral m e t h o d  15  2.2  T h e quadrature discretization method  28  2.2.1  C o n s t r u c t i o n of d e r i v a t i v e m a t r i c e s  28  2.2.2  G e n e r a t i o n of o r t h o g o n a l p o l y n o m i a l s , q u a d r a t u r e w e i g h t s a n d p o i n t s  33  3  1  2.3  T i m e integration  42  2.4  N u m e r i c a l tests a n d results  45  2.4.1  One-dimensional Poisson equation  45  2.4.2  A n analytic example  46  2.4.3  Three-dimensional Poisson problems  .  60  Applications of the Q D M to the Solution of Fokker-Planck Equation  66  3.1  66  Introduction  iv  4  3.2  The Q D M matrix representation of the Fokker-Planck operator  69  3.3  Eigenfunction expansion  72  3.4  Temporal discretization  73  3.5  Applications and results  74  3.6  Summary  A p p l i c a t i o n s of the Q D M to the S o l u t i o n of Schrodinger equation  103  4.1  103  4.2  5  102  Introduction The equivalence of the Fokker-Planck eigenvalue problem and the Schrodinger equation  106  4.3  The Q D M matrix representation of the Schrodinger equations  108  4.4  Applications and results  110  4.5  Summary  140  S u m m a r y and future work  142  Bibliography  148  v  L i s t of Tables  2.1  Errors for the numerical approximation of u" =  2.2  Maximum error in a  2.3  Maximum error in the numerical solution of a  n  18  sin(7rx)  and 8 calculated by the Stieltjes's procedure . . .  36  n  n  and (3 for the Hermite n  polynomials  38  2.4  Maximum error for the numerical approximation of points and weights  .  2.5  Standard error for the numerical approximation of u" =  2.6  Numerical error in the eigenvalues for the Hermite and Q D M methods . .  48  2.7  Numerical error in the eigenvalues for the Legendre method  49  2.8  Domain generated from the Q D M and Hermite weight function  52  2.9  Maximum time step in the F E time integration  55  sin(7rx)  40 45  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 C P U time for the numerical approximations for example A  62  2.13 Comparison of the standard Errors and C P U times for the three numerical approximations  63  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  3.3  potential  80  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) = y + J^I,  4.2  Convergence of the eigenvalues with V(y) = y +  4.3  Convergence of the eigenvalues with V(y) = y + \ Q 2, u>i(y) = exp(—ay ). 116  4.4  Convergence of the eigenvalues with V(y) = y + T+^2> 9  2  2  2  2  /  2  V )-  1 / 2  114  2  = exp( —ay ). 115 2  2  1) i( )  =  w  =  v  , • 117  )  P  , g = 10, w\(y) =  )  118  Convergence of the eigenvalues with V(y) = y + ^ 2, 9 = 100, vh{y) = 2  1  exp(-y (l + T ^ - ) ' / ) 2  4.7  — a  0  2  2  p(  V  Convergence of the eigenvalues with V(y) = y + ex (-y (l + T T ^ )  4.6  , \{y)  2  P  e x  w  2  1+  ex (-y (l + T ^ Y 4.5  \{y)  w  GY  ' . .'  2  119  Convergence of the eigenvalues with V(y) = y + T+^pi 9 2  =  1) i{u) — w  ex (-a )/(l+gy ) * 2  V  4.8  2  120  a  iy  Convergence of the eigenvalues with V(y) = y + j ^ ~ 2 , 9 = 10, w\(y) = 2  ex (-a )/(l+gy )^ 2  V  4.9  121  2  W  The convergence of eigenvalues with V(y) = y +  ; 9 = 100, wi(y) =  2  exp(-a  2 i y  )/(l + ^ ) ^  122  2  4.10 ( a i , a ) used for Tables 4.7-4.9  123  2  4.11 Comparison of the results of Ai with V(y) = y +  124  2  4.12 Convergence of the eigenvalues for SE with V(y) = y — 3y 6  130  2  4.13 Convergence of the eigenvalues of even parity with V(y) = \y  2  4.14 Convergence of the eigenvalues with V(y) = y + ey 2  vii  4  + 2y + |y .132 4  6  134  4.15 Convergence of the eigenvalues with V(y) = y + ey calculated by fitting 2  4  weight function to ground state eigenfunction  135  4.16 Convergence of the eigenvalues, A ^ , for the Henon-Heiles potential with n  Hermite points  137  4.17 Convergence of the eigenvalues, A ;, for the Henon-Heiles potential with n>  new points  139  4.18 Eigenvalues for the Henon-Heiles potential  vm  141  L i s t of Figures  2.1  Maximum error E°° in variation with N for the solution of u" =  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-  sin(7rx).  endre discretization 2.5  23  Largest eigenvalues of the 2nd derivative matrix with Chebyshev and Legendre discretization  2.6  25  Largest eigenvalues of the 1st and 2nd derivative matrices calculated by the KoslofF modified method  2.7  19  27  Largest eigenvalues of the 1st and 2nd derivative matrices calculated by the Hermite method  29  2.8  Errors in a  37  2.9  Minimum spacing of points versus number of quadrature points N.  n  and B calculated by the Stieltjes's procedure n  . . .  41  2.10 Numerical approximation of the eigenvalues for the analytic F P E problem.  51  2.11 Largest numerical eigenvalues for the analytic F P E example  54  2.12 Numerical solution of the time dependent F P E with N = 6  56  2.13 Numerical solution of the time dependent F P E 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  76  The first ten eigenfunctions for the quartic potential model  ix  3.2  The variation of the lowest nonzero eigenvalue, Ai versus e for the quartic potential  3.3  83  The potential U{x) in the equilibrium distribution for the optical bistability problem.  3.4  86  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 N P O potential  126  4.3  Ground state eigenfunction for the N P O potential with g = \ — 100, N = 25  127  4.4  Ground state eigenfunction for the N P O potential  128  4.5  Contour plots of the eigenfunctions of the Schrodinger equation for the Henon-FIeiles potential  138  x  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. T h e major numerical techniques widely used include finite difference method ( F D ) , finite element method ( F E M ) , and spectral methods. T h e m a i n objective of these numerical methods is to approximate derivatives by algebraic expressions involving the solution evaluated on a grid.  Ordinary differential  equations ( O D E ) are thus reduced to algebraic equations and partial differential equations ( P D E ) ( P D E ) are reduced to systems of, O D E s .  In this section, we shall confine our  discussion to linear differential equations. A differential equation can be written as  Lu where u(x)  =  -f{x)  is the solution of the equation, and L is the linear differential operator.  Discretizing the variable £ on a grid w i t h N points, E q . (1.0.1) can be approximated by algebraic equations  w i t h L N , the m a t r i x representative of the differential operator L, and UJV, the approximate solution at the grid points. W i t h different numerical techniques, the choice of the 1  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, -^-  = L  + f.  N  N  (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 - ^  = Lu. N  (1.0.5)  N  There are several ways to solve the O D E system Eq. (1.0.5). The solution of this time dependent P D E can be evaluated by an eigenfunction expansion. The formal solution (with f(x) = 0) is then of the form u(x,t) = ^ c e " V n ( x ) , A  n  (1.0.6)  n=0  where A and <f> (x) are the eigenvalues and eigenfunctions of the eigenvalue problem n  n  L(f) (x) - \ <f> (x), n  n  n  (1.0.7)  the coefficients c are determined from the initial condition, u(x,0), and are given by n  c = J u(x, 0)(j) (x)dx. n  n  (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 O D E system Eq. (1.0.5) in time, and thus integrate the O D E 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 RungeKutta 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 CkM )x  k=i  The transform coefficients  (1-0-9)  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 functionals. 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 approximation (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 F D 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 F E M is frequently used to solve partial differential equations that occur in engineering applications, especially, in the areas of solid mechanics, elasticity, and computational 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 matrix 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 F E M is therefore very useful in solving problems with complex geometry in multi-dimensional problems. The F E M 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 disvantage of the F E M 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) « u {x) = N  ak4>k{x)  (1.0.11)  k=0  in a set of orthogonal basis functions (trial functions), (j>k{x). When this series is substituted into Eq. (1.0.1), we can define a residual function R(x] a , a i , ...,a^) = LUN + / . 0  (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 F E M . The choice of the test functions distinguishes between the spectral Galerkin, tau and collocation (also referred to as pseudospectral) 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 methods 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 nonlinear 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 F D methods when the order of accuracy tends to infinity [24]. Unlike low order F D method and F E M 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 F E M is particularly well suited to problems in very complex geometries, whereas spectral methods can offer superior accuracies 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 convergence 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 decomposition or spectral element methods [25-32]. The techniques can apply to problems with complicated domain and provide superior accuracy and convergence over the F E M . The advantage of spectral methods is that they are more accurate than the F D method  Chapter 1.  Introduction  7  and the F E M 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 stability 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 advantages 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 accuracy 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 solution 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(N ) 2  condition num-  ber. For spectral methods, the largest eigenvalue of the first derivative operator grows as 0(N ), 2  and of the second derivative operator grows as 0(N ), 4  as N increases, leading to  Chapter 1.  Introduction  8  poorly conditioned matrices [18,19,29,36,39,52-58]. This effect upon the numerical computations 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 using 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 matrices 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 solution 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 expensive for large matrices. Iterative schemes are sometimes a practical necessity. Therefore, reduction of the condition number of spectral differential matrices is necessary for stability 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 F D [61,62], the F E M [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) (l a  + 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~ on [0, co] and for w(x) = e~* on [-co, oo], x  2  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 characteristics 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 differentiation matrices are 0(\/ N) and O(N), respectively [29,68,73]. Numerical results from r  our work also showed that the condition number for a second differentiation matrix can be improved to be between O(N) and 0(N ) 2  depending on the physical problem and  polynomial basis set chosen. This is better than second order finite difference which is 0(N ), 2  and far better than 0(N ), 4  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). In this thesis we introduce the quadrature discretization method (QDM), which allows nonclassical polynomials as basis functions in a spectral approximation. The idea of the Q D M 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 x exp(—x ) was used in this work. The method was later developed fur2  2  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 nonclassical 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 Q D M 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 Q D M . 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 F P E 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  Chapter 1.  Introduction  12  F P E . 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 F P E 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 F P E is a spectral method which usually 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. A n alternate approach involves the discretization of the P D F 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]. Fourier series or Chebyshev polynomials are almost exclusively chosen as basis functions in the application 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 F P E 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 include variational and spectral techniques with Gaussian-type basis functions [29,101-111], perturbative approach [112-117], Pade approximation [118-122], F D methods [123-125],  Chapter 1.  Introduction  13  and some other methods [126-131]. A very similar approach to the Q D M 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 Q D M and D V R differ in the use of polynomial basis functions. The Q D M uses nonclassical polynomials and the D R V 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 connection 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 Q D M 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 F D methods, it is possible for the Q D M discretization to preserve the symmetry for the self-joint joint differential operator. Furthermore, the flexibility of the Q D M 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 Q D M are introduced. The Q D M is tested and compared with other numerical methods. The applications of the Q D M to the F P E s 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 Q D M , 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 Q D M 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 Q D M 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 Q D M 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 Q D M 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 Q D M , 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), E  2  and E°° (maximum error) defined by  14  (standard error)  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 p o l y n o m i a l s and the pseudospectral m e t h o d  Classical polynomials referred to in this thesis are the polynomials used by most researchers 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 / . A are the eigenvalues and u {x) are the corresponding eigenfunctions. n  n  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) (l a+1  + xY+\  b{x) = 0,  w{x) = ( l - x ) ( l + x)^ a  (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) =  x e~ a+1  b(x) = 0,  x  i  w(x) = x e~\  (2.1.3)  a  for \/x G (0, oo), a > — 1, whereas for Hermite polynomials, these are given by, a(x) = w(x) = e~ \  b(x) = 0,  x  for  (2.1.4)  6 ( — 00,00). The polynomials satisfy a three recurrence relation [39] U (x)  =  n+1  where u (x) 0  (p X n  + O )u (x) n  n  +  (2.1.5)  T U --L, n  n  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 = ———xL (x)  L (x) n+1  n —L _i(x),  n  77. — — I  Tl -\- 1  _L  n  (2.1.6)  with Lo(x) = 1, Li(x) = x; Chebyshev T {x) n+1  = 2xT {x)-T ^(x),  (2.1.7)  = 2xHn(x)-2nHn-1,  (2.1.8)  n  n  with To(x) = 1, T\(x) = x; Hermite H (x) n+1  with H (x) = 1, H^x) 0  = 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 J  I  m J2f(xj)w , 3  3=1  (2.1.9)  Chapter 2. The Quadrature Discretization Method (QDM)  where Wj = J lj(x)w(x)dx,  17  are the weights, and Xj are the quadrature points of the  I  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, lj(xi) = { ll,  iii^j (2.1.10) 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 function 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  (D f)( ) N  = £(Av), /(*-. ),  Xi  ?  (i = 0 , 1 , 2 , A O ,  ;  .  (2.1.11)  3=0  where DN is the derivative matrix with the entries (Dw)ij = l'j(xi). Higher order derivative 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 du 2  —  . ,  N  = sm(7rx)  .  .  (2.1.12)  with homogeneous Dirichlet boundary condition u( — 1) = u(l) = 0. The problem has an analytic solution u(x) = — ir sin(7rai). 2  We solve the problem with a Legendre spectral  Chapter 2. The Quadrature Discretization Method (QDM)  E  N  4 8 16 32 64 128 256 512  x  4.09(-02) 2.40(-04) 1.50(-11) 2.73(-14) 1.21(-14)  Legendre E 2  5.01(-02) 2.83(-04) 1.72(-11) 3.23(-14) 1.52(-14)  E°°  FD E  Ex  6.13(-02) 4.46(-04) 2.91(-11) 6.03(-14) 4.36(-14)  18  E°°  2  1.56(-01) 3.66(-02) 8.68(-03) 2.11(-03) 5.19(-04) 1.29(-04)  1.91(-01) 4.01(-02) 9.46(-03) 2.31(-03) 5.73(-04) 1.43(-04)  Table 2.1: Errors for the numerical approximation of u" = spectral method and a finite difference method.  2.34(-01) 5.30(-02) 1.30(-02) 3.22(-03) 8.04(-04) 2.01(-04) 5.01(-05) 1.26(-05)  sin(7rx)  with a Legendre  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~ ) within 32 points. Increasing TV beyond this 14  point would not reduce the errors significantly, and the errors remain close to machine accuracy. For the F D method, the error decays at a much slower rate. With 128 grid points, the error is O ( 1 0 ) and with 512 grid points, it is O ( 1 0 ) . This indicates that to obtain -4  -5  4 decimal places accuracy for the solution of Eq. (2.1.12), the F D 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 F D 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 polynomials 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 = »Mx)  ^j^  on (-1,1)  X  (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  eigenvalue grows like 0(N ) 2  is plotted in Fig. 2.4. The largest  as N —> oo. For the Le gendre method, |A/v| ~ 0.0S0N , and 2  for the Chebyshev method, \ X \ w 0.089N . 2  N  For the second derivative matrix, we consider the eigenvalue problem of the diffusion operator = X f ()  ^M  n( )n  x  o  n  (-1,1)  (2.1.15)  with the boundary conditions of Dirichlet type: </> (l) = </>„( —1) = 0. The problem has n  the analytical solution, X = — (nir/2)  2  N  and (j) (x) = siny/X^x, n = 0,1,2,.... n  Chapter 2. The Quadrature Discretization Method (QDM)  21  s> oh  128  I  1  1  1  1  1  32  256  512  512  1  N=32-  o  256  64 h  o  -  ro  16  -  Oh -  -  o  -256  -64  h  0 -128 -128  1  1 , 1 , -64 0 Real  -512  1  64  128  -512  -256  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 0.40 1  '  1  0.80  '  1 ' 1 1.20 1.60 log (N)  '  1 2.00  '  1  2.40  1 0  Figure 2.4: Eigenvalues of 1st derivative matrix with Chebyshev and Legendre discretization. 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(N ) 4  as N —» co, and |A | = n  0.026A and | A | = 0.047iV for the Legendre and Chebyshev method, respectively. r4  4  n  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 A t to ensure the stability of the time integration may be 0(N~ )  for a 2nd order time-dependent differential equation. To  4  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 A s  m m  between the adjacent grid  points. The largest eigenvalue (in modulus) of the first derivative matrix is of the order of  1 / A x  m  i , n  and for the second derivative matrix it is of the order of  1 / A x 2  m  i . n  Most  classical points in finite domains, for example, Chebyshev and Legendre points, have minimum spacing of 0(1/N ) 2  0(N ) 2  and 0(N ), 4  such that the 1st and 2nd derivative matrices grow like  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 discretization. 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) the largest eigenvalues of 1st and 2nd derivative matrices improve to 0(N)  such that and  0(N ), 2  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.2SN ^ for the first derivative 0  matrix and XN ~ 0.3SN '  1 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, F (x), n  orthogonal with respect to a weight function w(x),  that is, = 8  J w(x)F (x)F (x)dx n  m  n m  .  (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 Q D M 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 Q D M is based on the discrete representation of a function j\x) the set of N quadrature points Xi, that is, f(xi).  by its values at  We expand f(x) in the {F } orthogonal n  basis set, N  f{x) = J2a F (x), n  n  (2.2.2)  n=l  where a are the expansion coefficients of f(x) and are given by, n  a = jw(x)f(x)F (x)dx. n  n  (2.2.3)  The expansion coefficients can be evaluated numerically with the quadrature rule in Eq.(2.1.9), and N  a = Y, W3f(x,)F (xj). n  n  (2.2.4)  Chapter 2. The Quadrature Discretization Method (QDM)  0.00 0.40 1  1  0.80  1  1  1  1.20  log  1 0  29  '  1  1.60  1  2.00  (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 {a }. n  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 {a } n  to the physical space {/(a^)}, N  f{x ) = '£a F {x ). i  n  n  (2.2.5)  i  71=1  Substituting Eq. (2.2.4) into Eq. (2.2.2), we obtain f(x) = J2J:w f(x )F (x )F (x). l  l  n  1  n  (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 ' i W / W .  (--) 2  2  7  where the interpolation polynomial lj(x) is given by N  71=1  and satisfies (0, ifi^j Ijixi) = <| 1, i f i = j If Eq. (2.2.7) is differentiated, one obtains DF{XI)  -EA./fe),  (2.2.9)  (2.2.10)  where A,-j are the entries of the derivative matrix A and are defined by  = /;(.'.•,) = E 71 =  '^./^(x,)/.;:!,-;).  (2.2.11;  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 F (x) are a , then there is a unitary n  n  (or orthogonal) transformation between the physical and transform spaces, that is, N  f( *) = E V^Fn{  (2.2.13)  X  n=l  N  ^ ^ E ^ - ' N / W .  (2.2.14)  i=i  where the matrix elements of the symmetric transformation T are T; =  s/w^n^Xi).  n  If in Eq. (2.2.13), Xi is replaced with x, and Eq. (2.2.14) is used for a , one obtains n  a A^th-order interpolation, given by  / » = E4(*)/(*;),  (2-2-15)  3=1  where the interpolation polynomial Ij(x) is given by AT  I {x) = /w^J]^2F {x )F (x). j  y  n  j  (2.2.16)  n  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 N df{ Xi - E dx . .  (2-2.17)  .7 = 1  where the derivative operator, Dij = y/W{W T'j(xi)  is given by  3  N  D  = ^JT  i3  3  E Fn(x )F^{x ). 3  n=l  (2.2.18)  t  We refer to the derivative matrix D with the elements Dij as the Q D M modified derivative matrix. We have the following relation between the derivative matrices A,- and Dij 7  &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 d  nm  = j w(x)F {x)F' {x)dx. n  m  (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  d  E WkF (x )F' (x ).  =  nm  n  k  n  (2.2.21)  k  For some specific physical models, using the transform matrix have an attractive advantage 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 Q D M derivative matrix D = T d T , that is. T  N  v  D  N  N  =EE  B ( )^/u7 [J2 n Xl  z  n=l m=l  (2.2.22)  w F (x )F^(x )}F (x )^IJ]. k  n  k  k  m  3  k=l  The sum over n yields 6i and the sum over k can subsequently be performed and one k  obtains Dij given by Eq. (2.2.18) The second derivative matrix A ^ ' or D ^ ' can be obtained simply by multiplying it 2  2  by itself. It is easy to show that the matrix elements for A ^ ' and D^ ) satisfies 2  Ag> = y/wj/wM?-  2  (2-2.23)  The application to differential equations is based on the algorithm for numerical differentiation, defined by, df lf-U*,  N  =  E vf( ), A  i=i  x  3  (2-2.24)  or df  N  V^[^],=,, = E A , / ( ^ ) v % , where the matrix A „ and  (2.2,25)  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 Q D M 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 Q D M 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  G e n e r a t i o n of orthogonal p o l y n o m i a l s , 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 procedure, has been provided recently by Gautschi [134].  Basis set generation; T h e 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 {Q { )} x  n  with Qo = 1,  given by Q +i{x) n  = (x - a )Q (x) n  n  - 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  (2.2.27)  — ,  (2.2.28)  Pv1  Pn =  PnPn-l  anc  Q  n  =  — • Pn-\pn-1---P\  For the classical polynomials, the coefficients a example, a  and P are known and exact. For n  = 0, Pi = 1/3 and P = n / ( 4 n — 1) (n > 1) for the Legendre case; a 2  n  n  (2.2.29)  2  n  Pi = 1/2 and P = 1/4 (n > 1) for the Chebyshev case; and a n  n  n  = 0,  = 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 P ( x ) = <5n(^)/ 7n satisfy /  n  v  the recurrence relation xP (x) n  = ^/P iP (x) n+  + a P (x)  n+1  n  n  + ^JP P _,(x), n  n  (2.2.30)  where the normalization factors are given by w(x)Q  7n = /  (x)dx.  2 n  (2.2.31)  For the Q D M , we are interested in the non-classical sets of orthogonal polynomials and a , P are usually unknown and calculated numerically. The recurrence coefficients n  in Eq.  n  (2.2.26)  are given by _ a  n  J  w(x)xQ (x)dx n  f w(x)Ql(x)dx  ~  >  and f /  w(x)Ql(x)dx w(x)Q  2 n  _ ( )dy 1  X  K  L  L  i  l  )  Chapter 2. The Quadrature Discretization Method (QDM)  35  One of the practical approaches of numerical computation of a  and B in Eqs. (2.2.32)  n  n  and (2.2.33) is the Stieltjes's procedure introduced by Gautschi [134]. The method involves the accurate calculation of the integrals in Eqs. (2.2.32) and (2.2.33) by subdividing 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, Q = 1 and 0  o = A) = 0. Equation (2.2.26) can then be used to generate Qi and then ot\ and B\ with  a  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 discretized 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 a  n  and 3 for most non-classical recurrence n  relations, the accuracy of the numerical calculation of a  n  and B will be essential to the n  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 a  n  and 8 for up n  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 a and B computed by Stieltjes's procedure for several n  n  classical polynomials with TV quadrature points are shown in Table 2.2. The error curves for the a  n  and B with N = 50 is given in Fig. 2.8. n  As shown in the table, the  maximum errors in B for both Legendre and Hermite polynomials and in a for all the n  n  polynomials are very small and at most O(10~ ). The maximum error in B for the 13  n  Chebyshev polynomials are 7.76 x 10~ . It is unchanged with changing 5  because the  maximum error eventually occurs at the first coefficient B^. The reduced accuracy of  Chapter 2. The Quadrature Discretization Method (QDM)  36  N  5 10 20 50 100  Chebyshev 9 528058591081317e5 386806816078257e9 710341644410283e1 805638645361598e2 490836730411662e-  17 16 16 15 15  7 7 7 7 7  763988590681060e- 05 763988590681060e- 05 763988590681060e- 05 763988590681060e- 05 763988590681060e- 05  5 10 50 100  Legendre 7 070397677170780e7 070397677170780e4 679789671731484e6 021950962055382e-  16 16 15 15  3 3 3 3  053113317719181e053113317719181e053113317719181e164135620181696e-  15 15 15 15  5 10 50 70 100  4 4 2 5 5  15 15 14 14 14.  9 103828801930001e4 840572387370000e3 375077994860000e4 760636329590000e7 034373084020000e-  15 14 13 13 13  Hermite 996208974320000e996208974320000e609915956900000e530643324450000e530643324450000e-  Table 2.2: Maximum error in ct and f3 calculated by the Stieltjes's procedure. n  n  the recurrence coefficients /3 's for the Chebyshev polynomial is due to the less accurate n  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)  .20.00 0.00 1  ' 10.00  1  '  ' 20.00  1  1  30.00  37  1 ' 40.00  1  50.00  n  Figure 2.8: Errors in a and B 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. n  n  Chapter 2.  The Quadrature Discretization Method (QDM)  N  E°° 01 [-20, 20] 996208974320000e-15 996208974320000e-15 609915956900000e- 14 530643324450000e- 14 530643324450000e- 14 530643324450000e- 14  38  N  5 10 50 70 90 100  4 4 2 5 5 5  5 10 50 100  [-5, 5] 2 346651055977530e1 263572315381132e3 598600987376839e8 320267776295266e-  5 10 50 70 >91  [-50, 50] 658312592764631e6 1 382306109940420e1 895834758617558e3 565939958333915eoverflow  9.103828801930001e-15 4.840572387370000e-14 3.375077994860000e-13 4.760636329590000e-13 7.034373084020000e-13 7.034373084020000e-13  15 14 14 14  7.246219748680005e-06 2.919377614958041e-02 1.874888749856461e+01 4.374979589455547e+01  16 15 14 14  1.154631945610163e-14 3.907985046680551e-14 1.492139745096210e-13 2.060573933704291e-13  Table 2.3: Maximum error in the numerical solution of a nomials on three intervals [-20,20], [-5, 5] and [-50, 50].  n  and B for the Hermite polyn  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 a  n  and fi . If the domain is too large, the integrals can be very large n  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 a and /3 . The computer overflow in the integration of low n  n  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 a  n  and ft for the Hermite polynon  mials calculated by Stieltjes's procedure with several choices of integration intervals are given in Table 2.3.  C a l c u l a t i o n of quadrature points and weights If the a  n  and f3 for the recurrence relation are known, the quadrature points which n  are the roots of the Nth polynomial constructed from the recurrence can be calculated by diagonalizing the tri-diagonal matrix [136] /  J =  V  0  0  0  0  0  0  ayy-i  /SJV-I  \  0  y/p\  a  0  0  0  ...  0  0  0  ... PN-I CJV J  2  (2.2.34)  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 a  n  and /3 , we n  Chapter 2. The Quadrature Discretization Method (QDM)  poo  poo  5 20 50 100  Chebyshev 2.954894016793386e-05 7.757441346845083e-06 3.121430352637233e-06 1.569879214957659e-06  7.806354429684692e-05 2.321553559997414e-05 9.635848417718917e-06 4.898200317991069e-06  5 50 100  Legendre 2.220446049250313e- 15 1.415534356397075& 15 3.885780586188049e- 15  3.774758283725532e-15 5.384581669432010e-15 7.507883204027621e-15  5 50 70 100  Hermite [-20, 20] 4.440892098500626e- 15 5.773159728050815e- 14 7.460698725481052e- 14 9.592326932761353e- 14  1.032507412901395e-14 1.271205363195804e-14 1.068589661201713e-14 1.276756478318930e-14  5 50 100  [-5, 5] 1.495381356519943e-07 4.192398352706889e+00 8.408353932924946e+00  5.903613908841976e-08 5.962550119618734e-02 8.056316459982088e-02  5 50 70  [-50, 50] 2.664535259100376e-15 5.417888360170763e-14 4.618527782440652e-14  2.109423746787797e-15 5.828670879282073e-15 5.190292640122607e-15  N  Table 2.4: Maximum error for the numerical approximation of points and weights.  40  Chapter 2. The Quadrature Discretization Method (QDM)  0.80  1.20 log (N)  1.60  41  2.00  1 0  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)  compare the results calculated directly from the exact a  n  42  and 8 , as well as the exact n  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 a  n  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 A i x m n  ~ 9.87/N , A i x 2  m n  « 12.2/A^ , and A i x 2  m n  &  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/N ) 2  for the Chebyshev  and Legendre points and 0(1 /s/N) for the Hermite points as number of quadrature points N increases.  2.3  T i m e integration  In most applications to the P D E , the spatial discretization uses the Q D M but the temporal discretization uses conventional finite difference. Consider a time dependent P D E | ^ = /(u,i)  <>0,  (2.3.1)  where the general (linear or nonlinear) operator / contains the spatial part of the P D E . After discretizing the spatial operator by the Q D M , 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  =  L  U  '  <--» 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> (x) = -e <f) (x), n  n  n  (2.3.4)  where e are the eigenvalues and 4> (x) are the eigenfunctions. The linear time-dependent n  n  equation, admits a solution of the form  oo u(x,t) =  5> e- "Vn(x), e  n  (2.3.5)  71=0  where the coefficients c are determined from the initial condition, u(x, 0), and are given n  by c = J u(x,0)<f> (x)dx. n  n  (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 approximate 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 order 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  +  A t L  n  ^ g  For the backward Euler method, U  n+1  = U + AtL \ n  (2.3.9)  n+  and for the Crank-Nicolson method, Tjn+l  rjn  =z  }-&t[L  n+1  +  + L }. n  (2.3.10)  The F E method is an explicit method and the absolute stability region is |1 + A A i | < 1. If the eigenvalue A is real, then —2 < XAt < 0. Both F E and C N 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)  N 4 6 8 10 12 14 16 20 25  QDM1 7.08(-03) 1.59(-02) 5.25(-04) 2.66(-05) 4.43(-07) 8.77(-09) 8.59(-ll) 1.20(-14) 4.95(-14)  QDM2 3.11(-03) 4.47(-03) 1.99(-04) 5.43(-06) 9.86(-08) 1.29(-09) 1.28(-11) 4.44(-15) 7.35(-15)  QDM3 9.26(-02) 9.26(-03) 4.01(-04) 1.04(-05) 1.84(-07) 2.36(-09) 2.31(-11) 3.48(-15) 1.05(-14)  45  QDM4 9.62(-02) 1.18(-02) 6.10(-04) 1.85(-05) 3.75(-07) 5.43(-09) 5.90(-ll) 6.78(-15) 1.23(-14)  Table 2.5: Standard error for the numerical approximation of u" = sin(7ra;) for the four Q D M 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  N u m e r i c a l tests a n d results  In this section, we give three analytic examples to study the properties of the Q D M . 2.4.1  O n e - d i m e n s i o n a l 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 F D method. Since the intention of the Q D M is to use nonclassical base functions, in this section we applied four nonclassical polynomial sets based on four weight functions. They are w ( x ) = sin (7roj), 2  x  w>2(x)  — e , w z ( x ) = e~ , x2  x2  and 104(2) = 1 — x . We refer to the Q D M based on these 2  weight functions as Q D M 1 , 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. A l l 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 Q D M competes well with that for the classical Legendre method and the method is much faster than the F D method (see Table 2.1). We also study the largest eigenvalue of the derivative operators based on the nonclassical Q D M , the result is similar to the classical methods, that is, the largest eigenvalue of the first derivative matrix grows like 0(N ) 2  and of the second derivative matrix grows like  0{N ). 4  2.4.2  A n a n a l y t i c example  In order to demonstrate the applicability of the Q D M and some of its properties, we consider a time-dependent Fokker-Planck equation defined by  where (2.4.2) The physical interpretation of P(x,t) is the probability density function. It is defined in x £ (—00,00), and satisfies P ( ± o o , £ ) = 0 . With a ^-function initial condition at xo, Eq. (2.4.1) has an analytic solution  where z = e / . The equilibrium solution as t —>• 00 is given by 1  4  (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 solution  PQ{X).  1. Eigenvalue  problem.  First, we consider the eigenvalue problem of Eq. (2.4.5): -L<f> (x) n  =  [- i )-^)  + ( )^}M )  A x  B  X  = -KM ),  X  X  (2-4.6)  where A and <f) are the eigenvalues and eigenfunctions, respectively. The exact eigenn  n  values of Eq. (2.4.6) are A = n/4, n  7i = 0 , l , 2 , . . . .  (2.4.7)  Three set of polynomials are applied in the Q D M 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 polynomials with weight function  w(x) = P o ( x ) .  w(x) = e~ ' , 4x 2  and nonclassical  We refer to the methods based on these  polynomials as Legendre, Hermite, and Q D M , respectively. The discretization matrix L N of the differential operator in Eq. (2.4.4) is in the form - L  N  = - A A + BA  2  (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 L N of the operator L can be obtained in a simple form L  = BDD , ' T  N  (2.4.9)  where D is the modified Q D M 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)  N  Ai  A  Exact  0.25  1.25  5  48  Aio  Al5  A20  2.5  3.75  5.0  .319(+02) .616E+01) .654(+00) .919(+00) .574(+00) .285(+00)  .458(4-02) .103(4-02) .171(4-01) .168(4-01) .141(4-01)  .124(+02) .140(4-01) .230(4-00) .152(-01) .767(-06) .439(-12)  .213(4-02) .358(4-02) .100(4-01) .227(-01) .577(-05)  .808(4-01) .151(4-01) .402(-02) .315(-01) .312(-05) .276(-ll)  .153(4-02) .369(4-01) .559(4-00) .759(-02) .240(-05)  Error Hermite 4 6 10 15 20 25 30 40 50  .218(-01) .123(-02) .592(-06) .205(-12) .104(-13) .905(-14) .133(-13) .508(-13) .143(-13)  4 6 10 15 20 25 30 40 50  .205(-02) .467(-04) .530(-08) .722(-15) .999(-15) .000(4-00) .666(-15) .555(-15) .355(-14)  4 6 10 15 20 25 30 40 50  .154(-03) .257(-05) .347(-07) .158E-08) .140(-11) .934(-13) .433(-14) .278(-15) .472(-14)  .465(+01) .524(+00) .682(-02) .889(-04) .146(-07) .184(-08) .389(-12) .913(-13)  .186(+02) .278(+01) .126(+00) .504(-01) .192(-01) .681(-04) .777(-07) QDM  .204(+00) .242(-02) .308(-05) .212(-07) .206(-ll) .422(-14) .666(-14) .000(+00)  .576(+01) .282(+00) .104(-01) .506(-04) .453(-07) .151(-13) .444(-14) MQDM  .388(+00) .216(-01) .180(-06) .999(-09) .107(-11) .102(-13) .666(-15) .688(-14)  .306(+01) .408(+00) .447(-02) .193(-03) .163(-07) .431(-13) .160(-13)  Table 2.6: Numerical error in the eigenvalues of the analytic example for the Hermite and Q D M methods.  Chapter 2. The Quadrature Discretization Method (QDM)  N  Exact Error  Ai  A  0.25  1.25  5  Aio  Al5  A20  2.5  3.75  5.0  .207(4-02) .823(4-01) .453(4-01) .268(4-01) .672(4-00) .523(4-00)  .386(4-02) .174(4-02) .113(4-02) .466(4-02) .137(4-01)  .606(4-01) .128(4-01) .796(4-00) .845(-02) .109(-03) .102(-05)  .127(4-02) .118(4-02) .511(4-01) .440(-00) .728(-00)  .133(4-01) .257(4-01) .229(4-01) .194(4-01) .135(4-01) .841(4-00)  .424(4-01) .343(4-01) .317(4-01) .250(4-01) .183(4-01)  Legendre 4 6 10 15 20 25 30 40 50  .295(-02) .102(-03) .143(-06) .105(-11) .156(-13) .354(-10) .491(-08) .462(-04)  .845(4-01) .259(4-01) .978(4-00) .395(4-00) .854(-01) .128(-02) .414(+00) .341(+00) .583(4-00)  .250(-01) .262(-02) .468(-06) .200(-12) .247(-14) .164(-14) .239(-14) .711(-14) .799(-14)  b  .761(+00) .292(+00) .186(4-01) .928(4-00) .424(4-00) .656(-01) .199(-03) .128(-09) .695(-12)  .568(-01) .566(-05) .143(-07) .770(-ll) .134(-12) .222(-12)  Legendre 4 6 10 15 20 25 30 40 50  a  .500(+00) .230(+00) .192(-01) .140(-02) .125(-04) .242(-06) .691(-03)  Legendre 4 6 10 15 20 25 30 40 50  49  0  .105(+00) .313(-01) .297(-04) .137(-09) .927(-14) .255(-14) .555(-14) .389(-15) .178(-14)  .901(+00) .739(+00) .494(4-00) .314(4-00) .465(-01) .942(-05) .443(-ll) .254(-10)  Legendre points in  .255(4-00) .168(4-01) .139(4-01) .109(4-01) .916(4-00) .388(4-00) .216(-05) a  (-2,2); (-3,3); 6  c  (-4,4).  Table 2.7: N u m e r i c a l error in the eigenvalues of the analytic example for the Legendre method.  Chapter 2.  The Quadrature  Discretization  Method  (QDM)  50  discretization as modified Q D M ( M Q D M ) . The main advantage of using the M Q D M representation is that the differential matrix is simple to construct. The matrix representative 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 A  1 ;  A , 5  Aio, A 1 5 , and A o is shown in Table 2.6 and 2.7 for each method. The exact eigenvalues 2  are given in the top of the tables. Table 2.6 compares the errors of these eigenvalues calculated by the Hermite, Q D M and M Q D M . As we can see from the table, Q D M and M Q D M 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 M Q D M is slightly slower than that for the Q D M . To obtain 14 decimal places accuracy of Ai it requires A^ = 30 points for the M Q D M in comparison with N = 15 points for the Q D M . The Hermite method calculates small eigenvalues (e.g. A A 5 ) very well 1 ;  but the results for the larger eigenvalues (e.g. A i , A o) are not as good. The eigenvalue 5  2  A20 is only accurate to 0(1) for the Hermite method with A^ = 50 in comparison with O(10~ ) for the Q D M and M Q D M . 5  For the Q D M , 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  Hermite, N=8  (QDM)  51  Hermite, N=32  0  0.5  • oo  O O O O O O O O O  o o ooaasao  -0.5 -1 -100 Legendre, N=8 0.5  ca E  -0.5  -50 Legendre, N=32  o o oo o o  o  CO  c  0  -2  -1  o o  Real  Figure 2.10: Numerical approximation of the eigenvalues for the analytic F P E problem.  Chapter  2.  The Quadrature  Discretization  N 4  6 10 15 20 25 30 40 50  Method  QDM  (QDM)  52  Hemite  (-1.834, (-2.077, (-2.375, (-2.605, (-2.764,  1.834) 2.077) 2.375) 2.605) 2.764)  (-2.885, (-2.983, (-3.136, (-3.252,  2.885) 2.983) 3.136) 3.252)  (-1.010, (-1.326, (-1.834, (-2.344, (-2.775, (-3.155, (-3.498, (-4.106, (-4.632,  1.010) 1.326) 1.834) 2.344) 2.775) 3.155) 3.498) 4.106) 4.632)  T a b l e 2.8: D o m a i n g e n e r a t e d f r o m t h e Q D M a n d H e r m i t e w e i g h t f u n c t i o n for t h e a n a l y t i c FPE  example.  t h e s o l u t i o n . W e w i l l discuss m o r e a b o u t t h e choice of t h e w e i g h t f u n c t i o n i n C h a p t e r 3 a n d C h a p t e r 4. T a b l e 2.8 gives t h e d o m a i n s for t h e Q D M a n d H e r m i t e m e t h o d versus n u m b e r of grid. T o a p p l y t h e L e g e n d r e p o i n t s , one has t o t r u n c a t e t h e i n f i n i t e d o m a i n ( —oo, o o ) i n t o a f i n i t e d o m a i n (—c, c). T a b l e 2.7 gives t h e n u m e r i c a l e r r o r of these eigenvalues c a l c u l a t e d by  L e g e n d r e m e t h o d i n t h e i n t e r v a l s a:  (-2,2), b :  (-3,3), a n d c: (-4,4), r e s p e c t i v e l y .  W e refer t o t h e a b o v e t h r e e L e g e n d r e m e t h o d s as L e g e n d r e " , L e g e n d r e m e t h o d , r e s p e c t i v e l y . A s we c a n see f r o m t h e t a b l e , t h e L e g e n d r e (-3,3) is t h e best a m o n g t h e t h r e e L e g e n d r e m e t h o d s .  6  6  and Legendre  c  method on interval  A s N i n c r e a s e , t h e e r r o r of t h e  eigenvalues c a l c u l a t e d b y t h e L e g e n d r e " m e t h o d decays i n i t i a l l y , a n d t h e n i n c r e a s e . T h i s is b e c a u s e t h e t r u n c a t e d d o m a i n is t o o s m a l l , a n d w i t h N i n c r e a s i n g , t h e t r u n c a t i o n e r r o r c a u s e d b y t h e b o u n d a r y c u t off a c c u m u l a t e s a n d is g e t t i n g l a r g e r a n d l a r g e r . T h e Legendre  0  m e t h o d is i n g e n e r a l slower t h a n L e g e n d r e  6  m e t h o d , since t h e d o m a i n (-4,4)  is w i d e r t h a n i t requires a n d some p o i n t s near t h e b o u n d a r y are n o t c o n t r i b u t i n g . I n c o m p a r i s o n w i t h t h e Q D M a n d M Q D M , t h e convergence of t h e L e g e n d r e  6  m e t h o d 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 Q D M and M Q D M , all the eigenvalues are real. This property is obviously important with regard to C P U 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^TV - ) for the 1  Q D M , 0(A^ - ) for the M Q D M , C^/V - ) for the Hermite, and 0(N ) 1  85  1  2  2  76  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(N )), 4  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 x as initial condition. The 0  application of the eigenvalue expansion stated in section 2.3 in the solution the the timedependent problem with x chosen as one of the quadrature points is straight forward. In 0  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 F P E problem. We set x = 1. Since it is very hard to discretize the ^-function numerically based on 0  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 Q D M , 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 A t has to satisfy the condition A t <  2/\  m a x  ,  where A  m a x  is the largest eigenvalue of the differential  Figure 2.11: Largest numerical eigenvalues for the analytic F P E example. The results with the Q D M , M Q D M , Hermite and Legendre are denoted by solid lined with asterisks, crosses, solid circles and diamonds, respectively.  Chapter 2. The Quadrature Discretization Method (QDM)  N 10 15 20 25 30 35 40 50 60 70 80 90  QDM 0.24 0.12 0.075 0.051 0.037 0.028 0.022 0.015 0.011 0.0085 0.0067 0.0054  Hermite 0.094 0.056 0.039 0.030 0.024 0.020 0.017 0.013 0.010 0.0093 0.0079 0.0070  55  Legendre 0.45 0.20 0.11 0.071 0.048 0.035 0.026 0.016 0.011 0.0084 0.0063 0.0050  Table 2.9: Maximum time step in the F E time integration.  matrix. Table 2.9 gives the maximum time step for the forward Euler method for the Q D M , Hermite and Legendre differential matrix. The numerical solution of Eq.  (2.4.1) solved by the forward Euler method with  A i = 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 Q D M approximates the exact solution very well. Among the three discretization methods in spatial dimension, the Q D M approximation 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 solution 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 F P E with N = 6. The solid lines are the exact solution. The plus signs are the results with the Q D M , 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 F P E with TV = 15. The solid lines are the exact solution. The plus signs are the results with the Q D M , 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) dt  =  l,dP{y,t) 4 dy l  ldP\y,t) 2 dy 2  (2.4.10)  The equilibrium solution of this equation is  dP(y,t) dt  =  ldP(y,t)  4  dy  ldP (y,t) 2  8  dy  2  The eigenvalues of its eigen-problem are rc/4, (n = 0 , 1 , N ) and the eigenfunctions are the normalized Hermite polynomials H (y). n  We solve the eigen-problem numerically by  the Q D M and M Q D M with w(y) = Po(y) and Legendre methods in the interval (-6,6). Since the eigenfunctions are polynomials, we expect accurate results in numerical approximations. Table 2.10 gives the convergence of the eigenvalues for the three methods. For the Q D M and M Q D M , 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 Q D M and M Q D M 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.  N  The Quadrature Discretization Method (QDM)  domain  Ai  A  5  Aio  Al5  A20  3.75  5.0  .755(-14) .933(-14) .933(-14) .355(-14) .444(-14) .933(-14)  .109(-12)  Exact 0.25  1.25  2.5  Error Hermite 4  (-2.020, 2.020)  .555(-16)  10 20 30 40 50 70 90  (-3.668, 3.668) (-5.550, 5.550) (-6.996, 6.996)  .222(-15) .111(-14) .178(-14)  .222(-15) .266(-14) .155(-14)  (-8.213, (-9.284, (-11.14, (-12.74,  .888(-15) .355(-14) .160(-13)  .799(-14) .666(-15) .622(-14)  .738(-14)  .533(-14)  4 10  (-6, 6)  8.213) 9.284) 11.14) 12.74)  .111(-13) .311(-14) .666(-14) .933(-14) .266(-14) .311(-14) .191(-13)  .400(-13) .320(-13) .258(-13) .977(-14) .888(-15)  Legendre .111(-15) .222(-15)  .205(-12)  .129(-11)  20 30  .833(-15) .350(-14)  .566(-12)  .141(-09) .317(-08)  .768(-10) .549(-07)  40 50  .173(-13) .533(-14) .279(-13)  .183(-09) .133(-09) .223(-09) .127(-10)  .647(-08) .331(-07)  .250(4-00)  .496(4-00)  .449(4-01)  .100(4-01)  .342(-06) .128(-05) .994(4-00) .979(4-00)  70 90  .192(-13)  .169(-12) .187(-12) .571(-13) .400(-09) .240(+01)  Table 2.10: Comparison of the numerical error of the eigenvalues of E q . (2.4.12) for H e r m i t e 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(N )). 4  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(N ) 2  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 ( i V ) , and 2  the condition number of the second derivative matrix is 0(N ). 4  However, the condition  numbers of the differential matrix representatives for the above examples are 0(N ) 2  O(N),  2.4.3  and  respectively.  Three-dimensional Poisson problems  In one dimensional examples, we already show that the Q D M provides very high accuracy and fast convergence in the solution of differential equations. In comparison with the F D method, the Q D M usually requires much coarser grid of points for the solution to converge to the same accuracy. This indicates a great advantage of the Q D M 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)  Boundary Condition  f(x,y,z)  A  3g(^+J/+2)  B  -207T ^(x, y, z)* + sin(47ra;) sin(27ry)a(z)*  - 5  ^  u\ c B  - 1  u(x,y,z)  (x+y+z)  exact solution  2  * <j>(x, y, z) = sin(47ra;) sm(2%y)[l - e  61  e  =0  ) ] , 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 considered. Numerical solution of example A is calculated by a second order central differencing F D method and the Q D M on a I B M RS6000 computer.  Two set of polynomials are  applied in the Q D M . One is the classical Chebyshev polynomials. Since the main intention of the Q D M is to use nonclassical polynomial basis set, we also applied a set of nonclassical polynomials associated with an arbitrary weight function  w(x)  =  e  6 x 2  .  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 different grids in different dimension in order to obtain superior convergence of the solution. The standard errors (E2) of the solution and the C P U 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 Q D M with nonclassical polynomials, the first C P U time is for the calculation of the points and the Q D M derivative matrix associated with the weight function w(x) = e  6 x 2  .  The second C P U time is for solving the discretized algebraic equation. As seen from  the table, the error of the solution calculated with the F D method decays at the rate  Chapter 2. The Quadrature Discretization Method (QDM)  FD •N  E  2  4 10 20 30 40 50 70 90  .754(-2) .143(-2) .366(-3) .163(-3) .920(-4) .589(-4) .300(-4) .182(-4)  62  Chebyshev CPU 0.01s 0.03s 0.39s 1.99s 6.75s 16.5s 61.1s 182.s  N  4 6 8 10 12 16  CPU 0.01s 0.03s 0.05s 0.10s 0.21s 0.50s  E  2  .486(-3) .263(-5) .731(-8) .136(-10) .122(-12) .522(-13)  N  4 6 8 10 12 16  E  2  .113(-2) .882(-5) .254(-7) .478(-10) .145(-12) .631(-13)  QDM CPU1 1.72s 1.76s 1.77s 1.81s 1.83s 1.86s  CPU2 0.00s 0.00s 0.01s 0.03s 0.08s 0.15s  Table 2.12: Standard Error and C P U time for the numerical approximations of example A of the Poisson problems.  of 0(1/N )  as N increases, whereas the error of the solution calculated with the Q D M  2  decays at a much faster rate. With N = 10, or 10 x 10 x 10 grid, E2 = 0.136 x 10~ for the Chebyshev approximation, and E  10  = 0.478 x 1 0 ° for the Q D M . For the F D _1  2  method, E = 0.589 x 10~ with N = 50, or 50 x 50 x 50 grid. To achieve accuracy of 4  2  O(10  -12  ) , only a 12 x 12 x 12 grid is needed, and it costs less than a second for both  Q D M and Chebyshev method. However, for the F D method, with 90 x 90 x 90 grid, it takes about 3 minutes and one can only obtain accuracy of O(10~ ). Comparing the 4  C P U time, the Q D M is much faster than the F D method if the same accuracy of the solution is required. Although the F D method performs relatively poorly in the calculation of Example A in comparison with the Q D M , the result calculated by it is still acceptable. In some cases, the F D method may fail while the Q D M 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 — e ) 5  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)  TV  FD E .218(02) .732(01) .452(01) .261(01) .194(01) .150(01) .108(01) .688(00) .349(00) .211(00) 2  10 16 20 26 30 34 40 50 70 90  CPU 0.04s 0.19s 0.44s 1.19s 2.20s 3.67s 7.15s 17.4s 36.5s 118.s  Chebyshev E .154(02) .121(00) .144(-02) .449(-06) .947(-09) .137(-11) .267(-12)  CPU 0.06s 0.18s 0.44s 1.18s 2.07s 3.49s 6.73s  2  63  QDM E .272(02) .375(00) .492(-02) .186(-05) .460(-08) .719(-11) .849(-12) 2  CPU 1.77s 1.89s 2.33s 3.22s 4.19s 5.59s 8.94s  Table 2.13: Comparison of the standard Error E and C P U times for the three numerical approximations of example B of the Poisson problems. 2  standard error of the numerical solution versus log N w  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 F D method converges very slowly. Even with a 90 x 90 x 90 grid, the F D solution barely reaches 1 decimal place accuracy. While for the Chebyshev method and the Q D M , the numerical solution can achieve 11 decimal place accuracy with 40 x 40 x 40 grid. For both examples, it is evident that the Q D M 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 0.80 1  1  1  1  1  1.20  1  1.60 log  1 0  (N)  Figure 2.15: Standard error of the solution for example B .  1  2.00  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 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 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 F P E 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]. A general Fokker-Planck equation in one dimension is of the form dP(x,t) dt  where P(x,t)  dA(x)P(x,t) ~  dx  d B(x)P(x,t) 2  +  d*x  •  [  6  A  A  )  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 applications may be characterized by two or more stable states. 66  Examples of systems  Chapter  3. Applications  of the QDM to the Solution  of Fokker-Planck  Equation  67  studied with a F P E 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 applications [87,91], plasma physics [88], nuclear dynamics [90], and numerous other applications. 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 F P E usually involves the expansion of the probability density function in a suitable basis set, and the reduction of the differential 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]. A n alternate approach involves the discretization of the P D F 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 pseudospectral approach as described by Fornberg [19]. Fourier series or Chebyshev polynomials are almost exclusively chosen as basis functions in the application 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 F P E 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 Q D M to the solution of several  Chapter 3. Applications of the QDM to the Solution of Fokker-Planck Equation  Fokker-Planck equations. distribution, P(x,t),  68  The F P E describes the relaxation of some non-equilibrium  to a "steady distribution", Po(x), at infinite time which can be  written in the form, (3.1.2) with (3.1.3) 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 Q D M , 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 F P E , it is anticipated that the present methodology will be easily adapted to nonlinear and multidimensional problems. The distinction between the F P E 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 Q D M will provide an accurate and efficient approach to these problems.  Chapter 3. Applications  3.2  of the QDM to the Solution of Fokker-Planck  Equation  69  The Q D M matrix representation of the Fokker-Planck operator  Equation (3.1.1) can be rewritten as dP(x,t)  d  dt  A(x)P(x,t)  dx  dB{x)P{x,t)  +  (3.2.1)  dx  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)  - A W ^ ^ l dx  dt =  +  B i x f dx  ^  9 2  (3.2.2)  -Lg(x,t),  It is easily shown that the operator L is self-adjoint with respect to the equilibrium solution Po(x). Direct application of the Q D M derivative matrix to the "spatial" operator L gives the matrix "representative" N  N  N  (3.2.3) J2 J2 A g{x )j=l j=l k=l which is not symmetric and its eigenvalues may be imaginary. It is therefore important Lgi = -A(xi)  Aij9j + (xi) B  3k  k  to construct a symmetric representative of the operator L. Let's consider the matrix representative of operator L in the transform space, L  r  where F (x) n  J  w(x)F (x)LF (x)dx, n  m  (3.2.4)  is the set of polynomials orthonormal with respect to the weight function  w(x), that is J w(x)F (x)F (x)dx n  =S .  m  nm  (3.2.5)  For convenience, it will be useful to introduce a second set of orthogonal functions, defined by Qn —  to a;  \  Po(x  F {x) n  (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 P (x)Q (x)Q (x)dx 0  n  = cw  m  (3.2.7)  The matrix elements of the operator L in the basis set {Q (x)} are given by n  L  nm  =  J P (x)Q (x)LQ (x)dx,  =  J  0  n  (3.2.8)  m  P (x)Q (x)(A(x)Q' (x)-B(x)Q' (x))dx. 0  n  m  m  Integrating by parts on Eq. (3.2.8), L  nm  =  -B(x)P {x)Q {x)Q' (x) 0  n  (3.2.9)  m  + J A(x)P (x)Q (x)Q' (x) 0  =  n  +  m  (B(x)P (x)Q' (x)Q _(x)dx, ff  0  n  -B(x)P (x)Q {x)Q' 0  n  + J (A(x)P (x)  m  + (B(x)P (x))')Q (x)Q' {x)  0  0  n  m  +  B(x)P (x)Q (x)'Q' (x)dx, 0  and using the facts that Po(x) = 0 at the boundaries and A(x)P (x)  n  m  + ( ) °( ) dB x P  0  x  =  Q f  or  all x on the defined interval, we obtain that L  = J P (x)B(x)Q' (x)Q' (x)dx.  nm  0  n  In terms of the polynomial basis {F (x)}  (3.2.10)  m  and weight function w(x), upon which the  n  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)]F (x)[~ n  + q(x)]F (x)dx,  (3.2.11)  m  where HK  w'(x) P^{x) = ^ - ^ - T . ' 2w(x) 2P {x)  V  (3.2.12) ;  0  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  L  nm  ~ E B{x )w [F' {x ) k  k=i  k  n  k  + q(x )F (x )][F (x ) k  n  k  m  k  + q(x )F (x )], k  m  k  (3.2.13)  Chapter 3. Applications  of the QDM to the Solution of Fokker-Planck  Equation  71  where x , k = 1,2, ...,N are the quadrature points which are the zeros of FN(X). k  matrix representative L  The  is transformed back to the physical space with transformation  nm  matrix T defined in Chapter 2, that is, T,- = y/wlF [xi) n  N  and we obtain  n  N  Lij — ^ ' / ~] Ti L T j, n  nm  (3.2.1.4)  m  n=l m=l  Substituting E q (3.2.13) and transformation matrix T into Eq.(3.2.14), that is, N  N  N  v = E  E  \/w F (xi)Y^B(x )w [Fl (x )+q{x )F (x )][F^ k=\  L  l  71=1 771 = 1  n  k  k  l  k  k  n  k  (3.2.15)  and using the fact that N  =8  w F (x )F (x,) k  n  k  n  (3.2.16)  %k  71=1  and N  J2 kF {x )F (x ) k=i  =6 ,  w  n  k  m  k  (3.2.17)  NM  one finds that N  ij ~ zZ B(x )[D  L  k  kl  + q(x )6 i][D j K  k  k  + q(x )S j], k  k  (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(x )D D j k=\ k  ki  k  (3.2.19)  The great advantage of the Q D M is that the matrix representation of the FokkerPlanck 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 F P operator. It is essential for the stability in the time integration of the O D E system ^ ^ .  -Y L {x ,t).  =  t  ii9  i  (3.2.20)  .7 = 1  3.3  E i g e n f u n c t i o n expansion  The formal solution of the Fokker-Planck equation, Eq. (3.2.1), may be evaluated by eigenfunction expansion. If <j) (x) are the eigenfunctions of operator L in Eq. (3.2.3), the n  eigenvalue problem of the F P equation is defined by L<f> (x) = -e <j> (x), n  n  n  (3.3.1)  where e are the eigenvalues and <f> (x) are the eigenfunctions. The linear time-dependent n  n  equation, Eq. (3.2.1), admits a solution of the form oo  = P (x) £  P(x,t)  0  c e-^Vn(x), n  (3.3.2)  where the coefficients c are determined from the initial condition, P(x, 0), and are given by n  c = J P(x, 0)(f) (x)dx. n  n  (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 e are real and satisfy e > 0, for n > 0. n  n  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 F P 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 F P equation. Their method corresponds to a standard finite difference scheme for solving time-dependent F P 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 F P 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 F P 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 F P 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 Q D M . First, the F P 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 Q D M 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 Q D M , the differential matrix of the F P 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 Q D M differential matrix is very easy to construct and the method is very easy to implement. The finite difference schemes seem more complex to implement. A t 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 F P E with bistable states. Our main objective is to determine the optimum basis sets with the Q D M 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) — x — x and B(x) = e with x € [—00, 00] . The equilibrium 3  P D F , 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 F P E . Blackmore and Shizgal [71] employed the Q D M , with a quadrature based on the steady equilibrium P D F 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 P D F w (x) = Po(x) = exp[—(x /4 — x /2)/e], 4  2  a  (b) a Gaussian weight function centered at the  origin, Wb(x) = exp(—x /2e), (c) a second narrower Gaussian weight function centered at 2  the origin, w (x) = exp(—4x /e), and (d) the sum of the equilibrium P D F and a Gaussian 2  c  weight function peaked at the origin, Wd(x) = PQ(X) + exp(—x /2.37e). The motivation 2  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 w (x) was used. a  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 P D F ,  w (x), a  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 A , A and A are triply degenerate [71] 3  4  5  and tend towards 2. Similarly, the eigenvalues A 7 , A and A9 tend towards 4. Hence, 8  with decreasing e the calculation of the splitting of these nearly degenerate eigenvalues becomes more difficult, as seen for A —> A . Results for e = 0.01 with the equilibrium 3  5  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  N  Ai  A  2  of the QDM to the Solution of Fokker-Planck  A  3  A € = 1.0 11.065 10.829 10.827  A  4  A  5  6  A  Equation  A  7  8  A  78  9  a  10 15 20 30 40 10 15 20 30 40 50 60  0.79216 0.79209 0.79209  3.5537 3.5489 3.5489  6 8623 6 8441 6 8436 6 8435  e = 0.01 1.8659 1.8658  0.96787 0.96786  15.910 15.445 15.415 15.415  2.6791 2.6776 2.6768 2.6765 2.6772 2.6772  1.7728 1.8708 1.8670 1.8670 e = 0.001° 1.9960 1.9880 1.9879  4 6 8 10 15  30.185 26.644 26.161 26.153 26.153  46.074 33.051 32.291 32.219 32.219  42.023 38.846 38.709 38.709  3.3328 2.9702 3.3050 3.3045 3.3044  3.4024 3.3652 3.3441 3.3321 3.3689 3.3655 3.3655  4.0276 3.8922 3.7107 3.5456 3.4526 3.4515 3.4515  3.9999 3.9529 3.9512 3.9512  4,0268 3.9547 3.9513 3.9512  6  1 8642 1 8645  1.7206(-3) 6.9871(-6) 6.0771(-8)*  23.238 20.607 20.543 20.540  2.0060 1.9882 1.9879 1.9879  e = 0.001  6  4 6 8 10  0.99698  6 8 10 15 20 30 40  1.0096 0.82578 1.0016 0.99953 0.99698  1 9880 1 9879  2 0220 1 9897 1 9912 1 9879  2.9727 2.9726 e = 0.001 2.0383 2.0024 1.9933 1.9879  3.9524 3.9513 3.9512  d  3.6293 4.1423 1.9900 1.9885 1.9879  4.0967 4.1748 3.1069 2.9723 2.9727 2.9726  17.0192 3.9524 3.9511 3.9512  18.3206 3.9546 3.9512  Table 3.1: Convergence of the eigenvalues for the quartic potential. b  w(x)  = eT^IM;  d  w(x)  = P { x ) + -* K™^). 0  e  solution of the Fokker-Planck equation.  P ( ) Q  X  =  w(x)  a  4.3375 4.0025 3.9541 3.9512 =  Po[x);  - * t ) I * \ the equilibrium s  Chapter 3. Applications  of the QDM to the Solution of Fokker-Planck  Equation  79  the e q u i l i b r i u m distribution (w ), a Gaussian centered on the origin (wb), and the sum of a  the equilibrium and a Gaussian at the origin (w ). W i t h the equilibrium distribution, only c  the eigenvalues A , A , A and A , are calculated (see F i g . 3.1: solid curves), whereas for 4  5  8  9  the Gaussian weight function centered at the origin, only the eigenvalues A , A , A , A and 2  3  6  7  A , are calculated. In both cases, the convergence is very rapid. This demonstrates that 1 0  the choice of weight function which determines the distribution of grid points is extremely important. T h e eigenfunctions <f> , 4  <TS , 5  (f>$ and <T5 i n F i g . 3.1 have their greatest variation 9  near x = ± 1 and the equilibrium distribution provides the required distribution of grid points for the rapid convergence obtained. B y contrast, the eigenfunctions </>, </>, (f> and 2  3  6  4> have their greatest variation near x = 0 and the Gaussian distribution centered at the 7  origin provides the required distribution of grid points for the rapid convergence of these eigenfunctions.  T h e eigenfunction pairs cf> and 3  CT5 , 5  and </3 and <f> , are antisymmetric 9  7  (see F i g . 3.1) and the corresponding eigenvalues are (nearly) degenerate (see Table 3.1). T h e final entry i n Table 3.1 is for e = 0.001 and a weight function which is the s u m 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 i n Table 3.2. Here, A ^ , Ni, N and Nd are the number of quadrature points associated w i t h the four c  weight functions shown at the bottom of the table.  T h e entries i n the table are the  approximate numbers of points required to calculate the eigenvalues shown to 5 significant figures.  As can be seen i n the results, the rate of convergence of the eigenvalues is  extremely rapid for e = 1.0 with both weight function w (x) and w (x), for e = 0.01 a  c  w i t h 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  Ai  A  e = 1.0 A N N N c  0.79209 15 50 20  3.5489 15 70 20  6.8435 30 50 20  10.827 20 70 15  € = 0.01 A N N  6.16(-12) 50 80  0.96786 80 15  1.8645 80 50  1.2016(-109)(°)  0.99698 ** 4 20  1.9879 ** 6 15  a  b  a  b  e = 0.001 A N N N a  b  d  A  2  A  A  Equation  80  Aio  Al5  A20  A25  15.415 20 70 20  45.599 40 90 30  85.415 40 110 50  133.03 50 140 40  187.38 60 160 50  1.8658 40 20  1.8670 70 50  3.9435 90 60  5.9608 100 60  8.7931 110 70  12.269 100 60  1.9879 8 ** 15  1.9879 8 ** 30  4.9234 ** 10 40  7.8013 20 ** 40  9.6867 ** 20 50  11.545 30 ** 60  3  4  5  ( )The value is estimated by Kramers' approximation. **More than 200 quadrature points are required for convergence. a  Table 3.2: Comparison of the rate of convergence of the eigenvalues for the quartic potential. N ,'Nb, N , and Nj, are the number of quadrature points for weight function W{ ) = P (x), W (x) = e-* /( ), W (x) = - / ( 0 . 2 5 e ) (j p (j -x /(2.37e) respectively. a  c  2  a  X  0  b  2£  2  c  2  e  a  n  d  Wd  x  =  Q  x  +  e  )  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 centered near the origin (as long as a proper width is selected). For e = 0.01, the eigenvalues converge faster with the equilibrium function  P (x), 0  i.e.  w (x), a  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 w (x) and for others with a  Wb(x).  But both weight functions fail to provide accurate results for other eigenvalues.  Chapter 3. Applications  e  Table 3.3: w{x) = P {x) 0  of the QDM to the Solution of Fokker-Planck  Ai  A  50 100 200 300 400  9.041867 13.04778 18.71484 23.06393 26.73058  QDM 30.45742 43.49627 61.94008 76.09384 88.02646  56.83964 81.02501 115.2331 141.4834 163.6139  50 200 200 300 400  9.04 13.1 18.7 23.1 26.7  Wu and Kapral 30.5 43.5 61.9 76.1 88.0  56.8 81.0 115 142 164  2  A  Equation  81  3  Comparison of the eigenvalues for large e for the quartic potential. = e x p [ - ( £ - 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  e 0.01 0.02 0.04 0.05 0.06 0.08 0.1 0.15 0.2 0.25 0.3 0.4 0.45 0.5  yQDM  ^Kramers  6.1596(-12) 1.6232(-6) 8.0874(-4) 2.7761(-3) 6.3146(-3) 1.7774(-2) 3.3545(-2) 8.2136(-2) 1.3478(-1) 1.8717(-1) 2.3802(-l) 3.3407(-l) 3.7934(-l) 4.2294(-l)  6.2518(-12) 1.6776(-6) 8.6901(-4) 3.0331(-3) 6.9792(-3) 1.9779(-2) 3.6951(-2) 8.5024(-2) 1.2897(-1) 1.6560(-1) 1.9564(-1) 2.4095(-l) 2.5828(-l) 2.7303(-l)  ^Kramers  _  82  yQDMyyQUM  0.0150 0.0335 0.0752 0.0926 0.1052 0.1128 0.1015 0.0352 -0.0431 -0.1152 -0.1781 -0.2787 -0.3191 -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 Q D M results with the equilibrium distribution as weight function and N = 20. For these large values of e, the Q D M 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 equilibrium 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 P D F 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 Q D M . 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. A << A < A a  2  3  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  (3.5.2)  e  With the explicit expression for U(x) for the quartic potential, we get that, A « -e'^ .  (3.5.3)  4t  a  7T  A comparison of the converged Q D M 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 A with decreasing e consistent with the form given by Eq. (3.5.3). Departures T  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, I , is a nonlinear function of the incident light intensity, I;. The relationship t  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 I versus t  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 F P E with coefficients  2C  A{x) = y — x  1 4- x2  (3.5.4)  '  and  2qx  2  (3.5.5)  (1 + z ) ' 2  2  T h e dimensionless input and output amplitudes are denoted by y and x £ [0, oo], respectively. The parameter q = C/2N  S  where N  s  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, E q . (3.1.2), is shown i n F i g . 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 b i m o d a l PQ(X) as a sum of two Gaussians. In this work, we choose w (x) a  = 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. T h e first ten eigenfunctions for C = 8, q = 0.4, and y = 8 and y — 9 are shown i n F i g . 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 i n F i g . 3.4 are localized in two regions; one is near x = 0 and the other near x = 6.6. O n the other hand, the equilibrium distribution (the n = 0 eigenfunction in F i g . 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. T h e eigenfunctions for y = 9 are more localized in the two regions than for y = 8. A comparison of the numerical convergence of the eigenvalues, A , with weight funcn  tions w (x) a  and Wb(x) is shown i n 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  N  Ai  A  of the QDM to the Solution of Fokker-Planck  A  2  3  A  A  4  y = 7.5 10 20 30 40  0.7779 0.5809  3 4958 0 9297 0 9291  5 4983 1 3063 1 3033  10 20 30 40 45 50  1.609(-5) 5.494(-7) 6.653(-7) 4.642(-7) 4.641(-7)  0.3773 0.3771  0 0 0 0 0  7078 6692 5255 5210 5209  6  A  7  A  8  A  88  9  a  6.7301 1.7289 1.6991  y = 8.0  A  5  Equation  9 2306 2 1979 2 1135  13.2532 2.7500 2.5445 2.5440  19.6935 3.7418 2.9903 2.9886  29.7897 5.4299 3.4510 3.4455  46.2350 5.7352 3.9434 3.9132  1 4857 0 9843 0 8822 0 8729 0 8728  2 1 1 1 1  3 7191 1 5440 1 3013 1 2742 1 2738 1 2738  7 3766 1 9160 1 5612 1 5206 1 5195 1 5195  16.4200 2.3760 1.8527 1.7880 1.7860 1.7858  a  1.0355 0.7342 0.6980 0.6976 0.6976  1358 2403 0626 0468 0467  y = 9.0"  2 4 6 8 10  0.6606 0.6605  1.3178  1.9735 1.9718  2.6224 2.6223  3.2768 3.2693  3.9130 3.9126  y = 9.0  b  6 8 10 20 30 35 40  0.1721 0.04931 0.1482 0.2064 0.2063  0.6613 0.6606 0.6606 0.6605  1.3510 1.3219 1.3172 1.3178 1.3178  40.5944 2.0198 1.3946 1.4933 1.4849 1.4849  5.0508 2.0090 1.9718 1.9718  3.3048 2.6223 2.6223  21.3499 2.9576 2.9039 2.8388 2.8388  3 2698 3 2693 3 2693  3 9180 3 9126 3 9126  Table 3.5: Convergence of the eigenvalues for the optical bistability model with q = 0.4. w(x) = PQ(X). w(x) = Po(x) + 0.33Po(0.33x + 6.6). Po(x) is the normalized equilibrium solution of the Fokker-Planck equation. a  b  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  w (x), a  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 constructing the Q D M 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,  and  w (x) a  Wb(x),  the rate  of convergence of the eigenvalues is more than twice as fast as that for weight functions w (x) c  = x exp(—x ) and Wd(x) = x exp(—x ) + x exp[—(x — 6.6) /0.1)] used in the 2  2  2  2  2  2  previous work by Blackmore et al [70]. For q = 0.4 and y = 7.5 and 8, the convergence is the same with both w (x) c  w (x) a  and  Wb(x).  The so-called "speed" quadrature points with  gives poorer convergence. For y = 9, w (x) is superior for some of the eigenvalues a  and worse for others. For example, for y = 9 and q = 0.4, the eigenvalues A , A and A T  4  7  converge with 90 to 110 points, while eigenvalues A , A , A , A , A and A converge with 2  3  5  6  8  9  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 w (x) = Po(x) at x = 6.6 dominates and, with small N , most of the a  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  Ai A N N N  a  b  c  A N N N  a  b  c  A N Nb N N a  c  d  A N N N  a  b  c  A N N  a  b  N  c  A N N  a  b  N  c  A  of the QDM to the Solution of Fokker-Planck  A  A  0.5809 20 20 30  0.9291 30 30 40  1.3033 30 30 40  A y = 7.5 1.6991 30 30 50  0.3771 20 20 60  0.5209 50 50 40  2/= 8.0 0.6976 40 40 70  q = 0.4  4.642(-7) 40 40 60  y =9.0  q = 0.4  0.2063 90 ' 30 30 30  0.6605 4 20 100* 60  1.3178 4 20 100 70  1.4849 100 30 30 40  1.9718 6 20 100* 90*  0.6883 20 20 30  1.1146 20 20 40  y = 7.5 1.5723 20 20 40  q = 1.0  0.3037 20 20 30  0.3343 30 30 40  0.4922 30 30 40  y = 8.0 0.7460 30 30 50  q = 1.0  1.514(-3) 20 20 40  0.6592 4 20 80  1.3099 6 20 100  y = 9.0 1.8157 60 30 30  q = 1.0  0.3166 50 40 30  2  3  4  5  q = 0.4 2.1135 30 30 60  0.8728 50 50 70  2.0519 30 30 40  1.0443 30 30 50  1.9516 6 30 100  A  A  6  7  Equation  As  A  90  9  2.5440 40 40 70  2.9886 40 40 80  3.4455 40 40 80  3.9132 40 40 80  1.0467 50 . 50 70  1.2738 50 50 80  1.5195 50 50 80  1.7858 50 50 90  2.6223 6 20 100* 90*  2.8388 110 40 30 50  3.2693 8 30 100* 90*  3.9126 10 30 100* 90*  2.5446 30 30 40  3.0443 30 30 50  3.5542' 4.0872 30 30 30 30 50 60  1.3761 30 30 60  1.7367 30 30 60  2.1225 30 30 60  2.5324 30 30 70  2.5835 8 40 100  3.2045 10 50 100  3.2588 70 40 100  3.8136 10 50 100  * Converge to essentially 3 figures.  Table 3.6: Comparison of the rate of the convergence of the eigenvalues for the optical bistability model. w (x) = P (x); w (x) = P (^) + 0.33P (0.33x + 6.6); w (x) = x e~ \ w (x) = x e-* + 2 -(x-6.6)Vo.i_ a  b  a  d  2  d  x  e  0  c  b  o  o  2  c  x  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 P D F of the temperature is then given by Eq. (3.1.5) with B(T) = q /2 and 2  A{T) = Q[l-a(T)]-toT\  (3.5.6)  l-a(r) = i,r<T ,  (3.5.7)  with  7  l-a(T) = l-a(T)  =  1  7 o  7  + /?r,T < x < T ; 1  2  2,T>T . 2  The coefficient A(T) is the difference between the solar influx Q[l — a(T)] (Q being related to solar constant, taken to be 3 4 0 W m  -2  with a the albedo) and the infrared cooling rate,  ecrT , (e being the emissivity and o the Stefan constant). In Eq. (3.5.7), 70, 71 and 72 are 4  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~ , T = 297.0, 7, = 61.2, 8  2  72 = 0.75, and for the continuity of A(x), 70 and T\ are defined by 70 — 72 — /?T and 2  Ti - (71 - 7o)//3, respectively. For appropriate values of these parameters the system is bistable, and the corresponding climate potential has two minima at T = T± and a maximum at T = T which lies 0  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: P (T) 0  = exp(-- U(T)), 2  (3.5.8)  Chapter 3. Applications  10 20 30 40 50 60 70  Ai 1.3990(-4)  A 0.57262  3.6683(-5) 3.3949(-5) 3.3889(-5) 3.3888(-5)  0.51358 0.49091 0.48983 0.48981 0.48980  2  of the QDM to the Solution  A 1.0695 0.65862  A 1.1826 1.0430 0.94902 0.94177 0.94160 0.94160  3  4  0.61516 0.61350 0.61347 0.61347  A 1.8395 1.1819 1.1703 1.1687 1.1687  of Fokker-Planck  Aio  Al5  3.5430 2.9259 2.7497 2.7335 2.7331  7.6004 5.9582 5.2499 5.0674 5.0466 5.0458  5  80 90  Equation  A20  A25  10.019 8.6076 8.0389 7.8278 7.7945 7.7931  15.110 12.757 11.684 11.045 10.865 10.824 10.822  92  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) = — J  T  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 F P E numerically to determine this long time dependence. 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 y r K . With an increase in 8, the steady _ 1  state moves from one peak at T  +  2  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 Q D M 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 P (T) for the climate model, q = 7.13 y r K . 0  - 1  2  Chapter  3. Applications  of the QDM to the Solution  of Fokker-Planck  Equation  function. In Fig. 3.6, we show the first ten eigenfunctions for q = 7.13yr  _1  94  K and 0 2  = 0.007123, corresponding to the symmetric potential and symmetric stationary distribution 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 log (Ai) versus 0 for several q values is shown in Fig. 3.7. The decrease in 10  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 y r ^ K , and 0 = 0.007123. 2  Chapter  3. Applications  of the QDM to the Solution  of Fokker-Planck  Equation  96  Figure 3.7: The variation of the lowest nonzero eigenvalue, A for the climate model. The solid line is the converged result with the Q D M . The dashed line is the Kramers' approximation, Eq. (3.5.2). a  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 distribution 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 y r K , and 8 = 0.007123 with S function at x = 1.5 ( at high _ 1  2  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 distribution 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 Q D M and other numerical methods of solution of FokkerPlanck 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 y r K , and 8 = 0.007123. The initial distribution is a delta function at T = 301 i f corresponding to x = 1.5. _ 1  2  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 distribution 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 A and Ai. 0  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 T d . The solid lines represent the exact results. The circle symbols are the results with the Q D M 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 Q D M .  Chapter 3. Applications  3.6  of the QDM to the Solution of Fokker-Planck  Equation  102  Summary  In the present chapter, we have demonstrated that accurate eigen-solutions of different Fokker-Planck equations can be determined with the Q D M . 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 Q D M . We have demonstrated that the convergence of the eigenfunctions 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 quartic 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 Q D M . The slow approach to equilibrium is then characterized 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  A p p l i c a t i o n s of the Q D M to the S o l u t i o n 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 solution 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 distribution 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 x e x p ( — x ) on the interval [0, oo]. 2  2  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 Q D M 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) = y +  (4.1.2)  2  i  +gy  This potential is chacterized by the parameters A and g. It is close to a harmonic oscillator 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 unnecessary. 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 — 2g , 2  then E = l - 2 g and for A = -7g -6g±g^25g 2  1  2  - 12g + 4, then E = (9g + X)/g. These 2  results are useful for benchmarking different numerical methods. Hautot [105] reconsidered 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)  = y  6  - 3 y  (4.1.3)  2  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 S W K B results [114,115] and an exact calculation 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)= -y 1  2  (4.1.4)  + 2y + -y 4  i  6  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)  = y + ey 2  4  (4.1.5)  which has been studied by several workers. Banerjee et al [102] and Banerjee [103] employed 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] obtained 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 Q D M 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 Q D M 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  T h e equivalence of the F o k k e r - P l a n c k eigenvalue p r o b l e m and the Schrodinger equation  The eigenvalue problems of Fokker-Planck equation are equivalent to Schrodinger equations: Recall the eigenvalue problem of F P E presented in Eq. (3.3.1). If the independent variable, x, is transformed to a new variable, y, (4.2.1)  and we define ip (y) by n  My)  =  \/Po{x{y))Mx(y)),  (--) 4  2  2  Chapter 4.  Applications  of the QDM to the Solution of Schrodinger  with Po(y) = yB( (y))Po( (y))i x  equation  107  then the Fokker-Planck eigenvalue equation, Eq. (3.3.1),  x  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 eigenvalue 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 authors [160]. It can be shown that the time-dependent Fokker-Planck equation equivalent to this stationary Schrodinger equation is given by dP(y,t)  _  dt  d dy  W> y)P(y,t) -rm d  {  +  (4.2.6)  with stationary solution given by P ( y ) = e x ( - j W(y)dy). 0  P  The Fokker-Planck eigenvalue problem is given by ,d<j)n d cf) 2  u  /  (  n  (4.2.7)  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 {S (y)} n  H  nm  orthogonal with unit weight function, that is, = - J S (y)S' (y)dy n  + J S (y)V(y)S (y)dy.  m  n  (4.3.1)  m  With an integration by parts in the first integral, we have that H  = J S' (y)S' (y)dy  nm  where V  nm  n  = J S (y)V(y)S (y)dy. n  +V ,  m  (4.3.2)  nm  We now introduce a second polynomial set,  m  {F } n  orthogonal with weight function w(y), that is, S (y)  where w(y) = exp(— / W(y')dy').  n  Eq. (4.3.2) can be rewritten as in  in'  1  /  (4,3.3)  = y/^y)F {y),  n  + 7 T ^ n + ^ F \dy + V .  (4.3.4) If one of the cross terms in the integrand above is integrated by parts one gets that ™(y)[K  F  n  = J wF' F' dy  H  nm  F  n  m  + [Km -  Vnm],  nm  (4.3.5)  where V(y) = \w\y) If the matrix representative H  nm  -\w\y).  (4.3.6)  is transformed back to the discrete representation in  physical space with the transformation T, that is, N  N  Hij — } ] } ] Ti H T j, n  n=0 m=0  nm  m  (4.3.7)  Chapter 4.  Applications  of the QDM to the Solution of Schrodinger  equation  109  one finds that N  Hi, = E kiDkj  + [V{ ) - V(y ))S .  D  yi  t  (4.3.8)  i3  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, d  d  2  2  ~Q^~]py~  1 ^nm(x,y)  ( >)  +  V  X  y  (4.3.9)  = X ^n {x,y) nm  m  the eigenfunctions are represented by a two-dimensional grid constructed from the product space of orthogonal polynomials in x and y. The matrix representative of the twodimensional Hamiltonian for bases set {X (x),  orthogonal with unit weight func-  Y (y)}  n  m  tion is given by 8  d  2  = - JJ X ,{x)Y ,{y)(— n  2  +—  m  )X (x)Y (y)dxdy n  m  + 11 X (x)Y ,(y)V(x,y)X (x)Y (y)dxdy nl  m  n  (4.3.10)  m  W i t h an integration by parts, we obtain that  -hf ' i n  where  m  V 'm\nm n  nm  - If  Sm'mJ  X' ,(x)X' (x)dx n  n  +S > n  Y ,(y)Y (y)dy+  J  n  m  (4.3.11)  V i,nm,  m  nlm  X (x)Y '{y)V(x,y)X (x)Y {y)dxdy. nl  m  n  m  As for the one-dimensional case, consider polynomial sets {G (x)}, {E (y)} n  m  orthog-  onal with weight function u(x),v(y) respectively, that is, X (x)  = ^u(x)G (x)  Y (y)  = yJv~ty~)E (y)  n  m  (4.3.12)  n  m  (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' ,(x)G' (x)dx n  n  + $' J  v(x)E' ,(y)E' (lj)dlJ  n n  m  m  +[V  ],  n  (4.3.14)  where J(\u\x)^-U'{x))u(x)G' {x)G' {x)dx nl  +*n'n  and U(x),V(y)  n  /{\v {y) - \v\y))v{y)E' ,{y)E' {y)dy,  (4.3.15)  2  m  m  satisfy U'(x) = -}nu(x),  (4.3.16)  V'(y) = -\nv(y).  (4.3.17)  We obtain the Q D M representation of the Hamiltonian by transforming matrix to the discrete representation Hij i, tk  Hij,ki = hi E k'iD ' fc'=i fc'=i D  k 3  + Sij  i/r ' ,m'm n  n  that is, Dk'kDkn + [V(xi, y ) - V{x , y )]8ij8 i, k  t  k  k  (4.3.1.8)  where V(x, y) = \u\x) and N , and x  4.4  - \u\x)  + \v\y)  -  -V\y).  l  (4.3.19)  are the numbers of quadrature points chosen in actual applications.  Applications and results  The Q D M has been applied to several one-dimensional Schrodinger equations by Shizgal and Chen [51]. In this work, the Q D M 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 Q D M , 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, F (x), n  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 considerable 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 Q D M to the two-dimensional Henon-Heles potential. The first potential that we have chosen and which has been studied extensively [101][129] is the N P O model, Eq. (4.1.2), shown in Fig. 4.1 as the solid curves. The dashed curves are the harmonic potentials, V(y) = y + X/g, for A = g = 100 and A = g = 10 2  (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(—ay ), where 2  a is a scaling parameter. For this N P O 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) = y + ^ . A and g equal to (a) 10, 10 (b) 100, 100 and (c) 10, 100. The dash lines are the corresponding harmonic potential V(y) — y + X/g. 2  1  2  gy2  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 convergence. The Q D M 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 Q D M representative of the Hamiltonian, 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 empirically 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  N  Ai  A  2  A= 0 5  1.00000000  3.00000000 A= 1  3  A  4  A  5  g=_0 5.00000000  7.00000000  9.00000000  g=l 5.62484751 5.58979112 5.58977905 5.58977894 5.58977893 5.58977893  7.83801690 7.64836479 7.64820406 7.64820127 7.64820124 7.64820124  10.32756366 9.68550819 9.68407574 9.68404264 9.68404202 9.68404202  13.61840145 13.38898345 13.38834923 13.38832431 13.38832353 13.38832349 13.38832349  16.71237123 15.82571826 15.81924074 15.81888806 15.81887214 15.81887152 15.81887149 15.81887149  53.83672948 53.83909110 53.83909383 53.83909326 53.83909327 53.83909326 53.83909326  64.45752724 64.18782502 64.18745791 64.18744198 64.18744105 64.18744100 64.18744100  10 20 25 30 35 40  1.23249101 1.23235080 1.23235072 1.23235072  10 20 25 30 35 40 45 50  2.78256744 2.78233128 2.78233044 2.78233052 2.78233052  A = 10 7.41859167 7.41750446 7.41750609 7.41750588 7.41750590 7.41750590  g= 1 10.73364911 10.70106074 10.70102615 10.70102563 10.70102557 10.70102558 10.70102558  9.35966852 9.35941813 9.35941803 9.35941803  A = 100 26.70397902 26.70596477 26.70596566 26.70596563 26.70596563  g = 1 41.44872496 41.44110330 41.44109963 41.44109976 41.44109975 41.44109975  10 20 25 30 35 40 45  3.51099389 3.50738872 3.50738837 3.50738835 3.50738835  A  Table 4.1: Convergence of the eigenvalues with V(y) = y + ^ where a is chosen for the fatest convergence. 2  114  n • w\(y) =  e x  p(  — a  !/ ) 2  5  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  N  Ai  10 30 50 60 70 80 90 100  1.11770702 1.05932983 1.05929698 1.05929690 1.05929689 1.05929688 1.05929688  10 30 50 70 80 90 100 .110  10 50 70 90 100 110 120 130  A  2  A  A  3  A= 1 3.54168906 3.08883073 3.08809133 3.08809085 3.08809085  g = 10 6.69900727 5.09038673 5.08285715 5.08284796 5.08284769 5.08284768 5.08284768  1.65877686 1.58013523 1.58002278 1.58002235 1.58002233 1.58002233  A = 10 4.53929108 3.88195452 3.87904292 3.87903684 3.87903683 3.87903683  g = 10 8.04585051 5.85711306 5.83286153 5.83276776 5.83276755 5.83276753 5.83276753  g = 10  5.82541635 5.79394465 5.79394241 5.79394231 5.79394230 5.79394230  A = 100 12.16555870 11.57221790 11.57219684 11.57219677 11.57219678 11.57219678  15.97213490 13.62913696 13.62877371 13.62877143 13.62877142 13.62877142  T a b l e 4.2: C o n v e r g e n c e of t h e eigenvalues w i t h w h e r e cv is c h o s e n for t h e fatest c o n v e r g e n c e .  V(y)  As  4  11.03978681 7.13612019 7.09048160 7.09037430 7.09037053 7.09037041 7.09037041  16.69475171 9.26980877 9.08892124 9.08805809 9.08801960 9.08801815 9.08801810 9.08801810  13.12551490 8.03082593 7,90413992 7.90315755 7.90315433 7.90315417 7.90315416 7.90315416  19.66357857 10.30455803 9.88876928 9.88233330 9.88230079 9.88229884 9.88229873 9.88229873  22.13479362 15.99309324 15.98848089 15.98843454 15.98843423 15.98843421 15.98843421  29.83816388 17.99876164 17.97250413 17.97208972 17.97208598 17.97208565 ' 17.97208562 17.97208562  y  2  +  115  h(y)  v  = exp(-m/  2  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  N  10 30 50 100 150 170 180  10 50 100 150 170 180  10 50 100 150 170 180  A  Ai  A  1.74034331 1.04840358 1.01083691 1.00841233 1.00841061 1.00841060 1.00841060  A= 1 6.61289499 3.33049421 3.04192292 3.00987806 3.00983181 3.00983177 3.00983177  g = 100 14.45434406 6.20396393 5.19309999 5.00989049 5.00927636 5.00927557 5.00927553  25.73416848 9.96961039 7.64468189 7.01481472 7.00985617 7.00984578 7.00984517  40.57425497 14.75340147 10.57282857 9.03654508 9.00959190 9.00949511 9.00948856  2.12557689 1.09321568 1.08408954 1.08406343 1.08406335 1.08406335  A = 10 8.03895659 3.19606870 3.09891916 3.09831922 3.09831722 3.09831706  g = 100 17.66231586 5.54663516 5.09892856 5.09279892 5.09276616 5.09276332  31.55531547 8.42645412 7.13621223 7.09883559 7.09850083 7.09846755  49.82378424 11.98070336 9.24759634 9.09763231 9.09530236 9.09503285  2.92175390 1.84742726 1.83638157 1.83633621 1.83633594 1.83633590  A = 100 9.34160581 4.11049745 3.98422018 3.98310435 3.98309903 3.98309857  g = 100 19.45189904 6.47955464 5.93857806 5.92841712 5.92834037 5.92833282  34.28680807 9.57139193 8.04492347 7.98535022 7.98458485 7.98449794  53.64778578 13.28590654 10.17201242 9.95499695 9.95023642 9.94960676  2  A  3  A  4  Table 4.3: Convergence of the eigenvalues with V(y) = y + where a is chosen for the fatest convergence. 2  1 +  5  ^ . w\(y) 2  = exp  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  N  X  1  10  1.23272180  30 40 50 60 70  1.23235100 1.23235074 1.23235073 1.23235072 1.23235072  10  2.78258502  20  2.78233137  30 40 45 50  2.78233053 2.78233052 2.78233052  10 15 20 25  9.35945915 9.35941761 9.35941803 9.35941803  A  2  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  3  g= 1 5.59128149 5.58978006 5.58977901 5.58977894 5.58977893 5.58977893  A  A  4  5  7.64562064 7.64819920 7.64820110 7.64820123 7.64820124 7.64820124  9.69097632 9.68404558 9.68404226 9.68404204 9.68404202 9.68404202  13.60393371 13.38888589 13.38832411 13.38832349 13 38832349  16.67001724 15.82484591 15.81888465 15.81887151 15.81887149  g=l 10.73118613 10.70105591 10.70102563 10.70102558 10.70102558  15 81887149  30  A = 100 26.70572641 26.70596964 26.70596558 26.70596563 26.70596563  35  g= 1 41.44628014 41.44114465 41.44110119 41.44109978  53 91385775 53 84147035 53 83917850 53 83909702  41.44109975  53 83909346  41.44109975  53 83909328  64 86511926 64 23043100 64 19022920 64 18763807 64 18745616  40 45  53 83909327 53 83909326  64 18744228 64 18744111 64 18744101  50  53.83909326  64 18744100  55  64.18744100  Table 4.4: Convergence of the eigenvalues with V(y) = y + ^ 2  1  u>i{y) =  1  exp(- yi + y  T^)-  2  • Weight functi  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  N  Ai  A  10 30 50 100 150 170 180  1.06515662 1.06003407 1.05946708 1.05930820 1.05929829 1.05929756 1.05929736  10 30 50 100 150 170 180  10 30 50 100 130 140 150 160 170  As  A  A= 1 3.08692234 3.08794408 3.08805698 3.08808859 3.08809057 3.08809071 3.08809075  g = 10 5.08700592 5.08337022 5.08296831 5.08285570 5.08284868 5.08284816 5.08284802  7.08838586 7.09012096 7.09031284 7.09036658 7.09036993 7.09037018 7.09037025  9.09234744 9.08856179 9.08814362 9.08802645 9.08801914 9.08801860 9.08801845  1.61407526 1.58268033 1.58046212 1.58003812 1.58002355 1.58002282 1.58002265  A = 10 3.87252286 3.87852707 3.87895272 3.87903382 3.87903660 3.87903674 3.87903677  g = 10 5.85545056 5.83441119 5.83303946 5.83277730 5.83276829 5.83276784 5.83276773  7.90700483 7.90231956 7.90301643 7.90314922 7.90315378 7.90315400 7.90315406  9.98909626 9.88401555 9.88258285 9.88230893 9.88229952 9.88229905 9.88229894  5.89164179 5.79569188 5.79404301 5.79394280 5.79394234 5.79394232 5.79394231 5.79394230 5.79394230  A = 100 11.65464995 11.57183960 11.57217532 11.57219667 11.57219677 11.57219677 11.57219677 11.57219678 11.57219678  g = io 14.22630311 13.62953798 13.62879829 13.62877155 13.62877143 13.62877142 13.62877142  17.92322840 15.99205694 15.98841237 15.98843409 15.98843420 15.98843420 15.98843421 15.98843421  22.47689280 17.99649871 17.97213192 17.97208577 17.97208563 17.97208562 17.97208562  2  Table 4.5: Convergence of the eigenvalues with V(y)  (y) = e x p ( - y y i T ^ J ^ ) . 2  Wl  118  As  4  = y + 2  1  _^ • Weight function gy2  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  N  Ai  10 30 50 100 150 170 180  1.00943021 1.00902768 1.00883028 1.00860665 1.00851699 1.00849611 1.00848760  A= 1 3.00981139 3.00981943 3.00982338 3.00982785 3.00982964 3.00983006 3.00983023  = 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  1.09402833 1.08994543 1.08798805 1.08583151 1.08499665 1.08480631 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  1.92323022 1.87989757 1.86200758 1.84549528 1.84038615 1.83936741 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  10 30 50 100 150 170 180  10 30 50 100 150 170 180  A  3  g  Table 4.6: Convergence of the eigenvalues with V(y) ™ (y) = 1  exp(- yi y  + -^-). A  TT  119  = y -f ^ 2  1  y 2  • Weight function  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  120  N  5 10 15 20 25 30  1.39754248 1.23347542 1.23235072 1.23235072  A= 1 6.00256900 3.50334632 3.50738845 3.50738835 3.50738835  00  1.23235072 1.23235353  3.50738835 3.50739706  g = 1 11.66691170 5.91911961 5.58977876 5.58977894 5.58977893 5.58977893 5.58977892 5.58983355  2.78138892 2.78233156 2.78233052 2.78233052  A = 10 8.72184392 7.41816173 7.41750593 7.41750590 7.41750590  g = 1 14.67163572 10.81174060 10.70102881 10.70102558 10.70102558  29.92249451 13.48916964 13.38872711 13.38832349 13.38832349  2.78233054  7.41767206  10.70448059  13.39000325  5 10 15 20 25 30  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  (0  9.35941803  26.70596563  41.44109975  53.83909296  (c)  5 10 15 20 25 30  27.12138709 7.73245910 7.64821025 7.64820124 7.64820124  52.46889173 13.42453838 9.68403519 9.68404205 9.68404202 9.68404202 9.68404195  7.64820121 7.64906899  55.15913549 18.62437460 15.82253275 15.81887215 15.81887149 15.81887149  91.26562732 64.45670875 64.18766807 64.18744157 64.18744100 64.18744100  Table 4.7: Convergence of the eigenvalues with V(y) = y + ^ 2 • Weight function W\(y) = e x p ( — O i i y ) / ( l + gy )" , where ot\ and a are given in Table 10. Results from (b) Fack and Vanden Berghe (1985) [124]; (c) Lai and Lin (1982) [118]. 2  1  2  2  2  2  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  A  N  Ai  10 12 15 20 25 30 35  1.06078000 1.05934121 1.05929700 1.05929690 1.05929688 1.05929688  10 20 25 30 35 40 45 (b) 10 20 25 30 35 40 45 50  A  2  A  3  A= 1 3.10750137 3.08925905 3.08810693 3.08809085 3.08809085  g= 10 5.34002953 5.10521689 5.08309896 5.08284789 5.08284767 5.08284768 5.08284768  1.67530513 1.58002638 1.58002232 1.58002233 1.58002233  A = 10 4.31996615 3.87916818 3.87903807 3.87903684 3.87903683 3.87903683  g= 10 9.75230993 5.83491979 5.83278723 5.83276775 5.83276753 5.83276753  1.58002233  3.87903683  5.83276752  7.48981433 5.79394731 5.79394193 5.79394232 5.79394230 5.79394230  A = 100 8.03655640 11.57682425 11.57225713 11.57219704 11.57219678 11.57219678  g= 10 44.09323078 13.70481854 13.62958358 13.62878139 13.62877147 13.62877142 13.62877142  = exp(—aiy )/(l + 2  gy 2) 012,  where  Fack and Vanden Berghe (1985) [124].  ct\  and  a  2  A  5  7.84826441 7.20306392 7.09455928 7.09037308 7.09037041 7.09037041  12.39506550 9.76442476 9.11338782 9.08806709 9.08801812 9.08801810 9.08801810  10.43741586 7.91869454 7.90346720 7.90315728 7.90315420 7.90315416 7.90315416 7.90315413  29.54884661 9.97591517 9.88457979 9.88233966 9.88229915 9.88229873 9.88229873 9.88229866  49.71269815 149.70360371 16.27467472 19.37491150 16.00289112 18.03283044 15.98861076 17.97411182 15.98843650 17.97211015 15.98843422 17.97208593 15.98843421 17.97208562 15.98843421 17.97208562  Table 4.8: Convergence of the eigenvalues w i t h V(y) wi(y)  4  121  are  = y 2 + ^ 2 • Weight function given in Table 10. (b) results from  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  N  Ai  A  2  1.04129101 1.00841165 1.00841060 1.00841060  10 20 30 40 45 50 55  1.45426410 1.08416996 1.08406336 1.08406334 1.08406334  A  3  = 1 3.29082677 3.00986701 3.00983203 3.00983177 3.00983177  g= 100 6.43712676 5.00988241 5.00928092 5.00927556 5.00927551 5.00927551  A = 10 5.01654653 3.10035001 3.09831782 3.09831700 3.09831700  g = 100 29.97472996 5.11504606 5.09277943 5.09276191 5.09276189 5.09276190 5.09276190  A  10 20 25 30 35 40 45  A  A  4  122  5  10.95519923 7.01509885 7.00992522 7.00984573 7.00984496 7.00984495 7.00984495  55.86316463 9.04372955 9.01020204 9.00949652 9.00948602 9.00948596 9.00948596  72.98922478 5.81559287 7.09864757 7.09844919 7.09844907 7.09844907  243.24844347 .7.21912972 9.09661478 9.09486638 9.09486470 9.09486466 9.09486466  A = 100 g = 100 .25273328 58.55206092 526.35484897 565.24148666 1785.60704400 1.89936536 4.43121455 8.05177942 12.63918879 27.70284628 1.83635795 3.98374475 8.03086862 5.93607533 10.15714058 1.83633587 3.98309869 5.92833557 7.98453363 9.95000108 1.83633584 3.98309834 5.92832858 7.98444358 9.94916197 1.83633583 3.98309834 5.92832857 7.98444352 9.94916096 1.83633583 5.92832857 7.98444352 9.94916096 1.83633444 3.98309836 5.92832790 7.90315413 9.88229866  10 20 30 40 50 60 65 0>)  Table 4.9: Convergence of the eigenvalues with V(y) = y + + 2 • Weight function i(y) P ( ~ ~ i y ) / ( l + 9y ) 1 where a\ and a are given in Table 10. (b) Results from Fack and Vanden Berghe (1985) [124]. 2  Y  1  w  =  e x  a  2  2  a2  2  GY  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  g/A  1  1 10 10 100  (1,10) (1.4,6) (1.4,6) (2,6)  10 (1.2,10) (2,8) (2,8) (2.4,8)  123  100 (3,12) (2,14) (2,14) (2.5,16)  Table 4.10: ( a i , a ) used for Tables 4.7-4.9. 2  has suggested that a useful choice of weight function would be derived from the "superpotential" 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, w (y)=exp{-a y )/(l+gy ) . 2  1  1  2  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 a in the weight function. For all pairs of A and g, we obtain 2  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, A is converged to 9.950 with 180 quadrature points 5  whereas it is converged to 9.94916096 with 50 points in Table 4.9. This demonstrates the usefulness of the Q D M 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  QDM a b c d e f  500 100 =1 1.232350723406 2 782330515932 9.359418026324 21.65874769959 1.23235072 2 78233052 9.35941803 1.23235072 1.24213 2 78233054 1.23235353 9.35941803 2 782330 1.23237205 9.35941803 21.6587477 2 1.23235 78233 9.3594  QDM a b e f  1.059296880862 1 580022327392 1 58002233 1.05929688 1 58002233 1.05929700 1 5800249 1 58002 1.05929  QDM a b c e f  g = 100 1.008410597947 1 084063335494 1.836335833449 5.083683913501 1.00841060 1 08406334 1.83633583 1.83633444 1 08411 1.8411 1 0840543 1.0084106 1.8363850 5.0840857 1.00841 1 08406 1.8364  QDM  g = 500 1.001849154630 1 084063335494 1.18486023962  124  10  a c e  1.0018491  1 0184910  g= 10 5.793942300193 16.73274738223 5.79394230 5.793947 5.794  16.73919  1.92317625551  1.18451 1.1848632  1.92255 1.9232260  Table 4.11: Comparison of the results of Ai with V(y) = y + ^ • 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]. 2  1 +  2  Chapter  4. Applications  of the QDM to the Solution  of Schrodinger  equation  125  with values of c*i and a which are chosen empirically for different values of A and g. The 2  Q D M 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 Q D M results are far superior to the previous results. Figure 4.2 shows the variation of the error in Ai for the N P O model (g = A = 10) versus the number of quadrature points, n, for four different weight functions. A  e x a c i  is defined  as the eigenvalue converged to 14 significant figures calculated with the Q D M . 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. However, 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  0  '  1  40  1  '  1  80  1  1  120  1  160  n  Figure 4.2: Variation of the error in Ai, A A i = |Ai — A j | for the N P O potential versus the number of grid points n for different weight functions, (a) w(y) = exp(—y ), (b) w(y) = e x p ( - j / y i + (c) w(y) = exp(-5.8j/ ), (d) w(y) = f^; A = g = 10 r a c i  2  2  e  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  127  Figure 4.3: Ground state eigenfunction for the N P O potential with g = X — 100, N = 25, with different weight functions. (*) w(y) — exp(—y ), (+) w(y) — exp( —17?y ) and (o) (y) (i+g~y2)* • The solid curve is for the last weight function with N=140; (A) full scale; (B) small scale near the origin. 2  w  =  e  2  Figure 4.4: Ground state eigenfunction for the N P O 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 mechanics [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, w (y) = exp(- / /4)  (4.4.2)  4  2  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 Q D M 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 corresponds to the Hermite polynomials, W2(y) defined earlier with a = 5, and another given by, w (y) = e x p ( - y / 4 - 5y ) 4  2  2  (4.4.3)  The results with the three weight functions are shown in Table 4.12. The overall convergence 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 eigenvalue. 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.  N  Applications  of the QDM to the Solution of Schrodinger  Ai  A  A  3  equation  130  Aio  5  w (y) = e x p ( - 5 y ) 2  2  5 10 15 20 25 30 35 40 45 50  3.20578381 1.92166391 1.93541230 1.93548442 1.93548209 1.93548210  19.34421619 11.48428663 11.67877474 11.68098869 11.68097117 11.68097087 11.68097089  24.44773756 25.22960745 25.25435384 25.25461676 25.25460450 25.25460490 25.25460488  72 71 71 71 71 71 71 71  68872477 64137641 57368183 57923539 57902800 57903698 57903668 57903669  75 72 71 71 71 71 71  81549114 04071624 58445530 57920993 57903737 57903670 57903669  83 73 71 71 71 71 71  65936104 07062002 64923337 57992969 57904422 57903671 57903669  w (y) = exp(-2/ /4) 4  2  5 10 15 20 25 30 35 40 45  1.95003306 1.93549705 1.93548226 1.93548210  13.51720225 11.68815652 11.68108903 11.68097109 11.68097089  w (y) = e x p ( - y / 4 4  2  5 10 15 20 25 30 35 40 45  4.54466778 2.23089971 1.94701006 1.93570651 1.93548392 1.93548212 1.93548210  23.73100470 13.39054786 11.78606709 11.68371077 11.68099725 11.68097108 11.68097089 11.68097089  25.58769695 25.26571988 25.25463882 25.25460546 25.25460488  by ) 2  31.59035642 25.89623300 25.28051990 25.25493023 25.25460771 25.25460489 25.25460488  Table 4.12: Convergence of the eigenvalues for SE with V(y) = y — 3y . 6  2  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 orthogonalization which is subject to considerable round-off errors [75,143]. Kaluza avoids these numerical difficulties by using symbolic algebraic techniques in mathematics. For arbitrary weight functions, this analytic approach is not feasible, whereas the Gautschi algorithm is generally convergent. Braun et al [125] employed a spectral method of solution 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, w (y) = exp(-2y - y /2) 2  (4.4.4)  4  3  and determined the matrix representative of the Hamiltonian in the "polynomial basis" representation, Eq. (4.3.5) with V  nm  —V  nm  = 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 H  nm  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 H  nm  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  N  1 3 5 8 10 15 20 25 30 35 N  30 40 50 55 60 65 70 75 80 85 90 95 100  A  Ai  A30  ^20  310.4920471524 309.4993497820 309.4993484837 309.4993484837  848.8060217068 588.5806628599 566.4282265701 566.4026817440 566.4026355012 566.4026354734 566.4026354734  A48  MO  1346.579274312 947.4614543288 893.9968790569 872.0907745529 868.2562193165 868.1457422322 868.1452015357 868.1452006773 868.1452006767 868.1452006767  1597.421054106 1364.247596709 1248.445773964 1183.544197185 1149.943901457 1138.668487703 1137.541785229 1137.522672203 1137.522588690 1137.522588541 1137.522588541  2  2  y  7  1.00000000000000000 1.00000000000000000 15.832389169799 15.124216267224 40.6232236023546 15.118931530866 36.367167641896 66.261603950851 15.118929992544 36.343021051640 62.648395926012 15.118929986242 36.342716214160 62.356049424923 15.118929986242 36.342716212413 62.356028944861 36.342716212413 62.356028944603 62.356028944604 62.356028944604  Table 4.13: Convergence of the eigenvalues of even parity with V(y) = \y w() = exp(-2y -y /2). 3  132  4  + 2y + \y . 4  &  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 determined 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, A is converged to 8 3  significant figures with 20 points whereas 25 are required with scaled Hermite polynomials. The choice of weight function is clearly important for the rapid convergence of the eigenvalues. The Q D M is also applied to a two-dimensional Schrodinger equation with HenonHeiles potential, given by 1 d  2  2d  I ill+ 2d y  2  2  V(x,y)  where V(x,y) = | ( x + y ) + Xx(y — \x ) 2  of previous workers.  2  2  2  1pnl(x,lj)  =  \ l^nl(x,y), n  (4.4,5)  and A = A/0.0125, consistent with the choice  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 Q D M 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 Q D M matrix yields the eigenvalues and eigenfunctions. The convergence of the two basis sets in each variable is shown in Table 4.16.  N  x  and N are the numbers of quadrature points in x and y dimension, respectively. In the y  tables, the eigenvalues are labeled with the principle quantum number n and angular  Chapter 4.  Applications  N  of the QDM to the Solution of Schrodinger  Ai  A  equation  As  Aio  2.44917485 2.44917408 2.44917407 2.44917407  € = 10( ) 16.63595545 16.63591955 16.63592150 16.63592149 16.63592149  35.88068953 35.88506209 35.88517148 35.88517122 35.88517122  94.30085478 95.81165911 96.15949348 96.15623411 96.15626312 96.15626298 96.15626298  10 12 15 20 25 30 35 40  4.99945382 4.99941563 4.99941758 4.99941755 4.99941755  e = 100< ) 34.87447875 34.87402295 34.87398862 34.87398427 34.87398426 34.87398426  75.72914876 75.88739267 75.87689375 75.87700463 75.87700403 75.87700403  253.32604009 201.40793502 205.27637088 204.79428957 204.79476335 204.79477459 204.79477451 204.79477451  10 12 15 20 25 30 35 40 45  22.86146298 22.86161867 22.86160889 22.86160887 22.86160887  e = lOOOOW 160.68335404 160.68601691 160.68588347 160.68591272 160.68591261 160.68591261  350.84170426 350.38352262 350.43503532 350.43589703 350.43589621 350.43589622 350.43589622  1022.19210882 924.84691394 944.02953926 947.71986787 947.68562278 947.68596392 947.68596166 947.68596167 947.68596167  3  a  12 15 20 25 30 35 40  b  Table 4.14: Convergence of the eigenvalues with ^w {y) = exp(-l(h/ ); ^w {y) = exp(-6(h/ ). 2  4  2  A  V(y)  = y  2  + ey . 4  ^w (y) 4  = exp(  Chapter  4.  Applications  N  of the QDM to the Solution  of Schrodinger  Ai  equation  135  no  10 12 15 20 25 30 35  2.44917318 2.44917406 2.44917407 2.44917407  e = lQ(> 16.63603391 16.63593038 16.63592170 16.63592149 16.63592149  12 15 20 25 30 35 40  4.99941762 4.99941755 4.99941755  6 = lQQfc) 34.87397375 34.87398436 34.87398426 34.87398426  22.86160088 22.86160897 22.86160887 22.86160887  e = 10000<) 160.68728162 160.68596913 160.68591446 160.68591261 160.68591261  a  35.86694240 35.88380588 35.88516632 35.88517122 35.88517122  107.31938413 98.84824260 96.71828902 96.16096863 96.15625913 96.15626298 96.15626298  75.87733275 75.87701004 75.87700401 75.87700403 75.87700403  210.04422203 205.20313119 204.79819433 204.79477654 204.79477452 204.79477451 204.79477451  c  10 12 15 20 25 30 35 40  350.26068143 1055.94778633 350.42650209 973.50076873 350.43583997 952.12677503 350.43589612 947.72238259 350.43589622 947.68593841 350.43589622 947.68596166 947.68596167 947.68596167  Table 4.15: Convergence of the eigenvalues with V(y) = y + t y calculated by fitting weight function to ground state eigenfunction. ^w (y) = exp( — (y + 5y )/2); ^w (y) = exp(-(2y + 6y )); ^w (y) = exp(-(50j/ + 25y )). 2  4  4  4  4  4  2  4  4  2  2  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 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 , Ai 2, only 2>  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 functions. 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 Q D M 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, (4.4.6) 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 N = N = 50. With the new quadrature defined by the weight function given by Eq. x  y  (4.4.6), this convergence is achieved with N = N = 32. This demonstrated a significant x  y  saving of computational time in comparison with other researchers. As can be seen from  Chapter  4. Applications  Ny/N  X  4  of the QDM to the Solution  6  8  10  12  of Schrodinger  16  20  equation  24  28  A ,o (4th) 2  4 6 8  2.9692 2.9632 2.9631  2.9605 2.9566 2.9566  4 6 8 10 12 16  7.0182 5.9505 5.8784 5.8716 5.8713 5.8712  5.9903 5.8914 5.8311 5.8204 5.8198 5.8198  2.9601 2.9563 2.9563 A  13.1807 12.4574  11.0336 9.6935 9.0536 9.0499 9.0497  11.1240 10.0448 9.9003 9.8905 9  15.5873 14.9696  12.4479 11.8715 11.4840 11.4788  5.9092 5.8322 5.8242 5.8176 5.8171 5.8170  5.9089 5.8256 5.8192 5.8176 5.8170 5.8170  10.1869 9.1708 9.0347 9.0334 9.0330  13.2657 10.1688 9.1287 9.0264 9.0247 9.0243  12.6973 10.1686 9.1284 9.0241 9.0223 9.0218  12.2964 10.1686 9.1284 9.0240 9.0222 9.0217  12.2319 10.7096 10.1248 10.1240 10.1237 10.1237  12.1574 10.6113 10.0684 10.0441 10.0425 10.0425  15.4681 11.6740 10.4985 10.0643 10.0379 10.0361 10.0360  14.5918 11.5832 10.4785 10.0640 10.0375 10.0356 10.0355  14.5446 11.5749 10.4774 10.0639 10.0375 10.0355 10.0354  12.3695 11.5307 11.1145 11.1072 11.1055 11.1051  14.1919 12.1856 11.2580 11.0690 11.0642 11.0637 11.0636  13.5202 11.8637 11.1959 11.0547 11.0515 11.0512 11.0512  17.8647 13.2560 11.8527 11.1959 11.0532 11.0506 11.0499 11.0498  17.2548 13.1758 11.8502 11.1699 11.0531 11.0502 11.0498 11.0497  17.0406 13.1734 11.8489 11.1609 11.0531 11.0501 11.0497 11.0497  14.4322 13.0533 12.5890 12.3128 12.2177 12.2139 12.2138  16.2730 13.7976 12.8019 12.3091 12.1503 12.0693 12.0689 12.0689  26.0197 16.0380 13.4736 12.6324 12.2442 12.1477 12.0666 12.0660 12.0660  21.9893 15.8520 13.4400 12.6245 12.2176 12.1465 12.0655 12.0652 12.0651  21.0597 15.5594 13.4335 12.6207 12.2143 12.1463 12.0651 12.0651 12.0650  (44th)  M  A ,_ 4 6 8 10 12 16 20  (16th)  5 j l  5.9126 5.8677 5.8279 5.8179 5.8174 5.8174 A  4 6 8 10 12 16  2.9601 2.9563 2.9562  (54th)  9  11.2850 10.6346 10.5047 10.3447 10.3234 Aio.io (65th)  4 6 8 10 12 16 20 24  18.1812 17.7635  14.4481 13.4192 13.2653 13.2631  12.3214 11.9189 11.6239 11.6110 11.6095 Ai2, (80th) 2  4 6 8 10 12 16 20 24 28  25.3965 21.9153 21.3602  16.1730 15.9680 15.6562 15.5597  17.7483 14.3616 13.4822 13.3034 13.2932 13.2926  17.2961 13.5335 12.8665 12.4343 12.3997 12.3940 12.3932  Table 4.16: Convergence of the eigenvalues, A / , for the Henon-Heiles potential wi Hermite points". "Weight function, u(x) = e x p ( — x ) and v(y) = e x p ( — y ) . ni  2  2  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 N = N — 32. x  y  Chapter 4. Applications of the QDM to the Solution of Schrodinger equation  Ny/N  X  4  6  8  10  12  16  10.1687 9.1286 9.0250 9.0232 9.0227  12.7920 10.1686 9.1284 9.0241 9.0222 9.0218  12.2620 10.1686 9.1284 9.0240 9.0222 9.0217  12.2173 10.6143 10.0691 10.0420 10.0403 10.0403  11.7764 10.5439 10.0644 10.0379 10.0361 10.0360  14.7114 11.5794 10.4775 10.0639 10.0375 10.0355 10.0354  12.3709 11.5196 11.0714 11.0702 11.0697  14.1440 11.9700 11.1959 11.0553 11.0517 11.0514  13.2565 11.8518 11.1959 11.0531 11.0504 11.0498  17 13 11 11 11 11 11  1066 1733 8489 1608 0531 0501 0497  14.3895 12.9474 12.3391 12.1806 12.1705 12.1701  16.0422 13.4600 12.6293 12.2266 12.1469 12.0660 12.0654  23 15 13 12 12 12 12 12  8554 6786 4334 6209 2142 1463 0651 0651  20  139  24  A2,o (4th) 4 6 8  2.9616 2.9575 2.9575  2.9601 2.9563 2.9562  4 6 8 10 12 16  7.0082 5.9531 5.8413 5.8764 5.8761 5.8761  5.9439 5.8738 4.9878 5.8185 5.8179 5.8179  A ,i  (16th)  6  5.9089 5.8283 4.9865 5.8176 5.8171 5.8170 A ,8  (44th)  8  4 6 8 10 12 16  13.2.729 12.5442  11.1456 9.7190 9.0446 9.0436 9.0433  11.2597 10.0603 10.0034 9.9984 A  4 6 8 10 12 16 20  15.6786 15.0393  5.9089 5.8248 4.9865 5.8172 5.8170 5.8170  _  9  9  (54th)  11.3010 10.6767 10.5473 10.3703 10.3499  12.5191 11.9125 11.5413 11.5353  Aio.10 (65th) 4 6 8 10 12 16 20  18.2794  12.3221 11.9199 11.7027 11.6918  14.4451 13.4728 13.3615  A , 1 2  4 6 8 10 12 16 20 24  25.6830 21.9864  16.2757 15.9901 15.7582  2  (80th)  17.2781 13.6380 12.8445 12.4913 12.4590 12.4532  17.9498 14.4303 13.5128 13.3842 13.3790  20.5669 15.2330 13.4323 12.6180 12.2141 12.1462 12.0651 12.0650  Table 4.17: Convergence of eigenvalues, A /, for the Henon-Heiles potential with new points . ^Weight function, u(x) = e x p ( — x 4- §Ax ) and v(y) = e x p ( — y ) . n)  6  2  3  2  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 Q D M 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 HenonHeiles 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 Q D M 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  Feit et al. 3.9825 3.9859 5.8672 5.8816 6.9991 6.9996 7.6979 7.7371 8.8116 8.8154 10.0356 10.0359 10.5727 10.5907 11.1603 11.3253 11.7497 11.7525 12.3335 12.2771 12.7474 13.0310 13.0868 13.0800  n 3  1 3 -3 5. 3 -3 6 6 -6 7 3 -3 8 6 -6 9 9 -9 10 6 -6 11 3 -3 11 9 -9 12 6 -6 12 12 -12 13 3 -3  QDM 3.982 417 3.985 761 5.867 015 5.881 446 6.998 932 6.999 387 7.697 721 7.736 885 8.811 327 8.815 188 10.035 413 10.035 592 10.572 480 10.590 470 11.160 258 11.325 231 11.749 519 11.752 297 12.333 785 12.277 192 12.748 445 13.032 062 13.086 873 13.081 196  equation  141  QDM nc nc nc nc nc nc nc nc nc nc nc nc nc nc 11.160 259 nc nc nc 12.333 786 nc 12.748 520 13.032 065 nc 13.081 199  a  6  c  Table 4.18: Eigenvalues for the Henon-Heiles potential. Hermite qudrature points based on u(x) = exp(—x ) and v(y) = exp(—y ), with N = N = 50. New quadrature points based on u[x) = exp(—x + §A£ ) and v(y) = exp(—y ), with N = N = 32. nc indicates no change from the QDM° results. a  2  2  6  x  2  3  y  2  c  x  y  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 Q D M 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 F D methods, it is possible for the Q D M discretization to preserve the symmetry for the self-joint joint differential 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 Q D M in choosing weight function provides one the opportunity to achieve the best solution with the least computing work. The Q D M 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 advantage of the Q D M over the classical spectral method is that the Q D M 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 Q D M is equivalent 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 Q D M (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 M Q D M has similar properties as the Q D M . Furthermore, the M Q D M 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 threedimensional Poisson equations, the Q D M 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 Q D M is more accurate and converges at a much faster rate. In the solution of the three dimensional Poisson problem, the Q D M requires much smaller number of grid points and less C P U time than the F D method to converge to a given accuracy. In the solution of the model problem, an analytic time-dependent Fokker-Planck equations, the Q D M with the equilibrium solution as weight function converges faster than the classical Hermite and Legendre methods.  A l l the eigenvalues of the Fokker-  Planck discretization matrix calculated by the Q D M (and M Q D M ) are real. By contrast, spurious imaginary eigenvalues occur for the classical Hermite and Legendre methods. The condition number of the F P operator is less than or equal to 0(N ) 2  for the Q D M  and classical methods. The main objective of the Q D M 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  C h a p t e r 5.  S u m m a r y and future  work  144  matrix is symmetric. All the eigenvalues, except the zero eigenvalue, are real and negative. 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 Q D M  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 FokkerPlanck 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(N ), P  with 1 < p < 2, as N —> co.  C h a p t e r 5.  S u m m a r y and future  145  work  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 Q D M can work well for nonsmooth functions. The Q D M  is also applied to the one-dimensional Schrodinger equations with several  potentials. The Q D M  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 Q D M  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 Q D M 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 Q D M  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 Q D M is very efficient to the solution of the differential equations discussed in this thesis. The Q D M 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 Q D M the Q D M ,  gives one more flexibility in the solution of differential equations. With  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 Q D M  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 Q D M 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 Q D M , 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 Q D M 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 Q D M in the solution of the advection-diffusion equation is our main task in the near future. (v) A l l 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  C h a p t e r 5.  S u m m a r y and future  work  147  equation) has been carried out. Our interest in the future is to study how to apply the QDM (vi)  to the solution of the nonlinear differential equations. When compared with the F D 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 easily.  QDM  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, equations: 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: Kutta and general linear methods. (John Wiley & Sons, New York, 1987).  [6] E . Hairer, The numerical solution of differential-algebraic  systems by  Runge-  Runge-Kutta  methods. (Springer-Verlag, 1989).  [7] D. Gottlieb, M . Gunzburger, and E . Turkel, S I A M 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, S I A M 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 semiconductor device equations, [microform]. Thesis. University of British Columbia. (1997)  [14] M . A . Crisfield, Non-linear finite element analysis of solids and structures.  New York, 1997)  148  (Wiley,  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  its applications;  element analysis for composite structures. Solid mechanics  and  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 (Pineridge Press, Swansea, 1983).  transonics.  [21] E. Tadmor, S I A M 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, S I A M 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 C L . 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 A I A A 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  applications.  analysis  of spectral methods: theory and  (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, U K (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 Symposium 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, S I A M 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, S I A M 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, S I A M 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, I A M J. Numer. Analysis 6, 273 (1986). [66] D. Funaro, S I A M J. Numer. Anal. 24(5), 1024 (1987). [67] T. Tang, S I A M 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 . L i , 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 (Cambridge University Press, Cambridge, 1995).  Derivatives  [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, Oxford, 1963). [97 J. S. Chang and G. Cooper, J. Computat. Phys. 6, 1 (1970). [98 E. W . Larsen, C. D. Levermore, G. C. Pomraning and J. G. Sanderson, J. Computat. Phys. 61, 359 (1985). [99 E. M . Epperlein, J. Computat. Phys. 112, 291 (1994); Laser and Particle Beams 12, 257 (1994). [100  B. T. Park and V . Petrosian, Astrophys. J. Suppl. Series 103, 255 (1996); Astrophys. J . 446, 699 (1995).  [101 A . K . Mitra, J . Math. Phys. 19, 2018 (1978). [102 K . Banerjee, S.P. Bhatnagar, V.Choudhry and S.S.Kanwal, Proc. R. Soc. Lond. A 360, 575 (1978). [103 K . Banerjee, Proc. R. Soc. Lond. A 364 (1978) 265 [104 N . Bessis and G. Bessis, J. Math. Phys. 21, 2760, (1980). [105 A. Hautot, J . Computat. Phys. 39, 72 (1981). [106 F. M . Fernandez, G. A . Arteca and E. A . Castro, Physica A 122, 37 (1983). [107 R. N . Chaudhuri and B. Mukherjee, J . Phys. A : Math. Gen. 16, 4031 (1983). [108 G. A . Arteca, F. M . Fernandez and E. A . Castro, J. Math. Phys. 25, 932, 1984. [109 M . R. M . Witwit, Pramana, J . Phys. 43, 279 (1994). [110 M . R. M . Witwit, J. Computat. Phys. 129, 220 (1996). [111  M . R. M . Witwit, J . Math. Chem. 19, 75 (1996).  [112 R. S. Kaushal, J . Phys. A : Math. Gen. 12, 253 (1979). [113 N . Bessis, G. Bessis and G. Hadinger, J. Phys. A : Math. Gen. 16, 497 (1983). [114; ' R . Dutt, A . Khare and U . P. Sukhatme, A m . J. Phys. 56, 163 (1988).  [115 A. Comtet, A . D. Bandrauk and D. K . Campbell, Phys. Letters B 150, 159 (1985). [116 R. Dutt, A . Mukherjee and Y . P. Varshni, Phys. Rev. A 52, 1750 (1995). [117 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. B a c k 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, M A , 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, S I A M 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, Y u . 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, Princeton, N . J . , 1987); J . Binney and S. Tremaine, Galactic Dynamics, (Princeton University 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).  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo 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

Comment

Related Items