Lab 1 - Vectors and Matrices

This notebook demonstrates the use of vectors and matrices in IPython. Note that the basis is not explicit in any of these operations. You must keep track of the basis yourself (using variable names, or notes etc).



In [1]:

    
from numpy import array, dot, outer, sqrt, matrix
from numpy.linalg import eig, eigvals
from matplotlib.pyplot import hist



In [2]:

    
%matplotlib inline



In [3]:

    
rv = array([1,2])  # a row vector
rv









    Out[3]:





array([1, 2])



In [4]:

    
cv = array([[3],[4]])  # a column vector
cv









    Out[4]:





array([[3],
       [4]])

Two kinds of vector products we'll see: inner product (dot product) and outer product

1) Use the function dot(vector1, vector2) to find the dot product of rv and cv. Does the order of the arguments matter?



In [5]:

    
dot(rv,cv)









    Out[5]:





array([11])



In [6]:

    
dot(cv,rv)









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-c47554307e31> in <module>
----> 1 dot(cv,rv)

<__array_function__ internals> in dot(*args, **kwargs)

ValueError: shapes (2,1) and (2,) not aligned: 1 (dim 1) != 2 (dim 0)

2) Use the function outer(vector1, vector2) to find the outer product of rv and cv. Does the order of the arguments matter?



In [7]:

    
outer(rv,cv)









    Out[7]:





array([[3, 4],
       [6, 8]])



In [8]:

    
outer(cv,rv)









    Out[8]:





array([[3, 6],
       [4, 8]])

II. Complex vectors



In [9]:

    
# Complex numbers in python have a j term:
a = 1+2j



In [10]:

    
v1 = array([1+2j, 3+2j, 5+1j, 4+0j])

The complex conjugate changes the sign of the imaginary part:



In [11]:

    
v1.conjugate()









    Out[11]:





array([1.-2.j, 3.-2.j, 5.-1.j, 4.-0.j])

3) Use dot() and .conjugate() to find the dot product of v1 and it's own conjugate:



In [12]:

    
dot(v1.conjugate(),v1)









    Out[12]:





(60+0j)

III. Matrices



In [13]:

    
# a two-dimensional array
m1 = array([[2,1],[2,1]])
m1









    Out[13]:





array([[2, 1],
       [2, 1]])



In [14]:

    
# can find transpose with the T method:
m1.T









    Out[14]:





array([[2, 2],
       [1, 1]])



In [15]:

    
# find the eigenvalues and eigenvectors of a matrix:
eig(m1)









    Out[15]:





(array([3., 0.]), array([[ 0.70710678, -0.4472136 ],
        [ 0.70710678,  0.89442719]]))

Can also use the matrix type which is like array but restricts to 2D. Also, matrix adds .H and .I methods for hermitian and inverse, respectively. For more information, see Stack Overflow question #4151128



In [16]:

    
m2 = matrix( [[2,1],[2,1]])



In [17]:

    
m2.H









    Out[17]:





matrix([[2, 2],
        [1, 1]])



In [18]:

    
eig(m2)









    Out[18]:





(array([3., 0.]), matrix([[ 0.70710678, -0.4472136 ],
         [ 0.70710678,  0.89442719]]))



In [19]:

    
# use a question mark to get help on a command
eig?

Examples:

Example 1.4

Find the eigenvalues and eigenvectors of M = ([0,1],[-2,3]])



In [20]:

    
M14 = array([[0,1],[-2,3]])



In [21]:

    
eig(M14)









    Out[21]:





(array([1., 2.]), array([[-0.70710678, -0.4472136 ],
        [-0.70710678, -0.89442719]]))

Interpret this result: the two eigenvalues are 1 and 2 the eigenvectors are strange decimals, but we can check them against the stated solution:



In [22]:

    
1/sqrt(2)  # this is the value for both entries in the first eigenvector









    Out[22]:





0.7071067811865475



In [23]:

    
1/sqrt(5)  # this is the first value in the second eigenvector









    Out[23]:





0.4472135954999579



In [24]:

    
2/sqrt(5)  # this is the second value in the second eigenvector









    Out[24]:





0.8944271909999159



In [25]:

    
eigvals(M14)









    Out[25]:





array([1., 2.])

Signs are opposite compared to the book, but it turns out that (-) doesn't matter in the interpretation of eigenvectors: only "direction" matters (the relative size of the entries).

Example: Problem 1.16 using Python functions



In [26]:

    
M16 = array([[0,-1j],[1j,0]])



In [27]:

    
evals, evecs = eig(M16)



In [28]:

    
evecs









    Out[28]:





array([[-0.        -0.70710678j,  0.70710678+0.j        ],
       [ 0.70710678+0.j        ,  0.        -0.70710678j]])



In [29]:

    
evecs[:,0]









    Out[29]:





array([-0.        -0.70710678j,  0.70710678+0.j        ])



In [30]:

    
evecs[:,1]









    Out[30]:





array([0.70710678+0.j        , 0.        -0.70710678j])



In [31]:

    
dot(evecs[:,0].conjugate(),evecs[:,1])









    Out[31]:





-1.6653345369377348e-16j

Part 2: Using QuTiP

Keeping track of row and column vectors in Ipython is somewhat artificial and tedious. The QuTiP library is designed to take care of many of these headaches



In [32]:

    
from qutip import *



In [33]:

    
# Create a row vector:
qv = Qobj([[1,2]])
qv









    Out[33]:





Quantum object: dims = [[1], [2]], shape = (1, 2), type = bra\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 2.0\\\end{array}\right)\end{equation*}



In [34]:

    
# Find the corresponding column vector
qv.dag()









    Out[34]:





Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}1.0\\2.0\\\end{array}\right)\end{equation*}



In [35]:

    
qv2 = Qobj([[1+2j,4-1j]])
qv2









    Out[35]:





Quantum object: dims = [[1], [2]], shape = (1, 2), type = bra\begin{equation*}\left(\begin{array}{*{11}c}(1.0+2.0j) & (4.0-1.0j)\\\end{array}\right)\end{equation*}



In [36]:

    
qv2.dag()









    Out[36]:





Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}(1.0-2.0j)\\(4.0+1.0j)\\\end{array}\right)\end{equation*}

Vector products in QuTiP

Only need to know one operator: "*" The product will depend on the order, either inner or outer



In [37]:

    
qv2*qv2.dag()  # inner product (dot product)









    Out[37]:





Quantum object: dims = [[1], [1]], shape = (1, 1), type = bra\begin{equation*}\left(\begin{array}{*{11}c}22.0\\\end{array}\right)\end{equation*}



In [38]:

    
qv2.dag()*qv2  # outer product









    Out[38]:





Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}5.0 & (2.0-9.0j)\\(2.0+9.0j) & 17.0\\\end{array}\right)\end{equation*}

Matrix in QuTiP



In [39]:

    
qm = Qobj([[1,2],[2,1]])
qm









    Out[39]:





Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 2.0\\2.0 & 1.0\\\end{array}\right)\end{equation*}



In [40]:

    
qm.eigenenergies()  # in quantum (as we will learn) eigenvalues often correspond to energy levels









    Out[40]:





array([-1.,  3.])



In [41]:

    
evals, evecs = qm.eigenstates()



In [42]:

    
evecs









    Out[42]:





array([Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
Qobj data =
[[-0.70710678]
 [ 0.70710678]],
       Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
Qobj data =
[[0.70710678]
 [0.70710678]]], dtype=object)



In [43]:

    
evecs[1]









    Out[43]:





Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}0.707\\0.707\\\end{array}\right)\end{equation*}

Some additional examples we have already done by hand:

Problem 1.2 using the hist() function.

Create a histogram of the data and from it estimate the probability distribution. Use min and max values for the range of 5 and 14, and use 5 bins (5 and 6 go in the first bin).



In [44]:

    
# Solution
data = [10,13,14,14,6,8,7,9,12,14,13,11,10,7,7]
# Fill in the hist() function:
n, bins, patches = hist(data, bins=9, range=(5,14))

Problem 1.8

Find the constant $c$ that normalizes the probability density: $$p(x) = \begin{cases} ce^{-ax}, & x \geq 0 \\[2ex] 0 & x < 0 \end{cases} $$

Hint: using sympy, we can calculate the relevant integral. The conds='none' asks the solver to ignore any strange conditions on the variables in the integral. This is fine for most of our integrals. Usually the variables are real and well-behaved numbers.



In [51]:

    
# Partial Solution:
from sympy import *
c,a,x = symbols("c a x")
Q.positive((c,a))
first = integrate( c*exp(-a*x) ,(x,0,oo),conds='none')
first









    Out[51]:





$\displaystyle \frac{c}{a}$



In [55]:

    
check = integrate( a*exp(-a*x) ,(x,0,oo),conds='none')
check









    Out[55]:





$\displaystyle 1$



In [ ]: