Application Center - Maplesoft

App Preview:

Diatomic anharmonic oscillator

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


The Diatomic Anharmonic Oscillator according to Matrix Mechanics

by M. B. Monagan and J. F. Ogilvie, Simon Fraser University, Centre for Experimental and Constructive Mathematics, Vancouver, Canada

NOTE: The calculation undertaken here employs perturbation theory of multiple orders; the extent of the calculation, according to the order, is adjustable by a user within limitations of machine capacity and acceptable duration of calculation.


Matrix mechanics is the original form of quantum mechanics that Werner Heisenberg initiated in 1925 and elaborated in association with Max Born and Pascual Jordan. By this means one can solve physical problems, such as that of an anharmonic oscillator that usefully models vibrational properties of a diatomic molecule as we demonstrate here. The full calculation is practicable with Maple statements and procedures that follow, and an augmented description and explanation appears in a paper to be published during 1999. Alternative procedures in Maple to generate energies of vibration-rotational states of a diatomic molecule, based on BKW quasi-classical theory and hypervirial perturbation theory in wave mechanics, are available in separate worksheets.

The calculation undertaken here employs perturbation theory of multiple orders; the extent of the calculation, according to the order, is adjustable by a user within limitations of machine capacity and acceptable duration of calculation.

> restart:


To calculate terms of perturbative order K in series for matrices E and U , we proceed as follows. First we construct a matrix [Maple Math] to a given rank. Next we create a matrix equation for the coefficient of [Maple Math] by expanding an equation [Maple Math] symbolically as a power series in [Maple Math] . We insert all known matrices truncated to a specified rank and solve this equation for [Maple Math] and off-diagonal entries of [Maple Math] . Next, to solve for diagonal entries of [Maple Math] , we make use of an equation [Maple Math] . We expand it as a power series in [Maple Math] symbolically and select the coefficient of [Maple Math] . We insert all known quantities calculated up to this point and solve for diagonal entries of [Maple Math] . After solving for [Maple Math] and [Maple Math] , one can proceed to calculate terms of order [Maple Math] .

Setting up Maple

The following statement tells Maple that [Maple Math] is a scalar quantity.

> constants := constants, lambda;

The following Maple procedure is used to simplify matrices of formulae; it factors out any common factor present in each entry of a matrix.

> factormatrix := proc(A)
local multigcd, B, L, N, D, G;
multigcd := proc(a)
local i,g;
g := a[1];
for i from 2 to nops(a) while g <> 1 do
g := frontend(gcd,[g,a[i]]);
B := map(factor,A);
L := {op(map(op,convert(B,listlist)))} minus {0};
if nops(L) = 0 then RETURN( eval(B) ) fi;
N := multigcd(map(numer,L));
if N <> 0 then N := N / icontent(N) fi;
if N = 0 then RETURN(eval(B)) fi;
D := multigcd(map(denom,L));
if D <> 0 then D := D / icontent(D) fi;
if D = 0 then RETURN(eval(B)) fi;
G := N/D;
if G = 1 or G = -1 then RETURN(eval(B)) fi;
G * map( unapply( x/G, x ), B )

With the following command we enable to extract a submatrix from a matrix.

> with(linalg,submatrix);

Construction of matrix equations

In construction of these two matrix equations, in the following procedure we assign coefficients of [Maple Math] to variables [Maple Math] and [Maple Math] to a specified order [Maple Math] . This part of the computation is truly ``symbolic'' in that equations appear in terms of unknown matrices. We subsequently fix the rank of these matrices, fill in known entries and solve for unknown entries.

> makeEquations := proc(D)
local i,H,U,E,G,Ut,UC;
H := add(H.i*lambda^i,i=0..D);
U := 1+add(U.i*lambda^i,i=1..D);
E := H0 + add(E.i*lambda^i,i=1..D);
G := expand(U &* H - E &* U);
G := algsubs(lambda^(D+1)=0, G);
Ut := 1 + add( transpose(U.i)*lambda^i, i=1..D );
UC := expand( U &* Ut - 1 );
UC := algsubs(lambda^(D+1)=0, UC);
for i from 1 to D do
G.i := coeff(G,lambda,i);
UC.i := coeff(UC,lambda,i);

> makeEquations(3);

For instance here are linear coefficients as determined previously,

> G1; UC1;

and here are quadratic coefficients.

> G2; UC2;

Matrix [Maple Math] is diagonal; this condition makes the algorithm work so simply. In the linear coefficient above, a quantity [Maple Math] has diagonal entries that are zero; diagonal entries of [Maple Math] are thus equal to those in [Maple Math] . In the quadratic coefficient, the unknown quantities are [Maple Math] and [Maple Math] ; diagonal entries of the first expression are again zero. Diagonal entries of [Maple Math] are hence given by diagonal entries of a matrix [Maple Math] . The same observation holds for terms of greater order that have a general form [Maple Math] . Diagonal entries of [Maple Math] are simply determined and off-diagonal entries of [Maple Math] give rise to simple linear equations. Here is the cubic coefficient.

> G3; UC3;

Construction of matrices H

In the following procedure we construct matrix [Maple Math] for K>0 , bearing in mind that all matrices are formally of infinite rank. A parameter N specifies the rank of a truncated matrix. To obtain [Maple Math] correct to N rows and N columns we compute [Maple Math] ; we begin with N+K rows and columns for X .

> makeH := proc(K,N)
local m,X,H;
X := matrix(N+K,N+K,0);
for m from 0 to N+K-2 do
X[m+1,m+2] := A*sqrt(m+1);
X[m+2,m+1] := X[m+1,m+2];
H := evalm( k.(K+2) * X^(K+2) );

Here is matrix [Maple Math] truncated to rank 6.

> H1 := makeH(1,6): factormatrix(H1);

Formule for bands we obtain by interpolation:

[Maple Math] for [Maple Math]


[Maple Math] for [Maple Math] .

Likewise, here is matrix [Maple Math] truncated to rank 6.

> H2 := makeH(2,6): factormatrix(H2);

Diagonal entries are given by

[Maple Math] for [Maple Math]

and the next band is given by

[Maple Math] for [Maple Math]


[Maple Math] for [Maple Math] .

In general, matrices H are banded and symmetric with K+3 non-zero bands and with bands with zero entries between non-zero bands. Formulae for the bands can be determined by inspection and interpolation.

Solving for [Maple Math] and [Maple Math] .

The following programme in the form of a Maple procedure calculates matrices E and U to order K (entered as D ) truncated to rank N .

> AnHarmonic := proc(D,N)
global U0,H0,E0,Order;
local i,j,m,x,K,T,Z,U,E,H,eqns,vars,sols;
U0 := array(1..N,1..N,identity);
H0 := array(1..N,1..N,diagonal);
for m from 0 to N-1 do H0[m+1,m+1] := (m+1/2)*omega od:
E0 := H0;
for K from 1 to D do
print(`Solving for the coefficient of`, lambda^K);
H.K := makeH(K,N);
E.K := array(1..N,1..N,diagonal);
U.K := array(1..N,1..N);
U := evaln(U.K); E := evaln(E.K); H := evaln(H.K);
printf(`matrix %s\n`,H); print(factormatrix(H.K));
printf(`matrix equation to be solved\n`); print(G.K);
Z := evalm(G.K);
if K=1 then printf(`expanded matrix equation to be solved\n`); print(Z) fi;
eqns := normal(convert(Z,set));
vars := {seq(seq(U[i,j],j=1..N),i=1..N), seq(E[i,i],i=1..N)} minus {seq(U[i,i],i=1..N)};
sols := solve(eqns,vars); assign(sols);
printf(`matrix %s\n`, E); print(factormatrix(E));
T := evalm( UC.K );
printf(`condition for the next U matrix`); print(UC.K = 0);
eqns := normal(convert(T,set));
vars := {seq(U[i,i],i=1..N)};
sols := solve(eqns,vars); assign(sols);
printf(`matrix %s\n`,U); print(factormatrix(U));

We calculate and display [Maple Math] then [Maple Math] ; we display also the symbolic matrix equations. Like [Maple Math] , all matrices E of odd order are 0. Before executing this procedure it is advisable to activate the following solver that greatly accelerates these calculations.

Fast solve

> unprotect(solve):

> solve := proc(eqns,unks)
if type(eqns[1],`=`) then RETURN( solve( map(proc(x) lhs(x)-rhs(x) end,eqns) minus {0}, unks ) ) fi;
map(proc(e,u) local x; if e=0 then RETURN() fi; x := indets(e) intersect(u); if nops(x)=1 then x := x[1]; x = normal( -coeff(e,x,0)/coeff(e,x,1) ) else ERROR(e,u) fi end, eqns, unks ) end;

The following test ensures that this procedure operates correctly.

> solve( {3*x=2,4*y=a}, {x,y} );

> AnHarmonic(2,6);

To verify that all computations performed with programme AnHarmonic are correct, that is, that we have correctly diagonalized H to [Maple Math] , we test that equation [Maple Math] holds to [Maple Math] .

> eqn := (1 + add(lambda^i*U.i,i=1..2)) &* add(lambda^i*H.i,i=0..2) =
add(lambda^i*E.i,i=0..2) &* (1 + add(lambda^i*U.i,i=1..2));

> map(simplify@series, evalm(lhs(eqn)-rhs(eqn)), lambda, 3 );

We proceed to examine the contribution to energy of second order; it is the matrix

> map(expand,E2);

Because we truncated matrices H , not all entries in matrices E and U are correct: actually only the first three diagonal entries in [Maple Math] are correct. We seek to compare our results with those of Wu [Quantum Mechanics, World Scientific, 1985], who uses the following definition for H .

[Maple Math] .

We convert our results for [Maple Math] to conform to Wu's notation.

> W2 := subs(k3=k2/3!, k4=k3/4!, omega=omega[0]*h,

These diagonal entries are consistent with equation (II-140) in Wu's book . The correction to energy of the harmonic oscillator in Wu's equation we obtain by interpolation from three values in the above matrix. Using Maple, we obtain

> interp([0,1,2],[W2[1,1],W2[2,2],W2[3,3]],n);

Two alternative ways to write this expression are as a polynomial in A or in [Maple Math] . In Maple

> collect(%,A,factor);

> map(collect,series(%,n=-1/2),A);

Corrections to greater order are obtained by running programme AnHarmonic with matrices of sufficient rank. Five statements, as given below starting with "AnHarmonic(4,10):" for instance, suffice to generate corrections of a given order; the duration of calculation increases exponentially with order (according to this, or any alternative, approach).

We end this section with comments on computations. Firstly, with regard to construction of [Maple Math] for K>0 , when we first computed these matrices, we used [Maple Math] [Maple Math] and [Maple Math] [Maple Math] etc. This procedure yields an algorithm to compute [Maple Math] of which the running time is exponential in [Maple Math] . An algorithm based on computing the matrix power [Maple Math] is superior, and it is sufficiently quick that there is little point in using formulae for bands of [Maple Math] .

The linear system of equations that is constructed and solved is a diagonal system, but, because formulae involve several square roots of small integers, Maple runs much more slowly than if coefficients were simple integers: the reason is that Maple is careful to determine algebraic relationships between square roots, e.g. [Maple Math] . In our case, if one does not mind seeing [Maple Math] and [Maple Math] in the answer, one can speed calculations by a factor 10 or more if one solves equations ``blindly''; we devised for this purpose a solver presented in the paragraph above marked Fast solve.


We illustrate here how matrix mechanics can be effectively used to solve an important problem in quantum mechanics using perturbation theory. Because matrix [Maple Math] is diagonal for a diatomic anharmonic oscillator, the requisite matrix algebra is simple; a consequence is that computations can be carried out exactly to high order using a system such as Maple for symbolic algebra.

We include in the following appendix statements that generate formulae for terms in energies of fourth, sixth and eighth orders.

Appendix Correction terms of Fourth, Sixth and Eighth Orders

> AnHarmonic(4,10):

> factormatrix(submatrix(E4,1..4,1..4));

> interp([0,1,2,3],[seq(E4[i,i],i=1..4)],n):

> collect(%,[A,k3,k4,k5,k6,omega],factor);

> evaln( E4[n,n] ) = map(collect, taylor(%,n=-1/2), A, factor);

Succeeding calculations require progressively increasing duration and memory. The first argument in each call to AnHarmonic indicates the perturbative order of a calculation, and the second argument specifies the rank of matrices involved in the calculation; this rank must increase exponentially with order to achieve accurate results.

> AnHarmonic(6,18):

> interp([0,1,2,3,4],[seq(E6[i,i],i=1..5)],n):

> evaln( E6[n,n] ) = map(collect, taylor(%,n=-1/2), A, factor);

E8 eighth order

> AnHarmonic(8,28):

> interp([0,1,2,3,4,5],[seq(E8[i,i],i=1..6)],n):

> evaln( E8[n,n] ) = map(collect, taylor(%,n=-1/2), A, factor);