A MAPLE Approach to Wavelet Analysis
of the Morse Oscillator

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 for calculating a wavelet analysis of the Morse oscillator motion in an anharmonic region, i.e close to dissociation of a single bond. This Maple V 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].

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.

Section I: 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=POINT,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]);

Section II: 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.

Section III: Integration 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 N 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.

> N:=nops(morsedverk);

> f0:=morsedverk[1..N,2]:

Section IV: 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. 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.

NN is determine by dividing the natural log of the total number of data points in the data set by the natural log of 2 and rounding the resultant value down and then subtracting 3.

> 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.

> t0:=time():
for m from 1 to NN do
for i to N/2 do
bma:=0:bmb:=0:
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[i]:=bmb/2:
od:N:=N/2:
od:
(time()-t0)*seconds;
#N=512 M=7, calculation = 5.83 seconds
#N=1024 M=7, calculation = 24.33 seconds
#N=2048 M=7, calculation = 100.94 seconds
#N=4096 M=7, calculation = 432.66 seconds
#N=8192 M=7, calculation =1899.33 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`);
time() - st;

>

Section V: Conclusion

The information pulled from the slow moving signal of x, for the time between 10 and 40 seconds, has a sinusoidal shape and 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.

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

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.

>