有限差分方法

第四章 有限差分方法 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 + 21 2 1 3 2 Δ x , ..., x = x + h = b , 然 后求出 f(x)在这些点上的近似值。显然步长 h越小，近似解的精度就越好。与节点 相邻的节点有 和 n+1 n ix hxi − hxi + ，因此在 点可以构造如下形式的差值： ix ),()( ii xfhxf −+ 节点 的一阶向前差分 xi ),()( hxfxf ii −− 节点 的一阶向后差分 xi ).()( hxfhxf ii −−+ 节点 的一阶中心差分 xi 2 与 点相邻两点的泰勒展开式可以写为 xi f x h f x hf x h f x h f x h f xi i i i i i( ) ( ) ( ) ( ) ! ( ) ! ( ) ......' '' ''' ''''− = − + − + − 2 3 4 2 3 4 . (4.1.1) f x h f x hf x h f x h f x h f xi i i i i i( ) ( ) ( ) ( ) ! ( ) ! ( ) ...... ' '' ''' ''''+ = + + + + + 2 3 4 2 3 4 . (4.1.2) （4.1.1）-(4.1.2)，并忽略 h的平方和更高阶的项得到一阶微分的中心差商表示： f x f x h f x h hi i i' ( ) ( ) ( )≈ + − − 2 . (4.1.3) 利用(4.1.1)和(4.1.2)式我们还可以得到一阶微分的向前，向后一阶差商表示 f x f x h f xhi i i' ( ) ( ) ( )≈ + − , (4.1.4) f x f x f x h hi i i' ( ) ( ) ( )≈ − − . (4.1.5) 将(4.1.1)和(4.1.2)式相加，忽略 h的立方及更高阶的项得到二阶微分的中心差商表示 f x f x h f x f x h hi i i i'' ( ) ( ) ( ) ( )≈ + − + −2 2 . (4.1.6) 3 利用(4.1.3) ~ (4.1.6)式，我们就可以构造出微分方程的差分格式。这里要指出的是：在构造 差分格式时，究竟应该选择向前，向后还是中间差分或差商来代替微分方程中的微分或微商，应当根 据由此得到的差分方程解的稳定性和收敛性来考虑。同时兼顾到差分格式的简单和求解的方便。 上述差分步骤应用于偏微分： 例如，对于 的情况，拉普拉斯算符在 0点作用在此函数上的值f f x y= ( , ) 2 22 2 2 f ff x y ∂ ∂ ∂ ∂ ⎛ ⎞⎛ ⎞⎟ ，也可 以用临近的点上的函数值来表示出来。(见图 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) 3 0 1 h2 2 (i+1, j+1)(i, j+1) j+1 (i-1, j+1) 图 4.1.1 节点 0及邻近节点 4 对微分方程数值求解的误差的来源： （1）方法误差（或截断误差）。这是由于采用的计算方法所引起的误差。例如上面我们介绍的差商 表示中，采用的泰勒展开式展开到第 1+n 项时的截断误差阶数为 ( )1+nhO 。具体方法的误差阶数取决于 在离散化时的近似阶数。因此若改进算法就可以减小截断误差。 （2）舍入误差（或计算误差）。这是由于计算机的有限字长而造成数据在计算机中的表示出现误差。 在计算机运算的过程中，随着运算次数的增加舍入误差会积累得很大。如果在多次运算后，舍入误差 的精度影响是有限的，那么这个算法是稳定的，否则是不稳定的。不稳定的算法是不能用的。 本书中我们将略去对差分法稳定性和收敛性理论的讨论，尽管这方面的内容是相当重要的。以下 的讨论中所讲到的各种差分格式，我们均假定求解方法满足稳定性和收敛性的要求。 5 4.2 有限差分法和偏微分方程 利用上节所介绍的微分的差分表示，我们就很容易地将微分方程离散化为差分方程组的形式。但 是由差分方程所得的解完全取决于待求微分方程的特性。正如我们在物理上所知道的，边界条件的情 况变化将会引起差分方程组的不同。在求解微分方程中，我们会遇到两类问题：一类是初始值问题； 另一类是边值条件的问题。在初始值问题中，部分边界上的函数值和部分的函数偏导值是给定的。通 常在这类问题中的独立变量之一是时间 t。在边界值问题中，边界上的信息是给定的。本书中我们仅 讨论后一类问题。 假定某方程形式上可以写为： L qφ = . (4.2.1) 其中 L为含偏微商的算符． 它的边界条件一般可写为： ).(|)(| 21 sgnsg GG =+ ∂ ∂φφ (4.2.2) G 表示场域D的边界， 为边界上 s点的逐点函数。 )(),( 21 sgsg 6 三类边界条件： （1）第一类边界条件，或称为狄利克莱(Dirichlet)问题 ( , )g g1 20 0= ≠ 。 ).(| sgG=φ (4.2.3) （2）第二类边界条件，或称诺伊曼(Neumann)问题 ( , )g g1 20 0≠ = 。 ).(| sgn G =∂ ∂φ (4.2.4) （3）第三类边界条件，或称混合问题 ( , )g g1 20 0≠ ≠ 。 对于算符 L为斯杜-刘维尔(Sturm-Liouville)算符的特定情况，即 L p f≡ −∇ ∇ +( ) . (4.2.5) 公式中的 p 和 f 是给定的函数。我们将会得到一类很重要的微分方程。它是在流体力学、等离子物 理、天体物理…等学科中，势函数起关键作用的许多问题当中的基本方程。当 p=1, f=0 时，我们得 到（4.2.1）式的特殊情况—泊松(Poisson)方程。 我们现在考虑方程(4.2.1)中 p 为常数的二维的情况，我们可以得到下面的方程： ∂ φ ∂ ∂ φ ∂ φ 2 2 2 2x y f x y q x y+ + =( , ) ( , ). (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不相等的最一般情况。以 φ φ φ φ φ0 1 2 3 4 分别代表在节点 0，1，2，3，4处, , , , φ的函数值。如前所述，0点的一阶偏导数可以通过先 前或向后的差商，由 1和 3节点近似写出 ∂φ∂ φ φ x h ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − 0 1 0 1 . （4.2.7） 或 ∂φ∂ φ φ x h ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − 0 0 3 3 . （4.2.8） 显然这种单侧差商的误差较大。 如果要寻求更精确的差分格式，我们可以引入待定常数α β, ，由φ1 和φ3 的泰勒展开，构造出如下的 关系式： ....)( 2 1)()()( 23 2 1 0 2 2 31 0 0301 ++⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛+−⎟⎠ ⎞⎜⎝ ⎛=−+− hh x hh x βα∂ φ∂βα∂ ∂φφφβφφα （4.2.9） 令 ∂ φ∂ 2 2 0x ⎛ ⎝⎜ ⎞ ⎠⎟ 项的系数为零，则得到α和 β 之间应当满足 α β= − h h 3 2 1 2 . （4.2.10） 9 将公式（4.2.10）带入（4.2.9），并舍去高阶项，得到 ∂φ∂x ⎛ ⎝⎜ ⎞ ⎠⎟ 0 的另一个差分表达式 ∂φ∂ φ φ φ φ x h h h h h h ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − − − +0 3 2 1 0 1 2 3 0 1 3 1 3 ( ) ( ) ( ) . （4.2.11） 当选用等步距 时，上式成为 h h hx1 3= = ∂φ∂ φ φ x hx ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − 0 1 3 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 3 1 0 1 3 0 1 3 1 3 2 x h h h h h h ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − + − + ( ) ( ) ( ) . （4.2.14） 当用等步长 h h hx1 3= = 时，上式成为 10 ∂ φ∂ φ φ φ2 2 0 1 0 3 （4.2.15） 22x hx ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − + . 用完全相同的计算方法，我们可以推导出 ∂ φ∂ 2 2 0y ⎛ ⎝⎜ ⎞ ⎠⎟ 的差分表达式： ∂ φ ∂ φ φ φ φ2 2 0 4 2 0 2 4 0 2 4 2 4 2 y h h h h h h ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − + − + ( ) ( ) ( ) . （4.2.16） 当采用等步长 时， 有 h h hy2 4= = ∂ φ∂ φ φ φ2 2 0 2 0 4 2 2 y hy ⎛ ⎝⎜ ⎞ ⎠⎟ ≈ − + . （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方向的步长分别相等， 即 h h hx1 3= = 和 hh h y2 4 = 时，则上式化为 = φ φ φ φ φ φ φ1 0 32 2 0 42 0 0 02 2− + + − + + =h h f qx y ， （4.2.19） 一般可以用角标来表示节点的标记，将上式写为 11 1 2 1 22 1 1 2 1 1h h f q x i j i j i j y i j i j i j i j i j i j( ) ( ), , , , , , , , , ， (4.2.20） φ φ φ φ φ φ φ+ − + −− + + − + + = 特别是当h h hx y= = 时，我们得到： , , （4.2.21） φ φ φ φ φi j i j i j i j i j i j i jh f h q+ − + −+ + + + − =1 1 1 1 2 24, , , , , ,( ) f = 0的时候方程（4.2.6）为泊松方程，由（4.2.21）式得到 , ， （4.2.22） φ φ φ φ φi j i j i j i j i j i jh q+ − + −+ + + − =1 1 1 1 24, , , , , f q= = 0的时为拉普拉斯方程，从（4.2.21）式得 φ φ φ φ φi j i j i j i j i j+ − + −+ + + − =1 1 1 1 4 0, , , , , ， （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点的距离。这种方法的误差为1h 2h ( )2hO 。 （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 qhfh =+⎟⎟⎠ ⎞⎜⎜⎝ ⎛ +−+++++++ φφβαφβφαφββφαα , (4.2.26) 14 （2） 第二类和第三类边界条件 图(4.2.5)场域的混合边界条件 第二类和第三类边界条件可以统一写为 .| gn 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点的插值给出。 ( )2hOcb RVP ++= φφφ . (4.2.29) 将其代入公式(4.2.27)，并由于 ( )hO nn QO +⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂=⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ φφ ，则得到 ( ) ( )hOncbah QRV +⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂=−− φφφφ01 . (4.2.30) 从边界条件(4.2.27)，有关系式 ( ) ( )QgQn QQ +−=⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ φαφ . (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 ( ) ( ) ( )QgQh OVO =+− φαφφ 1 . (4.2.34) 如 n∂ ∂φ 与 方向平行，设图(4.2.5)中 点与y V R点重合，则（4.2.28）成为 hyn RO OO φφφφ −≈⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂=⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ . (4.2.35) 得到 点的差分计算公式为 O ( ) ( ) ( )QgQh 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 所示）。则用 )ig , ( )(N M+ +1 1 网格划分的单位矩形求解区间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 ijh f h q+ − + −+ + + − − =1 1 1 1 2 24, , , , ， 在区域 内， D φ i j i jg, ,= ， 在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+ +1 1 网格划分的（4.2.6）的求解范围。 18 我们现在考虑 f fi j= =0 0)的特殊情况，此时要求解的是狄里克莱边界问题的泊松方程，它可以写 成为 ( , ),(2 2 2 2 yxq yx =+ ∂ φ∂ ∂ φ∂ ， 在D内， φ| ( ).G g p= 在 的边界 G上 (4.3.2) D 微分方程问题(4.3.2)对应的差分方程组为（参见公式（4.3.1）） ijjijijijiij qh4)(4 1 2 1,1,,1,1 −=+++− −+−+ φφφφφ ， 在D内， φ i j i jg, ,= ， 在D的边界 G上。 (4.3.3) 引入 y方向的层向量（也可以取 x方向分层的层向量） ， φ φ φ φ j j j N j = ⎛ ⎝ ⎜⎜⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟⎟⎟− 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 Nb 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 ML 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），我们可以得到 y h= 上各个节点的差分方程有如下形式 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） 同样沿 hy 2 上的各节点可以列出差分方程为 = 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） Φ = −K B1 在实际运用中，由于矩阵K的维数通常都很大，且计算机计算的舍入误差会引起数值结果的不稳定， 因而在实际应用中还存在较大困难。但对象泊松方程这类问题的求解，就可以采用直接法。泊松方程 用富里叶变换改写后，采用直接求解法计算会非常快。这部分内容我们将在本章第四节中介绍。 （2） 用随机游动：我们已在§2.6 节中对此做了详细介绍。 22 （3） 迭代求解法： 在采用有限差分法求解微分方程的实践中得到广泛应用的方法，其中尤以超松弛迭代法的使用效 果最佳。迭代求解法实际上是一种极限方法，它被用来求方程组的近似解。本节我们仅详细介绍迭代 求解法。 差分方程组（4.3.4）中的系数矩阵 K具有如下特征： 系数矩阵K是一个仅有少数不为零的元素的大型稀疏矩阵。K矩阵的每一行元素中只有少数几个不 为零的事实可以从上面我们给出的五点格式中看出：它的每一行元素中非零元素的个数不超过 5个。 因而在程序中只需记存系数矩阵K中的非零元素及每个非零元素在此一维数组中地址码的信息便足 够了。这样就可以极大节省计算机的存贮空间。当边界与网格节点重合时，矩阵K是个对称正定矩阵， 其非零元素都是实数。可以证明这时的矩阵K有 ( ) 个正交本征向量，其对应的本征值为 21−N ( )hqhpKpq ππλ coscos2 11)( +−= ， )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 的本征值为 ( )hqhpRpq ππλ coscos2 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 iE φ−Φ= )()( iiiE φ−Φ= )0()0( )1()1()1()()( )()( −−− =Φ−Φ=+Φ−+Φ=Φ−Φ= kkkkk RERBRBRE ， (4.3.16) 即 )0()( ERE kk = . kR 应当在 ∞→k 时收敛到零矩阵，迭代得到的值才收敛到满足微分方程(4.3.4)的解，并且与初始选 择的 无关。数学上， 的条件是)0(Φ 0→kR R 矩阵的谱半径应满足不等式 1)(