(* Mathematica program for Problems 9.1 to 9.6

  date : 28/10/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt
*)

(* Definitions *)
A={{0,1,0,0,0,0},
   {0,0,Sqrt[2],0,0,0},
   {0,0,0,Sqrt[3],0,0},
   {0,0,0,0,Sqrt[4],0},
   {0,0,0,0,0,Sqrt[5]},
   {0,0,0,0,0,0}};

AP=Transpose[A];
H=hbar w DiagonalMatrix[{1/2,3/2,5/2,7/2,9/2,11/2}] ;
x= Sqrt[hbar/(2 m w)] ( A + AP) ;
p= -I Sqrt[hbar  m w /2] ( A - AP) ;
u0={1,0,0,0,0,0} ;

(* Problem 9.1 *)
AP2= AP . AP;
AP3=AP2 . AP;
AP4=AP3 . AP;
Print["(A^+)^2 = ",AP2]
Print["(A^+)^3 = ",AP3]
Print["(A^+)^4 = ",AP4]

(* Problem 9.2 *)
Print["x = ",x]
Print["x^2 = ",x . x]

(* Problem 9.3 *)
Print["p = ",p]
Print["p^2 = ",p . p]

(* Problems 9.4 e 9.5 *)
u1=1/Sqrt[1] AP . u0;
u2=1/Sqrt[2!] AP2 . u0;
u3=1/Sqrt[3!] AP3 . u0;
u4=1/Sqrt[4!] AP4 . u0;
Print["u0 = ",u0]
Print["u1 = ",u1]
Print["u2 = ",u2]
Print["u3 = ",u3]
Print["u4 = ",u4]

(* Problem 9.6 *)
psi=1/Sqrt[6] {1,2,1,0,0,0};
Hmed = psi . H . psi;
xmed = Simplify[ psi . x . psi];
x2med = Simplify[psi . x . x . psi];
pmed = Simplify[ psi . p . psi];
p2med = Simplify[ psi . p . p . psi];
dx2 = Simplify[x2med-xmed^2];
dp2 = Simplify[p2med-pmed^2];
heisenberg=Simplify[Sqrt[dx2 dp2],Assumptions->hbar>0];
heisenbergN=N[heisenberg];
Print["<H> = ",Hmed]
Print["<x> = ",xmed]
Print["<x^2> = ",x2med]
Print["<p> = ",pmed]
Print["<p^2> = ",p2med]
Print["(Δx)^2 = ",dx2]
Print["(Δp)^2 = ",dp2]
Print["Δx Δp = ",heisenberg," = ",
      heisenbergN," > ", 1/2 hbar]