b=4; N=5000; M=N+1; Dx=(b-a)/N; for m=1:M
x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;
F(m,m)=2/((Dx)^2)+x(m)^2; end
for m=1:M-1;
F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end [v,d]=eig(F); for n=1:1 E5(n)=d(n,n)/2
A=sqrt(trapz(x,v(:,n).*v(:,n))); hold on
plot(x,v(:,n)/A,'-k', 'LineWidth',2) end clear a=-4; b=4; N=6000; M=N+1; Dx=(b-a)/N; for m=1:M
x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;
F(m,m)=2/((Dx)^2)+x(m)^2; end
for m=1:M-1;
F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2);
图7计算区间取为??4,4?, 计算点数分别为1000、2000、3000、4000、5000和6000.所有的波函数曲线重合为一条曲线.
17
end [v,d]=eig(F); for n=1:1 E6(n)=d(n,n)/2
A=sqrt(trapz(x,v(:,n).*v(:,n))); hold on
plot(x,v(:,n)/A,'-k', 'LineWidth',2) end
计算得到的本征值为:
E1 =0.5000, E2 =0.5000, E3 =0.5000, E4 =0.5000, E5 =0.5000, E6 =0.5000 计算所得的波函数如图7所示, 所有波函数的曲线均重合为一条曲线.
下面我们设计的程序计算区间分别取??4,4?,??5,5?,??6,6?,??7,7?,??8,8?和??9,9?, 计算点分别取2000、2500、3000、3500、4000和4500, 这样可保证计算精度相同. 计算程序:
clear a=-4; b=4; N=2000; M=N+1; Dx=(b-a)/N; for m=1:M
x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;
F(m,m)=2/((Dx)^2)+x(m)^2; end
for m=1:M-1;
F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end [v,d]=eig(F); for n=1:1 E1(n)=d(n,n)/2 figure
A=sqrt(trapz(x,v(:,n).*v(:,n))); plot(x,v(:,n)/A,'-k', 'LineWidth',2) axis([-9 9 0 0.8]) set(gca,'XTick',-9:1:9)
18
set(gca,'YTick',0:0.1:0.8)
t=['\\fontsize{14}\\rm\\Psi','_{',int2str(n-1),'}','(\\xi)']; grid on
xlabel ('\\fontsize{14}\\rm\\xi', 'Color','k') ylabel (t, 'Color','k')
set(gca,'ygrid','on','GridLineStyle','-') set(gca,'xgrid','on','GridLineStyle','-') end clear a=-5; b=5; N=2500; M=N+1; Dx=(b-a)/N; for m=1:M
x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;
F(m,m)=2/((Dx)^2)+x(m)^2; end
for m=1:M-1;
F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end [v,d]=eig(F); for n=1:1 E2(n)=d(n,n)/2
A=sqrt(trapz(x,v(:,n).*v(:,n))); hold on
plot(x,v(:,n)/A,'-k', 'LineWidth',2) end clear a=-6; b=6; N=3000; M=N+1; Dx=(b-a)/N;
19
for m=1:M
x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;
F(m,m)=2/((Dx)^2)+x(m)^2; end
for m=1:M-1;
F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end [v,d]=eig(F); for n=1:1 E3(n)=d(n,n)/2
A=sqrt(trapz(x,v(:,n).*v(:,n))); hold on
plot(x,v(:,n)/A,'-k', 'LineWidth',2) end clear a=-7; b=7; N=3500; M=N+1; Dx=(b-a)/N; for m=1:M
x(m)=a+(m-1)*Dx; end F=zeros(M); for m=1:M;
F(m,m)=2/((Dx)^2)+x(m)^2; end
for m=1:M-1;
F(m+1,m)=-1/((Dx)^2); F(m,m+1)=-1/((Dx)^2); end [v,d]=eig(F); for n=1:1 E4(n)=d(n,n)/2
20
