1、将 G 剖分为网格区域，见图 9-1 。h1，h2分别称为 x 方向和 y 方向的剖分步长，网格 交点（xi，yi）称为剖分节点 （区域内节点集合记为 Gh=（xi，yi）; （xi，yi）G ）, 网格线与边界 的交点称为边界点，边界点集合记为 h。现在将微分方程（ 9.1 ）在每一个内节点（ xi ，yi）上进行离散。在节点（ xi， yi）处，方程（ 9.1 ）为 u2（xi,yi） u2（xi,yi） f （xi , yi ）, （xi,yi） Gh （9.5）需进一步离散（ 9.5 ）中的二阶偏导数。为简化记号，简记节点（ xi ， yi）=（ i, j） ，节2u2 （i, j）
2、x2 （i, j） y点函数值 u（xi，yi）=u（i, j ） 。利用一元函数的 Taylor 展开公式，推得二阶偏导数的 差商表达式12 u（i 1, j） 2u（i, j） u（i 1, j） 0（h12） h1212 u（i, j 1） 2u（i, j） u（i, j 1） 0（h22） h2代入（ 9.5 ）式中，得到方程（ 9.1 ）在节点 （ i , j ）处的离散形式1112u（i 1, j） 2u（i,j） u（i 1,j） 12u（i,j 1） 2u（i,j） u（i,j 1） h1 h2fi,j 0（h12 h22）, （i, j） Gh其中 fi,j f （xi, y
3、i ） 。舍去高阶小项 0（h12 h22 ） ，就导出了 u（i,j）的近似值 ui,j所满足的差分方程1h12ui 1,j 2ui,j ui 1,j h122ui,j 1 2ui,j ui,j 1 fi,j, （i,j） Gh （9.6）在节点 （i,j）处方程（ 9.6 ）逼近偏微分方程（ 9.1 ）的误差为 O（h12 h22） ，它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的 截断误差 ，它的大小将 影响近似解的精度。在差分方程（ 9.6 ）中，每一个节点 （i, j）处的方程仅涉及五个节点未知量 ui,j， ui+1,j，ui-1,j，ui,j+1，ui,j-1，因此
4、通常称（ 9.6 ）式为五点差分格式 ，当 h1= h2=h时， 它简化为h2 ui 1,j ui 1,j ui,j 1 ui,j 1 4ui,j fi,j, （i, j） Gh差分方程（ 9.6 ）中，方程个数等于内节点总数，但未知量除内节点值 ui,j ,（ i, j） Gh外，还包括边界点值。例如，点 （1, j ）处方程就含有边界点未知量 u0, j。因此，还要利用给定的边值条件补充上边界点未知量的方程。对于第一边值条件式（ 9.2 ），可直接取 对于第三（ k=0 时为第二）边值条件式（ 以左边界点 （1, j） 为例，见图 9-2, 利用一阶差商公式u u（0, j） u（1, j
5、）（0, j） O（h1）n h1则得到边界点 （0, j） 处的差分方程u0,j u1,jk0,j u0,j r0,jh1ui, j=（xi，yi）， （ i, j）h （9.7 ）9.4），9.8 ）联立差分方程（ 9.6 ）与（ 9.7 ）或（ 9.8 ）就形成了求解 Poisson 方程边值问题 的差分方程组， 它实质上是一个关于未知量 ui, j的线性代数方程组， 可采用第 2，3 章介绍的方法进行求解。这个方程组的解就称为偏微分方程的 差分近似解 ，简称 差 分解 。考虑更一般形式的二阶椭圆型方程u u u u （A ） （B ） C D Eu f （x,y）, （x,y） G （
6、9.9 ） x x y y x y0。引进半节点 x 1 xi 1h1,i 12 i 2 1,其中 A（ x, y） Amin0, B（x, y） Bmin 0, E（ x,y）u 1 u 1 u 1 2x（A xu）（i,j） h11 （A ux）（i 12, j） （A ux）（i 12, j） O（h12）u（i 1,j） u（i, j） A u（i, j） u（i 1, j） O（h2） h1 i 2,j h1O（h12）ux（i,j） u（i 1,j）2hu（i 1, j）x 2h1uu对 （B ）, 类似处理，就可推得求解方程（ 9.9 ）的差分方程 y y y9.10 ）ai 1
7、,jui 1,j ai 1,jui 1,j ai,j 1ui,j 1 ai, j 1ui,j 1 ai,jui, j f（i, j）, （i, j） Gh其中显然，当系数函数 A（x,y）=B（x,y）=1, C（ x, y）= D （ x, y）= E（ x, y）=0 时，椭圆型方程（9.9） 就成为 Poisson 方程（ 9.1 ），而差分方程（ 9.10 ）就成为差分方程（ 9.6 ）。容易看 出，差分方程（ 9.10 ）的截断误差为 O（h12 h22） 阶。9.1.2 一般区域的边界条件处理前面已假设 G 为矩形区域， 现在考虑 G 为一般区域情形， 这里主要涉及边界条 件的处理
8、。考虑 Poisson 方程第一边值问题9.12 ）u f（x,y）, （x,y） G u （x,y）, （x,y）其中 G 可为平面上一般区域，例 如为曲边区域。仍然用两组平行直线： x=x0+ih1, y=y0+jh2, i, j=0, 1, , 对区域 G 进行矩形网格剖分，见图 9-3 。如果一个内节点（ i , j）的四个相邻节点（ i+1,j），（i-1, j ），（ i , j +1）和（ i , j-1 ） 属于 G G ，则称其为 正则内点 ，见图 9-3 中打“。”号者；如果一个节点 （i, j）属于 G 且不为正则内点，则称其为 非正则内点 ，见图 9-3 中打“ . ”
9、号者。记正则 内点集合为 Gh ，非正则内点集合为 h 。显然，当 G 为矩形区域时，Gh Gh, h h 成立。在正则内点（ i, j ）处，完全同矩形区域情形，可建立五点差分格式u（B） xC xB u（A） xB xA u（C） O（h12） xC xA xC xA则得到点 B 处的方程线性插值法精度较高，为二阶近似。对每一个非正则内点进行上述处理，将所得到的方程与（ 9.13 ）式联立，就组 成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般区域上边值问题（ 9.12 ）的差分近似解。对于一般区域上二阶椭圆型方程（ 9.9 ）的第一边值问题，可完全类似处理。 第二
10、、三边值条件的处理较为复杂，这里不再讨论。9.2 抛物型方程的差分方法本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。9.2.1 一维问题作为模型，考虑一维热传导的初边值问题9.15 ）9.16 ）u（x,0） （x）, 0 x lu（0,t） g1（t）, u（l,t） g2（t）, 0 t T其中 a 是正常数，f（x,t）, （x）, g1（t）和 g2（t）都是已知的连续的函数。现在讨论求解问题（ 9.14 ）-（9.18） 的差分方法。首先对求解区域 G=0 x l, 0t T进行网格剖分。取空间步长 h=l/ N,时间步长 =T/M, 其中 N,M 是正整数， 作
11、两族平行直线x x j jh, j 0,1, ,Nt tk k , k 0,1, M将区域 G 剖分成矩形网格，见图 9-5 ，网格交点（ xj， tk）称为节点。用差分方法求解初边值问题（ 9.14 ）-（ 9.16 ）就是要求出精确解 u（x, t）在每 k个节点（ xj，t k）处的近似值 ujk u（xj,tk）。为简化记号， 简记节点（xj，tk）=u（j, k）。利用一元函数的 Taylor 展开公式，可推出下列差商表达式u（j,k） u（ j,k 1） u（ j,k） O（ ）tu（j,k） u（ j,k） u（j,k 1） O（ ） u（j,k） u（ j,k 1） u（ j,
12、k 1） O（ 2） t22u（j,k） u（j 1,k） 2u（j,k） u（ j 1,k） O（h2） x2 h21.9.17 ）9.18 ）9.19 ）（ 9.20 ）古典显格式在区域 G的内节点 （j,k）处，利用公式（9.17）和（9.20），可将偏微分方程 （9.14 ） 离散为u（j,k 1） u（j,k） au（j 1,k） 2u（j,k） u（j 1,k） f k O（ h2）h2 j其中 fjk f （xi,tk ）。舍去高阶小项 O（ h 2 ） ，就得到节点近似值 （差分解） ukj 所 满足的差分方程显然，在节点 （ j, k）处，差分方程（ 9.21 ）逼近偏微分方
13、程 （9.14 ）的误差为 O（ h2），这个误差称为 截断误差 ，它反映了差分方程逼近偏微分方程的精度。现将（ 9.21 ）式改写为便于计算的形式，并利用初边值条件（ 9.15 ）与（ 9.16 ）补充上初始值和 边界点方程，则得到k 1 k k k kuj ru j 1 （1 2r）uj ru j 1 f j9.22 ）j 1,2, ,N 1, k 0,1, ,M 1 u0j （xj ）, j 1,2, , N 1kku0k g1（tk ）, uNk g2（tk ）, k 0,1, M与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程k k 19.22 ），当第 k 层节
14、点值 ukj已知时，可直接计算出第 k+1层节点值 ukj 1。这样，从第 0 层已知值 u0j （xi ）开始，就可逐层求出各时间层的节点值。差分方程（9.22 ）的求解计算是显式的，无须求解方程组，故称为 古典显格式 。此外，在式（9.22 ）中，每个内节点处方程仅涉及 k和 k+1两层节点值，称这样的差分格式为双层格式 。差分方程（ 9.22 ）可表示为矩阵形式1 2r r r 1 2r rArr 1 2rk k k Tu （u1 , ,uN 1）Fk （f1k rg1（tk）, f2k, , fNk 2, fNk1 rg2（tk）T （ （x1）, , （xN 1）T2.古典隐格式在区
15、域 G的内节点 （ j, k）处，利用公式（ 9.18 ）和（9.20 ），可将偏微分方程 （9.14 ） 离散为u（j,k） u（j,k 1） u（j 1,k） 2u（ j,k） u（j 1,k） k 2a 2 f j O（ h ） h2舍去高阶小项 O（ h2 ），则得到如下差分方程它的截断误差为 O（ h2） ，逼近精度与古典显格式相同。改写（ 9.24 ）式为便于计 算的形式，并补充上初始值与边界点方程，则得到k k k k 1 k ruj 1 （1 2r）uj ru j 1 uj f j9.25 ）j 1,2, ,N 1, k 0,1, ,Mu0j （xj ）, j 1,2, , N
16、 1u0 g1（tk）, uN g2 （tk）, k 0,1, M与古典显格式不同，在差分方程（ 9.25 ）的求解中，当第 k-1 层值 ukj 1 已知时，必须通过求解一个线性方程组才能求出第 k层值ukj ，所以称（ 9.25 ）式为古典隐 格式 ，它也是双层格式。差分方程（ 9.25 ）的矩阵形式为9.26 ）Buk uk 1 F k, k 1,2 ,MBr r 1 2r向量 uk,Fk, 同（ 9.23 ）式中定义。从（ 9.26 ）式看到，古典隐格式在每一层计算 时，都需求解一个三对角形线性方程组，可采用追赶法求解。3.Crank-Nicolson 格式（六点对称格式）利用一元函数
17、 Taylor 展开公式可得到如下等式ut （j,k 12） u（j,k 1） u（j,k）其截断误差为 O（ 2 h2） ，在时间方向的逼近阶较显格式和隐格式高出一阶。这个 差分格式称为 Crank-Nicolson 格式 ，有时也称为 六点对称格式 ，它显然是双层隐式 格式。改写（ 9.27 ）式，并补充初始值和边界点方程得到方程这是一个三层显式差分格式。在逐层计算时，需用到 ukj 1和 ukj 两层值才能得到 k+1层值 ukj 1 。这样，从第 0层已知值 u0j （xj ）开始，还须补充上第一层值 u1j ，才能逐层计算下去。可采用前述的双层格式计算 u1j 。除上述四种差分格式外
18、，还可构造出许多逼近偏微分方程（ 9.14 ）的差分格式， 但并不是每个差分格式都是可用的。一个有实用价值的差分格式应具有如下性质：（ 1）收敛性 。对任意固定的节点（ xj，tk）, 当剖分步长 ,h 0时，差分解 ukj 应收敛到精确解 u（ xj，tk）。（ 2）稳定性 。当某一时间层计算产生误差时，在以后各层的计算中，这些误差 的传播积累是可控制而不是无限增长的。理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。因此，这里 主要研究差分格式的稳定性。作为例子，先考查 Richardson 格式的稳定性。设 ujk 是当计算过程中带有误差时，按 Richardson 格式（9.3
19、0 ）得到的实际计算值。 ukj 是理论值，误差 ejk ujk ukj 。假定右端项 f jk的计算是精确的，网比 r ，则 ejk满足j 2 jk 1 k1 k k kej ej （ej 1 2ej ej 1） （9.31 ）设前 k-1 层计算时精确的，误差只是在第 k 层 j0 点发生，即k 1 k k ej 0,ej0 ,ej 0, j j 0。则利用（ 9.31 ）式可得到误差 的传播情况，见表 9-1 。表 9-1 r=1/2 时 Richardson 格式的误差传播kjj0-4j0-3j0-2j0-1j0j0+1j0+2j0+3j0+4kk+1-2k+2-4 7k+3-6 17
20、-24 -6k+4-831-68 89k+5-10 49-144 273-388 k+671-260 641-1096 1311-260从表中看出，误差是逐层无限增长的。表中的计算虽然是就网比 r 12 进行的，实际上对任何 r0 都会产生类似现象，所以 Richardson 格式是不稳定的。 利用误差传播图表方法考查差分格式的稳定性虽然直观明了，但只能就具体取 定的 r 值进行，并且也不适用于隐式差分格式。9.2.2 差分格式的稳定性前节构造的几种双层差分格式都可以表示为如下的矩阵方程形式其中 H 称为 传播矩阵 。对于显格式 H=A, 隐格式 H=B-1, 六点对称格式 H=（I+B） （
21、I+ A） 。一般的三层格式也可以转化为双层格式。为了讨论方便，设在初始层产生误差 0 ，且假定右端项 Fk 的计算是精确的。用 uk表示当初始层存在误差 0时，由差分格式（ 9.32 ）得到的计算解，则k k 1 ku Hu F00u记误差向量 k uk uk，则 k 满足方程k H k 1, k 1,2,0为初始误差定义 9.1 称差分格式（ 9.32 ）是稳定的，如果对任意初始误差 0 在某种范数下满足k C 0 , k 0, 0 0其中 C 为与 h, 无关的常数。 这个定义表明，当差分格式稳定时，它的误差传播是可控制的。从（ 9.34 ）式递推得到k k 0 TH , 0 k 因此，
22、差分格式稳定的充分必要条件是定理 9.3 （稳定性必要条件） 差分格式稳定的必要条件是存在与 常数 M, 使谱半径（H ） 1 M定理 9.4 （稳定性充分条件） 设 H 为正规矩阵，即 HH * H * 式也是差分格式稳定的充分条件。下面讨论几种差分格式的稳定性。为便于讨论，引进 N-1 阶矩阵 01 1 0 1S1 0 1uk 满足方程（9.33 ）（ 9.34 ），误差向量 k（9.35 ）（ 9.36 ）h, 无关的（ 9.37 ）H , 则（ 9.37 ）这个特殊矩阵的特征值为9.38 ）Sj 2cos j , j 1,2, ,N 1 jN例 9-1 古典显格式 此时 H=A=（1-
23、2r）I+rS 。 利用（9.38 ）式和三角函数公式， 可求得 H 的特征值为2jj 1 4r sin2 （2jN ）, j 1,2, ,N 1j1 4r sin2（ ） 1, j 1,2, ,N 12Nr0 ，条件（ 9.37 ）成立。注意， H =B-1仍为实对称矩阵，所以古典隐例 9-2 古典隐格式 此时 H=B -1 征值为显然，对任意 格式对任何网比 r0 都是稳定的，称为 绝对稳定 。例 9-3 六点对称格式 此时 H=（I+B） -1（I+A） ，利用矩阵 A 和 B 的特征值可得 到矩阵 H 的特征值为2 4rsin2 （ j ）j 2N , j 1,2, ,N 1 j 2
24、j2 4rsin2 （ ） 2N则对任意 r0, 条件（ 9.37 ）成立。由于 A 和 B 均为实对称矩阵，且 AB=BA，则可 验证 H 也是实对称矩阵。所以六点对称格式是绝对稳定的。习题九9-1 试用五点差分格式求解 Poisson 方程的边值问题2 2 16, （x,y） G xyu| 0, （x,y）其中 G=-1 x, y1 。取步长 h=0.5 求解。9-2 试写出求解 Laplace 方程边值问题222 2 16, 0 x 4, 0 y 3x2 y2u（0,y） y（3 y）, u（4,y） 0, 0 y 3u（x,0） sin x, u（x,3） 0, 0 x 44的五点差分
25、格式，取步长 h=1, 并写出差分方程的矩阵形式。9-3 试用古典显格式求解热传导方程定解问题x20 x 1, 0 t Tu（x,0） 4x（1 x）, 0 x 1 u（0,t） u（1,t） 0, 0 t T只计算 k=1,2 两层上的差分解，取网比 r=1/6, h=0.29-4 用古典隐格式求解热传导方程定解问题u（x,0） sin x, 0 x 1 u（0,t） u（1,t） 0, 0 t 0.3取 0.1, h 0.2 精确解为 u（x,t） e t sin x。9-5 导出求解热传导方程的差分格式k 1 k k k 1uj u j u j u j 1 （ k 2 k k ）（1 ） 2 （uj 1 2uj u j 1）h2的截断误差，并选取 使其达到二阶。
