Commutators and Maxima

Students who can not afford Mathematica and/or Maple should consider Maxima. This "free" resource will here be used to illustrate commutators in quantum chemistry.


We start with the definition of momentum:

which we use to form the commutator of the Hamiltonian with a ladder operator.

Consider the harmonic oscillator, whose Hamiltonian is:

and the ladder operator (up/+):

where: . We incorporate this to eliminate "k" from the Hamiltonian, i.e., , so



The commutator of these two is  but we wish to use Maxima to help us understand what this means. Thus, intially using Maxima along with a function f(x), we enter:

k : mu*omega^2$ H : (-%i * hbar * diff(-%i * hbar * diff(f(x),x),x))/(2*mu) + (((mu*omega^2)/2)*x^2)*f(x)$ aplus : -%i * hbar * diff(f(x),x)+(%i *mu*omega*x)*f(x)$ Haplus : (-%i * hbar * diff(-%i * hbar * diff( f(x) ,x),x))/(2*mu) + (((k)/2)*x^2)* f(x) $ Haplus : (-%i * hbar * diff(-%i * hbar * diff( (aplus) ,x),x))/(2*mu) + (((k)/2)*x^2)* (aplus) $ aplusH:-%i * hbar * diff( (H) ,x)+(%i *mu*omega*x)* (H) $ comm : expand(Haplus-aplusH); comm2 : expand(subst(diff(f(x),x)=-p_op/(%i*hbar),comm)); tex(comm2);

which, when sent to Maxima becomes

which is the correct result.

Since this article deals with commutators, which are not algebraic in the traditional sense of computational calculus, we have to use workarounds to implement what we're interested in. First, we write "H" as the first part of "Haplus" in a form isolating "f(x)" and then we (while editing) substitute "(aplus)" for "f(x)" twice. Then we write "aplus" as the first part of "aplusH" in a form isolating "f(x)" and the nwe (while editing) substitute "(H)" for "f(x)" twice.

With these two object obtained, "Haplus" and 'aplusH" we form the commutator by subtraction and clean up the resulting expression.

Once we're sure that all is correct, we change the (normally used line terminating) semicolons into dollar signs so that we can obtain a clean output of the last line.
This output needs editing and interpretation. Specifically, we need to transform   (really ) back into a momentum operator.

The result is  which is known to be correct.

We can alter the code to:
k : mu*omega^2; H : (-%i * hbar * diff(-%i * hbar * diff(f(x),x),x))/(2*mu) + (((mu*omega^2)/2)*x^2)*f(x)$ aplus : -%i * hbar * diff(f(x),x)+(%i *mu*omega*x)*f(x); Haplus : (-%i * hbar * diff(-%i * hbar * diff( aplus ,x),x))/(2*mu) + (((k)/2)*x^2)* aplus ; aplusH:-%i * hbar * diff( (H) ,x)+(%i *mu*omega*x)* (H) ; comm : expand(Haplus-aplusH); comm2 : expand(subst(diff(f(x),x)=-p_op/(%i*hbar),comm)); tex(comm2);
which works as well.

Here's my last attempt:

k : mu*omega^2;
H : (-%i * hbar * diff(-%i * hbar * diff(f(x),x),x))/(2*mu) + (((mu*omega^2)/2)*x^2)*f(x);
aplus : -%i * hbar * diff(f(x),x)+(%i *mu*omega*x)*f(x);
Haplus : expand(subst(f(x)=aplus,H));
aplusH:expand(subst(f(x)=H,aplus));

comm : ratexpand(Haplus-aplusH);
comm:comm,diff,expand,factor;/*this line, from Daniel Volinsk, made this code work*/ 

comm2 : expand(subst(diff(f(x),x)=-p_op/(%i*hbar),comm));