1. sympy
  2. examples
  3. notebooks

Example Notebook for sho1d.py

Import the sho1d.py file as well as the test_sho1d.py file

In [1]:
%load_ext sympy.interactive.ipythonprinting 
from sympy import *
from IPython.display import display_pretty
from sympy.physics.quantum import *
from sympy.physics.quantum.sho1d import *
from sympy.physics.quantum.tests.test_sho1d import *

The sympy.interactive.ipythonprinting extension is already loaded. To reload it, use:
  %reload_ext sympy.interactive.ipythonprinting

Printing Of Operators

Create a raising and lowering operator and make sure they print correctly

In [2]:
ad = RaisingOp('a')
a = LoweringOp('a')
In [3]:
ad
Out[3]:
RaisingOp(a)
In [4]:
a
Out[4]:
a
In [5]:
print latex(ad)
print latex(a)

a^{\dag}
a
In [6]:
display_pretty(ad)
display_pretty(a)

RaisingOp(a)
a
In [7]:
print srepr(ad)
print srepr(a)

RaisingOp(Symbol('a'))
LoweringOp(Symbol('a'))
In [8]:
print repr(ad)
print repr(a)

RaisingOp(a)
a

Printing of States

Create a simple harmonic state and check its printing

In [9]:
k = SHOKet('k')
b = SHOBra('b')
In [10]:
k
Out[10]:
|k>
In [11]:
b
Out[11]:
<b|
In [12]:
print pretty(k)
print pretty(b)

❘k⟩
⟨b❘
In [13]:
print latex(k)
print latex(b)

{\left|k\right\rangle }
{\left\langle b\right|}
In [14]:
print srepr(k)
print srepr(b)

SHOKet(Symbol('k'))
SHOBra(Symbol('b'))

Properties

Take the dagger of the raising and lowering operators. They should return eachother.

In [15]:
Dagger(ad)
Out[15]:
a
In [16]:
Dagger(a)
Out[16]:
RaisingOp(a)

Check Commutators of the raising and lowering operators

In [17]:
Commutator(ad,a).doit()
Out[17]:
-1
In [18]:
Commutator(a,ad).doit()
Out[18]:
1

Take a look at the dual states of the bra and ket

In [19]:
k.dual
Out[19]:
<k|
In [20]:
b.dual
Out[20]:
|b>

Taking the InnerProduct of the bra and ket will return the KroneckerDelta function

In [21]:
InnerProduct(b,k).doit()
Out[21]:
KroneckerDelta(k, b)

Take a look at how the raising and lowering operators act on states. We use qapply to apply an operator to a state

In [22]:
qapply(ad*k)
Out[22]:
sqrt(k + 1)*|k + 1>
In [23]:
qapply(a*k)
Out[23]:
sqrt(k)*|k - 1>

But the states may have an explicit energy level. Let's look at the ground and first excited states

In [24]:
kg = SHOKet(0)
kf = SHOKet(1)
In [25]:
qapply(ad*kg)
Out[25]:
|1>
In [26]:
qapply(ad*kf)
Out[26]:
sqrt(2)*|2>
In [27]:
qapply(a*kg)
Out[27]:
0
In [28]:
qapply(a*kf)
Out[28]:
|0>

Notice that akg is 0 and akf is the |0> the ground state.

NumberOp & Hamiltonian

Let's look at the Number Operator and Hamiltonian Operator

In [29]:
k = SHOKet('k')
ad = RaisingOp('a')
a = LoweringOp('a')
N = NumberOp('N')
H = Hamiltonian('H')

The number operator is simply expressed as ad*a

In [30]:
N().rewrite('a').doit()
Out[30]:
RaisingOp(a)*a

The number operator expressed in terms of the position and momentum operators

In [31]:
N().rewrite('xp').doit()
Out[31]:
-1/2 + (m**2*omega**2*X**2 + Px**2)/(2*hbar*m*omega)

It can also be expressed in terms of the Hamiltonian operator

In [32]:
N().rewrite('H').doit()
Out[32]:
-1/2 + H/(hbar*omega)

The Hamiltonian operator can be expressed in terms of the raising and lowering operators, position and momentum operators, and the number operator

In [33]:
H().rewrite('a').doit()
Out[33]:
hbar*omega*(1/2 + RaisingOp(a)*a)
In [34]:
H().rewrite('xp').doit()
Out[34]:
(m**2*omega**2*X**2 + Px**2)/(2*m)
In [35]:
H().rewrite('N').doit()
Out[35]:
hbar*omega*(1/2 + N)

The raising and lowering operators can also be expressed in terms of the position and momentum operators

In [36]:
ad().rewrite('xp').doit()
Out[36]:
sqrt(2)*(m*omega*X - I*Px)/(2*sqrt(hbar)*sqrt(m*omega))
In [37]:
a().rewrite('xp').doit()
Out[37]:
sqrt(2)*(m*omega*X + I*Px)/(2*sqrt(hbar)*sqrt(m*omega))

Properties

Let's take a look at how the NumberOp and Hamiltonian act on states

In [38]:
qapply(N*k)
Out[38]:
k*|k>

Apply the Number operator to a state returns the state times the ket

In [39]:
ks = SHOKet(2)
qapply(N*ks)
Out[39]:
2*|2>
In [40]:
qapply(H*k)
Out[40]:
hbar*k*omega*|k> + hbar*omega*|k>/2

Let's see how the operators commute with each other

In [41]:
Commutator(N,ad).doit()
Out[41]:
RaisingOp(a)
In [42]:
Commutator(N,a).doit()
Out[42]:
-a
In [43]:
Commutator(N,H).doit()
Out[43]:
0

Representation

We can express the operators in NumberOp basis. There are different ways to create a matrix in Python, we will use 3 different ways.

Sympy

In [44]:
represent(ad, basis=N, ndim=4, format='sympy')
Out[44]:
[0,       0,       0, 0]
[1,       0,       0, 0]
[0, sqrt(2),       0, 0]
[0,       0, sqrt(3), 0]

Numpy

In [45]:
represent(ad, basis=N, ndim=5, format='numpy')
Out[45]:
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  1.41421356,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  1.73205081,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  2.        ,  0.        ]])

Scipy.Sparse

In [46]:
represent(ad, basis=N, ndim=4, format='scipy.sparse', spmatrix='lil')
Out[46]:
<4x4 sparse matrix of type '<type 'numpy.float64'>'
	with 3 stored elements in Compressed Sparse Row format>
In [47]:
print represent(ad, basis=N, ndim=4, format='scipy.sparse', spmatrix='lil')

  (1, 0)	1.0
  (2, 1)	1.41421356237
  (3, 2)	1.73205080757

The same can be done for the other operators

In [48]:
represent(a, basis=N, ndim=4, format='sympy')
Out[48]:
[0, 1,       0,       0]
[0, 0, sqrt(2),       0]
[0, 0,       0, sqrt(3)]
[0, 0,       0,       0]
In [49]:
represent(N, basis=N, ndim=4, format='sympy')
Out[49]:
[0, 0, 0, 0]
[0, 1, 0, 0]
[0, 0, 2, 0]
[0, 0, 0, 3]
In [50]:
represent(H, basis=N, ndim=4, format='sympy')
Out[50]:
[hbar*omega/2,              0,              0,              0]
[           0, 3*hbar*omega/2,              0,              0]
[           0,              0, 5*hbar*omega/2,              0]
[           0,              0,              0, 7*hbar*omega/2]

Ket and Bra Representation

In [51]:
k0 = SHOKet(0)
k1 = SHOKet(1)
b0 = SHOBra(0)
b1 = SHOBra(1)
In [52]:
print represent(k0, basis=N, ndim=5, format='sympy')

[1]
[0]
[0]
[0]
[0]
In [53]:
print represent(k1, basis=N, ndim=5, format='sympy')

[0]
[1]
[0]
[0]
[0]
In [54]:
print represent(b0, basis=N, ndim=5, format='sympy')

[1, 0, 0, 0, 0]
In [55]:
print represent(b1, basis=N, ndim=5, format='sympy')

[0, 1, 0, 0, 0]