by Delmar N. S. Permann,
Ottawa, Ontario, Canada,
delmar.permann@home.com,
2000 Delmar N. S. Permann

keywords: Morse, oscillator, wavelet, analysis, anharmonic, harmonic, chemistry, physics

NOTE: This worksheet demonstrates the use of Maple V for calculating a wavelet analysis of the Morse oscillator motion in an anharmonic region, i.e close to dissociation of a single bond. This worksheet illustrates how Maple can now do calculations that were only in the realm of FORTRAN or C as little as a few years ago, such as in Permann[1]. Further, this worksheet demonstrates how breaking the information into smaller subsets decreases the calculation time required dramtically, i.e. in this case by a factor of five, with some changes introduced into the analyized data.

The author expects that this worksheet will only be used for teaching and educational purposes and not for commerical profit without contacting the author for a licenced agreement.

Introduction:

• Equations of Motion:
• Integration of Equations of Motion
• Wavelet Analysis
• Conclusion
• References:

Introduction:

The Morse [2] oscillator is used as a model for a diatomic molecule, to represent anharmonic motion, in both classical and quantum mechanical calculations. In this work only the classical model will be presented. The reader is referred to the references [1..3] for further information. Morse designated a potential energy using the equation . The x value represents the bond length, x=0 is the equilibrium bond length and a is a spectroscopic constant specific to each individual diatomic molecule or bond. A diatomic molecule, exhibiting behaviour determined by the Morse potential, dissociates at D_e which occurs at extreme bond expansion (x>>0) while the compression of the bond past the equilibrium point (x<<0) needs a large energy, both of which are illustrated in blue in the diagram below. The red lines indicate an energy of 80% and 99% of the dissociation energy, respectively. In contrast to the anharmonic Morse potential a harmonic potential ( ) is shown using a green line in the diagram below. In this work, D_e is to 1.0, as is the spectroscopic constant a.

> restart;

> with(DEtools):with(plots):

> Digits:=12;a:=1.0;D_e:=1.0;

> p1:=plots[pointplot](
[seq([j/100,D_e*(1-exp(-a*j/100))^2],j=-70..600)],
### WARNING: the definition of the type `symbol` has changed'; see help page for details
color=blue,style=LINE,symbol=POINT):

> p2:=plots[pointplot](
[seq([j/100,D_e*(1-exp(-a*j/100))^2],j=-70..600)],
### WARNING: the definition of the type `symbol` has changed'; see help page for details
color=blue,style=LINE,symbol=POINT):

> p3:=plots[pointplot](
[seq([j/100,D_e*(j/100)^2],j=-110..110)],
### WARNING: the definition of the type `symbol` has changed'; see help page for details
color=green,style=LINE,symbol=POINT):

> l1:=plots[pointplot]({[-0.7,0.99],[5,0.99]},color=red,style=line):l2:=plots[pointplot]({[-0.65,0.8],[2.3,0.8]},color=red,style=line):

> display({p1,p2,p3,l1,l2},view=[-2.0..8.0,0..1.1],
axes=BOXED,title=`Morse vs a Harmonic Potential`,labels=[x,E]);

Equations of Motion:

The equations of motion of a bond using the Morse potenial are obtained by differentiating the Morse potential with respect to x the bond length. This is easily done using Maple and yields and where the forcing frequency is set to 1.0 rad/sec., the forcing amplitude F is set to , C is a dampening coefficient that has been added to include the dissipation of energy by the surroundings of a system and here it is set to weak dampening at . The reader is referred to reference [1] for a deeper explanation of the Morse anharmonic behaviour and reasons for choosing the above mentioned values.

I ntegration of Equations of Motio n

The original work was done using IMSL FORTRAN routines to do numerical integration of the equations of motion of a bond under the influence of the Morse potential. The equations of motion were set out using FORTRAN and a GEAR numerical integration was done by selecting such in the FORTRAN program. A listing of the routines used is given in reference [1].

This work uses the Maple DEplot routine with the numerical integration method set to dverk78 which is a 7th-8th order continuous Runge-Kutta method (setting the tolerance to 0.00001) and it reproduces the data done using a fore mentioned FORTRAN routines. The initial value for x is 0 and the initial value for y is calculated using where E_int is the initial energy with respect to the dissociation energy D_e and in this work it is set to 0.99 of D_e.

> C:=1e-6;F:=6e-3;

> st := time():
morsedverk204:=DEplot({D(x)(t)=y(t),D(y)(t)=F*sin(1.0*t)\-C*y(t)-a*D_e*2*(1-exp(-a*x(t)))*exp(-a*x(t))},{y(t),\
x(t)},t=0..204.75,[[x(0)=0,y(0)=-(2.0*0.99)^(1/2)]],\
stepsize=.05,scene=[t,x(t)],method=dverk78,\
tolerance=0.00001):time() - st;

> morsedverk:=op(1,op(1,morsedverk204)):

> morsedverk[1],morsedverk[2],morsedverk[4096];

The data points are generated and saved into a variable instead of plotting directly. This way the data can then be plotted using the plots package pointplot and is also available to be processed with the wavelet analysis routine. The reader may also want to read data in using the dataread statement which can then be plotted using the plots pointplot package as well.

>

> plots[pointplot](morsedverk,title=`Data E=0.99, Digits=12, F=0.006, dverk78, tol=0.00001`,color=blue, view=[0..200,-1.0..7.0]);

>

There is a need to prepare the data set in a specific format to be fed to the wavelet analysis below. So using nops the total number of data points is determined and kept in the variable Ntotal for later use. The wavelet code expects that the first signal (data) is in f0 (f zero) and we read the data in the variable as noted in the do loop below.

W avelet Analysi s

The reader is referred to Daubechies [4] seminal work for a complete explanation of wavelets or the reader can go to the share library share/numerics/daub for a descriptive explanation given by Gelinas[5]. Wavelets are a way to analyse data that is not as nicely behaved as a Fourier analysis needs. A Fourier analysis works best on a set of frequencies that are repeated (over sampled) and have an average value that is close to zero, i.e. a sinusoidal waveform or a summation of sinusoidal waveforms. A windowing function is used to alleviate some of these concerns. In this Maple worksheet a slowly changing background (x value of the integration) is analysed in order to find the forcing frequency used in the equation of motion. For this work the wavelet library of Gelinas is used to get an array of Daubechies coefficient values to use in the analysis of the data.

Reading the share file directly is the only way I could get the Gelinas library daub to load correctly and it is done using the read statement and in back quotes having the file name preceded with the correct drive/directory/share/numerics/daub/daub.mpl then a back quote and a full colon. For example,
read `c:,d:e:,whichever is correct/yourmapledir/share/numerics/daub/daub.mpl`:

This would usually be read `c:/maple/share/numerics/daub/daub.mpl`:

> with (share); with (daub);

See ?share and ?share,contents for information about the share library

Share Library:  daub

Author: Gelinas, Jacques.

Description:  Package which approximates the coefficients of Daubechies low pass filters and the Daubechies Minium Phase filter

Next is the selection of the number of Daubechies coefficients to be used. For this work 8 coefficients are used by setting p=4 (so there will be 2*p-1=7 coefficients plus the zeroth coefficient for a total of 8 (c0..c7)) and then calling the daubc routine of Gelinas' daub share library package.

> p := 4; c:=array( 0..2*p-1, daubc(p) );

Now control constants, to be used in the calculation of the wavelet transform, are set.

Ntotal is the total number of data points from the integration routine used above. In this worksheet the data subset size has been choosen to be 512 data points and that is what N is set to while the number of subsets needed in Ntotal is calculated and set equal to ll. Then kk and mm are set initially to 1 and 512, respectively and each are incremented by 512 for each pass through the do loop for the data subsets.

> Ntotal:=nops(morsedverk);

> N:=512;ll:=Ntotal/N;kk:=1;mm:=512;

M is the largest value of the Daubechies coefficients, and NN is the number of passes through the detail, and new signal (degraded signal) do loop. M is used in the do loop to determine which parts of the summation to keep. Only data points that would fall into the range of the Daubechies coefficients (c0 to c7 in this work) are multiplied times the correct coefficient and summed to produce the analysed data. The code is currently corrected for wrap around of the Daubechies coefficients and the signal (the first few data points and last few data points) and has now been corrected to set the wrap for a M value up to and including 19 which is set when p=10 in the code listed above that calls daubc(p).

NN is determine by dividing the natural log of the number of data points in the data subset by the natural log of 2 and rounding the resultant value down and then subtracting 3 thus leaving 8 data points in the final file.

> M:=2*p-1;NN:=round(eval(ln(N)/ln(2)))-3;m:=1;

The wavelet analysis produces a degraded signal, f(1..NN) and details, D(1,2,..NN) each with half the data points of the previous signal. The new signal, for example f1 is a degraded signal produced from f0, and f1 is used in the do loop to produce the next degraded signal (f2) and the next detail (D2). This process is continued until the point where there is only 8 data points ( ) left in the signal (f9). The values in brackets are specific to this work.

If the wavelet routine is run without real numbers in the original signal f0, it takes much more time to do the calculations than with numbers and if the data set is doubled in size the time increases by a factor of four.

The first few lines of code listed below are for reading data into the smaller subsets. The next lines of code do a check to see if the data points be analyzed are within the Daubechies coefficient range of values and corrects the first few data points and last few data points depending upon the value of M. Currently data points being anlyzed with Daubechies coefficients up to and including M=19 (daub20) are corrected. The code in not elegant but runs quickly. The latter lines of code perform the wavelet analysis.

> t0:=time():
for lk from 1 to ll do
f0:=morsedverk[kk..mm,2]:

N:=nops(f0):NN:=round(eval(ln(N)/ln(2)))-3:

for m from 1 to NN do
for i to N/2 do
if i=1 then
bma:=add(c[kl]*f.(m-1)[N-(kl-3)],kl=3..M)+0: bmb:=-c[0]*f.(m-1)[N]:
elif i=2 then bma:=add(c[kl]*f.(m-1)[N-(kl-5)],kl=5..M)+0:bmb:=0:
elif i=3 then bma:=add(c[kl]*f.(m-1)[N-(kl-7)],kl=7..M)+0:bmb:=0:
elif i=4 then bma:=add(c[kl]*f.(m-1)[N-(kl-9)],kl=9..M)+0:bmb:=0:
elif i=5 then bma:=add(c[kl]*f.(m-1)[N-(kl-11)],kl=11..M)+0:bmb:=0:
elif i=6 then bma:=add(c[kl]*f.(m-1)[N-(kl-13)],kl=13..M)+0:bmb:=0:
elif i=7 then bma:=add(c[kl]*f.(m-1)[N-(kl-15)],kl=15..M)+0:bmb:=0:
elif i=8 then bma:=add(c[kl]*f.(m-1)[N-(kl-17)],kl=17..M)+0:bmb:=0:
elif i=9 then bma:=add(c[kl]*f.(m-1)[N-(kl-19)],kl=19..M)+0:bmb:=0:
elif i=(N/2)-8 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-18],kl=19..M)+0:bma:=0:
elif i=(N/2)-7 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-16],kl=17..M)+0:bma:=0:
elif i=(N/2)-6 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-14],kl=15..M)+0:bma:=0:
elif i=(N/2)-5 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-12],kl=13..M)+0:bma:=0:
elif i=(N/2)-4 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-10],kl=11..M)+0:bma:=0:
elif i=(N/2)-3 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-8],kl=9..M)+0:bma:=0:elif i=(N/2)-2 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-6],kl=7..M)+0:bma:=0: elif i=(N/2)-1 then bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-4],kl=5..M)+0:bma:=0
elif i=(N/2) then
bmb:=add((-1)^(kl+1)*c[kl]*f.(m-1)[kl-2],kl=3..M):bma:=c[0]*f.(m-1)[1]:
else
bma:=0:bmb:=0:fi:
for j to N do
if (2*i-j+1)<0 then bma:=bma:
elif (2*i-j+1)>M then bma:=bma:
else bma:=bma+c[2*i-j+1]*f.(m-1)[j] fi:
if (j+2-2*i)<0 then bmb:=bmb:
elif (j+2-2*i)>M then bmb:=bmb:
else bmb:=bmb+(-1)^(j+1)*c[j+2-2*i]*f.(m-1)[j] fi:
od:
f.m[i]:=bma/2:detail.m[N/2*(lk-1)+i]:=bmb/2:
od:N:=N/2:
od:kk:=kk+512:mm:=mm+512:
od:
(time()-t0)*seconds;
#Ntotal=4096 N=512 M=7, calculation = 50.42 sec
# compared to the total set of 4096 taking 432.66 seconds

When plotting "Details" keep in mind that the number of data points decreases by a factor of two for each increasing detail number. f0 started with 4096 therefore Detail_1 has 2048, Detail_2 has 1024 and Detail_3 (plotted below) only has 512 data points. Therefore the maxiumu value for j must be increased and the muliplier for time data point must be decreased. For example, in order to plot Detail_1 the sequence of data points is produced by using [seq([j*2*204.75/4096, detail1[j]], j=1..2048)] instead of the set of instructions below. The reader may want to plot the intervening signals as well which are in files f.(0..NN). In this work they are in f0, f1, f2, .. f9.

> st:=time():
### WARNING: the definition of the type `symbol` has changed'; see help page for details
p1:=plots[pointplot]([seq([j*8*204.75/4096, detail3[j]], j=1..512)],color=magenta, style=point,symbol=circle):
### WARNING: the definition of the type `symbol` has changed'; see help page for details
p2:=plots[pointplot]([seq([j*8*204.75/4096, detail3[j]], j=1..512)],color=magenta, style=line,symbol=circle):
display(p1,p2,view=[10..40.0,-0.00001..0.00001],axes=BOXED,labels=[time,`D `],title=`Detail_3 vs Time (512 subset)`);
time() - st;

The reader will notice that there is a spike in the waveform presented Detail3 compared to the same figure in the previous work and this difference is due to the 512 data point data subsets being used and the slight wavelet analysis weakness in begining and ending data points of a data set[6].

Conclusion

The information pulled from the slow moving signal of x, for the time between 10 and 40 seconds, has a sinusoidal shape with a period of approximately as expected for the forcing frequency being 1 rad/sec. The reader may want to decrease the forcing amplitude and observe that the amplitude of the waveform in the detail decreases in value.

There is a eight times decrease in the amount of time need to calculate the wavelet analysis using the 512 data point subsets and an overall decrease of five times for the whole worksheet which may be acceptable if the errors introduced into the data are easily recognizable as above. The reader may want to change the subset size and is referred to the parts of the code that are incremented by 512 for those changes. Large time savings can be accomplished but errors are introduced at the leading and trailing edges of each data subset. Caution is advised.

The reader could recalculate the results using different integration methods and see how this effects the results. For example, try setting method=mgear[adamspc] and removing the tolerance=0.0001 information from the routine and running the integration again and the wavelet analysis routine again.

The reader may also wish to view the other details that have been calculated. The details are from 1 to M (9 for the example above). The reader is cautioned that if the Digits setting is increased from 12 to 15 the integration time can increase by a factor of 10. If a data set is read in, to be evaluated by the wavelet routine, the reader is again cautioned in that the doubling of the data set size can increase the wavelet calculation time by a factor of four.

References:

[1] D. N. S. Permann, "Nonlinear Effects in Chemical Dynamics and Chemical Kinetics: Chaos in Physical Chemistry", (Ph. D., University of Ottawa, 1993)

[2] P. M. Morse, Phys. Rev., 34, 57 (1929)

[3] C. N. Banwell, "Fundamentals of Spectroscopy" (3rd edition, McGraw-Hill, London, 1983)

[4] I. Daubechies, "Ten Lectures on Wavelets" (Society for Industrial and Applied Mathematics, Philadelphia, 1992)

[5] J. Gelinas, Maple share library numerics/daub, included in MAPLE V 4,5

[6] R. Timothy Edwards, at John Hopkins University "Discrete Wavelet Transforms: Theory and Implementation" and the paper is in Tim's research section under "wavelet chip project" that he did at Stanford. The web page is currently at http://bach.ece.jhu.edu/~tim/research/wavelets/wavelets.ps.

Disclaimer: While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.