Wavelet based in eigenvalue problems in quantum mechanics

8 March 1996

CHEMICAL PHYSICS LETTERS ELSEVIER

Chemical Physics Letters 250 (1996) 485-494

Wavelet bases in eigenvalue problems in quantum mechanics Jason P. Modisette

a,b, Peter Nordlander a,b, James L. Kinsey a,c, Bruce R. Johnson

a.c

a Rice Quantum Institute, Rice University, Houston, TX 77251-1892, USA b Department of Physics, Rice University. Houston, TX 77251-1892. USA c Department of Chemistry, Rice University, tlouston, TX 77251-1892. USA Received 22 November 1995

Abstract The use of Daubechies' compact support wavelets in quantum mechanical eigenvalue problems is investigated. It is shown that these orthogonal multiresolution functions provide an efficient basis for systems in which the potentials vary strongly in different regions. 1. Introduction

The applications of wavelet analysis have grown steadily in number over the last few years, including areas such as seismology, still and motion picture image compression, signal analysis, and several others [I-9]. Much of this activity has come from the promise of wavelets as a natural tool in the analysis of phenomena where different levels of detail are required in different places. Wavelet analysis, like Fourier analysis, uses expansion functions with different characteristic scales. Wavelets also possess an additional degree of freedom over sines and cosines because they are localized in space. One may then economize in, for instance, the description of a digital image by expanding the pixel information in a wavelet basis and discarding those basis functions that are found to be unimportant based on the regions in which they reside. Similar considerations apply in solving a differential equation by expansion in a basis, where it is important for the sake of efficiency to use the smallest possible basis that will represent the solution with the desired accuracy. One frequently finds that the solution varies strongly in Elsevier Science B.V. Pil S0009-2614(96)00060-7

some places and slowly in others. This information, which can often be estimated in advance+ is exactly what is needed for the selection of a minimal wavelet basis set. In the case of the time-independent Schr~dinger equation, the regions where the wavefunction is appreciable and where it possesses nodes are strongly determined by the details of the potential function V, which may vary widely from case to case. Wavelet bases appear attractive as a general approach to calculating eigenvalues and eigenfunctions for a wide range of potentials. There have already appeared a few investigations into the use of wavelets in quantum contexts [10-14]. Particularly interesting from the viewpoint of electronic structure calculations is the work by Cho et al. [ 15], exploring the use of the 'Mexican hat' wavelet basis in density functional calculations. This is in the same spirit as the current work, though the bases used are different. In this Letter we present an initial investigation of the utility of Daubechies' wavelets [1,5] as general and flexible bases for the variational solution of the time-dependent SchrSdinger equation. Properties specific to this sort of compact support wavelet are exploited to

486

J.P. Modisette et ai. / Chemical Physics Letters 250 f!996) 485-494

allow rapid evaluation of Hamiltonian matrix elements. This is, however, a first exploratory calculation, and we do not attempt to compare computationai speeds with other state-of-the-art methods for solving this sort of PDE. Our ultimate interests lie primarily in their use for electronic structure and large-amplitude nuclear motion wave functions, although the results we obtain are not limited to such applications. There are several related advantages to using wavelets as expansion bases in quantum mechanics. As mentioned, the simultaneous localization of wavelets in both coordinates and momenta allows one to customize the basis set to provide resolution locally where it is required. This localization also imparts sparse and banded characters to the kinetic and potential energy matrices. Orthonormal wawflets can be used, eliminating the need to solve a generalized eigenvalue problem. Also, these basis sets are easily extensible to multiple dimensions by taking the tensor products of one-dimensional wavelets. It is remarkable that all of these characteristics may in principle be embodied in a single set of functions. There are many types of wavelets O(x), each derived from specific choices of a mother wavelet or scaling function ~b(x) (see below). For the family of compactly supported wavelets discovered by Daubechies [1,5], both ~b(x) and ~ ( x ) rigorously vanish outside of a finite interval. This particular choice leads to orthonormal bases, which we regard as a significant advantage in avoiding complications of linear dependency. Latto et al. [16] and Beylkin [17] have given methods to obtain the required matrix elements of differential operators in these bases, enabling exact construction of the kinetic energy matrix. Sweldens and Piessens [18] have also derived a quadrature method by which expansion coefficients of the potential function in the Daubechies bases can be numerically obtained. With the aid of this expansion, matrix elements of the potential between two basis functions are re-expressed as a sum of integrals over the product of three basis functions. These may be calculated by the methods of Latto et al. [16] and Dahmen and Micchelli [19]. Thus one has available the means to deal with quite general potentials. The simple harmonic oscillator is examined as a first example to demonstrate that quantitative accu-

racy can be obtained and to quantitatively determine the importance of smoothness of the functions ¢,(x) and 0(x). A more significant system which possesses both nearly singular (softened Coulomb) and multicenter (double well) behavior is then investigated. These characteristics are shared by problems of physical interest such as the electronic structure of molecules. In Section 2 the theory and definition of the wavelet bases is given, followed in Section 3 by a description of the methods required for construction of the Hamiltonian matrices. The harmonic oscillator and steep double well examples are demonstrated in Section 4. A discussion of future extensions of the present calculations and further applications of wavelet bases in quantum mechanical problems is given in Section 5.

2. Properties of wavelets The mathematical theory underlying wavelet analysis is far more involved than needs to be considered here. Our concern is only to describe the essential ideas of the decomposition into orthogonal functions of different scales and locations. Greater detail can be found in several places [1-3,5,7,8]. For our purposes, a 'wavelet basis' consists of a set of translated copies of a single scaling function, ¢~(x), and a set of scaled and translated copies of a wavelet, O(x). The fundamelital characteristic of the present wavelet basis is that its members can be represented as linear combinations of finite sets of scaled, translated, copies of the scaling function ~(x). Stating this mathematically, 0 ( x ) is required to obey L-|

,/,(x/2) =

E h,¢,(x- k).

(l)

k-0

and the wavelet 0 ( x ) is defined in terms of these translated scaling functions, L-!

¢,(x/2) = ¢T E g,¢,( x - k).

(2)

k=6

These two equations can be applied recursively to relate all the wavelets and scaling functions that we will be using, and will be referred to as the recursion relations.

J.P. Modisette et a l . / Chemical Physics Letters 250 (1996) 485-494

Eqs. (1) and (2) do not completely specify the functions ck(3c) and ~b(x). We also require that ~(x) be orthogonal to copies of itself translated by integer steps, o k ( x - k ) , that ~(x) be similarly orthogonal to translated copies of itself, and that ~k(x) be orthogonal to ~ ( x - k). Requiring that all functions are unit normalized, we have that f._+==~k(x)~b( x - k) d x ffi 6k,0, f+f~b( x ) ~ ( x

k) d x

8k.0,

(4)

(5)

Since these functions are of compact support, a finite number of them will overlap spatially even though their overlap integrals vanish. This means that ~b(x) and qJ(x) must possess nodes as opposed to, say, a simple (3aussian function. The compactly supported wavelets of Danbechies [1] obey Eqs. (1) to (5), but these requirements still do not completely specify the functions. Daubechies imposed the additional constraints that the lowest moment of ~ ( x ) be unity, and that the lowest M moments of ~b(x) vanish.

1,

dx~-0,

(6) O
1

(7)

The smoothness or regularity of O(x) increases with M, which is taken to be an integer greater than one. If all of the conditions discussed here are introduced into the recursion relations (1) and (2), it can be seen that L ffi 2M, where M is the number of vanishing moments of O(x). Furthermore, h k and gk axe not independent, but are related by g k = (--1) L-k hL_l_ k

!.

j

i

L=6

x)

".,,. ~..~/'~ ...~ .-

(3)

k) dxffio.

dx

487

(8)

To complete the picture, the h~ can be obtained as the solutions of simple linear equations. A tabulation of these coefficients up to L = 20 may be found in the book by Daubechies [5]. Although there is no analytic representation of the scaling functions and wavelets, they may be evaluated by the method of Daubechies [5] and their

I IA si+<,> L=10

,

io

,

i/ij-. -

Fig. I. Scaling function ~b(x) (solid line) and wavelet ~b(x) (dotted line) as a function of x for L = 6 (upper panel) and L = 10 (lower panel)

values tabulated on a very fine mesh. Since they are all translations and scalings of a single wavelet and a single scaling function, this has negligible effect on overall memory requirements. Therefore, reconstruction of the wavefunctions is no more difficult than with a Gaussian or STO basis, just substituting a lookup table for the analytic evaluation of the basis functions used in the latter methods. The reconstruction of the wavefunctions is unneccessary for the evaluation of Hamiltonian matrix elements, which can be performed explicitly in terms of only the wavelet expansion coefficients. Examples of tk(x) and O(x) are shown for the cases L = 6 and L = 10 in Fig. 1. As stated above, these functions are smoother (more highly differentiable) for higher L; L =-6 is the lowest order for u,hich tk(x) and O(x) possess continuous first derivatives, and L ffi 10 is the lowest order with continuous second derivatives. Fig. 2 illustrates for L = 6 both the scaling function 2-1/~ck(x/2) and the six functions & ( x - k ) from which it is composed by Eq. (l). The wavelet 2 - 1 / ~ ( x / 2 ) , though not shown, is constructed from these same functions using Eq. (2).

488

J.P. Modisette et aL / Chemical Physics Letters 250 (1996) 485-494

a basis consisting of wavelets at each level (the 'detail' functions) until the maximum level j = J is reached. For this level, one also wants to include the scaling functions ~k~. The final orthonormal set to be used is then

012345k

//~/,..,= 2-'/2

~(x/Z)

X

10

{{:},

0

Fig. 2. Spatial dependence of the scaling function 2 - t / ~ ( x / 2 ) and the six scaling functions functions ~b(x- k), k = I-6. from which it is constructed.

One may continue to recur, generating wider and wider scaling functions. This leads to the set of functions

i(x)

(2-Jx- k),

=

(9)

where j represents the number of recursion steps, k the number of translations from the origin by 2 j, and the factor 2 -j/2 maintains unit normalization. All of these functions have exactly the same shape, but differ in width, height, and location. From Eq. (1), they are related by

.....

Although the scaling function at the coursest level must be included in the basis, it will be referred to as the wavelet basis. In practical calculations, both the number of levels of resolution J and the range of translations k within each level must be finite. If we are only interested in a certain region of space, then k can be truncated at each level to include only functions that are non-zero in that region.

3. Hamiltonian matrix elements We now examine the solution of the time-independent Schr~dinger equation,

(12) I d2

H f T-I- V( x) =

2 ~dx -l- V( x),

(13)

L-l k'-O

At each level of resolution, we also define scaled and translated wavelets,

using a finite wavelet basis like the one described above. We need to calculate matrix elements of this Hamiltonian, which can be written as integrals of the form

L-I =

E

(ll)

k"O

An orthonormai basis may be generated [2] based on the fact that the set of scaling functions { ~ - '(x)} on level j - 1 (with k miming from - = to + =) is exactly equivalent to the combined set of scaling functions and wavelets {(~k~(x)), {0/(x)}} on level j. At the next level, 4~ +l and ¢,/+1 are both constructed only. from the subspace {0~/(x)}, and therefore do not contain components in the orthogonal subspace {0/(x)}. Thus {0/(x)} and {0/+ t(x)} are orthogonal sets. Similarly, the wavelets on level j + 2 are orthogonai to those on both level j and j + I. Starting from j - - 0, then, one can accumulate

which break into kineticand potentialparts. Similar considerations were encountered and described by Schult and Wyld [20] in their use of the Daubechies wavelets to solve the Burgers equation with shockwave-like phenomena (see also Refs. [21,22]). The additional complication we face here is that matrix elements of general potentials must be calculated. Considering first the kinetic energy matrix elements, we see by Eqs. (10) and (11) that all wavefunctions here can be obtained by recursion fron, the set of scaling functions on any level I below the minimum level of functions used in the basis. Thus, it suffices to calculate matrix elements of T between two such functions. By a simple change of variable,

J.P. Modisette et a l . / Chemical Physics Letters 250 (1996) 485-494

however, these may be related to matrix elements of T between scaling functions at level j = 0,

:: ~p(2-J'x-k) IdOl

(~b[ITI~:)---

_

x ~,(2-~' x - E) dx ,J.+=

I d2 I

489

Here ~- is a simple shift and the oJi are the weights for the quadrature of order r. This quadrature is exact for polynomials of order r, and efficient in general by virtue of being able to use overlapping ranges of samples of f ( x ) for a range of k. For j different from 0, we may use a simple change of variable again,

=

+.~

•

X t ~ ( y - k') d y = 2 - / ( ~ ° l T I alP,). (14) The evaluation of the latter matrix elements can be reduced to the solution of coupled systems of linear algebraic equations by the methods of Latto et al. [16] and Beylkin [17]. We have implemented this scheme and performed the necessary rccursions in Mathematica using an adaptation of a package by Cohen et al. [23] to extended precision. In this way, we have obtained all desired second derivative integrals in the wavelet bases for L = 6,8 . . . . . 20 to an accuracy of 14-15 decimal places. These integrals are stored in binary form and read from disk as required. Evaluation of the potential energy matrix elements is more complicated. With other common basis sets, the integrals equivalent to (14) are often evaluated by direct numerical quadrature. The existence of efficient quadratures for evaluation of the coefficients in an expansion of a function in terms of wavelets, and the fact that three-function product integrals can be calculated accurately, allow for a more efficient method of generating these matrix elements. This method requires an accurate expansion of the potential in the wavelet basis J

Z Z(

x) + Z(

j=O k

k

(is) The expansion coefficients here may all be evaluated by recursion from scaling function integrals. Sweldens and Piessens [18] have developed a quadrature scheme for evaluation of the projections onto the scaling functions of sufficiently smooth functions f(x). In our applications, this takes the form r

(~,°1f)= EtoJ(k+i-l-':). i=I

(16)

=f

y- k):<2'y) dy r

i=l

(17)

In practice, we frequently subsample by performing this quadrature on a level j that is 2-4 orders below our lowest scale (~.e. negative j) and then recur upwards in j to get the desired coefficients accurately. Sweldens and Piessens [18] have examined ways to reduce the subsampling in order to increase efficiency. The r = 1 case of Eq. (17) approximates the overlap of a scaling function with f by a single sample of the latter. Higher-j overlaps computed by recursion from this starting point have an error proportional to the cube of the spacing o~ the finest scale [18,24]. This is only accomplished by use of the appropriate shift ~" away from the grid points on that scale. The standard sampling approximation without such a shift (e.g. as used in numerical recipes [25], to name only one commonly used reference) is equivalent to a single-point quadrature without such a shift. The price to be paid is drastic, i.e. the error scales linearly with the smallness of the finest spacing rather than cubically. Armed with the ability to obtain the expansion coefficients of general potentials, we m,st then select which levels of resolution a~d which k-values for each level should be included in the expansion. This is very straightforward; one only needs to choose a numerical threshold and to eliminate candidate basis functions whose expansion coefficients fall below that threshold. The candidate functions can be pre-selected on the basis of their locality and scale for the physical problem at hand.

490

J.P. Modisette et aL / Chemical Physics Letters 250 (1996) 485-494

Once the potential expansion has been obtained to the desired accuracy, the quantum mechanical matrix elements of the potential are reduced to linear combinations of integrals over a product of three wavelets and/or scaling functions. Latto et al. [16] and Dabmen and Micchelli [19] give methods for obtaining integrals between three scaling functions (or their derivatives) on a common scale, called 'connection coefficients', from which all the other integrals may be obtained through the recursion relations. All but a finite number of these connection coefficients are zero due to the compact support of the scaling functions. The number of nonzero integrals increases rapidly with L, however, because of the increasing interval of support and the number of spatially overlapping functions. We have evaluated all necessary connection coefficients across seven scales of resolution up to L ffi 12 using Mathematica. We note that due to the large number of computational steps involved in the recursion, extended precision is required. It is computationally more expedient to calculate the integrals only once; while the number of such matrix elements that could possibly be required in all problems to be investigated is large, it can be greatly reduced to a minimal unique set by appropriate changes of variables and use of permutational symmetry in the integrals. Matrix elements of other operators (for instance the dipole operator) can be easily calculated following the same proceOure as used for the potential.

By the vanishing of the moments of 0 ( x ) (cf. Eq. (7)), we can see that, for each j and k,

( ~b~lV)= O.

(20)

Thus, the expansion in Eq. (14) reduces to just the contributions from the scaling functions. That being the case, we may take Y = 0 with no loss of generality. Using y as the variable and appropriately normalizing the scaling function, we have V:

olY-Yo ~ Zc~A -'/~ ~ * [ - - ' V - ) '

(21)

k

+®V( y)A -'/2 #~k[ T[ Y0- Yo ,) dy. c , " f_®

(22)

Eq. (21) is an exact expansion. This is a result of Strang [4], who pointed out that the scaling function exactly interpolates polynomials up to order M - I = L / 2 - i, i.e low-order powers of x can be expressed without error as a linear combination of the scaling functions whose support includes the point x. Furthermore, the coefficients c k are evaluated exactly by the quadrature method of Sweldens and Piessens [18], making the oscillator example atypi. cally easy to study. We use as a basis all #~t° for i k i < kmax and diagonalize the resulting Hamiltonian matrix. Minimizing the lowest eigenvalue (1/2) with respect to Yo and A for different values of kmax and L yields the results shown in Fig. 3. It is clear that quantitative accuracy can be obtained (approximately 14

4. Simple harmonic oscillator The simpi~ harmonic oscillator is a good benchmark for any new basis set since it has analytical solutions and its smooth potential is similar to those encountered in many types of common quantum mechanical calculations. Choosing unit mass and unit frequency, 1 d2

H=

_~ !o -,o

1

2 -'-'~ d y +.~y2,

"~"~.~-~..--".--..

(18)

the eigenvalues are simply n + I / 2 , n = 0, I . . . . . . . We take the argument x of ~ ( x ) and ~b(x) to be related to y by possible changes of origin and scale,

Y f Yo + Ax.

~ - - ~ L= 12

(19)

L=20

I0-" 0

I0

20

30

Max/mum Ikl Fig. 3. Error as a function of kmax for the calculated energy of the ground state of the harmonic oscillator for different L.

J.P. Modisette et a l l Chemical Physics Letters 250 (1996) 485-494

significant figures), providing a partial check on codes and integrals. However, it is also clear that higher L (smoothness of the basis functions) is very important for quicker convergence of the eigenvalue with respect to basis size for this problem. This is a potential disadvantage that must be balanced in smooth problems if one wishes to use sparse matrix techniques for Hamiltonian diagonalization (not done here) since higher L also means wider bands in the Hamiltonian matrix.

t

-2 I

-1 I

491

0

1

2

3

I

I

I

I

/

-5

-10

5. Steep double well To strongly test the flexibility of the Daubechies basis, we now examine a problem which simultaneously exhibits both nearly singular behavior and a multicenter character. This is provided by a double well potential with softened Coulomb attraction at

Fig. 4. Plot of Ihe double well potential, ~b(x), F..q. (23) as a function of x. The lowest six energy levels are indicated with horizontal bars. The fifth and sixth stales lie very close together. and have turning points well outside the range shown.

y = :i: I/2 V(y)=

-

]/( Y -- 1/2) 2 + a 2 1

~/(y+ 1/2) 2 +a:' '

(23)

where a = 0.1 (atomic units are used here). This represents an approximation to the one-dimensional one-electron two-center Coulomb problem which re~ quires no special considerations with regard to placement of quadrature points, yet which is still more singular than would ever be encountered in, e.g., pseudopotentials in electronic calculations or Morse potentials in nuclear calculations. The potential is shown in Fig. 4 along with the first six energy levels calculated by Numerovooley integration (an inefficient method for this application, but one capable of providing accurate energies). It is seen that the levels occur in pairs, with the lower pairs undergoing the largest splitting. The lower level in each case is even around the origin, the upper odd. One finds that the wavefunctions concentrate probability density away from the origin with increasing excitation. In the limit that there is very little density around the origin and the wells, the left and right portions of the wavefunction decouple so that the even and odd linear combinations

become degenerate. This potential represents a numerical challenge because of the combination of its localized detail around the wells and its slow asymptotic falloff. The variational method applied to this potential will obtain all of the lower eigenvaiues simultaneously, but the very lowest states are more localized ~nd thus of less interest. Therefore we focus on the convergence of the fifth excited state, which requires taking into account both the localized and extended aspects of the potential. We found, as expected, that the lower states all converged more rapidly than the fifth excited state. For the calculations, the wavelet basis with L ffi 10 was cho~n. The first step in the expansion of the potential is choice of an origin and spacing for the coarsest scale functions in the basis. The origin was taken to be at y = 0 without further consideration. The expansion coefficients of the scaling functions (but not the wavelets) on this level undergo a slow l/k decrease, so it is necessary to choose them to provide adequate coverage of the spatial region of interest. The spatial boundaries of the coarsest level of scaling functions were chosen to be +25 au, which were the limits of the mesh required for a separate numerical calculation of the fifth excited state with an accuracy oi l0 -s aa. A scaling function spacing of 0.5 au was then found to provide an adequate representation of the potential in the

492

J.P. Modisette et a l . / Chemical Physics Letters 250 (1996) 485-494

Table I Ranges of k used for the different levels of resolution in the calculation of the expansion of the double well model potential. Spacing on coarsest level is 0.5 au Type of function

,~

Resolution level j

Minimum k

5 5 4 3 2 2

-50 -!! -13 -15 -18 -3 -26 6 - 37 27

I

I 0 0

Maximum k

50 3 5 7 -5 10 -14

18 - 36 28

asymptotic regions (except at the very edges where the truncation of the expansion is importan0. Thus, a total of 101 scaling functions were symmetrically spread around the origin (Yo " 0) as the first step in the expansion of the potential. The wavelets, on the other hand, represent the detail between adjacent levels of scaling function expansions; consequently, the expansion coefficients decrease quickly in magnitude once the wavelets move out of the region where the potential exhibits significant detail. A numerical threshold of 10 -5 was chosen as a minimum for the expansion coefficients in order for the corresponding wavelets to be included in the potential expansion. This led to six different levels of resolution, from j ~ 0 (fines0 to j=5 (coarsest), with spacings 0 . 5 × 2 J-s au. The k-ranges are listed in Table I. It can readily be seen

that the finer level wavelets accumulate in the regions of the wells. The asymmetry between positive and negative k in the included ranges is a direct result of the asymmetry of the compact support wavelets. Using the potential expansion of Table 1, a series of variational calculations with different wavelet bases representing the wavefunctions were performed. The calculations consisted of generating and diagonalizing the Hamiltonian matrix for each choice of wavelet basis, as was described in Section 3. The convergence of the fifth eigenvalue E 5 was monitore~l with respect to inclusion of additional basis functions at each step, and each candidate basis function was retained in subsequent runs only if the change in the eigenvalue was approximately 10 -s or higher. The initial calculations included only scaling functions with spacing 0.5 an, but it was possible to recur to spacing 2.0 au, retain the smaller number of coarser scaling functions, and eliminate several of the wavelet functions generated. Thus, the coarsest spacing of the basis functions for the quantum wavefunction is 2.0 au, and the initial runs of Table 2 consist only of scaling functions on this level. Next, wavelets at the same spacing were then added until new members were far enough to either side of the origin that they failed to change the eigenvalue by the approximate convergence threshold. This was repeated at finer and finer levels of resolution until finer basis functions also fell beneath the convergence criterion. The results are summarized in Table 2. We see quick convergence of the accuracy with increasing

Table 2 Calculated energy of the fifth exciled state of the double well model ootential for increasingly detailed basis sets. Spacing on coarsest level

is 2.0 au

Type of function

Resolution level j

Minimum k

Maximum k

Cumulative basis size

Calculated energy

+ + + + + +

¢~ ~ 0 0 ~ 0 ~

7 7 6 5 4 3 2

-

13 !! 12 12 !I 11 14

!I 3 3 2 2 3 - 10

25 40 56 71 85 100

-0.17052602 -0.22122345 -0.24563155 -0.24798580 -0.24824769 -0.24825948

+

0

2

2

6

110

- 0.24825966

exact

- 0.24825962

.

J.P. Modisette et a L / Chemical Physics Letters 250 (1996) 485-494

resolution, culminating in a total of only 110 basis functions. It was not necessary to go to finer resolution than j = 2 with this threshold. This corresponds to a minimum spacing of 0.0625 au. Thus, a total of six resolution scales is required for the wavefunction, although these are not exactly the same scales as used for the potential expansion. This is to be expected since, e.g., the s-wave radial functions for a hydrogen atom are quite smooth compared to the Coulomb potential. Conversely, for smoother potentials, very excited wavefunctions may require relatively higher resolution due to the crowding of nodes. In any case, the current problem shows around an order of magnitude savings over the complete basis consisting of wavelets and scaling functions at the included resolutions and spanning the 5:25 au region. The converged wavelet result is lower than the exact result by a few digits in the eighth decimal place. This minor discrepancy is not surprising given the effects of truncation of the potential expansion, but is systematically improvable by tightening of the numerical threshold used. Given the finite nature of the potential and wavefunction expansions, one also expects a residual sensitivity of the calculated results to the scale and origin of the finest-level grid. The choice of optima~ ~rids is clearly an important question which needs to be addressed in a future investigation. We nevertheless have demo~ ",trated quantitative variational accuracy using a modest basis (for the challenging problem at hand) of orthonormal, compact support wavelet basis functions.

6. Conclusion This exploratory calculation has shown that the compact support wavelet families found by Daubechies are promising as quantum rrechanical bases for potentials with one or more centers of strong variation. One is able to place detailed basis functions where needed without the need for similar detail in other locations. This multiresolution advantage is obtained without the need to sacrifice orthonormality of the basis functions. The use of a quadrature scheme and methods for evaluation of wavelet integrals allows the construction of Hamiltonian matrices for general potentials. This has been

493

shown to allow quantitative accuracy for potentials ranging from the simple harmonic oscillator to a nearly singular double well problem. There are many other types of wavelets that are available for use. For instance, the Coifman wavelet family (calculated by Daubechies) [26] is also an orthogonal, compact support, set; however, it allows a one-point quadrature (instead of r points as used here) at the price of approximately 50% greater support length. Several further examples are mentioned in the recent review by Jawerth and Sweldens [27]. While applications such as image processing have driven much of these developments, they may also be of use in quantum mechanical problems. The determination of optimal wavelet families for such contexts will be an issue to be explored in future work. An important aspect of the wavelet approach is that it extends to higher dimensions. In Cartesian coordinates, one simply needs tensor products of the wavelets and scaling functions in the individual degrees of freedom. It is actually the promise of use in higher dimensions that has motivated the present investigation. In almost any problem where motion in two or more degrees of freedom are highly correlated, we expect the localized resolution of the wavelets to provide an important efficiency/. Oue of the intended areas of applications is density fuactional theory (within the local density approximation), for which an investigation using non-orthogohal wavelets has appeared [15]. Another area is in the large-amplitude variational calculation of molecular vibrations, which is pushing the envelope of today's computing power in its applications to four-atom systems (see, for example, Bramley and Carrington [28] and Antikainen et al. [29], and references therein). Much of the wavelet approach extends to (and is expected [20] to be, I~-~.~lelmtqort,~t in) higher dimensions. Yet another application of importance is the solution of time-dependent quantum problems [30]. There has recently appeared a first attack on propagation in time of (a nonlinear version of) the Schr'6dinger equation [14]. While many earlier propagation algorithms have utilized uniform grids in order to take advantage of the fast Fourier transform, the same multiresolution advantages offered by wavelets in principle apply here as in time-independent problems.

494

J J~. Modisette et a L / Chemical Physics Letters 250 (1996) 485-494

Acknowledgement This work was supported by the Robert A. Welch Foundation and National Science Foundation Grant CHE-9220278. Grateful acknowledgment is given to P.O. Wells, X. Zhou, C.S. Burrus and R.A. Gopinath for very helpful conversations.

References [l] !, Daubechies, Commun. Pure Appl. Math. 49 (1988) 909. [2] S,G. Mallat, IEEE Trans. Pattern Anal. Machine lntell. I I (1989) ~74. [3] J,M. Combes, A. Grossman and P. Tchamitchian eds., Wivelets: time-frequency methods and phase space, (Springer, Berlin, 1989) p. 315. [4] G. Strtmg, SIAM Rev. 31 (1989) 614. [5] I. Daubecl~i~. in: CBMS-NSF Series in Applied Mathematics No. 61. T~n lectures on wavelets (SIAM, Philadelphia, 1992) p. 357. [6] G. Sl~ng, Bull. (New Set.) AMS 28 (1993) 288. [7] R.A. Gopinath and C.S. Bunis, Technical Report No. TR9301, Dept. of Elec. and Comp. Eng., Rice Univ., Houston, "IX 77251 (unpublished). [8] JJ. Benedetto and M.W. Prazier, eds., Wavelets: ma~ematics and applications, studies in advanced mathematics (CRC Press, Boca Raton, 1994) p. 575. [9] R.O. Wells Jr., Technical Report No. TR94.12, Computational Mathematics Laboratory, Rice Univ., Houston, TX 77251 (unpublished). [10] T. Paul, in: Wavelets: time.-frequency methods and phase space, ads, J.M. Combes, A. Grossman and P. TehamitchinJt (Springer, Berlin, 1989) pp. 204-208. [11] P. Fischer and M. Defranceschi, Int. J. Quant. Chem. 45

(199~; ¢~t9.

[12] D. Permann and I. Hamilton, J. Chem. Phys. 100 (1994) 379. [13] Z. Li, A. Borrmann and C. Martens, Chem. Phys. Letters 214 (1993) 362. [14] L. Gagnon and J.M. Lina, J. Pilys. A 27 (1994) 8207. [15] K. Cho, T.A. Arias, J.D. Joannopolous and P.K. Lam, Phys. Roy. Letters 71 (1993) 1808. [16] A. Latto, H.L Rcsnikoff and E. Tenenbanm. Technical Report No. AD91078, Aware, Inc., 1 Memorial Drive, Cambridge, MA 02142, USA (unpublished). [17] G. Beylkin, SIAM J. Numer. Anal. 6 (1992) 1716. [18] W. Sweldens and R. Piessens, SIAM J. Numer. Anal. 31 (1994) 1240. [19] W. Dahmen and C.A. Mieehelli, SIAM J. Numer. Anal. 30 (1993) 507. [20] R.L. Schult and H.W. Wyld, Phys. Rev. A 46 0992) 7953. [21] A. Latto and E. Tenenbanm, Compt. Rend. Acad. Sci. ! 311 (1990) 9O3. [22] R. Glowinski, W. Lawton, M. Ravacbol and E. Tenenbaum, in: Computing methods in applied sciences end engineering, eds. R. Glowinski and A. Lichnewsky (Society for Industrial and Applied Mathematics, Philadelphia, 1990) pp. 55-120, proceedings of the Ninth International Conference on Cornporing Methods in Applied Sciences and Engineering. [23] 1". Cohen, K. Jack and Chen and M. Xu, ftp://hilbert. mines`coloredo.edu:pob/wavelets (unpublished). [24] R.A. Gopinath and C.S. Bun'is, in: Proceedings of the 1SCAS-92, San Diego (1992). [25] W.H. Press, S.A. Teukolsky, W.T. Wetterling and B.P. Finnnery, Numerical recipes (Cambridge Univ. Press, Cambridge, 1992). [26] G. Beylkin, R. Coifman and V. Rokhlin, Commun. Pure Appl. Math. 44 (1991) 141. [27] B. Jawerth and W. Sweldens, SIAM Roy. 36 (1994) 377. [28] MJ. Bramley and J. Tucker Canington, J. Chem. Phys. 99 (1993) 8519. [29] J. Antikainen, R. Friesner and C. Leforestier, J. Chem. Phys. 102 (1995) 1270. [30] R. Kosloff, Ann. Rev. Phys. Chem. 45 (1994) 145.