CH02抛物方程差分法CH21-24,25-28

1、第第2 2章章抛物型方程的差分方法抛物型方程的差分方法 2.1 2.1 差分格式建立的基础差分格式建立的基础 2.2 2.2 显显式式差分格式差分格式 2.3 2.3 隐式差分格式隐式差分格式 2.4 2.4 解三对角形方程的追赶法解三对角形方程的追赶法 2.5 2.5 差分格式的稳定性和收敛性差分格式的稳定性和收敛性 2.6 2.6 非线性抛物型方程的差分解法举例非线性抛物型方程的差分解法举例 2.7 2.7 二维抛物型方程的差分格式二维抛物型方程的差分格式 2.8 2.8 交替方向的隐式差分格式交替方向的隐式差分格式( ( ADIADI 格式格式) ) 本章，我们研究线性抛物型方程的差分解

2、法，主要讨论差分方程的构造方法和有关的理论问题以及研究方法等，重点在于一维线性抛物型方程的差分方法，对于非线性以及多维抛物型方程的差分解法也进行了研究。 其中, 为 平面上某一区域。,),( , 0),(, 0),(),(txtxctxatxxtutxcxutxbxutxaxtutx),(),(),(),(2.1)众所周知，一维线性抛物型方程的一般形式为 (2) 初边值问题(或称混合问题) 通常考虑的定解问题有：(1) 初值问题(或称Cauchy问题) 在区域 上求函数，使满足Ttxtx0 ,| ),(xxxutx)()0 ,(),() 1 . 2(方程(2.2) 为给定的初始函数。)(xTt

3、ttuttuxxxutx0)(), 1 (),(), 0(10)()0 ,(),() 1 . 2(21方程(2.3)(2.4) 在区域上 ，使满足Ttxtx0 , 10| ),(),(txu边值条件初值条件求函数 为了构造微分方程(2.1)的有限差分逼近,首先将求解区域 用二组平行 轴和 轴的直线构成的网格覆盖，网格边长在方向 为 ，在 方向为 (如图2.1所示)。 分别称为空间方向和时间方向的步长，网格线的交点称为网格的结点。对初值问题来说，网格是xtxhx tkt kh,kTNNnnktn;, 2 , 1 , 0, 2, 1, 0mmhxm2.1 2.1 差分格式建立的基础差分格式建立的基

4、础在 上的结点称为边界结点,属于 内的结点 0t称为内部结点。对于初边值问题，设 ，则网格是Ttxtx0 , 10| ),(kTNNnnktn;, 2 , 1 , 01;, 2 , 1 , 0MhMmmhxm研究导数的差商近似表达式。对二元函数 定义 ，且假定 具有我们需要的有界偏导数。),(txuu ),(nmnmtxuu),(txuu 在 上的结点称为边界结点，属于 内的结点称为内部结点。1, 0, 0 xxt 差分方程就是在网格点上求出微分方程解的近似值的一种方法，因此又称为网格法。 构造逼近微分方程的差分方程的方法。构造逼近微分方程的差分方程的方法。由Taylor展开，有 nmnmnm

5、nmnmxuhxuhxuhtxutxu)(! 3)(! 2)(! 1),(),(3332221nmnmnmnmnmxuhxuhxuhtxutxu)(! 3)(! 2)(! 1),(),(3332221则 在 处对 的一阶偏导数有三个可能的近似:u),(nmtxxhuuhtxutxuxunmnmnmnmnm11),(),()(huuhtxutxuxunmnmnmnmnm11),(),()(huuhtxutxuxunmnmnmnmnm22),(),()(1111(2.5)(2.6)(2.7)向前差商向后差商中心差商 显然，用差商近似导数存在误差，令huuxuEnmnmnmnm1)(2.8)则1,2

6、2)(2mmtxnmxxxxuhEn 关于导数的近似差商表达式，也可以通过线性算子作为推导工具得到，定义： 截断误差,阶为)(hO)(hO用向后差商近似导数的截断误差阶也为)(2hO而中心差商近似导数的截断误差阶为xDx为 方向偏导数算子xxTnmnmxnmnmxuuTuuT111,为为 方向位移算子方向位移算子，xx)(212121nmnmnmxuuu为为 方向平均算子方向平均算子，x),2(21nmnmthxuu其中: 方向的差分算子方向的差分算子:xnmnmnmxuuu1(2.9)前差算子前差算子:,xxnmnmnmxuuu1(2. 10)后差算子后差算子:,中心差算子中心差算子: :(

7、2.11)nmnmnmxuuu2121x,nmnmxuuT2121,nmnmxuuT2121 建立差分算子和导数算子之间的关系，由建立差分算子和导数算子之间的关系，由Talyor 展开，有展开，有nmnmnmnmnmxuhxuhxuhuu)(! 3)(! 2)(! 13332221nmxxuDhDhI22! 2! 1为恒等算子IuhDnmx)exp(由nmxnmuTu1得)exp(xxhDT (2.12)或者xxThDln(2.13)同理有)exp(1xxhDT1lnxxThD因为ITITxxxx,故323121)ln(xxxxxIhD(2.14)同理323121)ln(xxxxxIhD(2.

8、15)因为 2121xxxTT)21exp()21exp(xxxhDhD)21sinh(2xxhD(2.16)则54232! 523! 321)21sinh(2xxxxxarhD(2.17) 式(2.14)，(2.15)，(2.17)分别给出了偏导数算子关于前差、后差、中心差的级数表达式双曲正弦3246nmxxxxxxnmxxxnmxxxnmuuuxuh533232)(403)(6131213121)(2.18.1)(2.18.2)(2.18.3) 利用这些关系式就可给出偏导数的差分表达式返回又由222)ln(xxIDh222)ln(xxIDh222)21sinh(2xxarDh可得二阶偏导数

9、的差分表达式nmxxxnmxxxnmxxxnmuuuxuh64233243222290112112111211)(2.19.1)(2.19.2)(2.19.3)返回返回4235nmxxxxxxnmxxxnmxxxnmuuuxuh753543543333)(12037)(21)(47234723)(2.20.1)(2.20.2)(2.20.3)nmxxxnmxxxnmxxxnmuuuxuh86465465444424076161726172)(2.21.1)(2.21.2)(2.21.3)对于三阶、四阶偏导数的差分表达式为 从以上这些偏导数的差分表达式，我们可以得到偏导数的各种精度的近似表达式。

10、nmnmnmxnmuuuxuh1)(且nmxnmxnmnmnmuuuuxuh3213121)()( 又由二阶导数的前差表达式(2.19.1)，得nmxnmuxuh2222)(因此)()(1)(1hOuuhxuEnmnmnmnm 在 的前差表达式中取第一项，则有nmxuh)(即截断误差阶 为。)(hO 现在研究构造微分方程(2.1)的差分方程的方法，为此记微分方程(2.1)为uDDtxLtuxx),(2(2.22) L 是关于 的线性算子， 。包括二个相邻时间层的网格结点的差分方程可以从Talor 展开式推出2,xxDDxDx),()! 3! 2! 11 (),(333222txutktktkk

11、txu),()exp(txutk返回设 ，于是) 1( ,(,1knmhuunktmhxnmnmnmutku)exp(1(2.23)如果算子L不依赖于t，即 ，则),(2xxDDxLL nmxxnmuarharhmkkLu)21sinh(2(),21sinh(2,(exp(21(2.25)21sinh(2xxarhD将式(2.17)， ,代入算子L中，即在L中用中心差分算子 代替了微分算子 ，于是有 xxD(2.24)nmnmukLu)exp(138 目前通常用于解方程(2.1)的各种差分方程，都是方程(2.25)的近似表达式。下面各节，我们将以式(2.25)为基础，对简单的抛物型方程，推导一

12、些常用差分格式。 对于用差分方法求偏导数方程的数值解来说，设计差分方程，用之作为微分方程的近似，仅仅是第一步。本章除致力于这一研究外，特别着重讨论了诸如差分格式的稳定性、收敛性等基本问题，它们也是本书研究的主要内容之一。2.2 2.2 显式差分格式显式差分格式 现在，对抛物型方程(2.1)的几种特殊情况，从方程(2.25)出发，构造微分方程的有限差分近似。2.2.1 2.2.1 一维常系数热传导方程的古典显示格式一维常系数热传导方程的古典显示格式 首先考虑一维热传导方程22xutu(2.26)的差分近似。差分方程的构造由 ，方程(2.24)为2xDL nmxnmukDu)exp(21nmxxu

13、DkkD2222)(211代入式(2.19.3)，得 )901121(164222xxxxhD算子之间的关系则nmxxxnmurrrrrru62421)15121(61)61(211(2.27)其中 为步长比。2hkr 返回在上式中，如果仅仅保留二阶中心差分，且设 为相应差分方程解在结点(mh,nk) 上的值，则nmUnmxnmUrU)1 (21(2.28)代入 的表达式，则得差分方程2xnmnmnmnmrUUrrUU111)21 (2.29)xxxuTtxxutu)()0 ,(0 ,022将格式(2.29)应用于解初值问题(初边值问题)古典显式差分格式图2.2差分格式(2.29)也可简单地由

14、导数的差商近似表达式得到)(1)(1nmnmnmuuktu)2(1)(11222nmnmnmnmuuuhxu代入微分方程(2.26)，并令差分方程解为 即可。虽然在边界结点上，差分方程和微分方程具有相同的初值或者初边值条件，但是，一般而言，结点 上微分方程的精确解 和古典显式差分格式(2.29)的精确解 不相等。nmU) 1,(nm1nmu1nmU111nmnmnmUuz(2.30)记 假定 具有下面推导中所需要的有界偏导数，则由 展开，有 ),(txuTaylornmnmnmnmnmtuktuktukuu)(6)(2)(3332221nmnmnmnmnmxuhxuhxuhuu)(6)(2)(

15、3332221nmnmxuhxuh)(120)(24555444nmnmnmnmnmxuhxuhxuhuu)(6)(2)(3332221nmnmxuhxuh)(120)(24555444截断误差截断误差42nmnmnmnmnmxutukuururu)()()21 (22111nmxurtuk)61(244222(2.31)则由式(2.26)，(2.29)，(2.30)，(2.31)得nmnmnmnmnmxurtukzzrzrz)61(2)()21 (44222111(2.32)从式(2.31)有nmnmnmnmnmxutukuururu)()()21 (22111nmxurtuk)61(214

16、422或nmnmnmnmnmnmxutuhuuukuu)(222211nmxurtuk)61(214422(2.33)从而，上式右边量描写了古典显式差分格式(2.29)在 点对微分方程的近似程度，将其定义为差分格式在点 的截断误差，记为 ，即),(nm),(nmnmRnmnmxurtukR)61(24422(2.34) 假定假定 在所考虑的区域保持有界，则古典显在所考虑的区域保持有界，则古典显式差分格式的截断误差阶为式差分格式的截断误差阶为 。4422,xutu)(2hkO442222,xutuxutu61r从式(2.33)又可见到，如令 ，因为故截断误差 的阶可以提高，这时 )(42hkOR

17、nmnmRnmxxnmUrrrU421)61(211(2.35.1)或者)(32(32)652(211121nmnmnmnmUUrrUrrU)(61 (12122nmnmUUrr(2.35.2)相应的截断误差阶为 。通常，格式可用图2.3表示。 )(42hkO 为了提高截断误差的阶，我们也可用在式(2.27)中保留四阶中心差分项的办法达到，这时有差分格式(2.27)m,n+1m-2,nm-1,nm,nm+1,nm+2,n图图2.32.3m,n+1m-1,nm,nm+1,n图图2.22.2返回2.2.2 2.2.2 系数依赖于系数依赖于 的一维热传导方程的显式格式的一维热传导方程的显式格式x0)

18、()(22xaxuxatu(2.36)这时， 。2)(xDxaL L保留右边前二项，由 ，则有差分方程2221xxhD)()(21 (111nmnmmnmmnmUUxraUxaU(2.37)nmxnmuDxkau)(exp(21nmxxxuaDaDkkaD)(2112222nmxxxxuaDDaDaakkaD )2(21143222则 这一差分格式可用图2.4表示，其中 ，这是一个显式差分格式，其截断误差阶为 。)(mhaa )(2hkOm,n+1m-1,nm,nm+1,n图图2.42.4 由方程右边22)()()(xuxaxuxaxuxaxuDxaDxaxx)()(2xxDxaDxaL)()

19、(2nmxxnmuDaaDku)(exp(21nmxxuDaaDk)(12 进一步，考虑热传导方程0)()(xaxuxaxtu(2.38)的差分近似。12 在上式中保留前二项，并且 和 分别用 和 代替，则得差分方程)(2111nmnmuuh)2(1112nmnmnmuuuhnmnmnmnmUaharUaharUraU111)21()21()21 (2.39)nmxuDnmxuD2 也可通过直接用中心差分算子 代替微分算子 的办法获得方程(2.38)的差分近似 xh1xDx)()(1)(121nmxmxnmnmUxahUUknmmmnmUxaxarU)()(121211hxxUxaUxarmm

20、nmmnmm21,)()(21121121(2.40)这也是一个显式差分格式。 格式(2.39)和(2.40)的截断误差阶都是 。易见，由)(2hkOaa,mhx 注：注： 均在 处计算。Deltaahxaxamm21)()(21ahxaxamm21)()(21,2xxDD代入格式(2.40)即为格式(2.39)，差分格式(2.40)的推导方法,即在微分方程中直接用差分算子代替 正如前面已经指出的是推导差分格式的一个常用方法。),(nmtxaa 显然，微分方程(2.36)，(2.38)中的 如果为 ，即其自变量包括空间变量和时间变量，这时差分格式(2.37)，(2.39)，(2.40)同样是微

21、分方程的具有截断误差阶 的差分近似，这时格式(2.37)，(2.39)中 和 ,格式(2.40)中 和 分别换成 , 。)(xa),( txa)(2hkO),(nmxtxaa)(21mxa)(21mxa),(21nmtxa),(21nmtxa2.3 2.3 隐式差分格式隐式差分格式 隐式差分格式特点： 1. 具有二个或二个以上结点处的值未知； 2. 计算工作量较大； 3. 稳定性较好。nmxnmukDu)exp(21得 nmnmxuukD12)exp(nmnmxxuuDkkD1422)211 (由22xutu推导其最简单的隐式差分逼近古典隐式格式。 现在对热传导方程2.3.1 古典隐式格式古典

22、隐式格式1715格式用图2.5表示，其截断误差阶为 ，与古典显式差分格式相同。 )(2hkO或者nmnmnmnmUrUUrrU11111)21 (2.41)nmnmxUUhk122)1 (保留二阶导数项，且以 替代 ，则得差分格式221xh2xD 我们也可通过直接用差分算子代替 的方法，即2,xxDDkuutunmnmnm11)(2111111222)(huuuxunmnmnmnm代入微分方程，得到格式(2.41)。古典隐式差分格式m,n+1m-1,n+1m+1,n+1m,n图图2.5 隐式差分格式是解热传导方程(2.26)的常用的差分格式，由式(2.24)，有 NicolsonCranknm

23、nmukLukL)21exp()21exp(1NicolsonCrank2.3.2 隐式格式隐式格式由2xDL 得1222)21(21211nmxxukDkDnmxxukDkD222)21(21211(2.42)42两边仅保留二项，用 代替 ，则得差分格式221xh2xDnmxnmxUrUr)211 ()211 (212(2.43)这是一个隐式差分格式，称为 差分格式，截断误差阶为 。NicolsonCrank)(22hkO)(21)1 (11111nmnmnmUUrUr)(21)1 (11nmnmnmUUrUr(2.44) 由于格式(2.44)中包括六个结点，故也称为六点格式(如图2.6所示

24、)。m,n+1m-1,n+1m+1,n+1m,n图图2.6m-1,nm+1,n44也可将kuutunmnmnm121)(21121111121222221)(huuuhuuuxunmnmnmnmnmnmnm代入微分方程(2.26)，得到 格式。NicolsonCrank nmxnmxukDukD)211 ()211 (212由式(2.19.3)，可令 nmxxnmxuhuD)1211 (12222则可得12222)1211 (1xxxhD 另一精度较高的六点差分格式，如前在式(2.42)中仅保留直到 的项，即有2xD13代入上式，则有如下差分格式：nmxnmxUrUr212)61(211)61

25、(211(2.45) 称为 差分格式。Douglas38)(42hkO截断误差阶 2323)()(121)()(2124214422122122khOxuhxuuuhnmnmnmnmx因为 )()()(1242144212khOxuhuuknmnmnmx48 前面，我们已经推导了热传导方程(2.26)的古典显式格式。古典隐式格式及 格式等。实际上，它们都可以作为本节推导的加权六点隐式格式的特殊情形。NicolsonCrank 2.3.3 加权六点隐式格式加权六点隐式格式由nmxnmukDu)exp(2110)1exp()exp(212nmxnmxukDukD得到14222211nmxxuDkk

26、DnmxxuDkkD4222)1 (21)1 (1即用 代替 ，则得差分格式221xh2xDnmxnmxUrUr212)1 (1)1 (或者)()(1 ()21 (1111111nmnmnmnmnmUUrUUrUr10)1 (21nmUr(2.46) 这是一个六点差分格式(如图2.7所示)，称为加权加权六点差分格式。400时, 为古典显式格式；1时, 为古典隐式格式；21时, 为 格式；NicolsonCranknmxnmxnmnmUhUhkUU2212211)1 (1加权六点格式亦可直接由差商代替导数得到 2.3.4 2.3.4 系数依赖于系数依赖于 的一维热传导方程的一个的一维热传导方程的

27、一个隐式格式的推导隐式格式的推导 由其 展开式可得Taylor)(12644222hODhDhxxx22),(xutxatu(2.47)的差分逼近。 考虑方程)21sinh(2xxhD已知x,t11令)()(111)1(22121222122khOuukahtuaxnmnmnmxnm代入式(2.48)，则)(1)(21121122nmnmnmnmnmxuukauuh)()(111212412122khOuuahrnmnmnmx)()(121)()(2124214422122122khOxuhxuuuhnmnmnmnmx)()1(12)1(242122221khOtuaxhtuanmnm(2.4

28、8)因此 43格式(2.49.1)具有截断误差阶 。)(24khO这是一个隐式差分格式(如图2.8所示)。1212121)611 (21)(1nmnmxnmnmnmUraUUranmnmxUra)611 (21212(2.49.1)因此得差分方程122112122121)(1211nmxnmnmxnmUraaanmxnmnmxnmUraaa22112122121)(1211(2.49.2)可写成形式m,n+1m-1,n+1m+1,n+1m,n图图2.8m-1,nm+1,n 前节引进的隐式差分方程，在要求解未知函数值的时间层 上包括三个未知函数值 。因此，这些隐式差分格式仅仅适合于解如图 中所示

29、的边值问题。在每一时间层，需要求解的隐式差分方程形成了一个线性代数方程组，它的系数矩阵是三对角形矩阵，即仅在主对角线及其相邻二条对角线上有非零元素。方程组写成一般形式是kntn) 1(111111,nmnmnmUUU)( 1 . 2b2.4 解三对角形方程的追赶法解三对角形方程的追赶法m,n+1m-1,n+1m+1,n+1m,nm,n+1m-1,n+1m+1,n+1m,nm-1,nm+1,n111213433323232221212111MMMMMdUUdUUUdUUUdUU(2.50)这一类方程可用追赶法求解。由方程组(2.50)中的第一个方程解出 ，得1U112111dUU将此式代入方程组

30、(2.50)的第二个方程，得到232221212)(dUUgU即2322gUU令 ，则上式可写为111111,dg1211gUU其中122122212222,gdg完全类似地，可以推出下面的公式22111MmgUUmmmm(2.51)其中21,111Mmgdgmmmmmmmmmmmm注意当 时, 。1m111111,dg即2112111MMMMMMMgdU1112121)(MMMMMMMdUgU 将关系式 代入式(2.50)中最后一个方程，得到2122MMMMgUU若令2112111MMMMMMMgdg则有 。11MMgU 如果 已经算出，那么解向量 的最后一个分量 就已求得，为了求得 的所有

31、分量，只有利用方程(2.51)即可逐步求出 ,因此，整个求解过程分为两大步： 1MgU1MUU1232,UUUUMM1222211,MMMgggg,第一步 依次确定12,11111Mmgdgdgmmmmmmm22,1111Mmmmmmm21,111MmUgUgUmmmmMM计算公式可归结为1221,UUUUMM第二步 依相反次序确定 通常，第1步称为“追”的过程，第2步称为“赶”的过程，整个求解过程称为追赶法。(2)0; 1, 2 , 1011MmmmMm定义(3)0; 1, 2 , 111MmmmMm定义则上述追赶法过程是稳定的。1, 3 , 20Mmm1, 2 , 10Mmm2, 2 ,

32、10Mmm(1)可以论证，如果例例 2.22.2 说明用 方法数值解如下定解问题的过程： NicolsonCrankTtttuttuxxuTtxxutut0)(), 1 (),(), 0(10)(|0 ; 1021022 由前已知 格式为NicolsonCrank)(21)1 ()(21)1 (1111111nmnmnmnmnmnmUUrUrUUrUr 如果选择 ，则 ，要解的方程组写成矩阵形式是81h8M171615141312111212112121121211212112121121211nnnnnnnUUUUUUUrrrrrrrrrrrrrrrrrrr18810076543212121

33、0000021211212112121121211212112121121211nnnnnnnnnnnrUrUrUrUUUUUUUUrrrrrrrrrrrrrrrrrrrkTNNn; 1, 2 , 1 , 0(2.52)相应于上述定解问题的差分方程组为011, 1 , 0UNneUBUAnnnnn其中， 为七阶方阵， 为列向量，它们的表达式从式(2.52)可知。因为在求第 层 时， 已计算得， (它们在 中出现)由边值条件已知，故方程组右边已知，且nnBA ,nnneUU,1) 1( n1nmUnnnUUU721,181080,nnnnUUUUne7 , 3 , 2021mrm7 , 2 ,

34、101mrm021rm又111mmmr7 , 11mrmmm,1g因此可用追赶法求解方程组(2.52)，由方程组右边值及 可求出 ，然后顺次，可求出 。mmm,717,g111617,nnnUUU计算结果。以及并以问题方法数值解下面的定解用例略解三对角方程组追赶法解三对角方程组追赶法20)()(),1 ()(0)(), 1 (),(), 0(10)(|0 ; 102 . 2. 2)(. 14 . 22121022TttxxxTtttuttuxxuTtxxutuNCt形式如下：维列向量。它们的具体为的矩阵，为其中，要解的方程组可写成对于确定的则比如格式为解：已知) 1(,) 1() 1(,)()255, 1 , 0(,255, 1 , 0;)7 , 6 , 5 , 4 , 3 , 2 , 1 (),(, 5 . 0/,128/1256/2, 8/1, 7 ,