Dao động tử điều hoà – Phần 2

function wave_packet_parabol_potential_well
% Created by Tran Hai Cat
% 2019.04.15
clc;
clear variables
close all

global h m k k_parabol E

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

%% INPUT PARAMETRS FOR USER
omega0 = sqrt(k_parabol/m);
E0 = h*omega0*(7-1/2);

N = 30;
k_parabol = 400;
x0 = 0;
delta_x = 0.7*Angstrom;

Emin = -20*q_element;
Emax = 200*q_element;
xmin = -6*Angstrom;
xmax = 6*Angstrom;
x1 = -10*Angstrom;
x2 = 10*Angstrom;
x_inter = 0.5769*Angstrom;

%% CALCULATION
omega_array = zeros(1,N);
E_array = zeros(1,N);
psi = 1e-10;
psi1 = 0;
for i = 1:N
omega_classic = sqrt(k_parabol/m);
E_array(i) = h*omega_classic*(i-1/2);
omega_array(i) = E_array(i)/h;

E = E_array(i);
[x_left,psi_left] = ode45(@Schrodinger_equation,linspace(x1,x_inter,1500),[psi psi1]);
[x_right,psi_right] = ode45(@Schrodinger_equation,linspace(x2,x_inter,1500),[psi psi1]);
psi_right = sign(psi_left(end,2)*psi_right(end,2))*psi_right;
X = [x_left;flipud(x_right)];
Psi_E{i} = [psi_left(:,1);flipud(psi_right(:,1))];
S_norm = trapz(X,Psi_E{i}.*Psi_E{i});
Psi_E{i} = Psi_E{i}/sqrt(S_norm);
end

C = packet_analys(E0,x0,delta_x,N,X,Psi_E);

Psi_synthes = wavepacket_synthes(C,omega_array,Psi_E,0);
Psi_synthes_max = abs(max(Psi_synthes));
Psi_draw = 0.25*Psi_synthes/Psi_synthes_max*Emax/q_element;
psi_psi = 0.3*(conj(Psi_synthes).*Psi_synthes)...
/Psi_synthes_max/Psi_synthes_max*Emax/q_element;

Amp = sqrt(2*E0/k_parabol);
x_particle = 0;
y_particle = U_parabol(k_parabol,x_particle);

%% FIGURE
figure('name','Spectrum of Stationary States',...
'color','w','numbertitle','off');
subplot(2,1,1);
stem(abs(C),'fill','-','color','blue');
ylabel('Amplitude');

subplot(2,1,2);
stem(angle(C),'fill','-','color','blue');
xlabel('n (harmonic)');
ylabel('Phase');

x_plot = linspace(xmin,xmax,1000);
U_plot = U_parabol(k_parabol,x_plot);

figure('name','Harmonic oscillator',...
'color','black','numbertitle','off');
hold on
plot(x_plot/Angstrom,U_plot/q_element,'linewidth',2,'color','w');
line([xmin/Angstrom,xmax/Angstrom],[E0/q_element,E0/q_element],...
'linewidth',2,'color','b');
line_psi = plot(X/Angstrom,E0/q_element+real(Psi_draw),'linewidth',1,'color','y');
line_psi_psi = plot(X/Angstrom,E0/q_element+psi_psi,'linewidth',2);
draw_partile = plot(x_particle/Angstrom,y_particle/q_element,...
'ro','MarkerSize',10,'markerfacecolor','r');
title(sprintf('E0 = %0.3f [eV]',E0/q_element),'color','w');
axis([xmin/Angstrom xmax/Angstrom Emin/q_element Emax/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('x [Angstrom]');
ylabel('E, U [eV]');

t = 0;
dt = 1e-19;
while 1
t = t+dt;
Psi_synthes = wavepacket_synthes(C,omega_array,Psi_E,t);
Psi_draw = 0.25*Psi_synthes/Psi_synthes_max*Emax/q_element;
psi_psi = 0.3*(conj(Psi_synthes).*Psi_synthes)...
/Psi_synthes_max/Psi_synthes_max*Emax/q_element;
x_particle = Amp*sin(omega0*t);
y_particle = U_parabol(k_parabol,x_particle);

set(line_psi,'ydata',E0/q_element+real(Psi_draw));
set(line_psi_psi,'ydata',E0/q_element+real(psi_psi));
set(draw_partile,'xdata',x_particle/Angstrom);
set(draw_partile,'ydata',y_particle/q_element);
pause(0.001);
end

function dy = Schrodinger_equation(x,y)
%% Schrodinger Equation
global E k k_parabol
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k*(E-U_parabol(k_parabol,x))*y(1);

function y = U_parabol(k_parabol,x)
%% Potential Energy
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
y(i) = 0.5*k_parabol*x(i)*x(i);
end

function psi = wavepacket_Gauss(E,x0,delta_x,x)
%% wave packet with Gauss's distribution (normal distribution)
global h m
k0 = sqrt(2*m*E)/h;
psi = 300*exp(-(x-x0).^2/4/delta_x.^2).*exp(1i*k0*(x-x0));

function C = packet_analys(E,x0,delta_x,N,X,Psi_E)
%% Calculate koefficients of Stationary State
psi0 = wavepacket_Gauss(E,x0,delta_x,X);
C = zeros(1,N);
for i = 1:N
C(i) = trapz(X,conj(Psi_E{i}).*psi0);
end

function psi = wavepacket_synthes(C,omega_array,Psi_E,t)
%% Synthesys initial wave packet from koefficients of Stationary State
N = length(C);
psi = 0;
for i = 1:N
psi = psi+C(i)*Psi_E{i}*exp(-1i*omega_array(i)*t);
end