In [2]:

    

from pylab import *
from copy import copy
%matplotlib inline


    

In [5]:

    

E1, E2, E3 = 0., 20., 0.
V12, V23 = 1., 1.

psi0 = array([1, 0, 0], dtype='complex')

Nt = 100000
psi = np.zeros((Nt, 3), dtype='complex')
psi[0, :] = psi0

times = np.linspace(0, 20, Nt)
for E2 in np.arange(4) * 20:
    H = array([[E1, V12, 0],
                [V12, E2, V23],
                [0, V23, E3]])
    lambd, Q = eigh(H)
    Q_inv = Q.T.conj()
    for it in range(1, Nt):
        psi[it, :] = Q_inv.dot(psi0)
        psi[it, :] = np.diag(np.exp(-1j * lambd * times[it])).dot(psi[it, :])
        psi[it, :] = Q.dot(psi[it, :])
    plt.plot(times, abs(psi) ** 2)
    plt.figure()
    times = np.linspace(0, 200, Nt)


    

    













    













    













    













    






<matplotlib.figure.Figure at 0x7ff107afb2b0>

In [7]:

    

plt.plot(times, psi[:, 0].imag ** 2 + psi[:, 0].real ** 2)
plt.plot(times, psi[:, 2].imag ** 2 + psi[:, 2].real ** 2)


y = cos(V12 ** 2 / E2 * times) ** 2
plt.plot(times, y)
y = sin(V12 ** 2 / E2 * times) ** 2
plt.plot(times, y)


    

    
Out[7]:





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






    

In [16]:

    

y = cos(V12 ** 2 / E2 * times) ** 2
plt.plot(times, y)
plt.plot(times, psi[:, 0].imag ** 2 + psi[:, 0].real ** 2)
plt.xlim([0, 5])
plt.ylim([.975, 1])


    

    
Out[16]:





(0.975, 1)






    

In [9]:

    

plt.plot((psi[:, 1].imag ** 2 + psi[:, 1].real ** 2)[:100])


    

    
Out[9]:





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






    

In [36]:

    

Q_inv.dot(D.dot(Q))


    

    
Out[36]:





array([[ 0.33333333,  1.22474487,  0.23570226],
       [ 1.22474487,  0.5       ,  0.8660254 ],
       [ 0.23570226,  0.8660254 ,  0.16666667]])

In [78]:

    

np.linspace(0, 10, Nt)[1]


    

    
Out[78]:





0.0010001000100010001

In [87]:

    

1 / 64


    

    
Out[87]:





0.015625

In [109]:

    

Q


    

    
Out[109]:





array([[  5.08756637e-01,  -7.07106781e-01,   4.91087247e-01],
       [ -6.94502245e-01,  -2.08166817e-17,   7.19490536e-01],
       [  5.08756637e-01,   7.07106781e-01,   4.91087247e-01]])

In [110]:

    

Q_inv


    

    
Out[110]:





array([[  5.08756637e-01,  -6.94502245e-01,   5.08756637e-01],
       [ -7.07106781e-01,  -2.08166817e-17,   7.07106781e-01],
       [  4.91087247e-01,   7.19490536e-01,   4.91087247e-01]])

In [60]:

    

E1, E2, E3 = 0., 1., 0.
V12, V23 = 100., 100.

H = array([[E1, V12, 0],
            [V12, E2, V23],
            [0, V23, E3]])
Nt  = 10000
dt = times[1]
H = H * dt
Lambda , Q= eigh (H)
psi0 = array([1, 0, 0], dtype='complex')

Q_inv =Q. conj (). transpose () # inverse of Q
# Gaussian wave packet with mean momentum one as initial condition
# propagate 512 time steps
t =0
psi = np.zeros((Nt, 3), dtype='complex')
for k in range (0 , Nt):
    psi[k, :] = copy(psi0)
    psi0 = dot ( Q_inv , psi0 )
    # expand into eigen - functions
    psi0 *= exp ( -1j* Lambda ) # propagate eigen - functions
    psi0 = dot (Q , psi0 )
    # superpose eigen - functions
    t += dt
    psi[k, :] = copy(psi0)


    

In [61]:

    

plt.plot(np.absolute(psi))


    

    
Out[61]:





[<matplotlib.lines.Line2D at 0x7f95d74a29b0>,
 <matplotlib.lines.Line2D at 0x7f95d74a2b70>,
 <matplotlib.lines.Line2D at 0x7f95d74a2d68>]






    

In [15]:

    

t = np.linspace(0, 1, Nt)
y = cos(V12**2 / * t) ** 2
plt.plot(t, y)


    

    
Out[15]:





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






    

In [152]:

    

psi0 = dot ( Q_inv , psi0 )
# expand into eigen - functions
psi0 *= exp ( -1j* Lambda ) # propagate eigen - functions
psi0 = dot (Q , psi0 )


    

In [153]:

    

dot ( Q_inv , psi0 )* exp ( -1j* Lambda )


    

    
Out[153]:





array([ 0.50875664 +1.08646273e-24j, -0.70710678 -7.99576416e-37j,
        0.49108725 -1.12555382e-24j])

In [73]:

    

from sympy import *
init_printing(False)


    

In [41]:

    

E, V = symbols('E V')


    

In [48]:

    

H =Matrix([[0, V, 0],
           [V, E, V],
           [0, V, 0]]).evalf(subs={E:1, V:100})


    

In [49]:

    

Q, lambd = H.diagonalize()


    

In [74]:

    

lambd


    

    
Out[74]:





array([ -9.95049384e-02,   1.94543417e-15,   2.00995049e+01])

In [79]:

    

7 + sqrt(7**2 + 8 * 9)


    

    
Out[79]:





$$18$$

In [80]:

    

x = symbols('x')


    

In [81]:

    

solve(x** 2 - x  * E - 2* V ** 2, x)


    

    
Out[81]:





$$\left [ \frac{E}{2} - \frac{1}{2} \sqrt{E^{2} + 8 V^{2}}, \quad \frac{E}{2} + \frac{1}{2} \sqrt{E^{2} + 8 V^{2}}\right ]$$

In [86]:

    

t = symbols('t')


    

In [127]:

    

D =Matrix([[0, 0, 0],
        [0, exp(-1j *lambd[1, 1] * t), 0],
        [0, 0, exp(-1j *lambd[2, 2] * t)]])


    

In [128]:

    

psit = Q * D * Q.T * np.array([1, 0, 0])


    

In [131]:

    

psit[2].evalf(subs={t:0})


    

    
Out[131]:





$$2.0$$

In [99]:

    

x.evalf(subs={x:2})


    

    
Out[99]:





$$2.0$$

In [113]:

    

Q.evalf(subs={E:1, V:10}) * Q.T.evalf(subs={E:1, V:10})  * array([1, 0,0])


    

    
Out[113]:





$$\left[\begin{matrix}9.0\\0.399999999999999\\7.0\end{matrix}\right]$$

In [15]:

    

E1, E2 = 1., 1.
V = 1.

H = array([[E1, V],
            [V, E2]])

psi0 = array([1, 0], dtype='complex')
#dt = np.linspace(0, 10, Nt)[1]
lambd, Q = eigh(H)
Q_inv = Q.T.conj()

Nt = 10000
psi = np.zeros((Nt, 2), dtype='complex')
psi[0, :] = psi0

times = np.linspace(0, 5, Nt)
for it in range(1, Nt):
    psi[it, :] = Q_inv.dot(psi0)
    psi[it, :] = np.diag(np.exp(-1j * lambd * times[it])).dot(psi[it, :])
    psi[it, :] = Q.dot(psi[it, :])


    

In [18]:

    

plt.plot(times, psi[:, 0].imag ** 2 + psi[:, 0].real ** 2)
plt.plot(times, cos(V * times) ** 2)


    

    
Out[18]:





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






    

In [ ]: