第四章 有限差分方法 4.1 引言 有限差分法：数值求解常微分方程或偏微分方程的方法。 物理学和其他学科领域的许多问题在被分析研究之后, 往往可以归结为常微分方程或偏微分方 程的求解问题。一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同 时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。 有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。 在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离 散取值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的 连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值插值计算来近似得 到。这种方法是随着计算机的诞生和应用而发展起来的。其计算格式和程序的设计都比较直观和简单， 因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分。 1 有限差分法的具体操作分为两个部分： (1)用差分代替微分方程中的微分， 将连续变化的变量离散化，从而得到差分方程组的数学形式； (2)求解差分方程组。 在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。通 常采用的是规则的分割方式。这样可以便于计算机自动实现和减少计算的复杂性。网络线划分的交点 称为节点。若与某个节点P相邻的节点都是定义在场域内的节点，则P点称为正则节点；反之，若节 点P有处在定义域外的相邻节点，则P点称为非正则节点。 在第二步中，数值求解的关键就是要应 用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。 有限差分法的差分格式： 一个函数在x点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。 如对一个单变量函数f(x)，x为定义在区间[a,b]的连续变量。以步长h= Δx将[a,b]区间离散化，我 们得到一系列节点x = a , x = x + h , x = x + h = a + 2 1 2 1 3 2 Δx , ..., x = x + h = b , 然 后求出 f(x)在这些点上的近似值。显然步长h越小，近似解的精度就越好。与节点 相邻的节点有 和 n+1 n i x hx i ? hx i + ，因此在 点可以构造如下形式的差值： i x ),()( ii xfhxf ?+ 节点 的一阶向前差分 x i ),()( hxfxf ii ?? 节点 的一阶向后差分 x i ).()( hxfhxf ii ??+ 节点 的一阶中心差分 x i 2 与 点相邻两点的泰勒展开式可以写为 x i fx h fx hf x h fx h fx h fx iii i i i ( ) () () () ! () ! ( ) ...... ' '' ''' '''' ?=?+?+? 234 23 4 . (4.1.1) fx h fx hf x h fx h fx h fx iiii i i ( ) () () () ! () ! ( ) ...... ' '' ''' '''' += + + + + + 23 4 . (4.1.2) （4.1.1）-(4.1.2)，并忽略h的平方和更高阶的项得到一阶微分的中心差商表示： fx fx h fx h h i ii' () ()() ≈ +? ? 2 . (4.1.3) 利用(4.1.1)和(4.1.2)式我们还可以得到一阶微分的向前，向后一阶差商表示 fx fx h fx h i ii' () ()() ≈ +? , (4.1.4) fx fx fx h h i ii' () () ( ) ≈ ?? . (4.1.5) 将(4.1.1)和(4.1.2)式相加，忽略h的立方及更高阶的项得到二阶微分的中心差商表示 fx fx h fx fx h h i iii'' () ()()() ≈ +? + ?2 2 . (4.1.6) 3 利用(4.1.3) ~ (4.1.6)式，我们就可以构造出微分方程的差分格式。这里要指出的是：在构造 差分格式时，究竟应该选择向前，向后还是中间差分或差商来代替微分方程中的微分或微商，应当根 据由此得到的差分方程解的稳定性和收敛性来考虑。同时兼顾到差分格式的简单和求解的方便。 上述差分步骤应用于偏微分： 例如，对于 的情况，拉普拉斯算符在0点作用在此函数上的值 ffxy= (,) 22 2 22 f f f xy ?? ?? ???? ? ，也可 以用临近的点上的函数值来表示出来。(见图4.1.1, 且 ?= + ??? ???? hhhhh ==== 4321 时) ).( !4 24 4 4 4 42 2 043212 y f x fh h fffff f ? ? ? ? +? ?+++ ≈? (4.1.7) j-1 j i+1 j-1 (i-1, j-1) (i, j-1) h1 h3 4 h4 (i+1, j-1) (i+1, j) (i, j) j(i-1, j) 30 h2 2 (i+1, j+1)(i, j+1) j+1 (i-1, j+1) 图4.1.1节点0及邻近节点 4 对微分方程数值求解的误差的来源： （1）方法误差（或截断误差）。这是由于采用的计算方法所引起的误差。例如上面我们介绍的差商 表示中，采用的泰勒展开式展开到第 1+n 项时的截断误差阶数为 ( ) 1+n hO 。具体方法的误差阶数取决于 在离散化时的近似阶数。因此若改进算法就可以减小截断误差。 （2）舍入误差（或计算误差）。这是由于计算机的有限字长而造成数据在计算机中的表示出现误差。 在计算机运算的过程中，随着运算次数的增加舍入误差会积累得很大。如果在多次运算后，舍入误差 的精度影响是有限的，那么这个算法是稳定的，否则是不稳定的。不稳定的算法是不能用的。 本书中我们将略去对差分法稳定性和收敛性理论的讨论，尽管这方面的内容是相当重要的。以下 的讨论中所讲到的各种差分格式，我们均假定求解方法满足稳定性和收敛性的要求。 5 4.2 有限差分法和偏微分方程 利用上节所介绍的微分的差分表示，我们就很容易地将微分方程离散化为差分方程组的形式。但 是由差分方程所得的解完全取决于待求微分方程的特性。正如我们在物理上所知道的，边界条件的情 况变化将会引起差分方程组的不同。在求解微分方程中，我们会遇到两类问题：一类是初始值问题； 另一类是边值条件的问题。 在初始值问题中，部分边界上的函数值和部分的函数偏导值是给定的。通 常在这类问题中的独立变量之一是时间t。在边界值问题中，边界上的信息是给定的。本书中我们仅 讨论后一类问题。 假定某方程形式上可以写为： Lqφ = . (4.2.1) 其中 L为含偏微商的算符． 它的边界条件一般可写为： ).(|)(| 21 sg n sg GG =+ ? ?φ φ (4.2.2) G 表示场域 D的边界， 为边界上s点的逐点函数。 )(),( 21 sgsg 6 三类边界条件： （1）第一类边界条件，或称为狄利克莱(Dirichlet)问题 (, )gg 12 00= ≠ 。 ).(| sg G =φ (4.2.3) （2）第二类边界条件，或称诺伊曼(Neumann)问题 (, )gg 12 00≠ = 。 ).(| sg n G = ? ?φ (4.2.4) （3）第三类边界条件，或称混合问题 (, )gg 12 00≠ ≠ 。 对于算符 L为斯杜-刘维尔(Sturm-Liouville)算符的特定情况，即 Lpf≡?? ? +() . (4.2.5) 公式中的p 和f是给定的函数。我们将会得到一类很重要的微分方程。它是在流体力学、等离子物 理、天体物理…等学科中，势函数起关键作用的许多问题当中的基本方程。当p=1, f=0时，我们得 到（4.2.1）式的特殊情况—泊松(Poisson)方程。 我们现在考虑方程(4.2.1)中p为常数的二维的情况，我们可以得到下面的方程： ?φ ? ?φ ? φ 2 2 2 2 xy fxy qxy++ =(,) (,). (4.2.6) 设函数 φ 在区域 D内满足方程（4.2.6）式(区域 的边界为G)。 D 7 区域 的离散化： D 即通过任意的网络划分方法把区域 离散为许许多多的小单元。原则上讲这种网格分割是可以任 意的，但是在实际应用中，常常是根据边界G的形状，采用最简单，最有规律，和边界的拟合程度最 佳的方法来分割。常用的有正方形分割法和矩形分割法（如图4.2.1）。有时也用三角形分割法（见 图4.2.2）。对圆形区域，应用图（4.2.3）所示的极网络格式也许更方便些。这些网络单元通常称 为元素，网络点称为节点。 D D 0 x y D 0 x y 4.2.1 求解区域的矩形分割。 4.2.2 求解区域的三角形分割。 4.2.3 求解区域的极网络分割。 8 用节点上的函数值来表示节点上偏导的数值。 设区域内部某节点0附近的各节点如图4.1.1所示。这里我们取步长h不相等的最一般情况。以 φ φ φ φ φ 01234 分别代表在节点0，1，2，3，4处 ,,,, φ 的函数值。如前所述，0点的一阶偏导数可以通过先 前或向后的差商，由1和3节点近似写出 ?φ ? φφ xh ? ? ? ? ? ? ≈ ? 0 10 1 . （4.2.7） 或 ?φ ? φφ xh ? ? ? ? ? ? ≈ ? 0 03 3 . （4.2.8） 显然这种单侧差商的误差较大。 如果要寻求更精确的差分格式，我们可以引入待定常数 α β, ，由 φ 1 和 φ 3 的泰勒展开，构造出如下的 关系式： ....)( 2 1 )()()( 2 3 2 1 0 2 2 31 0 0301 ++ ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? =?+? hh x hh x βα ? φ? βα ? ?φ φφβφφα （4.2.9） 令 ?φ ? 2 2 0 x ? ? ? ? ? ? 项的系数为零，则得到 α和 β 之间应当满足 αβ=? h h 3 2 1 2 . （4.2.10） 9 将公式（4.2.10）带入（4.2.9），并舍去高阶项，得到 ?φ ?x ? ? ? ? ? ? 0 的另一个差分表达式 ?φ ? φφ φφ x hh hh h h ? ? ? ? ? ? ≈ ?? ? + 0 3 2 10 1 2 30 13 1 3 ()() () . （4.2.11） 当选用等步距 时，上式成为 hhh x13 == ?φ ? φφ xh x ? ? ? ? ? ? ≈ ? 0 13 2 . （4.2.12） （此为中心差商表达式。） 二阶偏导数的差分表达式（“五点格式”或“菱形格式”） 在（4.2.9）式中，如果令 ?φ ?x ? ? ? ? ? ? 0 的系数为零，则有 α 和 β 间存在关系式： αβ= h h 3 1 . （4.2.13） 将上式（4.2.13）代入（4.2.9）式中，并忽略h三阶以上的高次项，则得到表达式： ?φ ? φφ φφ 2 2 0 31 0 13 0 13 1 3 2 x hh hh h h ? ? ? ? ? ? ≈ ?+ ? + ()() () . （4.2.14） 当用等步长 hhh x13 = = 时，上式成为 10 ?φ ? φφφ 2 2 0 103 （4.2.15） 2 2 x h x ? ? ? ? ? ? ≈ ?+ . 用完全相同的计算方法，我们可以推导出 ?φ ? 2 2 0 y ? ? ? ? ? ? 的差分表达式： ?φ ? φφ φφ 2 2 0 42 0 24 0 24 2 4 2 y hh hh h h ? ? ? ? ? ? ≈ ?+ ? + ()() () . （4.2.16） 当采用等步长 时， 有 hhh y24 == ?φ ? φφφ 2 2 0 204 2 2 y h y ? ? ? ? ? ? ≈ ?+ . （4.2.17） 将公式（4.2.14）和（4.2.16）两式代入方程（4.2.13），我们就得到该方程的差分表达式为 000 4242 042024 3131 031013 0 2 )( )()( )( )()( 2)( qf hhhh hh hhhh hh =+ ? ? ? ? ? ? + ?+? + + ?+? =? φ φφφφφφφφ φ （4.2.18） 如果在x和y方向的步长分别相等， 即 hhh x13 = = 和 hhh y24 = 时，则上式化为 = φφφφφφ φ 103 2 204 2 00 0 22?+ + ?+ += hh fq xy ， （4.2.19） 一般可以用角标来表示节点的标记，将上式写为 11 1 2 1 2 2 11 2 11 hh fq x ij ij ij y ij ij ij ij ij ij ()() ,, , , ,, ,,, ， (4.2.20） φφφ φφφ φ +?+? ?+ + ?+ + = 特别是当 hhh xy = = 时，我们得到： ,, （4.2.21） φφφφ φ i j i j ij ij ij ij ij hf hq +? +? ++++ ?= 11 11 22 4 ,,,, , , () f = 0 的时候方程（4.2.6）为泊松方程，由（4.2.21）式得到 , ， （4.2.22） φφφφ φ i j i j ij ij ij ij hq +? +? ++ = 11 11 2 4 ,,,, , fq==0的时为拉普拉斯方程，从（4.2.21）式得 φ φ φ φ φ i j i j ij ij ij+? +? + + + ? = 11 11 40 ,,,, , ， （4.2.23） 边界条件的离散化的处理： 若场域的网络节点都落在边界G上，则显然无需再做处理。但是在一般情况下，边界G是不规则 的。网络节点不可能全部都落在边界G上。对(4.2.3)式给出的第一类边界条件，通常有两种处理办 法。 一种是所谓的直接转移法，如果0节点靠近边界，则取最靠近0点的边界节点上的函数作为0 点的函数值。这是一种比较粗糙的近似。另一种方法是较为精确的线性插值法。对第二、三类边界条 件也可以用插值法求出临近边界节点上的函数值。 12 （1）第一类边界条件 在二维情况下，第一类边界条件如公式(4.2.3)所示。如果网格的边界节点恰好落在边界G上， 则显然无需再做近似处理，边界节点的函数值就等于边界条件(4.2.3)给出的函数值。但是一般情况 下网格边界节点不在边界上，我们通常用以下三种方法处理。 （a） 直接转移法 在图(4.2.4)中网格是按正方形分割，步长为 。0点为靠近边界G的一个网格节点，1和2为边 界节点。我们取最靠近0点的边界节点1上的函数值作为0点的函数值。即取 h 10 φφ ≈ 。这种方法称为 直接转移法，又称为零次插值法。 图(4.2.4)场域的第一类边界条件 13 （b） 线性插值法 如图(4.2.4)所示，先判断 x方向的边界节点1和 y 方向的边界节点2哪一个更靠近0点。如果1 更靠近0点，则可以用 x方向的线性插值给出0点的函数值： 1 311 0 hh hh + + = φφ φ (4.2.24) 如果边界节点2更靠近0点，则可以类似地用 y 方向2，4边界节点的函数值 2 φ 和 4 φ 的线性插值求出 0点函数值 0 φ 。 2 422 0 hh hh + + = φφ φ ， (4.2.25) 其中 和 分别为边界节点1，2到0点的距离。这种方法的误差为 1 h 2 h ( ) 2 hO 。 （c） 双向插值法（更为精确的方法） 如果 hh α= 1 ， hh β= 2 ，将它们代入公式(4.2.15)和(4.2.17)（注意：这里 hhh yx == ），由方程 (4.2.6)得到0点附近的差分计算格式： () () 0 2 00 2 04321 22 11 1 1 1 1 1 1 1 1 q h f h =+ ? ? ? ? ? ? ? ? +? + + + + + + + φφ βα φ β φ α φ ββ φ αα , (4.2.26) 14 （2） 第二类和第三类边界条件 图(4.2.5)场域的混合边界条件 第二类和第三类边界条件可以统一写为 .| g n G = ? ? ? ? ? ? +αφ ? ?φ (4.2.27) 其中 0=α 时，即为第二类边界条件； 0≠α 时，即为第三类边界条件。 一种比较简单的处理方法： 如图(4.2.5)所示，过 点向边界 做垂线 交边界于 点，交网线段 于O G PQ Q VR P ，假定 ，OP PR和 的长度分别为 ah ， 和 ，则对 点有 VP bh ch O 15 ()hO nah O P + ? ? ? ? ? ? ? ? = ? φφφ 0 . (4.2.28) 因为 点一般不是节点，其值应当以 点和 P V R点的插值给出。 ( ) 2 hOcb RVP ++= φφφ . (4.2.29) 将其代入公式(4.2.27)，并由于 ()hO nn QO + ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? φφ ，则得到 ()()hO n cb ah Q RV + ? ? ? ? ? ? ? ? =?? φ φφφ 0 1 . (4.2.30) 从边界条件(4.2.27)，有关系式 () ()QgQ n Q Q +?= ? ? ? ? ? ? ? ? φα φ . (4.2.31) 结合上面两个公式(4.2.30)和(4.2.31)，并将 Q φ 用 O φ 来近似，就得到 点的差分计算公式 O ()()()QgQcb ah ORVO =+?? φαφφφ 1 . (4.2.32) 16 若 n? ?φ 的方向与网格线平行： 如 n? ?φ 与 x 方向平行，设图(4.2.5)中 点与O R点重合，则（4.2.28）成为 hxn VO OO φφφφ ? ≈ ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? . (4.2.33) 得到 点的差分计算公式为 O ( ) () ()QgQ h OVO =+? φαφφ 1 . (4.2.34) 如 n? ?φ 与 方向平行，设图(4.2.5)中 点与 y V R点重合，则（4.2.28）成为 hyn RO O O φφφφ ? ≈ ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? . (4.2.35) 得到 点的差分计算公式为 O ()()()QgQ h ORO =+? φαφφ 1 . (4.2.36) 17 4.3有限差分方程组的迭代解法 求解微分方程(4.2.6)。假定该问题是个在边界G上的狄里克莱问题。其求解的区域 D是个单位矩 形区间（ )1,0 ≤≤ yx 。我们在平行于 x y, 轴的方向分别用 和 )1( +N )1( +M 个点以等步长作网络划分，边界 G上的节点函数值为 j （如图4.3.1所示）。则用 ) i g , ()(N M+ +11网格划分的单位矩形求解区间 D中， x y, 方向的步长分别是 h N= 1/ 和 k M=1/ 。对这样的问题利用差分计算格式（公式 （4.2.21） ），并取 （即 ），则方程（4.2.6）可以近似写为 NM = kh = [ ] φφφφ φ i j i j i j i j ij ij ij hf hq +? +? +++?? = 11 11 22 4 ,,,, ， 在区域 内， D φ ij ij g ,, = ， 在 D的边界G 上。 ( 4.3.1) x y (n-1,1) (n-1,2) (2,1) (1,1) (1,2) y=h y=2h y=(M+1)h x=h x=2h x=(n-1)h 图4.3.1用 )()(N M+ +11 网格划分的（4.2.6）的求解范围。 18 我们现在考虑 ff ij = =00) 的特殊情况，此时要求解的是狄里克莱边界问题的泊松方程，它可以写 成为 ( , ),( 2 2 2 2 yxq yx =+ ? φ? ? φ? ， 在 D内， φ|(). G gp= 在 的边界G 上 (4.3.2) D 微分方程问题(4.3.2)对应的差分方程组为（参见公式（4.3.1）） ijjijijijiij q h 4 )( 4 1 2 1,1,,1,1 ?=+++? ?+?+ φφφφφ ， 在 D内， φ ij ij g ,, = ， 在 D的边界G 上。 (4.3.3) 引入y方向的层向量（也可以取x方向分层的层向量） ， φ φ φ φ j j j Nj = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 2 1 , , , M ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= ? jN j j j qh qh qh q ,1 2 ,2 2 ,1 2 4 1 M . 19 并记 ， . Φ= ? ? ? ? ? ? ? ? ? ? ? ? ? φ φ φ 1 2 1 M N ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ?1 2 1 N b b b B M 则方程组（4.3.3）就可以写为 B=ΚΦ . （4.3.4） 其中K矩阵的形式如 . （4.3.5） ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? =Κ GI I IGI IG 4/0 4/ 4/4/ 04/ L OOM M M L L I是 阶的单位矩阵。G为 阶方阵，其具体表示为 2 )1( ?N 2 )1( ?N . （4.3.6） ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = 14/10 4/1 14/1 04/11 L OOM MO L G 20 从公式（4.3.3）—（4.3.6），我们可以得到 yh= 上各个节点的差分方程有如下形式 1 0,11,1,1 2 0,21,2 2 0,11,01,1 2 2,1 2,2 2,1 1,1 1,2 1,1 4 1 100 0 10 001 4 1 14/10 4/1 14/1 04/11 b ggqh gqh ggqh NNNNN = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ?? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? MM L OOM MO L MOOM MO L φ φ φ φ φ φ （4.3.7） 即 121 4 1 bIG =? φφ . （4.3.8） 同样沿 hy2 上的各节点可以列出差分方程为 = 2321 4 1 4 1 bIGI =?+? φφφ . （4.3.9） 其中 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? 2,2,1 2 2,2 2 2,02,1 2 2 4 1 NN gqh qh gqh b M . （4.3.10） 将求解微分方程（4.3.2）变成为求解（4.3.4）的线性代数方程。 21 计算结果的收敛性检验： 原则上，只要我们取网格间距足够小，这个方程可以得到精确解。但是我们注意到差分方程的解 不大可能与原来的偏微分方程（4.3.2）的解完全相同。两者间的偏差正是由于用差分公式代替偏微 分所带来的。在实践中，有必要在求解方程（4.3.4）时采用不同的网格间距h值来计算以检验结果 的收敛性。 求解方程（4.3.4）有各种各样的方法，下面我们将列出三种方法。 （1） 采用直接求解法来解方程（4.3.4）： 这是通过求系数矩阵K的逆矩阵来得到方程的解 . （4.3.11） Φ= ? KB 1 在实际运用中，由于矩阵 K的维数通常都很大，且计算机计算的舍入误差会引起数值结果的不稳定， 因而在实际应用中还存在较大困难。但对象泊松方程这类问题的求解，就可以采用直接法。泊松方程 用富里叶变换改写后，采用直接求解法计算会非常快。这部分内容我们将在本章第四节中介绍。 （2） 用随机游动：我们已在§2.6节中对此做了详细介绍。 22 （3） 迭代求解法： 在采用有限差分法求解微分方程的实践中得到广泛应用的方法，其中尤以超松弛迭代法的使用效 果最佳。迭代求解法实际上是一种极限方法，它被用来求方程组的近似解。本节我们仅详细介绍迭代 求解法。 差分方程组（4.3.4）中的系数矩阵K具有如下特征： 系数矩阵 K是一个仅有少数不为零的元素的大型稀疏矩阵。 K矩阵的每一行元素中只有少数几个不 为零的事实可以从上面我们给出的五点格式中看出：它的每一行元素中非零元素的个数不超过5个。 因而在程序中只需记存系数矩阵 K中的非零元素及每个非零元素在此一维数组中地址码的信息便足 够了。这样就可以极大节省计算机的存贮空间。当边界与网格节点重合时，矩阵 K是个对称正定矩阵， 其非零元素都是实数。可以证明这时的矩阵 K有 ( ) 个正交本征向量，其对应的本征值为 2 1?N ()hqhpK pq ππλ coscos 2 1 1)( +?= ， )1,...,2,1,( ?= Nqp . (4.3.12) 但是当边界与网格节点不重合时， K的对称性将被破坏。 K矩阵通常是不可约的，因而该方程组不能 由其中的某一部分单独求解。 23 求解方程组（4.3.4）的迭代求解法： 求解差分方程组(4.3.4)的最简单的办法是雅可比方法， 又称直接迭代法。 我们可以将公式(4.3.4) 等价地写成如下形式 BR +Φ=Φ , (4.3.13) 其中 KIR ?= ， I 为单位矩阵。 如果已经得到一组势函数估计值 )1( ? Φ k ，则我们可以通过如下公式得到“改进”后的估计值 )(k Φ 。即 BR kk +Φ=Φ ? )1()( ， i j k jij k i BR += ∑ ? )1()( φφ . (4.3.14) 如果 K矩阵为实对称矩阵，则 R 矩阵也应当是实对称矩阵。 R 也称为雅可比(Jacobi)迭代矩阵。这个 矩阵的一个重要特征是它有一个特定的谱半径 )(Rρ ，此谱半径等于该矩阵最大本征值的绝对值。 从 K矩阵的本征值表示式(4.3.12)可知矩阵 R 的本征值为 ()hqhpR pq ππλ coscos 2 1 )( += ， )1,...,2,1,( ?= Nqp . (4.3.15) 24 若我们开始时给出一组任意的猜测势函数值 ，我们可以证明在连续使用公式(4.3.13) 的迭代后，会得到的改进值收敛到满足差分方程组(4.3.3)或(4.3.4)的解 }{ )0()0( i φ=Φ }{ i φ=Φ 。设在第 次迭代后 的一组误差为 ，初始时误差为 ，我们有 k i k i k i E φ?Φ= )()( iii E φ?Φ= )0()0( )1()1()1()()( )()( ??? =Φ?Φ=+Φ?+Φ=Φ?Φ= kkkkk RERBRBRE ， (4.3.16) 即 )0()( ERE kk = . k R 应当在 ∞→k 时收敛到零矩阵，迭代得到的值才收敛到满足微分方程(4.3.4)的解，并且与初始选 择的 无关。数学上， 的条件是 )0( Φ 0→ k R R 矩阵的谱半径应满足不等式 1)( <Rρ 。由(4.3.5)和(4.3.6) 式可以看出 R 矩阵的每一行元素之和均小于1，而只要有一行元素之和小于1，则可保证 矩阵具 有满足不等式 R 1)( <Rρ 的特性。 因此在连续使用公式(4.3.13)的迭代后可以得到与初始估计值 )0( Φ 无 关的满足差分方程组(4.3.4)的解。 定义误差矢量和迭代矩阵的范数为 ? ? ? ? ? ? = ∑ = N i k i k eE 1 2 )()( ， ? ? ? ? ? ? ? ? = ∑ j ij i RR max . (4.3.18) 谱范数 R 表示当矩阵 R 作用在单位矢量上时被放大的最大因子。则从公式(4.3.16)可导出 )1()( ? ≤ kk ERE . (4.3.19) 25 当 足够大时,从(4.3.19)可以写为(参考文献[2]) k )1()( )( ~ ? < kk ERE ρ . (4.3.20) 因此矩阵 R 的谱半径越接近1，则迭代收敛速度越慢。 直接迭代方法的实际计算步骤： 我们将直接迭代式(4.3.14)具体写成为 )( 4 1 , 2)( 1, )( 1, )( ,1 )( ,1 )1( , ji k ji k ji k ji k ji k ji qh?+++= ?+?+ + φφφφφ . (4.3.21) （1） 首先任意给出各内节点处的初始函数值 ，然后代入上面方程（4.3.21）的右端，求出各内 节点的第一次迭代法的函数近似值 。 )0( , ji φ φ ij, ()1 （2） 然后依次循环下去，以第 次迭代法的近似值来求出k 1+k 次的近似值。 这种所谓松弛迭代法的缺点是： 它需要两套存贮单元，分别存贮两次相邻次数迭代的函数近似值，因 而需占用的内存较大。此外, 该方法的收敛速度也慢, 因此它也就没有什么实用价值。 26 高斯—赛德尔迭代法：一种改进后比较好的迭代方法。 1+k 它实际上是在 次迭代中，将已经得到的某些相关节点上的第 1+k 次迭代近似值代入进行计算。 具体地说，如果在沿 方向（或 y x 方向）求得了 ( ) yjh （或 =?1 ( )xi h=?1) 层的 1+k 次迭代值，则可在 （或 yjh= xih= ）层节点的计算中，将临近的点上已有的 1+k 次迭代值代入进行计算。这实际上 就是将直接迭代公式(4.3.21)做如下的改写： )( 4 1 , 2)1( 1, )1( ,1 )( 1, )( ,1 )1( , ji k ji k ji k ji k ji k ji qh?+++= + ? + ?++ + φφφφφ . （4.3.22） 这就是高斯—赛德尔迭代式。 与迭代式对应的矩阵形式， . (4.3.23) BUL kkk GS +Φ+Φ=Φ ++ )()1()1( 其中矩阵间存在关系式 RKIUL =?=+ 。 矩阵的所有对角线和对角线以上的元素为零，而 矩阵则 有所有对角线和对角线以下的元素为零。该迭代方法的误差矢量是以下面公式逐步减小的。 L U )1(2)( )( ? = kk ERE ρ . (4.3.24) 27 超松弛迭代法： 高斯—赛德尔迭代方法在起始阶段的收敛速度可能比简单的直接迭代法快些，但仍然不是很理 想。为了加快收敛速度，通常引入一个松弛因子 ω ，而把用（4.3.22）式迭代计算出的结果作为一 个中间结果，它表示为 )( 4 1 , 2)1( 1, )1( ,1 )( 1, )( ,1 )1( , ji k ji k ji k ji k ji k ji qh?+++= + ? + ?++ + φφφφφ . （4.3.25） 取 次迭代的最后近似值为 1+k )( , k ji φ 和 的加权平均，即 )( , k ji φ )( , )1( , )( , )1( , )( , )1( , )1()( k ji k ji k ji k ji k ji k ji φωφωφφωφφ ?+=?+= +++ )4( 4 )( ,, 2)1( 1, )1( ,1 )( 1, )( ,1 )( , k jiji k ji k ji k ji k ji k ji qh φφφφφ ω φ ??++++= + ? + ?++ . （4.3.26） 利用公式(4.3.23)我们得到 , （4.3.27） BLIR kkk GS kk GS kk 1)1()1()()1()()1()( )()1()( ????? ?+Φ=Φ?+Φ=Φ?Φ+Φ=Φ ωωωωω ω 其中矩阵 ， 的定义同前。迭代矩阵 的表示式为 L U ω R . （4.3.28） ])1([)( 1 IULIR ωωω ω ?+?= ? (4.3.27)式与（4.3.13）式相似。 28 类似与前面由公式(4.3.13)导出(4.3.16)式,我们可以由(4.3.26)推出 . (4.3.29) k )1()( ? = kk ERE ω 当 足够大时, 从(4.3.29)可以得到 )1()( )( ~ ? < kk ERE ω ρ . (4.3.30) )( ω ρ R 为矩阵 的谱半径。因此，超松弛迭代法的收敛速度决定于松弛因子 ω R ω 的选取。 松弛因子的取值： 由(4.3.30)知，松弛因子 ω 的值的选择标准应当是它能减小矩阵 的最大本征值的数值。 ω R ω 的取值 范围在 21 ≤≤ω 时，收敛速度较好。当 1=ω 时，这就是高斯—赛德尔迭代法。一般情况下确定 ω 的最 佳值 0 ω ，只能靠经验来选取。 对于正方形区域的第一类边值问题，最佳的 ω 可从理论上选为 ω π 0 2 1 = + sin( / )l . （4.3.31） l +1为每边的节点数。若是矩形区域，用正方形网格分割，每边的节点数分别为 l + 1和 1+m ，则可选取 ωπ 0 22 22 11 =? +() l m . （4.3.32) 一般地讲，只要超松弛因子 ω选得合适，就可以大大地加快收敛速度，可以做到有阶的改善。 29 迭代精度判断： 由于我们事先不可能知道满足方程(4.3.4)的函数值 {} i φ=Φ ，要决定迭代值是否满足精度要求 ii k i εφφφ <? )( ，我们可以采用判断不等式 ε< Φ Δ )( )( k k (4.3.33) 是否满足。 如果满足，则迭代可以结束。其中 ε 为设定的相对精度值，移位矢量 )(k Δ 的定义为 {} }{ )1()()()( ? ?≡Δ=Δ k i k i k i k φφ . (4.3.34) ε 一般将 值取得比要求的精度要小10倍以上，以保证实际所要求的精度。 30 4.4 求解泊松方程的直接法 利用迭代法来求解差分方程组往往计算量非常大，尤其会涉及大量的反复计算。当然，我们可以 采用上一节中介绍的一些有效的办法来加快有限差分法数值求解的收敛速度。但是，即使是这样，我 们仍然感觉到这样来加快求解速度的效果是有限的。 为避免采用前面介绍的反复计算，并加快有限差分法数值求解的收敛速度，我们将会用到直接法 求解。本节我们介绍Hockney于1970年提出的基于有限傅立叶级数展开和循环相消法的直接求解法 [3]。 我们以一个非常简单的第一类边界条件情况下的泊松方程的求解问题来说明。该问题的数学形式 为 () ?φ ? ?φ ? 2 2 2 2 xy qxy+=, ， Dyx ∈),( , 0| = G φ ， . (4.4.1) )( 的边界为 DG 其中场域 D为一个边长为 L 的正方形，在边界G上的势函数值均为零。 31 采用离散傅立叶级数展开的方法，将任意周期性函数 ( ()xf ( ) ( )mLxfxf += ， ,...),2,1( =m L为周期） 做离散傅立叶级数展开的形式为 () () () ∑ ∞ = ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? += 1 0 2 sin 2 cos 2 j jj L jx xB L jx xA A xf ππ . (4.4.2) 对正方形场域 D内任意确定的x值，将位势函数 ( )φ xy, 做以 为变量的有限傅立叶级数展开，则得到 与式(4.4.2)对应的有限傅立叶级数为 y () () () ∑ = ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? += 2/ 1 0 2 sin 2 cos 2 , n j jj L jy xB L jy xA A yx ππ φ . (4.4.3) 这里取 ynxnL Δ=Δ= 。 考虑到后面将采用快速傅立叶变换，我们取 （其中 为正整数）。又注意到边界上位势为零的 特殊情况，上式中傅立叶级数的系数为 k n 2= k ， ， () 0 0 =xA () 0=xA j () () dy L jy yx L xB L j ? ? ? ? ? ? = ∫ π φ 2 sin, 2 0 ， )2/,...2,1( nj = . (4.4.4) 则公式（4.4.3)就写为 () () ∑ = ? ? ? ? ? ? = 2/ 1 2 sin, n j j L jy xByx π φ . (4.4.5) 32 把(4.4.5)式代入泊松方程 (4.4.1)，并利用三角函数的正交特性，则得到系数 ( ) Bx j 满足的方程为 () () ()xCxB L j xB jjj =?′′ 2 22 4π )2/,...,2,1( nj = (4.4.6) C j 是电荷分布 对确定的 ()yxq , x ，以 y 为变量的傅立叶级数展开的系数。它等于 () () dy L jy yxq L xC L j ? ? ? ? ? ? = ∫ π2 sin, 2 0 . (4.4.7) 通过有限傅立叶级数展开，已经把求解满足泊松方程(4.4.1)的势函数 ( ) φ xy, 值的问题变成了计算与 之对应的系数 值（在此问题中 ()Bx j ( ) 0=xA j ）。 采用差分法来解微分方程(4.4.6)。 首先对该正方形场域 采用正方形分割法离散化，即选择 方向等步长 的网格划分，使得 D yx, h nhL= 。 然后再利用公式(4.2.15) 将方程(4.4.6)改变为差分方程。这就是将任意一个网络的内节点 xih i = 上 系数函数 ( )xB j 的微商值，用临近两节点的函数值来表示。 为表示方便，我们记系数 ( ) ij xB 为 。与(4.4.6)式对应的差分方程为 )(i j B ， )(2)1()()1( i j i j i jj i j ChBBB =+? +? λ )1(,...,2,1( ?= ni ， )2/,...,2,1 nj = . (4.4.8) 33 其中我们定义了 ? ? ? ? ? ? ? ? += 2 4 2 222 L hj j π λ , (4.4.9) ∑ = ? ? ? ? ? ? = n m mi i j n mj q n C 0 , )( 2 sin 2 π . (4.4.10) 由于边界上位势为零的条件要求(见(4.4.1)式)，在边界节点上的系数 和 都必须为零。采用直 接法对差分方程组(4.4.8)求解，也就是在确定的 层的网格节点上，计算 个系数值 （ ）。 B j ()0 B j n() jhy j = 2/n )(i j B 2/,...,2,1 nj = 通过这些系数值，利用公式(4.4.5)就可以构造出计算内节点上势函数值的公式 () ? ? ? ? ? ? = ∑ = n jk Byx n j i jki π φ 2 sin, 2/ 1 )( . (4.4.11) 因此，其求解步骤分为三步： (1)对电荷分布函数做离散傅立叶变换，见公式(4.4.10)； (2)求解方程组(4.4.8)中的n(n-1)/2个系数 ； )(i j B (3)利用 的数值做公式(4.4.11)下的傅立叶分析，就得到位势函数的解。 )(i j B 34 循环相消方法： 对于上述三步求解步骤，计算量一般还是很大的。但是我们注意到在这个问题中我们可以采用一些数 学技巧。由于我们取 （其中 为正整数），因而我们可以在第一和第三步中采用快速傅立叶变 换的方法。此外方程组(4.4.8)还可以采用循环相消的方法来求解。下面我们以一个例子来说明循环 相消方法的求解步骤。 k n 2= k 为了叙述方便，我们将方程(4.4.1)的场域 D作进一步的简化。我们在x和y方向上分别将正方 形场域 D划分为仅只有 个小区间，即 nhL 88× = ， )8( =n 。对任意的 )8,...,2,1( =jj ， 我们可以将方程组(4.4.8) 中的三个含i点的函数值的方程写出 )1(2)()1()2( ??? =+? i j i j i jj i j ChBBB λ )(2)1()()1( i j i j i jj i j ChBBB =+? +? λ ， )1(2)2()1()( +++ =+? i j i j i jj i j ChBBB λ )8,...,2,1( =i （4.4.12） 35 公式(4.4.12)已将x方向节点上的系数值联系在一起。将 j λ 乘上上面方程组中的第二式，再与另外 两个公式相加，我们就得到 ′ =+′? +? )(2)2()()2( i j i j i jj i j ChBBB λ ， （ i=2,4,6） (4.4.13) 其中我们定义 . ,2 2 ?≡ ′ jj λλ )1()()1()( +? ++≡ ′ i j i jj i j i j CCCC λ 按照上面一样的处理办法，我们又可以得到 ″ =+ ″ ? +? )(2)4()()4( i j i j i jj i j ChBBB λ ， （i=4） (4.4.14) 即 ″ =+ ″ ? )4(2)8()4()0( jjjjj ChBBB λ . (4.4.15) 公式(4.4.12)和(4.4.13)中又定义 ,2)( 2 ? ′ ≡ ″ jj λλ ′ + ′ ′ + ′ ≡ ″ +? )1()()1()( i j i jj i j i j CCCC λ 36 (4.4.15)式说明： （1） 可以用已知的边界节点系数值 ， 和由公式(4.4.10)所示的系数 来表示，从而得到 ； )4( j B )0( j B )8( j B )(i j C )4( j B （2）再取公式(4.4.11)中i的分别为2和6，则 和 就也可以用已知的数值表示。 )2( j B )6( j B （3）最后，我们再应用所有已经得到的数值代入方程(4.4.12)，就可以计算出剩下的系数 ， ， 和 。 )1( j B )3( j B )5( j B )7( j B 这样对于边界上势函数为零的泊松方程求解问题就完全解决了。 循环相消方法的操作就是对差分方程组(4.4.8)做连续相消，直到只剩下少量可以求解的方程组。 当然，这种方法的应用范围也是有限的，它只适合于具有特殊类型的边界条件情况，并且在实际应用 中，其运算往往也比上面所举的简单例子中要复杂得多。但是在某些特殊情况下，尽管计算复杂些， 但是仍然比超松弛迭代法要快。 37