In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as LA
from scipy.linalg import dft
import pandas as pd
from scipy import integrate

%matplotlib inline


    

In [31]:

    

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

N = 2 ** 9            # 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 [27]:

    

E, psi = LA.eigh(H)


    

In [28]:

    

plt.plot(p, E)


    

    
Out[28]:





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






    

In [14]:

    

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 [36]:

    




    

In [37]:

    




    

In [38]:

    




    

In [39]:

    

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


    

    
Out[39]:





array([ 1.,  0.])

In [40]:

    

r = integrate.ode(df).set_integrator('zvode', method='bdf')
r.set_initial_value(X0, 0.)
dt = .01
E = 0.5
N = 150
x = np.zeros((N, 2))
i=1
print(r.integrate(r.t+dt))
while r.successful() and i < N:
    x[i, :] = r.integrate(r.t+dt)
    i+=1
plt.plot(np.absolute(x[:, 0]))


    

    




[ 1.01015427+0.j  0.00000000+0.j]






    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:10: ComplexWarning: Casting complex values to real discards the imaginary part
  # Remove the CWD from sys.path while we load stuff.
/home/lorenzo/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ode.py:869: UserWarning: zvode: Excess work done on this call. (Perhaps wrong MF.)
  'Unexpected istate=%s' % istate))






    
Out[40]:





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






    

In [ ]:

    

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

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

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

# initial condition
X0 = [20.,10.]
parameters = (.2, .01, .001, .1)
t = np.linspace(0,300,1000)
X, infodict = integrate.odeint(dX_dt, X0, t, full_output=True)
infodict['message']