(* Mathematica program for Problem 3.9

  date : 26/10/2009
  author: Jorge C. Romão
  email : jorge.romao@tecnica.ulisboa.pt
*)

(* Conventions

 We use as unit of length y=Sqrt[(m w)/hbar] x. We also take w=1, so time 
 is defined in units of 1/w.

*)

(* Normalized eigenfunctions *)

u0=Function[y, 1/Pi^(1/4) E^(-y^2/2)] ;
u1=Function[y, 1/Pi^(1/4) 1/Sqrt[2] 2 y E^(-y^2/2)] ;

(* Wave function *)

Psi=Function[{y,t},1/Sqrt[2] (u0[y] E^(-I 1/2 t) + u1[y] E^(- I 3/2 t))];

(* Fast definiton of the density of probability *)

Prob=Function[{y,t}, 1/2(u0[y] + u1[y] Cos[t])^2 + 1/2 u1[y]^2 Sin[t]^2];

(* Verify Normalization *)

Testu0 = Integrate[u0[y]^2,{y,-Infinity,Infinity}];
Testu1 = Integrate[u1[y]^2,{y,-Infinity,Infinity}];
TestPsi = Integrate[Conjugate[Psi[y,t]] Psi[y,t],{y,-Infinity,Infinity},
		    Assumptions->t>0];

Print["Test u0=",Testu0]
Print["Test u1=",Testu1]
Print["Test Psi=",TestPsi]

(* Probability Density *)

Animate[Plot[Abs[Psi[y,t]]^2, {y,-3,3},PlotRange->{0,0.8}, 
  AxesLabel -> {"x","|ψ|^2"}],{t,0,2 Pi,Pi/4}, 
  AnimationRunning -> False]

(* Average <x>, <x^2> as function of time *)

xaverage=Integrate[y Prob[y,t],{y,-Infinity,Infinity},Assumptions->t>0];
x2average=Integrate[y^2 Prob[y,t],{y,-Infinity,Infinity},Assumptions->t>0];

Plot[xaverage,{t,0,2 Pi},AxesLabel -> {"t","<x>"}]

(* Average of <p>, <p^2> as function of time *)

paverage=-I Integrate[Conjugate[Psi[y, t]] D[Psi[y, t], y], 
                      {y, -Infinity, Infinity},Assumptions -> t > 0];
p2average=- Integrate[Conjugate[Psi[y, t]] D[Psi[y, t], {y,2}], 
		      {y, -Infinity, Infinity},Assumptions -> t > 0];

(* Heisenberg Relation *)
(* In all cases DxDp > 1/2 (hbar=1) *)

DxDp=Sqrt[x2average-xaverage^2] Sqrt[p2average-paverage^2]

Plot[DxDp,{t,0,Pi},AxesLabel -> {"t","ΔxΔp"},
     PlotRange->{0.5,0.8}]

(* Virial Theorem *)
 (*  The virial theorem states that for the harmonic oscillator

     < T > = < V > => < E > = < T > + < V > = 2 < T > = 2 < V >

with, in our conventions,

T= p^2/2 and V=x^2/2 
 
  and ( in units of hbar w)

< E > = 1/2 1/2 + 1/2 3/2 = 1

 *)

Print["< T > = ",1/2 p2average]
Print["< V > = ",1/2 x2average]
Print["< E > = < T > + < V >= ",1/2 p2average+1/2 x2average]