基于matlab的薛定谔方程数值计算应用（源码）.zip

%**************************************************************** % Program 4:本程序完全从数值的角度出发,把求解薛定谔方程的过程看成是一个 % 求解哈密顿数值矩阵的特征值以及对应特征向量的一个过程 %**************************************************************** % 确定参数|求解薛定谔方程中使用到的一些参数: hbar=1;m=1; %在原子单位体系里(hbar= 1,质量m= 1) L = 5; % 区间长度为2L:-L < x < L N = 1000;x = linspace(-L,L,N)'; % 总共在区间内均匀地取N个点 dx = x(2) - x(1); % 步长大小 % 选定U-x关系|在这里可以任意以多项式的形式进行组合 (这里首先假定的是一维谐振子的情况) U = 1/2*100*x.^2; %一维谐振子满足的U-x关系 %U = 1/2*100*x.^2+1/2*x.^(4); %经过多项式组合的U-x关系(可以把上面U-x注释掉，替换成这个U-x关系) %读者也可以自己以多项式的形式组合 %写出哈密顿的数值表达式 e = ones(N,1); Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2; %二阶导数的数值表达形式 H = -1/2*(hbar^2/m)*Lap + spdiags(U.*e,0,N,N); % 完整哈密顿矩阵 %计算能量本征值与能量本征函数 nmodes = 3; options.disp = 0; [V,E] = eigs(H,nmodes,'sa',options); %利用求解特征值与特征向量的方式求出本征值与本征函数 [E,ind] = sort(diag(E));% 将E按数值大小从低到高排序 V = V(:,ind); %取出前三级的本证函数 Usc = U*max(abs(V(:)))/max(abs(U)); %为了便于观察，将U的范围按比例缩小U-x关系不变 %展示结果(如果读者认为展现的图片看不清其中的曲线，请调整显示的曲线数目) plot(x,V,x,Usc,'--k'); %展现前三级在一定的U-x关系下本征函数V-x的关系 lgnd_str = [repmat('E=',nmodes,1),num2str(E)]; %添加标注 legend(lgnd_str) shg %显示图形窗口 %%读者可以取消下面几行程序的注释,验证求得的能量本征函数是否满足归一化 % SumV = 0; % for numV=1:1000 % SumV = SumV+V(numV,1)^2; % end % SumV