In [1]:

    

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


    

In [107]:

    

# p=diag(-ones(N-1), -1) + diag(ones(N-1), +1)

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

H = - 1 / 12 *(np.diag( np.ones(N - 2), -2) + np.diag( np.ones(N - 2), 2)) \
    + 4 / 3 * (np.diag(np.ones(N - 1), -1) + np.diag(np.ones(N - 1), 1)) \
    - 5 / 2 * np.diag(np.ones(N))

H[0, -1] = H[-1, 0] = 4 / 3
H[0, -2] = H[-2, 0] = H[1, -1] = H[-1, 1] = -1 / 12
H /= -dx**2
dp = 2 * np.pi / N / dx
p = np.arange(N) * 2 * np.pi / N / dx

#pd.DataFrame(H)


    

In [108]:

    

E, psi = LA.eigh(H)


    

In [109]:

    

plt.plot(p, E)


    

    
Out[109]:





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






    

In [110]:

    

H1=np.diag(np.ones(N-1), -1) - 2*np.diag(np.ones(N)) + np.diag(np.ones(N-1), 1)
H1[0, N-1]=H1[N-1, 0]=1
H1*=-1/(2*dx**2)


    

In [15]:

    

E1, psi1 = LA.eigh(H1)


    

In [16]:

    

plt.plot(E)
plt.plot(E1)


    

    
Out[16]:





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






    

In [8]:

    

for i in range(5):
    plt.plot((psi[:, i ]/np.sqrt(dx))**2)


    

In [9]:

    

np.exp(np.arange(3)).reshape((1, 3)) ** np.arange(3).reshape((3, 1))


    

    
Out[9]:





array([[  1.        ,   1.        ,   1.        ],
       [  1.        ,   2.71828183,   7.3890561 ],
       [  1.        ,   7.3890561 ,  54.59815003]])

In [20]:

    

E_ana = np.sum(H[:, 0].reshape(((1, N))) * np.exp(((2 * np.pi * 1j / N)) * np.arange(N).reshape((N, 1)))** (N - np.arange(N)). reshape((1, N)), axis=1)
E_ana1 = np.sum(H1[:, 0].reshape(((1, N))) * np.exp(((2 * np.pi * 1j / N)) * np.arange(N).reshape((N, 1)))** (N - np.arange(N)). reshape((1, N)), axis=1)


    

In [47]:

    

plt.plot(p / (N/2 * dp), E_ana.flatten().real/( dp * N /2)**2, '--')
plt.plot(p / (N/2 * dp), E_ana1.flatten().real/( dp * N /2)**2, '--')
plt.plot(p / (N/2 * dp), E/( dp * N /2)**2)
plt.plot(p / (N/2 * dp), E1/( dp * N /2)**2)
plt.plot (p/( N/2*dp ), 0.5*p**2/(dp*N/2)**2, linestyle='dotted')
plt.xlabel(r'$p/(\Delta p \, N /2) $')
plt.ylabel(r'$E/((\Delta p \, N /2)^2/ m_0)$')
plt.ylim([0, .7])


    

    
Out[47]:





(0, 0.7)






    

In [43]:

    

#0.5*p**2/(dp*N/2)**2


    

    
Out[43]:





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






    

In [117]:

    

E.shape


    

    
Out[117]:





(512,)

In [118]:

    

psi.shape


    

    
Out[118]:





(512, 512)

In [113]:

    




    

In [130]:

    




    

In [133]:

    

m_0 = 1.
hbar = 1.
omega = 1.
alpha = 0
coef = -(2 * m_0) / hbar ** 2

def V(x):
    return 0.5 * omega**2  * x ** 2 + alpha * x ** 4

def df(t, X):
    psi, psi_p = X
    return np.array([psi_p, (-coef * V(t) + coef * E) * psi])


    

In [186]:

    

X0 = np.array([0, 1.])
df(0, X0)


    

    
Out[186]:





array([ 0., -4.])

In [187]:

    

r = integrate.ode(df).set_integrator('zvode', method='bdf')
r.set_initial_value(X0, 0.)
dt = .001
E = .5 + 1
N = 5000
x = np.zeros((N, 2), dtype='complex')
i=1
t = 0
while r.successful() and i < N:
    x[i, :] = r.integrate(r.t+dt)
    i+=1
plt.plot(np.arange(N) * dt, (x[:, 0]).real)


    

    
Out[187]:





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






    

In [189]:

    

M = 11
y = np.zeros(M)
X0 = np.array([0, 1.])
for j in np.arange(M):
    E = j * 0.1 + 1
    r = integrate.ode(df).set_integrator('zvode', method='bdf')
    r.set_initial_value(X0, 0.)
    dt = .001
    N = 5000
    x = np.zeros((N, 2), dtype='complex')
    i=1
    t = 0
    while r.successful() and i < N:
        x[i, :] = r.integrate(r.t+dt)
        i+=1
    #print (E)
    y[j] = x[-1, 0]
    print(x[-1, 0])
    #plt.plot(np.arange(N) * dt, (x[:, 1]).real)
plt.plot(np.linspace(0, 1, 11) + 1, y)


    

    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:17: ComplexWarning: Casting complex values to real discards the imaginary part






    




(6078.43408422+0j)
(4105.23601103+0j)
(2591.67085899+0j)
(1449.81831301+0j)
(606.250228922+0j)
(-0.00944749706659+0j)
(-419.393426388+0j)
(-693.369954945+0j)
(-855.831823149+0j)
(-934.287726521+0j)
(-950.983353841+0j)






    
Out[189]:





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






    

In [170]:

    

plt.plot(y)


    

    
Out[170]:





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






    

In [ ]: