In [13]:

    

import numpy as np
import matplotlib.pyplot as plt
from numpy import linalg as LA
%matplotlib inline


    

In [14]:

    

def V(x, w=.0001, l=10., V_0=200.):
    return -V_0 / 2. * (np.tanh((x + l / 2) / w) \
                      - np.tanh((x - l / 2) / w)) + V_0


    

In [15]:

    

x = np.linspace(-20, 20, 1000)
for k in - np.arange(-1, 2):
    plt.plot(x, V(x, 10. ** k, 10, 10),\
             label='w = 10^'+str(k))
plt.legend()

plt.figure()
for k in np.arange(1, 10, 3):
    plt.plot(x, V(x, 1, 10, k),\
             label='V_0 = '+str(k))
plt.legend()

plt.figure()
for k in np.arange(1, 10, 3):
    plt.plot(x, V(x, 1, k, 4),\
             label='l = '+str(k))
plt.legend()


    

    
Out[15]:





<matplotlib.legend.Legend at 0x7f1f48233828>






    













    













    

In [16]:

    

# define the grid
N = 2 ** 10 # number of grid points
x0 , x1 = -20. , 20. # grid boundaries
dx =(x1 - x0) / (N -1) # grid spacing
x = np.linspace (x0 , x1 , N) # grid points
# setup the Hamiltonian matrix
H = -1 / (2 * dx ** 2) * (np.diag(np.ones(N - 1) , -1) - 2 * np.diag(np.ones(N)) + \
    np.diag(np.ones(N - 1), 1)) + np.diag(V(x))
# compute eigenvalues
E = LA.eigvalsh(H)
# plot first eigenvalues
n =100
ns =  np.array (range (0 , n ))
plt.plot( ns, np.sort (E)[0: n ] ,'+', label ='eigenvalues of discretized Hamiltonian ')
plt.plot(ns[E_teo(ns) < 200], E_teo(ns)[E_teo(ns) < 200])
plt.grid()
sum(E < 200)


    

    
Out[16]:





66






    

In [17]:

    

E , Psi_E = LA.eigh(H)
Psi_E /= np.sqrt (dx)
for n in np.arange(2):
    c = np.zeros ( n +1)
    c [n]=1
    plt.plot (x ,- Psi_E [: , n ] , label ='$n =% i$ ' % n )
    
plt.legend()
def E_teo(n):
    return n** 2 * np.pi ** 2 / (2 * (10**2))
print (E[:5])
print (E_teo(np.arange(5))/78)


    

    




[ 0.04821121  0.19283743  0.43385644  0.77123123  1.20490998]
[ 0.          0.00063267  0.00253067  0.005694    0.01012267]






    

In [18]:

    

plt.plot(x, V(x))


    

    
Out[18]:





[<matplotlib.lines.Line2D at 0x7f1f472fee48>]






    

In [30]:

    

c=137.0359998

N = 2 ** 9 # number of grid points
x0 , x1 = -20. , 20. # grid boundaries
dx =(x1 - x0)/(N -1) # grid spacing
x = np.linspace (x0 , x1 , N) # grid points
# setup the Hamiltonian matrix
sigma_1=np.array([[0, 1], [1, 0]])
sigma_3=np.array([[1, 0], [0, -1]])
p=np.diag(-np.ones(N-1), -1) + np.diag(np.ones(N-1), +1)
p[0, N-1]=-1
p[N-1, 0]=1
p/=2*dx
H_free=c*(-1j)*np.kron(p, sigma_1) + np.kron(np.eye(N), sigma_3*c**2)
H=H_free+np.kron(np.diag(V(x)), np.eye(2))

# compute eigenvalues
E = LA.eigvalsh(H)
# plot first eigenvalues
n =100
ns =  np.array (range (0 , n ))
plt.plot( ns, E[:n] ,'+', label ='eigenvalues of discretized Hamiltonian ')
#plt.plot(ns[E_teo(ns) < 200], E_teo(ns)[E_teo(ns) < 200])
plt.grid()
sum(E < 200)


    

    
Out[30]:





512






    

In [31]:

    

E


    

    
Out[31]:





array([-18860.24184989, -18860.24184989, -18860.10041793, ...,
        19060.26754386,  19060.28365803,  19060.28365803])

In [21]:

    

E , Psi_E = LA.eigh(H)
Psi_E /= np.sqrt (dx)
for n in np.arange(2):
    c = np.zeros ( n +1)
    c [n]=1
    plt.plot (x ,- Psi_E [: , n ] , label ='$n =% i$ ' % n )


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-21-3e8512c33d07> in <module>()
      4     c = np.zeros ( n +1)
      5     c [n]=1
----> 6     plt.plot (x ,- Psi_E [: , n ] , label ='$n =% i$ ' % n )

/home/lorenzo/anaconda3/lib/python3.6/site-packages/matplotlib/pyplot.py in plot(*args, **kwargs)
   3315                       mplDeprecation)
   3316     try:
-> 3317         ret = ax.plot(*args, **kwargs)
   3318     finally:
   3319         ax._hold = washold

/home/lorenzo/anaconda3/lib/python3.6/site-packages/matplotlib/__init__.py in inner(ax, *args, **kwargs)
   1896                     warnings.warn(msg % (label_namer, func.__name__),
   1897                                   RuntimeWarning, stacklevel=2)
-> 1898             return func(ax, *args, **kwargs)
   1899         pre_doc = inner.__doc__
   1900         if pre_doc is None:

/home/lorenzo/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py in plot(self, *args, **kwargs)
   1404         kwargs = cbook.normalize_kwargs(kwargs, _alias_map)
   1405 
-> 1406         for line in self._get_lines(*args, **kwargs):
   1407             self.add_line(line)
   1408             lines.append(line)

/home/lorenzo/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py in _grab_next_args(self, *args, **kwargs)
    405                 return
    406             if len(remaining) <= 3:
--> 407                 for seg in self._plot_args(remaining, kwargs):
    408                     yield seg
    409                 return

/home/lorenzo/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py in _plot_args(self, tup, kwargs)
    383             x, y = index_of(tup[-1])
    384 
--> 385         x, y = self._xy_from_xy(x, y)
    386 
    387         if self.command == 'plot':

/home/lorenzo/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py in _xy_from_xy(self, x, y)
    242         if x.shape[0] != y.shape[0]:
    243             raise ValueError("x and y must have same first dimension, but "
--> 244                              "have shapes {} and {}".format(x.shape, y.shape))
    245         if x.ndim > 2 or y.ndim > 2:
    246             raise ValueError("x and y can be no greater than 2-D, but have "

ValueError: x and y must have same first dimension, but have shapes (512,) and (1024,)





    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [36]:

    

from pylab import *

close('all')

def V(x, V0, l, w):
    return V0/2*(tanh((x+l/2)/w) - tanh((x-l/2)/w))

c=137.0359998           # speed of light

# define the grid
N=512                   # number of grid points
x0, x1=-0.5, 0.5        # grid boundaries
dx=(x1-x0)/(N-1)        # grid spacing
x=linspace(x0, x1, N)   # grid points

l=3.2/c
w=0.3/c

# setup the Hamiltonian matrix
sigma_1=array([[0, 1], [1, 0]])
sigma_3=array([[1, 0], [0, -1]])
p=diag(-ones(N-1), -1) + diag(ones(N-1), +1)
p[0, N-1]=-1
p[N-1, 0]=1
p/=2*dx
H_free=c*(-1j)*kron(p, sigma_1) + kron(eye(N), sigma_3*c**2)


H=H_free+kron(diag(V(x, -1000, l, w)), eye(2))
E=eigvalsh(H)
plt.plot(E)
xlabel('$V_0/(m_0c^2)$')
ylabel('$E/(m_0c^2)$')
gca().set_yticks(arange(-8, 4))
tight_layout()
show()


    

In [ ]:

    




    

In [ ]: