有限差分法求一维谐振子本征值
clear all
clc
%令h/(2pi)=w=1
L=100; %范围
a=0.05;%精度调节
Np=L/a;
X=[0:a:L-a];
t0=1/(2*a^2);
Vn=zeros(Np);
for i=1:Np
Vn(i,i)=0.5*(i*a-L/2)^2; %原点移到中心
end
T1=(2*t0*diag(ones(1,Np)))-(t0*diag(ones(1,Np-1),1))-(t0*diag(ones(1,Np-1),-1));
T=T1+Vn;
[V,D]=eig(T);%求特征值
D=diag(D);
[Enum,ind]=sort(D); %排序
E1=D(ind(1));
psi1=abs(V(:,ind(1)));
P1=psi1.*conj(psi1);
E2=D(ind(2));
psi2=abs(V(:,ind(2)));
P2=psi2.*conj(psi2);
E3=D(ind(3));
psi2=abs(V(:,ind(3)));
P3=psi2.*conj(psi2);
E4=D(ind(4));
psi2=abs(V(:,ind(4)));
P4=psi2.*conj(psi2);
E5=D(ind(5));
psi2=abs(V(:,ind(5)));
P5=psi2.*conj(psi2);
figure
Ean=[1:Np]+0.5; %本征值理论计算
Ean=Ean';
hold on
plot(Ean,Enum);% 对比
figure
hold on
h1=plot(X,P1);%
h2=plot(X,P2);%
h3=plot(X,P3);%
h4=plot(X,P4);%
h5=plot(X,P5);%
grid on