> # This spread sheet computes the answers to the example problem
> # of Homework 1.
> 
> #  Please do not use the Maple worksheet QC6_Lib.mws for this example
> #  This worksheet works without it.
> 
> with(linalg):
> 
Warning, the protected names norm and trace have been redefined and
unprotected

> # We wish to observe the state psi with respect to the observable A
> 
> psi:=evalm( (1/sqrt(3)) * matrix(4,1,[1,I,0,-1]) );
> A:=matrix(4,4, [0, 0,1,-I,
>                 0, 0,I,-1,
>                 1,-I,0, 0,
>                 I,-1,0, 0  ]  );

                               [ 1/3 sqrt(3) ]
                               [             ]
                               [1/3 I sqrt(3)]
                        psi := [             ]
                               [      0      ]
                               [             ]
                               [- 1/3 sqrt(3)]


                           [0    0     1    -I]
                           [                  ]
                           [0    0     I    -1]
                      A := [                  ]
                           [1    -I    0    0 ]
                           [                  ]
                           [I    -1    0    0 ]

> # We now compute the two eingenvals of A along with their eigenvectors
> 
> List:=[eigenvects(A)];
> a1:=op(1, op(1,List) );
> v1:=op(1, op(3, op(1,List) )  );
> v2:=op(2, op(3, op(1,List) )  );
> 
> a2:=op(1, op(2,List) );
> v3:=op(1, op(3, op(2,List) )  );
> v4:=op(2, op(3, op(2,List) )  );

  List := [[sqrt(2), 2, {[0, -sqrt(2), I, 1], [1, I, sqrt(2), 0]}],

        [-sqrt(2), 2, {[0, sqrt(2), I, 1], [1, I, -sqrt(2), 0]}]]


                            a1 := sqrt(2)


                       v1 := [1, I, sqrt(2), 0]


                      v2 := [0, -sqrt(2), I, 1]


                            a2 := -sqrt(2)


                       v3 := [0, sqrt(2), I, 1]


                      v4 := [1, I, -sqrt(2), 0]

> # We now you the GramSchmidt to create orthogonal bases for the two
> eigenspaces V+ and V-.
> 
> VectorSet1:=GramSchmidt({v1,v2});
> ww1:=op(1,VectorSet1);   ww2:=op(2,VectorSet1);
> VectorSet2:=GramSchmidt({v3,v4});
> ww3:=op(1,VectorSet2);   ww4:=op(2,VectorSet2);

  VectorSet1 :=

        {[1, I, sqrt(2), 0], [-1/2 I sqrt(2), - 1/2 sqrt(2), 0, 1]}


                      ww1 := [1, I, sqrt(2), 0]


             ww2 := [-1/2 I sqrt(2), - 1/2 sqrt(2), 0, 1]


  VectorSet2 :=

        {[0, sqrt(2), I, 1], [1, 0, - 1/2 sqrt(2), -1/2 I sqrt(2)]}


                      ww3 := [0, sqrt(2), I, 1]


             ww4 := [1, 0, - 1/2 sqrt(2), -1/2 I sqrt(2)]

> # We now normalize the eigenkets
> 
> for i from 1 to 4 do
>    w||i:=convert(evalm( (1/norm(ww||i, 2)) * ww||i ), matrix):
> od;

                               [    1/2    ]
                               [           ]
                               [   1/2 I   ]
                         w1 := [           ]
                               [1/2 sqrt(2)]
                               [           ]
                               [     0     ]


                               [  -1/2 I   ]
                               [           ]
                               [   -1/2    ]
                         w2 := [           ]
                               [     0     ]
                               [           ]
                               [1/2 sqrt(2)]


                               [     0     ]
                               [           ]
                               [1/2 sqrt(2)]
                         w3 := [           ]
                               [   1/2 I   ]
                               [           ]
                               [    1/2    ]


                               [1/2 sqrt(2)]
                               [           ]
                               [     0     ]
                         w4 := [           ]
                               [   -1/2    ]
                               [           ]
                               [  -1/2 I   ]

> # Thus, the eigenvectors for the eigenvalue a1 are w1 and w2.
> # They form an orthonormal basis of the eigenspace  V+
> # Also, the eigenvectors for the eigenvalue a2 are w3 and w4.
> # They form an orthonormal basis of the eigenspace  V- 
> # Moreover, all w1 & w2 are each mutually orthogonal to w3 & w4
> # because eigenvectors corresponding to distinct eigenvectors
> # must be orthogonal.
> 
> # We can now verify that the spectral decomposition theorem is true,
> # which of course it is.
> 
> Proj1:=evalm( w1 &* htranspose(evalm(w1)) );
> Proj2:=evalm( w2 &* htranspose(evalm(w2)) );
> Proj3:=evalm( w3 &* htranspose(evalm(w3)) );
> Proj4:=evalm( w4 &* htranspose(evalm(w4)) );
> 
> 

             [    1/4            -1/4 I         1/4 sqrt(2)     0]
             [                                                   ]
             [   1/4 I            1/4          1/4 I sqrt(2)    0]
    Proj1 := [                                                   ]
             [1/4 sqrt(2)    -1/4 I sqrt(2)         1/2         0]
             [                                                   ]
             [     0               0                 0          0]


            [     1/4             1/4 I        0    -1/4 I sqrt(2)]
            [                                                     ]
            [   -1/4 I             1/4         0    - 1/4 sqrt(2) ]
   Proj2 := [                                                     ]
            [      0                0          0          0       ]
            [                                                     ]
            [1/4 I sqrt(2)    - 1/4 sqrt(2)    0         1/2      ]


             [0          0                0                0     ]
             [                                                   ]
             [0         1/2         -1/4 I sqrt(2)    1/4 sqrt(2)]
    Proj3 := [                                                   ]
             [0    1/4 I sqrt(2)         1/4             1/4 I   ]
             [                                                   ]
             [0     1/4 sqrt(2)         -1/4 I            1/4    ]


            [     1/2          0    - 1/4 sqrt(2)    1/4 I sqrt(2)]
            [                                                     ]
            [      0           0          0                0      ]
   Proj4 := [                                                     ]
            [- 1/4 sqrt(2)     0         1/4            -1/4 I    ]
            [                                                     ]
            [-1/4 I sqrt(2)    0        1/4 I             1/4     ]

> # So the Projection operators for the eigenspace V+ and V- are
> respectively
> ProjPlus:= evalm(Proj1 + Proj2);
> ProjMinus:=evalm(Proj3 + Proj4);

                    [1/2 ,  0 , 1/4 sqrt(2) , -1/4 I sqrt(2)]
                    [                                       ]
                    [0 , 1/2 , 1/4 I sqrt(2) , - 1/4 sqrt(2)]
        ProjPlus := [                                       ]
                    [1/4 sqrt(2) ,  -1/4 I sqrt(2) , 1/2 , 0]
                    [                                       ]
                    [1/4 I sqrt(2) , - 1/4 sqrt(2) , 0 , 1/2]


                     [1/2 , 0 , - 1/4 sqrt(2) , 1/4 I sqrt(2)]
                     [                                       ]
                     [0 ,  1/2 , -1/4 I sqrt(2) , 1/4 sqrt(2)]
        ProjMinus := [                                       ]
                     [- 1/4 sqrt(2) , 1/4 I sqrt(2) , 1/2 , 0]
                     [                                       ]
                     [-1/4 I sqrt(2) ,  1/4 sqrt(2) , 0 , 1/2]

> # Since by the spectral decomposition theorem, A = a1* ProjPlus +
> a2*ProjMinus,
> # we should get A with the following computation
> 
> evalm(  evalm(a1*ProjPlus) + evalm(a2*ProjMinus)  );

                         [0    0     1    -I]
                         [                  ]
                         [0    0     I    -1]
                         [                  ]
                         [1    -I    0    0 ]
                         [                  ]
                         [I    -1    0    0 ]

> # We now compute the coefficients of the expansion of psi in
> # the eigenbasis w1, w2, w3, w4
> 
> alpha||1:=evalm(htranspose(w1) &* psi)[1,1];
> alpha||2:=evalm(htranspose(w2) &* psi)[1,1];
> alpha||3:=evalm(htranspose(w3) &* psi)[1,1];
> alpha||4:=evalm(htranspose(w4) &* psi)[1,1]; 

                        alpha1 := 1/3 sqrt(3)


                   alpha2 := - 1/6 sqrt(2) sqrt(3)


            alpha3 := 1/6 I sqrt(2) sqrt(3) - 1/6 sqrt(3)


            alpha4 := 1/6 sqrt(2) sqrt(3) - 1/6 I sqrt(3)

> # We now verify that these coefficients are correct.
> 
> map(simplify, evalm(evalm(alpha1*w1 + alpha2*w2) + evalm(alpha3*w3 +
> alpha4*w4)));

                           [ 1/3 sqrt(3) ]
                           [             ]
                           [1/3 I sqrt(3)]
                           [             ]
                           [      0      ]
                           [             ]
                           [- 1/3 sqrt(3)]

> # Thus, the probability that eigenvalue a1 is read
> # as a result of the measurement is:
> 
> ProbPlus:= abs(alpha1)^2 + abs(alpha2)^2;
> ProbMinus:=abs(alpha3)^2 + abs(alpha4)^2;

                           ProbPlus := 1/2


                           ProbMinus := 1/2

> # The respective resulting states are:
> 
> StatePlus:= map(simplify, 
>           evalm( (1/sqrt(abs(alpha1)^2 + abs(alpha2)^2)) *
> evalm(alpha1*w1 + alpha2*w2) )
>                        );
> StateMinus:=map(simplify, 
>           evalm( (1/sqrt(abs(alpha3)^2 + abs(alpha4)^2)) *
> evalm(alpha3*w3 + alpha4*w4) )
>                        );
> 

                      [1/6 sqrt(2) sqrt(3) + 1/6 I sqrt(3) ]
                      [                                    ]
                      [1/12 sqrt(2) sqrt(3) (2 I + sqrt(2))]
         StatePlus := [                                    ]
                      [            1/3 sqrt(3)             ]
                      [                                    ]
                      [       - 1/6 sqrt(2) sqrt(3)        ]


                       [1/6 sqrt(2) sqrt(3) - 1/6 I sqrt(3)]
                       [                                   ]
                       [1/6 I sqrt(2) sqrt(3) - 1/6 sqrt(3)]
         StateMinus := [                                   ]
                       [           - 1/3 sqrt(3)           ]
                       [                                   ]
                       [       - 1/6 sqrt(2) sqrt(3)       ]

>