1. 基于matlab－PDE Toolbox的泊松（拉普拉斯）方程求解
在二维电磁场的有限元法计算中，用矩阵方程编制的计算程序长、大，而又复杂，且输入数据要化费很大的劳动。而MATLAB是一种以矩阵运算为基础的交互式语言，它是采用有限元法来求解偏微分方程的。因此在计算中，我们选用MATLAB提供的一个用户图形界面（GUI）的偏微分方程工具箱（PDE Toolbox）进行数值求解，采用的是三角形网格自动剖分。下面举例说明。
【例1－1】 横截面为矩形的无限长槽由3块接地导体板构成，如图3－3所示，槽的盖板接直流电压100V，求矩形槽的电位分布。
解：这是二维平面场问题。由于电位函数
利用麦克斯韦方程
式中，为介电常数，
和关系式
，得到泊松方程
和电场强度
之间关系为
为体电荷密度。由于所求区域内体电荷密度
，得到拉普拉斯方程：
其边界满足狄里赫利（Didchlet）条件：
， ，
本题运用MATLAB的偏微分方程工具箱（PDE Toolbox）进行数值求解。在命令窗口中输入命令pdetool，打开PDE图形用户界面，计算步骤为：
（1）网格设置：选择菜单Options下的Grid和Grid Spacing…，将X-axis linear Spacings设置为［-1.5：0.2：1.5］，Y-axis linear Spacings取Auto。 （2）区域设置：选择菜单Draw下的Rectangle／Square或按应用模式。
（4）输入边界条件：进入Boundary Mode或按
，输入：
，画矩形。
（3）应用模式设置：在工具条中单击Generic Scalar下拉列表框，选Electrostatics（静电学）
1、左边界：狄里赫利（Diriehlet）条件：h＝1，r＝0。 2、右边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 3、上边界：狄里赫利（Dirichlet）条件：h＝1，r＝100。 4、下边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 （5）方程参数设定：点击密度rho为0。
打开PDE Spacification对话框，设介电常数epsilon为1，体电荷
（6）网格剖分：点击按钮或选择菜单Initializc Mesh，加密网格再单击。
（7）图形解显示参数设置：点击菜单Plot下的Parameter或按设置等位线条数，在Colormap 选择不同的色图后，点击出电位的等位线和电场线分布图。 （8）解方程：点击按钮
，在Plot Selection对话框
中选择Color、Height（3-D plot）、Plot in x-y grid和 Show mesh四项，并在Contour plot levels
按钮画出电位的三维曲面图；选择
按钮后画
Color、Contour和Arrows三项，设置等位线条数和选择不同的色图参数后，点击
或选择Solve PDE菜单，结果如图3－4所示。
本题可由分离变量法直接求解，但其理论解是无穷级数的和，不易直观理解。当然也可以与上一章一样用有限差分法求解。为了分析场域内的电位分布，最简单的方法就是用MATLAB的偏微分方程工具箱进行有限元计算，只要在PDE图形用户界面按计算步骤一步步操作即可。读者可与上一章有限差分法的结果比较，它们的图形结果是一样的。
【例1－2】 设一个长直接接地金属槽，如图3－7所示，其侧壁和低面电位为零，顶盖电位为
，求槽内的电位分布。
解：金属槽中的电位函数
满足拉普拉斯方程
其边界条件满足混合型边值问题
，，
本题运用MATLAB的偏微分方程工具箱（PDE Toolbox）进行数值求解。在命令窗口中输入命令pdetool，打开PDE图形用户界面，计算步骤为：
（1）网格设置：选择菜单Options下的Grid和Grid Spacing…，将X-axis linear Spacings设置为［-1.5：0.2：1.5］，Y-axis linear Spacings取Auto。 （2）区域设置：选择菜单Draw下的Rectangle／Square或按学）应用模式。
（4）输入边界条件：进入Boundary Mode或按
，输入：
，画矩形。
（3）应用模式设置：在工具条中单击Generic Scalar下拉列表框，选择Electrostatics（静电
1、左边界：狄里赫利（Diriehlet）条件：h＝1，r＝0。 2、右边界：纽曼（Neumann）条件：g＝0，q＝0。
3、上边界：狄里赫利（Dirichlet）条件：h＝1，r＝ 4、下边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 （5）方程参数设定：点击密度rho为0。
（6）网格剖分：点击
。
打开PDE Spacification对话框，设介电常数epsilon为1，体电荷按钮或选择菜单Initializc Mesh，加密网格再单击
。
（7）图形解显示参数设置：点击菜单Plot下的Parameter或按设置等位线条数，在Colormap 选择不同的色图后，点击线和电场线分布图。 （8）解方程：点击按钮
或选择Solve PDE菜单，结果如图3－8所示
，在Plot Selection对话框中
选择Color、Height（3-D plot）、Plot in x-y grid和 Show mesh四项，并在Contour plot levels
按钮画出电位的三维曲面图；选择Color、
按钮后画出电位的等位
Contour和Arrows三项，设置等位线条数和选择不同的色图参数后，点击
本题是混合边值问题，只要用MATLAB的偏微分方程工具箱按相应的边值条件输入就能得到结果图。而对有限差分法来说，则需要根据不同的边界条件，对程序作相应的修改。而对理论解而言，混合边值问题可能使求解过程变得很繁。显然，边界条件的变换，并不会使MATLAB的偏微分方程工具进行有限元计算变难，这是它的优点所在。
【例1－3】 设两个同轴矩形金属槽如图3－5所示，外金属槽电位为零，内金属槽电位为求槽内的电位分布。
解：同轴矩形金属槽中的电位函数
满足拉普拉斯方程
，
其边界条件满足第一类边界条件（Dirichlet边界条件）问题
本题运用MATLAB的偏微分方程工具箱（PDE Toolbox）进行数值求解。在命令窗口中输入命令pdetool，打开PDE图形用户界面，计算步骤为：
（1）网格设置：选择菜单Options下的Grid和Grid Spacing…，将X-axis linear Spacings设置为［-1.5：0.2：1.5］，Y-axis linear Spacings取Auto。 （2）区域设置：选择菜单Draw下的Rectangle／Square或按为SQ1，第二个正方形名为SQ2，则所求区域为SQ1－SQ2。
（3）应用模式设置：在工具条中单击Generic Scalar下拉列表框，选择Electrostatics（静电学）应用模式
（4）输入边界条件：进入Boundary Mode或按
，输入：
，画正方形。若第一个正方形名
1、内左边界：狄里赫利（Diriehlet）条件：h＝1，r＝100。 2、内右边界：狄里赫利（Dirichlet）条件：h＝1，r＝100。 3、内上边界：狄里赫利（Dirichlet）条件：h＝1，r＝100。 4、内下边界：狄里赫利（Dirichlet）条件：h＝1，r＝100。 5、外左边界：狄里赫利（Diriehlet）条件：h＝1，r＝0。 6、外右边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 7、外上边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 8、外下边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 （5）方程参数设定：点击密度rho为0。
（6）网格剖分：点击
按钮或选择菜单Initializc Mesh，加密网格再单击
。
（7）图形解显示参数设置：点击菜单Plot下的Parameter或按在Colormap 选择不同的色图后，点击 （8）解方程：点击按钮
，在Plot Selection对话框
打开PDE Spacification对话框，设介电常数epsilon为1，体电荷
中选择Color、Height（3-D plot）和 Show mesh三项，并在Contour plot levels设置等位线条数，
按钮画出电位的三维曲面图；选择Color、Contour和Arrows
按钮后画出电位的等位线和电场线分布图。
三项，设置等位线条数和选择不同的色图参数后，点击
或选择Solve PDE菜单，结果如图3－6所示。
【例1－4】 在导电模拟的实验研究中，制备了如图3－15所示的二维电场模型，其内两种导电媒质的电导率分别为
和
，它们在场域的对角线
上接合。电极间外施电10V。求由此两种媒质构成的二维
电流场内的电位分布。
解：由于存在两种不同的导电媒质，所以应分别定义相应的电位函数它们在直角坐标系中，二维电流场内的电位函数满足拉普拉斯方程：
和
，
其边界条件满足混合型边值问题：
，
，
，
（此条件由MATLAB自动处理）
本题运用MATLAB的偏微分方程工具箱（PDE Toolbox）进行数值求解。在命令窗口中输入命令pdetool，打开PDE图形用户界面，计算步骤为：
（1）网格设置：选择菜单Options下的Grid和Grid Spacing…，将X-axis linear Spacings设置为［-1.5：0.2：1.5］，Y-axis linear Spacings取Auto。 （2）区域设置：选择菜单Draw下的Polygon或按学）应用模式。
，画上三角形、下三角形。
（3）应用模式设置：在工具条中单击Generic Scalar下拉列表框，选择Electrostatics（静电
（4）输入边界条件：进入Boundary Mode或按，输入：
1、左边界：纽曼（Neumann）条件：g＝0，q＝0。 2、右边界：纽曼（Neumann）条件：g＝0，q＝0。 3、上边界：狄里赫利（Dirichlet）条件：h＝1，r＝10。 4、下边界：狄里赫利（Dirichlet）条件：h＝1，r＝0。 （5）方程参数设定：点击
在PDE Spacification对话框内设定拉普拉斯方程
在上三角形选择Elliptic：c＝1；a=0；f＝0。 在下三角形选择Elliptic：c＝1；a=0；f＝0。 （6）网格剖分：点击
按钮或选择菜单Initializc Mesh，加密网格再单击
。
（7）图形解显示参数设置：点击菜单Plot下的Parameter或按择不同的色图后，点击 （8）解方程：点击按钮
按钮画出电位的等位线和电场线分布图。
或选择Solve PDE菜单，结果如图3－9所示。
，在Plot Selection对话框
中选择Color、Contour和Arrows三项，并在Contour plot levels设置等位线条数，在Colormap 选
由本题看到，根据不同的情况，可以把性质不同的区域分别用MATLAB的偏微分方程工具箱进行有限单元法的计算，而每个区域的计算与前面的类似，但要注意不同区域相邻边界的值的连续性或一致性，如
2.1 拉普拉斯方程的有限差分形式
二维拉普拉斯方程可以用有限差分法进行近似计算。首先把求解的区域划分成网格，把求解区域内连续的场分布用求网格节点上的离散的数值解代替。当然，把网格分得充分细，才能达到足够的精度。网格的划分有不同的方法，本书只讨论正方形和三角形网格划分。
如图2－1所示，在二维平面场中，我们将区域划分为许多正方形网格，每个网格的边长为（称为步长），两组平行线的交点称为网格节点，设网格节点位分别为
，
，
，
的电位为
，其上下左右4个节点的电为基点进行泰勒级数展开如下：
。
。在充分小的情况下，可以
把以上4式相加，得
（2.18）
及其以上
显然，在4式相加的过程中，的所有奇次方项都抵消了。所以，（2.18）式是在略去的项所得，其精度为的二次项。 由于场中任何点
都满足泊松方程
其中
为场源，则（2.18）式变为
，则二维拉普拉斯方程的有限差分形式：
（2.19）
对于无源场，
上式表示任意点的电位等于围绕它的四个等距离点的电位的平均值，距离越小则结果越精确。用方程（2.19）可以近似求解二维拉普拉斯方程，这样二阶偏微分方程就可以用差分代数方程来近似替代。
2.2 二维场域的边界条件
由于描述电磁场的偏微分方程是空间坐标的函数，只有在一组特定的边界条件下才能获得唯一解。大部分电磁场问题涉及到三种类型的边界条件：狄里赫利型、纽曼型和混合型边界条件。狄里赫利和纽曼两类型的边界条件通常分别称为第一类和第二类边界条件。 如图2－2所示，区域
的边界为曲线所包围，若二维场域边界上的电位
，
是一个事先已知的
电位函数，这种边界称为第一类边界条件（Dirichlet边界条件）。
如图2－3所示，若二维场域边界上电位法向的导数值满足这种边界称为第二类边界条件（Neumann边界条件）。
，
为已知数或一种连续函数，
如图2－4所示，若二维场域的边界条件为第一类边界条件和第二类边界条件的线性组合，这种边界条件称为混合型边界条件。
2.3 简单迭代法
简单迭代法的步骤是如下：
（1）先对某一网格点设一初值，这个初值完全可以任意给定，称为初值电位。虽然，问题的最终结果与初值无关，但初值选择估计得当，则计算步骤会得到简化。（当利用计算机来实现迭代计算时，为了简化程序初值电位一般可取为零值）。
（2）初值电位给定后，按一个固定顺序（点的顺序是从左到右，从下到上）依次计算每点的电位，即利用（2.19）式，用围绕它的四个点的电位的平均值作为它的新值，当所有的点计算完后，用它们的新值代替旧值，即完成了一次迭代。然后再进行下一次迭代，直到每一点计算的新值和旧值之差小于指定的范围为止。
简单迭代法的特点是用前一次迭代得到的网络点电位作为下一次迭代时的初值。例如，图2－5中的
点在
次迭代时计算公式为
（2.20）
2.4 超松弛法
简单迭代法在解决问题时收敛速度比较慢，一般来说，实用价值不大。实际中常采用超松驰法， 相比之下它有两点重大的改进。 第一是计算每一网格点时，把刚才计算得到的邻近点的电位新值代入，即在计算点即
（2.21）
上式称为松驰法或赛德尔法（relaxation method）。由于提前使用了新值，使得收敛速度加快。 第二，再把（2.21）式写成增量形式
这时每次的增量（即上式右边的第二项）就是要求方程局部达到平衡时应补充的量。为了加快收敛，我们引进一个松驰因子
式中的松驰因子
最佳值为
，将上式改写为
的电位时，把它左边的点
和下面的点
的电位用刚才算得的新值代入，
（2.22）
式中
、 为、
方向的网格数。不同的
值，可以有不同的收敛速度，其值范围一般为1
有一个
与2之间。即我们给予每点的增量超过使方程达到局部平衡时所需的值，这将加速解的收敛。最优值，如果选择比较恰当，收敛速度还将加快。
2.5 应用举例与计算程序
【例2－1】 横截面为矩形的无限长槽由3块接地导体板构成，如图2－6所示，槽的盖板接直流电压100V，求矩形槽的电位分布。
解：在直角坐标系中，矩形槽中的电位函数
满足拉普拉斯方程，
其边界条件满足第一类边界条件（Dirichlet边界条件）：
取步长
，，， 、
方向的网格数为
，
，共有16×10＝160个网孔，17×11＝187
个节点，其中槽内节点（电位待求点）15×9＝135个，边界节点（电位已知点）187－135＝52个。设迭代精度为10－6
，利用MATLAB编制的计算程序如下： % 例4-1：用简单迭代法求矩形槽内的电位分布 hx=17;hy=11; %设置网格节点数 v1=ones(hy,hx); %设置行列二维数组 %上下两行的Dirichlet边界条件值 v1(hy,:)=ones(1,hx)*100; v1(1,:)=zeros(1,hx);
%左右两列的Dirichlet边界条件值 for i=1:hy v1(i,1)=0; v1(i,hx)=0; end
v2=v1;maxt=1;t=0; %初始化 k=0
while(maxt>1e-6) %由v1迭代，算出v2，迭代精度为0.000001 k=k+1 %计算迭代次数 maxt=0;
for i=2:hy-1 %从2到hy-1行循环 for j=2:hx-1 %从2到hx-1列循环
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4; t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2 end
subplot(1,2,1),mesh(v2) %画三维曲面图 axis([0,17,0,11,0,100])
subplot(1,2,2),contour(v2,15) %画等电位线图 hold on
x=1:1:hx;y=1:1:hy
[xx,yy]=meshgrid(x,y); %形成栅格 [Gx,Gy]=gradient(v2,0.6,0.6); %计算梯度
quiver(xx,yy,Gx,Gy,-0.8,'r') %根据梯度数据画箭头 axis([-1.5,hx+2.5,-2,13]) %设置坐标边框
plot([1,1,hx,hx,1],[1,hy,hy,1,1],'k') %画导体边框 text(hx/2,0.3,'0V','fontsize',11); %下标注
%拉普拉斯方程差分式
text(hx/2-0.5,hy+0.5,'100V','fontsize',11) ;%上标注 text(-0.5,hy/2,'0V','fontsize',11); %左标注 text(hx+0.3,hy/2,'0V','fontsize',11); %右标注 hold off
这个程序运行的结果计算是：迭代次数K＝222，即导体矩形槽的节点经过222次迭代后，电位数值解收敛于某一固定值。场域内所划分的网格点电位的计算结果如表2－1所示： 表2－1 场域内网格点电位的计算结果
矩形槽内电位分布三维曲面图及槽内等位线、电场线分布图如图2－7所示：
本题是属于规则形状的第一类边值问题，原则上可以通过分离变量法求解域中的解析解，在此作为数值计算的例子可便于读者与解析解作比较，以验证该差分法的精度，请读者自行完成它们的比较。另外，查看本题的数值解容易发现，矩形场域中的电位分布是左右对称，说明我们计算的场域范围还可缩小一倍，即取矩形域左边或右边的一半进行分析和计算即可，但要学会处理第二类齐次边界的差分问题才行。 【例2－3】 设一个长直接接地金属槽，如图2－10所示，其侧壁和低面电位为零，顶盖电位为
，求槽内的电位分布。
解：在直角坐标系中，金属槽中的电位函数
满足拉普拉斯方程
其边界条件满足混合型边值问题
，
，
取步长， 、
－6
方向的网格数为
，
，共有16×10＝160个网孔，17×11＝187
个节点，其中槽内节点（电位待求点）15×9＝135个，边界节点（电位已知点）187－135＝52个。设迭代精度为10，利用MATLAB编制的计算程序如下： %例2－3：用简单迭代法求二维静电场的电位分布 hx=17;hy=11; %设置网格节点数 v1=ones(hy,hx); %设置行列二维数组 %上下两行Dirichlet边界条件 for j=1:hx
v1(hy,j)=100*sin(pi*(2*(j-1)/(hx-1)-1)); v1(1,j)=0; end
这个程序运行的结果计算是：迭代次数K＝267，即导体金属槽的节点经过267次迭代后，电位数值解收敛于某一固定值。场域内所划分的网格点电位的计算结果如表2－3所示：
表2－3 场域内网格点电位的计算结果
二维混合边界静电场域电位分布三维曲面图及域内等位线、电场线分布图如图2－11所示：
本题的区域形状虽然与前二例类似，但边界条件有关于的函数形式，比上述二例的边界条件复杂得多。另外，还含有第二类边界条件，这给编程带来了麻烦，但比起求解析解还是要简单一点。 【例2－4】 设两个同轴矩形金属槽如图2－12所示，外金属槽电位为零，内金属槽电位为求槽内的电位分布。
解：在直角坐标系中，同轴矩形金属槽中的电位函数
满足拉普拉斯方程
，
其边界条件满足第一类边界条件（Dirichlet边界 条件）问题
取步长
， ，
共有23×23＝529个网孔，设迭代精度为10，利用MATLAB编制的计算程序如下： %例2－4：用迭代法求两个同轴矩形金属槽的电位分布 FI=100;
hx=23;hy=23; %定外框网格:x轴1--23,y轴1--23
CX1=9;CX2=15;CY1=9;CY2=15; %定内框位置:x轴9--15,y轴9—15（刚好在正中间) %以上数据可以按要求更改 DX=1+CX2-CX1;DY=1+CY2-CY1;
v=ones(hy,hx)*50; %每个网格点的初值
mesh1=ones(hy,hx)*2; %每个网格点的计算标志值 %左右Dirichlet条件边界值：
v(1,:)=zeros(1,hx); %第一行的初值（下边界）
mesh1(1,:)=zeros(1,hx); %第一行标志值为零（下边界） v(hy,:)=zeros(1,hx);% 上边界初值
mesh1(hy,:)=zeros(1,hx); %上边界标志值为零
－6
， 、方向的网格数为
v(:,1)=zeros(hy,1); %左边界初值
mesh1(:,1)=zeros(hy,1); %左边界标志值为零 v(:,hx)=zeros(hy,1); %右边界初值
mesh1(:,hx)=zeros(hy,1); %右边界标志值为零 v(CY1:CY2,CX1:CX2)=ones(DX,DY)*FI; %内框区初值 mesh1(CY1:CY2,CX1:CX2)=ones(DX,DY); %内框区标志值 k=0;
difmax=1.0; %最大误差 while (difmax>1e-6) k=k+1; %计算迭代次数 difmax=0.0;
for i=2:hy-1 %从2到hy-1行循环 for j=2:hx-1 %从2到hx-1列循环 m=mesh1(i,j); %取（i,j）点的标志值
if (m>=2) %标志判断（大于等于2需要计算） vold=v(i,j); %取该点的原值
v(i,j)=(1/4)*(v(i-1,j)+v(i,j-1)+v(i+1,j)+v(i,j+1)); %拉普拉斯方程差分式 dif=v(i,j)-vold; %前后两次迭代值的差 dif=abs(dif); %取差值的绝对值
if(dif>difmax) difmax=dif; end %在所有网格中取最大误差值 end end end end
subplot(1,2,1),mesh(v) %画三维曲面图 axis([-2,hx+3,-2,hy+3,0,100])
subplot(1,2,2),contour(v,13) %画等电位线图 hold on
x=1:1:hx;y=1:1:hy;
[xx,yy]=meshgrid(x,y); %形成栅格 [Gx,Gy]=gradient(v,0.6,0.6); %计算梯度
quiver(xx,yy,Gx,Gy,-0.5,'r') %根据梯度数据画箭头 axis([-2,hx+3,-2,hy+3])
plot([1,1,hx,hx,1],[1,hy,hy,1,1],'k') %画外框边界线
plot([CX1,CX1,CX2,CX2,CX1],[CY1,CY2,CY2,CY1,CY1],'k' %画内框边界线 text(CX1+0.6,CY1+(CY2-CY1)/2,'U=100','fontsize',10); %内框文字标注 text(hx/2,hy+1,'U=0','fontsize',10); %外框上边界标注 text(hx/2,0,'U=0','fontsize',10); %外框下边界标注 text(-1.7,hy/2,'U=0','fontsize',10); %外框左边界标注 text(hx+0.4,hy/2,'U=0','fontsize',10); %外框右边界标注 hold off
这个程序运行的结果计算是：迭代次数K＝216，即同轴矩形金属槽的节点经过216次迭代后，电位数值解收敛于某一固定值。场域内所划分的网格点电位的计算结果如表2－4所示：
表2－4 场域内网格点电位的计算结果
同轴矩形金属槽内电位分布三维曲面图及槽内等位线、电场线分布图如图2－13所示：
本题属于复杂边界问题，很难求出解析解，要分析域中的电位分布，只能借助数值解。对这样的问题用有限差分法很容易处理计算，这也是数值解的优点所在。另外我们从数值结果分析可见，它存在冗余数据，读者思考如何消除这冗余数据，减小内存占用？
3.1 电磁场微分方程的泛函变分原理
物理和力学中的许多问题通常归结为满足一定边界条件或初始条件的偏微分方程。但是，在一般情况下解这些偏微分方程是非常困难的。但是对某些偏微分方程，可以构造相应的泛函，使求泛函的极小值问题与求解偏微分方程的问题等价，这就是求解物理问题中的变分原理。变分原理具有重要的理论和
实际意义，也是构造微分方程某些数值解法的基础。
泛函就是函数的函数，如函数 泛函，其中
是
函数，
，而
则称是的
是的函数。关于某物理问题变成求泛函极小
，
值问题的理解，可以从下述命题得到说明。如图3－1所示，设给定的边界条件为
，
，按自由下滑的前提，问选择什么
样的下滑曲线，即什么样的函数：
（3.1）
可使质点从就是
点下滑到
点所需时间最短。如果以
，因此，点到
表示时间，则显然
而变化的，
是不同的轨迹函数
有不同的时间
是随轨迹函数
的泛函。根据要求，质点由点所花时间最小，有：
即
（3.2）
这就是变分原理。当运动曲线满足下列参数方程时，可得到最速下降的结果：
（3.3）
其中
，
是圆滚线（运动曲线）的圆滚半径，这就是著名的伯努利命题。
一般变分法的几个主要步骤是： 1、从物理问题上建立泛函及其条件；
2、通过泛函变分，利用变分法定理得到尤拉方程； 3、求解尤拉方程。
在电磁场问题中，由于微分方程已知（除某些情况因边界条件复杂而无法求解的微分方程外），并且其微分方程可以通过对某个泛函的变分而得到。因此，这就将求解电磁场的微分方程问题化为泛函求极值的变分—积分问题，然后可通过数学离散方法得到一组代数方程，解这组代数方程，可得到近似解，这是有限元法的基本思想。 对于有多个自变量的函数
的泛函
来说，如果具有如下一般形式：
（3.4）
其中
，，。则在泛函极小的条件下，令，可得如下尤拉方程
（3.5）
对于由静电场电位
构成的函数
及
，即
（3.6）
（3.7）
（3.8）
代入（3.5）式，可得如下方程：
（3.9）
（3.10） （3.11）
（3.9）、（3.10）及（3.11）式正是电磁场的泊松方程及第二类和第一类边界条件。可见，求解（3.9）、（3.10）及（3.11）式的微分方程问题与下面（3.12）式的泛函极值等价，即求电磁场的微分方程问题可化为下述泛函求极小值的问题：
（3.12）
即
3.2 二维电磁场有限元法的数学离散形式 我们知道，如果
为电位函数，则有：
（3.12）式就是电磁场微分方程泛函变分表达式，也是电磁场数值计算的基础。
，
所以，势。因此，原理。
中的头三项乘以
， （3.13）
中的最后一项
代表单位体积的电荷的
代表静电场的能量密度，
具有能量的含义，即相当于力学中的哈密顿量，对它的体积求变分，即有如力学中的虚功
个单元（离散单元），其中，每一个单元可以有个
用坐标
的一次式近似表示（也称一阶有限元
对（3.12）式的求积分问题，可以将体积分成
节点。对于二维问题，每一个单元实际上就是一个面元，因此，它一般由三个节点组成，如图3－2所示，把场域按三角形分割，把各三角形单元内的电位法）如下：
（3.14）
因为电场为不变的。系数
，
，故（3.14）式的假定是说当面元取得足够小时，在单元内的电场可以看成是可从三角形单元的3个顶点（把它叫节点）坐标
与电位
（它们是
未知量），用下式给出：
(3.15）
从上式求出
，再代入（3.8）式则有：
（3.16）
式中，
为三角形
（
只要顺次改变
，
）
即可得到）。如用矩阵表示（3.16）式，则有 的面积。
（3.17）
式中
，
这样，单元的电位函数
可表示为：
，
（3.18）
，
，
也称为单元三个节点的形函数。对于整个区域
，将以上原则推广到整个区域，每
一个单元上的位函数都依照线性插值的原理用形函数和节点电位表示，就能解出电位在整个区域上的分布函数。总之，电位在一个单元上的分布由相应的三个节点电位和形函数的线性叠加而得，由于各个有限元互不重叠，因此，整个区域上的电位分布可近似由每个单元的电位组合而成，即
（3.19）
以上介绍了一阶二维有限元法的离散基本思路，下面将用变分法导出有限元计算的矩阵方程，即二维电磁场有限元方法的数学离散形式。对于二维电磁场的泛函表达式为
（3.20）
于是将 (3.19) 式代入 (3.20) 式可以写出相应于泊松方程的泛函：
（3.21）
这里的待定系数量的解必须使泛函系数
，
，
，
，
等即为节点上的电位值。由于
，
，
等是待求的变量，而这些变
取极小值，因此，对泛函求关于各个节点电位的导数，并使之为零，便得到待定
等变量所满足的线性方程组：
（3.22）
（3.23）
（3.24）
将（3.22）、（3.23）、（3.24）式与 (3.21) 式相比可以发现，在以上几个式子中，除了导数运算之外，积分与求和运算的次序也做了调整。上式中分别在每一个单元上对泛函求导，然后令算按每个有限元逐一进行，从而能更有效地利用计算机的数值计算和处理功能。 （3.22）、（3.23）、（3.24）式又可用矩阵方程表示：
（3.25）
式中
为
阶系数矩阵，
为
阶节点势函数矩阵，
为阶激励矩阵，为节点数。系数
个单元
的泛函导数之和为零。表面上看，这只是一个次序的调换，从实际数值计算的角度来看，这样做能使计
矩阵和激励矩阵的各元素依每个单元逐一计算，即整体矩阵的每一元素都由每个单元的贡献叠加而成，如下式所示：
（3.26） （3.27）
具有顶点
的一个任意元对矩阵系数的贡献为
（3.28）
其中各元素所在行与列由顶点由下式求出：
在整个区域中的次序所决定，式中局部系数矩阵的每一元素可
（3.29）
（3.30）
该单元对激励矩阵的贡献为
（3.31）
其中各元素所在的行由三角形的顶点在区域的次序决定，式中局部激励矩阵的每一元素可由下式求出：
（3.32）
（3.32）式可以进一步简化为
（3.33）
在有限元计算过程中，将以上局部系数矩阵和局部激励矩阵的各元素填充到整体矩阵（即系数矩阵和激励矩阵）中，其行和列的位置则由节点在分割整个区域的序号决定。最后，由（3.25）式联立（3.26）、（3.26）、（3.27）、（3.28）、（3.29）、（3.30）、（3.31）、（3.32）、（3.33）式求出场域内各点的电位。
本节在泛函变分原理基础上，介绍了二维场域有限单元法的基本概念及数值计算离散的数学建模。由于其矩阵方程编制的计算程序长、大，而又复杂，因此有限单元法的编程计算是浪费劳力的。
用MATLAB的偏微分方程工具箱进行有限元计算的平台式操作，这种方法既简单又有规可循，在二维电磁场的计算和分析中，而且其结果是比较准确的，因此可用来验证解析解或有限差分数值解的正确性。对于电磁场工程技术人员或电磁场理论学习者，掌握用MATLAB的偏微分方程工具箱进行有限元计算是非常有必要的。