数值法求解静态电磁场边值问题

2015-03-27电子科技大学物理电子学院

电子世界 2015年18期

电子科技大学物理电子学院薛冰

1 引言

在科学和工程中的许多问题，常常归结为求解方程或方程组。由于计算量大，大多数这类问题都不能采用手工计算，而需要利用计算机，采用数值方法进行计算[1-2]。而其中线性方程组求解不但在科学和工程中常常涉及到，而且数值计算方法其它分支的一些研究也往往归结为此类问题[3]。求解线性方程组的数值方法主要有两类：直接法和迭代法。在没有舍入误差的时候，采用直接法可以精确求得方程组的解。这类方法主要有高斯消去法、三角分解法、平方根法和追赶法等[3]。迭代法采用极限过程来逐步逼近准确解，包括雅克比迭代法、高斯-赛德尔迭代法及超松弛迭代法等。它对计算机内存要求低，更适合于求解高阶方程组。

本文采用三角分解法、高斯-赛德尔迭代法及超松弛迭代法，针对用有限差分法求解静态电磁场边值问题得到的差分方程组进行求解。采用不同网格步长离散化待求区域，得到规模不同的差分方程组。采用不同数值方法求解方程组，可以发现不同方法有不同的运算效率及适用范围。

2 基于有限差分法的静态电磁场边值问题

有限差分法是一种求解微分方程的数值方法，由于其逻辑清晰，方法简单，自上世纪五十年代以来得到了广泛应用[4]。为求解由微分方程定解问题所构造的数学模型，有限差分法将定解区域（场区）离散化，并以各离散点上函数的差商来近似该点的导数，使待求的微分方程转化为差分方程组。求解差分方程组，即可得到各离散点处的待求函数值。

下面我们首先采用一个简单的静态电磁场边值问题，来说明有限差分法的求解过程[5-6]。然后用数值法求解得到的方程组。

图1所示的一个长直接地金属矩形槽，其横截面为正方形D：，其侧壁和地板均接地，电位为0，顶盖与侧壁绝缘，电位为V0，求该区域中的电位分布。

由于该槽沿长度方向为均匀分布，故可以将其简化成二维电磁场问题。区域中无电荷源，故该区域中的电位符合二维拉普拉斯分布。设电位为，则有：

图1 静态电磁场边值问题网格划分

要采用有限差分法求解，首先将该区域离散化，沿x和y方向均匀划分网格。因为求解区域为正方形，故两个方向的网格步长都设为h，每边的网格数均为N+1个。区域内的电位用网格点上的电位来表示，其中分别表示网格在x和y方向的序号，如图1所示。

(1)式可转化成如下方程组[7]：

3 三角分解法

3.1 计算原理

三角分解法，又称为LU分解法。对线性方程组，其中系数矩阵A可以分解为两个矩阵的乘积，其中L是单位下三角矩阵，U是上三角矩阵。

根据矩阵乘法，可以依次求得矩阵U的第k行和L的第k列。原方程组就转化三角方程组：

最后通过回代过程可以方便地求出方程组的解。

3.2 编程思路和计算结果

采用MATLAB编程求解此问题[8-9]。设步长为h，则方程组的系数矩阵为N×N(N=L/h-1)阶的方阵。当h较小时，网格数很多，矩阵规模就会很庞大。若采用普通的方法来存贮系数矩阵，则需要N×N阶的数组空间。本问题中，(2)式左边的系数矩阵是一个稀疏阵，非零元素很少且为带状排列。为节省内存资源，我们采用对角存储法中的等带宽存贮法来存放矩阵元素[10-11]。因为对带状稀疏矩阵进行LU分解，不会在带状结构之外引入非零填入元，因此求解算法不会有任何改变，而且非零元集中存储，还可以提高计算效率[10]。

对角线存储法将各非零对角线作为数组单元进行存储。对于具有d条非零对角线的N×N稀疏矩阵A，可用两个数组表示：一个是N×d的值数组E，它的每一列存储了矩阵A的一条对角线元素，列数d为矩阵A的非零对角线的数量；另一个是1×d的偏移数组t=[l1,l2,…,ld]，它依次存储值数组E中每一列所对应的对角线相对于主对角线的偏移量，其中在主对角线下方的为负值，上方的为正值。

假设此问题中，正方形槽的边长L=1m，采用笔记本电脑（联想电脑，CPU为酷睿i 5，主频2.6GHz，内存2GB）进行计算。选取3种步长，分别是h=0.1m，0.05m和0.01m，将计算结果列于表1。由于每次计算的运行时间有少量偏差，故所列数值为5次计算的均值。而当网格步长为h=0.01m时，此时步长最小，网格数量为99×99个，方程组阶数较高，此电脑无法运算，故没有得到结果。

图2给出了网格步长为0.1m和0.05m时的电位分布图。很显然，此结果符合槽中的电位分布情况。

表1 直接法运行时间

图2 采用直接法计算得到的网格电位分布图，

4 迭代法

对于线性方程组Ax=b，设其同解方程组为x=Bx+f，若初始解为x(0)，其迭代格式可以表示为：

4.1 高斯-赛德尔迭代法

对于线性方程组Ax=b，令A=L+D+U，其中，L是矩阵A对角线以下元素构成的严格下三角矩阵，U是矩阵A对角线以上元素构成的严格上三角矩阵，D是A的对角元构成的矩阵。则可构造迭代格式为：

其中k为迭代次数，此即高斯-赛德尔迭代格式。

对于式(2)，可构造如下的高斯-赛德尔迭代格式：

进行迭代，就可计算出待求的网格电位。

4.2 超松弛尔迭代法

把高斯－赛德尔法的迭代值作为中间值与加权求平均，就得到超松弛迭代格式：

4.3 迭代法计算结果

应用MATLAB编程，用与直接法相同的电脑求解。选择同样的3种步长，并设迭代收敛的精度要求为。采用超松弛迭代时，根据计算得到的松弛因子为。用高斯-赛德尔迭代法和超松弛迭代法计算的迭代次数和时间（5次运行时间的均值）如表2所示。当网格步长为0.01m时，采用迭代法能够比较快速地完成计算。采用细网格，可以更准确地模拟空间的电位分布情况，只是运算时间会更长。

为了研究松弛因子对迭代收敛速度的影响，在步长为h=0.1m时，还比较了时的计算时间迭代收敛时间分别为112.3315(s)，11.07(s)和68.0965(s)，可见松弛因子的取值对收敛速度影响很大。同时也可发现，在的最优值附近取值时，超松弛迭代法都比高斯-赛德尔迭代法的收敛速度更快。

图4 步长为0.01m时的电位分布图

5 比较与结论

本文采用直接法和迭代法对一个静态电磁场边值问题进行求解。从求解情况可知，当矩阵规模较小时，采用直接法比迭代法的运算效率更高。但当矩阵规模超过102时，直接法的计算速度明显降低，且如果矩阵进一步增大，在普通电脑上甚至可能无法计算。而两种迭代法的求解结果表明，采用超松弛迭代法改进了高斯-赛德尔迭代法的收敛速度，选取合适的松弛因子更有利于提高收敛速度。