第 1章数字伪像 1.1引言 事实上，解决几乎所有电磁学问题都需要一台计算机。一个理想问题即使有准确的解析解或闭式解，对于给定的一组参数，我们也必须使用计算机将解析解转变为数值解。由于数据在计算机内存储方式的固有限制，计算结果总是不可避免地具有一定的误差，虽然这些误差通常都很小，但是我们应该知道这一现象。首先，我们将讨论一些关于有限精度产生的后果。 稍后，我们将讨论基于时域有限差分（ FDTD）方法的电磁问题数值解。 FDTD方法采用的近似技术使其得到的解也是近似的，也就是说这个方法本质上是近似方法，即使我们使用计算机可以获得无限精度的数值解，用 FDTD方法得到的结果仍是一个近似解。在后面的章节中将讨论在 FDTD方法中的近似问题。 使用数值方法时需要注意一点，如果数值方法在求解过程中不发生颠覆性错误，计算机总是会给出一个结果。我们可能经常发现运行程序很简单，调试程序的过程却非常复杂。当程序能够运行时，通常情况下假设所有的错误都已被纠正了，可惜的是，事情常常并非如此。让程序运行是一回事，得到正确的结果是另一回事，事实上，得到精确的结果又是另外一回事 ——解可能只对给定的情况是正确的，但是求解过程可能得不到足够精确的结果。因此，最好是采用更多的方式检验实现过程和结果。例如，可以通过一种分辨率得到一个解，使用更高的分辨率又可能得到另一个不同的解，如果这两个解不是很接近，至少解还没有收敛到“正确的”解，就必须使用更高的分辨率，或许在解的过程中有系统错误存在。基本事实是：计算机给出一个结果，但并不意味着这个结果一定正确。 1.2有限精确度 众所周知， 3个 13之和等于 1，即 1/3+1/3+1/3=1。但是，在计算机上做这样的运算正确吗？考虑程序 1.1所示的 C程序。 在这个程序中单精度变量 a设置为 31 。在第 9行，3个 a之和与实数 1.0进行比较，如 果它们相等，程序输出 “Equal”，否则输出“Not equal”。这个程序的输出是 “Not equal”，因此，对于计算机（至少可以运行一种常用的求解电磁学问题的语言）来说， 1/3+1/3+1/3不等于 1。值得注意的是如果把第 9行写成 a=1/3，因为使用整数运算法则求解除法的值，a将会被置零。当用 a=1.0/3.0表示时，计算机就会使用浮点型法则运算。 程序 1.1 oneThird.c:测试 1/3 + 1/3 + 1/3是否等于 1。 1 /* 1.0/3.0 + 1.0/3.0 + 1.0/3.0 == 1.0? */ 2 #include 3 4 int main() { 5 float a; 6 7 a = 1.0 / 3.0; 8 9 if (a + a + a == 1.0) 10 printf("Equal.\n"); 11 else 12 printf("Not equal.\n"); 13 14 return 0; 15 } 在 C语言或 FORTRAN语言中浮点型数据类型只能存储一个数的有限位。在大部分的机器中，单精度数用 4个字节（ 32个二进制位或比特）存储，双精度数用 8个字节（ 64位）存储。再次考虑 3个 13求和问题，举一个极端的例子，假设计算机只能存储两位小数， 13约 等于 0.33。把这个数求和 3次得到 0.33 + 0.33 + 0.33 = 0.99很明显结果不等于 1。如果存储数的位数更多，结果将更接近于 1，但永远不等于 1。十进 求和 3次的结果不等 制和二进制浮点型所表示的 1 3 都有无限的位数。因此，当我们试图把 1 1 1 3 存储到计算机时 数字会被截断，所以计算机只存储 3的近似值。因此对被截断的 3 于 1。 1 如果我们对 10求和 10次，结果是多少呢？因为 1/10等于 0.1，似乎这个数可以用有限数量的位数来存储。计算机认为下面的等式成立吗？ 1111111111 +++++++++ =1? 1010 10 10 10 10 10 10 1010 可以用同样的方法在程序 1.1中进行验证。奇怪的是，结果仍然不等于 1！虽然用十进1 制表示 10时具有有限的位数，但是当用二进制表示时就有无限的位数。在浮点型十进制数中每一位都表示 10的一个特定幂数。令每一个空格表示一位，一个数可以用下面的方式理解： ... . ... 103 102 101 100 10.1 10.2 10.3 10.4 每一位可以告诉我们它对应 10的幂数是多少，小数点是负指数和非负指数的分界线。与十进制相似，二进制中的一位表示 2的幂数： ... . ... 23 22 21 20 2.1 2.2 2.3 2.4 十进制数 0.110就是 1 × 10.1。为了利用二进制数获得相同的值，必须取 2.4 +2.5 + 2.8 +2.9 + ···，即无限位数的二进制数。另一种表示方法是 0.110 =0.0001100110011001100110011 ...21 如前所述，当这个数存储在计算机中时，数字被截断，所存储的值不再正好等于 10，对这个值 10次求和就不等于 1（尽管差别非常小）。如何将浮点型值存储在计算机中不是我们关心的主要问题，然而，了解比特如何分配是有益的。数字以指数形式存储，比特的标准分配是： 总比特数 符号 尾数 指数 单精度 32 1 23 8 双精度 64 1 52 11 本质上，指数给出一个数的量级，尾数给出一个数的位数 ——尾数决定精度。尾数的位数越多，表示的数越精确。双精度的位数是单精度的两倍，它的尾数有 52位而单精度只有 23位，因此，双精度数的精度是单精度的两倍多。尾数为 23位的二进制数比 7位十进制数略小，这是因为 223是 8,388,608，因此 23位二进制数可以表示 0和 8,388,607之间的数。换句话说，尾数为 52位的二进制数所表示的值对应于十进制数 15和 16位之间 (252 =4,503,599,627,370,496)。对于指数，一个双精度数的指数位数比单精度数多 3位。看起来似乎双精度的指数变化不大，因为它的位数不是单精度的指数位数的两倍。然而，指数表示的是一个数的大小，每多一比特从本质上就表示数值加倍。如果指数的位数为 9，它可以表示两倍于单精度数的数。一个双精度数所拥有的 3位额外位数表示的指数比单精度数大 8倍，即可以转变为比单精度数大（或小） 256倍的数。考虑下面的方程： a + b = a 从数学的角度分析，我们知道当 b等于 0时这个方程才满足。然而，当使用计算机时这个方程可以成立，也就是说， b对 a没有任何贡献，即使 b为非零值。 当一些数相加或相减时，它们的尾数就会改变直到它们的指数相等。此时，尾数可以直接相加或相减。然而，如果指数间的差别大于尾数的长度，当加上或从较大的数中减去较小的数时不会有任何影响。程序 1.2中的代码说明了这一现象。 这里 a的初始值为 1.0，b设定为 0.5。变量 c为 a和 b之和。 While循环开始于第 5行，只要 c不等于 a，将一直循环下去。在循环体中， b被 2除，c再次等于 a加 b。如果计算机有无限的精度，这将会是一个无限循环。 b的值将趋向于无穷小，但是永远不会等于 0，因此 a加 b永远不可能等于 a。然而，循环会终止，在第 10行printf()打印输出为 1 5.96046e-08 1 这说明当 b的值为 5.96046 × 10.8时，a和 c是一致的。注意到 b的这个值对应于 1 × 2.24，当 b等于这个值时，a和 b指数的差别将大于 23 (a是 1 × 20)。 程序 1.2测试一个非零数 b是否满足方程 a + b = a。 1 float a = 1.0, b = 0.5, c; 2 3 c = a + b; 4 5 while(c != a) { 6 b = b / 2.0; 7 c = a + b; 8 } 9 10 printf("%12g %12g %12g\n",a,b,c); 再举一个例子来说明有限精度不以显著的方式使计算结果变差。假设变量 a等于 2，对 a取平方根，再对其值进行平方运算，得到的结果应该接近于 2（理想情况下是 2，但是由于 √ 2有无限个位数，精度上将会有一些丢失）。然而，如果对数字 2取 23次平方根后再对结果进行 23次平方会发生什么呢？我们希望得到的结果接近于 2，但事实并非如此。程序 1.3所示的程序将帮助我们验证这种情况。 程序 1.3 rootTest.c:平方根测试。 1 /*平方根测试 */ 2 #include //为了sqrt()函数 3 #include 4 5 #define COUNT 23 6 7 int main() { 8 float a = 2.0; 9 int i; 10 11 for (i = 0; i < COUNT; i++) 12 a = sqrt(a); // a的平方根 13 14 for (i = 0; i < COUNT; i++) 15 a = a * a; // a的平方 16 17 printf("%12g\n",a); 18 19 return 0; 20 } 程序的输出是 1，即结果是 a=1.0。每进行一次平方根，得到的值越来越接近 1。最后，由于截断误差，计算机认为值是 1。这时，即使对它进行再多的平方也不会改变其大小。 1.3符号处理 使用常用的数值分析语言（如 C、C++、FORTRAN或者 MATLAB）时，截断误差是不可避免的。圆周长与半径的比值为 π =3.141592 ···，这是一个无限循环的无理数，这个数不可能在计算机内存储 π的准确数值，相反地，我们必须存储一个具有有限位的近似数。然而，有些软件（例如 Mathematica）允许对数学符号运算。在 Mathematica软件中，你可以将 π写成 Pi，Mathematica软件知道这个数学符号代表什么。例如， 10000000001*Pi的余弦值等于 .1。相似地，也可以写成 Sqrt[2]，Mathematica软件知道它的平方恒等于 2。只可惜这种符号运算在计算资源方面是非常昂贵的，许多电磁学的前沿问题包含成千上万甚至上百万的未知量。为了处理这些大量的数据，数据运算在内存和时间上必须尽可能地高效。 Mathematica软件虽然是非常好的工具，但它并不是解决大型数值问题的有效工具。 在 MATLAB中我们可以把pi看作 π的缩写。然而，这里所表示的 π和 Mathematica软件中用到的并不一样。在 MATLAB中，pi在本质上只是数值表示 ——它只是比写成完整的 pi方便。在 C语言中，将会提供包含 math.h的头文件，可以用 M PI作为 π的缩写。在 math.h中可以看到下面的语句： # define M_PI 3.14159265358979323846 /* pi */ 这恰好和 MATLAB中的一样， MATLAB只知道 pi的数值，而且这只是正确值的一个截断值。因此，将会用 10000000001*pi的余弦值等于 .0.99999999999954来代替准确值 .1（当然在这种情况下误差是微不足道的）。 图 2.1电荷之间的电场力电荷 Q2受到电荷 Q1的力沿两电荷间的直线方向，力的方向由两电荷的正负决定。电场指向径向远离正电荷的方向，如用指向远离 Q1的直线表示（这里假设 Q1是正电荷） 电场力与电荷量成正比，与电荷之间距离的平方成反比。电荷 Q2受到电荷 Q1的作用力可以表示为 1 Q1Q2F12 = a.12 R2 (2.1). 4πE0 12 *本书的符号用法保持英文原书的风格。 式中 a.12为从电荷 Q1指向电荷 Q2的单位矢量 , R12是电荷间的距离， 1/4πE0是常数。 E0是真空的介电常数，大小约为 8.854 × 10.12 F/m。电荷的单位是库仑（ C），电荷有正负之分。当两电荷电性相同时，两电荷间相互排斥： F12的方向与 a.12的方向相同。当两电荷电性相反时，F12的方向与 a.12的方向相反。 方程 (2.1)有一个不足，因为它意味着力作用于一段距离。这个方程似乎告诉我们相互作用力 F12的建立是即时的，从这个方程可以假定电荷之间距离 R12的改变将导致作用力 F12的瞬间改变，但事实并非如此。一个电荷位置的变化传达到另一个位置的电荷需要有限的时间间隔（类似地，电荷量的变化传达到另一个电荷也要花有限的时间）。为了克服这个不足，利用场的概念解释是很方便的。我们不说 Q1直接作用力于 Q2上，而是理解为 Q1产生了一个场，这个场产生了一个作用于 Q2的力。 Q1产生的场不依赖于 Q2 ——无论 Q2是否存在并感受到它，Q1的场都一直存在。 在静态情况下，场的方法并不比直接利用库仑定律更有优势。这是因为对于静电荷，库仑定律是正确的。为了能够区别，场是时变的。不过，与时变情况下一致，场同样适用于静态情况。点电荷 Q1产生的电场强度为 ! = .(2.2) E1 ar Q12 4πE0r其中 a.r为径向远离点电荷的单位向量， r为空间中任意一点到点电荷的距离。电场强度的单位为伏特每米 (V/m)。 为了得到作用于 Q2上的力，仅需要将其电荷乘以电场强度： F12 = Q2E1。一般地，任意电荷所受的力为该电荷与该点电场强度的乘积，也就是 F = QE。 2.3电位移矢量 所有物质都是由带电粒子组成的，物质可能不显示电性是因为它带有等量的正负电荷。然而，在包括电场影响在内的多种情况下，物质内部的正负电荷可能发生轻微偏移，由此造成的电荷分离将对整个电场产生影响。正因为如此，引入电位移矢量 D是非常方便的，电位移矢量 D的单位是库仑每平方米 (C/m2)，本质上 D忽略了约束在物质内部电荷局部影响。在自由空间中，电场强度和电位移矢量的关系是 D = E0E (2.3)高斯定律表明 D在闭合曲面的积分等于该闭合曲面包围的自由电荷量，即 D · ds = Qenc (2.4) S 其中 S是闭合曲面 , ds是一个面元增量且它的法线垂直向外， Qenc是闭合曲面包围电荷的电荷量。例如，考虑方程 (2.2)给定的电场强度，假定 S是一个球面且电荷在其中心，在方程 (2.4)中进行简单积分： π2πj j ! Q1 D · ds = E0 a.r · a.rr 2 sin θ dφ dθ = Q1 (2.5) 4πE0r2 Sθ=0 φ=0 实际上其结果并不依赖于所选择的积分曲面（假设它包围所有的电荷），但是对于球面来说，积分运算特别简单。 我们希望方程 (2.4)中的积分结果和在自由空间中一样总是等于包围的电荷量。然而，当物体存在时，事情就变得比较复杂。如图 2.2所示，考虑具有等量且均匀电荷分布的两个大平行板，但电荷的符号相反。 图 2.2自由空间中的带电平行板虚线表示积分表面 S 虚线表示一个距离平板边缘足够远的积分表面 S，使得电场在表面 S上均匀分布。设电场强度恒等于 E0，平行板外的电场强度为 0，平行板内的电场强度和 S的侧边相切。因此对积分有贡献的只是 S的顶部。积分的结果 :S E0E · ds是由表面 S包围的负电荷（也就是说，负电荷位于表面 S内部平行板的底部）。 现在考虑同样的平行板，带有同样的电荷，但是在平行板内充满介质。假设这个介质是“可极化的”，介质内的正负电荷可以略微偏移，但电荷不能完全自由移动 ——它们是被束缚的电荷。正电荷被上面的板排斥，被下面的板吸引。相反地，负电荷被下面的板排斥，被上面的板吸引，如图 2.3所示。 图 2.3在带电平行板之间有一个可极化介质 细长的椭圆表示分子，它们的定向电荷在上板的下面产生一个负电荷层，在下板的上面产生一个正电荷层。在介 质内部，正负束缚电荷相互抵消。在物体的表面必须考虑束缚电荷。因此，在整个图中，分子没有全部画出来。如 图 2.3右边所示，仅仅束缚电荷层被表示出来。平行板上的自由电荷形成的电场强度为 E0。束缚电荷形成了和 E0反向的电场强度 Em，因此使总电场强度变小。虚线仍然表示积分表面 S 尽管介质的存在，平行板上电荷所产生的电场强度仍然是 E0，也就是和图 2.2中的电场强度一样。然而，在两平行板之间由于极化介质上束缚电荷的偏移存在另外一个电场。 极化材料在与底板相邻处建立了一个正电荷层，在与顶板相邻处建立了一个负电荷层。由这些电荷层产生的电场仍然是均匀的，只是方向与板上的自由电荷产生的电场方向相反。在图 2.3中由束缚电荷产生的电场强度为 Em。则总的电场强度是束缚电荷和自由电荷所产生电场强度之和，也就是， E = E0 + Em。由于 E0和 Em是反向平行的，所以总电场强度 E将小于 E0。 由于介质不显电性，我们希望在面 S上电通量的积分等于底板上被包围的电荷量 ——不是介质内的束缚电荷。在某种意义上，这就意味着积分表面不能把任何一个分子内的正负束缚电荷分开。每一个分子要么完全在积分表面的内部，要么完全在积分表面的外部。由于每一个分子都不带电，所以对积分的贡献来自于板上的自由电荷。 由于介质的存在， :S E0E · ds的积分等于非常小的电荷量。这是因为，如上所述，总电场强度 E比没有介质存在时要小。为了纠正减少的电场强度并获得理想的结果，重新定义电位移矢量以便它能够把介质的影响考虑在内。电位移矢量更一般的表达式为 D = ErE0E = EE (2.6) 其中 Er是相对介电常数， E称为介电常数。通过把介质的介电常数考虑在内，高斯定律就能总是成立。在式 (2.6)中，D和 E通过一个标量常数相互联系在一起，这就意味着 D和 E对于所有频率、方向、场强，通过一个简单的比例常数联系在一起，可惜的是现实世界不是这么简单。很明显，如果电场强度足够大，就可以把束缚的正负电荷分开。因为电荷有一定的质量，它们对所有的频率并不具有同样的反应方式。此外，很多介质有一定的结构，比如晶体，在不同方向上的性质不一样。虽然如此，高斯定律既然是定律，它总是成立的。当事情变得复杂时，我们应该舍弃简单的介电常数，换一种适当的介电常数形式确保高斯定律的成立，例如用一个方向相关的张量介电常数。然而，除了与频率相关的特性（即色散媒质）外，我们不会考虑那些复杂的情况，一个标量介电常数就足够了。 2.4静电场 忽略材料可能的非线性特征，叠加原理适用于电磁场。因此，可以认为电荷分布是点电荷的集合。我们可以通过求得所有电荷作用之和来得到总场（如果电荷是连续分布的，其和可以用积分形式获得）。从方程 (2.2)中可以看到点电荷的电场放射性地从点电荷指向四周，没有场的“涡旋”。当存在不止一个电荷时，总场可能会弯曲，但是它不会涡旋。想象一个非常小的圆盘，正电荷沿它的圆周分布。圆盘的中心固定在一个位置，但是圆盘可以自由地以中心旋转。对于静电场，无论我们把圆盘放在什么地方，都没有合力作用于它使其旋转。在一个特定的方向上如果有一个合力作用于整个圆盘（一个平移的力），但是推动圆盘顺时针旋转的力与推动圆盘逆时针旋转的力相平衡。静电场的另一个性质是电位移矢量只开始或终止于自由电荷。如果没有电荷存在，那么电通量线将延续下去。没有涡旋电场和电位移矢量源是十分简单的概念。然而，为了能够正确分析场，我们需要一个关于这些概念的数学陈述。正确的表达式为 飞× E = 0 (2.7) 和 飞· D = ρv (2.8) 其中飞读作 “del”，是哈密顿算子， ρv是电荷密度（单位是 C/m3）。方程 (2.7)表示电场的旋度，方程 (2.8)表示电位移矢量的散度。这两个方程会在以后的部分做进一步讨论。 2.5梯度、散度和旋度 哈密顿算子独立于使用的坐标系 ——自然地，场的特征不取决于描述场的坐标系。然而，哈密顿算子可以表示在不同的坐标系中。在直角坐标系中飞表示为 飞≡ a.x . + a.y . + a.z . (2.9) .x .y .z 其中符号 ≡的意思是“定义为”。飞作用在标量场上产生该场的梯度。假设 f是一个标量场，飞f产生的矢量场表示为.f .f .f 飞f = a.x .x + a.y .y + a.z .z (2.10) 梯度 f的方向沿数值变化率最大的方向，其大小和变化率成正比。假定我们希望求解 f当 dx在 x方向上一个小变化时的变数量，可以通过飞f · a.xdx获得，即 飞f · a.xdx = .f dx =(x方向上的变化率) × (x方向上的移动) (2.11) .x 这可以推广到任意方向上的变化。定义一个非常小的长度增量 d. = a.xdx + a.ydy + a.zdz (2.12) 场的变化通过移动一个小量 d.来实现 飞f · d. = .f dx + .f dy + .f dz (2.13) .x .y .z 回到方程 (2.8)，当飞点乘一个矢量场时，就得到这个场的散度。散度可以看作为在给定一位置上点“源”或“汇点”场强的测量。一个矢量场的散度是一个标量，即 飞· D = .Dx + .Dy + .Dz (2.14) .x .y .z 在 xy面上考虑这一散度的有限差分近似，如图 2.4所示。 这里的散度是通过一个小盒子估算的，并且假定场在盒子每条边的值为常数。导数用中心差分近似： .Dx .Dy 2 到 . Dy 到(x + Δx ,y) . Dx (x . Δx ,y) Dy (x, y + Δ2y (x, y . Δ2y Dx 2 + ≈ + .x .y Δx Δy (2.15)