本文主要是介绍MATLAB定态氢原子波函数可视化，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
在知乎看到了
[王欢]的[用 Python+SciPy 可视化定态氢原子波函数（一）]
[王大可]的[利用Mathematica将定态氢原子波函数可视化]
手痒，想用MATLAB也实现一下，David Griffiths 所著的 Introduction to Quantum Mechanics 中:
ψ n , l , m ( r , θ , ϕ ) = R n l ( r ) Y l m ( θ , ϕ ) \psi_{n, l, m}(r,\theta, \phi)=R_{nl}(r)Y_{l}^{m}(\theta, \phi) ψn,l,m(r,θ,ϕ)=Rnl(r)Ylm(θ,ϕ)
即，函数由径向波函数及球形谐函数两部分组成，展开一下径向波函数是这样的：
ψ n , l , m ( r , θ , ϕ ) = ( 2 n a ) 3 ( n − l − 1 ) ! 2 n ( n + l ) ! e − r / n a ( 2 r n a ) l [ L n − l − 1 2 l + 1 ( 2 r n a ) ] Y l m ( θ , ϕ ) \psi_{n, l, m}(r,\theta, \phi)=\sqrt{\left(\frac{2}{n a}\right)^{3} \frac{(n-l-1) !}{2 n(n+l) !}} e^{-r / n a}\left(\frac{2 r}{n a}\right)^{l}\left[L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n a}\right)\right] Y_{l}^{m}(\theta, \phi) ψn,l,m(r,θ,ϕ)=(na2)32n(n+l)!(n−l−1)!e−r/na(na2r)l[Ln−l−12l+1(na2r)]Ylm(θ,ϕ)
其中 a a a 为玻尔半径，数值约为5.2917721092(17)×10^-11米，本文仅做可视化，方便起见，将 a a a 数值设置为1，并省略归一化系数，因此我们实际上本文绘制的是如下函数的图像：
ψ n , l , m ( r , θ , ϕ ) = e − r / n ( 2 r n ) l [ L n − l − 1 2 l + 1 ( 2 r n ) ] Y l m ( θ , ϕ ) \psi_{n, l, m}(r,\theta, \phi)=e^{-r / n}\left(\frac{2 r}{n}\right)^{l}\left[L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n}\right)\right] Y_{l}^{m}(\theta, \phi) ψn,l,m(r,θ,ϕ)=e−r/n(n2r)l[Ln−l−12l+1(n2r)]Ylm(θ,ϕ)
其中非常棒的是 l a g u e r r e laguerre laguerre 多项式 L n − l − 1 2 l + 1 ( 2 r n ) L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n}\right) Ln−l−12l+1(n2r)，MATLAB中有自带的函数 laguerreL可以使用，但是球谐函数(SphericalHarmonicY) Y l m ( θ , ϕ ) Y_{l}^{m}(\theta, \phi) Ylm(θ,ϕ)，并不是MATLAB的自带函数，不过好在，当 m ≥ 0 m\ge0 m≥0时球谐函数的展开式：
Y l m ( θ , ϕ ) = ( − 1 ) m ( 2 l + 1 ) 4 π ( l − ∣ m ∣ ) ! ( l + ∣ m ∣ ) ! e i m ϕ P l m ( cos θ ) Y_{l}^{m}(\theta, \phi)=(-1)^m \sqrt{\frac{(2 l+1)}{4 \pi} \frac{(l-|m|) !}{(l+|m|) !}} e^{i m \phi} P_{l}^{m}(\cos \theta) Ylm(θ,ϕ)=(−1)m4π(2l+1)(l+∣m∣)!(l−∣m∣)!eimϕPlm(cosθ)
其中 L e g e n d r e Legendre Legendre 函数 P l m ( cos θ ) P_{l}^{m}(\cos \theta) Plm(cosθ):
虽然，也很复杂。。但是耐不住MATLAB自带这个函数啊，甚至MATLAB还给出了如何通过该函数生成球谐函数的实例，简直太贴心了，爱了爱了：
dx=pi/60;
col=0:dx:pi;
az=0:dx:2*pi;
[phi,theta]=meshgrid(az,col);
% 计算 l=3 的网格上的 P_{l}^{m}(\cos \theta)
l=3;
Plm=legendre(l,cos(theta));
% 由于 legendre 为 m 的所有值计算答案，因此 Plm 会包含一些额外的函数值。
% 提取 m=2 的值并丢弃其余值。
% 使用 reshape 函数将结果定向为与 phi 和 theta 具有相同大小的矩阵。
m=2;
if l~=0Plm=reshape(Plm(m+1,:,:),size(phi));
end
% 计算Y_3^2的球谐函数值
a=(2*l+1)*factorial(l-m);
b=4*pi*factorial(l+m);
C=sqrt(a/b);
Ylm=C.*Plm.*exp(1i*m*phi);
% 球面坐标转换为笛卡尔坐标并绘图
[Xm,Ym,Zm]=sph2cart(phi,pi/2-theta,abs(real(Ylm)));
surf(Xm,Ym,Zm)
title('$Y_3^2$ spherical harmonic','interpreter','latex')
万事俱备开始编程：
工具函数
function psi=HydWave(n,l,m,x,y,z)
% @author:slandarer% 由坐标点计算向量模长及角度
r=vecnorm([x(:),y(:),z(:)]')';
theta=atan2(vecnorm([x(:),y(:)]')',z(:));
phi=atan2(y(:),x(:));% 恢复矩阵型状
r=reshape(r,size(x));
theta=reshape(theta,size(x));
phi=reshape(phi,size(x));% 利用MATLAB自带legendre函数计算球谐函数
Plm=legendre(l,cos(theta));
if l~= 0Plm=reshape(Plm(m+1,:,:),size(phi));
end
C=sqrt(((2*l+1)*factorial(l-m))/(4*pi*factorial(l+m)));
Ylm=C.*Plm.*exp(1i*m*phi);% laguerreL函数部分计算
Lag=laguerreL(n-l-1,2*l+1,2.*r./n);% 整合起来
psi=exp(-r./n).*(2.*r./n).^l.*Lag.*Ylm;
psi=real(conj(psi).*psi);
end
基本使用
[X,Z]=meshgrid(linspace(-30,30,120));
Y=zeros(size(X));psi=HydWave(4,2,0,X,Y,Z);
surf(X,Z,psi,'EdgeColor','none')
view(2);caxis([0,2])
封面图制作
function HydWaveDemo
% @author:slandarer% 构造一个好看的颜色映射实属不易
map=[0.1294 0.0549 0.1725;0.2196 0.1608 0.2902;0.3882 0.1804 0.4941;0.4392 0.1922 0.4706;0.5333 0.2235 0.4392;0.6471 0.2588 0.3686;0.7137 0.2745 0.3294;0.7725 0.3059 0.2902;0.8510 0.3725 0.2275;0.9137 0.4196 0.1804;0.9608 0.5020 0.2000;0.9765 0.5529 0.2078;0.9804 0.6431 0.2549;0.9843 0.6627 0.2706;0.9765 0.7176 0.3412;0.9765 0.7686 0.4000;0.9765 0.8118 0.4902;0.9725 0.8510 0.5961;0.9882 0.9020 0.6667;1.0000 0.9451 0.8431;1.0000 0.9961 0.9804;1.0000 1.0000 1.0000];
Xi=1:size(map,1);Xq=linspace(1,size(map,1),800);
map=[interp1(Xi,map(:,1),Xq,'linear')',...interp1(Xi,map(:,2),Xq,'linear')',...interp1(Xi,map(:,3),Xq,'linear')'];% 情况信息列表
condList=[1,2,0,0,-10,10,0.01;2,3,0,0,-20,20,0.01;6,2,1,0,-12,12,0.03;7,3,1,0,-20,20,0.08;8,3,1,1,-20,20,0.08;11,2,1,1,-12,12,0.03;12,3,2,0,-23,23,0.35;13,3,2,1,-23,23,0.35;14,3,2,2,-23,23,0.35;16,4,0,0,-36,36,0.008;17,4,1,0,-30,30,0.2;18,4,1,1,-30,30,0.1;19,4,2,0,-30,30,0.95;20,4,2,1,-30,30,0.95;21,4,2,2,-30,30,0.95;22,4,3,0,-35,35,6;23,4,3,1,-35,35,6;24,4,3,2,-35,35,6;25,4,3,3,-35,35,1.8];fig=gcf;
fig.Color=[0,0,0];
set(fig,'InvertHardCopy','off');% 文本信息
axh=subplot(5,5,[3,4,5]);
axh.XLim=[0,3];
axh.YLim=[0,1];
axh.XTick=[];
axh.YTick=[];
axh.XColor=[0,0,0];
axh.YColor=[0,0,0];
axh.Color=[0 0 0];
text(axh,.04,0.8,'Hydrogen Electron Orbita','Color',[1,1,1].*.98,...'FontName','cambria','FontWeight','bold','FontSize',15)
text(axh,.04,0.2,'$\psi_{n, l, m}=e^{-r / n}\left(\frac{2 r}{n}\right)^{l}\left[L_{n-l-1}^{2 l+1}\left(\frac{2 r}{n}\right)\right] Y_{l}^{m}(\theta, \phi)$',...'Color',[1,1,1].*.98,'FontSize',12,'Interpreter','latex')% 文本信息2
axh2=subplot(5,5,[9,10]);
axh2.XLim=[0,2];
axh2.YLim=[0,1];
axh2.XTick=[];
axh2.YTick=[];
axh2.XColor=[0,0,0];
axh2.YColor=[0,0,0];
axh2.Color=[0 0 0];
text(axh2,.04,0.8,'Not normalized by:','Color',[1,1,1].*.98,...'FontName','cambria','FontSize',13)
text(axh2,.04,0.2,'$\sqrt{\left(\frac{2}{n a}\right)^{3} \frac{(n-l-1) !}{2 n(n+l) !}}$','Color',[1,1,1].*.98,...'FontName','cambria','FontSize',12,'Interpreter','latex')% 循环绘图
for i=1:size(condList,1)ax=subplot(5,5,condList(i,1));% 绘制密度分布[X,Z]=meshgrid(linspace(condList(i,5),condList(i,6),120));Y=zeros(size(X));psi=HydWave(condList(i,2),condList(i,3),condList(i,4),X,Y,Z);surf(X,Z,psi,'EdgeColor','none')caxis(ax,[0,condList(i,7)]);colormap(map);view(2)axis tight% 坐标区域修饰ax.XLim=[condList(i,5),condList(i,6)];ax.YLim=[condList(i,5),condList(i,6)];ax.DataAspectRatio=[1,1,1];ax.XTick=[];ax.YTick=[];ax.XColor=[1,1,1].*.4;ax.YColor=[1,1,1].*.4;ax.LineWidth=1.5;ax.Box='on';ax.FontName='cambria';ax.XLabel.String=['(',num2str(condList(i,2)),',',...num2str(condList(i,3)),',',...num2str(condList(i,4)),')'];drawnow
end
saveas(fig,'result.png')
end
这篇关于MATLAB定态氢原子波函数可视化的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
