数学物理方程三维可视化仿真作者：江萍杨华军何文森罗志华来源：《教育教学论坛》2013年第03期摘要：数学物理方程三维可视化仿真及创新实践训练是《数学物理方法》教学模式改革中的重要内容。

本文通过MATLAB程序求解二维菱形晶格光子晶体的电磁场本征值方程，绘制出二维能带曲线，并将结果三维可视化，体现出复杂数学物理问题的物理图像，解决大学生在课程学习过程中理解困难的教学问题，加强大学生编程实践能力和创新能力的培养。

关键词：本征值问题；三维可视化仿真；光子晶体；平面波展开法；能带结构中图分类号：G642.0 文献标志码：A 文章编号：1674-9324（2013）03-0247-03一、课程背景《数学物理方法》是理工科学生的基础课程之一，也是科研中常用的基本方法。

数学物理方法课程的内容繁多，公式推导繁杂，尽管教材中的例题通常具有明确的物理意义，但是从眼花缭乱的数学表达式中看出其中所表达的物理图像，不仅学生会觉得困惑、枯燥，教师也难免觉得棘手。

探索数学物理方法数值化教学的新方法，是数学物理方法课程教学中的一项重要工作，也是数学物理方法教学改革中的重要内容。

利用MATLAB数值求解数学物理方程，将传统教学手段与计算机仿真教学相结合，改变只用公式符号教学的模式[1]，令学生对复杂、抽象、烦琐的数学物理问题具有更深刻的理解。

本论文旨在进行数学物理方程仿真求解实践训练，着力培养大学生应用数学物理思想解决实际问题的能力。

本着“重理论、强实践、突创新”的教育理念，结合科技前沿，以光子晶体的电磁场理论作为实践内容，利用MATLAB对复杂的电磁场本征值问题进行计算机仿真求解，将结果三维可视化，以此来展现复杂电磁场问题的物理图像，对培养大学生创新能力具有重要意义。

二、光子晶体电磁理论基础在利用分离变量法求解数学物理方程时，最后都归结到求解本征值问题。

在用本征函数系展开法解数学物理方程时，也要对所用的本征函数系有较好的理解[2]。

所以，各种本征函数系在数学物理方程课程的学习中有非常重要的地位。

周期结构对电磁波的调控是物理学领域的基础问题。

光子晶体是由介电常数周期排列形成的一种合成材料，是非均匀介质中少数可以严格遵循电磁理论的新型人工材料。

在一定的晶格常数和介电常数条件下，布拉格散射使在光子晶体中传播的电磁波受到调制形成类似于电子的能带结构[3]。

利用计算机仿真求解光子晶体中的复杂本征值问题，可以帮助学生熟悉并更好地掌握本征函数系的性质和求解方法。

1.理想二维光子晶体的结构。

假设介电常数为εa，半径为r的介质柱平行于z轴，背景介质的介电常数为εb，在（x，y）平面内的晶格常数为a，θ为相邻基矢a1和a2之间的锐角，当θ=90°和60°时，分别为正方形晶格和正三角形晶格。

（x，y）平面的傅里叶变换空间为倒易空间（如图1所示），对应于由波矢k定义的频谱。

根据几何关系，（x，y）平面内的基矢（a1，a2）和倒易空间的基矢（b1，b2）分别为a1=aex，a2=a（cosθex+sinθey）；b1=■ex-■ctgθey，；b2=■ey。

在倒易空间中，Γ、T、N、X和M等点为布里渊区的高对称点，它们所构成的多边形区域（深灰色部分）称为不可约布里渊区。

不可约布里渊区是倒易空间中最小的、可重复的区域，可以映射出电磁波在整个光子晶体中的传输特性[4]。

对称点坐标分别为：Γ=（0，0），T=■（1，■），N=■（1，■ctgθ，■）；X=■（0，■），M=■（■ctgθ-1，■）。

2.理想二维光子晶体的本征值问题。

平面波的指数形式表示为H（r，t）=H（r）e-iωt （1）E（r，t）=H（r）e-iωt （2）联立无源Maxwell方程组，分别得到电场和磁场的传播方程？塄×■？塄×H（r）=■■H（r）（3）■？塄×？塄×E（r）=■■E（r）（4）ε（r）是光子晶体介质分布的周期函数，本征值方程（3）和（4）式与电子材料中的周期性势场问题的Schr？觟dinger方程类似，称为光子晶体的支配方程[5]。

本征场H（r）和E （r）分别对应于理想二维光子晶体中横磁（TM）模式和横电（TE）模式的空间形态，通过求解本征值（ω/c）2，可以得到频率ω与波矢k之间的色散关系，即光子晶体的能带结构。

3.光子晶体中的平面波展开。

根据Bloch理论，将光子晶体本征场用平面波展开为HK（r）=■H■（G）e■ （5）EK（r）=■E■（G）e■ （6）G为倒格矢，将（5）和（6）式分别代入本征值方程（3）和（4）式，利用平面波基{G，exp[i（k+G）gr]，…}的正交性[6]，得到如下关于电磁场展开系数的本征值方程■k%（G-G'）（k+G）g（k+G'）Hk（G'）=■Hk（G）（7）■k%（G-G'）（k+G'）g（k+G'）Ek（G'）=■Ek（G）（8）矩阵k%（G）是k（r）=1/ε（r）=k%（G）e■的傅立叶展开系数k%（G）=■■k%（r）e■dr■=■■■e■dr2+■（■-■）■e■dr2 （9）u表示一个周期单元，Au为周期单元横截面的面积。

c表示一个散射单元横截面上的积分边界。

（9）式右边包含了G≠0和G=0的项k%（G）=■+（■-■）fr，G=02fr（■-■）■，G≠0 （10）其中J1（GR）为第一类贝塞尔函数，fr为填充比。

三、仿真求解电磁场本征值问题我们通过计算机仿真求解TM模式电磁场本征值方程（7）式，获得二维菱形晶格光子晶体的本征频率ωk与波矢k之间的色散关系，绘制出能带曲线。

1.光子晶体的数学建模。

对于θ=70°的二维菱形晶格光子晶体，背景介质的介电常数为εb=12，空气柱的半径r=0.4a。

仿真步骤和MATLAB程序如下：①定义光子晶体的结构参数。

ea=1；eb=12；R=0.4；sita=70；a=1；a1=a*[1，0]；a2=a*[cos（sita*pi/180），sin（sita*pi/180）]；b1=2*pi/a*[1，-1/tan（sita*pi/180）]；b2=2*pi/a*[0，1/sin（sita*pi/180）]；fr=pi*R*R/abs（a1（1）*a2（2）-a1（2）*a2（1））；②定义倒易空间中对称点的坐标。

Point（1，：）=[0，0]；Point（2，：）=pi/a*[1，（1-cos（sita））/sin（sita）]；Point（3，：）=pi/a*[1-（1-cos （sita））/sin（sita）*1/tan（sita），1/sin（sita）]；Point（4，：）=pi/a*[0，1/sin（sita）]；Point（5，：）=pi/a*[（1-cos（sita））/sin（sita）*1/tan（sita）-1，1/sin（sita）]；Point（6，：）=[0，0]；③产生一个20×20的矩阵，确定平面波的波数NPW，定义倒格矢G。

DimForG=20+1；NPW=DimForG*DimForG；gtemp=-10：10；gtemp1=repmat（gtemp，DimForG，1）；Gx=b1（1）*gtemp1+b2（1）*gtemp1'；Gy=b1（2）*gtemp1+b2（2）*gtemp1'；Gx=Gx（：）'；Gy=Gy（：）'；Gx_m=repmat（Gx，NPW，1）；Gx_n=Gx_m'；Gy_m=repmat（Gy，NPW，1）；Gy_n=Gy_m'；G=sqrt（（Gx_m-Gx_n）.*（Gx_m-Gx_n）+（Gy_m-Gy_n）.*（Gy_m-Gy_n））；④确定κ（r）=1/ε（r）的傅里叶展开系数。

ek0=fr/ea+（1-fr）/eb；ekc=（1/ea-1/eb）*fr*2；GR=G*R；na=find（GR==0）；GR（na）=1；ek=ekc*besselj（1，GR）./GR；ek（na）=ek0；2.仿真计算光子晶体TM模式能带曲线。

①定义倒易空间波矢路径。

用Keach代表波矢路径上的取值密度，Ktype为对称点的数目，第一布里渊区内沿波矢路径Γ→T→N→M→Γ的仿真程序为：Ktype=5；Keach=6；x=[]；y=[]；for n=1：Ktype x1=linspace（Point（n，1），Point（n+1，1），Keach+1）；y1=linspace （Point（n，2），Point（n+1，2），Keach+1）；x=[x，x1（1：Keach）]；y=[y，y1（1：Keach）]；end②求解本征值方程。

eigval=[]；for m=1：Ktype*Keachkx=x（m）；ky=y（m）；KGmn=sqrt（（kx-Gx_m）.^2+（ky-Gy_m）.^2）.*sqrt（（kx-Gx_n）.^2+（ky-Gy_n）.^2）；H=KGmn.*ek；eigvalue=sort（eig（H））；eigval=[eigval，eigvalue（1：20）]；endeigval=[eigval，eigval（：，1）]；eigval=real（sqrt（eigval）*a*0.5/pi）；③绘制二维能带曲线。

for m=1：KtypeD（m）=sqrt（（Point（m+1，1）-Point（m，1））^2+（Point（m+1，2）-Point（m，2））^2）；xtemp（m，：）=linspace（0，D（m），Keach+1）；endx=xtemp（1，1：Keach）；Dtotal=0；for m=2：KtypeDtotal=Dtotal+D（m-1）；x=[x，xtemp（m，1：Keach）+Dtotal]；endx=[x，xtemp（Ktype，Keach+1）+Dtotal]；x=x/max（x）；plot（x，eigval）；修饰过后的二维菱形晶格光子晶体TM偏振模式能带曲线如图2所示。

四、本征值函数的三维可视化仿真绘制三维等频面，关键是建立波矢平面（kx，ky）内二维点阵的坐标，再求解出每个点对应特征值，仿真步骤和MATLAB程序为：1.定义波矢（kx，ky）平面内点阵的坐标Keach=36；x=linspace（Point（1，1），Point（2，1），Keach+1）；y=linspace（Point（1，2），Point（4，2），Keach+1）；x1=[-x（Keach+1：-1：2），x（1：Keach+1）]；y1=[-y（Keach+1：-1：2），y（1：Keach+1）]；2.求解本征值方程。