工程数值方法上(2015-3-17讲稿)

1、工程数值方法学习内容：Chapter 1 线性代数方程组的数值解法Chapter 2 插值问题与数值微分Chapter 3 数值积分方法Chapter 4 常微分方程(组)初值问题的数值方法Chapter 5 常微分方程(组)边值问题的数值方法Chapter 6 椭圆型偏微分方程的数值方法Chapter 7 加权残值方法参考书目：1 武汉大学、山东大学合编，计算方法，高教版，19792 林成森编，数值计算方法（上、下），科学出版社，20003 中科院研究生数学丛书，工程中的数值方法，科学出版社，20004 曾绍林编，工程数学基础（研究生数学丛书），科学出版社，20015 李庆扬编，数值分析基础

2、教程，高等教育出版社，20026 李庆扬编，数值分析（第4版），清华版，20037 关治编，数值计算方法，清华版，20048 李岳生、黄有谦编，数值逼近，高教版，19789 李荣华编，微分方程数值解法，人教版，198010 邱吉宝编著，加权残值法的理论与应用，宇航版，1992Chapter 1 线性代数方程组的数值解法 线性代数方程组的求解是工程实践中最常遇到的问题。据不完全统计，在工程实践中提出的计算问题中，有近一半涉及到求解线性方程组。例如：结构有限元分析问题，大地测量问题，气象预报问题，电力传输网分析问题，各种电路分析问题，数据拟合问题，以及非线性方程组与微分方程的数值求解问题等等。因此

3、，学习并掌握线性代数方程组求解的基本理论与方法无疑是十分必需的。 本章将介绍目前一些利用计算机求解线性代数方程组常用的、且简单有效的数值方法。 求解线性方程组的数值方法尽管很多，但归并起来可分为两大类：（1） 直接法（精确法）凡经有限次的四则运算，若运算中没有舍入误差即可求得方程组精确解的方法。如：克莱姆（Cramer）法则方法、消元法、分解法、分解法等等。（2） 迭代法（近似法）将求解方程组的问题转化为构造一个无限迭代的序列，在实现该序列过程中的每一步计算结果，均是把前一步所得的结果施行相同的计算步骤进行修正而获得的，而这一无限序列的极限就是原方程组的精确解答。如：简单迭代法、赛德尔迭代法、

4、牛顿法、共轭斜量法等等。需要指出的，在一般情况下，我们使用直接法和迭代法两类方法都不可能完全获得原方程组的精确解答。原因很显然：（1）实际中在使用直接法时不可能没有数值计算的舍入误差，故此时所谓精确方法的解并不是绝对精确的；（2）实际中在使用迭代法时，不可能将极限过程无限进行到底，而只能进行有限次的迭代，故获得是满足精度要求的近似解答。 关于这两类方法求解的误差分析，我们将在每类方法的介绍之后进行简要讨论。§1.1 直接法Cramer法则与求逆方法设n元n个非齐次线性代数方程组为： (1.1)利用矩阵和向量符号，这个方程组可表为: (1.2)其中,方程组（1.1）的系数矩阵（阶方阵）

5、； 方程组（1.1）的右端已知向量（阶列阵）； 方程组（1.1）待求的解向量（阶列阵）。§1.1.1 Cramer法则由线性代数中关于线性方程组解的定理可知：若A的行列式，则方程组（1.1）有唯一解。此时，根据Cramer法则，其解的表达式为： (1.3)其中，即将中的第j列元素依次换为右端已知向量的元素所构成的阶方阵的行列式。易见，利用Cramer法则给出的（1.3）式求解n阶线性方程组（1.2）式，需要计算个阶行列式，而每个n阶行列式将有项，其中每一项又含有n个因子，故展开每一个n阶行列式仅乘法运算就需次（忽略加法运算次数），而对个n阶行列式，其乘法运算的次数为：对式（1.3）其

6、除法次数为：n故求解方程组（1.2）其乘法和除法总的运算次数为：次例如，若阶，则次这是一个十分惊人的数字，即使利用超高速的电子计算机能够胜任此计算次数，仅由于多个数的连乘亦有可能造成溢出而无法继续运算。因此，Cramer法则这个在理论上完善且精确的求解方法，仅在理论上和一些特殊情况下可以发挥作用，而对高阶线性代数方程组的实际求解中几乎没有多少实用价值。为此，人们不得不研究其它一些计算简单、且行之有效的求解方法。§1.1.2 求逆方法若，则A非奇异，即存在，则方程组（1.2）的解向量可表为： （1.4）常用的矩阵求逆方法有：（1）伴随矩阵法；（2）初等变换法；但当方阵A的阶数n较大时，

7、求逆非常麻烦且计算量很大。故对高阶线性代数方程组来说，求逆方法也是一种中看不中用的方法。§1.2 直接法Gauss（高斯）消去法尽管这是一种较为古老的方法，但至今仍不失为最常用和最有效的方法之一。基本思想：通过逐次消元处理，将原方程组化为等价的三角形方程组进行求解。§1.2.1 三角形方程组所谓三角形方程组无非是以下两种形式的方程组：(1)下三角形式(2.1)矩阵记为 (2.1')其中系数矩阵的元素满足关系：主对角线以上的元素均为零。（2）上三角形式 （2.2）矩阵记为 （2.2'）其中系数矩阵中主对角线以下的元素均为零，即：对三角形方程组的求解是十分简单的

8、。显然对于下三角方程组（2.1），其求解步骤如下： 从第一个方程中解得：； 将代入第二个方程中，从中解得：； 将代入第三个方程中，从中解得：；······如此逐个方程求解的过程向前递推下去，直到第n步。 将代入第n个方程中，从中解得：。对上述过程其完整的求解计算格式可归结为： （2.3）这一求解过程称为前推过程。同理，对于上三角方程组（2.2），其求解步骤亦可如法泡制，即： 从第n个方程中解得：； 将代入第n-1个方程中，从中解得：； 将、代入第n-2个方程中，从中解得：；·····&

9、#183; 将代入第1个方程中，从中解得。对上述过程其完整的求解计算格式可归纳为： （2.4）这一求解过程称为回代过程。通常是将（2.4）式合并为一式表为： （2.4a）式中约定：当足标j的取值大于其上界n时，和式。§1.2.2 Gauss消去法高斯消去法的求解过程就是首先利用矩阵的初等行变换方法将原方程组逐次消元，使之化为等价的具有上三角形式的方程组，然后再按上三角方程求解的计算格式（2.4a）式求出原方程组的解。整个求解过程可分为消元和回代两个过程。1. 简例 为了便于说明高斯消去法的求解过程，以如下4阶方程组的求解为例： (2.5)其展开式为： （2.5）其增广矩阵为： （2.

10、5.0）（1） 消元过程利用若干轮的初等行变换处理，设法将中的A部分化为上三角形式。若主元素，则可保留方程组中的第一个方程，并利用将其余三个方程中的第一个未知量消去，具体做法是取数：再对增广矩阵（2.5.0）中的第i行进行如下初等变换： （其中r表示行或排row）这样原方程的增广矩阵（2.5.0）被变换为如下形式： （2.5.1）其中被改变的元素为：至此，完成了第1轮初等行变换处理。若增广矩阵（2.5.1）中的主元素，则可保留第1个和第2个方程，利用同上的处理方法消去第3个和第4个方程中的第二个未知量，即取数：再对增广矩阵（2.5.1）进行如下初等行变换：得到新的增广矩阵如下： （2.5.2）

11、其中的被改变的元素为：至此，完成了第2轮初等行变换处理。若增广矩阵（2.5.2)中主元素，则保留第13个方程，将第4个方程中的第三个未知量消去，即取数：再对增广矩阵（2.5.2）进行如下初等行变换：得到新的增广矩阵为： （2.5.3）其中被改变的元素为：至此，完成了第3轮初等行变换处理。此时，原4阶方程组通过3轮初等行变换之后，将原方程组（2.5）的增广矩阵中的A部分已化为上三角形式。这一演变过程即为消元过程。（2） 回代过程由线性代数方程组理论可知：经初等行变换之后得到的新的增广阵（2.5.3）所对应的方程组与原方程组（2.5）为等价方程组，即（2.5.3）的解即为原方程组（2.5）的解。由

12、前关于上三角方程组求解的计算格式（2.4a），可依次逆序求得原方程组的解为：逐步求解出的过程即为回代过程。2. Gauss法的求解步骤与计算格式通过以上4阶方程组的求解过程，我们不难归纳出对于n阶方程组Gauss消去法的求解步骤与计算格式。设n阶线性方程组为： （1） 消元过程：通过（n-1）轮初等行变换处理将原方程组变为如下等价的上三角方程组，即： （为上三角方阵）具体做法是：将记为 在第k轮的消元过程中，增广阵中的前k行元素与增广阵中前k行元素保持不变，对其他改变元素和的计算格式为： （2.6）（2） 回代过程：求解由消元过程形成的上三角方程组，即由（2.4a）式的合并表达式，可得解的计算

13、格式为： （2.7）式中约定：当足标j的取值大于其上界n时，和式。3. n阶方程组高斯消去法的计算量（乘除法次数） 在消元过程中（公式（2.6）的计算中），即在第次消元执行过程中需作（）次除法，次乘法，故在整个消元过程中共需乘、除法次数为：。在回代过程中（即公式（2.7）的计算中），共需作次乘、除法。则高斯消元法总的乘除法次数为：其中是二次方以下的余项，可见主要运算量是花费在消去过程中。 例对阶，高斯消元法的乘除法运算次数为，而Cramer法的乘除法运算次数为。显见，高斯消元法的计算量减少的非常显著。§1.3 主元素消去法§1.3.1主元素选取的必要性 在前述的Gauss消

14、元法中，未知量是按其在方程式中的自然顺序逐次消去的；用来消去其他方程式中未知量的方程亦是按自然顺序选取的，故这种方法又称为自然顺序消元法。但实际计算结果表明：这种按自然顺序消元法的方法在某些情况下可能具有如下两个严重的缺陷：（1） 对，在消元过程中有可能出现主元素为零的情况，假设在第k主步的消元过程中，主元素，则有数（除法溢出），尽管，但此时消元过程将无法进行下去，即导致消元过程就此失败。（2） 若有主元素，但其绝对值已很小，此时消元过程仍可以继续下去，但用绝对值很小的主元素作为除数，将会导致其他元素数量级的严重增长和舍入误差的扩散，使得最终获得的解答极不准确。这里给出一个简单的例子。例1.

15、二元线性方程组其精确解为：现取三位浮点十进制数求解（即计算中保留小数后三位有效数字）：1 按自然顺序消元，即有：由此得方程的近似解为：（其结果已完全失真）2 行序交换后再消元（即原两个方程的位置对调），即有：由此得方程的近似解为：（与精确解相当接近）。现分析考察计算过程1产生解失真结果的原因： 令：（解1的误差）将和分别代入1中的第一个方程中，则分别有： 上两式相减得：即有：，这表明的误差传播给误差时被反向扩散了倍。上述例子表明：在消元过程中，适当选取主元素是十分必要的。大量实践表明，若A为对称正定阵，顺序消元法可以保证计算过程对于舍入误差的数值稳定性；但对于一般方阵A，则必须进行选主元处理后

16、，消元法才能得到满意的解答。§1.3.2 主元素消去法 该法仅是在高斯消元法中的每一消去主步之前增加一个主元素的选取过程，其余过程则完全相同。 最常用和最有效的主元素选取方法有两种：1. 列主元素消去法 未知数仍按自然顺序消去，但从各方程要消去的那个未知量的系数中选出模最大者作为下一步消元过程的主元素。 如拟消去未知数时，则在要消去的各方程中找出系数中模最大者，假如它是第（）个方程中的系数，则将第个方程与第k个方程交换位置，则使所选的主元素处于原来主元素的位置上，再按顺序消去法计算公式进行处理。例2 方程组 其增广阵为：回代解得 2. 全主元素选取法 未知数不再按自然顺序消去，而是在

17、每轮消去步之前，从所要消去未知数的所有方程中寻找出未知量系数模最大者作为当前的主元素，并再通过行或列的调换将其置换至原来主元素所在位置处。 在第k步主步消元之前，须先从第k至n个方程中的第k至n个未知数的系数中找出模最大者，假如它是，则将第k个与第个方程交换位置，再将各方程中与的两项位置交换（相当于系数阵A及右端项中的第k行与行交换，A中的第k列与列交换）。此时所选主元素处于原来主元素的位置，然后再按顺序消元法公式进行处理。 例： 两种主元消去法的优缺点：（1） 列主元消去法程序简单，在消去过程中可保持矩阵A的某些特征（带状，特征块型），并有足够的计算精度。（2） 全主元消去法计算精度高，但程

18、序稍复杂，在消去过程中将无法保持矩阵A的某些特征。§1.4 矩阵的分解法前面介绍了高斯消元法，其本质就是将线性方程组的系数矩阵A分解为下三角形方阵和上三角形方阵的乘积。本节首先从矩阵运算的角度来分析消元法，然后给出利用矩阵分解法求解方程组的步骤和计算公式。§1.4.1 消元法与矩阵分解设n阶方程组有解，即（矩阵A非奇异）。 由前知，Gauss消去法的求解过程为：记为 ，or 其中：为上三角形n阶方阵。即通过n-1轮消元主步处理将原方程组变成等价的上三角形方程组，即，而每一轮消元主步相当于对原方程组的增广阵进行了若干次初等行变换。按初等行变换的作用，即为用相应的初等方阵去左乘

19、被变换的矩阵。由此，容易验证： 第1轮消元主步过程相当于对增广阵进行了如下的矩阵乘法运算： ， ，其中，初等变换矩阵，是第1列对角线以下元素非零，其它元素与单位阵相同的下三角方阵。类似地，第2轮消元主步过程可等价表示为如下矩阵乘法：其中，初等变换矩阵是第2列对角线以下元素非零，其它元素与单位阵相同的下三角方阵。······第（n-1）轮消元主步等价为如下矩阵乘法：其中，初等变换矩阵是第(n-1)列对角线以下元素非零，其它元素与单位阵相同的下三角方阵。此时，即为上三角形方阵，为右端项由上述递推关系，显然有： （4.1）其中为第i轮消元主

20、步对应的初等变换矩阵，其表达式为： （4.2）它是第i列对角线以下元素非零，其它元素与单位阵相同的下三角形矩阵。按逆阵的定义，容易得到下三角形矩阵的逆阵为： （4.3）所以，从(4.1)中第一式中可解得： （4.4）其中 为三角方阵 （由矩阵乘法规则容易验证） （4.5）即L为下三角方阵。由（4.4）式可知： （4.4'）这表明原n阶方程组的系数阵A经n-1次高斯消元处理之后即被分解为下三角矩阵L与上三角矩阵的乘积。这一结果揭示了高斯消元过程与矩阵分解之间的本质关系，即为矩阵的三角形分解定理。定理4.1(三角形分解定理) 若n阶方阵A的各阶主子式，则A定可分解为下三角阵与上三角阵的乘积

21、，即：。将，代入方程组中，则原方程组可表为： （4.5）若令： （上三角方程组） （4.6）从而（4.5）可表为： （下三角方程组） （4.7）即原方程组的求解问题亦被相应分解为相继求解（4.7）下三角形系数的方程组和（4.6）上三角形系数的方程组问题。关于（4.7）和（4.6）的计算格式见前公式（2.3）和（2.4）。为表明其正确性，过程结果如下：即先从（4.7）中求出，即，代入（4.6）中，则有，从中解出：。这样处理的好处是显然的，它的突出优点是：在求解具有相同系数矩阵A而右端项不同的方程组时（在结构有限元中，对应同一结构受到不同的荷载工况作用），对系数矩阵A无需再作重新分解处理。如对另一

22、方程组：，顺序求解 即可。§1.4.2 直接求解法由定理4.1，对方程组的系数矩阵进行如下三角形分解：按矩阵乘法规则，A中的元素等于L中的第i行与中的第j列相乘的结果，即：当时（下三角元素），有： （4.8.1）其中约定：当足标大于上界时，上和式恒为零。当时（上三角元素），有： （4.8.2）分别从（4.8.1）式和（4.8.2）式中可获得分解的计算公式为： （紧凑格式） （4.9） 在上式中，由矩阵乘法规则显然有： 1° L的第1列元素即等于A中的第1列元素，即有： （4.9.1） 2° 的第1行元素为： （4.9.2） 不难获得在分解公式（4.9）中计算和的各

23、元素的顺序为： 第1步，由（4.9.1）式求出中的第1列元素； 第2步，由（4.9.2）式求出中的第1行元素； 第3步，由（4.9）第1式求出中的第2列元素； 第4步，由（4.9）第2式求出中的第2行元素； 如此交替进行，最终可求出三角阵和中的所有元素，即实现了对矩阵A的三角形分解。在三角形矩阵的分解公式(4.9)中，无须求出顺序消去法中A阵的中间结果，而可直接由原矩阵A的元素算出和U的元素，故该法称之为直接分解法。相应的分解公式(4.9)亦称为紧凑格式。当的三角形分解完成后，则原方程组可转化为逐次求解如下两三角形方程组，即 (4.7) (4.6) 由计算格式（2.3）式，对下三角形方程组（4

24、.7）其解的计算格式为： （4.10）当求出之后，由计算格式（2.4），对上三角形方程组（4.6）其解的计算格式为： （4.11）可以证明：上述的直接求解算法与高斯消元法是等价的，计算公式（4.10）和（4.11）分别对应于消元法中的消元过程和回代过程。§1.5 对称矩阵的分解法§1.5.1 平方根法（即法）设n阶方程组，若系数矩阵A为对称正定阵，即（对称），（正定），则由矩阵的三角形分解定理4.1，A可被分解为： （5.1） （5.2）为下三角阵，则为上三角阵。将按矩阵乘法规则展开，由出发，按直接分解法类似的步骤，（5.1）式两端元素之间关系为：其中约定：当时，和式； 当

25、时，和式。再由以上两式，很容易获得矩阵元素的计算式为： （5.3）此式即为对称正定矩阵的平方根分解法。在此法中，由于下三角阵的第k列元素即为上三角阵的第k行元素，因此，我们只需顺序地计算出阵的第1列至第n列元素即完成了对A阵的三角形分解，其分解的计算量仅为直接分解法的一半左右。若分解，则对原方程组的求解可分解为如下两三角形方程组的顺序求解:它们的求解公式为： 即将（4.11）式中的元素换为即可。由式（5.3）可见，在平方根分解法中，需要进行n次开方运算，而在计算机上作开方运算需要通过调用函数利用子程序来实现，这一则增大了运算量，二则可能降低了计算精度。为了避免开方运算的出现，人们在法的基础上又

26、研究提出了一种改进的分解法，即Cholesky（乔累斯基）法。§1.5.2 Cholesky（乔累斯基）法（法） 对，若A对称正定，则A可被分解为： （5.4）其中 （5.5）为下三角阵，其主对角元素均为1，其余元素亦与前（5.2）式L阵略有不同（见后）；对角阵 （5.6）事实上，我们只需在前（5.1）的分解表达式中，对阵L稍作数学上的处理即可获得（5.4）式结果。具体作法如下： 从（5.2）表出的下三角阵L的第j列元素中提出因子，由矩阵乘法，则有： （5.7）其中： （5.8）将（5.7）的关系式代入（5.1）式中，则有：（注意：此式L阵与（5.2）中L阵有所不同）仿照前面的处理，

27、可推得D阵和L阵元素的分解计算公式为： （5.9）此处约定：当，时，和式，。对称方阵的这种分解即为乔累斯基分解方法，显见在此分解公式中不再含有开方运算。利用乔累斯基分解法解方程的过程如下： （5.10） （5.11）对三角形方程（5.10）和（5.11）的求解过程和计算格式同前。§1.6 求解三对角方程组的追赶法在利用样条函数进行插值和利用有限差分法求解场问题的数值计算中，常常面临着求解下列形式的三对角方程组：用矩阵表示 （6.1）其中，系数方阵：称为三对角阵。其特点是：所有非零元素都集中在主对角及其相邻的两条次对角线上，而其余元素全为零。针对方程组系数阵的这种带状特点，人们给出高斯

28、消去法的一种简化形式追赶法，其求解步骤：消元过程：第一步：将中的第1个方程中的系数化为1，即有：其中，第二步：注意到在的剩下的方程中，只有第2个方程含有，则利用式将方程中第2个方程中的消去，并将的系数化为1，即有其中，如此一步步顺序加工方程组中的每个方程，设经第步后，方程中的第个方程被化为：其中，当经第步之后，方程中的第个方程被化为：其中，将与中第个方程（最后一个方程）联立，可解得变量为：其中，至此，通过上述消元过程，使原方程组化为如下形式更为简单的等价方程组：即 回代过程：此方程组为二对角形的，其系数阵中的非零元素均集中分布在主对角线以上的两条对角线上，且主对角元素均为1。显然，通过自下而上

29、的逐步回代，即可从方程中依次求得解，其求解公式为：其中，按式计算，按式计算。上述算法即为追赶法，它的消元过程和回代过程分别为追的过程和赶的过程。具体地说，追的过程就是按算式和顺序地求出系数和；赶的过程就是按算式逆序地求出解。可见，追赶法的消元原理与高斯法是相同的，但它针对三对角方程组的具体特点进行了相应的处理，即将系数阵中大量的零元素撇开，从而大大节省了计算量和存贮量，提高了算法的有效性。§1.7 简单迭代方法（Jacobi雅可比法）迭代法的基本思想：将求解方程组归结为重复计算一组彼此独立的线性表达式，产生一个解向量序列，并使该向量序列收敛至某个极限向量，则即为方程组的准确解。迭代法

30、的核心是建立相应的迭代计算格式。对于求解线性方程组常用的方法有：简单迭代法、Seidel（塞德尔）迭代法和松弛法等等。它们之间的差异仅在于迭代格式有所不同。本节介绍简单迭代方法。§1.7.1求解方法和迭代格式以3阶方程组为例说明迭代法的计算过程：设对角主元，将第个方程的第个未知量用其余项表出，则被化为等价形式：其中：式即为原方程组的等价方程，同时它又是简单迭代的计算格式。现任取一组方程组的近似解（通常称为初值或零次近似解）：代入式的右端，即完成了一次迭代，可求出一组新的近似解（一次近似解），即有：将一次近似解再代入式的右端，即完成了第二次迭代，求得二次近似解为：。如法重复上述过程，经

31、次迭代后，得到方程组的次近似解为：式即为简单迭代法的迭代计算公式，若引入矩阵符号：令 则迭代公式可以矩阵表为：（7.4）式中：B迭代矩阵；常数列向量。（7.4）式属于一阶线性定常迭代形式。一阶：新解仅是其前面一次近似解的函数；线性：新解为其前面近似解的线性函数；定常：无论为何值，计算的规则不变（即均为定常矩阵），且与迭代次数无关。按照3阶方程组简单迭代法求解过程和计算格式，容易写出n阶方程组的简单迭代公式为：等价变形为（从中得的形式解表达式）：其中：将缩写为：或统一表为：从上式可直接写出简单迭代法的迭代计算格式（表达式）为： 的矩阵形式完全同（7.4）式，即：其中，例1.（华中计算方法P208

32、）用简单迭代法求解：依次从三个方程中分离出，建立迭代计算格式：（a）取初值，迭代结果如下：000010.720.830.8420.9711.071.1531.0571.15711.248281.09981.199411.2997891.099941.199941.29992已趋向于本问题的准确解：。可见由迭代格式（a）所获得的解序列最终将收敛于准确解。§1.7.2 迭代的收敛条件简单迭代法虽然求解过程简单，但并不是随便构造一个迭代计算格式就能实现的，其中有一个非常重要的理论问题在迭代中必须予以关注，这就是所谓的收敛性问题。现再一次考察例1，与前处理不同，先将分别从原方程组（a）的第2

33、，3和1个方程中分离出来，并建立相应的迭代格式，即：（a）仍取初值迭代结果： 若继续迭代下去，所得解向量的值会越来越大，即迭代格式（a）是发散的（不收敛）。这里首先给出迭代收敛性的有关数学描述。收敛性定义：对迭代格式（7.10）当序列极限存在，则称该迭代格式为收敛的，且即为原方程组的解。在以上两例中，为何对同一方程组，构造的两个迭代计算表达式有所不同，会导致迭代结果收敛和发散两种截然不同的结果？其原因何在？进一步，欲确保迭代格式具有收敛性的条件是什么？收敛条件讨论：将准确解代入迭代公式中，必有两端相等，即由式-，得：显然有：令：即将矩阵中各行向量一阶范数中的最大者记为。（注：假如向量为：，则其

34、一阶范数为：） （7.14a） （7.14b）即准确解与第次近似解和第次近似解之间的最大误差。由（7.14a）和（7.14b）式，则式可表为： 此式即为第次迭代误差与第次迭代误差之间的关系不等式。由此可知，若，即使，必有当，即迭代解向量的计算误差将随迭代进程趋于零。这表明此迭代过程（格式）是收敛的。由此，可得迭代收敛性的判别定理如下：定理7.1 在迭代格式中，若，则该迭代格式对任意给定的初值均是收敛的。这一定理在实际中的应用是特别简单和方便的。这里，我们用此定理判别例1中两种迭代格式的收敛性：对迭代式（a）：（收敛）对迭代式（a）：（不收敛）这一判别结果与前面的数值计算结果是完全一致的。

35、67;1.8 Gauss-Seidel（高斯-塞德尔）迭代法与松弛法简单迭代法虽然具有迭代格式特别简单的优点，但其迭代的收敛速度较慢。为了加快迭代收敛速度，人们在简单迭代法的基础上研究提出了几种改进的迭代方案，主要有Gauss-Seidel迭代法和松弛法。§1.8.1 Gauss-Seidel法在简单迭代格式式中将第一式求得的阶近似解代替第二式右端的第阶近似解去求得近似解；再将、代替第三式右端的、去求出，逐次这样处理直到第n个方程。第轮迭代的计算式为：它是在简单迭代法的计算过程中，逐次以新值代替老值再进行迭代。式即是所谓的Gauss-Seidel迭代法的计算公式。为了写出其一般迭代公

36、式，先不妨详细写出中的第个表达式如下：将其中新值部分和老值部分分别用和式形式表为，则有： 若将的关系式代入，则式可表为： （8.2）式或（8.2）式即为Gauss-Seidel迭代法的计算公式。由于新值比其老值更准确一些，显然用这些已知的新值去代替老值进行迭代，自然迭代收敛于其准确解的速度要比简单迭代法来的快。为了进行比较，仍以上节例1为例。例3 方程组同上节例1，将例1的迭代格式改写为Gauss-Seidel迭代格式如下：其中迭代矩阵与简单迭代法中的完全相同。取初值，迭代结果如下：000010.720000.902001.1644021.043081.167191.2820561.09999

37、1.199991.30000此结果与前简单迭代结果相比可见，Seidel法的第6轮迭代结果的精度已高出简单迭代法第9次迭代结果的精度，即Seidel法的迭代收敛速率明显加快。显然，在Seidel迭代公式中，其迭代矩阵与简单迭代法中的完全相同。故关于简单迭代法的收敛性及其判别定理7.1亦完全适用于Seidel迭代法。§1.8.2 松弛法松弛法是另一种线性加速的方法。其基本作法：将前一步的结果与Seidel方法的迭代值适当进行线性组合，构造一个收敛速度较快的近似解序列，其迭代计算过程的每一步分为迭代和加速两个步骤。其迭计算公式为：迭代：加速：其中，系数称为松弛因子（即新解与老解的线性组合

38、系数）。我们可将迭代公式和加速公式合并在一起，则松弛法的迭代计算格式为：这就是松弛法的计算公式。显然，当取松弛因子时，式就变成了Gauss-Seidel法迭代公式。因此，适当选取因子可得到高于Gauss-Seidel法的收敛速度。其加速作用的简单解释为：若迭代收敛，则式的计算结果定比更精确些，因此在加速公式中应加大在线性组合中的比重，尽可能扩大它的影响效果。为了达到这一目的，通常应取松弛因子。但的值并不是越大越好，在实际计算中，可根据系数矩阵的性质，结合经验通过多次试算来选定最佳的松弛因子之值。数学上可以证明：为保证迭代格式式收敛，须要求：。当松弛因子，称为超松弛法；当松弛因子，即为Gauss

39、-Seidel法；当松弛因子，称为低松弛法。在松弛法中，关于“松弛”一词的数学作用，以及松弛因子取值范围的证明可参见计算方法（武大、山大P69）。因涉及到矩阵谱半径和谱范数等概念及其有关的定理，故在此不展开讨论。Chapter 2插值与数值微分§2.1 插值问题概述§2.1.1 背景工程上：在机械制造、产品几何造型中，常常给定一批离散的样点，要求利用数学的方法，自动构造一条光滑连续的曲线（面），使该曲线（面）通过这些样点，以满足几何造型设计要求，或者据此进行机械加工或放样。例如：赋形抛物面天线（非标准的抛物面），汽车、飞机的外壳等等。数学上：在利用数值方法求解有关连续函数问

40、题时，往往需要根据给定离散样点的数据插补中间点。如：用差商代替微商；结构有限元中求非节点处的位移和单元非形心处的应力等等。与此相关的数学问题就是插值问题。简单通俗地讲，所谓插值就是在给定的样点间再插进一些所需的中间值。§2.1.2 数学描述给定一组样点：，设法构造一个光滑连续的函数y=F（x），使该函数曲线（面）通过所有样点，即满足： （*）如图示： 满足插值条件（*）的函数F(x)即为所求的插值函数。构造F（x）的基本要求：（1）函数形式应为便于计算的初等函数，如多项式函数、分段多项式函数、样条函数等等；（2）F(x)中的自由度（即待定参数的个数）应与插值条件（*）的个数相等。为以

41、下讨论方便起见，统一约定如下：样点的节点；插值点；某个原函数（通常未知）；函数在节点处的函数值和一阶导数值（已知）；待求的插值函数（待求）；原函数与插值函数在点处的插值余项（绝对误差）；为含有点（）的最大区间。若插值点，则该插值问题称为内插；若插值点，则该插值问题称为外推。§2.2 Lagrange（拉格朗日）插值若给定样点，使（过点），则称为Lagrange（拉格朗日）插值函数。下面首先介绍拉氏插值中两种最简单的情况，即线性插值和抛物线插值，然后在此基础上给出高次拉氏插值函数的表达式。§2.2.1 线性插值（两点一次插值）。求：过两样点（x0，f0），（x1，f1）的线性

42、插值函数y=F（x）。由解析几何关系，过、两点的线性函数（直线方程）可表为： (点斜式) （2.1a）上式经整理，则有： (对称式) （2.1b）为今后表述方便且统一起见，令： (2.2)这是两个的一次多项式函数。显然有： （2.3），在节点、处构成了正交规范基，故称之为线性插值函数的基函数。从而（2.1b）式给出的线性插值函数可由基函数表为： （2.1c）此式表明：插值函数可表为基函数的线性组合，即曲线是由基函数张起构成的。利用微分中值定理，可推得线性插值函数的插值余项为： （2.4）式中：为微分中值点；=为的二次多项式。可见：线性插值函数对于原函数具有二阶精度，即对于原函数的一阶导数值亦有

43、逼近性。当原函数为线性函数，则插值余项，即此时插值函数就是原函数。§2.2.2 抛物线插值（三点二次插值）求过三样点（x0，f0），（x1，f1），（x2，f2）的二次多项式（抛物线）插值函数。由前知，所求插值函数可以表为其对应基函数的线性组合形式，故不妨假设： （2.5）其中，为待求的基函数，它们应满足基函数的正交规范化条件，即有： （2.6）这样的基函数不难直接构造，现以为例说明之。欲使满足，显然在的表达式中应含有和两个因子，故令： （2.7）其中为待求系数。由条件，解得：代入（2.7），则有 （2.8a）这是一个二次的多项式函数。同理可得： （2.8b） （2.8c）将（2.8

44、）式统一表为： （2.8）将（2.8ac）代入（2.5）式中，即得所求插值函数F(x)。易验证F(x)满足插值条件：利用微分中值定理，二次插值函数的插值余项为： （2.9）即具有三阶精度，直到二阶导数都有逼近性。 例1 利用100,121,144的平方根求（精确值=10.7238*）之值。解：已知样点为：（x0=100，f0=10），（x1=121，f1=11），（x2=144，f2=12），插值点：x=115，求插值函数值F（115）。代入插值公式（2.5）中，即有：可见误差甚小。§2.2.3 n次多项式插值求过（n+1）个样点（xi，fi）（i=0,1,2,n）的n次多项式插值函

45、数。设： （2.10）由前讨论可知，该n次多项式函数共用（n+1）项，它可由其（n+1）个基函数的线性组合形式表出，即有： （2.11）其中（n+1）个基函数必须满足正交规范条件，即： （2.12）由此条件仿前处理，可推得基函数的表达式为： （2.13）将（2.13）式代入（2.11）式中，得n次多项式插值函数的计算表达式为： （2.14）可见，拉氏插值公式在逻辑结构上表现为二重循环，其中内循环是由累积获得基函数，外循环由节点的函数值与基函数之积的累加获得插值点处的函数值。利用微分中值定理，可得n次拉氏插值的余项为： （2.15）该插值函数具有（n+1）阶精度。例2 已知原函数的四个样点值为：

46、0.10.150.250.30.9048730.8607080.7788010.740818求插值点处的f（x）的近似值。（精确值）解：因为给定四个样点，故采用四点三次拉氏插值公式。由（2.11）式，即有： （a）由（2.14）式，其中基函数为： （b）将和x=0.2分别代入（b）中，得基函数在x=0.2处的值为：由（a）式x=0.2处的插值函数值为： 由插值余项公式（2.15）式，则有：从而在插值点处的余项的估计为：由于为单减函数，故。从而得的估计区间为：在处插值余项的精确值为：落入误差的估计区间内。§ 2.3 Hermite（埃尔米特）插值在许多情况下，不仅要求插值函数过所有的样

47、点，甚至还要求与在样点处的导数值彼此相等，即要求既“过点”还要求“相切”。这就是所谓的埃尔米特（Hermite）插值问题。其数学表述为：若使：， 且使： 则称为Hermite插值函数。本节首先介绍埃尔米特插值中两种最简单情况，即一点一次带导数插值和两点三次带导数插值，然后再给出一般高阶埃氏插值的计算公式。§2.3.1 一点一次带导数插值（切线插值）给定一样点为，其中。求一次插值多项式函数F(x)，满足：， 根据以上两个插值条件，由解析几何，构造为： （3.1）显然有： 利用微分中值定理，其插值余项为： （3.2）该插值函数具有一阶导数的逼近性。§2.3.2 两点三次带导数插

48、值给定两样点：，求既过点又相切的多项式插值函数。即插值条件为： （*）共有四条，故设F(x)为三次多项式函数（共有四项），即 （3.3）将（*）式代入上式，得以a0，a1，a2，a3为未知量的线性方程组为：联立解得a0，a1，a2，a3，代入插值函数（3.3）式中，并经整理得的统一表达式为： （3.4）其中：即为基函数，它们表达式为： （3.5）其中：即为两点线性插值的基函数；为节点步长。易验证： 这是一组四维正交规范基，这表明：所求插值函数（3.4）式即为基函数的线性组合。这与拉氏插值的结果是一致的，差别仅在于基函数有所不同。由微分中值定理，（3.4）式的插值余项为： （3.7）=该插值函数

49、具有三阶导数的逼近性。对于基函数（3.5）的另一种表示因: x的二次多项式 (3.6)则有： 从而基函数亦可表为： （3.5a）此种形式便于推广到高次埃尔米特插值多项式函数的基函数表达式之中。§2.3.3 高阶埃氏插值给定(n+1)个样点，求既过点又相切的多项式插值函数。插值条件为： （3.8）共有个，故插值函数应为次多项式函数。仿照前两点三次带导数的埃氏函数的产生方法，这个次多项式插值函数可表为： （3.9）其基函数，（），可按（3.5a）式的格式表为： （3.10）其中 （3.11） （3.12）其插值余项为： （3.13）具有阶导数的逼近性。埃氏插值的几何意义是：插值函数()与原函数(x)在节点处即“等值”又“公切”。埃氏插值还可推广到包含更高阶导数以及各节点包含的导数阶不均等的情况。§1.4 关于多项式插值的稳定性根据拉氏和埃氏插值的余项公式(2.15)和(3.13)，似乎表明插值函数的次数愈高则精度愈高，与原函数的贴近效果就愈好。其实并非如此，现分析如下：在插值过程中总的误差来自于以下两个方面：1. 原函数f（x）被插值函数F（x）代替所引起的原理性误差，即余项E(x)=f(x)-F(x)（截断误差），其误差分析即为余项公式。2. 节点数据fi本身的误差（实验误差、或计