1、一、差商与微商,第二章：有限差分法初步,1 有限差分法基本概念,（i）、有限差分的数学基础是用差商代替微商。,有如下两种数学形式：,（i）微商（导数）的定义,是连续函数，则它的导数为：,若,（2.1）,式（2.1）右边,是有限的差商。,与,都不为零，,而式（2.1）左边,是,当,趋于零时极限情形下的差商，称之微商。,在,没有到达零之前，,只是,的近似。,趋于,的过程认为是近似向精确过渡，,用,代替,就是精确向近似过渡。,两者的差值,表示差商代替微商的偏差。,（ii）偏差-Taylor级数展开,（2.2）,稍加整理后可写成：,（2.3）,可见 与 只能是近似相等。,偏差为：,（iii）微商与差商
2、的几何意义,图2-1 差商与微商的比较,图2.1表示了差商与微商之间的关系。应当指出，用不同方法得到的差商去代替微商，它们带来的偏差是不同的。,向右（前）差商：,（iV）差商的几种表示,（2.4）,向左（后）差商：,（2.5）,中心差商，取向右差商与向左差商的平均值：,（2.6）,偏差分析：,将Taylor级数写成：,（2.7）,Taylor级数还可写成：,（2.8）,由式（2.7）可得,由式（2.8）可得,（2.9）,（2.10）,（2.9）+（2.10），得到,（2.11）,比较式（2.9）、（2.10）、（2.11）可看到，用不同的差商形式去代替微商，所带来的偏差是不同的。这些偏差都是截
3、去了Taylor级数展开式中的高阶项而引起的，常称“截断误差”。,用向右差商与向左差商代替微商，其截断误差为与,同量级的小量 ；,同量级的小量；,中心差商的截断误差小于向右差商或向左差商。,而用中心差商代替微商，其截断误差是与,讨论：,上述一阶差商一般仍是x的函数，对它们还可以求差商。这种一阶差商的差商称为二阶差商，它是二阶微商的近似，常用向右差商的向左差商来近似二阶微商，即：,（V）二阶差商,根据式（2.7）+（2.8），可得,（2.12）,由式（2.12）知，二阶差商的截断误差也为与同阶的小量。,结论：由于用差商代替微商必然带来截断误差，相应地用差分方程代替微分方程也必然带来截断误差。这是
4、有限差分法固有的。因此，在应用有限差分法进行数值解时，必须对差分的构成及其对方程造成的误差引起注意。,二、从微分形式出发的差分格式,图2.2给出了一个简单边界值问题。,图2.2 矩形区域离散化,问题是求图2.2所示的边值问题的解，其数学表达如下，方程：,（2.13）,边界条件：,（2.14）,（2.15）,（2.15）,（2.16）,式（2.13）、（2.14）、（2.15）、（2.16） 所示定解问题解法。,在问题的提法已经明白之后，差分格式的构成可通过以下几步来实现：（i） 区域离散法,下面分别予以说明:,（vi） 构成差分格式,（ii） 建立区域内差分方程,（iii） 边界条件的差分形式
5、,1区域离散化所谓离散化，就是把几何上连续的区域用一系列网格线把它划分开。一般说来，网格形式应视几何区域的不同而不同，对于矩形区域而言，用矩形的网格，如图2.2，用五条水平网线与五条垂直网线把矩形区域离散掉。网线与网线的交点称之为“节点”，节点与节点的距离称之为步长，x方向的步长表示为 ，y方向的步长表示为 。,节点编号：为便于计算，需对节点逐个编号。常用（i，j） 表示节点位置，其中，i、j是与网线相对应的正整数。,i，j 的排列：可有不同的方式。 习惯上，与x、y 轴相一致，i 由左而右逐个增长， j 由下而上逐个增长。,但也有，考虑到与矩阵的格式相一致， i 表示行数， 由上而下逐个增长
6、，j 表示列数，由左而右逐个增长。 这种从上到下，从左到右的编排与一般书写 习惯也是一致的，因此，在计算机上算题也常被采用。 在本章中，大都采用与坐标相一致的编排方法。,在区域内的节点称“内节点”，在边界上的节点称“边界节点”。图2.2所示边界是规则的,则节点或在区域内，或正好落在边界上。,步长 或 可以是不变的常量，即等步长，也可以在区域内的不同处是不同的，即变步长。如果区域内各处的温度梯度变化很大，则在温度变化剧烈处，网格布得密些；在温度变化不剧烈处，网格布得疏些。至于网格布置多少，步长取多大为宜，要根据具体问题，兼顾到计算的精确度与计算的工作量等因素而定。,步长:,从物理方面对区域离散化
7、可作这样的理解，即认为区域内离散的每个节点，都集中着它周围区域（尺度为步长）的热容，或者说，区域内连续分布的热容都被分别地集中到离散的节点上去了。这样，节点的温度代表着它周围区域的某种平均温度。一系列离散的节点温度值代表着连续区域内的温度分布。,区域离散化物理理解:,节点（i，j）处的温度表示成 。,2差分方程代替微分方程在上节我们已对有限差分法的数学基础作了简要 的介绍，说明了如何用差商代替微商，以及由此带来的误差。这里介绍用差商代替微商的办法来处理导热方程（2.13），得到相应的差分方程。,方程（2.13）,对区域内各个点都成立的，当然对任意一个内节点（i，j）也成立。,或者说，在（i，j
8、）处存在二阶偏微商 与，这些二阶偏微商所对应的差商可表示成：,i=2，3，4；j=2，3，4,（2.17）,（2.18）,i=2，3，4；j=2，3，4,其中， 与 表示相应的二阶差商 与二阶偏微商的差别为 与 的数量级。,将式（2.17）与（2.18）代入方程（2.13），得,（2.19）,式（2.19）中去掉 项，得到,（2.20）,i=2，3，4；j=2，3，4,式（2.20）被称为对应于方程（2.13）的差分方程。方程（2.20）被改写成：,（2.21）,若 ， ，则式（2.21）又被 改写成：,或,（2.22）,物理意义： 一点（i, j）处的温度是它周围4点温度的平均值。,由于差分
9、方程（2.20）是从式（2.19）中去掉项 得来的，称去掉的项 为差分方程（2.19）的截断误差。当 与 趋于零时，差分方程的截断误差也趋于零，即差分方程逼近微分方程。我们称这种逼近的差分方程与相应的微分方程为“相容”。,3边界条件的差分形式,对流换热边界条件：,（2.14）,-用T 对 x 的向前差商代替式（2.14）中的T 对 x 的一阶偏微商，使式（2.14）变成为如下差分形式 ：,这里介绍用差商代替微商的办法把定解问题中 的各种边界条件表示成差分的形式.,或,（2.23）,i=1；j=2，3，4,热流边界条件：,（2.15）,用T 对 y 的向前差商代替式（2.15）中T 对 y的一阶
10、偏微商，使式（2.15）变成为如下差分形式:,i=1，2，3，4，5；j=1,（2.24）,或,绝热边界条件：,变成为：,i=5；j=2，3，4,（2.25）,（2.16）,给定温度边界条件：,i=1，2，3，4，5；j=5,（2.26）,至此，我们对全部节点，包括内节点与边界节点，都用差分形式代替了原来的函数形式。对于内节点上差分形式，我们通称差分方程，因为内节点上温度都是未知的。对于边界节点的差分形式，在边界节点的温度为未知量时，它是差分方程。而对于边界节点为给定的温度时，得到的就不是差分方程了。但在实际应用中，人们往往习惯地把由内节点与边界节点建立起来的差分形式，都统称为差分方程。笼统地
11、讲，一个节点对应一个差分方程。,在边界节点的处理方面还有几点需要强调：,（i） 每一边界节点只应属于一种边界条件。如图2.2中，i=1，j=5的节点只属于边界条件式（2.16）。,（iii）若边界节点不正好落在区域的边界上，则需对它们进行特殊处理。,（ii）对应不同边界的差分方程式（2.23）、2.24）、（2.25）都是用一阶向前差商代替一阶微商得到的，也即它们的截断误差为,或,量级，与内节点差分方程的截断误差相比， 低了一个量级。这一点也是从微分形式出发 建立差分格式的弱点。,4差分格式的构成,由于式（2.13）、（2.14）、（2.15）、（2.16）、（2.17）所表示的方程式与边界条
12、件都是线性的，由此而得到的内节点与边界节点的差分方程也都是线性代数方程。由全部节点的差分方程构成一个线性代数方程组。在这个方程组里，方程式的个数等于节点的个数。针对前面讨论的例子，我们可以看到，这个有25个节点所构成的方程组，只需要5个式子，即式（2.21）、（2.23）、（2.24）、（2.25）、（2.26）就可表示，这里每个式子都表示几个节点方程。,也就是说，在组成方程组时，不必把每个节点方程都写出来，而只要写出几个规格化了的方程就可以了。因此，人们把规格化了的，由内节点与边界节点全部差分方程所构成的线性代数方程组，称之为“差分格式”。,一般地说，差分格式被写成如下的形式：,（2.27）
13、,其中，n是节点数，也即方程个数，每个方程对应一个节点。 和bi（i=1，2，n；j=1，2，n）都是常数 ，,方程组（2.27）可被进一步写成矩阵的形式 :,（2.28）,其中,如果在方程组中去掉其中已知温度节点的那些方程，由此构成的线性代数方程组中，方程的个数等于温度未知的节点个数，也即方程组的未知数。这样的方程组也是差分格式。也可写成式（2.28）。,综上所述，用有限差分法对式（2.13）、（2.14）、（2.15）、（2.16）、（2.17）所组成的边值问题的数值处理，最终归结成求解线性代数方程组（2.28），方程组的解即各节点的温度。如果整个区域的节点足够多，那么，离散节点的温度分布
14、就近似代替了区域内的连续温度分布。,为便于讨论各种差分格式的优缺点，最好把方程组（2.28）中系数矩阵具体地写出来。但当我们着手书写由式（2.21）、（2.23）、（2.24）、（2.25）、（2.26）组成的代数方程组时，发现它所占的版面太大，造成印刷的困难。,所以，为便于书写，采用图2.3所示的网格,并假定,得到由全部节点组成的线性代数方程组为：,（2.29）,表示成矩阵形式：,（2.30）,式（2.29）或（2.30）构成差分格式。若在方程组中去掉已知温度节点所对应的方程，即第1至第3个方程，则式（4.2.24）被改写成：,（2.31）,将式（2.31）与式（4.2.28）进行对照，即可
15、得到矩阵 、 、 的各个元素。这里特别提醒读者注意，在式（2.31）与（2.28）中温度T下角码的对应关系，在讨论二维稳定导热问题时，人们习惯用二个角码（i，j）来表示节点位置及对应的曾度 。,（二）解线性代数方程组的直接法,这里介绍计算机上常用的解线性代数方程组的直接法。 大家知道，线性方程组,（2.32）,中只要矩阵的行列式 ，方程组（2.32）就有唯一解，其表达式为:,（2.33）,其中 为用右端向量 替换行列式 的第j列而得的，这一公式为著名的Carmer法则。,显然，按Cramer法则求解方程组（2.32）需要计算n+1个n阶行列式。每个n阶行列式按直接展开办法来算，需作（n-1）n
16、！次乘法和n!次加法运算。当n=30时，共约需完成 次乘法和加法运算，这是一个十分惊人的数字，即使在一台每秒作一亿次运算的计算机上完成这一计算也是不可能的。所以，尽管这种办法也是一种直接法，并且理论上可行，但实际上是无法进行求解的。即使采用其它办法来计算行列式，按Cramer法则求解的工作量也比通常的直接法大得多。因而，Cramer法则对于数值计算来说是没有什么用处的，仅在一些特殊场合才有用。,另外，矩阵求逆也是人们熟悉的求解线性方程的一种方法。,(2.34),但求逆矩阵 时，要计算 个( n -1)阶行列式，和一个n阶行列式，它的计算工作量也是相当可观的，对于n较大的情况也没有什么现实意义。
17、,现在计算机上常用的直接解法大多数是以系数矩阵的三角形化为基础的。就是说说，先对方程组进行变换，使其化为等价的（即具有相同解的）三角形方程组。由于三角形方程组的求解十分容易，原方程的求解问题即告解决。下面我们简要地讨论这一类型的方法。为讨论方便起见，我们首先叙述三角形方程的解法，然后叙述将原方程化为等价三角形方程组的方法。由于计算机的字长是有限的，每次运算之后还要对结果进行舍入，所以，虽然理论上直接法在有限步内可以得到精确解，但计算机上实际得到的只是近似解。,1三角形方程组的解法所谓三解形方程组是指下面两种形式的方程组：,(2.35),或,(2.36),方程组(2.35)叫作下三角形方程组，方
18、程组 (2.36)叫作上三角形方程组。,若用矩阵符号可分别写为：,其中L为方程组(2.35)的系数所构成的下三角形矩阵，其元素满足关系：,U为方程组(2.36)的系数所构成的上三角形矩阵，其元素满足关系：,三角形方程组的求解是很简单的。方程组(2.35)的计算公式可归结为：,(2.37),这个计算过程通常也只作前推过程。,对于方程组(2.36)，其计算公式可归结为：,(2.38),这个计算过程通常也叫作回代方程。,由以上分析可以看到，只要把方程组化成了等价的三角形方程组，求解就容易了。,2Gauss消去法,Gauss消去法（简称消去法）的提出已有相当长的时间了，是一种古老的方法。然而，近年来在
19、计算机上求解线性代数方程组的实践表明，它仍是直接法中最常用的一种方法，也是最有效的方法之一。其基本思想是，用逐次消去一个未知数的办法把原来的方程组化为等价的（具有相同解）三角形方程组。这样，求解就很容易了。,假定把要求的n阶线性的方程组(2.32)改写成如下形式：,(2.39),用矩阵符号记为,其中 为 方阵， 为 向量，它们分别为：,分别从原方程组的第二个方程减去第一个方程乘以 ，第三个方程减去第一个方程乘以 ，如此等等，即可消去后面n-1个方程中的未知量 。,这时方程组(2.39)就变为如下等价方程组：,表示成矩阵形式为：,其中,若令 ，则系数 和 的计算公式应为：,类似地，分别从上述等价
20、方程组的第三个方程减去第二个方程乘以 ，第四个方程减去第二个方程乘以 ，如此等等，即可进一步消去后面n-2个方程中的未知量T2，而将方程(2.39)变为如下等价形式：,表示成矩阵形式为：,其中,若令 ，则系数 和 的计算公式应为：,上述的消去步骤还可以进行下去。如此继续之，重复上述步骤(n -1)次以后，我们即可得到如下等价三角形方程组：,(2.40),表示为矩阵形式：,其中 为如下上三角形矩阵， 为 向量；,(2.41),三角形方程组 很容易用前述的回代过程(2.38)求解，这样就完成了消去法求解n阶线性代数方程组的过程。从原来方程组(2.39)得出等价三角形方程组(2.40)的过程称之为消
21、去过程。采用前面的记号，我们可将消去过程的计算公式归结为对于 ，递推地计算如下各量：,(2.42),用Gauss消去法解线性代数方程组时，为了能求到最后的解（尽管已经具备了 的条件），并使解尽可能的精确，应注意如下两点：,（i）系数矩阵中对角线上的元素 都不应为零，因为在消元的过程中，不断用作为除数，倘若有一个 为零，计算就无法进行下去。（ii）在系数矩阵A每一行的元素中， 的绝对值最好比同一行的其他元素都大，这 样在作除法运算时，引起的舍入误差就比较小。,对照上节中差分格式，不难看到，由有限差分法得到的系数矩阵是能够满足以上两个条件的。因为每个节点方程（第i个方程）都代表着一个单元（第i个单
22、元）与周围单元或外界环境的热量交换关系。在这些热量交换中，无疑都与该单元的温度(Ti)有关，Ti的系数aii当然不能为零。,而且在第i个方程中，Ti的地位比其周围单元温度更为突出，表现在系数上，它的绝对值总是最大的。对于给定温度的节点方程而言，这种性质更明显。因为方程中除了该节点温度以外再也没有别的节点温度，它的系数当然也就最大了。,综上所述，由于对稳定导热问题用有限差分法得到的代数方程具有上述性质，因此在求解方程组时，可大胆放心使用Gauss消去法。（对于更一般的方程组，目前更多采用主元素消去法，而不用Gauss消去法。）,（三）解线性代数方程组的迭代法,前面介绍的解线性代数方程组的直接法对
23、于阶数不是很高的问题是非常有效的，这种场合一般不使用下面介绍的迭代法。然而对于阶数很高的稀疏矩阵，尽管提出了很多特殊的直接法来处理它们，在运算量和存储量的节省方面也取得了很大的进展，但仍然难于克服存储需要量大的缺点，特别在不具备大型计算机的条件下，采用下面介绍的迭代法更为合适。,迭代法的优点：由于不需要存储系数矩阵的零元素，所以占用的存储单元少。同时程序也比较简单，对于稳定导热用有限差分法所得到的方程组，求解收敛较快，因此广泛地被采用。迭代法的缺点是：它所得到的是一种近似解，在运算过程中需要进行多次迭代才能达到收敛指标的要求，而迭代次数事先是不知道的，这样，往往要耗费较多的时间。因此，一般地讲
24、，直接法与迭代法各有优缺点。,迭代法所要讨论的问题，仍是如下线性代数方程组：,(2.43),简写成：,(2.44),迭代法的基本思想是，构造一个由组成的向量序列，使其收敛于某个极限向量，并且 就是方程组(2.43)的精确解。根据构造向量序列的方法不同，常用的有简单迭代法，GaussSeidel迭代法与超松弛迭代法，下面分别予以介绍。,1简单迭代法最简单的迭代法称为简单迭代法，也称同步迭代法，Jakobi迭代法。迭代的最终目的是求解方程组(2.43)中的 。前面已经说到，用有限差分法（包括有限元法）得到的代数方程组能保证系数矩阵A对角线上元素不为零，即,则可将式(2.43)改写成：,(2.45)
25、,其中任一方程均可写成：,(2.46),任意给定各节点上的温度值 作为解的第零次近似，把它们代入式(2.46)的右端，由此算得的,作为解的第一次近似，把第一次近似得到的解再代入式(2.46)的右端，得到解的第二次近似。一般地讲，在已得到解的第k次近似 后，代入式(2.46)右端，得,为解的第k+1次的似。这样得到的序列为方程组的近似解。只要方程组(2.43)存在唯一解，则不论零次近似如何选取，当 时，此序列 必然收敛，且收敛于方程组的解 。实际计算中k不可能取 ，但可以说，当k充分大时，序列 已足够精确地接近方程组的。,通常，对充分大的k，其相邻两次迭代解之间的偏差小于预先给定的适当小量 （
26、大于零），即满足,就结束迭代过程，而取 作为方程组(2.43)的近似解。,2Gauss Seidel 迭代法简单迭代法虽然计算程序简单，但它计算每一个都要用到全部 的值，因此在计算机上，必须有两套工作单元来存放全部节点的旧值和新值。为了节省工作单元，并加快迭代收敛速度，对上述计算作如下修改。,方程组(2.45)中的第一式仍写成,或,方程组(2.45)中第二式，其中 换成第一次 算得 ，,(2.47),即,按此规律，可得到的一般关系式：,(2.48),这种计算方法称为Gauss-Seidel 迭代法，也称为异步迭代法。用Gauss-Seidel迭代法时，必须将节点或单元按顺序排列，并按顺序逐个进
27、行迭代。,由式(2.48)即可看到，用Gauss-Seidel迭代法进行迭代求解时，只需用一套工作单元存放近似值 或 ，这样节省了工作单元。同时在迭代过程中，有一半用迭代的新值，加速了收敛速度。因此，Gauss-Seidel迭代法是一种常用的方法。,3超松弛迭代法实际计算表明，当节点个数较多时，Gauss-Seidel迭代法的收敛速度仍然很慢，人们从改进Gauss-Seidel法的收敛速数出发，提出了超松弛法。要解的仍然是方程组(2.43)。先假定作为第零次近似的向量序列 。现在分两步来求 ，,第一步用Gauss-Seidel迭代法：,(2.49),求出一个近似解 ，第二步按下述公式(2.50)来改善 ，并得出第一次的近似解 。,或,(2.50),其中 是一个适当选择的参数，用它来改善 。,从式(2.49)与(2.50)中消去 ，可以得出由 直接求出 的公式。一般而言，由 算得 的公式为：,(2.51),这就是超松弛迭代法的计算公式， 为松弛因子。,由式(2.51)可以看到，当 时，式(2.51)简化成Gauss-Seidel代法的计算公式。由此可见，超松弛法是更为不般的迭代计算公式。对于不同的 值，超松驰法收敛的快慢也很不同，可以找出一个 值，使得超松驰法收敛得比取别的 值更快，此 称为最佳松驰因子。它的范围为 ，对于 的情况称为欠松驰，一般不被采用。对于 的情况称为超松驰。,
