import numpy as npy
import scipy.linalg as spyl
|
|
a = 1e-10
hbar = 1.0546e-34
q = 1.602e-19
m = 9.11e-31
|
|
N=1000 # (Roughly) number of steps
delta_x = a/N
|
|
# First try banded version!
|
|
A = npy.zeros((2,N-2))
|
|
A[0,:]=2
A[1,:]=-1
|
|
eigs, vectors = spyl.eig_banded(A, lower=True, eigvals_only=False, select='i', select_range=(0,10)) # Get the lowest 10 eigenvalues
|
|
Energies = eigs/(delta_x**2)*hbar**2/(2*m)
|
|
Energies
array([ 6.03664410e-18, 2.41465167e-17, 5.43294387e-17,
9.65851117e-17, 1.50913118e-16, 2.17312919e-16,
2.95783860e-16, 3.86325164e-16, 4.88935936e-16,
6.03615161e-16, 7.30361704e-16])
array([ 6.03664410e-18, 2.41465167e-17, 5.43294387e-17,
9.65851117e-17, 1.50913118e-16, 2.17312919e-16,
2.95783860e-16, 3.86325164e-16, 4.88935936e-16,
6.03615161e-16, 7.30361704e-16])
|
# Compare with theory
|
|
E_theory = ((npy.arange(1,12))**2 * hbar**2 * 3.1416**2)/(2*m*a**2)
|
|
E_theory
array([ 6.02460999e-18, 2.40984400e-17, 5.42214899e-17,
9.63937599e-17, 1.50615250e-16, 2.16885960e-16,
2.95205890e-16, 3.85575040e-16, 4.87993409e-16,
6.02460999e-16, 7.28977809e-16])
array([ 6.02460999e-18, 2.40984400e-17, 5.42214899e-17,
9.63937599e-17, 1.50615250e-16, 2.16885960e-16,
2.95205890e-16, 3.85575040e-16, 4.87993409e-16,
6.02460999e-16, 7.28977809e-16])
|
(Energies - E_theory)/E_theory*100
array([ 0.1997492 , 0.19950147, 0.1990886 , 0.19851057, 0.19776739,
0.19685907, 0.19578561, 0.19454702, 0.19314329, 0.19157444,
0.18984047])
array([ 0.1997492 , 0.19950147, 0.1990886 , 0.19851057, 0.19776739,
0.19685907, 0.19578561, 0.19454702, 0.19314329, 0.19157444,
0.18984047])
|
x = npy.arange(N-2)*delta_x
plots = [line2d(npy.column_stack((x, vectors[:,index]))) for index in range(3)]
|
|
"""I could have used "zip" in the above line for a neater syntax, but column_stack is probably more efficient."""
'I could have used "zip" in the above line for a neater syntax, but column_stack is probably more efficient.' 'I could have used "zip" in the above line for a neater syntax, but column_stack is probably more efficient.' |
all_plots = sum(plots)
plot(all_plots)
|
|
