# uncomment ONE line to choose matplotlib backend
# if using Jupyter Notebook, use interactive "notebook" backend for best results
# if using Jupyter Lab, use interactive "widget" backend for best results
# if both fail, use static "inline" backend
%matplotlib notebook
#%matplotlib widget
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as linalg
h_bar = 1.05e-34
m = 9.11e-31 #elektronmasse [kg]
start = -1
stop = 1
L = stop - start
N = 1000
x, dx = np.linspace(start, stop, N+1, retstep = True)
n = [1, 2, 3]
'''FUNKSJON SOM RETURNERER EGENVERDIER E OG TILHØRENDE, NORMERTE EGENVERDIFUNKSJONER PSI'''
def calc_E_and_Psi(V):
N_tot = len(V)
d = [v + h_bar**2/(m*dx**2) for v in V]
e = [-h_bar**2/(2*m*dx**2) for n in range(N_tot-1)]
E,Psi = linalg.eigh_tridiagonal(d,e)
return E, Psi
'''FUNKSJON SOM FRAMSTILLER POTENSIALET V, ENERGIEGENVERDIER E OG ENERGIEGENFUNKSJONER PSI'''
def plot_psi(analytic, numeric, analytic_2, numeric_2, x_list, start, stop, xlim):
'''
Denne funksjonen tar inn liste av lister med analytiske og numeriske verdier.
Hvert element i listene er en liste av psi for en bestemt tilstand nr. n.
'''
fig = plt.figure(figsize=(15,10))
for n in range(len(analytic)):
ax1 = plt.subplot("221")
ax1.plot(x_list, analytic[n][0], label="Ψ "+str(n+1))
ax1.set_xlabel("L", fontsize=13)
ax1.set_ylabel('Ψ', fontsize=13)
ax1.set_title("Analytisk", fontsize=20)
ax1.grid()
ax1.set_xlim(xlim)
plt.legend()
ax2 = plt.subplot("222")
ax2.plot(x_list, numeric[n][0], label="Ψ" +str(n+1))
ax2.set_xlabel("L", fontsize=13)
ax2.set_ylabel('Ψ', fontsize=13)
ax2.set_title("Numerisk", fontsize=20)
ax2.grid()
ax2.set_xlim(xlim)
plt.legend()
ax3 = plt.subplot("223")
ax3.plot(x_list, analytic_2[n][0], label="Ψ^2 "+str(n+1))
ax3.set_xlabel("L", fontsize=13)
ax3.set_ylabel('Ψ^2', fontsize=13)
ax3.set_title("Analytisk", fontsize=20)
ax3.grid()
ax3.set_xlim(xlim)
plt.legend()
ax4 = plt.subplot("224")
ax4.plot(x_list, numeric_2[n][0], label="Ψ^2 "+str(n+1))
ax4.set_xlabel("L", fontsize=13)
ax4.set_ylabel('Ψ^2', fontsize=13)
ax4.set_title("Numerisk", fontsize=20)
ax4.grid()
ax4.set_xlim(xlim)
plt.legend()
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
def plot_V_E(V_list, E_list_analytic, E_list_numeric, x_list, start, stop):
'''
Denne funksjonen tar inn liste for potensialet og energiegenverdiene
'''
#potensial
plt.figure(1,figsize=(10,5))
plt.plot(x_list, V_list, label = 'Potential' )
plt.xlabel("L", fontsize=13)
plt.ylabel('V', fontsize=13)
plt.xticks(np.arange(start, stop, step=stop/10))
plt.title("Potensial", fontsize=16)
plt.grid()
# analytic solutions
plt.figure(2, figsize=(10,5))
for i in range(len(n)):
plt.plot(x_list, E_list_analytic[i][0], label = 'Eigenvalue ' + str(i+1))
plt.xlabel("L", fontsize=13)
plt.ylabel('V, E', fontsize=13)
plt.xticks(np.arange(start, stop, step=stop/10))
plt.title("Analytisk løsning", fontsize=16)
plt.legend()
plt.grid()
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
#numeric solutions
plt.figure(3, figsize=(10,5))
for i in range(len(n)):
plt.plot(x_list, E_list_numeric[i][0], label = 'Eigenvalue ' + str(i+1))
plt.xlabel("L", fontsize=13)
plt.ylabel('V, E', fontsize=13)
plt.xticks(np.arange(start, stop, step=stop/10))
plt.title("Numerisk løsning", fontsize=16)
plt.legend()
plt.grid()
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
'''POTENSIALET'''
V_infwell = [0]*len(x)
'''ANALYTISK LØSNING'''
def psi_analytic(x,n):
return np.sqrt(2/L)*np.sin(n*np.pi*(x+1)/L)
def E_analytic(x_vals,n):
E_an = np.zeros(len(x_vals))
for i in range(len(x_vals)):
E_an[i] = n**2*np.pi**2*h_bar**2/(2*m*L**2)
return E_an
'''NUMERISK LØSNING'''
E_num_infwell_calc, Psi_num_infwell_calc = calc_E_and_Psi(V_infwell)
'''LISTENE SOM SKAL PLOTTES'''
E_an_infwell = [[] for _ in range(len(n))]
Psi_an_infwell = [[] for _ in range(len(n))]
E_num_infwell = [[] for _ in range(len(n))]
Psi_num_infwell = [[] for _ in range(len(n))]
for i in range(len(n)):
E_an_infwell[i].append(E_analytic(x, n[i]))
Psi_an_infwell[i].append(psi_analytic(x, n[i]))
Psi_num_infwell[i].append(np.array([0] + Psi_num_infwell_calc[:, n[i]-1] + [0]))
E_num_infwell[i] = np.full((1, len(x)), E_num_infwell_calc[n[i]-1])
xlim_infwell = (-1,1)
'''PLOTTER RESULTATENE'''
print(plot_V_E(V_infwell, np.array(E_an_infwell), np.array(E_num_infwell), x, start, stop))
print(plot_psi(np.array(Psi_an_infwell), np.array(Psi_num_infwell), np.abs(np.array(Psi_an_infwell))**2, np.abs(np.array(Psi_num_infwell))**2, x, start, stop, xlim_infwell))
from scipy import special
''' gammel kode '''
#w = 1
#V0=1.6E-19
#N = 100
#start_h = -0.1
#stop_h = 0.1
''' ny kode '''
w = 1
N = 1000
start_h = -1
stop_h = 1
'''Diskretiserer potensialet for harmonisk ocsillator, gammel kode'''
#V_harmosc = [V0]*4*N + [V0*((n-N)/(N*1.0))**2 for n in range(2*N+1)] + [V0]*4*N
#x_harmosc = np.linspace(start_h, stop_h, len(V_harmosc))
''' Diskretiserer potensialet, ny kode '''
# definerer x-akse for hele plottet (konstant + varierende potensial)
x_harmosc, dx = np.linspace(start_h, stop_h, N*3 , retstep = True)
''' Potensialet der det ikke er konstant'''
def V_middle_calc(m, w, x_list):
V = [0]*len(x_list)
for i in range(len(x_list)):
V[i] = 1/2 * m * w**2 * x_list[i]**2
return V
# potensialet i midten skal variere etter HO formel
V_middle = V_middle_calc(m, w, x_harmosc[N: 2* N])
# konstant potensial før og etter det varierende
V_const = [V_middle[0]]*N
# Samler de til en array
V_harmosc = np.array(V_const + V_middle + V_const)
'''Analytisk løsning'''
def psi_harmosc_an(x,n):
'''x: en array, n: et heltall. returnerer en array med psi-verdier for denne n-verdien'''
psi = np.array(1/(np.sqrt((2**n)*np.math.factorial(n)))*(m*w/(np.pi*h_bar))**(1/4)*np.exp(-m*w*x**2/(2*h_bar))*special.hermite(n)((x*(m*w/h_bar)**(1/2))))
return psi
def E_harmosc_an(x, n,w):
E_an = np.zeros(len(x))
for i in range(len(x)):
E_an[i] = np.array((n+1/2)*h_bar*w)
return E_an
'''Numerisk løsning'''
d = [v + h_bar**2/(m*dx**2) for v in V_harmosc]
e = [-h_bar**2/(2*m*dx**2) for n in range(len(V_harmosc)-1)]
E,Psi = linalg.eigh_tridiagonal(d,e)
'''LISTENE SOM SKAL PLOTTES'''
E_an_harmosc = [[] for _ in range(len(n))]
Psi_an_harmosc = [[] for _ in range(len(n))]
E_num_harmosc = [[] for _ in range(len(n))]
Psi_num_harmosc = [[] for _ in range(len(n))]
for i in range(len(n)):
E_an_harmosc[i].append(E_harmosc_an(x_harmosc, n[i] - 1, w))
Psi_an_harmosc[i].append(psi_harmosc_an(x_harmosc, n[i] - 1))
E_num_harmosc[i] = np.full((1, len(x_harmosc)), E[n[i] - 1])
Psi_num_harmosc[i].append(Psi[:, n[i] - 1])
xlim_harmosc = (-0.1, 0.1)
'''PLOTTER RESULTATENE'''
print(plot_V_E(V_harmosc, np.array(E_an_harmosc), np.array(E_num_harmosc), x_harmosc, start_h, stop_h))
print(plot_psi(np.array(Psi_an_harmosc), np.array(Psi_num_harmosc), np.abs(np.array(Psi_an_harmosc))**2, np.abs(np.array(Psi_num_harmosc))**2, x_harmosc, start_h, stop_h, xlim_harmosc))
V0=4.0*1.6E-19
dx=0.1E-10
N = 100
'''Diskretiserer potensialet for endelig brønn'''
V = np.asarray([V0]*10*N + [0]*N + [V0]*10*N)
x_finwell = np.linspace(start_h, stop_h, len(V))*1.0e9
'''Numerisk løsning'''
d = [v + h_bar**2/(m*dx**2) for v in V]
e = [-h_bar**2/(2*m*dx**2) for n in range(len(V)-1)]
E,Psi = linalg.eigh_tridiagonal(d,e)
'''LISTENE SOM SKAL PLOTTES'''
E_num_finwell = [[] for _ in range(len(n))]
Psi_num_finwell = [[] for _ in range(len(n))]
for i in range(len(n)):
E_num_finwell[i] = np.full((1, len(x_finwell)), E[n[i]])
Psi_num_finwell[i].append(Psi[:, n[i]])
'''PLOTTING'''
plt.figure(1, (15,5))
for i in range(len(n)):
plt.plot(x_finwell, E_num_finwell[i][0])
plt.title('Energiverdier',fontsize=20)
plt.xlabel('$x$ (nm)',fontsize=20)
plt.ylabel('$E(x)$',fontsize=20)
plt.grid()
plt.show()
plt.figure(2,figsize=(15,5))
plt.plot(x_finwell, Psi[:,0],x_finwell,Psi[:,1],x_finwell,Psi[:,2],x_finwell,Psi[:,3])
plt.title('Bundne tilstander',fontsize=20)
plt.xlabel('$x$ (nm)',fontsize=20)
plt.ylabel('$\psi(x)$',fontsize=20)
plt.grid()
plt.show()
plt.figure(3,figsize=(15,5))
plt.plot(x_finwell, Psi[:,15],x_finwell,Psi[:,16],x_finwell,Psi[:,17],x_finwell,Psi[:,18])
plt.title('Ubundne tilstander',fontsize=20)
plt.xlabel('$x$ (nm)',fontsize=20)
plt.ylabel('$\psi(x)$',fontsize=20)
plt.grid()
plt.show()