抛物型方程的差分方法

抛物型方程的差分方法第一页，共八十三页，编辑于2023年，星期六第2章抛物型方程的差分方法

2.1差分格式建立的基础2.2显式差分格式2.3隐式差分格式2.4解三对角形方程的追赶法2.5差分格式的稳定性和收敛性2.6非线性抛物型方程的差分解法举例2.7二维抛物型方程的差分格式2.8交替方向的隐式差分格式(ADI格式)第二页，共八十三页，编辑于2023年，星期六

本章，我们研究线性抛物型方程的差分解法，主要讨论差分方程的构造方法和有关的理论问题以及研究方法等，重点在于一维线性抛物型方程的差分方法，对于非线性以及多维抛物型方程的差分解法也进行了研究。其中,为平面上某一区域。(2.1)众所周知，一维线性抛物型方程的一般形式为第三页，共八十三页，编辑于2023年，星期六(2)初边值问题(或称混合问题)通常考虑的定解问题有：(1)初值问题(或称Cauchy问题)在区域上求函数，使满足(2.2)为给定的初始函数。(2.3)(2.4)

在区域上求函数，使满足边值条件初值条件第四页，共八十三页，编辑于2023年，星期六为了构造微分方程(2.1)的有限差分逼近，首先将求解区域用二组平行于轴和轴的直线构成的网格覆盖，网格边长在方向为，在方向为(如图2.1所示)。分别称为空间方向和时间方向的步长，网格线的交点称为网格的结点。对初值问题来说，网格是2.1差分格式建立的基础在上的结点称为边界结点，属于内的结点

称为内部结点。第五页，共八十三页，编辑于2023年，星期六对于初边值问题，设，则网格是研究导数的差商近似表达式。为此对二元函数定义，且假定具有我们需要的有界偏导数。在上的结点称为边界结点，属于内的结点称为内部结点。差分方程就是在网格点上求出微分方程解的近似值的一种方法，因此又称为网格法。构造逼近微分方程的差分方程的方法。第六页，共八十三页，编辑于2023年，星期六由Taylor展开，有

则在处对的一阶偏导数有三个可能的近似:(2.5)(2.6)(2.7)向前差商向后差商中心差商第七页，共八十三页，编辑于2023年，星期六

显然，用差商近似导数存在误差，令(2.8)则关于导数的近似差商表达式，也可以通过线性算子作为推导工具得到，定义：截断误差,阶为用向后差商近似导数的截断误差阶也为而中心差商近似导数的截断误差阶为为方向偏导数算子第八页，共八十三页，编辑于2023年，星期六为方向位移算子，为方向平均算子，其中:

方向的差分算子:(2.9)前差算子:,(2.10)后差算子:,中心差算子:(2.11),,第九页，共八十三页，编辑于2023年，星期六

建立差分算子和导数算子之间的关系，由Talyor展开，有由得(2.12)或者(2.13)同理有第十页，共八十三页，编辑于2023年，星期六因为故(2.14)同理(2.15)因为

(2.16)则(2.17)

式(2.14)，(2.15)，(2.17)分别给出了偏导数算子关于前差、后差、中心差的级数表达式双曲正弦3246第十一页，共八十三页，编辑于2023年，星期六(2.18.1)(2.18.2)(2.18.3)

利用这些关系式就可给出偏导数的差分表达式返回第十二页，共八十三页，编辑于2023年，星期六又由可得二阶偏导数的差分表达式(2.19.1)(2.19.2)(2.19.3)返回返回4235第十三页，共八十三页，编辑于2023年，星期六(2.20.1)(2.20.2)(2.20.3)(2.21.1)(2.21.2)(2.21.3)对于三阶、四阶偏导数的差分表达式为第十四页，共八十三页，编辑于2023年，星期六

从以上这些偏导数的差分表达式，我们可以得到偏导数的各种精度的近似表达式。且

又由二阶导数的前差表达式(2.19.1)，得因此在的前差表达式中取第一项，则有即截断误差阶为。第十五页，共八十三页，编辑于2023年，星期六现在研究构造微分方程(2.1)的差分方程的方法，为此记微分方程(2.1)为(2.22)

L是关于的线性算子，。包括二个相邻时间层的网格结点的差分方程可以从Talor展开式推出返回第十六页，共八十三页，编辑于2023年，星期六设，于是(2.23)如果算子L不依赖于t，即，则(2.25)将式(2.17)，,代入算子L中，即在L中用中心差分算子代替了微分算子，于是有(2.24)返回3835第十七页，共八十三页，编辑于2023年，星期六目前通常用于解方程(2.1)的各种差分方程，都是方程(2.25)的近似表达式。下面各节，我们将以式(2.25)为基础，对简单的抛物型方程，推导一些常用差分格式。对于用差分方法求偏导数方程的数值解来说，设计差分方程，用之作为微分方程的近似，仅仅是第一步。本章除致力于这一研究外，特别着重讨论了诸如差分格式的稳定性、收敛性等基本问题，它们也是本书研究的主要内容之一。第十八页，共八十三页，编辑于2023年，星期六2.2显式差分格式

现在，对抛物型方程(2.1)的几种特殊情况，从方程(2.25)出发，构造微分方程的有限差分近似。2.2.1一维常系数热传导方程的古典显示格式

首先考虑一维热传导方程(2.26)的差分近似。差分方程的构造第十九页，共八十三页，编辑于2023年，星期六由，方程(2.24)为代入式(2.19.3)，得

算子之间的关系则(2.27)其中为步长比。返回第二十页，共八十三页，编辑于2023年，星期六在上式中，如果仅仅保留二阶中心差分，且设为相应差分方程解在结点(mh,nk)上的值，则(2.28)代入的表达式，则得差分方程(2.29)将格式(2.29)应用于解初值问题(初边值问题)古典显式差分格式图2.2第二十一页，共八十三页，编辑于2023年，星期六差分格式(2.29)也可简单地由导数的差商近似表达式得到代入微分方程(2.26)，并令差分方程解为即可。虽然在边界结点上，差分方程和微分方程具有相同的初值或者初边值条件，但是，一般而言，结点上微分方程的精确解和古典显式差分格式(2.29)的精确解不相等。(2.30)记第二十二页，共八十三页，编辑于2023年，星期六假定具有下面推导中所需要的有界偏导数，则由展开，有截断误差42第二十三页，共八十三页，编辑于2023年，星期六(2.31)则由式(2.26)，(2.29)，(2.30)，(2.31)得(2.32)从式(2.31)有第二十四页，共八十三页，编辑于2023年，星期六或(2.33)第二十五页，共八十三页，编辑于2023年，星期六从而，上式右边量描写了古典显式差分格式(2.29)在点对微分方程的近似程度，将其定义为差分格式在点的截断误差，记为，即(2.34)假定在所考虑的区域保持有界，则古典显式差分格式的截断误差阶为。从式(2.33)又可见到，如令，因为故截断误差的阶可以提高，这时。第二十六页，共八十三页，编辑于2023年，星期六(2.35.1)或者(2.35.2)相应的截断误差阶为。通常，格式可用图2.3表示。为了提高截断误差的阶，我们也可用在式(2.27)中保留四阶中心差分项的办法达到，这时有差分格式(2.27)第二十七页，共八十三页，编辑于2023年，星期六m,n+1m-2,nm-1,nm,nm+1,nm+2,n图2.3m,n+1m-1,nm,nm+1,n图2.2返回第二十八页，共八十三页，编辑于2023年，星期六2.2.2系数依赖于的一维热传导方程的显式格式(2.36)这时，。L保留右边前二项，由，则有差分方程(2.37)则第二十九页，共八十三页，编辑于2023年，星期六这一差分格式可用图2.4表示，其中，这是一个显式差分格式，其截断误差阶为。m,n+1m-1,nm,nm+1,n图2.4第三十页，共八十三页，编辑于2023年，星期六由方程右边进一步，考虑热传导方程(2.38)的差分近似。12在上式中保留前二项，并且和分别用和代替，则得差分方程(2.39)第三十一页，共八十三页，编辑于2023年，星期六也可通过直接用中心差分算子代替微分算子的办法获得方程(2.38)的差分近似(2.40)这也是一个显式差分格式。格式(2.39)和(2.40)的截断误差阶都是。易见，由注：均在处计算。Delta第三十二页，共八十三页，编辑于2023年，星期六显然，微分方程(2.36)，(2.38)中的如果为，即其自变量包括空间变量和时间变量，这时差分格式(2.37)，(2.39)，(2.40)同样是微分方程的具有截断误差阶的差分近似，这时格式(2.37)，(2.39)中和,格式(2.40)中和分别换成,。代入格式(2.40)即为格式(2.39)，差分格式(2.40)的推导方法,即在微分方程中直接用差分算子代替正如前面已经指出的是推导差分格式的一个常用方法。第三十三页，共八十三页，编辑于2023年，星期六2.3隐式差分格式

隐式差分格式特点：1.具有二个或二个以上结点处的值未知；2.计算工作量较大；3.稳定性较好。第三十四页，共八十三页，编辑于2023年，星期六得由推导其最简单的隐式差分逼近─古典隐式格式。现在对热传导方程2.3.1古典隐式格式1715第三十五页，共八十三页，编辑于2023年，星期六格式用图2.5表示，其截断误差阶为，与古典显式差分格式相同。或者(2.41)保留二阶导数项，且以替代，则得差分格式我们也可通过直接用差分算子代替的方法，即代入微分方程，得到格式(2.41)。古典隐式差分格式第三十六页，共八十三页，编辑于2023年，星期六m,n+1m-1,n+1m+1,n+1m,n图2.5第三十七页，共八十三页，编辑于2023年，星期六隐式差分格式是解热传导方程(2.26)的常用的差分格式，由式(2.24)，有2.3.2隐式格式由得(2.42)42第三十八页，共八十三页，编辑于2023年，星期六两边仅保留二项，用代替，则得差分格式(2.43)这是一个隐式差分格式，称为差分格式，截断误差阶为。(2.44)由于格式(2.44)中包括六个结点，故也称为六点格式(如图2.6所示)。第三十九页，共八十三页，编辑于2023年，星期六m,n+1m-1,n+1m+1,n+1m,n图2.6m-1,nm+1,n44第四十页，共八十三页，编辑于2023年，星期六也可将代入微分方程(2.26)，得到格式。第四十一页，共八十三页，编辑于2023年，星期六由式(2.19.3)，可令则可得另一精度较高的六点差分格式，如前在式(2.42)中仅保留直到的项，即有13代入上式，则有如下差分格式：(2.45)称为差分格式。38截断误差阶

23第四十二页，共八十三页，编辑于2023年，星期六23因为48第四十三页，共八十三页，编辑于2023年，星期六前面，我们已经推导了热传导方程(2.26)的古典显式格式。古典隐式格式及格式等。实际上，它们都可以作为本节推导的加权六点隐式格式的特殊情形。2.3.3加权六点隐式格式由得到即第四十四页，共八十三页，编辑于2023年，星期六用代替，则得差分格式或者(2.46)这是一个六点差分格式(如图2.7所示)，称为加权六点差分格式。40时,为古典显式格式；时,为古典隐式格式；时,为格式；第四十五页，共八十三页，编辑于2023年，星期六加权六点格式亦可直接由差商代替导数得到第四十六页，共八十三页，编辑于2023年，星期六2.3.4系数依赖于的一维热传导方程的一个隐式格式的推导

由其展开式可得(2.47)的差分逼近。考虑方程已知11第四十七页，共八十三页，编辑于2023年，星期六令代入式(2.48)，则(2.48)因此43第四十八页，共八十三页，编辑于2023年，星期六格式(2.49.1)具有截断误差阶。这是一个隐式差分格式(如图2.8所示)。(2.49.1)因此得差分方程(2.49.2)可写成形式第四十九页，共八十三页，编辑于2023年，星期六m,n+1m-1,n+1m+1,n+1m,n图2.8m-1,nm+1,n第五十页，共八十三页，编辑于2023年，星期六前节引进的隐式差分方程，在要求解未知函数值的时间层上包括三个未知函数值。因此，这些隐式差分格式仅仅适合于解如图中所示的边值问题。在每一时间层，需要求解的隐式差分方程形成了一个线性代数方程组，它的系数矩阵是三对角形矩阵，即仅在主对角线及其相邻二条对角线上有非零元素。方程组写成一般形式是2.4解三对角形方程的追赶法第五十一页，共八十三页，编辑于2023年，星期六m,n+1m-1,n+1m+1,n+1m,nm,n+1m-1,n+1m+1,n+1m,nm-1,nm+1,n第五十二页，共八十三页，编辑于2023年，星期六(2.50)第五十三页，共八十三页，编辑于2023年，星期六这一类方程可用追赶法求解。由方程组(2.50)中的第一个方程解出，得将此式代入方程组(2.50)的第二个方程，得到即令，则上式可写为其中第五十四页，共八十三页，编辑于2023年，星期六完全类似地，可以推出下面的公式(2.51)其中注意当时,。即将关系式代入式(2.50)中最后一个方程，得到若令第五十五页，共八十三页，编辑于2023年，星期六则有。如果已经算出，那么解向量的最后一个分量就已求得，为了求得的所有分量，只有利用方程(2.51)即可逐步求出,因此，整个求解过程分为两大步：,第一步依次确定计算公式可归结为第二步依相反次序确定第五十六页，共八十三页，编辑于2023年，星期六通常，第1步称为“追”的过程，第2步称为“赶”的过程，整个求解过程称为追赶法。(2)(3)则上述追赶法过程是稳定的。(1)可以论证，如果第五十七页，共八十三页，编辑于2023年，星期六例2.2说明用方法数值解如下定解问题的过程：由前已知格式为如果选择，则，要解的方程组写成矩阵形式是第五十八页，共八十三页，编辑于2023年，星期六第五十九页，共八十三页，编辑于2023年，星期六(2.52)第六十页，共八十三页，编辑于2023年，星期六相应于上述定解问题的差分方程组为其中，为七阶方阵，为列向量，它们的表达式从式(2.52)可知。因为在求第层时，已计算得，(它们在中出现)由边值条件已知，故方程组右边已知，且又第六十一页，共八十三页，编辑于2023年，星期六因此可用追赶法求解方程组(2.52)，由方程组右边值及可求出，然后顺次，可求出。第六十二页，共八十三页，编辑于2023年，星期六我们先看一个数值例子，考虑初边值问题(2.53)其中2.5差分格式的稳定性和收敛性2.5.1问题的提出第六十三页，共八十三页，编辑于2023年，星期六利用显式差分格式(2.29)，即式中。连同初值条件边值条件逐层解出结点处的值。现在对，取二种，使与。图2.9和图2.10中的曲线表示不同时刻微分方程的精确解，图中“·”表示差分方程的解。图2.9所示时的计算结果是曲线自上而下依次为微分方程第六十四页，共八十三页，编辑于2023年，星期六的精确解。黑点是用差分格式在时算出的相应各层上的近似值。二者符合得很好，由于对称性我们只给出一半图形。图2.10是当时差分方程解和微分方程精确解的图示，黑点仍表示差分方程解，其中分别为在时的计算结果。从图中看出，随着的增大，差分方程的解越来越远离微分方程的解。由此可见，值的不同，得出的结果有很大的差别，如的结果是可用的，但是时的结果就完全没有用。当然上面各种情况所得的差分方程解是由计算机得到的，不可能是差分方程理论上的准确解第六十五页，共八十三页，编辑于2023年，星期六，而是差分方程的近似解，我们用表示。显然与之间存在着差别，差分方程的准确解与微分方程的解之间，如前所述,也是有差别的。因而从计算机上解得的差分方程近似解与微分方程解之间的差别实质上包括两方面的差别，即(2.54)下面我们先研究上式右边第二项，即差分方程的理论解与计算机上解得的近似解之间的差别是随着的增大而无限增加还是有所控制。如果这种差别是无限增加，则称差分格式不稳定，显然不稳定的格式是不能使用的,因为误差的无限增加淹没了真解。上例中时就是差分方程不稳定的情况。从差分方程比如格式(2.29)可知，在求第六十六页，共八十三页，编辑于2023年，星期六第一层的差分方程解时，用到第0层上的值，也就是初始值。由于计算机存储数据为二进制数位的限制,不可能完全精确地存储在机器中，也就是计算用到的是带有误差的初始值。一般来说，在计算时又出现了误差，因此中包括了由于参加运算而出现的误差，即初始误差的传递，以及本身计算过程中出现的误差。这样，在第时间层计算时得到的是由于前面的误差传递和本身计算中出现的误差引起的。下面我们给出研究差分格式稳定性的最直接的方法，就是在第0层的一个结点上给出一个误差，然后研究这个误差的发展情况，即图方法。第六十七页，共八十三页，编辑于2023年，星期六假定在固定的某个结点引入一个误差，即把改成了，而在这一层的其他结点上的初值还是，假定用带有初始误差的初值按差分格式去计算以后各排结点上的值，且假定计算时没有引入其他误差，我们把得到的值记做，这样满足原来的差分格式。假如我们使用差分格式(2.29)，于是2.5.2图方法第六十八页，共八十三页，编辑于2023年，星期六显然两解之差满足(2.55)(2.56)以下分析当和时，随着增加而变化的情况。先看的情况，由式(2.55)得由此利用条件(2.56)即可算出的值(见表2.3)。第六十九页，共八十三页，编辑于2023年，星期六表2.3第七十页，共八十三页，编辑于2023年，星期六由表2.3可知，用显式差分格式(2.29)()计算时，由初始数据的误差，在以后各层所引起的误差是逐层减小的，这说明差分格式(2.29)当时是稳定的。再看的情形，由(2.55)得由此利用条件(2.56)即可得出的值(见表2.4)。表2.4第七十一页，共八十三页，编辑于2023年，星期六可见，用显式差分格式(2.29)()计算时，由初始数据的误差所引起的误差在以后各层的计算中逐层迅速增大，以致不能控制，因此差分格式(2.29)在时是不稳定的。用图方法讨论格式的稳定性能直观地看到差分格式是稳定性，缺点是必先固定。第七十二页，共八十三页，编辑于2023年，星期六2.5.3稳定性定义、稳定性分析的矩阵方法

以下讨论求初边值问题差分方程写成矩阵形式为：(2.57)第七十三页，共八十三页，编辑于2023年，星期六其中为维列向量；；为已知向量；为包括边值条件的向量；为阶方阵，可以随而改变。如果差分方程为显式，则对所有的(2.58)如果，则隐式格式可写成显式形式第七十四页，共八十三页，编辑于2023年，星期六设是初始值引进的误差向量，而在边值以及其他各层计算中未引入其它任何误差。由于的引入，差分方程的解为。则我们说差分格式是稳定的，其中是某一向量范数。稳定性的定义：时，对于任何的，差分格式得到的解满足不等式对于任意给定的，存在与无关且依赖于的正数，使当第七十五页，共八十三页，编辑于2023年，星期六设向量，则常用的向量范数有：(1)(2)(3)它们分别称为2-范数，1-范数和无穷范数，其中2-范数亦称为欧氏范数。第七十六页，共八十三页，编辑于2023年，星期六(1)，其中，为的共轭