The package OscillatorPlot - 1

Bernd Thaller
Institute for Mathematics
University of Graz

Introduction

We want to develop a Mathematica package that represents any given function in the energy-representation of the harmonic oscillator.

The eigenstates (Subscript[ψ, n])(x) of a harmonic oscillator form an orthonormal basis in the Hilbert space of square-integrable functions. Hence any given (complex-valued) function f(x) can be expanded in the basis of oscillator eigenfunctions:

f(x) =(∑_(n = 0)^∞) 
                      (Subscript[c, n]) (Subscript[ψ, n])(
                    x), where (Subscript[c, n])=∫(Subscript[ψ, n])(x) f(
                  x) dx

|(Subscript[c, n])\( | \^2\) is the probability to find the state f(x) in the n-th oscillator state (Subscript[ψ, n])(x). The sequence c of expansion coefficients {(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),...} describes the function f in the oscillator basis. c is sometimes called the energy-representation of f.

We want to define a command which plots the energy representation of a function. We want to plot the expansion coefficient (Subscript[c, n]) as a vertical line on the x-axis, with the x-coordinate giving the energy (Subscript[E, n]) of the state (Subscript[ψ, n]) and the height giving the absolute value |(Subscript[c, n])|. We also plan to draw each line with a color corresponding to its complex phase.

Several steps are necessary to accomplish this task. We plan to design a Mathematica package providing a command for this type of plot. First of all, we have to find a way how to do this plot.

First of all, we will need to define the normalized oscillator eigenfunctions (whose form can be found in any book on quantum mechanics):

OscillatorFunction[n_, x_] := 
        HermiteH[n, x]*Exp[-x^2/2]/Sqrt[2^n*n!]/
            Pi^(1/4);

This function is an eigenfunction of the energy operator H

H f(x)=-(1/2) (((∂^2) f(x))/(∂(x^2)))+(1/2) (x^2) f(x)

and belongs to the eigenvalue

(Subscript[E, n]) = n + (1/2)

In order to be specific, assume that we want to plot the energy representation of

func[x_] := (1/Pi)^(1/4) E^(-1/2*(x - 2)^2 - I 2 x);

Designing the plot command

We first have to determine the expansion coefficients by a numerical integration. We decide to evaluate just the first few coefficients. This takes some time:

Do[coeff[n] = 
        NIntegrate[
          OscillatorFunction[n, y] func[y], {y, -Infinity, 
            Infinity}], {n, 0, 10}]

Now, the coefficients have been pre-evaluated, and now it only remains to plot them.

Table[coeff[n], {n, 0, 10}]

We see that the coefficients are complex numbers. We will define a colored vertical line at the x-coordinate given by the energy of the corresponding state of the harmonic oscillator.

coloredLine[n_] := 
{Thickness[0.01], 
        Hue[Arg[coeff[n]]/2/Pi], 
        Line[{{n  +  1/2, Abs[coeff[n]]}, {n  +  1/2, 0. }}]}

Show[Graphics[coloredLine[5]]]

But we don't want to plot just one line, we want to plot all coefficients:

lineTable = Graphics[Table[coloredLine[n], {n, 0, 10}]];

Show[lineTable]

Now we just have to add coordinate lines:

Show[lineTable, Frame -> True]

Perhaps its better to make some reasonable assumptions about the PlotRange: We take the maximal absolute value of the coefficients and add a little bit above and below the lines:

up = 
        Max[Table[Abs[coeff[n]], {n, 0, 10}]];
upper = up + up/10;
lower = -up/10;

Show[lineTable, Frame -> True, 
      PlotRange -> {{0, 11}, {lower, upper}}]

Combining commands into a procedure

This looks fairly nice. We can now proceed to make a useful Mathematica command. In principle, we could just combine the commands above in a single module: For the purpose of testing we only take 6 expansion coefficients

Remove[OscillatorPlot]

OscillatorPlot[f_] := 
  Module[{coeff, coloredLine, lineTable, up, upper, lower},
    Do[coeff[n] = NIntegrate[
       OscillatorFunction[n, x]*f, 
       {x, -Infinity, Infinity}], {n, 0, 5}]; 
    coloredLine[n_] := {Thickness[0.01], 
      Hue[Arg[coeff[n]]/2/Pi], 
      Line[{{n + 1/2, Abs[coeff[n]]}, 
        {n + 1/2, 0.}}]};
    lineTable = Graphics[Table[coloredLine[n], {n, 0, 5}]];
    up=Max[Table[Abs[coeff[n]],{n,0,5}]];
    upper=up+up/10;
    lower=-up/10; 
    Show[lineTable, Frame -> True, PlotRange->{{0,6},{lower,upper}}]]

OscillatorPlot[func[x]]

However, the command defined above is not without problems. Try, for example,

OscillatorPlot[func[y]]

So we need a more appropriate way to pass the function argument to the Module. Moreover, we want to be able to change the borders of the region for the numerical integration, and we want to be able to change the number of coefficients in the final image.

Remove[OscillatorPlot]

OscillatorPlot[f_, {x_Symbol, x1_:-Infinity, x2_:Infinity},
    n1_:0, n2_:10] := 
  Module[{coeff, coloredLine, lineTable, up, upper, 
    lower, func},
    func = Function[{x}, f]; 
    Do[coeff[n] = NIntegrate[OscillatorFunction[n, x]*
        func[x], {x, x1, x2}], {n, n1, n2}]; 
    coloredLine[n_] := {Thickness[0.01], 
      Hue[Arg[coeff[n]]/(2*Pi)], 
      Line[{{n + 1/2, Abs[coeff[n]]}, 
        {n + 1/2, 0.}}]};
    lineTable = 
     Graphics[Table[coloredLine[n], {n, n1, n2}]]; 
    up = Max[Table[Abs[coeff[n]], {n, n1, n2}]]; 
    upper = up + up/10; lower = -(up/10); 
    Show[lineTable, Frame -> True,
        PlotRange -> {{n1, n2 + 1}, {lower, upper}}]]

OscillatorPlot[func[y], {y, -6, 10}]

OscillatorPlot[func[x], {x}, 3, 8]

Putting symbols in a private context

We can now put all commands in a file. We place the commands OscillatorFunction and OscillatorPlot in a private context and export OscillatorPlot for the main session. The file looks as follows

OscillatorPlot::usage=""

Begin["Private`"]

OscillatorFunction[n_, x_] := (HermiteH[n, x]*Exp[-(x^2/2)])/
    (Sqrt[2^n*n!]*Pi^(1/4))

OscillatorPlot[f_,{x_Symbol,x1_:-Infinity,x2_:Infinity},
               n1_:0,n2_:10]:=
   Module[{coeff,coloredLine,lineTable,up,upper,lower,func},
      func=Function[{x},f];
      Do[coeff[n]=
         NIntegrate[OscillatorFunction[n,x]*func[x],{x,x1,x2}],{n,n1,n2}];
      coloredLine[n_]:={Thickness[0.01],Hue[Arg[coeff[n]]/(2*Pi)],
                     Line[{{n+1/2,Abs[coeff[n]]},{n+1/2,0.}}]};
      lineTable=Graphics[Table[coloredLine[n],{n,n1,n2}]];
      up=Max[Table[Abs[coeff[n]],{n,n1,n2}]];
      upper=up+up/10; lower=-(up/10);
      Show[lineTable, Frame->True, PlotRange->{{n1,n2+1}, {lower,upper}}]
   ]
   
End[]

We name this file "OscillatorPlot1.m" and save it as a pure text file.
If it is on the Mathematica search path, it can be read with the Get command.

<< OscillatorPlot1.m

OscillatorPlot[Exp[-x^2] Sin[x - 1], {x}] // Timing