61 lines
2.0 KiB
Matlab
61 lines
2.0 KiB
Matlab
clear all;
|
|
%physical constants in MKS units
|
|
|
|
hbar = 1.054e-34;
|
|
q = 1.602e-19;
|
|
m = 9.110e-31;
|
|
|
|
%generate lattice
|
|
|
|
N = 100; %number of lattice points
|
|
n = [1:N]; %lattice points
|
|
a = 1e-10; %lattice constant
|
|
x = a * n; %x-coordinates
|
|
t0 = (hbar^2)/(2*m*a^2)/q; %encapsulating factor
|
|
L = a * (N+1); %total length of consideration
|
|
|
|
%set up Hamiltonian matrix
|
|
|
|
U = 0*x; %0 potential at all x
|
|
main_diag = diag(2*t0*ones(1,N)+U,0); %create main diagonal matrix
|
|
lower_diag = diag(-t0*ones(1,N-1),-1); %create lower diagonal matrix
|
|
upper_diag = diag(-t0*ones(1,N-1),+1); %create upper diagonal matrix
|
|
|
|
H = main_diag + lower_diag + upper_diag; %sum to get Hamiltonian matrix
|
|
|
|
% Modify hamiltonian for circular boundary conditions
|
|
H(1, N) = -t0;
|
|
H(N, 1) = -t0;
|
|
|
|
[eigenvectors,E_diag] = eig(H); %"eigenvectors" is a matrix wherein each column is an eigenvector
|
|
%"E_diag" is a diagonal matrix where the
|
|
%corresponding eigenvalues are on the
|
|
%diagonal.
|
|
|
|
E_col = diag(E_diag); %folds E_diag into a column vector of eigenvalues
|
|
|
|
% return eigenvectors for the 1st and 50th eigenvalues
|
|
|
|
phi_4 = eigenvectors(:,4);
|
|
phi_5 = eigenvectors(:,5);
|
|
|
|
% find the probability densities of position for 1st and 50th eigenvectors
|
|
|
|
P_4 = phi_4 .* conj(phi_4);
|
|
P_5 = phi_5 .* conj(phi_5);
|
|
|
|
% Find first N analytic eigenvalues
|
|
E_col_analytic = (1/q) * (hbar^2 * pi^2 * n.*n) / (2*m*L^2);
|
|
|
|
% Plot the probability densities for 1st and 50th eigenvectors
|
|
|
|
figure(1); clf; h = plot(x,P_4,'kx',x,P_5,'k-');
|
|
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
|
|
xlabel('POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
|
|
legend('n=4','n=5');
|
|
|
|
% Plot numerical eigenvalues
|
|
figure(2); clf; h = plot(n,E_col,'kx'); grid on;
|
|
set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
|
|
xlabel('EIGENVALUE NUMBER'); ylabel('ENERGY [eV]');
|