Trạng thái dừng của hàm sóng trong hố thế parabol

  • Code chương trình

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
function standing_wave_parabol_potential_well
% Created by Tran Hai Cat
% 2018.05.15
clc;
clear variables
close all

global h m

%% CONSTANTS
h = 1.054e-34;
m = 9.1095e-31;
q_element = 1.6e-19;
Angstrom = 1e-10;
k = 2*m/h/h;

%% INPUT PARAMETRS
n = 2; % n<6

E_standing_wave = [-73.0982,-59.295,-45.49,-31.68,-17.86]*q_element;
U0 = -80*q_element;

Emin = -100*q_element;
Emax = 100*q_element;

xmin = -6*Angstrom;
xmax = 6*Angstrom;
xmin_cal = -6*Angstrom;
xmax_cal = 3*Angstrom;
dx = (xmax-xmin)/1000;

%% data processing
x = linspace(xmin,xmax,1000);
E = E_standing_wave(n);
U = U_parabol(U0,x);

omega = abs(E)/h;
T = 2*pi/omega;
dt = T/1000;
tmax = 100*T;
t = 0;
E_draw = E/q_element;

xcal = xmin_cal;
psi = 1;

%% FIGURE
fig_wave = figure('name','Standing_wave_potential_well','color','k','numbertitle','off');
hold on
line_psi = plot(xcal/Angstrom,E_draw+psi,'linewidth',1,'color','y');
plot(x/Angstrom,U/q_element,'linewidth',2,'color','w');
plot(x/Angstrom,Energy(E,x)/q_element,'b','linewidth',2);
ht = title(sprintf('E = %0.3f [eV]',E/q_element));
axis([xmin/Angstrom xmax/Angstrom Emin/q_element Emax/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')

%% CALCULATION
psi_array = psi;
U = U_parabol(U0,xcal);
psi1 = psi*sqrt(U-E);
psi_max = psi;
x_array = xcal;

psi_old2 = psi;
psi_old1 = psi;
dautien = 1;
while xcal<xmax_cal
xcal = xcal+dx;
x_array = [x_array xcal];
U = U_parabol(U0,xcal);
psi2 = -k*(E-U)*psi;
psi1 = psi1+psi2*dx;
psi = psi+psi1*dx;
psi_array = [psi_array psi];

if (dautien==1)&&(psi_old2<psi_old1)&&(psi<psi_old1)
psi_max = 2.5*abs(psi_old1);
psi_min = - psi_max;
dautien = 0;
end
psi_old2 = psi_old1;
psi_old1 = psi;

end

S = trapz(x_array,psi_array.^2);
% psi_array = psi_array/sqrt(S);

psi_psi = (psi_array/psi_max*Emax/q_element).^2*2e-2;
psi_draw = psi_array/psi_max*Emax/q_element;

figure(fig_wave)
plot(x_array/Angstrom,E_draw+psi_psi);
set(line_psi,'xdata',x_array/Angstrom,'ydata',E_draw+psi_draw);

j = sqrt(-1);
while t<tmax
t = t+dt;
Psi_time = psi_draw*exp(-j*omega*t);
set(line_psi,'ydata',E_draw+real(Psi_time));
pause(0.001);
end

function y = U_parabol(U0,x)
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
y(i) = U0+2e2*x(i)*x(i);
end

function y = Energy(E,x)
y = E*ones(size(x));