【薛定谔方程求解】基于MATLAB的薛定谔方程求解仿真

1.软件版本

matlab2013b

2.本算法理论知识

3.部分源码

clc;
clear;
close all;
figure
p     = 3*pi/2;
a     = 1;

alpha = -15:0.0005:15;
index = 1;

for i=1:length(alpha)
    if abs(alpha(i))<0.001
    y(i) = p + 1;
    else
    y(i) = p*sin(a*alpha(i))/(a*alpha(i)) + cos(a*alpha(i));    
    end     
end

plot(alpha,y,'b');hold on;
plot(-15:0.01:15,1,'r-','LineWidth',4);hold on;
plot(-15:0.01:15,-1,'r-','LineWidth',4);hold on;
title('y = p*sin(a*alpha)/(a*alpha) + cos(a*alpha)');
figure_s();
hold off;
%******注意通过以下代码求出图像中的填充区域*******
%四个区间：
%A:2.2540 - 3.1480
%B:4.7206 - 6.2965
%C:7.4290 - 9.4450
%D:10.3040- 13.2755 
% startpoint =  2.2540;
% endpoint   =  3.1480;
% N          =  10000;
% k          =  0;
% alpha2 = solvealpha(startpoint,endpoint,N,k);
% alpha2
%******注意通过以上代码求出图像中的填充区域*******

legend('p*sin(a*alpha)/(a*alpha) + cos(a*alpha)','k');
text(pi,  -1,'k=pi');
text(2*pi, 1,'k=2pi');
text(3*pi,-1,'k=3pi');
text(4*pi, 1,'k=4pi');
grid on;


figure
%使用四个区间的平方，作为图2的y轴区间
%使用四个区间的平方，作为图2的y轴区间
%使用四个区间的平方，作为图2的y轴区间
%AAAA
startpoint =  2.2540;
endpoint   =  3.1480;
N          =  10000;
k          =  -pi:(2*pi)/20:pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-pi:(2*pi)/80:pi,alpha2(length(alpha2)));hold on;

clear k
clear alpha2

%BBBB
startpoint =  4.7206;
endpoint   =  6.2965;
N          =  10000;
k          =  -2*pi:(2*pi)/20:-1*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
k          =  pi:(2*pi)/20:2*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;

plot(-1*pi:(2*pi)/80:1*pi,alpha2(1));hold on;
plot(-2*pi:(2*pi)/80:2*pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2


%CCCC
startpoint =  7.4290;
endpoint   =  9.4450;
N          =  10000;
k          =  -3*pi:(2*pi)/20:-2*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
k          =  2*pi:(2*pi)/20:3*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-2*pi:(2*pi)/80:2*pi,alpha2(1));hold on;
plot(-3*pi:(2*pi)/80:3*pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2
%DDDD
startpoint =  10.3040;
endpoint   =  13.2755;
N          =  10000;
k          =  -4*pi:(2*pi)/20:-3*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
k          =  3*pi:(2*pi)/20:4*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-3*pi:(2*pi)/80:3*pi,alpha2(1));hold on;
plot(-4*pi:(2*pi)/80:4*pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2

%画虚线
plot(pi,  0:4.7206^2,'r');hold on;
plot(2*pi,0:7.4290^2,'r');hold on;
plot(3*pi,0:10.3040^2,'r');hold on;
plot(4*pi,0:13.2755^2,'r');hold on;

text(-1, 145,'允带');
text(-1, 96,'禁带');
text(-1, 75,'允带');
text(-1, 45,'禁带');
text(-1, 30,'允带');
text(-1, 15,'禁带');

4.仿真结论

5.参考文献

[1]王忆锋, 唐利斌. Kronig-Penney能带模型的MATLAB分析与计算[J]. 激光与红外, 2010, 40(7):5.A16-07