数值计算教程

数值计算方法的对象与特点什么是数值计算方法现代的科学技术发展十分迅速，他们有一个共同的特点，就是都有大量的数据问题。比如，发射一颗探测宇宙奥秘的卫星，从卫星设计开始到发射、回收为止，科学家和工程技术人员、工人就要对卫星的总体、部件进行全面的设计和生产，要对选用的火箭进行设计和生产，这里面就有许许多多的数据要进行准确的计算。发射和回收的时候，又有关于发射角度、轨道、遥控、回收下落角度等等需要进行精确的计算。有如，在高能加速器里进行高能物理试验，研究具有很高能量的基本粒子的性质、它们之间的相互作用和转化规律，这里面也有大量的数据计算问题。计算问题可以数是现代社会各个领域普遍存在的共同问题，工业、农业、交通运输、医疗卫生、文化教育等等，各行各业都有许多数据需要计算，通过数据分析，以便掌握事物发展的规律。研究计算问题的解决方法和有关数学理论问题的一门学科就叫做计算方法。计算方法属于应用数学的范畴，它主要研究有关的数学和逻辑问题怎样由计算机加以有效解决。数值计算方法的内容数值计算方法也叫做计算数学或数值分析。数值计算方法主要内容包括非线性方程求根、线性代数方程组解法、微分方程的数值解法、插值问题、函数的数值逼近问题、概率统计计算问题等等，还要研究解的存在性、惟一性、收敛性和误差分析等理论问题。我们知道五次及五次以上的代数方程不存在求根公式，因此，要求出五次以上的高次代数方程的解，一般只能求它的近似解，求近似解的方法就是数值分析的方法。对于一般的超越方程，如对数方程、三角方程等等也只能采用数值分析的办法。怎样找出比较简洁、误差比较小、花费时间比较少的计算方法是数值分析的主要课题。在求解方程的办法中，常用的办法之一是迭代法，也叫做逐次逼近法。迭代法的计算是比较简单的，是比较容易进行的。迭代法还可以用来求解线性方程组的解。求方程组的近似解也要选择适当的迭代公式,使得收敛速度快，近似误差小。在线性代数方程组的解法中，常用的有塞德尔迭代法、共钝斜量法、超松弛迭代法等等。此外，一些比较古老的普通消去法，如高斯法、追赶法等等，在利用计算机的条件下也可以得到广泛的应用。在数值计算方法中，数值逼近也是常用的基本方法。数值逼近也叫近似代替，就是用简单的函数去代替比较复杂的函数，或者代替不能用解析表达式表示的函数。数值逼近的基本方法是插值法。初等数学里的三角函数表，对数表中的修正值，就是根据插值法制成的。在遇到求微分和积分的时候，如何利用简单的函数去近似代替所给的函数，以便容易求微分和求积分,也是计算方法的个主要内容。微分方程的数值解法也是近似解法。常微分方程的数值解法有欧拉法、预测校正法等。偏微分方程的初值问题或边值问题，目前常用的是有限差分法、有限元素法等。有限差分法的基本思想是用离散的、只含有限个未知数的差分方程去代替连续变量的微分方程和定解条件。求出差分方程的解法作为求偏微分方程的近似解。有限元素法是近代才发展起来的，它是以变分原理和剖分差值作为基础的方法。在解决椭圆形方程边值问题上得到了广泛的应用。目前，有许多人正在研究用有限元素法来解双曲形和抛物形的方程。数值计算方法的内容十分丰富，它在科学技术中正发挥着越来越大的作用。数值计算方法的特点第一，面向计算机，要根据计算机特点提供实际可行的有效算法。即算法只能包括加、减、乘、除和逻辑运算，是计算机能直接处理的。第二，有可靠的理论分析，能任意逼近并达到精度要求，对近似算法要保证收敛性和数值稳定性，还要对误差进行分析。第三，要有好的计算复杂性，时间复杂性好是指节省时间，空间复杂性好是指节省存储量，这也是建立算法要研究的问题，它关系到算法能否在计算机上实现。第四，要有数值实验，即任何一个算法除了从理论上要满足上述三点外，还要通过数值实验证明是行之有效的。根据“数值计算方法''的特点，学习时我们要注意掌握方法的基本原理和思想，要注意方法处理的技巧及其与计算机结合，要重视误差分析、收敛性及稳定性的基本理论：其次，要通过例子，学习使用各种数值方法解决实际计算问题；最后，为了掌握本课的内容，还应做一定数量的理论分析与计算问题练习。由于本课内容包括了微积分、代数、常微分方程的数值方法，以及高级程序设计的方法，读者必须掌握这几门课的基本内容才能学好这一课程。误差来源与误差分析误差的来源早在中学我们就接触过误差的概念，如在做热力学实验中，从温度计上读出的温度是23.4C就不是一个精确的值，而是含有误差的近似值。事实上，误差的在我们的生活中无处不在，无处不有。如量体裁衣,量与裁的结果都不是精确无误的，都有误差。人们可能会问：如果使用计算机来解决问题，结果还会有误差吗？下面我们通过考察数学方法解决实际问题的主要过程来思考这个问题。用数学方法解决一个具体的实际问题，首先要建立数学模型，这就要对实际问题进行抽象、简化，因而数学模型本身总含有误差，这种误差叫做模型误差。在数学模型中通常包含各种各样的参变量，如温度、长度、电压等，这些参数往往都是通过观测得到的，因此也带来了误差，这种误差叫做观测误差。当数学模型不能得到精确解时，通常要建立一套行之有效的数值方法求它的近似解，近似解与准确解之间的误差就称为截断误差或方法误差.由于在计算机中浮点数只能表示实数的近似值，因此用计算机进行实际计算时每一步都可能有误差，这种误差称为舍入误差。

例如，函数#幻用泰勒(Taylor)展开式"必平"家一铲近似代替，则数值方法的截断误差是(:介于o(:介于o与，之间)又如，在计算时用3.14159近似替代：产生的误差*■<-3.14159=0.0000026 就是舍入误差。上述种种意味着都会影响计算结果的准确性，因此需要了解与研究误差。在数值计算方法中将着重研究截断误差、舍入误差，并对它们的传播与积累作出分析。绝对误差、相对误差与有效数字.绝对误差与绝对误差限定义1.1设乂为准确值，：为x的一个近似值，称/・/一玄为近似值的绝对误差，或误差。通常我们无法知道准确值2,也不能算出误差的准确值一，只能根据测量或计算估计出误差的绝对值不超过某个正数鼠，则/称为绝对误差限。有了绝对误差限就可知道x的范围,--4*4即△落在区间［*•一『•内。例如用毫米测度尺测量一长度？，读出的长度为23mm,则有^一^“06mm。由此例也可以看到绝对误差是有量纲和单位的。.相对误差与相对误差限只用绝对误差还不能说明数的近似程度，例如甲打字时每百个字错一个，乙打字时平均每千个字错一个，他们的误差都是错一个，但显然乙要准确些。这就启发我们除了要看绝对误差大小外，还必须顾及量的本身。定义1.2把近似值的误差£与准确值工的比值称为近似值T的相对误差，记作实际计算时，由于真值工通常是不知道的，通常取'7。相对误差也可正可负，它的绝对值的上■•■8 1* '|t*| 4—,1%界叫做相对误差限。记作即II。根据定义，甲打字时的相对误差 100 ,乙打字时的r；5-L-01%相对误差 wooo易知相对误差是个无量纲量。.有效数字当下有多位时，常常接四舍五入的原则得到*的前几位近似值例如x-x-314159265*-取3位k；・3.®440002；取5位它们的误差都不超过末位数字的半个单位，即卜一生H|£Lx10、,卜-3J4iq4lxVT*2 2现在我们将四舍五入抽象成数学语言，并引入一个新名词“有效数字”来描述它。定义1.3若近似数T的误差限是某一位的半个单位，该位到丁的第一位非零数字共有q位，我们就说，•有七位有效数字。如取x•“殳14作丁的近似值，工•就有3位有效数字；取x'=3.M]6作：的近似值，丁就有5位有效数字。有效数字也可采用以下定义:X的近似数X*写成标准形式:其中,是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位有效数字的近似值是&0000而不是8,8只有一位有效数字。.有效数字与误差关系不难看出，若由（1.1）给出某近似数有V位有效数字，则可以从（1.2）求得这个数的绝对误差限2因此，在色相同的情况下，':越大则IL就越小，故有效数字越多，绝对误差限越小。定理1.1用（1.1）表示的近似数/，若工”具有,：位有效数字，则其相对误差限为阴. —■—*IOT*** .反之，若f的相对误差限 则T■至少具有工位有效数字。证明由式（1.1）得所当？‘有暗位有效数字时所当？‘有暗位有效数字时IH卜】户12•冒,球产・如。3<1--—心2，反之，若一的相对误差限 2gl+D有：卜"卜kka><(Wl卜"卜kka><(Wl>x1oWX2^由式（1.2）式知道，寸至少有三位有效数字，证毕这说明有效数字越多，相对误差限越小。例1.2要使亚的近似值的相对误差限小于om,要取儿个有效数字？解：因为屈.4.47…，所以q7,令—XI0*^01%

2"故取K-4即可满足。5.函数计算的误差估计一般情况，当自变量有误差时计算相应的函数值也会产生误差，其误差限可利用函数的泰勒展开式估计。设“X）是一元函数，工的近似值为味以/（,）近似,㈤，其误差界记作式/6»可用泰勒展开式/W 力-八戏”-力+婴Qr尸:介于了，丁之间，取绝对值得假定义电与/*（吟的比值不太大，可忽略贰”■）的高阶项，于是可得计算函数值的误差限加**力W4解|aT%#a.>例1.3设/的相对误差限为2%,则/"）.（*•『的相对误差限是多少？解：阳m”户浮Cx2M・“g»所以/（0的相对误差限为0.02no当J为多元函数时，如计算♦■/（%卬…，4）。如果、4…'*"的近似值为年*…H,则看的近似值尸心，于是函数值三*的误差心"）山泰勒展开式得/瓜,冷….4)-/(。卬….Q于是误差限为如疆*力例1.4已测得某场地长’的值为宽丁的值为已知卜小［«。**»,试求面积的绝对误差限马相对误差限。s-4/.— -/解：因也㈤，那么3牌卜）俘I"）其中令，■，-80w,琮），・r-110*1于是绝对误差限«（『）-»*02*110x01・2X6相对误差限1.3误差传播与若干防治办法由前述可知，在数值计算中每步都可能产生误差。而一个问题的解决，往往要经过成千上万次运算，我们不可能每步都有加以分析，下面，通过对误差的某些传播规律的简单分析，指出在数值计算中应注意的几个原则，它有助于鉴别计算结果的可靠性并防止误差危害现象的产生。.要避免两相近数相减在数值运算中两相近数相减会使有效数字严而损失，例如*'5父85.尸・53252都具有五位有效数字，但只有两位有效数字，所以要尽量避免这类运算。1g%-1g吗-1g—通常采用的方法是改变计算公式，例如当天与勺很接近时，由于 巧那么可用右端的公式代替左端的公式计算，有效数字就不会损失。当工很大时， 必也可用右端来代替左端。一般情况，当时，可用泰勒展开，）+，）+q：）（x-/+…取右端的有限项近似左端。如果计算公式不能改变，则可采用增加有效位数的方法。.要防止大数“吃掉”小数若参加运算的数的数量级相差很大，而计算机的位数有限，如不注意运算次序，就可能出现大数“吃掉“小数的现象，影响计算结果。例如在五位十进制计算机上，计算写成规格化形式MM $A-0.52492»10,+^0.000001x10由于计算时要对阶，在计算机中表示为o,计算出来A・0524S=1Q‘，结果严重失真!m小如果计算时，先将H计算出来，再与52492相加，就不会出现大数“吃掉”小数的现象了。.注意简化计算步骤，减少运算次数同样一个计算题，如果能减少运算次数，不但可节省计算机的计算时间，还能减少舍入误差。例如，计算x255的值，如果逐个相乘要用254次乘法，但若写成产-X//只要做14次乘法运算即可。.绝对值太小的数不宜作除数•X •X•十X*—设T与；分别有近似值 ，的近似值，，其绝对误差F®卜卜廿声，M

i/r1»严।显然,当g很小时，近似值工”的绝对误差,仁）有可能很大。因此，不宜把绝对值太小的数作除数。认.，■爷? ■631%KrI。W第2章非线性方程求根与线性方程相比，非线性方程问题无论是从理论上还是从计算公式上，都要复杂得多。对于一般的非线性方程/（力・0,计算方程的根既无一定章程可寻也无直接法可言。例如，求解高次方程组7/♦x-L5=0的根，求解含有指数和正弦函数的超越方程/■幻的零点。解非线性方程或非线性方程组也是计算方法中的•个主题。一般地，我们用符号/（*）来表示方程左端的函数，方程--般形式表示为〃工）・°,方程的解称为方程的根或函数的零点。通常，非线性方程的根不止一个，对于非线性方程，一般用迭代法求解。因此，在求解非线性方程时,要给定初始值或求解范围。实根的对分法使用对分法的条件对分法或称二分法是求方程近似解的一种简单直观的方法。设函数人力在〔。力1上连续，且/㈤则在［6网上至少有-零点，这是微积分中的介值定理，也是使用对分法的前题条件。计算中通过对分区间、缩小区间范围的步骤搜索零点的位置。例2.1用对分法求解,（x）*37.7/+19.2x-153在区间［词之间的根.解：/G）--2-8/o-Q3由介值定理可得有根区间1层四・［闭。（2）计算均-。+期2-1.“。9・《45有根区间［»］叩.切（3）计算三=（16+2）/2=1.75,/Q4=0.078125,有根区间—。工工了可。（4）•直做到y（/）K.（计算前给定的精度）或I。一坪£时停止。2.1.2对分法求根算法计算/（动■。的一般计算步骤如下：.输入求根区间SQ1和误差控制量凡定义函数〃力。IF THEN退出选用其他求求根方法.WH/LE・时（1）计算中点*.g*劫,2以及/（目的值；（2）分情况处理l/WK«:停止计算，・x,转向步骤4/（•VW（Q修正区间141卜>（«•*!/W/W<0修正区间［3卜》［@*1ENDWHILE.输出近似根在算法中，常用血贰/3）吨底/@>°代替八视/⑸>°的判断，以避免/1（«）・/（力数值溢出。对分法的算法简单，然而，若/（*）在S/1是有几个零点时，只能算出其中一个零点；另一方面，即使〃力在上有零点，也未必有/触怨）〈匕这就限制了对分法的使用范围。对分法只能计算方程的实根。2.2迭代法对给定的方程将它转达换成等价形式：给定初值外，由此来构造迭代序列尸，如果迭代法收敛，即出3,此武与）”，尸，如果迭代法收敛，即出3,此武与）”，有a・a°）,则现就是方程(2)--5(2)--52r・P-5,迭代格式 2人才■。的根。在计算中当I&*~匕1小于给定的精度控制量时，取@为方程的根。例如，代数方程？-2x-5・0的三种等价形式及其迭代格式如下:(1)/■&f+5.x■胱〃，5.迭代格式%.4匕*+5(1)(3)迭代格式(3)对于方程构造的多种迭代格式怎样判断构造的迭代格式是否收敛？收敛是否与迭代的初值有关？定理2.1 若定义在［。为|上，如果满足（D当有2］有《配（2）KG在上可导，并且存在正数"1,使对任意的问6村，有IrlWl"；则在［。.句上有惟一的点片满足,"»<**>,称寸为K©的不动点。而且迭代格式“_对任意的初值均收敛于♦<*）的不动点黑■,并有误差估计式证明：（1）先证明存在性： ）・x-a*）,则有/何…故有a 使得（2）再证明惟一性：设都是•<*）的不动点，且父**・，则有与假设矛盾，这表明*~xn,即不动点是惟一的。(3)当*0—同时，由于«x)w[a向可用归纳法证明，迭代序列QJu(«网，于是由微分中值定理加和得"|--，卜4川|勺-，| (2.2)因为4<1,所以当北~>9时，""70%"*’即迭代格式加・^^收敛。(4)误差估计式：设A固定，对于任意的正整数3有，1%,-qH7V-9.1Hl+1+-£(*4+tf*Z 卜注：定理2.1是判断迭代法收敛的充分条件，而非必要条件。要构造满足定理条件的等价形式一般是难于做到的。事实上，如果寸为JQ)的零点，若能构造等价形式“女工而"2KI,由。口)的连续性，一定存在£的邻域仁・艮八「)，其上有I"（才KL〈i,这时若初值迭代也就收敛了。由此构造收敛迭代式有两个要素，其一，等价形式*■•*）应满足i.s'ki；其二，初值必须取自寸的充分小邻域，这个邻域大小决定于函数/a）,及做出的等价形式“■«©。例2.2求代数方程？-2x-5-0,在4・2附近的实根。解：1）2・2/5―广取%+5、2 1 T’（2x+/.（杯】，当x€[l.*25]，构造的迭代序列收敛。取勺.？则=2.08008, 4=2.09235, =2.094217"*=2.094494,马=2.094543,三=2.094550准确的解是工=2.094551481502）将迭代格式写为亡-5,%a?-5

54--^-0任）・一^-力。；。）卜容卜产口」2S]’•迭代格式.仍（Q不能保证收敛。对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢,从而使计算量变得很大，因此迭代过程的加速是个重要的课题。以下介绍一种埃特金(Aitken)方法。对X.山)方程，构造加速过程，算法如下：(1)预测：(2)校正：(3)改进：有些不收敛的迭代法，经过埃特金方法处理后，变得收敛了。例2.3求方程J® 在勺附近的根一。解：若采用迭代公式："=匕7迭代法是发散的，我们现在以这种迭代公式为基础形成埃特金算法:取&■15计算结果如表2.1所示:表2.1计算结果t%举♦】01.512.375001339651.4162921.840925.238881.3556531.491402.317281.3289541.347101.444351.3248051.325181.327141.32472我们看到，将发散的迭代公式通过埃特金方法处理后，竞获得了相当好的收敛性。2.4牛顿迭代法对方程」(*)‘°可构造多种迭代格式，牛顿迭代法是借助于对函数,Q)■°作泰勒展开而构造的一种迭代格式。将八外・。在初始值工作泰勒展开：/w- a*弯.--专尸*••4■取展开式的线性部分作为的近似值，则有设fa/o,则勺/W令飞■/6)类似地，再将/(*)-0在均作泰勒展开并取其线性部分得到：一宜做下去得到牛顿迭代格式：■・4・ 上・L2…/⑷ , (2.3)牛顿迭代格式对应于/©)・°的等价方程是

(2.4)cr«)a(2.4)若M是/㈤的单根时，/8）・0./9）"°则仃*4°,只要初值三充分接近工，所以牛顿迭代收敛。当二为的五重根时，取下面迭代格式：比72“牛顿法的几何意义以小（勺）为斜率作过点的切线，即作了（外在士点的切线方程令*-Q,则得此切线与工轴的交点电，即再作在《》/（砧）处的切线，得交点与，逐步逼近方程的二根.如图2.2所示图2.2牛顿切线法示意图在区域［勺，4 的局部“以直代曲'’是处理非线性问题的常用手法。在泰勒展开中，截取函数展开的线性部分替代/㈤。例2.4用牛顿迭代法求方程/（x）-/-7.7/+19一2X-I5.3,在、附近的根。-7.7^*19.2^-15.3解：XA―乂754—192计算结果列于表2.2中。表2.2计算结果Kxkf(x)01.00-2.811.41176-0.72707121.62424-0.14549331.6923-0.013168241.69991-0.000151551.70比较表2.1和表2.2的数值，可以看到牛顿迭代法的收敛速度明显快于对分法。牛顿迭代法也有局限性。在牛顿迭代法中，选取适当迭代初始值三是求解的前题，当迭代的初始值壬在某根的附近时迭代才能收敛到这个根，有时会发生从一个根附近跳向另一个根附近的情况，尤其在导数数值很小时，如图2.3所示。图2.3失效的牛顿迭代法牛顿迭代算法.定义函数/(外，«(«),输入或定义迭代初始值和控制精度epsilon。.FORSOtoMAXREPTx_k\-g(x_JW) "g(才是牛慎迭代公式IF(\x_*1-x_JtO|<apsdoM)oRa/(x_JtQK^oK)THEN｛输出满足给定精度的近似解结束｝ENDFOR.输出：在J*0附近无根。注：在伪码中，"/r'表示注释语句。伪码又称过程设计语言，主要用于描述算法。它是某种高级语言和自然语言的混杂语言，它取某种高级语言中的一些关键字，用于描述算法的结构化构造和数据说明等。伪码的语句中嵌有自然语言的叙述，伪码易于2.5弦截法弦截法迭代格式j-X.2--TOC\o"1-5"\h\z在牛顿迭代格式中： 『g\o"CurrentDocument",一㈤・"32 〃、用差商 q-Ri代替导数并给定两个初始值七和工】，那么迭代格式可写成如下形式：〃必-〃2 （2.5）上式称为弦截法。用弦截法迭代求根，每次只需计算一次函数值，而用牛顿迭代法每次要计算一次函数值和一次导数值。但弦截法收敛速度稍慢于牛顿迭代法。弦截法的几何意义做过两点（、•/6》和的一条直线（弦），该直线与％轴的交点就是生成的迭代点为，再做过（%/（♦»和（马的一条直线，力是该直线与工轴的交点，继续做下去得到方程的根/依）-0,如图2.4所示图2.4弦截法示意图例2.5用弦截法求方程〃X）■/-7.7?*19.2x-15.3的根，取勺.1.5.-4.0。解：“‘

计算结果列于表2.3中。表2.3计算结果上二/W01.5-0.45142.321.909090.24883531.65543-0.080569241.717480.028745651.701160.0019590261.69997-0.000053924671.79.45910-8弦截法算法.定义函数输入或定义控制精度值epsilon和迭代初始值计算/匕・/。_初 表示％.勺.FOR:=0,1…MAXREPT/21・/«_北为2->x_A:-x_Jt2-/2(x_A2〃/_上表示4MTHEN｛输出满足给定精度的近似解结束｝2.4x_kl-x_k2x_k2:~x_k〃为下一次迭代准备数值ENDFOR.输出：在初始值、・劫*・上2附近无根。2.6非线性方程组的牛顿方法为了叙述简单，我们以解二阶非线性方程组为例演示解题的方法和步步骤，类似地可以得到解更高阶非线性方程组的方法和步骤。设二阶方程组以5・。 (2.6)其中X•，为自变量。为了方便起见，将方程组写成向量形式：“您做fl旧仁川，其中二»J将附近作二元泰勒展开，并取其线性部分，得到方程组4(勺㈤+(r)当型♦。-%)丝-oTOC\o"1-5"\h\zox dy〃勺㈤WR好必”卬逐衿・0(M .令X-(.&.y一“则有dx a&当也唱出rg川w qy

弧也生生如果dx分，解出&A,得再列出方程组：晔小巡守32•_〃”)n dy‘^2q --工的比M 夕解出,&7r 修・k田得出继续做下去，每一次迭代都是解一个方程组…脚£图力(2.7)&卬Ar.九*-%即“・".&①”・为.2,直到为止。例2.6求解非线性方程组4GM・4---y-0/式八力-y-0取初始值

画渔，〃、班打 -2y\解:玄笠-lJ5Sy解:"43(；：sa：：)(-2&+32--OLll-27l328Ar-^>-001828产.00042%]a(qJI-0.029849J解方程得'■国曾…斗尸力5J l-,7J(-0-029849J(-1.729849J继续做下去，直到mW&U》D<l°"时停止。将两个变量的非线性方程组推广到多变量的非线性方程组：人风,盯・…。-0((X*,-0 (28)记*…的向量形式:F(X)-0

的向量形式:F(X)-0用牛顿迭代法求方程组的解为:就⑺弥X) %⑺% I %跳⑶里(箱...秘⑶

血 与 改乜⑺孤就⑺弥X) %⑺% I %跳⑶里(箱...秘⑶

血 与 改乜⑺孤X)..^(X)"X）・其中：称J（x）为雅可比（jacobi）矩阵。计算中采用4犷乂/川-义与・化专NjMZ'-一或评 （29）一直做到小于给定精度为止。在上的领域中若”《也"AT,则迭代收敛。.7程序示例程序2.1用牛顿迭代法求/S）在工:附近的根。算法描述给定/（*）,从三开始，根据牛顿迭代公式计算」(x)在L附近的根。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII// Purpose:牛顿迭代法求根 〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h#include<math.h>//iterate。)为"#defineiterate(x)(x-((x*x-3)*x-1)/(3*x*x-3))#define/(力 ((x*x-3)*x-l)#definexO1.5 〃迭代初值xO#defineMAXREPT1000//最大迭代次数#defineepsilon0.0001 //精度voidmain(){inti;doublex_k=x0,x_kl=x0for(i=i<MAXREPT；i++)(printfCtGot...%f7n,\x_kl);x_k1=iterate(x_k); 〃迭代iffabs(x_kl-x_k)<epsilon||fabs(fi[x_kl)<epsilonprintfC'!Root:%f\n”,x_kl);〃满足精度，输出return;{x_k=x_kI{printfi(tfcAfter%drepeate,nosolvedAm'\MAXREPT)；// EndofFile 计算实例已知取4 用牛顿迭代法计算/㈤的根。程序输入输出由本程序预设的及小.L5,得到输出!Root:1.879385程序2.2用弦截法求/(“)在内通附近的根。算法描述给定/Q),从、•再开始，根据弦截法迭代公式求得/(*)在其附近的根。程序源码〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃///Purpose：弦截法求根〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<math.h>#definef(x)(x*x*x-7.7*x*x+19.2*x-15.3)//f函数#definexO0.0〃初值xO,x1#definexl1.0#defineMAXREPT1000 //最大迭代次数#defineepsilon0.00001〃求解精度voidmain(){inti;doublex_k=x0,x_k1=x1,x_k2=x1;fbr(I=O;i<MAXREPT;i++={printf("!Root:%f\n”,x_k2);x_k2=x_kl-(fi[x_kl)*x_k1・x_k))/(f(x_k1)・Rx_k));〃弦截求新x_niffabs(x_k2-x_k1)<epsilon||fabs(f(s_k2))<epsilon{printf("!Root:%f\n,,,x_k2);〃满足精度,输出return;(x_k=x_kl;x_kl=x_k2; 〃反复{printfi(t4After%drepdate,nosolvde.'n",MAXREPT);// EndofFile 计算实例用弦截法求方程〃x)・J-7.7/+19.2x-15.3的根，取.-I'q程序输入输出由本程序的,Q)及4•'得到输出!Root:1.700000第3章解线性方程组的直接法在近代数学数值计算和工程应用中，求解线性方程组是雨:要的课题。例如，样条插值中形成的“关系式，曲线拟合形成的法方程等，都落实到解一个v元线性方程组，尤其是大型方程组的求解，即求线性方程组(3.1)的未知量的数值。■再*aw+—>a.x.-4%4*31*-***-A与*""，-4 (31)其中Gjbi为常数。上式可写成矩阵形式4x=b,即其中，AF哂为系数矩阵("吃"，产为解向量…为常数向量。当deS=£)0时，由线性代数中的克莱姆法则，方程组及的解存在且惟•,且有q为系数矩阵三的第7列元素以◎代替的矩阵的行列式的值。克莱姆法则在建立线性方程组解的理论基础中功不可没，但是在实际计算中，我们难以承受它的计算量。例如，解一个100阶的线性方程组，乘除法次数约为(101T00!-99),即使以每秒】©的运算速度，也需要近产年的时间。在石油勘探、天气预报等问题中常常出现成百上千阶的方程组，也就产生了各种形式方程组数值解法的需求。研究大型方程组的解是H前计算数学中的一个重要方向和课题。解方程组的方法可归纳为直接解法和迭代解法。从理论上来说，直接法经过有限次四则运算，假定每一步运算过程中没有舍入误差，那么，最后得到方程组的解就是精确解。但是，这只是理想化的假定，在计算过程中，完全杜绝舍入误差是不可能的，只能控制和约束由有限位算术运算带来的舍入误差的增长和危害，这样直接法得到的解也不一定是绝对精确的。迭代法是将方程组的解看作某种极限过程的向量极限的值，像第2章中非线性方程求解一样，计算极限过程是用迭代过程完成的，只不过将迭代式中单变量X**】换成向量而已。在用迭代算法时，我们不可能将极限过程算到底，只能将迭代进行有限多次，得到满足一定精度要求的方程组的近似解。在数值计算历史上，直接解法和迭代解法交替生辉。一种解法的兴旺与计算机的硬件环境和问题规模是密切相关的。一般说来，对同等规模的线性方程组，直接法对计算机的要求高于迭代法。对于中等规模的线性方程组伊^^，由于直接法的准确性和可靠性高，一般都用直接法求解。对于高阶方程组和稀疏方程组（非零元素较少），一般用迭代法求解。消元法三角形方程组的解形如下面三种形式的线性方程组较容易求解。对角形方程组TOC\o"1-5"\h\z%^勺 ■与■■' %勺-X （3.3）设对每一个方程，/,取,TN…』。显然，求解〃阶对角方程的运算量为下三角方程组，马 -4* ■巧 "与+ 4勺*…♦ ・ 4 （3.4）按照方程组的顺序，从第一个方程至第K个方程，逐个解出由方程得看将』的值代入到第二个方程♦'gK,■幻得 勺-眄-’11砧〃必将的值代入到第？个方程。凝+…一x「孰■IN”-m得 P计算不需要,次乘法或除法运算，i-t2-,■**.因此，求解过程中的运算量为g・罕•*)上三角方程组%-*x* … *F.4* - *%。_%<…433* -如=% (3.5)与计算下三角方程组的次序相反，从第七个方程至第一个方程，逐个解出2j.从"-1•…2・L由第t个方程%。・.待。将三的值代入到第4-1个方程将、•/士…Eu的值代入到第;个方程F0-4得解的通式—部W.，…V.A1计算号需要“T+1次乘法或除法运算1**・"-L…・2・L因此求解过程中的运算量为7♦I・£/・心；D・吁）1 /-I 2消元法的基本思想就是通过对方程组做初等变换，把•般形式的方程组化为等价的具有上述形式的易解方程组。3.1.2高斯消元法与列主元消元法高斯消元法高斯消元法是我们熟悉的古老、简单而有效的解方程组的方法。下面是中学阶段解二元方程组（高斯消元法）的步骤：fx,-X,-1 （3.6）13X（*2x）-8 （3.7）方程（3.6）乘以-3加到第（3.7）个方程中得03-57代入（3.6）得“IN其方法相当于对方程组的增广矩阵做行的初等变换:

（3。：：）3』（；；；卜史）W已是上三角矩阵，而为原方程组的等价方程组，已化成易解的方程组形式。再用回代方法求解，得到：-2这就是高斯消元法解方程组的消元和回代过程。一般地，可对线性方程组（3.1）施行以下一系列变换:（1）对换某两个方程的次序;（2）对其中某个方程的两边同乘•个不为零的数;（3）把某•个方程两边同乘•个常数后加到另一个方程的两边。记变换后的方程组为:才“%*-♦!。-E(38)显然方程组（3.1）与（3.8）是等价方程组,显然方程组（3.1）与（3.8）是等价方程组,的增广矩阵为:可以看出，J实际上是由《4与按-系列初等换后得到的(1)对换(/杨某两行元素；(2)(九力中的某行乘一个不为零的数；(3)把(4年的某一行乘一个常数后加到另一行。高斯消元法就是通过以上(3)的变换，把《49化为等价的上三角形式。下面我们以"・彳为例演示消元过程。设方程组：ZlA -瓦0M巧*0^X2+〜占+0Mx.-6,♦Mi*<V3+.5F玉f%。+4用+3,-4 (39)其增广矩阵为：(1)若,则将第一行乘以-%/%i加到第二行上；将第乘以-加到第三行上；将第一行乘以一加到第四行上得到a*-呼其中:a,..-egj.2%*f°叱-5%/,/-2,X4⑵若词“a则将第二行乘以一造,右加到第三行上；将第一行乘以-马‘啰加到第四行上,得到all

00(3.11)0(3.11)其中:谭■碎-啰-3.4

中-if9-岬噌3,4⑶铲ff"则将第三行乘以一点/蟾*加到第四行上，得到(3.12)其中喑•虐-噌喈础＞--已是上三角矩阵，于是得到了与原方程等价的易解形式的方程组：再对方程组（3.13）依次回代解出

从式(3.12)可以得到系数矩阵＜的行列式的值为彳的对角元素的乘积。即这也正是在计算机上计算方阵W的行列式的常规方法。要将上述解方程步骤推广到三阶方程组，只需将对控制量“4”的操作改成对控制量V的操作即可。t元方程组高斯消元法的步骤如下：对于士-L2”.〃-l约定《■。有6・铲，铲/邛)中6・铲，铲/邛)中**l<i<m经过以上消元，我们已得到与出等价的方程组后.A,其中之已是一个上三角矩阵。为简单起见，仍记2的元素为“宕的兀素力当月归即己得到原方程组的解。高斯消元法算法1.输入：方程组阶数〃，方程组系数矩阵4和常数向量入2.FORk:=lTOn-1〃消元过程2.FORk:=lTOn-1〃消元过程FORi:=k+lTOn[1■/aFORj:=k+lTOn

{■产}〃endforj//ENDFORI//ENDFORK3.FORi:=nTOI〃回代求解3.FORi:=nTOI〃回代求解{s:=biFORj:=i+lTOnDO4.输出方程组的解勺J-12…。高斯消元法的运算量整个消元过程即式（3.14）的乘法和除法的运算量为JU4 Al 4回代过程即式（3.15）的乘除运算量为网”.1）-2~故高斯消元法的运算量为(3.16)3+- ■g,）

(3.16)高斯消元法的可行性在上面的消元法中，未知量是按照在方程组中的自然顺序消去的，也叫顺序消元法。在消元过程中假定对然元素小 消元步骤才能顺利进行，由于顺序消元不改变£的主子式值,故高斯消元法可行的充分必要条件为W的各阶主子式不为零。但是，实际上只要修才~“方程组就有解。故高斯消元法本身具有局限性。另一方面，即使高斯消元法可行，如果簿很小，在运算中用它作为除法的分母会导致其它元素数量级的严重增长和舍入误差的扩散。这是高斯消元法的另一缺陷。(0.0003X.*10000x.-20001 (3J7)L000011*10OOOx,-1.0000 (3.18)的精确解为:2-3在高斯消元法计算中取5的精确解为:2-3在高斯消元法计算中取5位有效数字。解：方程(3.17)x(-1)/0.0003+方程(3.18)得:09003不♦39000*,-200019999.0»s-6666.0,-.Xj-0.M67代入方程(3.17)得“】•°。由此得到的解完全失真，如果交换两个方程的顺序,得到等价方程组I.OOOOAj*1.0000^-1.00W

10.0003jh+39000啊-20001经高斯消元后有(10000、+10000x2-100002劳的》・1诩8.-.Xj-06667,«-0.3333由此可看到，在有些情况下，调换方程组的次序对方程组的解是有影响的，对消元法中抑制舍入误差的增长是有效的。如果不调换方程组的次序，取6位有效数字计算方程组的解，得到取9位有效数字计算方程组的解，得到力・06666£7.X1-0333333由此可见有效数字在数值计算中的作用。列主元消元法如果在一列中选取按模最大的元素，将其调到主干方程位置再做消元，则称为列主元消元法。调换方程组的次序是为了使运算中做分母量的绝对值尽量地大，减少舍入误差的影响。用列主元方法可以克服高斯消元法的额外限制，只要方程组有解，列主元消元法就能畅通无阻地顺利求解，同时又提高了解的精确度。max{kD・kl更具体地，第一步在第一列元素中选出绝对值最大的元素,交换第一行和第噌行的所有元素，再做化简吗为零的操作。对于每个七狂12…-T在做消元前，选出闾…中绝对值最大的元素对既行和航行交换后，再做消元操作，这就是列主元消元法的操作步骤。由于修//0,可证<«射，蜀震…中至少有一个元素不为零，从理论上保证了列主元消元法的可行性。列主元消元法与高斯消元法相比，只增加了选列主元和交换两个方程组（即两行元素）的过程。列主元消元法算法.输入：方程组阶数用，方程组系数矩阵W和常数向量项£。.FORt=lTOd-1 〃选主元的消元过程{〃选择. 小K*lFOR«=k+!TOnIFk|〉6THBH{胸・《；6・卜JFORk^TOq //交换第七行和第z行9f；@a：・7IFORr=k*lTOn（匕・。/小FORj:=kHTOn{,・％-«*/}ffENDPORJW4} //ENDFORI} //ENDFORK3.FORi:=nTO1 //回代求解4.输出方程组的解如果对于第克步，从全行至七行和从袅列至y列中选取按模最大的wW5记月卜蚓1aH对第—行交换,对第i列交换,这就是全主元消元法。在《列和第V列交换时，还要记录下v的序号，以便恢复未知量M和XV的位置。3.1.3高斯一若尔当(Gauss—Jordan)消元法高斯消元法将系数矩阵化为上三角矩阵，再进行回代求解；高斯一若尔当消元法是将系数矩阵化为对角矩阵，再进行求解，即对高斯消元法或列主元消元法得到的等价增广矩阵：alt%flU *10图4第b?0。电审Ap)00 0a曾即用第4行乘右加到第3行上,用第4行乘如曾加到第2行上,用第4行乘布里加到第1行上，得到〜0°裁0aSe0守0。嫉。守用第3行乘号'卡加到第2行上，用第3行乘-媾加到第1行上，再用第2行乘一//姆加到第1行上，得到'.0 0 06fo0喈0Q中0。点。理0 0 。眼,用 (3,19)为方便起见，仍记式(3.19)系数矩阵元素为气,常数项元素为当。那么用初等变换化系数矩阵为对角矩阵的方法称为高斯一若尔当消元法。从形式上看对角矩阵（3.15）比上三角矩阵（3.11）更为简单，易于计算三但是将上三角矩阵（3.11）化为对角矩阵（3.15）的工作量略大于上三角矩阵回代的工作量。高斯一若尔当消元法求逆矩阵设a为非奇异矩阵，方程组/**/・的增广矩阵为,-（工〃）。如果对。应用高斯-若尔当消元法化为则小，123,/-245例3.2用高斯-若尔当消元法求P56J的逆矩阵TOC\o"1-5"\h\z1 0 -1/2fo -5/2 2->0 1 3/20 3/2 -I0 0 1/21 -172 03.2直接三角分解法

相当于用了三个初等矩阵左乘山和土。记仍以为例，在高斯消元法中，从对方程组进行初等变换的角度观察方程组系数矩阵的演变过程。相当于用了三个初等矩阵左乘山和土。记第一次消元步骤将方程组（3.9）化为方程组（3.10）,〃.. .2*4,容易验证,1000W10 1000W10 100 0由卒”■Q•得到■6吗其中将方程组（3.10）化为方程组（3.11）,相当于用了两个初等矩阵左乘士和E。记Q•咽，渡八切.有I0 I0 00,0 1000 0 100-3014 0 0 0、0 10 00 -Zn I 00 0 0 10 0同理，将方程组(3.11)化为方程组(0 0同理，将方程组(3.11)化为方程组(3.12),相当于用一个初等矩阵左乘完成了消元过程，有亦有所有消元步骤表示为：（4与左乘一系列下三角初等矩阵。容易验证，这些下三角矩阵的乘积工仍为下三角矩阵，并有「I000'*11 *12 ■ °/•IGiLq彳彳/■a于是有X-7X或这里『‘仍为下三角矩阵，其对角元素为1,称为单位下三角矩阵，而三已是上三角矩阵。记L・£,]■",则有A-LU结果表明若消元过程可行，可以将三分解为单位下三角矩阵上与上三角矩阵力的乘积。由此派生出解方程组的宜接分解法。由高斯消元法得到启发，对E消元的过程相当于将工分解为一个上三角矩阵和一个下三角矩阵的过程。如果直接分解W得到上和二'，A-LUQ这时方程痴化为,令由。解出了；再由改解出工。这就是直接分解法。将方阵工分解为4.AU,当A是单位下三角矩阵，二.是上三角矩阵时，称为多利特尔（Doclittle）分解；当上是下三角矩阵，二是单位上三角矩阵时，称为库朗（Courant）分解。分解的条件是若方阵三的各阶主子式不为零，则多利特尔分解或库朗分解是可行和惟一的。多利特尔分解多利特尔分解步骤若W的各阶主子式不为零，三可分解为A其中心为单位下三角矩阵，二为上三角矩阵，即矩阵上和下共有『个未知元素，按照二的行元素上的列元素的顺序，对每个,列出式（3.16）两边对应的元素关系式，一个关系式解出一个工或二’的元素。下面是计算上和二’的元素的步骤：（1）计算二的第一行元素要计算则列出式（3.20）等号两边的第1行第1列元素的关系式:.1・各也＜■(10 0故一般地，由了的第一行元素的关系式Jav・■(10 0)得到

Vv-OqJ-1,2.--.*(2)计算上的第一列元素♦Ju要计算‘二】，则列出式(3.20)等号两边的第2行第1列元素的关系式:.1。…可,T故一般地，由£的第1列元素的关系式得到ifl-Ofl/Ujpt-N….“(3)计算-的第2行元素(4)计算上的第2列元素若已算出二”的前£一】行，上的前左-1列的元素，则(5)计算二的第2二行元素*里*由二的第2行元素的关系式:

、■江，7皿41…,3是上三角矩阵，"列标二二行标上。* A 1-4%N以0 •2以f•-1 ，■】 ，」得到口(3.21)气■、•L・・9m

z(3.21)(6)计算L的第C列元素j匕尸‘电由£的第2列元素的关系式：、・金5■4…41Tla…⑼蓝o,:L是下三角矩阵，二行标标：-行标:二o得到一直做到上的第,7列，二的第工行为止。用4U直接分解方法求解方程组所需要的计算量仍为3 ,和用高斯消元法的计算量基本相同。可以看到在分解中的每个元素只在式(3.21)或(3.22)中做而且仅做一次贡献，如果需要节省空间，可将二以及上的元素直接放在矩阵占相应元素的位置上。用直接分解法解方程4*首先作出分解/■幺“令解方程组“由于上是单位下三角矩阵，容易得到Mk (3.23)再解方程组及■，，其中时工进行4U分解时，并不涉及常数项因此，若需要解具有相同系数的•系列线性方程组时，使用直接分解法可以达到事半功倍的效果。多利特尔直接分解算法1.输入：方程组阶数3系数矩阵-二和常数项2FORk=lTOn{FORj4TOn//计算:的第七行元素FORi=kMTO«〃计算£的第士列元素完成分解////解方程组0-6////解方程组0-6//解方程组3.FORI-ITOn4FORi-ITOn5.输出方程组的解巧/例3.2用多利特尔分解求解方程组:对*-1,有«v-«vJ-U3«n-Z«R-l«u-14・•，2・23.-1/2-0.S^m-1/2-0.5对比-2,有«» -3-05-15«»"«»-的口-2-05-1.5外-圆-WiI)g-(2-05)/25-0.6对*-3,有♦fT—・Q6.』・4力・，力・06解皿・，，即-Lq-l3.2.2追赶法很多科学计算问题中，常常所要求解的方程组为三对角方程组其中(3.26)并且满足条件:Ml*#。lARk1*1*11^，。八1.2.…l4H，po称a为对角占优的三对角线矩阵。其特征是非零元素仅分布在对角线及对角线两侧的位置。在样条插值函数的“关系式中就出现过这类矩阵。事实上，许多连续问题经离散化后得到的线性方程组，其系数矩阵就是三对角或五对角形式的对角矩阵。利用矩阵直接分解法，求解具有三对角线系数矩阵的线性方程组十分简单而有效。现将工分解为下三角矩阵£,及单位上三角矩阵）的乘积。即‘其中利用矩阵乘法公式，比较人.工•“两边元素，可得到u,G=N"=2.3「避-1于是有!4.巧.吗.令，吗吟一q,_1/-2,W…,“(3.28)H叫/-L2,…/-I

(3.28)由些可见将三分解为£及了，只需计算值｝及｛呼｝两组数，然后解37,计算公式为:(3.29)M.工"Of-<W.i)/-2.3.--.M(3.29)再解“■’则得x.-XA-x.-XA-Xf4**m-L…21(3.30)整个求解过程是先由(3.28)及(3.29)求MLM)及E),这时i=12“工是，，追，，的过程，再由(3.24)求出｛&L这时1=®，…/是往回赶，故求解方程组(3.25)的整个过程称为追赶法。3.2.3对称矩阵的皿f分解对称正定矩阵也是很多物理问题产生的•类矩阵，正定矩阵的各阶主子式大于零。可以证明，若三正定，则存在下三角矩阵工，使/・££「，直接分解/的分解方法，称为平方根法。由于在平方根法中含有计算平方根的步骤，因此很少采用平方根法的分解形式。对于对称矩阵，常用的是人・£。上1'分解。

由不对称正定，可得% 令由d的对•称性和分解的惟一性可证£是对角元素为1的单位下三角矩阵。对矩阵工作多利特尔或库朗分解，共计算炉个矩阵元素；对称矩阵的UE分解，只需计算2个元素，减少了近一半的工作量。借助于多利特尔或库朗分解计算公式，容易得到分解计算公式。设t有多利特尔分解形式：]Pi也-小•山-?；•.心；呼其中4在分解可套用多利特尔分解公式，只要计算下三角矩阵上和总的对角元素d*。计算中只需保存,"3)的元素，三的「行/列的元素用工的L表示。由于对称正定矩阵的各阶主子式大于零，直接调用多利特尔或库朗分解公式可完3d分解计算，而不必借助于列主元的分解算法。对于1.1»2・一・4.有_J_J4・%-2初F-IM由UDlT^x~b，解方程组AxF可分三步完成：(1)解方程组&.，其中z…-加―(2)解方程组。其中,.dx。尤・〃4, 1-1,2(3)解方程局再**_习j ….1(3.32)JL4…,〃(3.33)(3.34)…” (3.35)(3.36)对称矩阵的LDU分解算法.输入：方程组阶数七，系数矩阵上和常数项X.FORk=1TOnd*：=oSdd：r-1FORi.HcMTOn.略去解方程组步骤从矩阵分解角度看，直接分解法与消元法本质上没有多大区别，但实际计算时它们各有所长。一般来说，如果仅用单字长进行计算，列主元消元法具有运算量较少、精度高的优点，故是常用的方法。但是，为了提高精度往往采取单字长数双倍内积的办法（即作向量内积计算时，采用双倍加法，最终结果再舍入成单字长数），这时直接分解法能获得较高的精度。3.4用LDC分解求解方程组：J28t.r1 o o|p'TOC\o"1-5"\h\z-1 I 0.0- 21-0.5I 39 \由有Lx.仇z-«-4冉?1,xmytx•-\2Li3.3范数向量范数在一维空间中，实轴上任意两点41*距离用两点差的绝对值|*b|表示。绝对值是一种度量形式的定义。范数是对函数、向量和矩阵定义的一种度量形式。任何对象的范数值都是一个非负实数。使用范数可以测量两个函数、向量或矩阵之间的距离。向量范数是度量向量长度的一种定义形式。范数有多种定义形式，只要满足下面的三个条件即可定义为一个范数。同一向量，采用不同的范数定义，可得到不同的范数值。定义3.1对任一向量/按照一个规则确定一个实数与它对应，记该实数记为帆1,若满足卜面上个性质：⑴有P1W当且仅当*・。时，（非负性）⑵Uek,aw虐，有（齐次性）」（3.37）⑶“心肥，有产♦平阳中（三角不等式）那么称该头数L■为向量上的氾数。几个常用向量范数向量的J范数定义为巴•序if其中，经常使用的是P・L2.三种向量范数。双■力讣附+NI++M»】wl,标+x：或写成m•国WL•靖{Id・g{kL同'>kO例3.5计算向量*・a・M_5f.P.L28的三种范数。设卜1*，*5-979161|X|L-max{l.X|-5|-5向量范数的等价性上两种不同的范数定有限维线性空间片中任意向量范数的定义都是等价的。若4⑶出⑶是上两种不同的范数定义，则必存在0<*<*《。,使¥2€肥均有

（证明略）向量的极限有了向量范数的定义，也就有了度量向量距离的标准，即可定义向量的极限和收敛概念了。设为X上向量序列，若存在向量使!m-十°,则称向量列心是收敛的（H是某种向量范数），仪称为该向量序列的极限。由向量范数的等价知，向量序列是否收敛与选取哪种范数无关。向量序列m.（物‘守""产），m=12…收敛的充分必要条件为其序列的每个分量收敛,即出者B存在。若妈/■角，则” …J就是向量序列（榕卜｛（产，近—y）的极限。解：算出每个向量分量的极限后得取格,为极限在计算方法中，计算的向量序列都是数据序列，当pL-ei小于给定精度时,向量*二取格,为极限3.3.2矩阵范数矩阵范数定义定义3.2如果矩阵,.（哂电广的某个非负实函数*3,记作风，满足条件：（1）Ml?。,当且仅当人・。时，°（非负性）⑵M_HM«tfe*（齐次性）⑶对于任意两个阶数相同的矩阵a»b有（三角不等式性）（4）46矩阵为同阶矩阵（相容性）则称M为矩阵范数。矩阵的算子范数常用的矩阵范数是矩阵的算子范数，可用向量范数定义：设AcL.Xe**,记方阵上的范数为阿，那么皆或 3（3.38）称为矩阵的算子范数或从属范数。这样定义的矩阵范数满足矩阵范数的所有性质外，还满足相容性:W为华阶矩阵，丫£wR■.恒有卬BMai （3.39）根据定义，对任一种从属范数有即单位矩阵的范数是1。常用矩阵范数向量有三种常用范数，相对应的矩阵范数的三种形式为：同,陶加《的行范数）（3.40）

（正的列范数）(3.41)（正的列范数）(3.41)（为是"4的最大特征值）（4的2范数）（3.42）证明：既然矩阵的算子范数是配上满足牍I.〕向量范数wi的上确界，那么，找到这个上确界也就找到了矩阵的范数。⑴任取万名必.设|・1,则血卜:MI/-4血卜:MI/-4喀快胴《卜加）副小踹加即14,踹加另一方面设极大值在上列达到，即端UI岁喇取入7。,0「・.0」.0「-，09,¥除第2分量为1外，其余分量均为0,于是有W卜隆.冬引r[・gjoj■能却J由于H",故因此有

⑵任取二…则皿■二皿工2卜播热旧卜即ML，福第J另一方面设极大值在上行达到，取”・•…产叫了这里{l.akO-l.a<0于是M办I故心㈱割(3)，)为对称非负矩阵，具有非负特征值，并具有v个相互正交的单位特征向量。设*'A的特征值为4 2— 相应的特征向量为足，其中聘为相互正交的单位向量。设XrgfAT*并旦Ml即不"，则j^AX■ +乙巧当♦…+斗!4"I"(山幽・，找㈤■之常j为4’-A1 M即对任意均有卬L"JZ,故于是如果A是对称矩阵，那么,4・44,设W的特征值是4,则有14・°«口-皿=卜1还有一种与向量范数1rL相容的矩阵范数，称为欧几里得（Euclid）范数或舒尔（Schur）范数，用Ml表示，其定义为因为欧几里得范数易于计算，在实用中是一种十分有用的范数。但它不能从属于任何一种范数，因为与向量范数的等价性质类似，矩阵范数之间也是等价的。?KKKM例3.7I,I,分别有 。解：-mae(卜1|♦3,2•17)■9pyL-皿吠卜1|+2»3*7)70的特征值4■依电4・281Pl-^019Ml-（1+4+9+49）%-屈矩阵范数等价性定理对l上的任两种范数Hi及H,存在常数。了三使矩阵特征值与范数关系若W是矩阵三的特征值（即存在非零向量者使得：AK=OX）,对任一算子范数有（相容性）.,MNMII即矩阵特征值的模不大于矩阵的任一范数。谱半径若3有特征值卬•记。（&-慵"1"则称。（'）为上的谱半径。有了谱半径的定义，矩阵的2范数可记为：风■山（，'）。谱半径与矩阵范数关系由矩阵谱半径定义，可得到矩阵范数的另一重要性质，。（⑷“ML在解方程组时，我们总是假定系数矩阵片和常数项上是准确的，而在实际问题中，系数矩阵A和常数项士往往是从前面的近似计算所得，元素的误差是不可避免的。这些误差会对方程组小.5的解2产生多大的影响？矩阵的条件数将给出一种粗略的衡量尺度。定义3.3若三非奇异称为三的条件数。其中■表示矩阵的某种范数。用矩阵V及其逆矩阵的范数的乘积表示矩阵的条件数，由于矩阵范数的定义不同，因而其条件数也不同，但是由于矩阵范数的等价性，故在不同范数下的条件数也是等价的。矩阵条件数的大小是衡量矩阵“坏”或“好”的标志。对于线性方程组若系数矩阵有小扰动W这时方程组的解也有扰动，于是\-Cottd"+3A）（X+6工）•8。34对，*的影响可表示为：\-Cottd(3.44)若常数上有小扰动.，其对$方的影响表示为:(3.45)(3.45)因此，C*（内大的矩阵称为“坏矩阵''或"病态矩阵对于式4）大的矩阵，小的误差可能会引起解的失真。一般说来，若上的按模最大特征值与按模最小特征值之比值较大时，矩阵就会呈病态。特别地，当修4很小时，/总是病态的。例3.8方程组2160^+01441X,-01440H2969jq+0耐呢-08M2方程组的解为对常数项上引入扰动热9SOOEMOgT,则解为化…]串”1]司[-0.48701即取・(-1.009,1.5130);虽然.很小，但解已完全失真。detX--Idr,^-1^(°8648…][-1,29690.2161)14,=21617卜、L=L5l3xl(f故.&2TO02lx|”,4的病态是显而易见的。.5病态方程组的解法如果系数矩阵W的条件数远大于1,则&T 为病态方程组，对病态方程组求解可采用以下措施：(1)采用高精度运算，减轻病态影响，例如用双倍字长运算；⑵用预处理方法改善后的条件数，即选择非奇异阵RQeL,使尸等价，而".侬的条件数比工改善，而求无■阳的解£・0"x•即才为原方程组的解，计算时可选择P・Q为对角阵或三角阵:

(3)平衡方法，当三中元素的数量级相差很大，可采用行平衡或列平衡的方法改善■的条件数。设加廿l非奇异,计门移令。-蛔于是求…等价于求&检.。6,或法■瓦这时孤-”的条件数可得到改善，这就是行均衡法。例3.9给定方程组为A的条件数Cg〃).7()f,若用行均衡法可取D"箍趴旷，1),则平衡后的方程九为3/而.*4,用三位有效数字的列主元消元法求解得nmimV3.6程序示例3.6程序示例程序3.1用列主元高斯消元法求解方程组:算法描述输入：方程组阶数n,矩阵a及列向量b。分解过程：对充・L2…选择q*・g{l%U%uLT%iQ交换第上行和第产;行方程对….熊jn….〃对f-A+l…回代过程：对…」•解得■夫•色-石仙)―程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII〃RT/KW•则主元离裔箱元相隔口»〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<stdio.h>#defineMAX-N20〃(xi,yi)的最大维数intmain(){intn;intij,k;intmi;doublemx,tmp;staticoublea[MAXN],[MAXN],b[MAXN],x[MAX_N];printf(44\ninputnvalue(dimofAx=b):");〃输入Ax=b的维数scanR"％d'',&n);ifi(n>MAXN){printf(4*TheinputnislargerthanMAXN,pleaseredefinetheMAXN.\n");return1;}if(n<=0){printfi(tfcPleaseinputanumberbetween1and%d.\n",MAX_N);return1;(〃输入Ax=b的A矩阵printf(4fcNowinputthematrixa(i,j),i,j=0,...,%d:\n”,n-1);for(i=0,i<n;j++)fbr(j=Ou<n;j-H-scanC%lf,&a[i][j])〃输入b向量printR"Nowinputthematrixb(i),i=0,…,%d:\n”,n-1)；fbr(i=O；i<n;i++)scanf(u%If”,&b[i]);fbr(i=O；i<n-l;i++)fbr(j=i+l,mi=i,mx=fabs(a[i][i]);j<n；j++)〃找主元素mi=j;mx=fabs(a[j][i]);{if(i<mi) 〃交换两行{tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=l;j<n;j++)(tmp=a[i][j];a[i][j]=a[mi][j]；a[mi][j]=tmp;))for(j=i+l;j<n;j++) 〃高斯消元{tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for=(k=i,k<n;k-H-)a|j][k]+=a[i][k]*tmp;})x[n-l]=b[n-l]/a[n-l][n-l]; 〃求解方程fbr(i=n-2;i>0;i—){x[i]=b[i];

x[i]-=a[i][j]*xU]；x[i]/=a[i][i];〃输出)〃输出printfpSlove…xj=\n")；fbr(i=0;i<n++)printf("％f\n,x[i「)；return0;)// EndofFile 计算实例用列主元高斯消元法水求解方程组2、ff-4

<勺+34Zr*■6r*2j^*2xi-5程序输入输出Inputnvalue(dimofAx=b):3Nowinputthematrixa(i,j),i,j=0,...211132122Nowinputthematrixb(i),i=0, 2:465Solve...x_i=1.0000001.0000001.000000程序3.2用库朗分解公式求解方程组:算法描述输入方程组阶数巴矩阵工及列向量入将矩阵工分解为£口，其中Ai£-勺?,4】。对土■对，■**7…上j，・%-2kJ%-®r£3NJH对i=l2…4对i=nja-I,…AJ程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII//Purpose:库朗分解解Ax=b//IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<stdio.h>#defineMAXN20 〃(xi,y_i)的最大维数intmain(){intn;inti,j,k;intmi;doublemx,tmp;statedoudlea[MAXN][MAXN],b[MAXN],x[MAXN],y[MAXN];statedoubleI[MAXN][MAXN],u[MAXN][MAXN];printfC'nlnputnvalue(dimofAx=b):");〃输入加小的维数scanfT%d”,&n);ifi(n>MAXN){printf(fc*TheinputnislargerthenMAXN,pleadefinetheMAXN.\n");return1;if(n<=0){printfi[fcTleaseinputanumberbetween1and%d.\n",MAX_N);retun1;〃输入AxW>的才矩阵printff'Nowinputthematrixa(ij),ij=O...，%d.\n,,,n-l);fbr(i=O;i<n;i-H-)for(j=O;j<n;j-H-)scanC4%If;&a[i][j]);〃输入b矩阵printfit^Nowinputthematrixb(i),i=0,...d:\n,,,n-l);fdr(i=O;i<n;i++)scanfi[44%If;&b[i]);for(i=0;i<n;i++)u[i][i]=l;//U矩阵对角元素为1for(k=0;k<n;k++){for(i=k;i<n;i++)〃计算上矩阵的第三列元素(l[i][k]=a[i][k];For(j=0;j<=k-lJ-H-,l[i][k]-=(l[i][j]*u]j][k])；)for(j=k+l；j<n;j++) 〃计算二•矩阵的第匕行元素(u[k][j]=a[k][j];fbr(i=O;i<=k-l;i++)u[k]U]-=(l[k]U]*u[i]UD；for(i=0;i<n;i++)//Ly=b计算y{y[i]=b[i]；for(j=0;i<=i-l;j++)y[i]-=(i[i]U]*yU]);y[i]/=l[i][i];(fbr(i=n-1;i>=0;i-)//Ux=y计算x{x[i]=y[i];for(j=i+l;j<n;j++)x[i]-=(u[i][j]*x[j]);)printf("Solve...x」=\n")；〃输出for(i=0;i<=n;I++)printf(11%return0;}计算实例用库朗分解公式求解方程组：程序输入输出Inputnvalue(dimofAx=b):3Nowinputthematrixa(ij)ij=O,...,2：121-2-1-50-16Nowinputthematrixb(i),i=0,...,2：24-6350solve…x-i=7.0000004.0000009.000000程序3.3用Sf分解求解方程组:算法描述输入矩阵/及列向量5。将矩阵工分解为人・皿,，其中解小・A,先解£Z•凯再由0■”得I-12,-,M最后由/x•尸得a4.乂・就"%，'rx-1…21程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIPurpose:LD6分解解方程组〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<math.h>#defineMAX-N20〃(x・i,y・i)的最大维数intmain(){intn;inti,j,k;ilntmi;doublemx,tmp;staticdoublea[MAX-N][MAX-N],B[MAX-N],x[MAX-N],y[MAX-N]