第第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.1)众所周知，一维线性抛物型方程的一般形式为 (2) 初边值问题(或称混合问题) 通常考虑的定解问题有：(1) 初值问题(或称Cauchy问题) 在区域 上求函 数，使满足(2.2) 为给定的初始函数。(2.3)(2.4) 在区域上 求函数 ，使满足边值条件初值条件 为了构造微分方程(2.1)的有限差分逼近，首先将求解区域 用二组平行于 轴和 轴的直线构成的网格覆盖，网格边长在方向 为 ，在 方向为 (如图2.1所示)。 分别称为空间方向和时间方向的步长，网格线的交点称为网格的结点。对初值问题来说，网格是2.1 2.1 差分格式建立的基础差分格式建立的基础在 上的结点称为边界结点，属于 内的结点 称为内部结点。对于初边值问题，设 ，则网格是 研究导数的差商近似表达式。为此对二元函数 定义 ，且假定 具有我们需要的有界偏导数。在 上的结点称为边界结点，属于 内的结点称为内部结点。 差分方程就是在网格点上求出微分方程解的近似值的一种方法，因此又称为网格法。 构造逼近微分方程的差分方程的方法。构造逼近微分方程的差分方程的方法。由Taylor展开，有 则 在 处对 的一阶偏导数有三个可能的近似:(2.5)(2.6)(2.7)向前差商向后差商中心差商 显然，用差商近似导数存在误差，令(2.8)则 关于导数的近似差商表达式，也可以通过线性算子作为推导工具得到，定义： 截断误差,阶为用向后差商近似导数的截断误差阶也为而中心差商近似导数的截断误差阶为为 方向偏导数算子为为 方向位移算子方向位移算子，为为 方向平均算子方向平均算子，其中: 方向的差分算子方向的差分算子:(2.9)前差算子前差算子:,(2. 10)后差算子后差算子:,中心差算子中心差算子: :(2.11), 建立差分算子和导数算子之间的关系，由建立差分算子和导数算子之间的关系，由Talyor 展开，有展开，有由得(2.12)或者(2.13)同理有因为故(2.14)同理(2.15)因为 (2.16)则(2.17) 式(2.14)，(2.15)，(2.17)分别给出了偏导数算子关于前差、后差、中心差的级数表达式双曲正弦3246(2.18.1)(2.18.2)(2.18.3) 利用这些关系式就可给出偏导数的差分表达式返回又由可得二阶偏导数的差分表达式(2.19.1)(2.19.2)(2.19.3)返回返回4235(2.20.1)(2.20.2)(2.20.3)(2.21.1)(2.21.2)(2.21.3)对于三阶、四阶偏导数的差分表达式为 从以上这些偏导数的差分表达式，我们可以得到偏导数的各种精度的近似表达式。且 又由二阶导数的前差表达式(2.19.1)，得因此 在 的前差表达式中取第一项，则有即截断误差阶 为。 现在研究构造微分方程(2.1)的差分方程的方法，为此记微分方程(2.1)为(2.22) L 是关于 的线性算子， 。包括二个相邻时间层的网格结点的差分方程可以从Talor 展开式推出返回设 ，于是(2.23)如果算子L不依赖于t，即 ，则(2.25)将式(2.17)， ,代入算子L中，即在L中用中心差分算子 代替了微分算子 ，于是有 (2.24)返回3835 目前通常用于解方程(2.1)的各种差分方程，都是方程(2.25)的近似表达式。下面各节，我们将以式(2.25)为基础，对简单的抛物型方程，推导一些常用差分格式。 对于用差分方法求偏导数方程的数值解来说，设计差分方程，用之作为微分方程的近似，仅仅是第一步。本章除致力于这一研究外，特别着重讨论了诸如差分格式的稳定性、收敛性等基本问题，它们也是本书研究的主要内容之一。2.2 2.2 显式差分格式显式差分格式 现在，对抛物型方程(2.1)的几种特殊情况，从方程(2.25)出发，构造微分方程的有限差分近似。2.2.1 2.2.1 一维常系数热传导方程的古典显示格式一维常系数热传导方程的古典显示格式 首先考虑一维热传导方程(2.26)的差分近似。差分方程的构造由 ，方程(2.24)为代入式(2.19.3)，得 算子之间的关系则(2.27)其中 为步长比。返回在上式中，如果仅仅保留二阶中心差分，且设 为相应差分方程解在结点(mh,nk) 上的值，则(2.28)代入 的表达式，则得差分方程(2.29)将格式(2.29)应用于解初值问题(初边值问题)古典显式差分格式图2.2差分格式(2.29)也可简单地由导数的差商近似表达式得到代入微分方程(2.26)，并令差分方程解为 即可。虽然在边界结点上，差分方程和微分方程具有相同的初值或者初边值条件，但是，一般而言，结点 上微分方程的精确解 和古典显式差分格式(2.29)的精确解 不相等。(2.30)记 假定 具有下面推导中所需要的有界偏导数，则由 展开，有 截断误差截断误差42(2.31)则由式(2.26)，(2.29)，(2.30)，(2.31)得(2.32)从式(2.31)有或(2.33)从而，上式右边量描写了古典显式差分格式(2.29)在 点对微分方程的近似程度，将其定义为差分格式在点 的截断误差，记为 ，即(2.34) 假定假定 在所考虑的区域保持在所考虑的区域保持有界，则古典显式差分格式的截断误差阶为有界，则古典显式差分格式的截断误差阶为 。从式(2.33)又可见到，如令 ，因为故截断误差 的阶可以提高，这时 。 (2.35.1)或者(2.35.2)相应的截断误差阶为 。通常，格式可用图2.3表示。 为了提高截断误差的阶，我们也可用在式(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 系数依赖于系数依赖于 的一维热传导方程的显式的一维热传导方程的显式格式格式(2.36)这时， 。L保留右边前二项，由 ，则有差分方程(2.37)则 这一差分格式可用图2.4表示，其中 ，这是一个显式差分格式，其截断误差阶为 。m,n+1m-1,nm,nm+1,n图图2.42.4 由方程右边 进一步，考虑热传导方程(2.38)的差分近似。12 在上式中保留前二项，并且 和 分别用 和 代替，则得差分方程(2.39) 也可通过直接用中心差分算子 代替微分算子 的办法获得方程(2.38)的差分近似 (2.40)这也是一个显式差分格式。 格式(2.39)和(2.40)的截断误差阶都是 。易见，由注：注： 均在 处计算。Delta 显然，微分方程(2.36)，(2.38)中的 如果为 ，即其自变量包括空间变量和时间变量，这时差分格式(2.37)，(2.39)，(2.40)同样是微分方程的具有截断误差阶 的差分近似，这时格式(2.37)，(2.39)中 和 ,格式(2.40)中 和 分别换成 , 。代入格式(2.40)即为格式(2.39)，差分格式(2.40)的推导方法,即在微分方程中直接用差分算子代替 正如前面已经指出的是推导差分格式的一个常用方法。2.3 2.3 隐式差分格式隐式差分格式 隐式差分格式特点： 1. 具有二个或二个以上结点处的值未知； 2. 计算工作量较大； 3. 稳定性较好。得 由推导其最简单的隐式差分逼近古典隐式格式。 现在对热传导方程2.3.1 古典隐式格式古典隐式格式1715格式用图2.5表示，其截断误差阶为 ，与古典显式差分格式相同。 或者(2.41)保留二阶导数项，且以 替代 ，则得差分格式 我们也可通过直接用差分算子代替 的方法，即代入微分方程，得到格式(2.41)。古典隐式差分格式m,n+1m-1,n+1m+1,n+1m,n图图2.5 隐式差分格式是解热传导方程(2.26)的常用的差分格式，由式(2.24)，有 2.3.2 隐式格式隐式格式由得(2.42)42两边仅保留二项，用 代替 ，则得差分格式(2.43)这是一个隐式差分格式，称为 差分格式，截断误差阶为 。(2.44) 由于格式(2.44)中包括六个结点，故也称为六点格式(如图2.6所示)。m,n+1m-1,n+1m+1,n+1m,n图图2.6m-1,nm+1,n44也可将代入微分方程(2.26)，得到 格式。由式(2.19.3)，可令 则可得 另一精度较高的六点差分格式，如前在式(2.42)中仅保留直到 的项，即有13代入上式，则有如下差分格式：(2.45) 称为 差分格式。38截断误差阶 2323因为 48 前面，我们已经推导了热传导方程(2.26)的古典显式格式。古典隐式格式及 格式等。实际上，它们都可以作为本节推导的加权六点隐式格式的特殊情形。2.3.3 加权六点隐式格式加权六点隐式格式由得到即用 代替 ，则得差分格式或者(2.46) 这是一个六点差分格式(如图2.7所示)，称为加权加权六点差分格式。40时, 为古典显式格式；时, 为古典隐式格式；时, 为 格式；加权六点格式亦可直接由差商代替导数得到 2.3.4 2.3.4 系数依赖于系数依赖于 的一维热传导方程的一个的一维热传导方程的一个隐式格式的推导隐式格式的推导 由其 展开式可得(2.47)的差分逼近。 考虑方程已知11令代入式(2.48)，则(2.48)因此 43格式(2.49.1)具有截断误差阶 。这是一个隐式差分格式(如图2.8所示)。(2.49.1)因此得差分方程(2.49.2)可写成形式m,n+1m-1,n+1m+1,n+1m,n图图2.8m-1,nm+1,n 前节引进的隐式差分方程，在要求解未知函数值的时间层 上包括三个未知函数值 。因此，这些隐式差分格式仅仅适合于解如图 中所示的边值问题。在每一时间层，需要求解的隐式差分方程形成了一个线性代数方程组，它的系数矩阵是三对角形矩阵，即仅在主对角线及其相邻二条对角线上有非零元素。方程组写成一般形式是2.4 解三对角形方程的追赶法解三对角形方程的追赶法m,n+1m-1,n+1m+1,n+1m,nm,n+1m-1,n+1m+1,n+1m,nm-1,nm+1,n(2.50)这一类方程可用追赶法求解。由方程组(2.50)中的第一个方程解出 ，得将此式代入方程组(2.50)的第二个方程，得到即令 ，则上式可写为其中完全类似地，可以推出下面的公式(2.51)其中注意当 时, 。即 将关系式 代入式(2.50)中最后一个方程，得到若令则有 。 如果 已经算出，那么解向量 的最后一个分量 就已求得，为了求得 的所有分量，只有利用方程(2.51)即可逐步求出 ,因此，整个求解过程分为两大步： ,第一步 依次确定计算公式可归结为第二步 依相反次序确定 通常，第1步称为“追”的过程，第2步称为“赶”的过程，整个求解过程称为追赶法。(2)(3)则上述追赶法过程是稳定的。(1)可以论证，如果例例 2.2 2.2 说明用 方法数值解如下定解问题的过程： 由前已知 格式为 如果选择 ，则 ，要解的方程组写成矩阵形式是(2.52)相应于上述定解问题的差分方程组为其中， 为七阶方阵， 为列向量，它们的表达式从式(2.52)可知。因为在求第 层 时， 已计算得， (它们在 中出现)由边值条件已知，故方程组右边已知，且又因此可用追赶法求解方程组(2.52)，由方程组右边值及 可求出 ，然后顺次，可求出 。我们先看一个数值例子，考虑初边值问题(2.53)其中2.5 差分格式的稳定性和收敛性差分格式的稳定性和收敛性2.5.1 问题的提出问题的提出利用显式差分格式(2.29)，即式中 。连同初值条件边值条件逐层解出结点处的 值。 现在对 ，取二种 ，使 与 。图2.9和图2.10中的曲线表示不同时刻微分方程的精确解，图中“”表示差分方程的解。 图2.9所示 时的计算结果是曲线自上而下依次为 微分方程的精确解。黑点是用差分格式在 时算出的相应各层上的近似值。二者符合得很好，由于对称性我们只给出一半图形。 图2.10是当 时差分方程解和微分方程精确解的图示，黑点仍表示差分方程解，其中 分别为在 时的计算结果。从图中看出，随着 的增大，差分方程的解越来越远离微分方程的解。 由此可见， 值的不同，得出的结果有很大的差别，如 的结果是可用的，但是 时的结果就完全没有用。 当然上面各种情况所得的差分方程解是由计算机得到的，不可能是差分方程理论上的准确解 ，而是差分方程的近似解，我们用 表示。显然 与 之间存在着差别，差分方程的准确解 与微分方程的解 之间，如前所述,也是有差别的。因而从计算机上解得的差分方程近似解 与微分方程解 之间的差别实质上包括两方面的差别，即(2.54) 下面我们先研究上式右边第二项，即差分方程的理论解与计算机上解得的近似解之间的差别是随着的增大而无限增加还是有所控制。如果这种差别是无限增加，则称差分格式不稳定，显然不稳定的格式是不能使用的,因为误差的无限增加淹没了真解。上例中 时就是差分方程不稳定的情况。从差分方程比如格式(2.29)可知，在求 第一层的差分方程解 时，用到第0层上的 值，也就是初始值。由于计算机存储数据为二进制数位的限制, 不可能完全精确地存储在机器中，也就是计算 用到的是带有误差的初始值 。一般来说，在计算 时又出现了误差，因此 中包括了由于 参加运算而出现的误差，即初始误差的传递，以及本身计算过程中出现的误差。这样，在第 时间层计算 时得到的 是由于前面的误差传递和本身计算中出现的误差引起的。下面我们给出研究差分格式稳定性的最直接的方法，就是在第0层的一个结点上给出一个误差 ，然后研究这个误差的发展情况，即 图方法。 假定在固定的某个结点 引入一个误差 ，即把 改成了 ，而在这一层的其他结点上的初值还是 ，假定用带有初始误差的初值 按差分格式去计算以后各排结点上的 值，且假定计算时没有引入其他误差，我们把得到的值记做 ，这样 满足原来的差分格式。假如我们使用差分格式(2.29)，于是2.5.2 2.5.2 图方法图方法显然两解之差 满足(2.55)(2.56) 以下分析当 和 时， 随着 增加而变化的情况。先看 的情况，由式(2.55)得由此利用条件(2.56)即可算出 的值(见表2.3)。 表2.3 由表2.3可知，用显式差分格式(2.29)( )计算时，由初始数据的误差，在以后各层所引起的误差是逐层减小的，这说明差分格式(2.29)当 时是稳定的。再看 的情形，由(2.55)得由此利用条件(2.56)即可得出的值(见表2.4)。表2.4 可见，用显式差分格式(2.29)( )计算时，由初始数据的误差所引起的误差在以后各层的计算中逐层迅速增大，以致不能控制，因此差分格式(2.29)在 时是不稳定的。 用 图方法讨论格式的稳定性能直观地看到差分格式是稳定性，缺点是必先固定 。 2.5.3 2.5.3 稳定性定义、稳定性分析的矩阵方法稳定性定义、稳定性分析的矩阵方法 以下讨论求初边值问题差分方程写成矩阵形式为：(2.57)其中为 维列向量； ； 为已知向量； 为包括边值条件的向量； 为 阶方阵，可以随而改变。如果差分方程为显式，则对所有的(2.58)如果，则隐式格式可写成显式形式 设 是初始值引进的误差向量，而在边值以及其他各层计算中未引入其它任何误差。由于的引入，差分方程的解为 。则我们说差分格式是稳定的，其中是某一向量范数。 稳定性的定义：时，对于任何的 ，差分格式得到的解 满足不等式 对于任意给定的 ，存在与 无关且依赖于 的正数 ，使当 设向量 ，则常用的向量范数有： (1) (2) (3) 它们分别称为2-范数，1- 范数和无穷范数，其中2-范数亦称为欧氏范数。 (1) ，其中 ，为 的共轭 转置矩阵， 为 的最大特征值； (2)(3) 设矩阵 为A的元素，则相应的矩阵范数为：分别称为矩阵 A 的2-范数，1-范数和无穷范数。对所有的范数都有 其中 为矩阵 的谱半径， 为矩阵A的特征值。 故满足 上面定义的稳定性，只考虑初始值引进的误差的传播，称为差分格式关于初始值的稳定性。因为 满足如下方程因此，可推得如果存在一正常数 ，使在一定范数下满足(2.59)则差分格式(2.57)，(2.58)稳定。 通过对矩阵C的直接估计探求差分格式稳定性条件的方法称为稳定性分析的直接法(矩阵法)。解抛物型方程初边值问题的差分格式常利用矩阵法求得稳定性条件。 设 为 的特征值，用 表示 的最大值，即 的谱半径，则有定理定理 2.1 2.1 差分格式(2.60)稳定的必要条件是，存在与 无关的常数 ，使矩阵 的谱半径满足(2.62) 以下仅讨论差分方程系数不依赖于时间层数，即 ，故 ，差分方程为(2.60)存在常数 ，使(2.61)稳定性条件为：充要条件定理定理 2.3 若在差分格式(2.60)中， 为正规矩阵，即其满足 ，则条件(2.62)是差分格式(2.60)按欧几里德范数稳定的充分条件。 稳定性的必要条件(2.62)十分重要，在很多情况下，它也式充分条件。应用矩阵的欧几里德范数，则我们有以下定理。定理定理 2.2 若A为正规矩阵，则。有关矩阵特征值计算的几个结论： (1) 阶三对角线方阵 的特征值为这里a, b, c 可以是实数或复数。 (3) 设 为 的多项式， 非奇异，则 为矩阵 和 的特征值，相应的特征向量为 。 (2) 如果 是A的具有系数 的多项式，于是表明 为 的特征值，相应的特征向量为x。