In [3]:

    

from pylab import *

close('all')

def V(x, Z):
    return -Z / sqrt(2 / Z ** 2 + x ** 2)

# define the grid
N= 2 ** 9                   # number of grid points
x0, x1=-10., 10.        # grid boundaries
dx=(x1-x0)/(N-1)        # grid spacing
x=linspace(x0, x1, N)   # grid points

Z = 1
plot(x, V(x, Z), label='$V(x, Z)$')

ylabel('$V(x)$')
legend()
tight_layout()
show()

figure()  # add new figure


# setup the Hamiltonian matrix
H=diag(ones(N-1), -1) - 2*diag(ones(N)) + diag(ones(N-1), 1)
H[0, N-1]=H[N-1, 0]=1
H*=-1/(2*dx**2)
H+=diag(V(x, Z))

# compute eigenvalues
E=eigvalsh(H)
plot(E)
xlabel('index')
ylabel('$E$')
tight_layout()
show()

figure()  # add new figure
H=diag(ones(N-1), -1) - 2*diag(ones(N)) + diag(ones(N-1), 1)
H[0, N-1]=H[N-1, 0]=1
H*=-1/(2*dx**2)
H+=diag(V(x, Z))
E=eigvalsh(H)
E_bound=E[E<0]
plot(ones_like(E_bound), E_bound, 'k.')
xlabel('$V_0$')
ylabel('$E$')
tight_layout()
show()

figure()  # add new figure

V0=-2  # potential depth 
l=5    # width of well
w=0.5  # width of step-like function

# setup the Hamiltonian matrix
H=diag(ones(N-1), -1) - 2*diag(ones(N)) + diag(ones(N-1), 1)
#H[0, N-1]=H[N-1, 0]=1
H*=-1/(2*dx**2)
H+=diag(V(x, Z))

# compute eigenvectors
E, Psi=eigh(H)
plot(x, Psi[:, 0]/sqrt(dx), label=r'$\Psi_0(x)$')
plot(x, Psi[:, 1]/sqrt(dx), label=r'$\Psi_1(x)$')
plot(x, Psi[:, 2]/sqrt(dx), label=r'$\Psi_2(x)$')
plot(x, Psi[:, 3]/sqrt(dx), label=r'$\Psi_3(x)$')
#plot(x, Psi[:, 40]/sqrt(dx), label=r'$\Psi_{40}(x)$')
xlabel('$x$')
legend(ncol=2)
tight_layout()
show()


    

In [2]:

    

def E(t, E0, omega, n):
    return E0 * sin(omega * t) * sin(omega * t / (2 * n))

E0 = 1.
n = 5
t = 1.
omega = np.linspace(-np.pi, 2 * np.pi * (n + .5), 10000)
plt.plot(omega, E(t, E0, omega, 3 / 2))
show()
plt.plot(omega, E(t, E0, omega, 5))
show()


    

In [3]:

    

def V_e(x, t, E0, omega, n):
    return - x * E(t, E0, omega, n)


    

In [4]:

    

omega = .02
n = 3 / 2
E0 = .05

for t in linspace(0, 50 * 4 * pi, 9):
    plt.figure()
    plt.plot(V(x, 1) + V_e(x,t, E0, omega, 5.))
plt.show()


    

In [26]:

    

plt.cla()

N =2 ** 7
x0 , x1 = -7.5 , 7.5
x = linspace (x0 , x1 , N)
dx =( x1 - x0 )/( N -1)
# size of spatial grid spacing
dt =1. / 8
# temporal step size
# Gaussian wave packet with mean momentum 2 as initial condition

# constuct momentum grid
p = fftfreq (N , d= dx /(2* pi ))


# setup the Hamiltonian matrix
H=diag(ones(N-1), -1) - 2*diag(ones(N)) + diag(ones(N-1), 1)
H[0, N-1]=H[N-1, 0]=1
H*=-1/(2*dx**2)
H+=diag(V(x, Z))

# compute eigenvectors
_, Psi_tot=eigh(H)
Psi = Psi_tot[:, 0].astype('complex')

#prop_cycle = plt.rcParam['axes.prop_cycle']
#colors = prop_cycle.by_key()['color']


def V_tot(x, t):
    return V(x, 1) + V_e(x,t, E0, omega, 5.)

#U_1 = exp (-0.5*1j*V(x)* dt )
U_2 = exp (-1j * 0.5* p **2* dt )
for k in range (0 , 1000):

    #plot(x, Psi.real , '--')
    #plot(x, Psi.imag, ':')
    plot(x, abs(Psi))
    xlim(x0 , x1 )
    #ylim( -1 , 1)
    xlabel (r' $x$ ')
    #legend ( loc = ' upper left ')
    tight_layout ()
    U_1 = exp (-0.5*1j*V_tot(x, k * dt)* dt )
    Psi *= U_1
    # apply U_1 in real space
    Psi = fft(Psi)
    # go to Fourier space
    Psi *= U_2
    # apply U_2 in Fourier space
    Psi = ifft ( Psi ) # go to real space
    Psi *= U_1
plt.show()


    

In [27]:

    

omega = .02
n = 3 / 2
E0 = .05

-x * E(t, E0, omega, n)


    

    
Out[27]:





array([  1.59086286e-16,   1.56580990e-16,   1.54075694e-16,
         1.51570398e-16,   1.49065102e-16,   1.46559807e-16,
         1.44054511e-16,   1.41549215e-16,   1.39043919e-16,
         1.36538623e-16,   1.34033327e-16,   1.31528032e-16,
         1.29022736e-16,   1.26517440e-16,   1.24012144e-16,
         1.21506848e-16,   1.19001552e-16,   1.16496257e-16,
         1.13990961e-16,   1.11485665e-16,   1.08980369e-16,
         1.06475073e-16,   1.03969777e-16,   1.01464482e-16,
         9.89591857e-17,   9.64538898e-17,   9.39485940e-17,
         9.14432981e-17,   8.89380023e-17,   8.64327065e-17,
         8.39274106e-17,   8.14221148e-17,   7.89168189e-17,
         7.64115231e-17,   7.39062273e-17,   7.14009314e-17,
         6.88956356e-17,   6.63903397e-17,   6.38850439e-17,
         6.13797481e-17,   5.88744522e-17,   5.63691564e-17,
         5.38638605e-17,   5.13585647e-17,   4.88532689e-17,
         4.63479730e-17,   4.38426772e-17,   4.13373814e-17,
         3.88320855e-17,   3.63267897e-17,   3.38214938e-17,
         3.13161980e-17,   2.88109022e-17,   2.63056063e-17,
         2.38003105e-17,   2.12950146e-17,   1.87897188e-17,
         1.62844230e-17,   1.37791271e-17,   1.12738313e-17,
         8.76853544e-18,   6.26323960e-18,   3.75794376e-18,
         1.25264792e-18,  -1.25264792e-18,  -3.75794376e-18,
        -6.26323960e-18,  -8.76853544e-18,  -1.12738313e-17,
        -1.37791271e-17,  -1.62844230e-17,  -1.87897188e-17,
        -2.12950146e-17,  -2.38003105e-17,  -2.63056063e-17,
        -2.88109022e-17,  -3.13161980e-17,  -3.38214938e-17,
        -3.63267897e-17,  -3.88320855e-17,  -4.13373814e-17,
        -4.38426772e-17,  -4.63479730e-17,  -4.88532689e-17,
        -5.13585647e-17,  -5.38638605e-17,  -5.63691564e-17,
        -5.88744522e-17,  -6.13797481e-17,  -6.38850439e-17,
        -6.63903397e-17,  -6.88956356e-17,  -7.14009314e-17,
        -7.39062273e-17,  -7.64115231e-17,  -7.89168189e-17,
        -8.14221148e-17,  -8.39274106e-17,  -8.64327065e-17,
        -8.89380023e-17,  -9.14432981e-17,  -9.39485940e-17,
        -9.64538898e-17,  -9.89591857e-17,  -1.01464482e-16,
        -1.03969777e-16,  -1.06475073e-16,  -1.08980369e-16,
        -1.11485665e-16,  -1.13990961e-16,  -1.16496257e-16,
        -1.19001552e-16,  -1.21506848e-16,  -1.24012144e-16,
        -1.26517440e-16,  -1.29022736e-16,  -1.31528032e-16,
        -1.34033327e-16,  -1.36538623e-16,  -1.39043919e-16,
        -1.41549215e-16,  -1.44054511e-16,  -1.46559807e-16,
        -1.49065102e-16,  -1.51570398e-16,  -1.54075694e-16,
        -1.56580990e-16,  -1.59086286e-16])

In [41]:

    




    

    
Out[41]:





0.0

In [20]:

    

plt.plot(x, abs(Psi))
plot(x, Psi.imag, ':')
plt.show()


    

In [32]:

    

#_, Psi_tot=eigh(H)
sum(abs(Psi_tot[:, 2] * dx))


    

    
Out[32]:





1.2214055110259801

In [ ]: