In [1]:

    

from pylab import *
# the harmonic oscillator potential with hbar = m = omega =1
def V(x):
    return 0.5 * x ** 2
# computational grid running from -5 to 5 with 512 grid points
N = 512
x0, x1 = -5., 5.
x = linspace(x0, x1, N)
dx = (x1 - x0) / (N -1) # size of spartial grid spacing
dt = 1. / 64
# temporal step size
# create Hamiltonian matrix and calculate its eigenvectors and
# eigenvalues
# use finite differences to approximate derivatives to 2 nd order
H =(-0.5 * (diag(ones(N - 1), -1) - diag(2 * ones( N )) +
    diag(ones(N - 1), + 1)) / (dx ** 2) + diag(V(x))) * dt
Lambda, Q = eigh(H)
Q_inv = Q.conj().transpose() # inverse of Q
# Gaussian wave packet with mean momentum one as initial condition
Psi0 = 1. / (2 * pi) ** (0.25) * exp(-x ** 2 / 4 + 1 j * x)
Psi = Psi0.copy()
# propagate 512 time steps
t = 0
for k in range (0, 512):
    clf()
    plot (x, Psi.real, '--' , color = 'r ', label=r'$\mathrm{Re}\, \Psi(x)$')
    plot (x, Psi.imag, ':', color ='r ', label =r'$\mathrm{Im}\, \Psi(x)$')
    plot (x , Psi . real **2+ Psi . imag **2 ,
    color = ' #266 bbd ', label = r '$ |\ Psi ( x )|^2 $ ')
    gca (). set_xlim ( x0 , x1 )
    gca (). set_ylim ( -1 , 1)
    xlabel (r ' $x$ ')
    legend ( loc = ' upper left ')
    show ()
    pause (0.01)
    Psi = dot ( Q_inv , Psi )
    # expand into eigen - functions
    Psi *= exp ( -1 j* Lambda ) # propagate eigen - functions
    Psi = dot (Q , Psi )
    # superpose eigen - functions
    t += dt
figure ()
plot (x , abs ( Psi0 )**2 , label = r '$t =0 $ ')
plot (x , abs ( Psi )**2 , ' -- ', label =( r '$t ={: g } $ '. format ( t )))
xlabel (r ' $x$ ')
ylabel (r '$ |\ Psi (x , t )|^2 $ ')
legend ()
tight_layout ()
show ()
savefig ( ' figures / harmonic_os


    

    




  File "<ipython-input-1-7a0bb375a9b6>", line 20
    Psi0 = 1. / (2 * pi) ** (0.25) * exp(-x ** 2 / 4 + 1 j * x)
                                                         ^
SyntaxError: invalid syntax

In [2]:

    

from pylab import *
%matplotlib inline
# the harmonic oscillator potential with hbar = m = omega =1

def V(x):
    return x**2
# computational grid running from -5 to 5 with 512 grid points
V0 = 1.
N = 512
x0, x1 = -5., 5.
x = linspace(x0, x1, N)
dx = (x1 - x0) / (N -1) # size of spartial grid spacing
dt = 1. / 64 * 8
# temporal step size
# create Hamiltonian matrix and calculate its eigenvectors and
# eigenvalues
# use finite differences to approximate derivatives to 2 nd order
H =(-0.5 * (diag(ones(N - 1), -1) - diag(2 * ones( N )) +
    diag(ones(N - 1), + 1)) / (dx ** 2) + diag(V(x))) * dt
Lambda, Q = eigh(H)
Q_inv = Q.conj().transpose() # inverse of Q
# Gaussian wave packet with mean momentum one as initial condition
Psi0 = 1. / (2 * pi) ** (0.25) * exp(-x ** 2 / 4 + 1j * x)
Psi = Psi0.copy()

t=0
dt = 100
for k in range (0, 25):
    #plot(x, Psi.real, '--' , color = 'r', label=r'$\mathrm{Re}\, \Psi(x)$')
    #plot(x, Psi.imag, ':', color ='r', label =r'$\mathrm{Im}\, \Psi(x)$')
    plot(x ,Psi.real ** 2 + Psi.imag ** 2, color='#266bbd', label=r'$|\Psi(x)|^2$')
    Psi = dot(Q_inv, Psi)
    # expand into eigen - functions
    Psi *= exp(-1j*Lambda) # propagate eigen - functions
    Psi = dot(Q, Psi)
    # superpose eigen - functions
    t += dt


    

In [14]:

    

Q.shape
plt.plot(Q[, :]*Psi)


    

    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)






    
Out[14]:





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






    

In [4]:

    

cond(H)**100


    

    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: overflow encountered in double_scalars
  """Entry point for launching an IPython kernel.






    
Out[4]:





inf

In [ ]: