［收稿时间］2019-12-18
［基金项目］哈工大（威海）研究生教育教学改革研究项目（WH2019014）；哈工大研究生教改研究项目（JGYJ-2019036）。［作者简介］周洪娟（1980-），女，山东烟台人，博士，副教授，主要从事电磁理论方面的教学和研究工作。
［摘
要］电磁场边值问题的求解是电磁理论教学中的难点和重点。课题组以简单的静态二维电场边值问题为例，同时采用 解析法和数值法求解，基于Matlab 仿真平台编程实现，从解析法和数值法的结论互相呼应的角度来逐层次地设计实验，使学生对电磁场边值问题求解方法、抽象复杂的数学结论以及唯一性定理产生感性认识。 ［关键词］电磁场边值问题；唯一性定理；解析法；数值法
［中图分类号］O411.1［文献标识码］A ［文章编号］2095-3437（2021）02-0004-04
2021年2University Education
“电磁场理论”或“电磁场与电磁波”是工科院校电子信息、无线电技术类专业的一门重要基础课，其涉及的矢量微积分公式繁多、概念抽象，需要学生具备较为扎实的数学和物理基础，学生普遍反映难度大。“电磁场理论”这门课的教学虽然要侧重电磁场、电磁波的基础理论，但也要注重与工程实践的结合，为工科院校的学生在相关课程以及方向的学习研究提供较为直接的理论指导。这其中，电磁场边值问题的求解就是联系电磁场麦克斯韦方程等基础理论与各种复杂工程实践，如天线设计、电磁干扰与电磁兼容以及雷达散射截面积等相关应用的桥梁[1-7]，但由于其涉及数理方程等复杂数学理论，使之成为本科教学中的难点。
整形归来2
电磁场边值问题指的是满足特定偏微分方程和边值条件的数理方程，静态电、磁场的边值问题的求解是指满足泊松方程或拉普拉斯方程和指定边值条件的数理方程的求解。电磁场边值问题的求解方法主要分为解析法和数值法两大类。解析法是指能从电磁理论出发通过公式推导可直接得到所求解问题的精确表达式的方法，该类方法通常只适合一些边界形状简单的特殊边值问题的求解，如边界形状为规则的平面状、球状或圆柱状。数值法是将求解区域划分成离散的网格、待求解的连续变量离散成求解区域中一系列离散点处的求解的一系列方法，常见的有限差分法[8-10]、有限元法[11]和矩量法[11]。由于该类方法可灵活处理各种复杂的边界问题，在实践工程计算中被广泛采用。现阶段本科教学中，一般只强调一些特殊边值问题的解析求法，如分离变量法和镜像法，而由于授课学时等的限制，对更通用的数值法讲解不透，甚至一带而过，这样就难以达到服务于工程实践的真实目的。
本文以二维静电场边值问题的求解为例，设计结合
解析法和有限差分数值法的自我验证式实验，提高学生对边值问题求解方法和对边值问题求解的唯一性定理的理解。
一、待求解的边值问题及其解析解
待求解的二维静电场边值问题如式（1）所示，满足第一类边值条件，求解区域内面电荷密度为ρ，x 轴和y 轴的长度分别为a 和b ，U 0为一个常数，求解区域内的电位满足的偏微分方程和边值条件写成：
ìí
îïïïï
ttv
∂2φ∂x
2+
∂2φ
∂y 2=-ρ(x ,y )εφ(0,y )=0,φ(a ,y )=0
(0≤y ≤b )φ(x ,0)=0,φ(x ,b )=U 0(0≤x ≤a )
（1）
根据静态电磁场问题的唯一性定理，只要求出的位函数满足相应的泊松方程或拉普拉斯方程，又满足给定的边界条件，则此位函数就是所求的唯一正确解。该边值问题可以由叠加法分解成两个边值问题，如式（2）和
（3）所示，则有φ=φ1+φ2。
ìí
îïïïï
∂2φ1∂x
2+
∂2φ1
∂y 2=0φ1(0,y )=0,φ1(a ,y )=0
(0≤y ≤b )φ1(x ,0)=0,φ1(x ,b )=U 0(0≤x ≤a )（2）
ìí
îïïïï
∂2φ2∂x
2+
∂2φ2∂y 2=-ρ(x ,y )
εφ2(0,y )=0,φ2(a ,y )=0
(0≤y ≤b )φ2(x ,0)=0,φ2(x ,b )=0(0≤x ≤a )
（3）
对边值问题（2），由于电位满足拉普拉斯方程，可以采用分离变量法直接求解[5]，即将电位分解为两个独立变量的函数相乘的形式，即φ1(x ,y )=f (x )g (y )，由x 向
我的校长生涯边界处的边值条件得f (0)=f (a )=0，因此有f (x )=
电磁场数值计算仿真实验设计
周洪娟
金
涛
陈
铖
（哈尔滨工业大学（威海）信息科学与工程学院，山东
威海
264209）
4
University Education
∑n
C n sin(
n πx
a )，由y 向边界处的边值条件得g (0)=0，g (
b )≠0，
因此分离函数g (y )具有sin h (n πy
a
)的形式，进一步由边界条件和三角函数的正交性得到最终解析结
论如式（4）所示。
φ1(x ,y )=
∑n =1,3,5,⋯
∞4U 0n πsinh(n πb /a )sin(n πx a )sinh(n πy
a )（4）
对边值问题（3），电位满足泊松方程，可以采用级数
展开法求解[12]。考虑四个齐次边界条件，令解为二重正弦级数，即
φ2(x ,y )=∑m =1∞
∑n =1∞
A mn sin ()mπx a sin
()
nπy b （5）
将式（5）带入泊松方程，并由三角函数的正交性，可得系数表达式为：A mn =
4
ab
éëêêùûúú()mπa 2+()
nπb 2
-1
∫0a ∫
0b ρ(x ,y )εsin ()mπx a sin ()
nπy b dxdy （6）至此，边值问题（1）的解析表达式推导完毕。以上
求解在唯一性定理的统一框架下，采用分离变量法、叠加法和级数展开法联合求解，系数求解中需要用到三角函数的正交性原理，推导过程比较复杂，公式繁多，导致学生的接受难度大，因此电磁场边值问题的求解历来是本科教学中的难点。
二、数值法原理
本次实验采用有限差分法求解，该方法代码实现较为简单，比较适合本科阶段的学生自行编写代码实现，既可以体会数值计算方法的精要，在实现上又不占用太多的精力。另外，其求解联立方程的不同解法也可以体现数值计算中保证计算精度的同时加快收敛速度的需求，让学生对数值求解方法有一个比较全面的认识。
首先将求解区域划分成正方形网格，下标i 和j 分别代表节点的x 向和y 向的序号，网格长度为h ，则泊松方程的差分方程如式（7）所示[9]
。
φi ,j =
14
()
φi -1,j +φi ,j -1+φi +1,j +φi ,j +1+h 2
ρ
ε（7）
若求解区域划分的节点数为N ，则得到N 个联立方程组，进一步考虑已知边界条件，主要的求解方法包括矩阵求逆、高斯消元和迭代法。迭代法是先对每个待求节点赋初值，按照一定顺序和公式（7）逐步更新每个节点的电位，直至每个节点前后更新值误差在给定的精度范围内。由于迭代法更适合求解大规模边值问题，更具有通用性，同时实现简单，因此本次实验采用迭代法求解。迭代法分简单迭代法和超松弛迭代法两种，简单迭代法的更新公式如式（8）所示，即直接根据泊松方程的差分方程（7）用上一次迭代的值来更新本次迭代的值：
φi ,j n +1=
1
4
(
)
φi -1,j n +φi ,j -1n +φi +1,j n +φi ,j +1n +h 2
ρi ,j ε（8）
其中，上标n 和n+1代表迭代次数的序号。若更新
顺序为从左到右、从下到上，由于简单迭代法中每次更新时，左边和下边的节点已经更新过，为提高收敛速度引入了超松弛迭代法。超松弛迭代法通过引入超松弛因子α（1<α<2），充分利用已经得到的新值，在得到每点新值时立即把旧值冲掉，以提高迭代效率。设每次迭代时，节点更新顺序为从左到右、从下而上，则更新公式如式（9）和（10）所示。
φ(n +1)i ,j =φ(n )i ,j +a [φ
(n +1)i ,j -φ(n )i ,j ]（9）
无差异曲线
φ
i,j n +1=1
4
(
)
φi -1,j n +1+φi ,j -1n +1+φi +1,j n +φi ,j +1n +h 2
ρi ,j ε（10）
当α=1时φ(n +1)i ,j =φ(n +1)i ,j ，此时称为松弛迭代。超松
弛因子α的大小决定着收敛速度，迭代次数最小的超松
弛因子称为最优超松弛因子。若计算区域为长方形，设划分成正方形网格时两边的节点数分别为Nx 和Ny ，则最优松弛因子αopt 的经验值如式（11）和（12）所示[12]。
αopt =
8-4
4-t 2
t 2
（11）
t =cos
(
)πNx -1+cos
()
π
Ny -1
（12）
三、实验设计
本次实验要求学生自行编写代码来实现解析结论分析、数值法求解等。MATALB 提供了一种高效的基于矩阵运算的仿真环境，具有丰富的函数库和画图功能，代码编写简单，在理论研究中被广泛采用。本次实验基于MATLAB 平台进行有限差分代码编写及后期数据的画图分析，实验设计从学生可以自我验证计算结论的角度出发，实验过程及要求设计如下所示：1.基于公式（4）（5）（6）分析基于叠加法、分离变量法和级数展开法的解析解；
2.利用简单迭代有限差分法求解，采用不同的网格尺寸，数值解与解析法结论对比，分析求解结果的精确度；
3.利用超松弛迭代法求解，确定不同网格尺寸下收
敛速度最快的最优值松弛因子，并与公式（11）（12）经验值进行对比。
本次实验是以边值问题（1）的求解为例展开的，可以根据实际问题需要调整待求解的边值问题。通过调整网格尺寸，让学生体会数值计算中网格离散度对计算结论的精确度的影响，以及对计算资源方面的要求；超松弛因子的求解又让学生体会到数值计算中保证求解精度提高收敛速度方面的切实要求。同时，该实验无论从采用叠加法的解析求解过程，还是解析解与数值解呼应的验证过程，都遵循着电磁场边值问题求解的唯一性定理，即同一边值问题可以采用不同的方法求解，但所
5
仪征市第二中学得的解必须满足指定的偏微分方程和边值条件。唯一性定理在本科课堂教学中通常只是一掠而过，学生难以体会其在电磁场边值问题求解中的重要地位，通过该实验可以让学生深刻理解唯一性定理的意义和重要性，真正起到了实验辅助教学的目的。
四、实验结果分析
设求解区域尺寸a=10m，b=8m，U0=100V，电荷面密度ρ=10x（y-1）/ε0pC/m2；数值计算中采用三种网格尺寸h=0.1m、h=0.5m和h=1m，要求两次迭代差值最大为10-10V，以下为对应的分析结论。图1中的（a）（b）（c）（d）分别为解析法，网格尺寸h=1m、h=0.5m和h=0.1m 的有限差分法得到的整个计算区域内的电位分布彩图，其中有限差分法采用简单迭代法。图2中的（a）和（b）分别为沿着y=1m和y=5m两条线上三种网格尺寸下得到的数值解相对于解析解的相对误差。可见随着网格尺寸降低，数值解逐渐趋近于解析解，表明网格分得越细计算精确度越高，同时也证明了在严格一致的边界条件和偏微分方程下，同一边值问题多种解法的结论是趋于一致的，
即解的唯一正确性。
图1带求解边值问题的电位分布彩图：（a）解析解；（b）h=1m
数值解；（c）h=0.5m数值解；（d）h=0.1m
数值解。
图2不同网格尺寸下数值解相对误差分布：（a）y=1m沿线；
（b）y=5m沿线
图1中的（b）（c）（d）的迭代次数分别为393、1508和
32771，可见求解精度的提高是以计算时长为代价的，这
对所有的数值计算方法都是不可回避的事实，因此需要
在保证求解精度的同时尽量提高收敛速度，降级计算时
间，提高计算效率。本实验中进一步采用超松弛迭代法
加速收敛速度，设要求的求解精度与简单迭代法相同。
实验结果表明，超松弛迭代法在保证计算结论精确度的
前提下，可以得到更快的收敛速度，且收敛速度与松弛
因子有关。图3中的（a）（b）（c）为三种网格尺寸h=0.1m、
h=0.5m和h=1m下不同松弛因子对应的迭代次数，在具
体实验中，首先采用0.1间隔来遍历搜索，得到最小迭代
次数分别为45、90和913，对应的松弛因子分别为1.5、
1.7和1.9。由式（11）（12）可得h=0.1m、h=0.5m和h=1m
三种网格尺寸下的最优松弛因子经验值分别为1.484、
掘金黑客
1.699和1.931，仿真结论与经验值基本吻合。若要得到
分辨率更高的αopt数值，可以进一步在最小值附近进行
细搜索，如图3中的（c）中所示，得到最优松弛因子αopt=
1.94，迭代次数为478。还有需要说明的是α=1.93时，迭
代次数为482，两者相差很小，说明经验公式是非常可信
的。还有需要说明的是，在松弛因子αopt=1时，h=0.1m、
h=0.5m和h=1m三种网格尺寸下的迭代次数分别为
16947、779和204，相对于简单迭代法的收敛速度已经有
很大改进。
最终简单迭代法和最优超松弛迭代的迭代次数对
比如表1所示，其中的比值代表采用超松弛迭代法的最
小迭代次数与简单迭代法的迭代次数的比值。可见，通
过调整超松弛因子，在保证计算精度的前提下，可以极
大地提高计算效率，并且节点数越多改善效果越明显，
这对大型的边值问题的数值求解具有非常重要的意义。6
University
Education
图3三种网格尺寸下迭代次数随松弛因子的变化：（a）h=1m；（b）h=0.5m；（c）h=0.1m。
表1超松弛迭代法和简单迭代法的迭代次数对比
五、结束语
电磁场边值问题的求解由于涉及的数学公式繁杂，方法多样化，历来是电磁场与电磁波本科教学中的难点。本文以一个二维静电场边值问题的求解为例，首先基于叠加法、分离变量法和级数展开法推导其解析解；然后进一步采用有限差分法进行数值求解，讨论了离散网格大小对数值计算精度的影响，并采用超松弛迭代法，让学生认识到数值算法中改善计算效率的问题。本实验内容的每一步结论都有参考结论，学生可以自我验证；基于MATLAB的算法代码比较简单，学生通过自己动手编写代码，加深了对算法的理解。该实验设计同时从解析法和数值法多个角度加深了学生对抽象的唯一性定理的理解及其在电磁场边值问题求解中的重要地位，帮助学生理解电磁场的数值计算方法，进一步促进学生
对电磁场理论中的基本方程、边界条件及各类二阶偏微分方程的活学活用。
［参考文献］
［1］周禹龙，曹彰玉，高军，等.双频频率选择表面及其在微带天线宽带RCS减缩中的应用［J］.电子与信息学报，
2017，39（6）：1446-1451.
［2］严雪飞，朱长青.线面结构天线的电磁建模与计算［J］.微波学报，2018，34（1）：6-10+20.
［3］刘战合，王菁，姬金祖，等.典型布局飞机电磁散射特性数值计算研究［J］.航空工程进展，2018，9（3）：341-347.［4］刘涛，曹祥玉，高军，等.宽带低RCS超表面天线阵设计［J］.电子与信息学报，2019，41（9）：2095-2102.
［5］韩书键，白锐锋，于涛.目标宽频带RCS的快速算法研究［J］.国防科技，2019，40（3）：7-12.
［6］罗晨，葛勇，王舒，等.基于矩量法的某轰炸机电磁散射特性研究［J］.空天防御，2019，2（3）：66-72.
［7］刘阳，周海京，郑宇腾，等.高性能计算在目标电磁散射特性分析中的应用［J］.电波科学学报，2019，34（1）：3-
11.
［8］葛德彪，闫玉波.电磁波时域有限差分方法（第三版）［M］.西安：西安电子科技大学出版社，2011.
［9］谢处方，饶克谨.电磁场与电磁波（第4版）［M］.北京：高等教育出版社，2006．
［10］陈伟军，龙世瑜，刘如军.电磁场课程设计中新的数值计算方法探索［J］.实验科学与技，2017，15（6）：76-79.［11］倪光正，杨仕友，邱捷，等.工程电磁场数值计算（第2版）［M］.北京：机械工业出版社，2010.
［12］Matthew N.O.Sadiku.MATLAB模拟的电磁学数值技术（第3版）［M］.喻志远，译.北京：国防工业出版社，2016.
［责任编辑
责任编辑：：钟岚］
7
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议