-
Notifications
You must be signed in to change notification settings - Fork 11
/
Schroed2D_FEM_f.m
123 lines (103 loc) · 7.04 KB
/
Schroed2D_FEM_f.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
function[E,psi]=Schroed2D_FEM_f(x,y,V0,Mass,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
e=1.602176487E-19; %% electron charge [C]
m0=9.10938188E-31; %% electron mass [kg]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nx=length(x);
Ny=length(y);
dx=x(2)-x(1);
dy=y(2)-y(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Building of the operators %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second derivative X %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DX2(Ny=5,Nx=4); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% -2 0 0 0 0 | 1 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 -2 0 0 0 | 0 1 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 -2 0 0 | 0 0 1 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 -2 0 | 0 0 0 1 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 0 -2 | 0 0 0 0 1 | 0 0 0 0 0 | 0 0 0 0 0 %
% ----------------------------------------------------------------------------- %
% 1 0 0 0 0 |-2 0 0 0 0 | 1 0 0 0 0 | 0 0 0 0 0 %
% 0 1 0 0 0 | 0 -2 0 0 0 | 0 1 0 0 0 | 0 0 0 0 0 %
% 0 0 1 0 0 | 0 0 -2 0 0 | 0 0 1 0 0 | 0 0 0 0 0 %
% 0 0 0 1 0 | 0 0 0 -2 0 | 0 0 0 1 0 | 0 0 0 0 0 %
% 0 0 0 0 1 | 0 0 0 0 -2 | 0 0 0 0 1 | 0 0 0 0 0 %
% ----------------------------------------------------------------------------- %
% 0 0 0 0 0 | 1 0 0 0 0 |-2 0 0 0 0 | 1 0 0 0 0 %
% 0 0 0 0 0 | 0 1 0 0 0 | 0 -2 0 0 0 | 0 1 0 0 0 %
% 0 0 0 0 0 | 0 0 1 0 0 | 0 0 -2 0 0 | 0 0 1 0 0 %
% 0 0 0 0 0 | 0 0 0 1 0 | 0 0 0 -2 0 | 0 0 0 1 0 %
% 0 0 0 0 0 | 0 0 0 0 1 | 0 0 0 0 -2 | 0 0 0 0 1 %
% ----------------------------------------------------------------------------- %
% 0 0 0 0 0 | 0 0 0 0 0 | 1 0 0 0 0 |-2 0 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 1 0 0 0 | 0 -2 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 1 0 0 | 0 0 -2 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 1 0 | 0 0 0 -2 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 1 | 0 0 0 0 -2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Axy = ones(1,(Nx-1)*Ny);
DX2 = (-2)*diag(ones(1,Ny*Nx)) + (1)*diag(Axy,-Ny) + (1)*diag(Axy,Ny);
%%%%%%%%%%%%%%%%%%%%%%%%%%% Second derivative Y %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DY2(Ny=5,Nx=4); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% -2 1 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 1 -2 1 0 0 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 1 -2 1 0 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 1 -2 1 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 1 -2 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% ----------------------------------------------------------------------------- %
% 0 0 0 0 0 |-2 1 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 1 -2 1 0 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 1 -2 1 0 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 0 1 -2 1 | 0 0 0 0 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 0 0 1 -2 | 0 0 0 0 0 | 0 0 0 0 0 %
% ----------------------------------------------------------------------------- %
% 0 0 0 0 0 | 0 0 0 0 0 |-2 1 0 0 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 1 -2 1 0 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 1 -2 1 0 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 1 -2 1 | 0 0 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 1 -2 | 0 0 0 0 0 %
% ----------------------------------------------------------------------------- %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 |-2 1 0 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 | 1 -2 1 0 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 | 0 1 -2 1 0 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 1 -2 1 %
% 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 | 0 0 0 1 -2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AA=ones(1,Nx*Ny);
BB=ones(1,Nx*Ny-1);
BB(Ny:Ny:end)=0;
DY2=(-2)*diag(AA) + (1)*diag(BB,-1) + (1)*diag(BB,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Building of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=(-hbar^2/(2*m0*Mass)) * ( DX2/dx^2 + DY2/dy^2 ) + diag(V0(:)*e) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Diagonalisation of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=sparse(H);
[PSI,Energy] = eigs(H,n,'SM');
E = diag(Energy)/e;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Normalization of the Wavefunctions %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n
psi_temp=reshape(PSI(:,i),Ny,Nx);
psi(:,:,i) = psi_temp / sqrt( trapz( y' , trapz(x,abs(psi_temp).^2 ,2) , 1 ) ); % normalisation of the wave function psi
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here is a small patch due to differences between Octave and Matlab
% Matlab order the eigen values while Octave reverse it
if E(1)>E(2)
psi=psi(:,:,end:-1:1);
E=E(end:-1:1);
end
end