数值计算方法(免费版)

1、数值计算方法第1章绪 论1.1数值计算方法的对象与特点1.1.1 什么是数值计算方法现代的科学技术发展十分迅速，他们有一个共同的特点，就是都有大量的数据问题。比如，发射一颗探测宇宙奥秘的卫星，从卫星设计开始到发射、回收为止，科学家和工程技术人员、工人就要对卫星的总体、部件进行全面的设计和生产，要对选用的火箭进行设计和生产，这里面就有许许多多的数据要进行准确的计算。发射和回收的时候，又有关于发射角度、轨道、遥控、回收下落角度等等需要进行精确的计算。有如，在高能加速器里进行高能物理试验，研究具有很高能量的基本粒子的性质、它们之间的相互作用和转化规律，这里面也有大量的数据计算问题。计算问题可以数是现

2、代社会各个领域普遍存在的共同问题，工业、农业、交通运输、医疗卫生、文化教育等等，各行各业都有许多数据需要计算，通过数据分析，以便掌握事物发展的规律。研究计算问题的解决方法和有关数学理论问题的一门学科就叫做计算方法。计算方法属于应用数学的范畴，它主要研究有关的数学和逻辑问题怎样由计算机加以有效解决。1.1.2 数值计算方法的内容数值计算方法也叫做计算数学或数值分析。数值计算方法主要内容包括非线性方程求根、线性代数方程组解法、微分方程的数值解法、插值问题、函数的数值逼近问题、概率统计计算问题等等，还要研究解的存在性、惟一性、收敛性和误差分析等理论问题。我们知道五次及五次以上的代数方程不存在求根公式

3、，因此，要求出五次以上的高次代数方程的解，一般只能求它的近似解，求近似解的方法就是数值分析的方法。对于一般的超越方程，如对数方程、三角方程等等也只能采用数值分析的办法。怎样找出比较简洁、误差比较小、花费时间比较少的计算方法是数值分析的主要课题。在求解方程的办法中，常用的办法之一是迭代法，也叫做逐次逼近法。迭代法的计算是比较简单的，是比较容易进行的。迭代法还可以用来求解线性方程组的解。求方程组的近似解也要选择适当的迭代公式，使得收敛速度快，近似误差小。在线性代数方程组的解法中，常用的有塞德尔迭代法、共轭斜量法、超松弛迭代法等等。此外，一些比较古老的普通消去法，如高斯法、追赶法等等，在利用计算机的

4、条件下也可以得到广泛的应用。在数值计算方法中，数值逼近也是常用的基本方法。数值逼近也叫近似代替，就是用简单的函数去代替比较复杂的函数，或者代替不能用解析表达式表示的函数。数值逼近的基本方法是插值法。初等数学里的三角函数表，对数表中的修正值，就是根据插值法制成的。在遇到求微分和积分的时候，如何利用简单的函数去近似代替所给的函数，以便容易求微分和求积分，也是计算方法的一个主要内容。微分方程的数值解法也是近似解法。常微分方程的数值解法有欧拉法、预测校正法等。偏微分方程的初值问题或边值问题，目前常用的是有限差分法、有限元素法等。有限差分法的基本思想是用离散的、只含有限个未知数的差分方程去代替连续变量的

5、微分方程和定解条件。求出差分方程的解法作为求偏微分方程的近似解。有限元素法是近代才发展起来的，它是以变分原理和剖分差值作为基础的方法。在解决椭圆形方程边值问题上得到了广泛的应用。目前，有许多人正在研究用有限元素法来解双曲形和抛物形的方程。数值计算方法的内容十分丰富，它在科学技术中正发挥着越来越大的作用。1.1.3 数值计算方法的特点 第一，面向计算机，要根据计算机特点提供实际可行的有效算法。即算法只能包括加、减、乘、除和逻辑运算，是计算机能直接处理的。第二，有可靠的理论分析，能任意逼近并达到精度要求，对近似算法要保证收敛性和数值稳定性，还要对误差进行分析。第三，要有好的计算复杂性，时间复杂性好

6、是指节省时间，空间复杂性好是指节省存储量，这也是建立算法要研究的问题，它关系到算法能否在计算机上实现。第四，要有数值实验，即任何一个算法除了从理论上要满足上述三点外，还要通过数值实验证明是行之有效的。根据“数值计算方法”的特点，学习时我们要注意掌握方法的基本原理和思想，要注意方法处理的技巧及其与计算机结合，要重视误差分析、收敛性及稳定性的基本理论；其次，要通过例子，学习使用各种数值方法解决实际计算问题；最后，为了掌握本课的内容，还应做一定数量的理论分析与计算问题练习。由于本课内容包括了微积分、代数、常微分方程的数值方法，以及高级程序设计的方法，读者必须掌握这几门课的基本内容才能学好这一课程。1

7、.2误差来源与误差分析1.2.1 误差的来源早在中学我们就接触过误差的概念，如在做热力学实验中，从温度计上读出的温度是23.4就不是一个精确的值，而是含有误差的近似值。事实上，误差的在我们的生活中无处不在，无处不有。如量体裁衣，量与裁的结果都不是精确无误的，都有误差。人们可能会问：如果使用计算机来解决问题，结果还会有误差吗？下面我们通过考察数学方法解决实际问题的主要过程来思考这个问题。用数学方法解决一个具体的实际问题，首先要建立数学模型，这就要对实际问题进行抽象、简化，因而数学模型本身总含有误差，这种误差叫做模型误差。在数学模型中通常包含各种各样的参变量，如温度、长度、电压等，这些参数往往都是

8、通过观测得到的，因此也带来了误差，这种误差叫做观测误差。当数学模型不能得到精确解时，通常要建立一套行之有效的数值方法求它的近似解，近似解与准确解之间的误差就称为截断误差或方法误差。由于在计算机中浮点数只能表示实数的近似值，因此用计算机进行实际计算时每一步都可能有误差，这种误差称为舍入误差。例如，函数f(x)用泰勒(Taylor)展开式近似代替，则数值方法的截断误差是 （介于0与之间）又如，在计算时用3.14159近似替代产生的误差-3.14159=0.0000026 就是舍入误差。上述种种意味着都会影响计算结果的准确性，因此需要了解与研究误差。在数值计算方法中将着重研究截断误差、舍入误差，并对

9、它们的传播与积累作出分析。1.2.2 绝对误差、相对误差与有效数字1绝对误差与绝对误差限定义1.1 设为准确值，为x的一个近似值，称为近似值的绝对误差，或误差。通常我们无法知道准确值，也不能算出误差的准确值，只能根据测量或计算估计出误差的绝对值不超过某个正数，则称为绝对误差限。有了绝对误差限就可知道x的范围，即落在区间内。例如用毫米测度尺测量一长度，读出的长度为23mm，则有mm。由此例也可以看到绝对误差是有量纲和单位的。2相对误差与相对误差限只用绝对误差还不能说明数的近似程度，例如甲打字时每百个字错一个，乙打字时平均每千个字错一个，他们的误差都是错一个，但显然乙要准确些。这就启发我们除了要看

10、绝对误差大小外，还必须顾及量的本身。定义1.2 把近似值的误差与准确值的比值称为近似值的相对误差，记作。实际计算时，由于真值通常是不知道的，通常取。相对误差也可正可负，它的绝对值的上界叫做相对误差限。记作。即。根据定义，甲打字时的相对误差，乙打字时的相对误差。易知相对误差是一个无量纲量。3有效数字当有多位时，常常接四舍五入的原则得到的前几位近似值，例如取3位；取5位；它们的误差都不超过末位数字的半个单位，即现在我们将四舍五入抽象成数学语言，并引入一个新名词“有效数字”来描述它。定义1.3 若近似数的误差限是某一位的半个单位，该位到的第一位非零数字共有位，我们就说有位有效数字。如取作的近似值，就

11、有3位有效数字；取作的近似值，就有5位有效数字。有效数字也可采用以下定义：x的近似数x*写成标准形式： （1.1）其中是0到9的一个数字，且不为0，为整数，若（1.2）则称有位有效数字。例1.1 依四舍五入原则写出下列各数具有5位有效数字的近似数，913.95872，39.1882，0.0143254，8.000033解: 按定义，上述的5位有效数字的近似数分别为 913.96，39.188，0.014325，8.0000注意，8.000033的5位有效数字的近似值是8.0000而不是8，8只有一位有效数字。4有效数字与误差关系不难看出，若由（1.1）给出某近似数有 位有效数字，则可以从（1.

12、2）求得这个数的绝对误差限因此，在相同的情况下，越大则就越小，故有效数字越多，绝对误差限越小。定理1.1 用（1.1）表示的近似数，若具有位有效数字，则其相对误差限为反之，若的相对误差限则至少具有位有效数字。证明 由式（1.1）得所当有位有效数字时反之，若的相对误差限有：由式（1.2）式知道，至少有位有效数字，证毕这说明有效数字越多，相对误差限越小。例1.2 要使的近似值的相对误差限小于，要取几个有效数字？解：因为，所以，令故取即可满足。5函数计算的误差估计一般情况，当自变量有误差时计算相应的函数值也会产生误差，其误差限可利用函数的泰勒展开式估计。设是一元函数，的近似值为，以近似，其误差界记作

13、，可用泰勒展开式介于，之间，取绝对值得假定与的比值不太大，可忽略的高阶项，于是可得计算函数值的误差限例1.3 设的相对误差限为2%，则的相对误差限是多少？解：所以的相对误差限为0.02n。当为多元函数时，如计算。如果的近似值为，则的近似值，于是函数值的误差由泰勒展开式得于是误差限为例1.4 已测得某场地长的值为，宽的值为，已知，试求面积的绝对误差限与相对误差限。解：因，那么其中于是绝对误差限相对误差限1.3误差传播与若干防治办法由前述可知，在数值计算中每步都可能产生误差。而一个问题的解决，往往要经过成千上万次运算，我们不可能每步都有加以分析，下面，通过对误差的某些传播规律的简单分析，指出在数值

14、计算中应注意的几个原则，它有助于鉴别计算结果的可靠性并防止误差危害现象的产生。1要避免两相近数相减在数值运算中两相近数相减会使有效数字严重损失，例如都具有五位有效数字，但只有两位有效数字，所以要尽量避免这类运算。通常采用的方法是改变计算公式，例如当与很接近时，由于那么可用右端的公式代替左端的公式计算，有效数字就不会损失。当很大时，也可用右端来代替左端。一般情况，当时，可用泰勒展开取右端的有限项近似左端。如果计算公式不能改变，则可采用增加有效位数的方法。2要防止大数“吃掉”小数若参加运算的数的数量级相差很大，而计算机的位数有限，如不注意运算次序，就可能出现大数“吃掉”小数的现象，影响计算结果。例

15、如在五位十进制计算机上，计算写成规格化形式由于计算时要对阶，在计算机中表示为0，计算出来，结果严重失真！如果计算时，先将计算出来，再与52492相加，就不会出现大数“吃掉”小数的现象了。3注意简化计算步骤，减少运算次数同样一个计算题，如果能减少运算次数，不但可节省计算机的计算时间，还能减少舍入误差。例如，计算x255的值，如果逐个相乘要用254次乘法，但若写成只要做14次乘法运算即可。4绝对值太小的数不宜作除数设与分别有近似值的近似值，其绝对误差显然，当很小时，近似值的绝对误差有可能很大。因此，不宜把绝对值太小的数作除数。第2章非线性方程求根与线性方程相比，非线性方程问题无论是从理论上还是从计

16、算公式上，都要复杂得多。对于一般的非线性方程，计算方程的根既无一定章程可寻也无直接法可言。例如，求解高次方程组的根，求解含有指数和正弦函数的超越方程的零点。解非线性方程或非线性方程组也是计算方法中的一个主题。一般地，我们用符号来表示方程左端的函数，方程一般形式表示 为 ，方程的解称为方程的根或函数的零点。通常，非线性方程的根不止一个，对于非线性方程，一般用迭代法求解。因此，在求解非线性方程时，要给定初始值或求解范围。2.1实根的对分法2.1.1 使用对分法的条件对分法或称二分法是求方程近似解的一种简单直观的方法。设函数在上连续，且，则在上至少有一零点，这是微积分中的介值定理，也是使用对分法的前

17、题条件。计算中通过对分区间、缩小区间范围的步骤搜索零点的位置。例2.1 用对分法求解在区间之间的根。解：（1），由介值定理可得有根区间。（2）计算，有根区间。（3）计算(1.5+2)/2=1.75，=0.078125，有根区间。（4）一直做到（计算前给定的精度）或时停止。2.1.2 对分法求根算法计算的一般计算步骤如下：1输入求根区间和误差控制量，定义函数。IF THEN 退出选用其他求求根方法2WHILE 时（1）计算中点以及的值；（2）分情况处理： 停止计算，转向步骤4修正区间 修正区间 ENDWHILE3。4输出近似根。图2.1给出对分法示意图图2.1 对分法示意图在算法中，常用代替的判

18、断，以避免数值溢出。对分法的算法简单，然而，若在是有几个零点时，只能算出其中一个零点；另一方面，即使在上有零点，也未必有。这就限制了对分法的使用范围。对分法只能计算方程的实根。2.2迭代法对给定的方程，将它转达换成等价形式：。给定初值，由此来构造迭代序列，如果迭代法收敛，即，有，则就是方程的根。在计算中当小于给定的精度控制量时，取为方程的根。例如，代数方程的三种等价形式及其迭代格式如下：（1）迭代格式 （2），迭代格式 （3）迭代格式 对于方程构造的多种迭代格式，怎样判断构造的迭代格式是否收敛？收敛是否与迭代的初值有关？定理2.1 若定义在 上，如果满足（1）当有有（2）在上可导，并且存在正数

19、，使对任意的，有；则在上有惟一的点满足，称为的不动点。而且迭代格式对任意的初值均收敛于的不动点，并有误差估计式(2.1)证明：（1）先证明存在性：令，则有故有，使得或（2）再证明惟一性：设都是的不动点，且，则有 与假设矛盾，这表明 ，即不动点是惟一的。（3）当时，由于可用归纳法证明，迭代序列，于是由微分中值定理，和，得 (2.2)因为，所以当时，即迭代格式收敛。（4）误差估计式：设固定，对于任意的正整数有，由于的任意性及，故有 注：定理2.1是判断迭代法收敛的充分条件，而非必要条件。要构造满足定理条件的等价形式一般是难于做到的。事实上，如果 为的零点，若能构造等价形式而，由的连续性，一定存在的

20、邻域，其上有，这时若初值迭代也就收敛了。由此构造收敛迭代式有两个要素，其一，等价形式应满足；其二，初值必须取自的充分小邻域，这个邻域大小决定于函数，及做出的等价形式。例2.2 求代数方程，在附近的实根。解：1），当构造的迭代序列收敛。取 则=2.08008, =2.09235, =2.094217=2.094494, =2.094543, =2.094550准确的解是=2.094551481502）将迭代格式写为，当迭代格式不能保证收敛。2.3迭代法的加工对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此迭代过程的加速是个重要的

21、课题。以下介绍一种埃特金（Aitken）方法。对方程，构造加速过程，算法如下：（1）预测：（2）校正：（3）改进：有些不收敛的迭代法，经过埃特金方法处理后，变得收敛了。例2.3 求方程在附近的根。解：若采用迭代公式： 迭代法是发散的，我们现在以这种迭代公式为基础形成埃特金算法：取计算结果如表2.1所示：表2.1 计算结果 0123452.375001.840921.491401.347101.3251813.39655.238882.317281.444351.327141.51.416291.355651.328951.324801.32472我们看到，将发散的迭代公式通过埃特金方法处理后，

22、竟获得了相当好的收敛性。2.4牛顿迭代法对方程可构造多种迭代格式，牛顿迭代法是借助于对函数作泰勒展开而构造的一种迭代格式。将在初始值作泰勒展开： 取展开式的线性部分作为的近似值，则有 设，则 令 类似地，再将在作泰勒展开并取其线性部分得到： 一直做下去得到牛顿迭代格式：,(2.3)牛顿迭代格式对应于的等价方程是 (2.4)若是的单根时，则有，只要初值充分接近，有，所以牛顿迭代收敛。当为的重根时，取下面迭代格式： 牛顿法的几何意义以为斜率作过点的切线，即作在点的切线方程 令，则得此切线与轴的交点，即 再作在处的切线，得交点，逐步逼近方程的根。如图2.2所示图2.2 牛顿切线法示意图在区域的局部“

23、以直代曲”是处理非线性问题的常用手法。在泰勒展开中，截取函数展开的线性部分替代。例2.4 用牛顿迭代法求方程，在附近的根。解： 计算结果列于表2.2中。表 2.2 计算结果Kxk f (x)0123451.00 -2.81.41176 -0.7270711.62424 -0.1454931.6923 -0.01316821.69991 -0.00015151.7 0比较表2.1和表2.2的数值，可以看到牛顿迭代法的收敛速度明显快于对分法。牛顿迭代法也有局限性。在牛顿迭代法中，选取适当迭代初始值是求解的前题，当迭代的初始值在某根的附近时迭代才能收敛到这个根，有时会发生从一个根附近跳向另一个根附近

24、的情况，尤其在导数数值很小时，如图2.3所示。图2.3 失效的牛顿迭代法牛顿迭代算法1定义函数，输入或定义迭代初始值和控制精度epsilon。2FOR to MAXREPT THEN输出满足给定精度的近似解，结束 ENDFOR3输出：在附近无根。注：在伪码中，“/”表示注释语句。伪码又称过程设计语言，主要用于描述算法。它是某种高级语言和自然语言的混杂语言，它取某种高级语言中的一些关键字，用于描述算法的结构化构造和数据说明等。伪码的语句中嵌有自然语言的叙述，伪码易于2.5弦截法弦截法迭代格式在牛顿迭代格式中：用差商代替导数，并给定两个初始值和，那么迭代格式可写成如下形式： (2.5)上式称为弦截

25、法。用弦截法迭代求根，每次只需计算一次函数值，而用牛顿迭代法每次要计算一次函数值和一次导数值。但弦截法收敛速度稍慢于牛顿迭代法。弦截法的几何意义做过两点和的一条直线（弦），该直线与轴的交点就是生成的迭代点，再做过和的一条直线，是该直线与轴的交点，继续做下去得到方程的根，如图2.4所示图2.4 弦截法示意图例2.5 用弦截法求方程的根，取。解： 计算结果列于表2.3中。表2.3 计算结果 012345671.5 -0.454 2.31.90909 0.2488351.65543 -0.08056921.71748 0.02874561.70116 0.001959021.69997 -0.000

26、05392461.7 9.45910-8弦截法算法1定义函数，输入或定义控制精度值epsilon和迭代初始值。计算 表示 2FOR:=0,1MAXREPT2.1 2.2 /表示2.3 THEN输出满足给定精度的近似解，结束2.4 /为下一次迭代准备数值ENDFOR3输出：在初始值附近无根。2.6非线性方程组的牛顿方法为了叙述简单，我们以解二阶非线性方程组为例演示解题的方法和步步骤，类似地可以得到解更高阶非线性方程组的方法和步骤。设二阶方程组 (2.6)其中为自变量。为了方便起见，将方程组写成向量形式： ，其中= 将附近作二元泰勒展开，并取其线性部分，得到方程组 令则有如果，解出,得 再列出方程

27、组：解出，得出 继续做下去，每一次迭代都是解一个方程组 (2.7) 即，直到为止。例2.6 求解非线性方程组取初始值。解： 解方程得 继续做下去，直到时停止。将两个变量的非线性方程组推广到多变量的非线性方程组： （2.8）记的向量形式： 用牛顿迭代法求方程组的解为： 其中： 称为雅可比（jacobi）矩阵。计算中采用 (2.9)一直做到小于给定精度为止。在的领域中若，则迭代收敛。.7程序示例程序2.1 用牛顿迭代法求在附近的根。算法描述给定，从开始，根据牛顿迭代公式 计算在附近的根。程序源码/ Purpose :牛顿迭代法求根 /#includestdio.h#includemath.h /i

28、terate(x)为 #define iterate(x) (x-(x*x-3)*x-1)/(3*x*x-3)#define (x*x-3)*x-1)#define x0 1.5 /迭代初值x0#define MAXREPT 1000 /最大迭代次数#define epsilon 0.0001 /精度void main( ) int i; double x_k=x0,x_k1=x0 for(i=iMAXREPT；i+) printf(“Got%f/n”,x_k1);x_k1=iterate(x_k); /迭代 if fabs(x_k1 x_k)epsilon|fabs(f(x_k1)epsil

29、on printf(“!Root:%fn”,x_k1);/满足精度，输出 return; x_k=x_k1 printf(“After % d repeate, no solved.m”,MAXREPT)；/-End of File-计算实例已知，取，用牛顿迭代法计算的根。程序输入输出由本程序预设的及，得到输出！ Root: 1.879385程序2.2 用弦截法求在附近的根。算法描述给定，从开始，根据弦截法迭代公式 求得在其附近的根。程序源码/ Purpose：弦截法求根 /#includestdio.h#includemath.h#define f(x) (x*x*x 7.7*x*x+19.

30、2*x 15.3)/f函数#define x0 0.0 /初值x0,x1#define x1 1.0 #define MAXREPT 1000 /最大迭代次数#define epsilon 0.00001 /求解精度void main( )int I; double x_k=x0,x_k1=x1,x_k2=x1; for(I=0;iMAXREPT;i+ printf(“!Root:%fn”,x_k2); x_k2=x_k1-(f(x_k1)*x_k1-x_k)/(f(x_k1)-f(x_k);/弦截求新x_n if fabs(x_k2-x_k1) epsilon|fabs(f(s_k2) ep

31、silon printf(“!Root:%fn”,x_k2); /满足精度，输出 return; x_k=x_k1; x_k1=x_k2; /反复 printf(“After %d repdate, no solvde.n”,MAXREPT); /-End of File-计算实例用弦截法求方程的根，取 程序输入输出由本程序的及得到输出!Root: 1.700000第3章解线性方程组的直接法在近代数学数值计算和工程应用中，求解线性方程组是重要的课题。例如，样条插值中形成的关系式，曲线拟合形成的法方程等，都落实到解一个元线性方程组，尤其是大型方程组的求解，即求线性方程组（3.1）的未知量的数值。

32、（3.1）其中ai j，bi为常数。上式可写成矩阵形式Ax = b，即 （3.2）其中，为系数矩阵，为解向量，为常数向量。当detA=D0时，由线性代数中的克莱姆法则，方程组的解存在且惟一，且有 为系数矩阵的第列元素以代替的矩阵的行列式的值。克莱姆法则在建立线性方程组解的理论基础中功不可没，但是在实际计算中，我们难以承受它的计算量。例如，解一个100阶的线性方程组，乘除法次数约为（101100！99），即使以每秒的运算速度，也需要近年的时间。在石油勘探、天气预报等问题中常常出现成百上千阶的方程组，也就产生了各种形式方程组数值解法的需求。研究大型方程组的解是目前计算数学中的一个重要方向和课题。解

33、方程组的方法可归纳为直接解法和迭代解法。从理论上来说，直接法经过有限次四则运算，假定每一步运算过程中没有舍入误差，那么，最后得到方程组的解就是精确解。但是，这只是理想化的假定，在计算过程中，完全杜绝舍入误差是不可能的，只能控制和约束由有限位算术运算带来的舍入误差的增长和危害，这样直接法得到的解也不一定是绝对精确的。迭代法是将方程组的解看作某种极限过程的向量极限的值，像第2章中非线性方程求解一样，计算极限过程是用迭代过程完成的，只不过将迭代式中单变量换成向量而已。在用迭代算法时，我们不可能将极限过程算到底，只能将迭代进行有限多次，得到满足一定精度要求的方程组的近似解。在数值计算历史上，直接解法和

34、迭代解法交替生辉。一种解法的兴旺与计算机的硬件环境和问题规模是密切相关的。一般说来，对同等规模的线性方程组，直接法对计算机的要求高于迭代法。对于中等规模的线性方程组，由于直接法的准确性和可靠性高，一般都用直接法求解。对于高阶方程组和稀疏方程组（非零元素较少），一般用迭代法求解。3.1 消元法3.1.1 三角形方程组的解形如下面三种形式的线性方程组较容易求解。对角形方程组 （3.3）设，对每一个方程，。显然，求解n阶对角方程的运算量为 。下三角方程组 （3.4）按照方程组的顺序，从第一个方程至第个方程，逐个解出。由方程，得。将的值代入到第二个方程得将的值代入到第个方程得计算需要次乘法或除法运算，

35、。因此，求解过程中的运算量为上三角方程组 （3.5）与计算下三角方程组的次序相反，从第个方程至第一个方程，逐个解出。由第个方程。将的值代入到第个方程得将的值代入到第个方程得解的通式 计算需要次乘法或除法运算。因此求解过程中的运算量为消元法的基本思想就是通过对方程组做初等变换，把一般形式的方程组化为等价的具有上述形式的易解方程组。3.1.2 高斯消元法与列主元消元法高斯消元法高斯消元法是我们熟悉的古老、简单而有效的解方程组的方法。下面是中学阶段解二元方程组（高斯消元法）的步骤：（3.6）（3.7）方程（3.6）乘以3加到第（3.7）个方程中得代入（3.6）得。其方法相当于对方程组的增广矩阵做行的

36、初等变换：已是上三角矩阵，而为原方程组的等价方程组，已化成易解的方程组形式。再用回代方法求解，得到：这就是高斯消元法解方程组的消元和回代过程。一般地，可对线性方程组（3.1）施行以下一系列变换；（1）对换某两个方程的次序；（2）对其中某个方程的两边同乘一个不为零的数；（3）把某一个方程两边同乘一个常数后加到另一个方程的两边。记变换后的方程组为： （3.8）显然方程组（3.1）与（3.8）是等价方程组，或者说它们有相同的解。分别记方程组（3.1）与（3.8）的增广矩阵为：可以看出，实际上是由按一系列初等换后得到的（1）对换某两行元素；（2）中的某行乘一个不为零的数；（3）把的某一行乘一个常数后加

37、到另一行。高斯消元法就是通过以上（3）的变换，把化为等价的上三角形式。下面我们以为例演示消元过程。设方程组： （3.9）其增广矩阵为：（1）若，则将第一行乘以加到第二行上；将第一乘以加到第三行上；将第一行乘以 加到第四行上得到 （3.10）即其中：（2）若则将第二行乘以加到第三行上；将第二行乘以加到第四行上，得到（3.11）其中：（3）若则将第三行乘以加到第四行上，得到 （3.12）其中：已是上三角矩阵，于是得到了与原方程等价的易解形式的方程组：（3.13）再对方程组（3.13）依次回代解出。从式（3.12）可以得到系数矩阵的行列式的值为的对角元素的乘积。即这也正是在计算机上计算方阵的行列式的

38、常规方法。要将上述解方程步骤推广到阶方程组，只需将对控制量“4”的操作改成对控制量的操作即可。元方程组高斯消元法的步骤如下：对于约定有（3.14）经过以上消元，我们已得到与等价的方程组，其中已是一个上三角矩阵。为简单起见，仍记的元素为 （3.15）即已得到原方程组的解。高斯消元法算法在算法中略去了变量，矩阵和向量的定义部分。在消元过程中，将仍放在元素的位置上。1输入：方程组阶数n，方程组系数矩阵A和常数向量b。2FOR k:=1 TO n-1 /消元过程 FOR i:=k+1 TO n / 假定 FOR j:=k+1 TO n / ENDFORJ / ENDFORI / ENDFORK3FOR

39、 i:=n TO1/回代求解 s:=bi FOR j:=i+1 TO n DO 4输出方程组的解。高斯消元法的运算量整个消元过程即式（3.14）的乘法和除法的运算量为回代过程即式（3.15）的乘除运算量为故高斯消元法的运算量为 （3.16）高斯消元法的可行性在上面的消元法中，未知量是按照在方程组中的自然顺序消去的，也叫顺序消元法。在消元过程中假定对然元素，消元步骤才能顺利进行，由于顺序消元不改变的主子式值，故高斯消元法可行的充分必要条件为的各阶主子式不为零。但是，实际上只要方程组就有解。故高斯消元法本身具有局限性。另一方面，即使高斯消元法可行，如果很小，在运算中用它作为除法的分母，会导致其它元

40、素数量级的严重增长和舍入误差的扩散。这是高斯消元法的另一缺陷。例3.1 方程组（3.17）（3.18）的精确解为：。在高斯消元法计算中取5位有效数字。解：方程（3.17）（1）/0.0003+方程（3.18）得：，代入方程（3.17）得。由此得到的解完全失真，如果交换两个方程的顺序，得到等价方程组经高斯消元后有由此可看到，在有些情况下，调换方程组的次序对方程组的解是有影响的，对消元法中抑制舍入误差的增长是有效的。如果不调换方程组的次序，取6位有效数字计算方程组的解，得到取9位有效数字计算方程组的解，得到由此可见有效数字在数值计算中的作用。列主元消元法如果在一列中选取按模最大的元素，将其调到主干

41、方程位置再做消元，则称为列主元消元法。调换方程组的次序是为了使运算中做分母量的绝对值尽量地大，减少舍入误差的影响。用列主元方法可以克服高斯消元法的额外限制，只要方程组有解，列主元消元法就能畅通无阻地顺利求解，同时又提高了解的精确度。更具体地，第一步在第一列元素中选出绝对值最大的元素，交换第一行和第行的所有元素，再做化简为零的操作。对于每个在做消元前，选出中绝对值最大的元素，对行和行交换后，再做消元操作，这就是列主元消元法的操作步骤。由于，可证中至少有一个元素不为零，从理论上保证了列主元消元法的可行性。列主元消元法与高斯消元法相比，只增加了选列主元和交换两个方程组（即两行元素）的过程。列主元消元

42、法算法1输入：方程组阶数，方程组系数矩阵和常数向量项。2 /选主元的消元过程 /选择 / 交换第行和第行 / ENDFORI / ENDFORK3FOR i:=n TO 1 / 回代求解4输出方程组的解 。如果对于第步，从行至行和从列至列中选取按模最大的，对第行和第行交换，对第列和第v列交换，这就是全主元消元法。在k列和第v列交换时，还要记录下v的序号，以便恢复未知量xk和xv的位置。3.1.3 高斯若尔当（GaussJordan）消元法高斯消元法将系数矩阵化为上三角矩阵，再进行回代求解；高斯若尔当消元法是将系数矩阵化为对角矩阵，再进行求解，即对高斯消元法或列主元消元法得到的等价增广矩阵：用第

43、4行乘加到第3行上，用第4行乘加到第2行上，用第4行乘加到第1行上，得到用第3行乘加到第2行上，用第3行乘加到第1行上，再用第2行乘加到第1行上，得到（3.19）为方便起见，仍记式（3.19）系数矩阵元素为 ，常数项元素为。那么用初等变换化系数矩阵为对角矩阵的方法称为高斯若尔当消元法。从形式上看对角矩阵（3.15）比上三角矩阵（3.11）更为简单，易于计算 ，但是将上三角矩阵（3.11）化为对角矩阵（3.15 ）的工作量略大于上三角矩阵回代的工作量。高斯若尔当消元法求逆矩阵设为非奇异矩阵，方程组的增广矩阵为。如果对应用高斯-若尔当消元法化为，则。例3.2 用高斯-若尔当消元法求的逆矩阵。解：解

44、得：3.2 直接三角分解法仍以为例，在高斯消元法中，从对方程组进行初等变换的角度观察方程组系数矩阵的演变过程。第一次消元步骤将方程组（3.9）化为方程组（3.10），相当于用了三个初等矩阵左乘 和。记，容易验证，由得到其中将方程组（3.10）化为方程组（3.11），相当于用了两个初等矩阵左乘 和。记有同理，将方程组（3.11）化为方程组（3.12），相当于用一个初等矩阵左乘 和。记，有完成了消元过程，有亦有所有消元步骤表示为：左乘一系列下三角初等矩阵。容易验证，这些下三角矩阵的乘积仍为下三角矩阵，并有于是有或这里仍为下三角矩阵，其对角元素为1，称为单位下三角矩阵，而已是上三角矩阵。记，则有结果表明若消元过程可行，可以将分解为单位下三角矩阵与上三角矩阵的乘积。由此派生出解方程组的直接分解法。由高斯消元法得到启发，对消元的过程相当于将分解为一个上三角矩阵和一个下三角矩阵的过程。如果直接分解得到和，。这时方程化为，令，由解出；再由，解出。这就是直接分解法。将方阵分解为，当是单位下三角矩阵，是上三角矩