《数值计算方法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
数值计方法1.1数值计算方法的对象与特点1.1.1什么是数值计算方法现代的科学技术发展十分迅速，他们有一个共同的特点，就是都有大量的数据问题。比如，发射一颗探测宇宙奥秘的卫星，从卫星设计开始到发射、回收为止，科学家和工程技术人员、工人就要对卫星的总体、部件进行全面的设计和生产，要对选用的火箭进行设计和生产，这里面就有许许多多的数据要进行
1准确的计算。发射和回收的时候，又有关于发射角度、轨道、遥控、回收下落角度等等需要进行精确的计算。有如，在高能加速器里进行高能物理试验，研究具有很高能量的基本粒子的性质、它们之间的相互作用和转化规律，这里面也有大量的数据计算问题。计算问题可以数是现代社会各个领域普遍存在的共同问题，工业、农业、交通运输、医疗卫生、文化教育等等，各行各业都有许多数据需要计算，通过数据分析，以便掌握事物发展的规律。研究计算问题的解决方法和有关数学理论问题的一门学科就叫做计算方法。计算方法属于应用数学的范畴，它主要研究有关的数学和逻辑问题怎样由计算机加以有效解决。1.1.1数值计算方法的内容数值计算方法也叫做计算数学或数值分析。数值计算方法主要内容包括非线性方程求根、线性代数方程组解法、微分方程的数值解法、插值问题、函数的数值逼近问题、概率统计计算问题等等，还要研究解的存在性、惟一性、收敛性和误差分析等理论问题。我们知道五次及五次以上的代数方程不存在求根公式，因此，要求出五次以上的高次代数方程的解，一般只能求它的近似解，求近似解的方法就是数值分析的方法。对于一般的超越方程，如对数方程、三角方程等等也只能采用数值分析的办法。怎样找出比较简洁、误差比较小、花费时间比较少的计算方法是数值分析的主要课题。在求解方程的办法中，常用的办法之一是迭代法，也叫做逐次逼近法。迭代法的计算是比较简单的，是比较容易进行的。迭代法还可以用来求解线性方程组的解。求方程组的近似解也要选择适当的迭代公式,使得收敛速度快，近似误差小。在线性代数方程组的解法中，常用的有塞德尔迭代法、共辗斜量法、超松弛迭代法等等。此外，一些比较古老的普通消去法，如高斯法、追赶法等等，在利用计算机的条件下也可以得到广泛的应用。在数值计算方法中，数值逼近也是常用的基本方法。数值逼近也叫近似代替，就是用简单的函数去代替比较复杂的函数，或者代替不能用解析表达式表示的函数。数值逼近的基本方法是插值法。初等数学里的三角函数表，对数表中的修正值，就是根据插值法制成的。在遇到求微分和积分的时候，如何利用简单的函数去近似代替所给的函数，以便容易求微分和求积分,也是计算方法的一个主要内容。微分方程的数值解法也是近似解法。常微分方程的数值解法有欧拉法、预测校正法等。偏微分方程的初值问题或边值问题，目前常用的是有限差分法、有限元素法等。有限差分法的基本思想是用离散的、只含有限个未知数的差分方程去代替连续变量的微分方程和定解条件。求出差分方程的解法作为求偏微分方程的近似解。有限元素法是近代才发展起来的，它是以变分原理和剖分差值作为基础的方法。在解决椭圆形方程边值问题上得到了广泛的应用。目前，有许多人正在研究用有限元素法来解双曲形和抛物形的方程。数值计算方法的内容十分丰富，它在科学技术中正发挥着越来越大的作用。1.1.2数值计算方法的特点第一，面向计算机，要根据计算机特点提供实际可行的有效算法。即算法只能包括加、减、乘、除和逻辑运算，是计算机能直接处理的。
2第二，有可靠的理论分析，能任意逼近并达到精度要求，对近似算法要保证收敛性和数值稳定性，还要对误差进行分析。第三，要有好的计算复杂性，时间复杂性好是指节省时间，空间复杂性好是指节省存储量，这也是建立算法要研究的问题，它关系到算法能否在计算机上实现。第四，要有数值实验，即任何一个算法除了从理论上要满足上述三点外，还要通过数值实验证明是行之有效的。根据“数值计算方法''的特点，学习时我们要注意掌握方法的基本原理和思想，要注意方法处理的技巧及其与计算机结合，要重视误差分析、收敛性及稳定性的基本理论；其次，要通过例子，学习使用各种数值方法解决实际计算问题；最后，为了掌握本课的内容，还应做一定数量的理论分析与计算问题练习。由于本课内容包括了微积分、代数、常微分方程的数值方法，以及高级程序设计的方法，读者必须掌握这几门课的基本内容才能学好这一课程。1.1误差来源与误差分析1.1.1误差的来源早在中学我们就接触过误差的概念，如在做热力学实验中，从温度计上读出的温度是23.4℃就不是一个精确的值，而是含有误差的近似值。事实上，误差的在我们的生活中无处不在，无处不有。如量体裁衣,量与裁的结果都不是精确无误的，都有误差。人们可能会问：如果使用计算机来解决问题，结果还会有误差吗？下面我们通过考察数学方法解决实际问题的主要过程来思考这个问题。用数学方法解决一个具体的实际问题，首先要建立数学模型，这就要对实际问题进行抽象、简化，因而数学模型本身总含有误差，这种误差叫做模型误差。在数学模型中通常包含各种各样的参变量，如温度、长度、电压等，这些参数往往都是通过观测得到的，因此也带来了误差，这种误差叫做观测误差。当数学模型不能得到精确解时，通常要建立一套行之有效的数值方法求它的近似解，近似解与准确解之间的误差就称为截断误差或方法误差。由于在计算机中浮点数只能表示实数的近似值，因此用计算机进行实际计算时每一步都可能有误差，这种误差称为舍入误差。
3例如，函数A划用泰勒(Taylor)展开式近似代替，则数值方法的截断误差是/*(*+1)m卷卅】(孑介于0与X之间)又如，在计算时用3.14159近似替代"产生的误差R=乃-3.14159=0.0000026…就是舍入误差。上述种种意味着都会影响计算结果的准确性，因此需要了解与研究误差。在数值计算方法中将着重研究截断误差、舍入误差，并对它们的传播与积累作出分析。1.1.1绝对误差、相对误差与有效数字1.绝对误差与绝对误差限定义1.1设x为准确值，丁为x的一个近似值，称G=x'-x为近似值的绝对误差，或误差。通常我们无法知道准确值x,也不能算出误差的准确值,，只能根据测量或计算估计出误差的绝对值不超过某个正数£二则一称为绝对误差限。有了绝对误差限就可知道x的范围/-f+£'，即x落在区间［X'-£'，/+内。例如用亳米测度尺测量一长度X,读出的长度为23mm,则有一405mm。由此例也可以看到绝对误差是有量纲和单位的。2.相对误差与相对误差限只用绝对误差还不能说明数的近似程度，例如甲打字时每百个字错一个，乙打字时平均每千个字错一个，他们的误差都是错一个，但显然乙要准确些。这就启发我们除了要看绝对误差大小外，还必须顾及量的本身。定义1.2把近似值的误差g与准确值x的比值称为近似值/的相对误差，记作e，gy=--实际计算时，由于真值X通常是不知道的，通常取'/。相对误差也可正可负，它的绝对值的上.xerS=Ivo界叫做相对误差限。记作号。即II。根据定义，甲打字时的相对误差10°,乙打字时的e；4—5—=0.1%相对误差1000o易知相对误差是一个无量纲量。
41.有效数字当x有多位时，常常接四舍五入的原则得到x的前几位近似值/，例如x==3.14159265--取3位4=3.12440.002；取5位x；=3.1416£«0.00005：它们的误差都不超过末位数字的半个单位，即11|^-3.14|i-xl0-2,|^-3,1416|^-xl0M现在我们将四舍五入抽象成数学语言，并引入一个新名词“有效数字”来描述它。定义L3若近似数/的误差限是某一位的半个单位，该位到丁的第一位非零数字共有力位，我们就说/有万位有效数字。如取x.=3.14作乃的近似值，,就有3位有效数字；取大=3.1416作兀的近似值，/就有5位有效数字。有效数字也可采用以下定义:X的近似数X*写成标准形式:x=+10*x(a1xl0-1+(32xl0-2+-+a],xio-1')其中…%是o到9的一个数字，且,不为0,化为整数，若|x-x*|<-xlOfc^112（1.2）则称:有力位有效数字。例1.1依四舍五入原则写出下列各数具有5位有效数字的近似数，913.95872,39.1882,0.0143254,8.000033解：按定义，上述的5位有效数字的近似数分别为913.96,39.188,0.014325,8.0000
5注意，8.000033的5位有效数字的近似值是8.0000而不是8,8只有一位有效数字。4.有效数字与误差关系不难看出，若由（1.1）给出某近似数有片位有效数字，则可以从（1.2）求得这个数的绝对误差限E=1x10*-"2因此，在上相同的情况下，附越大则10"-"就越小，故有效数字越多，绝对误差限越小。定理1.1用（1.1）表示的近似数个，若/具有*位有效数字，则其相对误差限为.——x.反之，若X’的相对误差限2（生+1）则Y至少具有〃位有效数字。证明由式（1.1）得
6,X10"】《卜］«(%+1)*1。1所当X有〃位有效数字时LlOi1Z——b=_Lxl0TT%xlO»i2al反之，若X的相对误差限2(的+1)xlO-'f有:xlO-"+1=2*10i2k-4=卜'W(/)431+1)X10"1X_由式(1.2)式知道，X’至少有片位有效数字，证毕这说明有效数字越多，相对误差限越小。例1.2要使而的近似值的相对误差限小于°」％,要取几个有效数字？解：因为痴=447…，所以，=4,令—xlQ-*+1<0.1%2x4故取〃=4即可满足。4.函数计算的误差估计一般情况，当自变量有误差时计算相应的函数值也会产生误差，其误差限可利用函数的泰勒展开式估计。设“X)是一元函数，x的近似值为以“X)近似/(X),其误差界记作,可用泰勒展开式/(x)-/(x)=/(x)(x-+(x-X*)3。介于X,父之间，取绝对值得
7|/（x）-4|7d）|£（x）+.‘Is2（x*）假定了'（x）与」“CO的比值不太大，可忽略£（x）的高阶项，于是可得计算函数值的误差限
8I"'I101W.试求面积的绝对误差限与相对误差限。s=ld,—=d,—=l解：因第3d,那么6糕轲）+部，）其中（§*=/=80/m,（—）**/*-110w£（d*）=0,lw,£（Z*）=0.2附于是绝对误差限£（3）、80x0.2+110x0.1=27（加2）相对误差限1.3误差传播与若干防治办法由前述可知，在数值计算中每步都可能产生误差。而一个问题的解决，往往要经过成千上万次运算，我们不可能每步都有加以分析，下面，通过对误差的某些传播规律的简单分析，指出在数值计算中应注意的几个原则，它有助于鉴别计算结果的可靠性并防止误差危害现象的产生。1.要避免两相近数相减在数值运算中两相近数相减会使有效数字严重损失，例如x=532.65,=532.52都具有五位有效数字，但x-y=013只有两位有效数字，所以要尽量避免这类运算。IgXj-lgXj=lg—通常采用的方法是改变计算公式，例如当X】与町很接近时，由于町那么可用右端的-V^=-r=-7=公式代替左端的公式计算，有效数字就不会损失。当X很大时，.x+1+yx也可用右端来代替左端。一般情况，当/(x)=/(x)时，可用泰勒展开/(X)-/(/)--X*)+-x*)a+- 2!取右端的有限项近似左端。
9如果计算公式不能改变，则可采用增加有效位数的方法。2.要防止大数“吃掉”小数若参加运算的数的数量级相差很大，而计算机的位数有限，如不注意运算次序，就可能出现大数“吃掉”小数的现象，影响计算结果。例如在五位十进制计算机上，计算1000工=52492+2012-1写成规格化形式10005A=0.52492*105+^0,000001x102-1由于计算时要对阶，000001*105在计算机中表示为0,计算出来工=0.52492*1()5,结果严重失真!1000殳】如果计算时，先将I计算出来，再与52492相加，就不会出现大数“吃掉”小数的现象了。3.注意简化计算步骤，减少运算次数同样一个计算题，如果能减少运算次数，不但可节省计算机的计算时间，还能减少舍入误差。例如，计算x255的值，如果逐个相乘要用254次乘法，但若写成X=X，X'X'X'X'XX只要做14次乘法运算即可。4.绝对值太小的数不宜作除数♦匚♦X♦XX同，Z=—Z=—设x与y分别有近似值沙的近似值y,其绝对误差显然，当1丁I很小时，近似值Z,的绝对误差e（Z）有可能很大。因此，不宜把绝对值太小的数作除数。
10第2章非线性方程求根与线性方程相比，非线性方程问题无论是从理论上还是从计算公式上，都要复杂得多。对于一般的非线性方程计算方程的根既无一定章程可寻也无直接法可言。例如，求解高次方程组7/-/+矛75=0的根，求解含有指数和正弦函数的超越方程/8$9乃=°的零点。解非线性方程或非线性方程组也是计算方法中的一个主题。一般地，我们用符号/（力来表示方程左端的函数，方程一般形式表示为,（x）=°,方程的解称为方程的根或函数的零点。通常，非线性方程的根不止一个，对于非线性方程，一般用迭代法求解。因此，在求解非线性方程时，要给定初始值或求解范围。2.1实根的对分法2.1.1使用对分法的条件对分法或称二分法是求方程近似解的一种简单直观的方法。设函数,（x）在句上连续，且/（a）/。）＜0,则/（x）在［a1］上至少有一零点，这是微积分中的介值定理，也是使用对分法的前题条件。计算中通过对分区间、缩小区间范围的步骤搜索零点的位置。例2.1用对分法求解/口）=#-7.7/+19.2钎15.3在区间［1,2］之间的根。解：⑴/⑴=-2.8,/⑵=0.3,由介值定理可得有根区间［。闻=口，2］。
11(2)计算为=(1+2)/2=1.5,/(1.》=-0.45,有根区间【。肉=口,5,2]。(3)计算叼=(1.5+2)/2=175,0=0.078125,有根区间1°向=(4)一直做到(计算前给定的精度)或1。一外〈£时停止。2.1.1对分法求根算法计算/(x)=°的一般计算步骤如下：1.输入求根区间SI1和误差控制量£,定义函数/(X)。IF>0)THEN退出选用其他求求根方法2.W///LEI。力1>£时(1)计算中点x=3+切/2以及/(x)的值；(2)分情况处理停止计算J=x,转向步骤4<0修正区间[a,x]~>[a,b]<0修正区间[x,句->[a,b]ENDWHILE&xs(a+/2J•o4.输出近似根亡。
12(2)2x=9-5,迭代格式为1怎-52在算法中，常用gg'c/gAag应/炽））＞0代替/（。）・/。）＞0的判断，以避免了3）,/（力数值溢出。对分法的算法简单，然而，若了5）在是有几个零点时，只能算出其中一个零点；另一方面，即使丁（X）在［虱加上有零点，也未必有,“）/。）＜°。这就限制了对分法的使用范围。对分法只能计算方程/⑶=0的实根.2.2迭代法对给定的方程,（*）=°,将它转达换成等价形式：x=W（x）。给定初值两，由此来构造迭代序列%「联）次=1,2,…，如果迭代法收敛,即明“1-原双凝）一:有&=岫），则a就是方程/（"）=0的根。在计算中当।-凝।小于给定的精度控制量时，取a=*切为方程的根。例如，代数方程芯3-2x-5=0的三种等价形式及其迭代格式如下:(1)xX迭代格式=2x+5,x=*2x+5,迭代格式Xfc+i=/凝+5
1324+5(3)3nu2x+5rx=2x+5,x=——-%
14对于方程」（力=0构造的多种迭代格式=飙凝），怎样判断构造的迭代格式是否收敛？收敛是否与迭代的初值有关？定理2.1P（x）若定义在【么句上，如果W（x）满足（1）当有xe［a,切有a4P（X）4九（2）P（x）在［凡句上可导，并且存在正数,使对任意的xe【a,可，有|p'（x）|4£：则在［。，句上有惟一的点/满足,=忒,），称/为W。）的不动点。而且迭代格式加1=°（秋）对任意的初值/丘［明3均收敛于P。）的不动点寸，并有误差估计式证明：（1）先证明存在性：令/（X）='-°（x）,则有/(a)=a-p(a)幺0-00)20/(/)=x-。(产)=0或x=0(高（2）再证明惟一性：设都是W（x）的不动点，且则有~x“|=|穴X*)-W(x“)I=|-0(/-x")|4£I，-x”|<|x-x**\Je[a,b]与假设矛盾，这表明f=x”,即不动点是惟一的。（3）当力]时，由于W（x）e【a，可可用归纳法证明，迭代序列SJu[a力],于是由微分中值定理耳+1-必=
15。（4）一«亡）=0（与（/一亡），^[a,b]和Id（x）1«£,得I一x'K川凝-x'卜£I/“）-W）I-户I加1-xI-严11两-xI（2.2）因为ZX1,所以当上一8时，乃—X即迭代格式马+1=44）收敛。（4）误差估计式：片+1一耳|=|奴/）一«4-]）区©/一4-」…“I玉-引设々固定，对于任意的正整数「有，I入+p-Xk1=14+P一4+;HI+I七+;H~Xk+p-2I幺d+户Z+…为|X「而|=Ix「与I1-£limk=x*由于。的任意性及，T9，故有.ilim%=x*Ix-&K---I々-%I"T9>[-L注：定理2.1是判断迭代法收敛的充分条件，而非必要条件。要构造满足定理条件的等价形式一般是难于做到的。事实上，如果X•为/（X）的零点，若能构造等价形式X=W（x）,而I|<1,由P'（x）的连续性，一定存在父的邻域（/-2，/+同，其上有
16|w'（x）|<£ 17（2）校正：才用=吊）▼_=一（赤一褊）2激+1_二一"二yz"~（3）改进：％1-/敌+1+醺有些不收敛的迭代法，经过埃特金方法处理后，变得收敛了。例2.3求方程/。）=工3-矛7=0在而=1.5附近的根父。解：若采用迭代公式：加1=北7迭代法是发散的，我们现在以这种迭代公式为基础形成埃特金算法：取勺=15计算结果如表2」所示:表2.1计算结果k■敢+1%01.512.3750013.39651.4162921.840925.238881.3556531.491402.317281.3289541.347101.444351.3248051.325181.327141.32472我们看到，将发散的迭代公式通过埃特金方法处理后，竟获得了相当好的收敛性。2.3牛顿迭代法对方程/（x）=0可构造多种迭代格式=°（醺），牛顿迭代法是借助于对函数/（x）=°作泰勒展开而构造的一种迭代格式。将/。）=°在初始值与作泰勒展开：/。）=/（而）+/'（/）。-工0）+，（X-飞）2+-- 18取展开式的线性部分作为/（X）"0的近似值，则有设/%"°,则L一」（%）A-An了也）x=yw令r（x0）类似地，再将/（X）=0在五作泰勒展开并取其线性部分得到：x-X」⑷21一直做下去得到牛顿迭代格式：牛顿迭代格式对应于/0）=0的等价方程是 19X= 20图2.2牛顿切线法示意图在区域口。，号+加的局部“以直代曲”是处理非线性问题的常用手法。在泰勒展开中，截取函数展开的线性部分替代“X）。例2.4用牛顿迭代法求方程7.7x2+19.2x-15.3,在5=1附近的根。x?-7.7xj+19.2xl-15.3x=》*.*解：“'3元75.4&+19.2计算结果列于表2.2中。表2.2计算结果Kxk/W01.00-2.811.41176-0.72707121.62424-0.14549331.6923-0.013168241.69991-0.000151551.70比较表2.1和表2.2的数值，可以看到牛顿迭代法的收敛速度明显快于对分法。牛顿迭代法也有局限性。在牛顿迭代法中，选取适当迭代初始值瓦是求解的前题，当迭代的初始值两在某根的附近时迭代才能收敛到这个根，有时会发生从一个根附近跳向另一个根附近的情况，尤其在导数/'（"）数值很小时，如图2.3所示。图2.3失效的牛顿迭代法 21牛顿迭代算法1.定义函数,(x),g(x),输入或定义迭代初始值X—旦)和控制精度epsilon。2.FORi二°toMAXREPTx_il:=g(x_^0)〃g(x)是牛顿迭代公式ZF(|x_k\-x_kO\ 22在牛顿迭代格式中:力”]=/（—）…用差商凝一X—代替导数/（演），并给定两个初始值X。和不，那么迭代格式可写成如下形式：X=x32…人JtyXl人1,R,，乙/（凝）-/（%）（2.5）上式称为弦截法。用弦截法迭代求根，每次只需计算一次函数值，而用牛顿迭代法每次要计算一次函数值和一次导数值。但弦截法收敛速度稍慢于牛顿迭代法。弦截法的几何意义做过两点SoJ（x。））和（占'/（々））的一条直线（弦），该直线与X轴的交点就是生成的迭代点*2,再做过（0/（々））和（叼，,（工2））的一条直线，与是该直线与X轴的交点，继续做下去得到方程的根/3）=0,如图2.4所示图2.4弦截法不意图例2.5用弦截法求方程"X）=--7.7,+19.2x75.3的根,取勒=1.5,々-4.0敌+1—敌解：/（xfc）-/（xw）计算结果列于表2.3中。表2.3计算结果k凝/（X） 2301.5-0.45142.321.909090.24883531.65543-0.080569241.717480.028745651.701160.0019590261.69997-0.000053924671.79.45910-8弦截法算法1.定义函数“X）,输入或定义控制精度值epsilon和迭代初始值J肛。计算/1：=」。_妇)x——x/2x表示X，凝2.FOR:=0,1...MAXREPT2.1/2:=/。_左2)22x_k:=x_k2-J2(x_k2-x_^l)/(/2-/I)/齐-｛表示0+123叩豆=OK(|/(x_熄|<叩面匕^THEN｛输出满足给定精度的近似解x-上，结束｝ 24/1：=/2x_kl-x_k2x_k2:-x_k〃为下一次迭代准备数值ENDFOR3.输出：在初始值X-妇，工一上2附近/（X）无根。2.6非线性方程组的牛顿方法为了叙述简单，我们以解二阶非线性方程组为例演示解题的方法和步步骤，类似地可以得到解更高阶非线性方程组的方法和步骤。设二阶方程组：工（“）=01东”）=0Q6）其中%丁为自变量。为了方便起见，将方程组写成向量形式:尸（W）/（X，"其中=将办（xj）,工f（xj）附近作二元泰勒展开，并取其线性部分，得到方程组k（x“0）+（x-x。）咨型+0»）吟®=。oxay、:、“2（XoJo）＞、历（为,北）「力（X。,为）+（x-X。）——+0-%）——-——=0dxdy令X-X。则有’及/（XoJo）+协通Jo）dxdy&羽（x0,既）上年电。0，%）dxdy=一力（丽jo）一力（而Jo） 25l（XoJo）l=如果皿笠dx血生生dxdy解出反，与，得&1%+8,△yj‘o+与」再列出方程组：'吟凶。F）+吟®8-%）一工（勺加dxay‘羽（须,％）,、_1.电（X1J1）/、，/、一3（X-X。+——7——⑶-M）=一力区,71）dxdy解出&=xfa=厂必，得出W2继续做下去，每一次迭代都是解一个方程组J(X"1-力（X"G(2.7)为止。即加「4+&,然a=4+与，直到max（|A4|®|）＜例2.6求解非线性方程组1（0）=4--/彳—（"）=]_"-y=’11用=取初始值（T7人 26、J-2y72p--flIk=.以一ayI打打一有菽明以小脩-23.4-2.71828-1F(w(j)XCWo))JO.ll『式Xo，yJ-0.01828L/LJ-2&+3.44=-0.11(-2,71828Zto-2>=0.01828△w='△x[f0,0042560.004256-0.029849」 27解方程得1004256-0.029849-1729849继续做下去,直到max（心|,|今|）＜10-5时停止。将两个变量的非线性方程组推广到多变量的非线性方程组:‘工(向，X2「・，勺=0fi(x1,x3,-xx=0(2.8)记x=（々,叼,…,x”）rF（X）=（工。1，项,…1ai，々，…,x*）,…,彳历程…,A））'的向量形式:斤⑶=0 28用牛顿迭代法求方程组的解为:那+1)=xw__巴?=X(E)-(J(X3))-%(1闺)F(X')见⑶巩⑻...巩⑺的dx28A电(X)电⑻电⑻J(X)=加的网姒5筑⑻.筑(制其中：遍妈四称J(X)为雅可比(jacobi)矩阵。计算中采用J(X⑻)(**】)-*对)=-F(X侬)“犷乂型=-尸(X产c°、即划=X(k)+LXW一直做到"△‘'△’L小于给定精度为止。在X的领域中若UJ(㈤1匕"£<1,则迭代收敛。.7程序示例程序2.1用牛顿迭代法求/(X)在X。附近的根。算法描述给定了(X),从勺开始，根据牛顿迭代公式k=0,1，…XM“xj 29计算/(")在X。附近的根。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII//Purpose:牛顿迭代法求根〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#inckide 30”,x_kl）;〃满足精度，输出return; 31{x_k=x_k1{printf（44After%drepeate,nosolved.\m",MAXREPT）；//EndofFile计算实例已知/（x）=x3_3x7,取£=勺=1.5,用牛顿迭代法计算/（x）的根。程序输入输出由本程序预设的/（X）及飞=15,得到输出!Root:1.879385程序2.2用弦截法求,（X）在飞，々附近的根。算法描述给定了（x）,从丽，X1开始,根据弦截法迭代公式X=x_/（凝）®-%）AJt+i4//、//、，&［，乙，…-，（k1）求得了（X）在其附近的根。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII//Purpose：弦截法求根〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 32#include 33M,x_k2);x_k2=x_kl-(f(x_kl)*x_kl-x_k))/(f(x_k1)・f(x_k));〃弦截求新x_niffabs(x_k2-x_kl) 34,\x_k2);〃满足精度，输出return;{x_k=x_k1;x_kl=x_k2;〃反复(printf(i4After%drepdate,nosolvde. 35",MAXREPT);//EndofFile计算实例 36用弦截法求方程“X）Rx?-7.7x?+19.2x75.3的根，取x°=1.5,Xj=4.0程序输入输出由本程序的/（X）及飞,々得到输出!Root:1.700000第3章解线性方程组的直接法在近代数学数值计算和工程应用中，求解线性方程组是重要的课题。例如，样条插值中形成的河关系式，曲线拟合形成的法方程等，都落实到解一个万元线性方程组，尤其是大型方程组的求解，即求线性方程组（3.1）的未知量为瓜2,…，々的数值。‘勺丙+a12x2+--+auxn=4白21再++•••+a2Mx=%d%为4/++am=bx其中mjbi为常数。上式可写成矩阵形式Ax=b,即其中，"=（%）为系数矩阵，x=（孙叼…，/）r为解向量，力=（"也…也）,为常数向量。当加丛=。工0时，由线性代数中的克莱姆法则，方程组月入=小的解存在且惟一，且有Di为系数矩阵月的第i列元素以小代替的矩阵的行列式的值。克莱姆法则在建立线性方程组解的理论基础中功不可没，但是在实际计算中，我们难以承受它的计算量。例如，解一个100阶的线性方程组，乘除法次数约为（101100!-99）,即使以每秒103的运算速度，也需要近100小年的时间。在石油勘探、天气预报等问题中常常出现成百上千阶的方程组，也就产生了各种形式方程组数值解法的需求。研究大型方程组的解是目前计算数学中的一个重要方向和课题。解方程组的方法可归纳为直接解法和迭代解法。从理论上来说，直接法经过有限次四则运算，假定每一步 37运算过程中没有舍入误差，那么，最后得到方程组的解就是精确解。但是，这只是理想化的假定，在计算过程中，完全杜绝舍入误差是不可能的，只能控制和约束由有限位算术运算带来的舍入误差的增长和危害，这样直接法得到的解也不一定是绝对精确的。迭代法是将方程组的解看作某种极限过程的向量极限的值，像第2章中非线性方程求解一样，计算极限过程是用迭代过程完成的，只不过将迭代式中单变量入3=°（底）换成向量而已。在用迭代算法时，我们不可能将极限过程算到底，只能将迭代进行有限多次，得到满足一定精度要求的方程组的近似解。在数值计算历史上，直接解法和迭代解法交替生辉。一种解法的兴旺与计算机的硬件环境和问题规模是密切相关的。一般说来，对同等规模的线性方程组，直接法对计算机的要求高于迭代法。对于中等规模的线性方程组5<200）,由于直接法的准确性和可靠性高，一般都用直接法求解。对于高阶方程组和稀疏方程组（非零元素较少），一般用迭代法求解。3.1消元法3.1.1三角形方程组的解形如下面三种形式的线性方程组较容易求解。对角形方程组’的内=瓦。22心=%、2*=b*（3.3）设由L0,对每一个方程，号=12…,”显然，求解〃阶对角方程的运算量为°（功。下三角方程组 38%再=瓦’21再+’22勺=&2+射叼+…+=bn(34)按照方程组的顺序，从第一个方程至第〃个方程，逐个解出玉/=12,…/由方程之为=4,得々=441。将々的值代入到第二个方程41X1+Z22X2=b2X?=02-JXi)〃22将々，盯，…,Xj_i的值代入到第i个方程—1+…=42-1玉=(4-£%Xj)4，i=1,2,…小得川计算为需要i次乘法或除法运算，,=1，2,….。因此，求解过程中的运算量为2?=-z-=6")2-1乙上三角方程组‘％1玉+uux2+•••+m=a“2/+…+2・=%4…4-“吊_1+%L*Xa=4_i,2,1"x*X*-4(35)与计算下三角方程组的次序相反，从第%个方程至第一个方程，逐个解出玉，由第力个方程"皿4=4彳可。将/的值代入到第力"I个方程 39U»-1>»-1X»-1+ux-l>xXx=耳-1得-（4-1-"*-1，*，*将\,勺_1，…，为+1的值代入到第j个方程%肉+%3k1+…+%*勺=4得解的通式■々=(4-Z%Xj)/%,/-i+1计算不需要非T+1次乘法或除法运算'=n'n-L…，2」。因此求解过程中的运算量为-4-1（，力（力+1）2、N力-1+1=2，=---=0（力）Ij-i乙消元法的基本思想就是通过对方程组做初等变换，把一般形式的方程组化为等价的具有上述形式的易解方程组。3.1.2高斯消元法与列主元消元法高斯消元法高斯消元法是我们熟悉的古老、简单而有效的解方程组的方法。下面是中学阶段解二元方程组（高斯消元法）的步骤：Xj-x2=1"6）’3々+2必=8（3.7）方程（3.6）乘以-3加到第（3.7）个方程中得%-Xa=1"0+5x2=5 40‘叼=L代入（3.6）得再=2。其方法相当于对方程组的增广矩阵做行的初等变换：/、fl-111（1-11）L、 41~瓦司»XM可当游一12*崂a~aa1AIX1*~的52瓦-~,6X可以看出，I'3实际上是由（A切按一系列初等换后得到的（1）对换（A"）某两行元素；（2）（A①中的某行乘一个不为零的数；（3）把（A与的某一行乘一个常数后加到另一行。高斯消元法就是通过以上（3）的变换，把（Aa化为等价的上三角形式。下面我们以％=4为例演示消元过程。设方程组：anXx+/2町+a13X3+，4勺=瓦+白22才2+a23X3+424工4=力20丙+。32灼+&3勺+。34勺=.产4丙+%2町+。43勺+。44乙=瓦（39：其增广矩阵为：an(3〜°31，2ana14瓦。22。23。24力2ana33。34力3a/。43。44“」（1）若，1工°,则将第一行乘以一。21ali加到第二行上；将第一乘以一。31/%1加到第三行上；将第一行乘以一41,/I加到第四行上得到 42(3.10)xr/\JZv)z勾Q2Q30414⑴24⑴34⑴“aaaa13⑴空⑴33⑴43aaaa^11000X⑴/即*=附一%%/。11£)=2,3,4防”=曲一瓦沏/。“/=2,3,4(2)若。黑"°，则将第二行乘以一%，/?加到第三行上；将第二行乘以-*/*加到第四行上,得到1/瓦培6iA»14(1)“⑵34⑵44aaaa旬艰瞪南小asoohlooo(3.11)。)Q)Q)Q)/Q),.oj白炉°炉"2/白口:。22，】，J$，48,)=可。)-珊a*/眼,i=3,4⑶若洲"°则将第三行乘以一滤"。翼加到第四行上，得到YJZ).).与Q2S04%⑴24⑵34⑶44aaaa小.(1)/⑵330aaa=(调(3.12)A已是上三角矩阵，于是得到了与原方程等价的易解形式的方程组: 43。11入1+anX2+al3x3+%4为=瓦吟+。氏+&4=.-%+思左4=理胆=“.444(3.13)再对方程组（3.13）依次回代解出々，/R2，々，。从式（3.12）可以得到系数矩阵工的行列式的值为彳的对角元素的乘积。即det（H）=/次gag）a3）这也正是在计算机上计算方阵/的行列式的常规方法。要将上述解方程步骤推广到”阶方程组，只需将对控制量“4”的操作改成对控制量力的操作即可。«元方程组高斯消元法的步骤如下：对于上=1,2,…,*7,约定=。步有*=4"一解"9）*岁）炉=母4-健、婕4））*6"k+l 44{讹/"依〃假定。依h°FORj:=k+lTOn{%：=%-X*%}〃ENDFORJ耳=*线}//ENDFORI}//ENDFORK「»、Xj=bi-3.FORi:=nTOI〃回代求解」{s:=biFORj:=i+lTOnDO5:=s一以歹*Xjxi:=s!ah}4.输出方程组的解々J=12…/。高斯消元法的运算量整个消元过程即式（3.14）的乘法和除法的运算量为职！要i/\M伍-2）2伽-左+D+-左）+—7~A-lJU1N回代过程即式（3.15）的乘除运算量为 45故高斯消元法的运算量为(3.16)高斯消元法的可行性在上面的消元法中，未知量是按照在方程组中的自然顺序消去的，也叫顺序消元法。在消元过程中假定对然元素“，消元步骤才能顺利进行，由于顺序消元不改变上的主子式值,故高斯消元法可行的充分必要条件为5的各阶主子式不为零。但是，实际上只要det'*0，方程组人入二小就有解。故高斯消元法本身具有局限性。另一方面，即使高斯消元法可行，如果।।很小，在运算中用它作为除法的分母(3.17)(3.18)(a,—a#-a诉会导致其它元素数量级的严重增长和舍入误差的扩散。这是高斯消元法的另一缺陷。…』何0003xi+3.0000%.=2.0001例3.1方程”0000%i+i.0000x2=1.0000=2=2的精确解为：/5'必3o在高斯消元法计算中取5位有效数字。解:方程(3.17)x(-1)/O.OOO3+方程(3.18)得：'0.0003再+3,0000x2=2.0001’9999.0与=6666.0=06667,代入方程(3.17)得々°。由此得到的解完全失真，如果交换两个方程的顺序,得到等价方程组'l.OOOOxj+1.0000x2=1.0000"0.0003%,+3.0000x2=2.0001经高斯消元后有10000玉+1.0000X2=10000[2.9997勺=1.9998,,.xa=0.6667,再=0.3333 46由此可看到，在有些情况下，调换方程组的次序对方程组的解是有影响的，对消元法中抑制舍入误差的增长是有效的。如果不调换方程组的次序，取6位有效数字计算方程组的解，得到x2=0.666667,X]=0,33取9位有效数字计算方程组的解，得到&=0.666667,句=0.333333由此可见有效数字在数值计算中的作用。列主元消元法如果在一列中选取按模最大的元素，将其调到主干方程位置再做消元，则称为列主元消元法。调换方程组的次序是为了使运算中做分母量的绝对值尽量地大，减少舍入误差的影响。用列主元方法可以克服高斯消元法的额外限制，只要方程组有解，列主元消元法就能畅通无阻地顺利求解，同时又提高了解的精确度。maxM="J更具体地，第一步在第一列元素中选出绝对值最大的元素日公，交换第一行和第加行的所有元素，再做化简…为零的操作。对于每个上土=1,2,…避-1在做消元前，选出（I碓讣到…旧?】中绝对值最大的元素脸"L对上行和m行交换后，再做消元操作，这就是列主元消元法的操作步骤。由于detnHO,可证_（^-d旗’*1次•'派中至少有一个元素不为零，从理论上保证了列主元消元法的可行性。列主元消元法与高斯消元法相比，只增加了选列主元和交换两个方程组（即两行元素）的过程。列主元消元法算法1.输入：方程组阶数附，方程组系数矩阵力和常数向量项5。 47{〃选择FORu:=k+lTOnlFM>sTHEN{w-u;s-|aM|}FORv:=kTOn//交换第上行和第次行“：=—=a-。tf“：=EA,；=tIFORi:=k+lTOn{。=a/a&FORj:=k+lTOn{%:=%-£*aJ//ENDFORJ}//ENDFORI}//ENDFORK3.FORi:=nTO1//回代求解 484.输出方程组的解玉」=1，2,…,”如果对于第k步，从此行至力行和从分列至万列中选取按模最大的叫k界记卜T)卜max{W|),,於舟&,对第发行和第〃行交换，对第上列和第V列交换，这就是全主元消元法。在k列和第v列交换时，还要记录下v的序号，以便恢复未知量H和刘的位置。3.1.3高斯一若尔当(Gauss—Jordan)消元法高斯消元法将系数矩阵化为上三角矩阵，再进行回代求解；高斯一若尔当消元法是将系数矩阵化为对角矩阵，再进行求解，即对高斯消元法或列主元消元法得到的等价增广矩阵：X/X/\JZ阳洲点谭如碌溜O如理oOhlooO用第4行乘一03：“44加到第3行上，用第4行乘一°9/比‘加到第2行上，用第4行乘-a»1加到第1行上，得到alla12fl130年”0(2)谭0欧,00a330欧000谭娟用第3行乘-a矍/%蓼加到第2行上，用第3行乘-/3翼加到第1行上，再用第2行乘一对/娟,加到第1行上，得到an000及少0碗00酸)00噌0理000溜用」(3.19)为方便起见，仍记式(3.19)系数矩阵元素为,常数项元素为“。那么 49=4/他,i=用初等变换化系数矩阵为对角矩阵的方法称为高斯一若尔当消元法。从形式上看对角矩阵(3.15)比上三角矩阵(3.11)更为简单，易于计算%,但是将上三角矩阵(3.11)化为对角矩阵(3.15)的工作量略大于上三角矩阵回代的工作量。高斯一若尔当消元法求逆矩阵设上为非奇异矩阵，方程组=4的增广矩阵为。如果对c应用高斯-若尔当消元法化为(S,则川=丁。3156J的逆矩阵才、rl2/=2435例3.2用高斯-若尔当消元法求560450231001p10720115/32001/32/3101-2/31/3110-1/31->0010-1/20-5/227013/203/2-1001/21-1/20,1001-32,010-33-1=（4团）0012-10-32,/=-33-1解得：2-103.2直接三角分解法 50相当于用了三个初等矩阵左乘力和方。记仍以〃=4为例，在高斯消元法中，从对方程组进行初等变换的角度观察方程组系数矩阵的演变过程。第一次消元步骤将方程组（3.9）化为方程组（3.10）,41=2,3,4,容易验证,ana12ai3ai4由&4X=布,得到Wx=W,其中将方程组（3.10）化为方程组（3.11）,相当于用了两个初等矩阵左乘月和小。记心如=3,4,有’100100J_/42ooWio】0001000o,001001 51(100010-z420L甲。)=卢=[同理，将方程组(3.11)化为方程组(3.12),相当于用一个初等矩阵左乘/和占。完成了消元过程，有 52亦有所有消元步骤表示为：（4与左乘一系列下三角初等矩阵。容易验证，这些下三角矩阵的乘积7仍为下三角矩阵，并有T-1=伉四看厂=写1石1写1101*32001’4301001于是有A=TAA=T-XA这里T"仍为下三角矩阵，其对角元素为1,称为单位下三角矩阵，而N已是上三角矩阵。记K一左W则有A=LU结果表明若消元过程可行，可以将总分解为单位下三角矩阵£与上三角矩阵U的乘积。由此派生出解方程组的直接分解法。由高斯消元法得到启发，对幺消元的过程相当于将幺分解为一个上三角矩阵和一个下三角矩阵的过程。如果直接分解工得到上和人A=LU。这时方程力x=6化为=令以=丁，由功=小解出丁；再由解出x。这就是直接分解法。 53将方阵力分解为工=£0,当£是单位下三角矩阵，是上三角矩阵时，称为多利特尔（Doclittle）分解；当Z是下三角矩阵，U是单位上三角矩阵时，称为库朗（Courant）分解。分解的条件是若方阵工的各阶主子式不为零，则多利特尔分解或库朗分解是可行和惟一的。3.2.1多利特尔分解多利特尔分解步骤若金的各阶主子式不为零，/可分解为工=£〃，其中上为单位下三角矩阵，为上三角矩阵，即矩阵£和。共有/个未知元素，按照U的行元素£的列元素的顺序，对每个°。列出式（3.16）两边对应的元素关系式，一个关系式解出一个上或^的元素。下面是计算£和。的元素的步骤：（1）计算。的第一行元素"11，%?，…，un0要计算,1,则列出式（3.20）等号两边的第1行第1列元素的关系式:■<211==(100)r-1故%1=/1。一般地，由U的第一行元素的关系式得到 54“n=aijJ=1，2，-，力(2)计算£的第一列元素JA1，…，41要计算"l,则列出式(3.20)等号两边的第2行第1列元素的关系式:。21=5?”％=fel10r-10)故’21=町1/%1。•般地，由上的第1列元素的关系式aa=Sirurl10r-10)U110Aian得到Ai=沏®/=2,3,…,％(3)计算的第2行元素(4)计算上的第2列元素若已算出U的前上-1行，£的前发-1列的元素，则(5)计算一的第发行元素…"心由。的第发行元素的关系式: 55afy=£｝口气=ik2…io…。)”r-l00"U是上三角矩阵，列标)-行标上.xki-1a/£%%=2?卢厂+%r«lr-1r>l得到fc-1=k,k+\,--,n(3.21)r-1(6)计算Z的第k列元素4+以如2,…‘由上的第上列元素的关系式：%3刃4-田…kJ,0,…,。)”；r-1U0i；£是下三角矩阵，.•.行标标》2行标上。xkIf州=多认k=工g=+4/玳r-1r-1得到-2?»"雄/"战，i=上+1，儿+2,…/，-1J一直做到£的第抬-1列，u的第%行为止。用直接分解方法求解方程组所需要的计算量仍为3,和用高斯消元法的计算量基本相同。 56可以看到在分解中上的每个元素只在式(3.21)或(3.22)中做而且仅做一次贡献，如果需要节省空间，可将以及£的元素直接放在矩阵4相应元素的位置上。用直接分解法解方程=3,首先作出分解为=£^令醒，解方程组切=力。由于七是单位下三角矩阵，容易得到i-1必=4一2%",/=1,2,…小"(3.23)再解方程组”=、，其中/=必-2〃1弓一,2,1i2(3.24)对力进行^分解时，并不涉及常数项右。因此，若需要解具有相同系数的一系列线性方程组时，使用直接分解法可以达到事半功倍的效果。多利特尔直接分解算法1.输入；方程组阶数％,系数矩阵力和常数项力。2FORk:=lTOn{FORj：=kTOn〃计算〃的第先行元素U&：=%r-1 57FORi:=k+lTOn〃计算£的第上列元素}//完成^分解3,FORi:=1TOn//解方程组功=84.FORi:=1TOn//解方程组”=丁'*'.=乂-2%勺/%5.输出方程组的解x”，=12…/例3.2用多利特尔分解求解方程组:2々+町+勺=4勺+3x2+2x3=6X]+2x2+2x3=5「21解：u?21,311’32O'"11u12%3'°0“22u231八00“33」对分=1，有%=a^.j=1,2,3«n=2,“12=1必3=1Al=an/un,i=2,3/21=1/2=0.5,Z31=1/2=05对分=2,有"22=a22=3-0.5=2.5 58“23=。23一，2/13=2-0.5=1.5,32=（。32-13声12）/〃22=（2~0.5）/2.5=0.6对上=3,有u33=&3_'31"13一’32“23=0-6r100W211'：.LU=0.51002.51.50.50,611000.6j解功=々即•>1=4,乃=4,乃=0.6x3=1,x2=1,X]=13.2.2追赶法很多科学计算问题中，常常所要求解的方程组为三对角方程组 59其中(3.26)并且满足条件:同>同>0,14阳4|+匕|，46=0」=1,2,---,»-1的1>1%1>0称/为对角占优的三对角线矩阵。其特征是非零元素仅分布在对角线及对角线两侧的位置。在样条插值函数的M关系式中就出现过这类矩阵。事实上，许多连续问题经离散化后得到的线性方程组，其系数矩阵就是三对角或五对角形式的对角矩阵。利用矩阵直接分解法，求解具有三对角线系数矩阵的线性方程组十分简单而有效。现将/分解为下三角矩阵£,及单位上三角矩阵u的乘积。即4=其中1V11,U=v»-l1(3.27)利用矩阵乘法公式，比较两边元素，可得到%=碑=2,3,…，力$=bL4=+吗,i=2,3,…,九ci=叫"，382,3,・・•，力-1于是有<=/，»=瓦，%=八一4¥_i，i=2,3,・•，力vi-1 60由些可见将上分解为z及u,只需计算｛吗｝及｛4｝两组数，然后解37，计算公式为：乃=工/叫Ji=(Ji=2,3,---,«(329)再解以=丁则得4=",Xj==冬%-1,…21(3.30)整个求解过程是先由(3.28)及(3.29)求｛用，｛吊｝及｛旬，谢日2…,n是“追,，的过程，再由(324)求出｛'J,这时1F，…」是往回赶，故求解方程组(3.25)的整个过程称为追赶法。3.2.3对称矩阵的£次5分解对称正定矩阵也是很多物理问题产生的一类矩阵，正定矩阵的各阶主子式大于零。可以证明，若金正定，则存在下三角矩阵£,使幺=££',直接分解幺=££7的分解方法，称为平方根法。由于在平方根法中含有计算平方根的步骤，因此很少采用平方根法的分解形式。对于对称矩阵，常用的是幺=£刀£「分解。对工作多利特尔分解工=£0,即 61由/对称正定，可得%,令由月的对称性和分解的惟一性可证即4=£戊。1%1…J1J(3.31)£是对角元素为1的单位下三角矩阵。n(n-1)对矩阵金作多利特尔或库朗分解，共计算M个矩阵元素；对称矩阵的£。尸分解，只需计算2个元素，减少了近一半的工作量。借助于多利特尔或库朗分解计算公式，容易得到分解计算公式。设/有多利特尔分解形式:…%*4其中Ju六%rv 62在分解可套用多利特尔分解公式，只要计算下三角矩阵£和。的对角元素"e。计算中只需保存L=%)的元素，尸的,行/列的元素用£的4表示。由于对称正定矩阵的各阶主子式大于零，直接调用多利特尔或库朗分解公式可完分解计算，而不必借助于列主元的分解算法。对于上=1,2,…，力，有〜比-1〜反-1八=%=%.-2.「求=akk-＞斤d//y-1r-1(3.32)i〜〜1A=(依1m=四一＞业小)"八r-17Jt-1L=(做-^见4)"Z=k=l,…/7(3.33)由L(DlF)x-b，解方程组=b可分三步完成：(1)解方程组友=b淇中z=DEx。j-1打(3.34)(2)解方程组°=z,其中1y=£4。X=zjd”j=1,2,--«(3)解方程"矛工丁毛=乃一£%%.,i1J-4+1对称矩阵的LDI；分解算法 631.输入：方程组阶数附，系数矩阵/和常数项合。2.FORk:=lTOnM=a*_£d心r-1FORi:=k+lTOn3.略去解方程组步骤从矩阵分解角度看，直接分解法与消元法本质上没有多大区别，但实际计算时它们各有所长。一般来说，如果仅用单字长进行计算，列主元消元法具有运算量较少、精度高的优点，故是常用的方法。但是，为了提高精度往往采取单字长数双倍内积的办法（即作向量内积计算时，采用双倍加法，最终结果再舍入成单字长数），这时直接分解法能获得较高的精度。例3.4用3广分解求解方程组:必==1k=1时，l21=叼/d]--1解：=%1*=1上_之时分=%一片4=2，32=(“32-自4/1)，d?=-0.5k=3时，dz=a33-耳id】-g2d2=3r100]p、Z=_110,Z)=21-0.513由LDlFx=b,有Lzsbfz(4,-4.6)rDy=zfy=(4-2.2/rl1 64ifx=y,x=-123.3范数3.3.1向量范数在一维空间中，实轴上任意两点a，b距离用两点差的绝对值区划表示。绝对值是一种度量形式的定义。范数是对函数、向量和矩阵定义的一种度量形式。任何对象的范数值都是一个非负实数。使用范数可以测量两个函数、向量或矩阵之间的距离。向量范数是度量向量长度的一种定义形式。范数有多种定义形式，只要满足下面的三个条件即可定义为一个范数。同一向量，采用不同的范数定义，可得到不同的范数值。定义3.1对任一向量XeR",按照一个规则确定一个实数与它对应，记该实数记为口”"，若满足下面三个性质：⑴VXd,有因|20,当且仅当。时，区卜0（非负性）⑵YXw*,QC&,有1a刈=1。惘I（齐次性）（3）v,有庐+”平I+IPI（三角不等式）那么称该实数为向量X的范数。几个常用向量范数向量工=（内，电，…，X"）的4范数定义为区「住同『””+8其中，经常使用的是P=l，2,8三种向量范数。区卜力引响+厉|+…+kl2-1或写成区卜J（X'X） 65I矶飞澧佃D=max{|xi|，kJ.，k|）例3.5计算向量X=（L3，-5），尸=1,2,8的三种范数。IM=l+3+5=9区卜（P+3?+5?产=5.9161|^|L=max{U|-5j-5向量范数的等价性有限维线性空间*中任意向量范数的定义都是等价的。若舄（*）'0（*）是我"上两种不同的范数定义，则必存在0<加<〃<8,使ae义■均有加国（与4鸟（X”初&（X）m<雪冷4M（Xh0）或为㈤ 66（证明略）向量的极限有了向量范数的定义，也就有了度量向量距离的标准，即可定义向量的极限和收敛概念了。则称向量列v（l）”…是收敛的〃是某种向量范数）二称为该向量序列的极限。设炉"，…，X®，…为？上向量序列，若存在向量ae使［史由向量范数的等价知，向量序列是否收敛与选取哪种范数无关。=f（»）w（«）V_1今向量序列<1'2,m=l，2,…收敛的充分必要条件为其序列的每个分量收敛,hm4a即「9存在。*尹F,贝产（L就是向量序列阳.图？物，'*））的极限。1r+i阴％归2k-\〃极限向量。例3.6求向量序列U解：算出每个向量分量的极限后得在计算方法中，计算的向量序列都是数据序列，当口।小于给定精度时,3.3.2矩阵范数矩阵范数定义定义3.2如果矩阵'的某个非负实函数V）,记作同，满足条件：⑴⑷20,当且仅当工=o时，|即=0（非负性）⑵1MH司⑷“6*（齐次性） 67⑶对于任意两个阶数相同的矩阵a,b有|a+b||«M|+|B||（三角不等式性）（4）工出矩阵为同阶矩阵1网14M眄（相容性）则称从（£*=MI为矩阵范数。矩阵的算子范数常用的矩阵范数是矩阵的算子范数，可用向量范数定义：VW,记方阵工的范数为网，那么网H冏六骑W或回=潮阳（338）称为矩阵的算子范数或从属范数。这样定义的矩阵范数满足矩阵范数的所有性质外，还满足相容性:上为"阶矩阵，▼“二,恒有।阳Kl'NM（3.39）根据定义，对任一种从属范数有M"=l,即单位矩阵的范数是1。常用矩阵范数向量有三种常用范数，相对应的矩阵范数的三种形式为：凤=幽2|%|4i-1（月的行范数）（3.40）LJ-*z(上的列范数)(3.41)(4是的最大特征值)(上的2范数)(3.42)证明：既然矩阵的算子范数是？上满足=1向量范数口⑷1的上确界，那么，找到这个上确界也就找到了矩阵的范数。(1)任取殷*,团！则 68网1=±1(4卜2支气勺卜七支同同2-12-1J-1I2-1=2冬kI)kl«(蟹热股"卜蟹热।atMil«喘Eki即修-“i-1xx另一方面设极大值在先列达到，即黎取X=e=(O,O,…,0,1,0,…,0)丁，e除第、个分量为1外，其余分量均为0,于是有阿=||(%,%力…%反升卜与⑷飞患与kl由于|*=1,故xX⑷2支力|蟆之同2.1J2-1因此有14▼kii=1l=ILr观一)卜腮勺l«腮£同1#电黑(2)任取X,U卜,则川，」川即K«果寿力J--J-1另一方面设极大值在此行达到，取X="1的%物知・•尸g-y这里0,a2O^=Vu<0II^L=iM 69于是.I胤2力4J-1凤=龈£同故■■川(3)为「j为对称非负矩阵，具有非负特征值，并具有附个相互正交的单位特征向量。设的特征值为4242.242°,相应的特征向量为丛，其中M为相互正交的单位向量。设.X=XM+jU2+...+、M,,并且区卜1,即含为，则|“『-(AX.AX)-[^AX,X)-力刎2<21sx:吗2-12-1即对任意优卜1均有"猊,故ML«网取£=巧,则有>闻2=4故MIL2M于是IMIL=x/a如果A是对称矩阵，那么A，工=A4,设j的特征值是4,则有IM]'=max5^=max|^|还有一种与向量范数l*l2相容的矩阵范数，称为欧几里得（Euclid）范数或舒尔（Schur）榛用I'L表示，其定义为 7011J”（3.43）因为欧几里得范数易于计算，在实用中是一种十分有用的范数。但它不能从属于任何一种范数，因为Mb与向量范数的等价性质类似，矩阵范数之间也是等价的。f-l2、公37W.MLML网g例3.7J人分别有。解：M卜max{|-1|+3,2+7}=9ML=max{H|+2,3+7）=10工7的特征值为=60.19,4=281MIL=疯而=（1+4+9+49）%=相矩阵范数等价性定理对上的任两种范数此及IHL,存在常数q£02>o,使巧帆=网心/风矩阵特征值与范数关系若兀是矩阵上的特征值（即存在非零向量x使得：ax，x）,对任一算子范数有 71又⑷因（相容性）邓||冏|即矩阵特征值的模不大于矩阵的任一范数。谱半径若上有特征值砧，-4记。⑷=腮明则称0⑷为上的谱半径。有了谱半径的定义，矩阵的2范数可记为：帆邛⑹4）。谱半径与矩阵范数关系由矩阵谱半径定义，可得到矩阵范数的另一重要性质， 72在解方程组时，我们总是假定系数矩阵/和常数项3是准确的，而在实际问题中，系数矩阵/和常数项b往往是从前面的近似计算所得，元素的误差是不可避免的。这些误差会对方程组人了=右的解X产生多大的影响？矩阵的条件数将给出一种粗略的衡量尺度。定义3.3若/非奇异称L为/的条件数。其中I比表示矩阵的某种范数。用矩阵力及其逆矩阵力-1的范数的乘积表示矩阵的条件数，由于矩阵范数的定义不同，因而其条件数也不同，但是由于矩阵范数的等价性，故在不同范数下的条件数也是等价的。矩阵条件数的大小是衡量矩阵“坏”或“好”的标志。对于线性方程组力x=5,若系数矩阵有小扰动人这时方程组的解也有扰动，于是(A+5A)(x+Jx)=i。JA对Jx的影响可表示为：同‘%")万明I(3.44)若常数3有小扰动，其对Jx的影响表示为：(3.45)因此，C。混(A)大的矩阵称为，，坏矩阵，，或，，病态矩阵，，。对于C"d(工)大的矩阵，小的误差可能会引起解的失真。一般说来，若金的按模最大特征值与按模最小特征值之比值较大时，矩阵就会呈病态。特别地，当det工很小时，/总是病态的。例3.8方程组'0.2160演+0.1441X2=0.1440[12969x1+0.8648x2=08642方程组的解为 732-2对常数项8引入扰动仍=(0.00000001,-0.00000001)1,则解为rx1+3x1[0,991,+3x2-0.4870俨=(7.009,1.5130)]虽然孙艮小，但解已完全失真。det=-10-8履=-10T'0.8648-1.2969-0.1441'0.2161M|j2.1617,|4*=1.513x108故C。叽(/)=3.2706521x108,/的病态是显而易见的。.5病态方程组的解法如果系数矩阵工的条件数%附“(工)远大于1,则■=合为病态方程组，对病态方程组求解可采用以下措施：(1)采用高精度运算，减轻病态影响，例如用双倍字长运算；(2)用预处理方法改善力的条件数，即选择非奇异阵尸，使融尸力与月x=5等价，而取。的条件数比月改善，而求4夕=否=母的解％=Q"x,即x=Q*为原方程组的解，计算时可选择P，Q为对角阵或三角阵; 74(3)平衡方法，当兑中元素的数量级相差很大，可采用行平衡或列平衡的方法改善总的条件数。设，=(a#)eK*非奇异，计算齿=蟹|%|(，一1，2,…㈤，令’sj,于是求4=B等价于求=或急=如这时急的条件数可得到改善，这就是行均衡法。例3.9给定方程组人工=小为1叫卜111叼,1042A的条件数3d⑷0M0,,若用行均衡法可取0=diag(10Y,l),则平衡后的方程卷=5,为iir^i[io4]~[io*i'弧[2，[11C(')9f用三位有效数字的列主元消元法求解得x=(1.00,1.00)T。3.6程序不例程序3.1用列主元高斯消元法求解方程组:算法描述输入：方程组阶数n,矩阵a及列向量b。分解过程：对上=L2,…选择交换第k行和第W行方程对i=k=1,…，此 75j=k+l,--,n阳=%一(4/%)*阳对，"+L…/4=4-3/%)*包回代过程：对，=乩"-1，…［解得X程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIItlpurpose:列主元高斯消元法解Ax二b〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 76inputnvalue(dimofAx=b):");〃输入Ax=b的维数scanf(u%dM,&n);if(n>MAXN)|printf(4tTheinputnislargerthanMAXN,pleaseredefinetheMAXN. 77^^); 78return1;)if(n<=0){printf(4tPleaseinputanumberbetween1and%d. 79",MAX_N);return1;|〃输入Ax=b的A矩阵printf(4*Nowinputthematrixa(i,j),i,j=0,…,%d: 80”,n-1);for(i=0,i 81”,n-1)；for(i=0；ivn;i++)scanf("％If^&bfi]);for(i=0；i 82a[mi][j]=tmp;}for(j=i+l;j 83M);//输出for(i=0;i 84,x[ir，);return0;)//EndofFile计算实例用列主元高斯消元法水求解方程组：'2X1+X2+X3=4<再+3x2+2x3=6A+2叼+2x3=5程序输入输出 85Inputnvalue(dimofAx=b):3Nowinputthematrixa(i,j),i,j=0,2:211132122Nowinputthematrixb(i),i=0,...»2:465Solve...x_i=1.0000001.0000001.000000程序3.2用库朗分解公式求解方程组: 86算法描述输入方程组阶数万，矩阵上及列向量占；将矩阵力分解为£沙，其中%T*21*22L=...Jul蜃**对上=1,2「一，力，对i=匕分+1,…/kk一，此一’"%息r-1对J=k+1,上=2,…/jt-i%=（%-E，b%）"心r-1对i—1,2,…,n必对i=n,n-1,…,2,1外=X-%E/Xj.J-2+1程序源码llllllllllllllllllllllllllllllll//Purpose:库朗分解解Ax=b//lllllllllllllllllllllllllllllll#include 87#include 88");return1;if(n<=0){printf(44Pleaseinputanumberbetween1and%d. 89,\MAX_N);retun1;〃输入Ax=b的力矩阵prinlf("Nowinputthematrixa(i,j),i,j=O...,%d. 90, 91-l);for(i=0;i 92, 93-l);for(i=0;i 94forti=0;i 95for(i=0;i 96M);〃输出for(i=0;i<=n;I++)printf(**%f 97”,x[i]);return0;}计算实例用库朗分解公式求解方程组：程序输入输出Inputnvalue(dimofAx=b):3Nowinputthematrixa(i,j)i,j=O,...»2：121-2-1-50-16Nowinputthematrixb(i),i=0,...»2： 9824-6350solve...x-i=7.0000004.0000009.000000程序3.3用分解求解方程组:算法描述输入矩阵4及列向量3.将矩阵/分解为兑=其中i=k+1,…，%7「Jl-1\k=%->44或/葭iJ解工x=B,先解£Z=s：再由a=z得乃=zJ4，i=最后由尸得 99x「yrZ%xki2,1程序源码iiiiiiiiiiiiiiiiiiiiiiiiiiiiiii〃Purpose：LDlF分解解方程组〃lllllllllllllllllllllllllllllll#include 100inputnvalue(dimofAx=b):M)；//输入Ax=b的维数scanf("％d”,&n)：If(n>MAX-N)pritf(4TheinputnislagerthenMAX-N,pleafineredefinetheMAX-N. 101M);return1;}if(n<=0){printf(44Pleaseinputanumbrbetween1and%d. 102,\MAX-N);retun1) 103//输入Ax=b的A矩阵printf(**Nowinputthematrixa(i,j),i,j=O,%d: 104”,n-1);for(i=0;i 105for(i=0;i 106M)〃输出for(i=0;i 107,\x|i|);return0;计算实例用尸分解求解方程组：程序输入输出Inputnvalue(dimofAx=6):3 108Nowinputthematrixa(i,j),i,j,=O,...，2:1-11-13-21-24.5Nowinputthematrixb(i),i=0,...,2:4-812Solve...x-i=1.000000-1.0000002.000000本章练习计算证明题1.分别用高斯消元法和列主元方法解方程(计算中取四位有效数字)：r0.002x1+87.13x2=87,15⑴[4.453^-7.26%2=37.27'001再-69.47x2=738.93(2)(2.01々+8.51叼=15.012.用高斯消元法编制程序解下列方程组：7.2%1+2.3x3-4,4x3+0,5x4=15.11.3xi+6.3x3-3.5x3+2.8x4=1.85.6x1+0.9勺+8.1%-1.3x,=16,61.5xj+0,4x2+3.7弓+5,9x4=36.93.用高斯消元法计算下列行列式的值:|-132A=21-236210-2-15=-210-1 109—1—25(2)第4章解线性方程组的迭代法用迭代法求解线性方程组工x='与第4章非线性方程求根的方法相似，对方程组M=y进行等价变换，构造同解方程组工=血"+8(对工乂=歹可构造各种等价方程组，如分解尸，"可逆，则由■=y得到公W'PX+旷、=MT+g),以此构造迭代关系式(41)任取初始向量=卜"耨)，…，x，)y，代入迭代式中，经计算得到迭代序列X”，*⑵，…。若迭代序列｛*("0｝收敛，设X®的极限为X、对迭代式两边取极限limX(*lim(W+g)(T9即x=mx+g,x是方程组工乂=丁的解，此时称迭代法收敛，否则称迭代法发散。我们将看到，不同于非线性方程的迭代方法，解线性方程组的迭代收敛与否完全决定于迭代矩阵的性质，与迭代初始值的选取无关。迭代法的优点是占有存储空间少，程序实现简单，尤其适用于大型稀疏矩阵；不尽人意之处是要面对判断迭代是否收敛和收敛速度的问题。。(河)=max|4|《i1可以证明迭代矩阵的与谱半径由公是迭代收敛的充分必要条件，其中乙是矩阵乱的特征根。事实上，若x*为方程组/乂=丁的解，则有x*rMX+g再由迭代式文=〃淤)+g，可得到X'- 110=V-炉D）=-那））limM”=o/»八/1由线性代数定理，…的充分必要条件°（故）<］。因此对迭代法（4.1）的收敛性有以下两个定理成立。定理4.1迭代法於+"=必"+g收敛的充要条件是十一0（—8）。定理4.2迭代法必""=g"+S收敛的充要条件是迭代矩阵M的谱半径。（⑼<1因此，称谱半径小于1的矩阵为收敛矩阵。计算矩阵的谱半径，需要求解矩阵的特征值才能得到，通常这是较为繁重的工作。但是可以通过计算矩阵的范数等方法简化判断收敛的工作。前面已经提到过，若।⑷ip矩阵上的范数，则总有“小2㈤。因此，若则乂必为收敛矩阵。计算矩阵的1范数和8范数的方法比较简单，其中HL=maxEkl1 1114.1.1雅可比迭代格式雅可比迭代计算«元线性方程组，如为+*+…+f（a2ix1+a23x2+-+a2],j,xj,=^2(4.2)，。杜玉+白也町+一+白皿々=居写成矩阵形式为*X=y。若的H°/=L2,…将式（4.2）中每个方程的他为留在方程左边，其余各项移到方程右边；方程两边除以内则得到下列同解方程组:记外=-%，％,&=yj%,构造迭代形式铲'摩+…+%炉+&,铲=坛承+J岁+…+%看）+g2,看叫=4看+1母）+…淄+g, 112迭代计算式（4.3）称为简单迭代或雅可比迭代。任取初始向量X'°）,由式（4.3）可得到迭代向量序列⑶巧雅可比迭代矩阵设力=L+D+U由■=»+力-方/=乙得到等价方程:DX=(D-A)X+y=Da(D-A)Xw+D~xy记B=D-\D-A)=l-£TlA,g=D-xy不难看出，B正是迭代式（4.3）的迭代矩阵，g是常数项向量。于是式（4.3）可写成矩阵形式:如D=B相Q+g其中：雅可比迭代算法下面描述解线性方程组"的雅可比迭代算法，为了简单起见，在算法中假定矩阵/满足雅可比迭代要求，即％=°，=1，2,…》，并设由系数矩阵/构造迭代矩阵8是收敛的。 1131.定义和输入系数矩阵/与常数项向量丫的元素。2.FORi：=192f.・.，n伊：=川%〃假定aijO，，72,…,”形成常数项向量FORj:=l,2,...,n%:=-'I/%=。]〃形成迭代矩阵元素3xl：={0,0,---,0)r,x2:={1,L…」产//赋初始值，xl和x2分别表示炉"和炉M4.WHILE卜1-叫1~>10司x1:=x2x2:=B*xl+g//FORu:=l,2,...,n//$：=g[u];//FORv:=l,2,...,n5；=s+Z>[u][v]*xl[v];//x2[w]:=5；ENDWHILE5.输出方程组的解x2”1=L2,例4.1用雅可比方法解下列方程组: 1142X]-x2-x3=-5 115代矩阵的谱半径也益时，迭代收敛，这是收敛的充分必要条件。迭代矩阵的某范数归八’时,迭代收敛。要注意的是范数小于1只是判断迭代矩阵收敛的充分条件，当迭代矩阵的一种范数1皿1>1,并不能确定迭代矩阵是收敛还是发散。例如，I02°用，贝"矶=1。悯1=1」，但它的特征值是0.9和0.8。°⑻ 116而0(1-£)%)”(1-h/),得0⑻=疝-£)切〈1由系数矩阵/构造的雅可比迭代矩阵收敛。'310"1-41_11-3一(如矩阵L」既是行对角占优阵，也是列对角占优阵)定理4.4若方程组系数矩阵兑为对称正定阵，并且20-力也为对称正定，则雅可比迭代收敛。4.2高斯-塞德尔(Gauss-Seidel)迭代法高斯-赛德尔迭代计算在雅可比迭代中，用卜,)，患)，…，的值代入方程(4.2)中计算出卷吗］=12…，起的值，xf"】)的计算公式是铲】)=£%摩+&J-1(JU1)(JU1)(JI+1)事实上，在计算不前，已经得到再再-1的值，不妨将已算出的分量直接代入迭代式中，及时使用最新计算出的分量值。因此玉的计算公式可改为：X尸尸)+支乐炉+&j-1;-i+l即用向量)计算出再的值，用向量Ml，X—,x*/计算出叼的值……，fx(M...X(Mx(e).•/勒『(M)用向量I1''i，।，"/计算出演的值，这种迭代格式称为高斯―塞德尔迭代。对于方程组4X=y,如果由它构造高斯-塞德尔迭代和雅可比迭代都收敛，那么，多数情况下高斯一塞德尔迭代比雅可比迭代的收敛效果要好，但是情况并非总是如此。 117构造方程组及的高斯廛德尔迭代格式的步骤与雅可比类似，设■产°,=L2,…,名将式（41）中每个方程的句号留在方程的左边，其余各项都移到方程的右边；方程两边除以“”，得到下列同解方程组：X]——A,—十anauanx2=-■再竺！馥+422222___%♦,一一上打X*———X]——叼…X*+—amtaxxa*.axx记为=-%/%，&=必/传，对方程组对角线以上的飞取第上步迭代的数值，对角线以下的不取第k+1步迭代的数值，构造高斯一塞德尔迭代形式：(4.5)例4.2用高斯-塞德尔方法解方程组：2xx-x2-x3=-5 118--0.2-(-0.93)+0.2-2.1+1.6=1.994x?)=-01-(-0.93)-0.111.994+1.1=0.9936计算结果如表4.2所示。表4.2计算结果k必兀3L00001-2.52.11.142.52一0.882.0040.98761.623-1.00421.99841.00060.12424-1.00052.00021.00000.0037高斯一塞德尔迭代矩阵设/=L+D+U写成等价矩阵表达式：AX=(D+L+U)X=(D+E)X+UX=y(D+L)X=-UX+y构造迭代形式：(D+£)月("1)=-UX8+y有1(川)一(。+£尸以⑶+»+£)+(46)则高斯-塞德尔迭代式(4.4)为 119炉川)=弦㈤+/(47)记G=-(。+Z)T。J=(Z)+£)>G称为高斯一塞德尔迭代矩阵。高斯-塞德尔迭代算法高斯一塞德尔迭代的程序实现与雅可比迭代步骤大致相同，对于方程组"=',在前面的雅可比算法中，假定雅可比迭代矩阵为B，xl表示M2X2表示x<*+]),其迭代核心部分是x2=B♦力+g。下面的算法给出由8和灯计算工2的过程，省略了形成迭代矩阵B和对xLx2初始化的部分。雅可比迭代的核心部分：while3-*2|L>1°"for(u:=1;u<=n,u++)xl[u]:=x2[u]for(i:=l;j<=n;j++){s:=g[i];for(j:=l;i<=n;i++){s:=s+b[i][力x2[j]}//注意x2|j]%2[i]:=sENDWHILE在高斯-赛德尔迭代计算中并不需要形成迭代矩阵G,由式(4.5)可知在计算中只要形成矩阵5=(%)=(-，/%)在程序中令向量=…x(ul)冽…ri=(xw/)…孙r工乙⑺51Mi5，毛，，小「它的核心部分是计算迭代式X2-B*X2+g计算中只需及时将x产”放到再⑸的位置上。高斯-塞德尔迭代的核心部分： 120WHILEllxl-z2L>1qH)for(u:=1;u<=n;u++)xl[u]:=x2[u]for(i:=l;j<=n;j++){s:=g[i];for(j:=l;j<=n;j4-+){s:=s+b[i][刃x2[j]}〃注意x2|j]x2[i]:=s}ENDWHILE上列算法是在假定迭代收敛的前提下，使用当型(WHILE)结构控制循环。更一般地，可将上列算法中WHILE循环改为FOR循环，通过控制循环次数和观测计算误差终止循环。届将上列算法中WHILE语句改为WHILE卜—ML〉：!。"以循环次数〈京这时在程序中要增加循环变量的设定和运算。判断高斯塞德尔迭代收敛的方法与判断雅可比迭代收敛类似，一方面从高斯-塞德尔迭代矩阵G获取信息，当0(G)<1或G的某种范数归Ik1时，迭代收敛：另一方面，直接根据方程组系数矩阵的特点作出判断。定理4.5若方程组系数矩阵4为列或行对角优时，则高斯塞德尔迭代收敛。定理4.6若方程组系数矩阵A为对称正定阵，则高斯塞德尔迭代收敛。例4.3方程组乂工=占中,证明当-U2〈a 121而△?=detA=1+2a3-3a2=(l-a)2(l+2a)>0得>-1/2.于是得到7/2 122X(K^=(1-(D)XW+a}(LX^+UXW+g)(48)这里田为修正因子，称为松弛因子，而式(4.8)称为松弛迭代。迭代的分量形式为'铲1;(1-°)染+或如蜡+如谈+…+'守+号)铲)=(1-0?)康)+峻1铲"+%镇)+…+%炉+g2)看入。-。对)+奴％白泮+&后9+…+1婿”+4)(4.9)式(4.9)称为松弛迭代的计算格式。例4.4给定方程组 123用SOR法求解，取卬=145解：用SOR迭代公式可得‘看川=（1-。）染+（0+2蜡+烧），刀用）=。_助老）+£「2+2看叫+X，）¥"=（1-。）#）+巴（3+x严+x严）取初始值:即：。1,球如果用高斯-赛德尔迭代法（0=1）迭代72次得：M科=39999952,0.9999945,1.999995）］用SOR迭代法3=145）,只须迭代25次即得：》到=（09999994,0.9999998,2.000000）T逐次超松弛迭代算法下列算法假定迭代矩阵收敛，否则要在WHILE循环中增加判断条件。1.定义和输入系数矩阵上与常数项向量y的元素，输入松弛因子。的值。2.FORi:=l,2,...,n回=%/%H假定,产0,=1,2,…,力,形成常数项向量forJ：=L2,…/3-%；阳=0I3xl:=（0,0-,0）r,x2：=（U--,l）r4.WHILEim-x2iL>ioYFOR(u=l,w<=nfu++) 124xl[w]:=x2[w],FOR(i=l;i<=nj++){s：=g「].FOR(J=\\j<-»J++)6=s+飒U]*x25}x2[i]:=(1-9ms)ENDWHILE5.输出吗'=12…,”松弛迭代矩阵将式(4.9)中含有太“"与》八的项分别放在方程两边:(1-OjL)X(M)=((1-3)1+Q}U)XW+QJg用八2乜力=2也代入上式，有(1+=((Z-(d)I-©Z尸Z7)X⑻+eug**】)=(/+皿匕)-1^-(d)I-(dD-xU)Xw+卬(/+£»。-1£)-也 125则松弛因子为田的迭代矩阵为W=S.Xg+/其中邑=(1+曲乜尸[(1-0})I-ojD-1U]/=oj（z+的乜严定理4.7逐次超松弛迭代收敛的必要条件0<©<2。定理4.8若上为正定矩阵，则当时，逐次超松弛迭代恒收敛。以上定理给出了逐次超松弛迭代因子的范围。对于每个给定的方程组，确定）究竟取多少为最佳，这是比较困难的问题，对某些特定的方程组，我们可以得到一些理论结果。通常，把©<1的迭代称为亚松弛迭代，把卬=1的迭代称为高斯-塞德尔迭代，而把的迭代称为松弛迭代。4.4逆矩阵计算在线性代数中逆矩阵是按其伴随矩阵定义的，若出产°则方阵人可逆，且,其中4为金的伴随矩阵。要计算笳个%-1阶的列式才能得到一个伴随矩阵，在数值计算中因其计算工作量大而不被采用。通常对（A，）做行的初等的效换，在将月化成/的过程中得到（"“）。在数值计算中，这仍然是一种行之有效的方法。由逆矩阵的定义4小7,令工=1=（脸,有化为万个方程组 126为是第J个分量为1,其余分量为0的〃维向量。或记为：'X'=',/=1,2,…储。用直接法或迭代法算出勺'」=L2…'乩也就完成了逆矩阵a"计算。如果依次对（4%），5品）,…,（'与）用高斯若尔当消元法，组合起来看有（当然也能组合起来做）：…,%）=（41）这正是在线性代数中用初等变换计算逆矩阵的方法。由此可见，计算一个万阶逆矩阵的工作量相当于解N个线性方程组。在数值计算中常常将计算矩阵逆的问题转化为解线性方程组的问题。例如，已知方阵总和向量》八有迭代关系式=对，在计算中不是先算出上"，再作a』与X⑸长的乘积得到如T;ffiB上作为线性方程组系数矩阵，求解方程组作为常驻数项解出玄叫4.4程序示例程序4.1用雅可比迭代求解方程组：算法描述输入系数矩阵4及常数项向量c。按雅可比迭代公式: 127‘铲"=(-如杯—%卓+田/%铲”=(F国对。2*炉+七)他2必")=(-%x，对4**1看1+G)/X、KR1it-XX，lift求解及=c。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII〃Purpose:雅可比迭代求解线性方程组〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 128lnputnvalue(dimofAx=c):〃);〃输入方程的维数scanf(ft%d〃,&n);if(n>MAX-N)printf("TheinputnislargerthenMAX-N,pleaseredefinetheMAX-N. 129〃);) 130if(n<=0){printf(〃pleaseinputanumberbetween1and%d. 131),f,MAX-N};retum1;}〃输入Ax=c的系数矩阵A.PRINTF(〃Nowinputthematrixa(I,j,)=O,...,%d: 132,f,n-l);for(i=0;i 134",MAXREPT);/输出return1;{计算实例用雅可比方法解下列方程组：64演-3x2-=14«2占-90x2+=-5西+町+40x3=20程序输入输出Inprtnvalue(dimofAx=c):3Nowinputthematrixa(j,j),i,j=0,...2:64-13-12-9011140Nowinputthematrixc(i),i=0,...,2: 135Solve...xi=0.2295470.0661300.492608程序4.2用逐次超松弛迭代（。作参数）求解方程组:输入矩阵工及列向量以按因子为。的超松弛迭代公式：'看。-。）染）-。3标/+-+a3*<#1）=（1-⑼章-承+与耳对…+%*看）f看+】）=（1-⑼看）-卬（％染）+…+。3出”工）小求解Ex=C。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII〃Purpose:超松弛迭代求解线性方程组〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 136#defineepsilon0.00001〃求解精度intmain(){intn;inti,j,k;doubleerr,w;staticdoublea[MAX_N][MAX_N],b[MAX_N][MAX_N],c[MAX_N],c[MAX_N],g[MAX_N];staticdoublea[MAX_N],nx[MAX_N];PRINTF(tfc 137lnputnvalue(dimofAx=c):");〃输入方程的维数ScanfC4%d”,&n);If(n>MAX_N){printf(kTheinputn-islargerthenMAX_N,pleaseredefinetheMAX_N. 138’‘)；return1;)if(n<=0){printf(44Pleaseinputanumberbetween1and%d. 139M,MAX_N);return1；}〃输入a=’的人矩阵printf(44Nowinputthematrixa(i,j),i,j=0,...%d: 140"，n-1);for(i=0;i 141”,n-1);for(i=0;i 142if(w 143’‘)；return1;(forti=0;i 144if(err 145");〃输出for(i=0;i 146,,,x[i]);return0;}(printf("After%drepeat,noresult... 147",MAXREPT);〃输出return1;计算实例解下列方程组：64演-3x2-=14«2占-90x2+=-5西+町+40x3=20程序输入输出Inputnvalue(dimofAx=c):3Nowinputthematrixa(i,j),i,j=0,...,2：64-3-12-9011140Nowinputthematrixc(i),i=0,...,2：14-520Nowinputthewvalue:1 1480.2295470.0661300.492608本章练习计算证明题1,计算下列矩阵的U八111，IIA|b，||A||8，三种范数:(511)00130162.计算下列矩阵的谱半径：r522、(31)3=260B=I⑴I131⑵1204)3.已知方程组：(1)写出解方程的雅可比迭代计算式,并以X(0)=(0,0,0,0)T为初值，迭代计算⑶。(2)写出解方程组的高斯-塞德尔迭代计算式，并以为初值，迭代计算 1492.以为初值，用雅可比迭代求解下列方程组的。2%!一叼+工3=73工］+3x2+9x3=03再+3x2+5的=45再一盯+弓=-13再+6x2+2x3=0x1-x2+2x3=43.用高斯・塞德尔迭代求解:lOxj_2/一弓=0-2再+10x2-弓=-21一再-2x2+5x3=-20因川)-X叫〈演自取初始值，当IL时迭代停止。5再-x2-x3=163再+60+2弓=11百一叼+2x3=-2闻||<10-3为初值，当“卜时迭代停止。4.已知方程组"1+2马f。(1)写出解方程组的雅可比迭代的迭代矩阵，并讨论迭代收敛条件。(2)写出解方程组的高斯-塞德尔迭代的迭代矩阵，并讨论迭代收敛条件。5.设有系数矩阵：'12-2)(2-1-/=1115=111证明(1)系数矩阵力对雅可比迭代收敛，而高斯-塞德尔迭代不收敛。(2)系数矩阵B对雅可比迭代不收敛，而高斯-塞德尔迭代不收敛。 1508*.选用一种直接方法或一种迭代方法，在计算机上编制程序计算下列矩阵的逆:r12or-10232315017333-21253-4B=-326112-135.1引言在实际问题中，有时只能给出函数,（x）在平面上的一些离散点的值｛（4/（玉）））,=°，L…*,而不能给出了（x）的具体解析表达式，或者/（X）的表达式过于复杂而难于运算。这时我们需要用近似函数P（x）来逼近函数/（X）,在数学上常用的函数逼近的方法有:•插值。•一致逼近。•均方逼近或称最小乘法。什么是插值？简单地说，用给定的未知函数/（X）的若干函数值的点构造了（X）的近似函数0。）要求P（X）与/（X）在给定点的函数值相等，则称函数伊（X）为插值函数。例如：在服装店订做风衣时，选择好风衣的样式后，服装师量出并记下你的胸围、衣长和袖长等几个尺寸，这几个尺寸就是风衣函数的插值点数值，在衣料上画出的裁剪线就是服装师构造的插值函数°（x）,裁剪水平的差别就在于量准插值点和构造合乎身材的插值函数。定义5.1,a）为定义在区间上的函数，，。，孙…，勺，为h⑸上力+i个互不相同的点，”为给定的某一函数类。若"上有函数P（x）,满足唉）=」（初i=0,l，…/则称伊（X）为了（X）关于节点X"在"上的插值函数；称点为，々，…，X”为插值节点：称｛（七，/（々））｝,=0,1,…为插值型值点，简称型值点或插点；/（X）称为被插函数。这样，对函数/（x）在区间卜”］上的各种计算，就用对插值函数0。）的计算取而代之。构造插值函数需要关心下列问题：・插值函数是否存在？・插值函数是否唯一？ 151・如何表示插值函数？如何估计被插函数与插值函数的误差？5.2拉格朗日（Lagrange）插值可对插值函数程。）选择多种不同的函数类型，由于代数多项式具有简单和一些良好的特性，例如,多项式是无穷光滑的，容易计算它的导数和积分，故常选用代数多项式作为插值函数。5.2.1线性插值问题5.1给定两个插值点（瓦沙。），（々J1）其中/x々，怎样做通过这两点的一次插值函数？过两点作一条直线，这条直线就是通过这两点的•次多项式插值函数，简称线性插值。如图5.1所示。图5.1线性插值函数4S）在初等数学中，可用两点式、点斜式或截距式构造通过两点的一条直线。下面先用待定系数法构造插值直线。设直线方程为=/+%*，将（/，M）），（/，为）分别代入直线方程A。）得: 152当X。时，因H0,所以方程组有解，而且解是唯一的。这也表明，平面上两个点，有且仅有一条直线通过。用待定系数法构造插值多项式的方法简单直观，容易看到解的存在性和惟一性，但要解一个方程组才能得到插值函数的系数，因工作量较大和不便向高阶推广，故这种构造方法通常不宜采用。当为=々时，若用两点式表示这条直线，则有:T/、X一再X一/WA——L^o+——勺一再再一而(5.1)这种形式称为拉格朗日插值多项式。XFoX)-Xo,4(x)虱，）称为插值基函数，计算'o（x）,氧幻的值，易见4(勺)=%(5.2)在拉格朗日插值多项式中可将A（x）看做两条直线瓦-不点的作用和地位都是平等的。演一号的叠加，并可看到两个插值拉格朗日插值多项式型式免除了解方程组的计算，易于向高次插值多项式型式推广。线性插值误差定理5.1记A。）为以（飞）。）＜々）1）为插值点的插值函数，这里%=/区）,%=/&）,设/（x）一阶连续可导，/"（x）在（。⑼上存在，则对任意给定的xe【a，句至少存在一点,使K(x)=/(*)-4(x)=:-^(x-x0)(x-x1),^e[a,b](5.3) 153证明令&⑶=/（x）-*x）,因氏（勺）=&⑸=0,%,々是&（X）的根，所以可设R（x）=上（力。-砧。-々）对任何一个固定的点X,引进辅助函数¥«）：¥«）=/«）-4。）~灰x）Q-xJQ-公）则平（为）=0，i=0,1。由定义可得+。）=°,这样+⑥至少有3个零点，不失一般性，假定为＜X"分别在［%,x］和〔冗'J上应用洛尔定理，可知+'⑴在每个区间至少存在一个零点，不妨记为县和易，即+'（«）=°和%©）=0,对+'©在［£易］上应用洛尔定理，得到P⑷在［£易］上至少有一个零点7,*©=00现在对+‘©求二次导数，其中4«）=°«i⑷为：的线性函数），故有+3=/"«）-2Mx）代入打得广©-2!稔）=0所以的此广©/2!R（x）=A^（x-Xo）（x-x）Je［a,»］即2!5.2.1二次插值问题5.2给定三个插值点（“"&））,，=°，1，2,其中不互不相等，怎样构造函数”X）的二次的（抛物线）插值多项式？平面上的三个点能确定一条次曲线，如图5.2所示。 154it式x)图5.2三个插值点的二次插值仿造线性插值的拉格朗日插值，即用插值基函数的方法构造插值多项式。设%(x)=4)+J1(x)/(^)+/2(x)/(x2),每个基函数4(x)是一个二次函数，对4(x)来说，要求修产2是它的零点，因此可设，o(x)=j(x-xj(x-x2)同理X"),4(x)也相对应的形式，得Z2(x)-^(x-x1)(x-x2)/(Ab)+5(X-Ab)(z-Xa)/(X1)+C(x-Ao)(x-x1)/(x2)将x=&代入&(X),得4(%)=4(4-石)(%-必)/(%)=〃％)A---(XO-X1)(XO-X2)"(x)=<(X-X,(X-X2)_("々)(.一盯)(x0-x1)(x0-x2) 155同理将X=X1，X="2代入乙u）得到B和C的值，以及4（x）和L（x）的表达式。（々-%）（公一刍）]（X2-XO）（X2-X1）AW-，2（X）（x-才0）（钎巧）卜「％）（々-电）（X一飞）（钎占）5一%）（々-々）.「仆_（xr）（x-X2）""尸许/F+（X"。）。-"（再一而）（再七2）一（X-两）（x-xj（勺-勺）（电一才1）NO/'）2-0fM〃西）（5.4）也容易验证：，0（%）=1，4（^1）=0.，0（勺）=04（%）=。，-=1，4（々）=04（X0）=0,，2（演）=0,，2（勺）=1插值基函数仍然满足:二次插值函数误差:鸟（X）=（X-X。）（X-再）（X-电），Je[min（x0,x1,x2,x},max（Ai）,x1,x2,x}]上式证明完全类似于线性插值误差的证明，故省略。 156插值作为函数逼近方法，常用来作函数的近似计算。当计算点落在插值点区间之内时叫做内插，否则叫做外插。内插的效果一般优于外插。例5.1给定应111'=0」90809,sinl20=0.207912。构造线性插值函数并用插值函数计算sinirs。'和sinl0030'解：构造线性插值函数：4⑺=--12）0.190809+）二0.20791211-1212-11分别将"=115,*=10・5代入上式，得kQi5)=0-199361准确值疝「11噂0'=0.1993684。05)=0.182258准确值Gn]0。30'=0.182236五（力=^^6"。）（工一4）=学一1）（12）=0.125例5.2给定sell0.190809snl20=0.207912,sui13'=0.224951。构造二次插值函数并计算sinmo'o解:4(*)=0.190809(x-12)(%-13)(11-12)(11-13)(11)(23)(12-11)(12-13)0.2079120.224951 1574"5)=0.199369,准确值sinll°30'=0.199368例5.3要制做三角函数的函数smx值表，已知表值有四位小数，要求用线性插值引起的截断误差不超过表值的舍入误差，试决定其最大允许步长。解：设最大允许步长h=hi=xi-xi_l1,I?1ooZh<0.02523%次拉格朗日插值多项式问题5.3给定平面上两个互不相同的插值点（0/a）），i=0」，有且仅有一条通过这两点的直线；给定平面上三个互不相同的插值点=有且仅有一条通过这三个点的二次曲线；给定平面上？j+i个互不相同的插值点a，/a）），i=o，i，2,…3，互不相同是指引互不相等，是否有且仅有一条不高于〃次的插值多项式曲线，如果曲线存在，那么如何简单地作出这条〃次插值多项式曲线？分析:万次多项式月（Q=/+。7+'”+%/，它完全由丁+1个系数两决定。若曲线只。）通过给定平面上%+i个互不相同的插值点aJa））'=°J2,…述，则々（x）应满足片阮）=上）八0,1,2「/,事实上一个插值点就是一个插值条件。将a=0,1,2,…,％依次代入々（x）中得到线性方程组： 158+白2片+•••+%X；=/（xo）,4+%再+叼X；+…+%X：=/（公）旬+为、+%*；+…+a*x：=/（x„）方程组的系数行列式是范德蒙（Vandermonde）行列式:na-弓）0 1594(x)=£4(x再)i-0(5.6)那么4/）不高于"次且满足'区）=/（%））=0,12…,”故4（x）就是关于插值点％,如…，々的插值多项式，这种插值形式称为拉格朗日插值多项式。&（X））称为关于节点｛玉｝的拉格朗日基函数。例5.4给出下列插值节点数据，做三次拉格朗日插值多项式，并计算,（0.6）»-2.000.001.002.00“X。17.001.002.0017.00解：拉格朗日插值基函数为:=上再乂钎动卜⑹°(%一再)(勺-电)(勺一弓)(x-0)(x-1,00)(x-2.00)"(-200-0)(-2.00-1.00)(-200-2.00)(x-O)(x-1,00)(x-2.00)"(-2.00-0)(-2.00-1.00)(-2.00-2.00)=一/57)(钎2)(钎.)(「一幻(”勺)1(x1-x0)(x1-x3)(x1-x3)=*+2)(x-1)(x-2),Zs=(…拼钎为玄一动“(X2-XO)(X2-X1)(X2-X3)»-l(x+2)x(x-2) 160(x-Xo)(x-xJ(x-弓)3(弓-工0)63-/)(弓一电)4(x+2)x(x-1)o三次拉格朗日插值多项式：1714⑶=-^x(x-1)(x-2)+w(x+2)(x-1)(x-2)-|(x+2)x(x-2)+?x+2)Q7)/(0.6)=4(x)=02560n次插值多项式的误差定理5.2设4㈤是［巴可上过aJ(xJ)，Xie【aj］，i=0,1，…*的抬次插值多项式，毛互不相等,当了eC^\a,b］时，则插值多项式的误差：43)=(、，。_/)(］一占)…(x-4),("1)!其中8［a,句(57)证明*:记名⑶=/(x)-4(x)。由于4出)=/(0i=0,1,…,”因而飞，如…,々是月⑶的根,于是可设均(x)=K(x)(x-%)(x-Xi)-(x-\)下面的目标是算出《(X),为此引入变量为£的函数¥")：中⑷=/W-4«)-K(x)("x())(£-X。…-xj(58)令£=包得¥(西)=0,J=0,l,-,«令£=x,由定义乎(不)=了⑶-4(，)-K(x)(x-x0)(x-X。…(x-4)=0即至少有％+2个零点，由于『wC"】［a,句，由洛尔定理，+'«)在+«)相邻的两个零点之间至少有一个零点，即+’(£) 161至少有"+i个零点。同理再对+“a）应用洛尔定理，即+”1）至少有附个零点，反复应用洛尔定理得到另一方面，对于（'）求%+1阶导数，有”+1)Q)=尸)0-K(x)5+1)10-3*+i)©=©-K(x)5+1)I得到（"1）!4(x)/叫©5+1)!(x_/)(x_Xi)…(x-X*)，e[a,b](5.9)若I,阳⑶|4正仅口k-小(5.10)由于¥（“叫（。的零点,与¥«）的零点与，电…，'有关，因而，为x的函数。F""(x)，Kxe[a,b],则4(x)可表示为由（5.9）式可以看到，当“X）是不高于附次的多项式时，4（x）=°,即4（x）=/（x）。对于函数/（x）h，，2=0」,….，关于节点X。，电…，X”的拉格朗日插值多项式就是其本身，故拉格朗日基函数4（x"满足4(X)=Zt(x)x；=x",k=0,1,一*2-0定理5.2给出了当被插函数充分光滑时的插值误差或称插值余项表达式，但是，在实际计算中，并不知道/（X）的具体表示，难以得到,的形式或较精确的界限河，因此也难以得到界在实际计算中，可对误差运用下面的事后估计方法。给出川+2个插值节点而，再，…，/+1,任选其中的力+1个插值节点，不妨取玉/=°，1，…构造一个万次插值多项式，记为43）。在力+2个插值节点中另选力+1个插值点，不妨取玉」=1，2,…,+1,构造一个万次插值多项式，记为 1624（“）。由定理2可得到(5.11)("1)1/（X）-4（X）=—_（X-演）（X-X?）…（X-%）（"1）!（5.12）设/"'（X）在插值区间内连续而且变化不大，有/叫刍A则」（x）-q（x）〜X-//⑶-4⑶X-%,+1从而可得到/㈤、工22±1_卬力+士^x。一、+1X»+1一%（5.13）“X）-4（x）“干2_&㈤-4㈤）XoF+1（5.14）拉格朗日插值多项式的算法下面用伪码描述拉格朗日插值多项式的算法。1：输入:插值节点控制数万，插值点序列a，乃“工…泮，要计算的函数点x,及变量抹="2：FORi：=0,1,n〃i控制拉格朗日基函数序列{tmp：=1；2.1FORj：=0,1,,i-1,i+1,…，n4。）=n--0ij<»X「Xj{〃对于给定X,计算拉格朗日基函数方tmp:={忤x（x-r）/（演-X,） 163）〃tmp表示拉格朗日基函数"（X）2.1力：=力+他xy}2.2出&（X）的计算结果亦。5.4埃特金逐步插值法我们已经知道，拉格朗日插值多项式的余项为&（♦）=f+1\;。一/）。一硝…（xrj,5+1）!其中句但在实际插值过程中，由于‘""（/一般无法计算，因而实际的误差也就无法得知。那么，如何根据所要求的精度来选取插值点的多少呢？这就是埃特金（Aitken）逐步插值所要解决的问题。首先介绍一个误差的事后估计法。设外-1,4X）为通过飞，孙…,点的分次插值多项式；外-i.川（X）为通过勺,々,…，私!/卬点的发次插值多项式。这两个上次插值多项式都通过通‘演'…''1,其中前者还通过乙，而后者还通过敌“，即它们所通过的插值线结点有上个是共同的，只有一个是不同的。若插值多项式用Aj（x）表示，则其中第一个下标i表示插值多项式通过；+1个插值结点x°'Xi'…'々一七；第二个下标,表示还通过的插值结点勺。 164现在考虑插值多项式At•式"与外-1.川（X）,由于它们所通过的插值结点中有一个不同的（敌”小，,因此其余项也不同。设它们的余项分别为(K\k川(x)=^2加f)X**1)(f!\M其中刍）与‘备）虽然不等，但由于在插值区间内这两个插值多项式所通过的插值结点大部分相同，只有一个不同，所以这两个分+1阶导数相差不多（假设费与易离得很近），可以近似看作相等，即产1）陷卜产D©）由此可以得到这两个插值多项式的余项比为〃X）-pyE（X）〜Xf•/（x）-Pnm（x）X-Xm从上式中解出了（X）得/（xAps（x）+''［pjuu（x）-Pe-u+1（x）］即〃X）-p川（力+X敌［p川（x）-味,N（x）］凝-X川LJ由上式可以看出，上次插值多项式Pet•式X）的误差可以由具有上个相同的结点与一个不同结点的两个上次插值多项式的差来估计。并且还可以看出，如果将这个误差加上，则可以得到更精确的,（'）的近似值。可以证明，这更精确的插值多项式为化+1次插值多项式，即 165P31(x)=Pm(x)+X醺[pt(x)-Pijui(x)]Xjif+iLJ由上式可以看出，两个k次插值多项式Pl•4x）与经过线性组合后，便得到k+1次插值多项式0瓦卬。）o这种方法称为埃特金逐步线性插值。它的优点在于可根据精度的要求逐步提高插值的阶,在插值过程中，由于都是线性运算，因而计算也方便。其计算格式如图5.3所示。xo/(xb)a凡】V叼V\\Poa\\\P\2\X3/将)\AP03—QPi3\\P23图5.3埃特金逐步线性插值的计算格式例5.8设函数J*t已知%=0.3(xx=0.4(x2=0.5弓=0.6y3=0.58813O=0.7乂=0.68122"_y0=0,29850=0,39646(乃=0.49311求“0・462）的近似值。解：根据埃特金逐步线性插值的计算格式，其计算结果如表5.2所示。表5.2Kxkf(xk)Po,MPi*(x)Pj,k(x) 16600.30.2985410.40.396460.45719520.50.493110.4561340.45653730.60.588130.4549000.4564840.45655740.70.681220.4535020.4564320.4565570.456557由表5.2可知，若取£=0.00001,则|P23（X）-P2.4（X）|《£所以有/（0.462）80.45656此时,73.4实际上已用不着计算了。若取£=0.001,则匹卜）-%（“<£所以有“0.462）皿56此时，以后的值也不必计算了。埃特金逐步线性插值的算法输入：插值结点（々'勺）（‘插值点£及精度要求输出：插值点£处的函数值严PROCEDUREGATKN{X,Y,n,t,£,F}Po-y。；eps<-1.0+e,i<-lWHILE（i 167FORj=lTOiDO尸<-%t+&-XjG（丹-i-尸）/St-否）Pi—F；塞s一3一Pi』」J+1）OUTPUTFRETURN在上述算法中，济・•分别表示/。=/（%），"」（£），…，乡」（£），…。当某相邻两个Pi-u0之差的绝对值小于£时，插值过程就停止。如果所有的插值结点都已用完，但还未满足精度要求，插值过程也停止。5.5*埃尔米特（Hermite）插值在构造插值时，如果不仅要求插值多项式节点的函数值与被插函数的函数值相同，还要求在节点处的插值函数与被插函数的一阶导数的值也相同，这样的插值称为埃尔米特插值或称密切插值。常用埃尔米特插值描述如下：设/（X）具有一阶连续导数，以及插值点不，,=°，1，…泮，若有至多为2%+1次的多项式函数巴22。）满足/+G）=/⑷用=/'（xji=0,1，…/则称“2x+l（X）为了（X）关于节点{4}皿的埃尔米特插值多项式。问题5.7：给定“X。）=%，/（々）=乃/（勺）=%/（々）%怎样构造给定两个节点的函数值和一阶导数值的埃尔米特插值多项式？分析：用4个条件，至多可确定三次多项式。设满足插值条件的三次埃尔米特插值多项式为：%（X）=a0+乎+4/2+守3 168将插值条件代入/（X）得到线性方程组：，旬+/勺+叼片+&君=30%+a1x1+a2Xi+a3xf=yxax+2a/o+3a3%=%%+2«洛+3%x；=w1因为方程组的系数行列式1/X；1々X1012x0012占所以方程组有解，可惟一解出旬，生，为，。3,即关于节点「，片的埃尔米特插值多项式存在惟一。类似于拉格朗日插值多项式样构造手法，也可通过插值基函数作出玛（X）O设&3CO=%⑶%+似力必+go⑶飞+gl（X）矽要修（力=&（乃必1=%可设&（%）=L4（%）=o,&（而）=0,&（勺）=o同理要月3（々）=似々）必=必有&G1）=0,似为）=14（硝=QgGi）=0要用（飞）=%有用（勺）=0,限而）=0a（%）=1保（勺）=0,要 169有用（再）=0,4（b）=O,go（xx）=0,或再）=1由上述要求，对为（乃来说，至多为三次多项式，”1是它的二重根，可设F_为（X）=（旬+如）X=（即+如）看（X）—。一再1（5.20）利用=（/+如濡（勺）=（即+加0）=1吊&）=%片（飞）+（旬+为而）4&）4（勺）=0由go(x())=O,g(Xi)=g'(Xi)=O,可设go(x)=a(x-xo)/：(x)由g：(xo)=l,算出a=1o所以go(x)=(x-x0)/q(x)同理gl(x)=(x-xj/；(x)(5.21)综上所述，得到以质，々为节点的埃尔米特插值为%(x)=%(x)/(x0)+4(力/(石)+SoW/Vo)+&(力/'(再) 170(5.22)容易证明，当时插值误差为R(x)==’『(，-丽J(xe[/加(5.23)如果要构造了(')关于用+1个节点=0'的2%+[米特插值多项式，手法与构造两个节点的方法类似。»»为2@)=工用。)“演)+£号(制/⑷02-0这里&(x)，&(x),z=°,L…,％分别为不高于2花+1次插值多项式，分别满足%[&(勺)=0(&3)=0及卜；(勺)=@由此可得到'11%(x)=—/；")、川国-'”(5.24) 171gKx)=(xfW(x)这里(段(x))为关于节点看"=0/，…*的拉格朗日基函数。容易证明，当/eC?叫a囱时，误差为:R(x)=/(x)-/n⑶二券1「力5止”工兵叫")例5.8给定/(T)=°J⑴=4JS)=2J'⑴=0,求埃尔米特插值多项式,并计算“05)。解：%⑸=&(x)・0+4(x)，4+go(x)・2+gi(xAO显然本题不必计算4。似=(2-x)(x+1)2/4g0(x)=(x-(-l))(琶)=(x+l)(x-1)2/4耳3(x)=%(x)・4+go(x)・2%(x)=(2-x)(x+1)2+；(x+1)(x7)2 172品")=3.5625利用构造基函数方法做插值多项式被广泛地应用在不同的插值条件中。例5.9给定/So)=必)，」'(入0)=%，/(甬)=必构造二次插值多项式函数。解：设』⑺=4。)尢+*x)M+〃(x>%这里£。(乃心(x)g(x)均为不高于二次的多项式，它们分别满足4(%)=17Go)=0=045)=04£；（而）=0«块硝=04(X1)-04（丸）=14（勺）=1于是与（X）可表示为X-XX-jjr舄（x）=（仆x+%）Ly0+a1—乃+a2（x-Xo）（x-XiMo勺一再々一%由45）=1,得％而+3T、n^o-1+（旬而+4）=0由4（而）=U，得勺一再而一再g二一—%所以而-演，而一々£0（x）=（2xoF-x）.（xf）所以瓦-西（瓦一公） 173他）=卢鼻同理（为一%）由4（%）一】,叼（/一占）=11a2=所以％一再/）=（—五）所以而一々 174用牛顿差商插值也能构造埃尔米特插值。对给定的插值点的函数值和一阶导数值（为=0，L…/定义序列Z（l==Xo，Z2=々/3=々,…，即z©=Z2“1=西/=0,1计算一阶差商时：加5引/区卜八味）Z2iZ2i-1由九"『和瞿力"】=〃而），取■/[zw，Zr+J=/'&）/=0,1，…述即构造差商表中用/（而）'/（再）=、,（/）代替/[^o，zjJ[Z2/3…,z2x+1]其余差商计算公式不变，得到差商型埃尔米特插值公式:2X4-1H2/1⑸=/%]+£/区,Z],…,ZJJt-1（5.26）（x-z。）…（x-Zj）其中Z以=z2川=4，例5.10用下列数据构造埃尔米特插值多项式，并计算/（1.36）。X1.2/W0.6广⑺0.51.41.60.91.10.70.6解：计算差商。 1751.20.61.20.60.51.40.91.551.40.90.7-4-451.61.11.01.5131451.61.10.6-2-17.5-76.25-553.125%(x)=0.6+05。-12)+5(x-1.2)2-45(x-1,2)2(x-1.4)+145(x-1,2)2(x-1.4)2-553.125(x-1.2)2(x-14)2(x-1.6)凡(1.36)=0.8655例5.11求尸（x）e%，使P（Xi）=/（Xi）G=°，l，2）及尸（再）=/'5）插值多项式及其余项表达式。解：这里给出了四个条件故可构造三次插值多项式尸（x）e“3,由「（不）=/（大刀"°，1，2）,可用Newton均差插值，令p（x）=/（%）+/［%+/［勺,xl,x2］（x-x0）（x-再）+a（x-7i））（x-x1）（x-x2）显然它满足条件P（xJ=/（%）（'=°，L2）,a为待定参数。由PG）=/'（々）可得-）=f［x0，再］+/［%,再,与］（々-而）+a（x1-x0）（x1-x3）7'a）a=11〃再）-/（飞，再】I解得：x「x。J于是就可得到插值多项式。它的余项为 176R3(x)=/(x)-p(x)=；/')(b(X-Xo)(X-Xi)2(X-X2)*T^在飞与X?之间，而a«勺《再4刍4小。5.6分段插值5.6.1龙格(Runge)现象在构造插值多项式时，根据误差表达式(5.9),你是否认为多取插值点总比少取插值点的效果好呢？答案是不一定。如果被插函数是高次多项式，那么多取插值点比少取插值点效果好；但对有些函数来说，有时点取的越多，效果越不近人意。请看下面的例子。j(x)—1给定函数1+25/,xe［-1,1］。构造io次插值多项式几⑶。对［-1,1］作等距分割，取10,■，构造Il+25xj,10次插值多项式4。(工)如图5.4所示。从图中可以看到，在零点附近，4o(x)对/(X)的逼近效果较好，在了=一0.90,-0.70,0.70,0.90时误差较大。 177图5.4Ao(x)和/(X)下面列出Ao(x)和/(x)的几个插值点的数值：X-0.90-0.70-0.50-0.30/wAo（x）1.578720.47060.07547-0.226200.137930.253760.307690.23535这个例子是由龙格（C.Runge）提出的，也称插值多项式在插值区间发生剧烈振荡的现象为龙格现象。龙格现象揭示了插值多项式的缺陷。它说明高次多项式的插值效果并不一定优于低次多项式的插值效果。在插值过程中，误差由截断误差和舍入误差组成。式（5.9）给出的是截断误差，它是插值函数°卜）与原函数/（“）的误差。另外由节点*计算产生的舍入误差，在插值计算过程中可能被扩散或放大，这就是插值的稳定性问题。而高次多项式的稳定性一般比较差，这从另一角度说明了高次插值多项式的缺陷。5.6.1分段线性插值既然增加插值点并不能提高插值函数的逼近效果，那么采用分段插值的效果又如何呢？若对给定区间［0'可作分割0=瓦＜々4=方，在每个小区间［不上作/（x）以0和1为节点的线性插值，记这个插值函数为P（x）=&（x）,则&（x）=X/（勺）+X4/（租1）%4X44+1石―石+1玉+1一&（5.27）把每个小区间的线性插值函数连接起来，就得到了,（X）以剖分（节点）a=/＜々＜公=3的分段线性函数P（x）。P*）在［辱上为一个不高于一次的多项式，事实上「（力是平面上以点aJa））为折点的折线。由线性插值误差公式，当xw［x”为+J时有fW-p（x）=/（x）-gj（x）=（x-个）（x-%）2!（5.28）|/（x）Mx）14争（x-xj 178因而«学・：（%即2=今（％-西）2其中圾=max|/"（x）|,a4x4小于是，当区间分割加密，恶生1"加不）—°时，分段线性插值收敛于/（X）。事实上，只要/（X）连续，分段线性插值序列就能收敛于/（x）。分段线性插值算法简单，只要区间充分小，就能保证它的误差要求。它的一个显著优点是它的局部性质，如果修改了某节点a的值，仅在相邻的两个区间为受到影响。分段线性插值的缺点是在插值节点处不光滑。图5.5给出分段线性插值P（x）（虚线表示）和/（X）的图形，可以看到分段线性插值的效果明显好于整体的拉格朗日插值的效果。图5.5分段线性插值P（X）和/（X）例5.12对下列数据作分段线性插值，并计算/（1.2）,/（3.3）o12121612P（x）=&（力为%/（西）+'//（F）,xe[&,kJ解：石一玉+ixM-V1.2€[-1,2] 1791.2-21.2+1P（1.2）=g】（x）=-l-2X5+2+1X1=2.06673.3-93.3-3(3.3)=g3(，)=-6x6+6x12=6.35.7三次样条函数在制造船体和汽车外形等工艺中传统的设计方法是，首先由设计人员按外形要求，给出外形曲线的一组离散点值万)1工…’，施工人员准备好有弹性的样条(一般用竹条或有弹性的钢条)和压铁，将压铁放在点｛X”“｝的位置上，调整竹条的形状，使其自然光滑，这时竹条表示一条插值曲线，我们称为样条函数。从数学上看，这一条近似于分段的三次多项式，在节点处具有一阶和二阶连续微商。样条函数的主要优点是它的光滑程度较高，保证了插值函数二阶导数的连续性，对于三阶导数的间断，人类的眼睛已难以辨认了。样条函数是一种隐式格式，最后需要解一个方程组，它的工作量大于多项式拉格朗日型式或牛顿型式等显式插值方法。定义给定区间卜'可上"1个节点a=为＜…＜五"和这些点上的函数值，了⑷=0,1,…,”若$(x)满足S(aj)=yifi=0.；S(x)在每个小区间【小小］上至多是一个三次多项式;£(x)在网上有连续的二阶导数，则称£("为了(X)关于剖分0=/＜々/=3的三次样条插值函数，称…'、为样条节点。要在每个子区间以‘玉+1］上构造三次多项式S(x)=S(须)=白/3+«/+xe［乐k皿=0,1,…*-1,共需要几个条件，由插值条件S(4)=X，i=°，L…，”提供了“+1个条件：用每个内点的关系建立条件S(.+0)=S(4-0)s'a+o)=s'a-o)S"&+0)=S"a-0),J=0,1,•••,«-!又得到%-3个条件；再附加两个边界条件，即可惟一确定样条函数了。用待定系数法确定了构造样条函数的存在性和惟一性。在具体构造样条函数时一般都不使用计算量大的待定系数法。下面给出构造三次样插值的"关系式和m关系式的方法。 1805.7.1三次样条插值的M关系式引入记号m=s'a),j…,修。用节点处二阶导数表示样条插值函数时称为大m关系式，用一阶导数表示样条插值函数时称为小w关系式。问题5.8给定插值点(0又＞""°工，怎样构造用二阶导数表示的样条插值函数，即怎样构造M关系式？假设S"⑷=M"=O,L"/.由于S"®)在【小小」上为线性函数，故在［不，上做”的分段线性插值函数：s"(xj=二1^±l此+士工=0i=-1令%=玉+1-玉,得到s"8)=①二土此+王1宝盟(5.29)=0i=1,2,..-1对S〃@)积分两次有%)=&W~M+若必"x+dM+M+i+C(x“I-x)+Q(x-Xj)(5.30)64将sa)=>y”sa+i)=%i代入式(527)可解出c=——"Nid="+i—々瞩+i 181一%丁’一了6故在［西，为+J上有S（X）=（%一-）3此+。一炉此+1+（Xi+i-x）M+（x-个）y”i々K4「x）M+-M+J6xe[xj（xi+1],j=l,2,...,?j,、（5.31）S（x）在每个小区间上具有不同的表达式，但由于8（力在整个区间[°力]上是二阶光滑的，故有卬（再+0）=£'（七-O）,i=O,l,--,«-l列出每一个关系式¥（为）=M,a）,再经计算得：“K-1+2此+4此+i=d”i=1,2,…,%-1（532）其中：\=-^—%+%丛=Fd=6x+「x_'%+%%%,=6乂0此再+1]由式（5.32）得到"+1个未知数的“7 182个方程组。现补充两个边界条件，使方程组只有惟一解。下面分三种情况讨论边界条件。（1）给定的值（%=°'M*=°时，称为自然边界条件），此时"-1阶方程组有力-1个未知量=L2,…,%-1,即 183，2A体-Ma1百2乙•--—4-22Ai-2M.2麓7“I12>(5.33)（2）给定S（x。）=S（x»）=%的值，它们分别代入S'（x）在［「，舷中的表达式，得到另外两个方程：2M）+a=「［Mx（），xJ-多卜色%M」+2%=二［外加小区］卜4如于是需要解%+1阶的方程组：’21以12%“224%L（3）被插函数以'一而为基本周期时，即3。=居，即S（'o）=£（x）S'（x。）=S（xJS《。）='"（/）；即帆0=阳，弧=%。此时化为力个变量、力个方程的方程组。样条插值构造的M关系式是对角占优的三对角带状矩阵，可用第3章中的追赶法求解。例5.13给出离散数值表： 184xi1.11.21.41.50.40000.80001.65001.800()取弧=a=0,构造三次样条插值的M关系式，并计算/(1.25)。解：由题中的数值，计算得%=0.1,%=0.2,%=0.1(4=0.6667,J4=0.3333=0,3333,(为=0.6667由=5,心=-55由弧=a=°的边界条件，得r206667)因)510.66672八弧「-55」解得弧=13.125,必=-31.875。因此，三次样条插值的分段表达式为r21.875a3-72.1875x2+83.1875x-32.875xe[l,1.1,2]S(x)=<-37.5x3+141.5625/-173.3125x+69.725xe[1.2,1,4]53.125x3-239.06259+359.5625x-178.95xe[1.4,1.5]特别地，/(1.25)^(1.25)=1.0336o详细的程序和算例请看5.8节。5.7.1三次样条插值的m关系式问题5.9给定插值点=°工…〃，怎样构造用节点处一阶导数表示的样条插值函数，即怎样构造加关系式? 185对给定的插值点（"⑻）/=°工…先假定已知S'。。飞,在每个小区间【不必］上做埃尔米特插值，那么在整个［飞，5］上是分段的埃尔米特插值，在［为4+J上£（力的表达式为4逸_1+2逸+m+i=c”i=1,2,",«-l其中：四=1-4q=3（4，［再_卜xj+必丁［演,玉+1］）再附加两个边界条件，即可解出幽i的值。附加的边界条件情况同M关系式中的讨论类似，不再详述。 186图5.6样条插值图示5.7程序示例程序5.1给定（为，M），'=0，1一，篦，构造牛顿插值多项式MS），*，互不相同。算法描述输入用值，及（a，x），i=o,L…/；记/缶）=%forj=0,!,•••,«了小用,…,凝卜八。.，…-Ha]计算差商敌一而其中/闻=/（均）对给定的入，由+（x-Ko）（x-再）/[而,和电卢…+（x_5）（x_占）一•（X_X»）/[Xo,勺,一X*]计算”（X）的值；输出M（x）。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII〃purpose：（x_i,y_i）的牛顿插值多项式〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 187{doublex；doubley；}POINT：intmain(){intn；inti,j；POINTpoints[MAX_N+1]；doublediff[MAX_N+1]；doublex»tmp,newton=0；printf(44 188lnputnvalue：M)；〃输入被插值点的数目scanf("％d”,&n)；if(n>MAX_N)(printf("TheinputnislargerthenMAX_N,pleaseredefinetheMAX_N. 189")；return1；)if(n<=0)(printf("Pleaseinputanumberbetween1and%d. 190",MAX_N)；return1：)//输入被插值点(x_i,y_i)printf(44Nowinputthe(x_i»y_i),i=0»...»%d: 191",n)；for(i=0；i<=n；i++)scanf("%If%IF',&points[i].x,&points[i].y)； 192printf("Nowinputthexvalue：")；//输入计算牛顿插值多项式的x值scanf(4t%ir,&x)；for(i=0：i<=n；i++)diff[i]=points[i].y；for(i=l；i<=n；i++){for(j=n；j>=i；j-)(diff[j]=(diff[j]—diff|j-l])/(points[)].x-points[j-l-i].x)；}//计算f(x_0,…，x_n)的差商}tmp=1；newton=diff[0]；for(i=0；i<=n；i++)(tmp=tmp*(x-points[i].x)；newton=newton+tmp*diff[i+1]；)printf(4inewton(%f)=%f 193M,x,newton)；〃输出return0；)计算实例给定sin11。=0.190809,sin12°=0.207912,sin13°=0.224951,构造牛顿插值函数并计算sinl1。30。程序输入输出Nowinputthe（x_i>y_i），i=O,…，2：110.190809120.207912130.224951NowInputthexvalue：11.5newton（11.500000）=0.199369 194程序5.2给定插值点…，"和二阶导数的端点值％，圾，用肠关系式构造三次样条插值多项式S（x）,求在给定点X处S（x）的值。算法描述1.输入忽值，及（不乃），'=0，1「，”如M,要计算的函数点X；2.for/=1ton-13段d——计算%+%,乌=1-4,d=6_MF；1%+%i%如f其中%=玉+1一玉：3.求解方程组对给定的x,由y、_（%-次跖+（Xr）3此+1S（x）=疏一（■+1-芯）乃+（工一个）乃+1国A一戈（大+1-x）M+（x-QM+J64.计算出S（x）的值；5.输出S（x）o程序源码 195lllllllllllllllllllllllllllllllllllllllllllllllllll//purpose：给定%,%值的三次样条插值多项式〃lllllllllllllllllllllllllllllllllllllllllllllllllll#include 196scanf(“%d"，&n)；if(n>MAX_N){printf(4TheinpulnislargerthanMAX_N,pleaseredefinetheMAX_N. 197")；return1；)if(n<=0){printf("Pleaseinputanumberbetween1and%d. 198",MAX_N)；return1}〃输入插值点(x_i,yj),MO值和“〃值printf(**Nowinputthe(x__i,y_i),i=O,...»%d： 199",n)；for(i=0：i 200”,points|0].x,points[0].x)；return1；}//计算M关系式中各参数的值h[0]=points[1].x-points[0].x： 201for(i=l；i<；i-H-){h[i]=points[i+l].x-points[i].x；b[i]=h[i]/(h[i]+h[i-l];c[i]=l-b[i]；d[i]=6*((points[i+l].y-points[i].y)/h[i]—(points[i].y-points[i—l].y)/h[i—l])/(h[i]+h[i-1])；}〃用追赶法计算Mi,…，n—\d[l]-=c[l]*M[0];d[n-l]~=b[n-l]*M[n]；b|n-11=0；c[l]=0；v[0]=0；for(i=l；i 202+(p*points[k],y+q*points[k+l].y)/h[k]—h[k]*(p*M[k]+q*M[k]+q*M[k+l])/6；printf(**S(%f)=%f\x,S)；//输出getchar()；return0；)计算实例给定离散点(1.1，0.4),(1.2,0.8),(1.4,1.65),(1.5,1.8),=圾=°,用河关系式构造三次样条插值多项式贺力，计算£.25)。程序输入输出Inputnvalue：3Nowinputthe(x_i,y_i),i=0,…，3：1.10.41.20.81.41.651.51.8NowInputtheMOvalue：0NowInputtheMnvalue：0NowInputthexvalue：1.25S(1.25000)=1.033594本章练习计算证明题1.作出插值点(-1.00.3.00),(2.00,5.00),(3.00,7.00)的二次拉格朗日插值多项式与炽)，并计算上2(°)。2.作出插值点(-2.00,0.00),(2.00,3.00),(5.00,6.00)的二次拉格朗日插值多项式'(入)，并计算£式-1.2),(£2(12)。 2031.作出下列插值点的三次拉格朗日插值多项式。(1)(一1,3),(0,一1/2),(1/2,0),(1,1)(2)(-1,2),(0,0),(2,1),(3,3)4.设/M",且/⑷=〃力=0,证明J"切圾，其中圾=max|/*(x)|5./CO=石在离散点有/(81)=9,/(100)=10,(121)=11,用插值方法计算芯=105的近似值，并由误差公式给出误差界，同时与实际误差作比较。6.给出函数表：X-1.002.003.004.00“X)3.005.007.005.00作出差商表，写出牛顿插值公式，用M.2)计算,(12)的近似值。7.要作三角函数l°gx的函数值表，已知表值有5位小数的近似值，要求用线性插值引起的截断误差不超过表值的舍入误差，试决定其最大允许步长。&/(%)=/-125/+237/-999,计算差商力2°,2”J[20,2],…,27]以及力2°,21,…,28]。9.给出函数表：X1.051.101.151.20“X)2.122.202.172.32构造分段线性函数，并计算/QBS)和/(1*5)的近似值。10.给定数据/（°），/（D，/'⑴，作出二次插值多项式，并写出插值余项。 2049.给定数据/⑶=5.00,/（5）=15.00/（5）=7.00,作出二次插值多项式，写出插值余项，并计算/07）的近似值。10.给定数据⑶/'（3）,作出三次插值多项式，并写出插值余项。11.给定数据,（°）=1°，/（1）=°75,/（3）=°.25*3）=°.56,作出三次插值多项式,并写出插值余项。12.给定数据,（°）=0，,（1）="（0）=0,*1）=1/（1）=°,构造四次插值多项式，并写出插值余项。13.给出下列数据：X-2.00-1.001.002.00一4.003.005.0012.00试求满足上列数据和S"L2）=0,S”（2）=0的三次样条函数，并计算S（0）的值。14.给出下列数据：X-1.000.001.003.0()“X)2.003.004.0029.00第6章曲线拟合的最小二乘法6.1拟合曲线通过观察或测量得到一组离散数据序列a’"）/=°/，…，冽，当所得数据比较准确时，可构造插值函数°（幻逼近客观存在的函数y=y（x）,构造的原则是要求插值函数通过这些数据点，即。（不）=必」=1,2,…，加。此时，序列。=（。51），。①），…，W（xG）r与y=31，乃，…，％）「是相等的。 205如果数据序列a‘x）/=°'L…，含有不可避免的误差（或称“噪音”），如图6」所示：如果数据序列无法同时满足某特定函数，如图6.2所示，那么，只能要求所做逼近函数程（X）最优地靠近样点，即向量Q=…,。（xQ/1与y=（m,也，…，%）,的误差或距离最小。按Q与y之间误差最小原则作为“最优”标准构造的逼近函数，称为拟合函数。图6.1含有“噪声”的数据图6.2一条直线公路与多个景点插值和拟合是构造逼近函数的两种方法。插值的目标是要插值函数尽量靠近离散点；拟合的目标是要离散点尽量靠近拟合函数。向量°与丫之间的误差或距离有各种不同的定义方法。例如：用各点误差绝对值的和表示：次舄川2-1用各点误差按模的最大值表示: 206凡=max|0（丹）-乃|用各点误差的平方和表示:”号?（唉）-方”地⑶_理其中R称为均方误差，由于计算均方误差的最小值的方法容易实现而被广泛采用。按均方误差达到极小构造拟合曲线的方法称为最小二乘法。本章主要讲述用最小二乘法构造拟合曲线的方法。在运筹学、统计学、逼近论和控制论中，最小二乘法都是很重要的求解方法。例如，它是统计学中估计回归参数的最基本方法。关于最小二乘法的发明权，在数学史的研究中尚未定论。有材料表明高斯和勒让德分别独立地提出这种方法。勒让德是在1805年第一次公开发表关于最小二乘法的论文，这时高斯指出，他早在1795年之前就使用了这种方法。但数学史研究者只找到了高斯约在1803年之前使用了这种方法的证据。在实际问题中，怎样由测量的数据设计和确定“最贴近”的拟合曲线？关键在选择适当的拟合曲线类型，有时根据专业知识和工作经验即可确定拟合曲线类型；在对拟合曲线一无所知的情况下，不妨先绘制数据的粗略图形，或许从中观测出拟合曲线的类型；更一般地，对数据进行多种曲线类型的拟合，并计算均方误差，用数学实验的方法找出在最小二乘法意义下的误差最小的拟合函数。例如，某风景区要在已有的景点之间修一条规格较高的主干路，景点与主干路之间由各具特色的支路联接。设景点的坐标为点列…，.；设主干路为一条直线°（x）=a+&x,即拟合函数是一条直线。通过计算均方误差0S力）最小值而确定直线方程（见图6.2）。Q3,8）=E（。5）-乃）2=Z（a+M-乃尸2-12-16.1线性拟合和二次拟合函数线性拟合给定一组数据（小乃））=°」，施，做拟合直线「。）=0+以，均方误差为W）=宜（P（Xj）-乃尸=支3+也-乃尸i-1i-1（6.2） 207OSA）是二元函数，°（a，与的极小值要满足8Q（a.b）da=2£（。+也-乂）=024-db~=2Z（a+^Xi-M）为=°2-1整理得到拟合曲线满足的方程:称式（6.3）为拟合曲线的法方程。用消元法或克莱姆法则解出方程:jn-jw加£不乃一工不工》2-12-12-1x1315y11io162111122212232513132930313612141617例6.1下表为P.Sale及R.Dybdall在某处作的鱼类抽样调查，表中入为鱼的数量，丁为鱼的种类。请用线性函数拟合鱼的数量和种类的函数关系。 2084042556062647072100130V13142214212124172334解：设拟合直线P（x）=°+5x,并计算得下表:编号Xyxy2X1131114316921510150225316111762564211225244152212264484*a一*■•21130344420169009563441891361640将数据代入法方程组（6.3）中，得到:2195634495661640b18913解方程得：a=8.2084,3=0.1795拟合直线为：P（x）=8.2084+0.1795X二次拟合函数给定数据序列=°，L…，加，用二次多项式函数拟合这组数据。设p（x）=a0+a^+a/2,作出拟合函数与数据序列的均方误差：(6.4)=23（再）一乃尸=£（%+alxi+a2x.-*）22-12-J 209由多元函数的极值原理，03。，生，。2）的极小值满足■=2£3()+。洛+a2m-乃)=02-1=22仇+aixi+a2xi=02-1融=2£3o+4玉+%/-»)x；=02-1整理得二次多项式函数拟合的法方程:次支不2-1滞2-1M2-1MMIX2-1M2-1M2-1湿少2-1源2-1j/f汽X也2-1(6.5)解此方程得到在均方误差最小意义下的拟合函数P（x）o方程组（6.5）称为多项式拟合的法方程，法方程的系数矩阵是对称的。当拟保多项式阶力＞5时，法方程的系数矩阵是病态的，在计算中要用双精度或一些特殊算法以保护解的准确性。例6.2给定一组数据，如下表。用二次多项式函数拟合的这组数据。解：设阿=/+守+呼2由计算得下表:-3-2-10123y4230-1-2-5Xy2A小3Ax4-34-12936-2781-22-448-816-13一313-110000000 210将数据代入式（6.5）,相应的法方程为：+0,+23a2=1«0aQ+281+0%=-3928ao+0%+196%=-7解方程得：=0.66667,01=-1.39286，°，=-0.13095二*（，）=0.66667-1.39286x—0.13095-77拟合曲线的均方误差：i-1i-1=3.09524结果见图6.3o6.1解矛盾方程组在6.2节中用最小二乘法构造拟合函数，本节中用最小二乘法求解线性矛盾方程组的方法构造拟合函数。给定数据序列，作拟合多项式P（x）+乎+”・+,如果要求外。）过这些点，那么有矛盾方程组： 211外(再)=%+a1xi+…+a*x；*(i=1,2,->»+1)即:■+好+.•.+"=%%+,勺+…+/X；=72(6.6)…+/x：我们作一辅助函数0防即,-4）=£3个）-珀224篇=工（。0+。丙+一+久邸-乂）2-1°对每这是自变量为/的多元函数，要使°达到最小值，采用多元函数求权值的方法,一个自变量的偏导数为0。哭=2之4+a洛+…+4胃-%)=03%字=22(%++…+%x：-R个=0此tr窘2热"+-"”=0整理成以。。n卜…为未知数的线性方程组 212JffJffffffff2-12-12-12-1M㈱MMZ、，2X—,*4.1、，大+,万营+…+%》玉=W不必1i-1i-1i-1Z斌JftMt㈱W：+，>+…+士m、2-12-12-12-1整理成矩阵形式3〃=C,其中：>x*+iI［小》〜»：EVS＞：»：+◎>>，.这是一个5+D*5+1）的对称方程组，称为法方程。只要8非奇异，就可以得出唯一解。这就是矛盾方程组的最小二乘解。有什么快捷的方法来求法方程的解？把矛盾方程组（6.6）写成矩阵形式忆=»,其中1玉^1…X：ao必 213容易验证B=/弘c，法方程就是/忆="九例如，拟合直线P。）=&+0户的矛盾方程组Wax=幺%的形式如下:化简得到与（6.3）相同的法方程:在线性代数中，我们知道，关于方程组月入=小，若秩（4切=秩（/），则方程组无解，这时方程组称为矛盾方程组。在数值代数中对矛盾方程计算的是在均方误差mm极小意义下的解，也就是在最小二乘法意义下的矛盾方程的解。定理6.1将证明，方程组的解就是矛盾方程组月入=小在最小二乘法意义下的解。定理6.11.』为次行〃列的矩陈，〃＞吃力为列向量，秩（工）=花，H⑷=工%称为矛盾方程的月入=占的法方程，法方程恒有解。2.X是mmII/X1的解，当且仅当X满足=,即X是法方程的解。证明：1）对上作行初等变换尸，使 214（产工）（产工）r=...秩（（日）（口力力而秩（（取）（尸⑷］）卬秩⑷］）=秩（力〜）=« 215方程组有解而且解惟一存在。2）设X痛足任取y=x+（y-x）=x+e,则\\AX-b\^=（AX-b+Ae/（AX-b+Ae）={AX-bf{AX-b）+2⑷「（AT-3）+（工e）11（㈤-\\AX-b\^+\\Ae|£+2er（^AX-Arb）=\\AX-b\^\\Ae\^>\\AX-b\^由于丫是任取的，故法方程组力=力的解为极小问题的解。事实上，对离散数据（小必）/=0」,…，搐作〃次多项式曲线拟合，要计算=。（曲,,,Jfl=£（&+。通+…a*x：-x）27的极小问题。这与解矛盾方程组阳+。丙+-+%*：=% 216 217故对离散数据"QL…，.所作的万次拟合曲线°。+/X+…+%/，可通过解下列方程组求得:例6.3给定一组数据，见下表，用二次多项式函数拟合这组数据。X-3-2-10123y4230-1-2-5解：记二次拟合曲线为•/@)=旬+。1，+叼,，形成法方程:2801960280 2181-39-7423oT-2T1oO1711~241-39必为不7S2.172-i-l/I=y1-39-7r1=X/012L:a-与。叫0280解方程得：°=0.66667,1=-1.39286,2=-0.13095:./⑺=0.66667-1.39286x-0.13095^例6.4给出--组数据，见下表。用最小二乘法求形如,（X）+的经验公式。XX-3-2-124yy14.38.34.78.322.7解：列出法方程:3再3r^3勺3勺3々e大5E-2-15S.Z%4954 219-27法方程为:-8-18111536364954bf\'58.3'1062解方程得到：a=10.675,B=0.137拟合曲线为：/（*）=10.678+0.137-有些非线性函数经过转换以后可化为线性函数计算。例如，'X令y例6.5求一个形如y=ae的经验函数公式，使它能够和下列数据相拟合。X12345678y15.320.527.436.649.165.687.8117.6解：化经验公式为线性形式，对经验公式的两边取自然对数有 220Iny=lna+云由矛盾方程组Ina+%= 221yxIna+bx2=In必Ina+b&=In/得到法方程组 222r836Win03.0197,36204)(b]63.9003,解方程组得Ma=2,3469,b=0.291"=*湖=11.4375,y=11.4375严2x拟合曲线的均方误差为:88珀2=0.02624012-1i-1计算结果见图6.4o图6.4拟合曲线图示例6.6解矛盾方程组 223\+x2+x3-2再+3x2-x3=-l2xj+5x2+2x3=13丑一弓+5/=-2解：写出法方程组/「凡*=/'丫，即rl113252135-1'1123解法方程得:々=1.5917,々=0.5899,r151119*3=0.7572«1-1图6.3拟合曲线与数据序6.4权有的实际问题中，已知数据不一定都是一次观测的结果，对于不同的a，x）可能观测次数不同，在矛盾方程组中，一组a，乂）确定一个方程，而最小二乘解对每个方程来说都还存在误差，这样，把每对a，x）都同等看待就不太合理，希望观测次数多的（即可靠性大些的）数据在矛盾方程组中占的比重大些。为了使问题的提法更有一般性，通常把最小二乘法中偏差平方和改为加权平方和，使I晓=冬（项八（再）2-1V.]2,、、C!最小，其中w（xj>°称作权。 224。（曲小…同样，可由误差函数（%）&2必为）西2以西）可D—4）=Wwa）（p*a）-必yh求极值，得到法方程3〃=C,其中£啾为）胃--工讽3-工以“彳…^w（Xj）z*+1＞尸（西）犬+2事实上，只需对矛盾方程组忆=、作一些修改,每个方程乘上仃㈤,得至产=尸其中一”⑷一Jw(q)s/w(aw)xmJ中(电电Jwa)x；•••屈石芯Jwg)考■■■打⑹专Jw(/)x：--Jw(/)x：这样，就使8=新'晒0=柄行，仍然可用前面的方法，对g=7，左乘柄'用&r柄^=柄「尸求得最小二乘解。6.5用正交多项式作最小二乘拟合我们不仅可用多项式来拟合函数，还可用一般的函数来拟合。定义6.1若。o（x），仍（力，…,W苫°[4句，如果40例⑶+%例⑶+…+%R(X)=0 225当且仅当/=白1=一.=%=0时成立，则称例（幻'例（力在〔兄切上线性无关，称。（.为［“力〕上线性无关族。最小二乘法的一般提法为：对给定的一组数据5Ji"，=L2,…,我），要求在函数类 226定义6.2：⑸%）,（%）%（&）%）称为叼3与伙（砌关于点集㈤荒的内积。.£勺•（6=8,伙）_这样，法方程式可简写为刀，记为B"=d，其中（仙，例）（仍，例）（%，00）30,何）（仍，例）（怯，仍）（00,也）（仍，俗）（例，。黑）6㈤.3,仍）Os）.（%,仍）（例,他）（%，。0）|用=（。0,何）（仍,例）（俗，程）（00，供）（仍，例）(Gw*)称为克莱姆行列式，记作5。定理6.2GJ。的充要条件是Po（x），研（X），…，On线性无关。证明：若存在做W长体=0，1，…/）使/例+。1例+…+/0*=0对此式两边分别取与例0），例（x）J[0”的内积得：’（仙，仙）/+（研，。0）的+…+（俗,例）%=0（例，仍）即+（仍，何）/+…+（外，例）％=0,（例,0*）。0+（例，。*）/+3+（。『供）勺=0 227这是一个以旬'生'…'怎为未知数的齐次方程组，有非零解的充要条件是系数矩阵行列式等于零，于是5*°的充要条件是方程有全零解，即全为0,所以例（X），何（力，…，0n线性无关。证毕。由于法方程有惟一解的充要条件是[*°,因而0。0），研（X），…，例（线性无关也是法方程3〃=d有惟一解的充要条件。特别当价（立研（X），…，％取为1，工，…，/时，由于｛/）岂是在中的线性无关函数族，因而必有最小二乘解。用印a双l,x,…上的多项式拟合，需要解一个万+1的线性方程组，且当〃取得大一此时，法方程的系数矩阵会出现病态。从系数矩阵B的形式看，里面的元素都是些内积，是否能取某些函数族，使对非对角元素全变为0?如果有这样的函数族，那么方程容易解，病态也得到改善。 228定义6.3函数族优（x），研（X），…，死如果在点集伪比1上满足潞（叼•，怯）=£州&）%（既）伍（西）2-1{1>0j*kj=k称仰（x）,例（x）,…,/为点集⑨）工带权w（"i）的正交函数族。例6.7三角函数族Lsm冗cosx,sin2x,cos2M…在［-n，句上是正交函数族（权.5）三］）(1,1)=fdx=2几实际上L(sinnxtsinmx)=Jsinnxsinmxdx0,whnit\k,m-n(cos«x,cos?«x)=fcos??xcoswxdx=〈n.m531,2，•一10,w(cosnxrsinmx)=Jcos阀xsinmxdx=0,n,m=1,2：如果拟合函数在印°”｛仙，研「：仰》上取,且｛伙｝to是正交函数族，则法方程式成为:一⑶㈤_季（钠小）趣-7；~直接可得到1,不用解线性方程组了。的d(a=ML团|max®,,依)0S尽*min(wg,保)且。**容易估算,是否病态也容易判断。 229完成以上工作的关键在于如何构造正交函数族。正交多项式是最简单的正交函数族。常用的正交多项式如：Chebyshev（切比雪夫）多项式、Legendre（勒让德）多项式等。现在我们根据给定结点及权函数w（x）＞°,造出带权w（x）正交的多项式族⑦⑶应。+用递推公式表示如下：,Po(x)=1 230Po=1㈱㈱3例）=工必演）回（%）=21=42-1MM4（砺o，Po）=£w（生）否p：（壬）=£为=102-1i-1„_M，A）_5/一^T548，Po）=»i=58i-i所以:3，为）=58=29（Po，p0）42Pi=x--（Pi，0）=a（研，介乳「I）q（啊,Pi）_5（P「Pi）2（P”Pi）-5（Po-Po'）W。，分）＜37所以：1（巧，/）55552外㈤二炽-叼加一凤％工工一汽一戏一厂x-5x+5（外#2）=Z（x；+5/=42-148，外）=£6-5号+5）制=214与=㈤=I所以：2»2#2）2 231得:2937.5、J穴cy=——+——(x--)+—(x3-5x+5)2522493+—x~—1026.6程序示例程序6.1用线性函数,（*）=a+5工拟合给定数据（'”£）「=L2,…，耀-1算法描述输入次值，及（和乂）,7=12…，耀7。解方程组:M-lJM-1输出p(x)=a+8x程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII//purpose：(x_i,y_i)线性拟合函数〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 232POINTpoints[MAX_N]；staticdoubleul1,ul2,u21,u22,cl,c2；doublea,b,tmp；printf(44 233lnputnvalue：")；〃输入点数mscanf("%d”,&m)；if(m>MAX_N){printf("TheinputnislargerthenMAX_N,pleaseredefinetheMAX_N.\rT)；return1；)if(m<=0){printf("Pleaseinputanumberbetween1and%d.\rT,MAX_N)；return1；}printf(44Nowinputthe(x_i,y_i),i=0,…，％m—1)；//输入点for(i=0；i<=m；i++){scanf(4i%lf\&tmp)；printf[i|.x=tmp；scanf&tmp)；printf[i].y=tmp；//scanf&printf[i].x,&printf[i].y)) 234for(i=0：i<=n；i++)〃列出方程U(a,b)T=c(u21+=points[i].x；u22+=points[i].xpoints[i].x；cl+=points[i].y；*c2+=points[i].ypoints[i].y；)u12=u21；ul1=m；//求解a=(cl*u22-c2*ul2)/(ull*u22-ul2*u21)；b=(cl*u21-c2ull)/(u21ul2-u22ul1)；printf(t4Solve：p(x)=%f+%fx 235*\a,b)；〃输出return0；}//EndofFile计算实例用线性函数拟合下列数据(例6.1)oX1315162122232529303136y1110II12121313121416174042556062647072100130y13142214212124172334程序输入输出inputmvalue：21 236Nowinputthe(x_i>y_i),i=O»...»20：13111510161121122212231325132912301431163617401342145522601462216421702472171002313034Solve：p(x)=8.208408+0.179522x程序6.2用形如P(R的函数拟合定数据(%,%),i=1,2,…，.-1算法描述输入沏值，及(XQi)J=L2,…盟7。解方程组：TM-lM-1'i-o1n«=i-o就一i——1I辰I•-i,i-0i-0i-Oj输出所求拟合函数P(x)=ae”。程序源码lllllllllllllllllllllllllllllllllllllllllllllllllllllllll//purpose：(x_i,y_i)形如y=aExp(bx)的拟合函数〃lllllllllllllllllllllllllllllllllllllllllllllllllllllllll#include 237doublex； 238doubley；}POINT；intmain(){intm；inti；POINTpoints[MAX_N]；staticdoubleul1,ul2,u21,u22,cl,c2；doubleA,B,tmp；printf(4< 239lnputnvalue：v)；〃输入点数mscanf("%cT,&m)；if(m>MAX_N){printf('TheinputnislargerthenMAX_N,pleaseredefinetheMAX_N. 240M)；return1；)if(m<=0)(printf(**Pleaseinputanumberbetween1and%d. 241",MAX_N)；return1；}printf(44Nowinputthe(x_i»y_i)，i=0,...»%d: 242M,m-1)；//输入点for(i=0；i<=m；i-H-){scanf(t4%lf,\&tmp)；printf[i].x=tmp；scanf("%lf”,&tmp)； 243printf[i].y=tmp：//scanf(44%lf%lf"，&printf|i].x»&printf[i].y)；}for(i=0；i<=n；i++)//列出方程X(a,b)T=Y{〃其中丫为对原方程取对数后的数据u21+=points[i].x；u22+=points[i].xpoints|i|.x；cl+=log(points[i].)y；*c2+=points[i].ylog(points[i].y)；)u!2=u21；ul1=m：//求解a=(cl*u22-c2*ul2)/(ull*u22-ul2*u21)；b=(clu21—c2ul1)/(u21ul2—u22ul1)；printf("Solve：p(x)=%Ifexp(%lfx) 244”,exp(A),B)；//输出return0；}//EndofFile计算实例求一个形如丁的经验函数公式，使它能够和下列数据相拟合(例6.5)。X12345678y15.320.527.436.649.165.687.8117.6 245程序输入输出inputmvalue：8Nowinputthe(x_i»y_i)，i=O,...»7：115.3220.5327.4436.6549.1665.6787.88117.6Solve：p(x)=11.437069exp(0.291215x)本章练习计算证明题1.下表中的资料说明各种等级汽车在碰撞时严重损伤的比率，求此数据的最小平方近似直经。车型车重(磅)严重损伤(％)美国制豪华型48003.1美国制中级型37004.0美国制经济型34005.2美国制轻便型28006.4其它美国制轻便型19009.62.下表给出一组数据:巧-1.00-0.500.000.250.75必0.22000.8002.0002.50003.8000分别用一次、二次多项式拟合这些数据，并给出最佳平方误差。3.给出下列数据：X-3-2-124y14.38.34.78.322.7用最小二乘法求形如y=°+力，的经验公式。 2464.给出下列数据:-0.70-0.500.250.75X0.991.212.574.23用最小二乘法求形如A的经验公式。5.用最小二乘法求解矛盾方程组。’公+2x2=52X1+勺=6为+叼=4XI-2x2=1Xi+5弓=13.12再+刍=7.6X1+弓=5.1第7章数值微分和数值积分7.1数值微分7.1.1差商与数值微分当函数/（x）是以离散点列给出时，当函数的表达式过于复杂时，常用数值微分近似计算了（x）的导数/'（X）。在微积分中，导数表示函数在某点上的瞬时变化率，它是平均变化率的极限；在几何上可解释为曲线的斜率；在物理上可解释为物体变化的速率。以下是导数/‘（，）的三种定义形式：r（x）=lin/（x+"）-/（x）'、/jo（2,lltn/（x）-/（x-^）7ft=lim2%（7.1）在微积分中，用差商的极限定义导数；在数值计算中返璞归真，导数取用差商（平均变化率）作为其近似值。 247最简单的计算数值微分的方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。向前差商用向前差商（平均变化率）近似导数有：(7.2)其中/+”的位置在X。的前面，因此称为向前差商。同理可得向后差商、中心差商的定义。由泰勒展开*/■（%+月）=/■（而）+V'（通）+天x修）,演4fW%+人得向前差商的截断误差：及"（%）*+?-/（*=-**©=0（〃）n2向后差商用向后差商近似导数有：一—八八"）一〃演7）k（7.3）与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差：&（归«）/■-华训n=。㈤x0-h<^ 248用中心差商（平均变化率）近似导数有:f'M/（号+〃）-/（勺-%）2h(7.4)由泰勒展开“Xo+〃）=/（而）+/（而）+4（刍）我2层/(^0-^)=/(^0)~¥/(^0)+—/*(X0)~—/*(^2)2!5\得中心差商的截断误差:&（X）=〃%）-麴储）丁一㈤=与二©）+“切］2月12=4/汽⑤=0(层)X。-h3«Xo+h0差商的几何意义.㈤咽.…-狗一、微积分中的极限定义1°力，表示/（X）在x一而处切线的斜率，即图7」中八而+»-/（勺）直线产的斜率；差商〃表示过（X"/（X。））和（%+九/6。+〃））两点直线。的斜率，是一条过X。的割线。可见数值微分是用近似值内接弦的斜率代替准确值切线的斜率。图7.1微商与差商示意图 249例7.1给出下列数据，计算八002)J,(0.06)J,(0.10)/(0.08),X-0.020.040.06().80.10f(x)5.065.075.0655.055.055解：［"02)”(5.07-5.06)/(0.04-0.02)=0.5八0.06)(5.05-5.07)/(0.08-0.04)=-0.5(5.05-5.055)/(0.08-0.10)=0.25/*(0,08)*(/(0/0)(0.06))/(0.10-0,06)=18,75设定最佳步长在计算数值导数时，它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断误差，由原始值X的数值近似产生舍入误差。在差商计算中，从截断误差的逼近值的角度看，1用越小，则误差也越小：但是太小的I川会带来较大的舍入误差。怎样选择最佳步长，使截断误差与舍入误差之和最小呢？一般对计算导数的近似公式进行分析可得到误差的表示式，以中心差商为例，截断误差不超过而舍入误差可用量1估计(证明略)，其中e是函数乂的原始值的绝对误差限，总误差为川»,e—Af,+ 6he\he4一弧+―=—M3--=0当I6h)3h时，总误差达到最小值，即(*)可以看到用误差的表达式确定步长，难度较大，难以实际操作。 250%),叫h,\通常用事后估计方法选取步长〃，例如，记I二1为步长等于2的差商计算公式，给定误差Z)(A)-Z)|-| 251设孙i=0,1…/为10力］上的节点，给定（个，/«）），？=QL…潭，以（西,丁（砌为插值点构造插值多项式4S）,以46）的各阶导数近似/（X）的相应阶的导数，即“X）=4Q）=（X）/（百）2-0X/⑶=4（X）=E4（x）/a）2-0八勺）⑶/（初八04，…/1（7.5）误差项为:云（X）4者!>F例7.3给定5J&）"=°工2,并有向-a=W〃，计算/（/）/&）J“2）解：作过aJ&））j=°12的插值多项式:3空零工（而）+（钎.，：一初〃再），(x)=W(x)=+，-电)一等以_/+工_町)h+xF 252将X=百代入了'(X)得三点端点公式和三点中点公式：，(Xo)=L(-3/(x0)+4/®)-y®))2%/U)=^-(-/(x0)+/(x2))2A八必)=工(〃而)-"(x1)+3」(盯))2n利用泰勒(Taylor)展开进行比较和分析，可得三点公式的截断误差是。色今。类似地，可得到五点中点公式和五点端点公式：/U)=-L[/(x0-2A)-8/(Ao-A)+8/(x0+A)-/(x0+2A)]12〃+Q⑸©Je[x0-2凡为+2〃]/Vo)=±[/(而)+48/(。+/)-36/(x0+2A)12«+16/(%+3〃)-3/(瓦+4〃)]+5/⑸©，火区,勺+4町 253把离散点按大小排列成a=勺<々々=葭用加关系式构造插值点a=°」'2,…/的样条函数8（力：当x=玉则/'（石）*n，当xe（4为+1）时，可用S'（x）=/'（x）计算导数。7.2数值积分在微积分中用牛顿-莱布尼兹（Newton-Leibniz）公式计算连续函数,（，）的定积分:（/@）以=尸©）一尸3）但是，当被积函数是以点列（斗-0,1,2,…/的形式给出时，当被积函数/（X）的原函数尸（力八inx2dk难以得到时，例如J1sm'则无法用牛顿-莱布尼兹积分公式计算。有时当被积函数的原函数过于复杂时，也不宜套用积分公式计算积分，而应采用数值积分公式计算定积分。在微积分中，定积分是黎曼（Rimann）和的极限，它是分割小区间长度趋于零时的极限，即不）J在数值积分公式中，只能用有限项的和近似上面的极限，通常由函数在离散点函数值的线性组合形式给出。 254/（力=3（力右/。）=£力（西）记w，在本章中，用表示精确积分值，用，表示近似积分值，｛'J称为积分节点，/称为积分系数。确定40）中积分系数巧的过程就是构造数值积分公式的过程。怎样判断数值积分公式的效果？代数精度是衡量数值积分公式优劣的重要标准之一。定义7.1（代数精度）记【凡句上以（小/（毛）尸=°工2「-/为积分节点的数值积分公式为4。）=汽%/5）2-0若满足与（/）=id）-4（/）=0*=。工…，掰,而纥陵吗*°,峭4。）具有次阶代数精度。由此可知当4。）具有〃阶代数精度时，对任意的次阶多项式都有1（,）7.2.1插值型数值积分对给定的被积函数在口»】上的点列aJa））J=°J2,…,力作拉格朗日插值多项式45）,以j4(x)dx=((x)_/(Xj)dxfZ（x）dx,“|'f（x）dxR八'近似计算J。，'',即j/（x）dx*.［也⑸可.动a,=f4（x）dx4⑺=14（x）dx=记'小，则有由i-。数值积分误差，也就是对插值误差的积分值纥（/）=（凡（x）dx 255玛(/)=工力Xo,X1，…,々㈤「|(x-xjdx对一般的函数玛(/"°,但若“X)是一个不高于〃次的多项式，由于//“(力=0,而有EH。因此，附阶插值多项式型式的数值积分公式至少有附阶代数精度。例7.4建立[02上节点为%"0，为=05用=2的数值积分公式。ba；=fL(x)dx解：由'k'得旬=i，o(x)dx=f(x-0.5)(x-2)(0-05)(0-2)dx=a「R(x)dx=。(x-°)S-2)办=竺Jo1(05-0)(05-2)9尸°)-dx=2(2-0)(2-0.5)9得以数值积分公式：.1-~、AC/)=-[-3/(0)+16/(0.5)+5/(2)]7.2.1牛顿-柯特斯(Newton-Cote's)积分卜=b-a把积分区间口力】分成力等分，记步长为«，取等分点玉=a+%G=°」，…/)作为数值积分节点，构造拉格朗日插值多项式4(X),取O(x)dx“"x)dx由此得到的数值积分称为牛顿―柯特斯积分。下面可以看到，牛顿―柯特斯积分系数和积分节点以及积分区间无直接关系，系数固定而易于计算。梯形积分以(a，/(a))和9，/(切)为插值节点构造线性函数A"),有 256j/(x)dx=j。⑶以那么，(。(x)dx=f(70(x)/(而)+4(x)/(再))办a0=fl0(x)dx=f——―dx=—(b-a)=(b-a^Q11J。Jaa-b2ax=f'i(x)dx=(^-^-dx=-(6-a)=(6-a)cf)J。6b-a2C(D=1C(D=1提取公因子9-a)后，得到牛顿―柯特斯积分的组合系数：°2,12,它们已与积分区间没有任何关系了。B5一ajy(x)dx=亭"⑷+/(/>)]b-aT(J)= [/(a)+/⑼记2(7.6)称丁(/)为梯形积分公式。它的几何意义是用梯形面积近似代替积分值(图7.2)。y 257怎样确定梯形积分公式的代数精度?我们可以取/(X)=LX,V…验证取了。)=1时，有/(力=£dx=b-a耿力=一C/(a)++1)=〜即：4/)=")取/(x)=x时，有/(力=£粒=匕；7(力=?(/(。)+/。))=((“+5)=即：«力=文/)取〃x)=，时，有/(7)=£0dx=H?。⑷+/(6))=7(力得梯形求积公式具有一阶代数精度。 258梯形积分公式的误差：/(x)-AW+华(X-0)(3),a«i由2! 2594(x)=ff(x-a)(x-b)dx得Ja2!因为(x-a)(x-8)在［a,句上不变号，由积分中值定理得到梯形求积公式的截断误差:&(x)=f(x-a)(x-b)dx(7.7)=^2^®aY，a&叮&b辛普森(Simpson)积分对区间口向作二等分，记5=。用=(a+妨/2,盯=b。以(aJ(a)),((a+b)/2J((a+与/2)和(瓦,9))为插值节点构造二次插值函数&(x),那么，有(卬乃办=(⑶〃/)+4⑶/(々)+。5)/5))心=j，o(x)dx_?(x-(a+6)/2)(x-6),-IaxJtt(a-(a+6)/2)(a-6)=?(5-a)=0-a)*)6="(x)dx=-(b-a')=(b-a)c|2)Ja6=(4⑶心=23-a)=(5-a)cf)a6计算得到积分组合系数：6,6 260O(x)dx~/“)=scn=?,(a)+”用")62J 2615(/)称为辛普森或抛物线积分公式。它的几何意义是用过三点的抛物线面积近似代替积分的曲边面积(图7.3)。分别将“X)"1,冗产,小代入到/(/)和SC/)中，可以得到SC/)=+”(审］+〃引=/C/)6I2J表明辛普森公式对于次数不超过三次的多项式准确成立，,(/)具有三阶代数精度。因此可设一个三次多项式满足条件：鸟⑷=/(a),鸟⑶=/0)对等M啜牛扑/愕)计算得到误差为：J(x)-月(x)=/，(x-a)(x-b).a<^ 2624(力=/(7)-W)757⑻)+(1⑻-SS) 263/㈤=S（第=彳8⑷+4川明+啊=sc/）故辛普森求积公式的截断误差：纥(力=10)7㈤J4)»)(x-a)「-(…)=/I2J钎等）（一）办2880/“⑺,a 264令工=。+戊，为=°+%,代入上式得-1)(Z-J+1)(2—I—1)一,(£-%),,I-3ndtJo(7.10)c；")='0ft(t-1),1,(/-j+1)(/-j-1),1,(Z-n)dt这里称"ST)!%”为牛顿-柯特斯系数-（»）可见在取等距节点时，积分系数q与积分节点和积分区间无直接关系，只与插值的节点总数有关，而在例7.3中的积分系数是待定系数，这就简化了数值积分公式，而不必对每一组插值节点xi都要计算一组相应的积分系数ai。在公式（7.10）中取力=1,可算出梯形积分系数；取力=2,可算出辛普森积分系数。在表7.1中列出力从1到6的牛顿-柯特斯系数。从表中可以看出牛顿-柯特斯系数具有对称性。7.2.3求积公式的收敛性与稳定性也看/⑷]=f/⑶以/"）=5/（石）定义7.2若,则称求积公式i-。是收敛的。定义中的〃78包含了0W&T源"，通常都要求计算积分的求积公式是收敛的。 265/式力=£%/（再）//\g[稳定性是研究计算公式皿当j有误差％时，的误差是否增长，现设/a）4片，误差记为a=|/&）-£|。=0」，…M。定义7.3对任给只要a=1/（初一片m蜕=°」，…，编，就有匕（力4⑦心是稳定的。3力则称求积公式.[*€/）=>0j（不）\m=n1x定理7.1若求积公式皿的系数%>U°=U,L…,），则求积公式是稳定的。证明由于角司r/a）T）蜕=o，l…㈤,故有I4（力-4（舟1=1B^C/（^）-Z）]|<占2>i=讥右-a）2-02-0于是对“£>°'第只要a=i/a）一邓£,就有I4S-4⑦|«呀"）«£故求积公式是稳定的。7.3复化数值积分由插值的龙格现象可知，高阶牛顿-柯特斯积分不能保证等距数值积分系列的收敛性，同时可证（略）高阶牛顿-柯特斯积分的计算是不稳定的。因此，实际计算中常用低阶复化梯形等积分公式。7.3.1复化梯形积分 266把积分区间分割成若干小区间，在每个小区间［不为+J上用梯形积分公式，再将这些小区间上的数值a积分累加起来，称为复化梯形公式。复化梯形公式用若干个小梯形面积逼近积分比用一个大梯形公式效果显然更好，如图7.4所示。这种作法使我们想起定积分定义，即它为被积函数无限分割的代数和。这也正是计算定积分最朴素的算法。图7.4复化梯形公式积分视图复化梯形积分计算公式h—b-a对［a,切作等距分割，有n，毛=a+淞i=O,l于是/（7）=『/（x）dx=三广/⑶办J1Mrxxi（初+/（%））-/"阊工在L不加J上，Jx,212,则有/（力=£（己"区）+/（演+1））-/”修）有,2-0乙"「1I11力31/2-1Ni-0"记〃等分的复化梯形公式为4（/）或丁㈤，有 2677⑶=看⑺=我卜/3)+2>3+谕+艰电122-12复化梯形公式截断误差由12j,根据均值定理，当/eC［a,3时，存在4厘。,即，有£?■"仁)=4©髭。，于是纥。）=-*"©=—0二⑶，屋"8由此看到复化梯形公式的截断误差按照/或者/的速度下降，事实上，可以证明，只要/（x）在（口力）b上有界并黎曼可积，当分点无限增多时.，复化梯形公式收敛到积分'（,）=工‘"）"'。记此飞剑⑻,则有对于任给的误差控制小量£>°,有就有IMC/）I〈£,式中［,］表示取其最大整数。7.3.1复化辛普森积分 268把积分区间分成偶数等分2加，记花=2冽，其中力+1是节点总数，加是积分子区间的总数。记«,Xj=a+%,i=O,l「/,在每个区间【叼“孙+2]上用辛普森数值积分公式计算，则得到复化辛普森公式，记为用（/）。/(力=f」(x)dx复化辛普森积分计算公式广”/（x）dx=乡0（0）+4/（孙Q+/（亏+1）-黑/”（约而加，62880,称M-l2卜用（力=+4/（0+1）+/D）2.0°1「JR—]M-1、=7/⑷+4\/（g+1）+22〃0）+/⑶〃百口1（7.13）为复化辛普森积分公式，它是/（X）在【叼”芯2“21上采用辛普森积分公式叠加而得。下面用图7.5显示2h复化辛普森积分计算公式中节点与系数的关系，取"=8,在每个积分区间上提出因子6后，三个节点的系数分别是1,4,1；将4个积分区间的系数按节点的位置累加，可以清楚地看到，首尾节点的系数是1,奇数点的系数是4,偶数点的系数是2。图7.5复化辛普森积分系数 269复化辛普森公式的截断误差 270设在［与，0+21上的误差为ZooU因此,3-s*C/)=-暮£*©)=-翳^ooU2ZooU~(b-a)52880“即卬…卷却C(7.14)与复化梯形公式类似，误差的截断误差按照/或者小的速度下降。可以证明，只要了（X）在3»）上有界并黎曼可积，当分点无限增多时，复化辛普森公式收敛到积分记圾=跚P"创，则有|纥（力壮^4小。目2880w）Jg-a)5监N2880c对任给的误差控制小量£>°,只要(b-a)5-——<£2880/就有国SK一 271/（力=je*dx例7.5求8，i博中要求有5位有效数字。用复化梯形和复化辛普森求积公式的分点应取多少？解：/（X）=/J"（x）=/"（x）=/irwi=i/（4）wi<^o 272=/（/）+?%）_a+b这里，22。可以看到，芍"）的值是看（/）与新增分点,（勺）的组合。取分点%=4,计算得Z（/）=一（与+午+/（药）+/（今）+/33））=9+?（/（再）+/53））西=<3+与），马=：8+8）这里，22。同理，计算方"）时只要在40）的基础上计算新增分点/。1）,了（勺）的值再做组合，如图7.6所示。图7.6石⑺与£⑺一般地，每次总对前一次的小区间分半，分点加密一倍，并可充分利用老分点上的函数值，每次只需计算新增分点的和。对［。囚上阀等分，“«,则有 273［1X-1］、看⑺=叫打（。）+£/（再）+力（与（222记［&&+J上的中点为才|+¥则L（1»-1M-11/（力吟-f（A）+£/⑷+£/&#+r@）NIN?-1i-0乙吟,八。）+£小）+9（斗?£〃如）/1/2-14）2-0&C7）=夕4（力+兄（力）凡（力=比%）其中皿纭=£+也汇/("(21)%)乙i.1我2n=二―其中2«。类似地，可得积分节点为“，2片的辛普森求积公式的关系式:$2<力=:凡（力+:（4即（7）-必⑶）玛”。）=d£u&h）+〃。，］2>0由误差公式:1（力-看（力=-与齐r©/（力工⑺ 274I»-112»-1/"©=—£/"4），/々）==，/"8）c由于«»-0m3，分别为%及2%个点上的均值，可视二O八*于是上式表明看（力的误差大约是乙°,要|7"）一心⑺K%只需心SYGOK3£即可,而用I&C/）-1（力I作为控制手段简单直接，序列北（力名式力，…在计算机上也不难实现。复化积分的算法描述从数值积分的误差公式可以看到，截断误差随分点%的增长而减少，控制计算的精度也就是确定分点数万。在计算中不用数值积分的误差公式确定分点数力的理论模式，而用IZ*。1）-北（力1〈3£作为控制,通过增加分点自动满足精度的方法称为数值积分公式的自动积分法。即在计算中构造序列《，石门4』…，凡一T双I（6直到|&<-41<£或因m时停止计算，由分点数自动控制积分值的误差，并取（力。下面描述复化数值积分公式的自动控制误差算法，详细程序和算例请看本章7.6节。1.输入：误差控制精度e=eps；初始分点值。2.计算力分点的复化梯形积分4，72=4 275Tl=T2+100〃迭代计算中T1和T2分别表示*和3.whilelTl-T2l>eT1=T2%(力=底/(3H=Hh〃计算新增节点的值MT2=(Tl+H)/2H=h/2,n=2n〃将区间一分为二endwhile4.输出积分值T2。在自动控制误差算法中初始分点值不宜过小，以防假收敛。7.3.4龙贝格(Romberg)积分由前面得到的关系式(7.17),可以将3作为的修正值补充到/即141W〜0(力+京4(7)YS)333(7.18)其结果是将梯形求积公式组合成辛普森求积公式，截断误差由。("2)提高到°(〃’)。这种手段称为外推算法。外推算法在不增加计算量的前题下提高了误差的精度，是计算方法中一种常用手法。不妨对号式力，工GO再做一次组合。由"力”-鬻眇e唱’4(7)-S")=-技-a八加ZooU得到 276/（力心式力4（巴式力一凡（力）IG）嗓±86）=551515（7.19）复化辛普森公式组成复化柯特斯公式，其截断误差是。（炉）。同理对柯特斯公式进行组合：1（/）-“/）=/W-%（/）=U得到具有7次代数精度和截断误差是°（必）的龙贝格公式：641凡（力=右弓式力-白”力03OJ还可以继续对0（力做上去。为了便于在计算机上实现龙贝格算法，将方看酎《，鸟，…统一用&7表示，列标j=1,2,3,…分别表示梯形、辛普森、柯特斯积分，行标上表示分点数%=2",或步长&=人/2"）。龙贝格计算公式：对每一个口/从2做到亡，一直做到I&人-I小于给定控制精度停止计算。龙贝格算法 277龙贝格算法按表7.2元素的行序进行运算，用.1，舄.1，立22「,在计算中每个元素只用到上一行和本行的元素。对上面的算法进一步优化，对每A行可将计算定义在两行元素之间，令凡力舄中衣1，为鸟.，；在每计算一行元素后，要将舄'j-Mj，J=1,2,…天。表7.2龙贝格算法计算元素顺序表R口国1出%--•■《2%31.输入区间端点小，精度控制值e,循环次数定义函数,(x),取N=L〃=&-a；2,%=6a)+/@)”2；3.fork=2to河R-+初+⑵-W)/2//hk=hf2^2-1forj=2tok{&J=&J4+(%』-%5)/(4*-D退出循环4.输出7.4重积分计算在微积分中计算二重积分是用化为累次积分的方法进行的。计算二重数值积分也是计算累次数值积分的过程。为了简化问题，我们仅讨论矩形域上的二重积分。有很多非矩形域上的二重积分可作变换将其转换到矩形域上。 278(7.20)其中：a,b,%d是常数，在。上连续。像在微积分中一样，将二重积分化为累次积分:£0（"）力"口O（xj）力0（721）二重积分的复化梯形公式对区间［名句和k，d］分别选取正整数礴Dn,在X轴和V轴上分别有步,b-a,d-ch=rk=mnf八x,向dy用复经梯形公式计算1'八；计算中将x当作常数，有F］1超7、（f^,y）dy"k5g>）+*“*）+沙科）出IZJ（7.22）再将刀当作常数，在x方向上计算式（7.23）中每一项的积分，有［力h（11次-15jy（“。）杰"万15〃两‘&）+y区'%卬1ahf11w-i'“了（““灿〜彳/（为/j+rk，乂）+2八/既）2a/1/乙2-1；乃）dx“塔+*（aM+*（3,））x-i/1i、=5与号，M）+5/（勺，乂）+AWa>11//7-1i-1"""）叱〜妓｛；（/（飞Jo）+/（x*Jo）+/（与/）+/（/,〃）） 279/M-l、+—|W（'”%）£/（4打）>/（,乃）2111+2-1+j-i+J"j?一1।mx+2X/（m4）｝处E>3（纣j）i-4j-lJ=i_oj-0积分区域的4个角点的系数是1/4,4个边界的系数是1/2,内部节点的系数是1。误差：四一*陛3唱词（%〃）和（万,"）在积分区间内。例子7.6用复化梯形公式计算二重积分JJs'"'+"）力公，取入=上=0.25。解：/（XJ）如表7.3所示：表7.3/（XJ）数值表X1.001.251.501.752.000.000.8414710.9489850.9974950.9839860.9092970.250.8735750.9668270.9999660.9709320.881530.500.9489850.9974950.9839860.9092970.7780730.750.9999660.9709320.881530.7373190.5472651.000.9092970.7780730.5984720.3816610.14112口的数值如表7.4所示:表7.4l-J数值表X0123401/41/21/21/21/411/21111/221/21111/231/21111/241/41/21/21/21A4 280.jjsin(x2+y)dydx=0.25x0.25-(0.841471+0.909297+0.909297+0.14112)+-(0.948985+0.997495+0.983986)+g(0.778073+0.598472+0.381661)+g(0.873575+0948985+0.999966)+A(0.88153+0.778073+0.547265)+(0.966827+0.999966+0.970932+0.997495+0.983986+0.909297+0,970932+0.88153+0.737319))=0.873601f1f2sin(x2+y)dydx“一JoJi的准确值是0.886176。二重复化辛普森求积公式对区间卜力］和［*”］分别加等分和力等分，在x轴和丁轴上分别有步长,b-a.d-cn=,k=,mfnm阀均为偶数类似于二重复化梯形公式推导，得>jKL.MM[[力dx=勺 281记Z7='=0,424,…,241产，={5a=(1424,・24,1产卬7=%叼误差：ECO「S-雷-叩斗3炉歌伉GE"）和（万M在积分区间内。按例7.6的化分区间，%的值如表7.5所示:表7.5X0123401424114168164228482341681644142417.5,高斯（Gauss）型积分公式介绍对插值型积分公式在牛顿・柯特斯积分公式中要求节点是等距的，其优点是计算积分系数的公式规则相同，缺点是制约了求积公式的代数精度。可以证明：当节点个数附为偶数时，求积公式具有N-1阶的代数精度：当节点个数«为奇数时，求积公式具有力阶的代数精度。如果我们不预先指定求积节点大的位置，升和权系数4都作为待定的常数，能否适当地确定它们， 282以提高积分公式的代数精度？回答是肯定的。2%个待定常参数，需要2〃个方程来确定，取一个函数组：Qx,/,…,这一组函数构成了2都-1次多项式的基，任一小于等于2"-1次的多项式，都可以用这组函数的线性组合来表示。如果某一积分公式，对这组函数都能精确积分，则此积分公式就有次代数精度。八cCXXf」(x)dx=。"(%)+。21/(工2)例7.7计算求积系数402和求积节点而，々，便导J-17''U"2J'T至少具有3阶代数精度。解：按照求积公式的代数精度定义，分别令/。)=1，冗,，了,得方程组：解方程组得：x=-4=4==0,5773503q=l,c2=1,的g求积公式：jj(x)dx=/(-0.57735)+/(0.57735)按例7.7的方式，构造更高阶的代数精度的求积公式，生成求积系数和求积节点的方程组并无困难，而求解该方程组则无一定的章法可循。一般地，通过计算正交多项式的零点作为求积节点。当取积分节点为正交多项式的零点时，则节点个数是管的求积公式具有2花一1阶的代数精度。并称积分节点为正交多项式的零点的数值积分公式为高斯型积分公式。为了一般性，考虑积分其中印。）（印⑶20）称为权函数。当印（x）=l时，即是普通的积分。 283对于不同的权函数取（X）选定的节点也不相同。如何构造高斯型积分公式呢？对给定的间及权函数次⑴，由施密特（Schmidt）正交化过程作出正交多项式4（x）/（x）,…,々（X）；解出正交多项式片（X）的找个零点加电,…,这个片个零点就是积分节点；以这些节点构造插值多项式，计算积分系数1,2,…，盟其中是拉格朗日插值基函数。高斯型求积公式为Jtic/)=2Xz"㈤i-l高斯型积分公式的优点是它的代数精度高，特别是对无穷区间或瑕积分更有效，但计算正交多项式的零点即积分节点有一定的工作量，好在数值学家们已算出一些特定的函数的积分节点和积分系数，在计算中我们可以查表直接得到这些数值。本章并不构造各种高斯型积分公式，有关的详细内容请参考有关的教材。下面给出口力卜17』上，取权函数的的高斯型积分。取［-1,1］上权函数/0）=1的正交多项式为勒让德（Legend优）多项式:14舅4(x)=rttK/-1力，Zn\ax高斯-勒让德〃=2,4,5相应的积分节点和积分系数表如下:n相2々=-0.5773503,弓=0.57735031.0000000,1.0000000 2844%1=-0.8611363,叼=一().3399810=0.3399810,々=0.86113630.3478548,0.65214520.6521452,0.34785485-^1=%5=0.9061798一式2="=0.538469勺=0.00.2369269,0.23692690.4786287,0.47862870.5688889b_a^b-a土要计算一般区间〔&可上的积分《只需作变量代换*22则有fy(x)dx=fg(t)dtg«)=/-r-+—T-J/(/)=r，(x)dxk2J“，其中I2二J,这样5L八J仍可用高斯积分求积，即b-ab-a[(力==-G*(g)= XaiS(zi)1=fcosudu例7.8应用两点高斯・勒让德积分公式计算J-io解：/=(-0.5773503ycos(-0.5773503)+(0.5773503)2cos(0.5773503)=0.558608I=jdx例7.9应用两点高斯-勒让德积分公式计算出1+一。解： 28511X=一+-Z令22,得到积分i=2—，、2=2(/1⑷+/③))=0786885口4+«+1>例7.10证明不存在氨，漏住二12…㈤使求积公式*WaJ(4)JI的代数精度超过2%一1次。证明：只要能找到一个2当次多项式，使求积公式两边不相等即可。用反证法，假定存在求积系数和节点燃及敌色=L2,…/)使求积公式对任何2%次多项式/(x)精确成立。现取」(x)=确》=n(x-xj■代入求积公式左端得Jw(K)/(x)dx=£w(x)4(x)dx>0X靛£福了(凝)=显®)=0而公式右端i5,故右端与左端不相等，与假设矛盾，说明不存在2〃次代数精度的求积公式，故高斯型求积公式是具有最高代数精度的求积公式。 2867.6程序示例(1»-111/⑷+£/(&+论)+“㈤2I21的自动控制误差算法。算法描述输入6,a的值，定义」(x)；n=m9h=(b-a)/〃T2=T/k++/)+：/“)1272人Tl=T2+100whileIT1-T2I>£T1=T2；(3由2)H=〃i-。!或用复化梯形公式计算T2T2=(T1+H)/2；h=h/2；n=2n；endwhile输出T2。程序源码lllllllllllllllllllllllllllllllllllllllllllllll//Purpose：梯形公式的自动控制误差算法〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 287#definea1.0#defineb2.0#defineepsilon0.00001voidmain(){intn=m；inti；doubleTn,H_n,T1,T2；doubleh=(b-a)/n；T_n=(f(a)+f(b))/2；for(i=1；i 288printfCT=%lf 289?,,T2)；//EndofFile计算实例对于/(x)=Sin(x),在区间［aj］上验证梯形公式的自动控制误差公式。程序输入输出对于不同/(X),修改程序#definef(x)项，本例/⑸=沏⑶，初始次=100,区间［a,占卜口,2］。)T=0.956447程序7.2龙贝格积分算法计算公式和算法描述第7.3.4节所述。程序源码IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII〃Purpose：龙贝格算法〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include 290sum+=f(aa+ih)；sum+=(f(aa)+f(bb))/2；*return(hsum)；)voidmain()(inti；longintn=N_H,m=0；doubleT[2][MAXREPT+1]；T[l][0]=computeT(a,b,n)；n=2；for(m=1；m 291M,T[l][m])；〃输出return；}printf("Returnnosolved... 292M)； 293//EndofFile计算实例1fsin(x)rfx利用龙贝格积分法计算.o程序输入输出初始(对于不同"X)，修改程序#define/")项，本例/⑶=皿⑶，区间［a,力卜［1,2］,h=(b-a,)fn=(b-a)/20)对于本程序的给定，输出结果。T=0.956447本章练习计算证明题1.设函数/(X)由下表给出，分别用向前、向后差商计算了'(°°2),f(0.°6)。X().0()0.020.040.06“X)2.001.931.881.822.设函数/(X)的数值由下表给出，用向后差商计算/'(020),f(040)。x0.100.200.300.400.050.100.150.203.*设函数,(x)的数值由下表给出:X“X)X“X)0.060.871475().110.8929950.070.8760220.120.896943 2940.080.8804460.130.900780.090.8847480.140.9045060.100.888930.150.90812取不同的步长用向前差商计算/,观察误差变化规律，从而确定最佳步长。1.写出下列两种矩形公式的余项：⑴。⑶dx~/(a)—a)/(a+f—(6-a)(2)\Z/5.构造积分‘5=L/(x)dx的数值积分公式。/C/)-«V(-A)+V(0)+V(2A)6.用梯形公式计算下列积分：⑴1.2⑵Jsin(cosx)ax7.用辛普森公式计算下列积分: 295f/(x)dx8.设函数由下表给出，分别用复化梯形和复化辛普森公式计算出6o9.六'取£=10',%=1,试用龙贝格公式计算积分直到1段厂&但1〈£时停止，并作出龙贝格积分数值表。10.用复化梯形公式计算二重积分：(1)U-^dxdy,取力=»4RJ1ff——dxdy(2)x+y，取力=k=5ii.用具有三阶代数精度的高斯-勒让德积分公式计算下列数值积分值：fV1-cos2x+sinx2dx⑴06一。2dx(2)小第8章常微分方程数值解在描述系统的动态演变时，例如，物种的增长和蜕变，物体的运动，电路的振动瞬变，化学反应过程等，都能将其表示为以时间,为变量的常微分方程或方程组。物体冷却过程的数学模型可用下式表示：dt-k(u-u0) 296du它含有自变量入未知函数〃以及它的一阶导数由，是一个常微分方程。在微分方程中我们称只有一个自变量函数的微分方程为常微分方程，自变量函数个数为两个或两个以上的微分方程为偏微分方程。给定微分方程及其初始条件，称为初值问题；给定微分方程及其边界条件，称为边值问题。本章主要讨论常微分方程的初值问题：,(a 297o用数值微商的方法，即用差商近似微商数值求解常微分方程。对于常微分方程初值问题[式(8.1)],在求解区间口，句上作等距分割的剖分，步长加，记x*=xnA+h,n=\,2,-,m用向前差商近似y（x）做出y（X）的在*=%处的一阶向前差商式：n又＞'（工0）=」（而丁（工0））,于是得到h而丁（々）的近似值必可按必;°=，%）或必=为+47（飞%）求得。类似地，由以及y（xj=/（“（/））得到计算近似值居一的向前欧拉公式：居+「％+%/（X*/）（82）由差商（差分）得到的方程（8.2）称为差分方程。由其直接算出乂+1值的计算格式称为显式格式，向前欧拉公式是显式格式。 298欧拉方法的几何意义以/（%,%）为斜率，通过点（X。Jo）做一条直线，它与直线*的交点就是必。依此类推，%+1是以/（X*，乂）为斜率过点（勺，居）的直线与直线x=x*+i的交点。欧拉法也称为欧拉折线法，如图&1所示。 299图8.1欧拉折线法例8.1假定某公司的净资产因资产本身产生了利息而以4%的年利率增长，同时，该公司以每年100万的数额支付职工工资。净资产的微分方程：dwdt=0.04w-100（f以年为单位）分别以初始值忒°）=1500万，v（°）=2500万，2（°）=3500万，用欧拉公式预测公司24年后的净资产趋势。解.qS+4（0.04吗-100）=1.04%一100,〃=1w0分别以%=1500,%=2500,Z。=3500代入，计算结果见表表8.1计算结果n乂4打X*乂4]14602500.354051283.352500.3716.6521418.42500.3581.661234.682500.3765.3231375.142500.3624.8671184.072500.3815.9341330.142500.3669.8681131.432500.3868.5791076.692500.3923.3117552.12500.4449.9 300101019.762500.3980.2418474.1832500.4525.8211960.5462500.4039.4519393.1512500.4606.8512898.9682500.4101.032()308.8772500.4691.1213834.9262500.4165.0721221.2322500.4778.7714768.3242500.4231.6822130.0812500.4869.9215699.0562500.4300.942335.28452500.4964.7216627.0192500.4372.982463.30422500.5063.3从表8.1可以看到当利息赢利低于工资的支出，公司的净资产逐年减少，以至净资产为负值；当利息赢利与工资的支出平衡时，公司的净资产每年保持不变：当利息赢利超过工资的支出，公司的净资产稳步增长。在图8.2中ZL£2,Z3分别表示初始值3500,2500和1500的三条净资产趋势曲线。L1.（3500）一£2.（2500）25图8.2三种初始值的净资产趋势用向后差商近似y'（X）做出丁（力在*=0处的一阶向后差商式：h又丁’（为）=，＞（西）），得到，（X1）的近似值乃的计算公式： 301乃=为+”（々仍）类似地，可得到计算丁（々+1）近似值居+1的计算公式:%+1=丁*+"（々+1，%+1）（8.3）公式（8.3）称为向后欧拉公式。通常产（XJ）为丁的非线性函数，因此式（8.3）是关于乂+1的非线性方程，称为隐式欧拉公式，需要通过迭代法求得居+1。其中，初始值乂方可由向前欧拉公式提供。从算法结构上看，显式公式比隐式公式简单：从方法的稳定性和精度上看，多数情况下隐式公式优于显式公式。最简单的迭代公式（皮卡（Picard）迭代）为：‘乂2=居+汇8.乂）I淄）=乂+"（3斓）’"'''…直到此一乂胃"给定精度。可以证明，人充分小时，以上迭代收敛。事实上，记"3）=居+"（/+1，＞）,则”8）=彷（4＞）〃充分小时，可以保证限（々Md股〈1,其中上为李普希兹（Lipschitz）常数。用中心差商近似y（X）做出A（x）的在*处的中心差商式： 302又y（xi）=/（为，（再）），可得到丁（马）的近似值乃的计算公式：乃=K+2”（々M类似地，可得到计算t（/+1）近似值居一的计算公式:(8.4)%+i=匕-1+2"公式（8.4）称为中心公式。按公式（8.4）,需要知道其4居的值才能求得居+1的值。因此，要先用其它公式算出乂，再用中心格式算出乃，丁3,…。需要指出的是，中心格式不是稳定格式，因此不能采用。8.1.1*欧拉公式的收敛性定义8.1（局部截断误差）在假设%=丁（々），即〃步是精确的前提下，考虑的截断误差方+i=y（/Q-%+i称为局部截断误差。对y（x*+i）在您做泰勒展开：,我27（xs+1）+&）可（勺）+如（Q+力丁（幻小"女"田（&5）欧拉公式（8.2）是由以上展开式中截断2!而得，故欧拉公式的局部截断误差为：2（8.6）若y（x）为线性函数，则丁“8）三°,这时局部截断误差为o,由欧拉公式得到的解为精确解，故欧拉公式是一阶方法。 303图8.3Euler公式的误差定义8.2如果给定方法的局部截断误差是看+1=。(3“)则称该方法是P阶的，或称具有P阶精度。以下考查欧拉向后公式的局部截断误差:Mx"-y(x")-汇(％•y@*+〃)-双)-如a+h)=加'心)+)+…-a)+如"(x*)+T-日丁气）+。尚）层故此方法的局部截断误差主项为一万,0口,也是一阶方法。整体截断误差和收敛性在计算乂+i的局部截断误差时，是假定在、的y（x）值是准确的前提下，即y（x*）=y*。实际上，从计算必开始，每个片估2,…/）都会产生截断误差，而且在先的误差会扩散到加1 304中，将这些 305前列点的误差累计到计算y（o+i）的误差中，称为整体截断误差。它将影响到方法的收敛性，我们将估计这一误差。由微分方程理论，为保证微分方程解的存在惟一性及稳定性，通常对丁应满足李普希兹条件，即对任意，乂了，存在满足于是由式（8.5）与（8.2）两式相减得到•d）一居“I?=＞（5）-匕+河/（x*+—1/©）*曾记2!,故有I-1«|4|+〃|」区j（Q）-/（/,乂）|+|与“国4|+近㈤|+区“|?=01产%|=。尚)记*,则有|%“|«(1+以)%|+7因此有反邛(1+例(1+8|=|+7]+--4(l+ZA)-i4|+T+(l+ZA)7+(l+ZA)a7,+--+(1+ZA)"T=(1+息严|%|」-(1+4产〉1-(1+ZA) 306〈Q+人力）\e0\+1Ln二（心啡吃）对z>0,由公式（l+z>«e”,最终可得：一仙+色“1）仙+叶卜Lh}\Lh}ECT/g）JIe-其中“屋。I由原始误差引起，当初始值为精确值时，这一项的值是0：£4仍是由截断引_2/…工70起，由于7二°（/），当〃7°时，Lh,即欧拉公式是收敛的。8.1.1基于数值积分的差分方法对常微分方程（8.1）即五两边在区间［x»/+J上积分得：双+1）-丁区）=广/（"）dx即=y（xj+『/（xJ）dx=y（Q+广外）山用矩形近似积分公式计算上取—（£）*/（々）=/（x,1y区）），有 307。y'it}dt=（%-x*»（x*）=则得向前欧拉公式：%1=居+"（/，乂）(8.7)取y«）*y'Sz）=/（xs+iM^+i））,有「y'Sdt"x*+i-x*）y（x*G="（uD）J。则得向后欧拉公式：%%+1）（8.8）用梯形近似积分公式计算上3r1ry⑥力〜F（y（%）+/（/））JX，ZhJ（/））+/（x*+iJ（x*+1））］则得到梯形公式：hy^i=y»+-办）+〃x*+i，/J）2（8.9）梯形公式也是隐式格式，可用皮卡迭代或牛顿迭代计算（8.9）的居+】。计算中为了保证一定的精度，又避免用迭代过程不菲的计算量，一般先用显式公式算出初始值，再用隐式公式进行一次修正。称为预估-校正过程。例如，下面是用显式的欧拉公式和隐式的梯形公式给出的预估-校正公式：{%+1=%+”"”,居）h_"+1=%+Jn）+/（。+1，区+1）］,（8.10）式（8.10）也称为改进的欧拉公式，它可合并写成h乂+1=片+/d，乂）+/（x*+i,居+¥（^»，%））） 308例8.2分别用梯形公式和预估-校正公式解初值问题：Idx,0.04x40.4卜（0）=1 309解:y0=1,A=0.1用下面的迭代公式，对每个迭代4次，上=L2,3,4y»1y°--该方程的精确是1-X,计算结果如表8.2所示。表&2计算结果nX*y»双）|居一人）110.11.11181.11110.000720.21.25201.25000.00203().31.43311.42360.009540.41.67631.66670.00048.1龙格-库塔方法8.1.1二阶龙格-库塔方法常微分方程初值问题:，'⑶=/(x，7)[(a)=乂>做y（x+〃）在x点的泰勒展开:y(.x+&)=y(x)+hy\x)+3/(x)+ 310+(x+网=丁⑶+hy'(x)++…+，⑺⑶+T2!p\这里”641,九。( )。取XF,就有>2京y(x*+D=>(々)+»(5)+7Tl/⑶+…+中’⑷+乙】2!P'-(8截断4+1可得到y(x*+D近似值居+1的计算公式，即欧拉公式:%+1=%+”(々，居)若取P=2,式(8.11)可写成：y(x*+i)=+»(/)++了=y(x*)+纭(x*j(Q)h2+5"(々J(xj)+力(X*J(X*)M(x*J1))]+4+1或双+i)=双)+处/61Mxit))}(8.12)+【〃”(/))+/(x*,(/))/(/j(x*))]+Zu 311*“什©，氏区公］5\截断4+i可得到y®+i)近似值%+i的计算公式:乂+1=居+刎（々，招）+引［〃\，乂）+力J*）］或居+1=〃+我/（""）+罚〃+力（""）/（"»）］）上式为二阶方法，一般优于一阶的欧拉公式（8.2）,但是在计算居+1时，需要计算工人，力在（々）J点的值，因此，此法不可取。龙格-库塔设想用了（XJ）在点（//（々））和（（五+办），》（々）+匕”（玉,丁（勺）））值的线性组合逼近式（8.12）的主体，即用力」（5，丁区））+七/（々+叫1y（Q+8”（x*j（4）））（8」3）逼近/（乙,〉（/））+/（/，双））+力（“（/））,（//（勺））］得到数值公式：%+1=%+彬1/（3（5））+cj（/+M-+W（x*，Mx"）））］（8.14）或更一般地写成 312‘%+i=乂+〃9由+。23）<左=/（%,〃）夕+a〃，乂+M占）对式（8.13）在（4J（/））点泰勒展开得到：cj（x*J（X*））J（X*）+B优(X*j(x*))+缎〃々,7(々))) 313(cl+c2)/(x*,y(x*))+32cMr(/MQ)+勿力对人工",y(x*))1/(x*j(x*))+。(层)]将上式与式（8.12）比较，知当。卜勺，。力满足q+C?=12c^a=12c2b-11q=一若取2时有最好的逼近效果，此时式（8.13）一式（8.14）=0这是4个未知数的3个方程，显然方程组有无数组解。c2=-,a=U=l2,则有二阶龙格-库塔公式，也称为改进欧拉公式:乂+1=%+g(左+多)左=/("”)*2=/(/+〃,居+磕)(8.15)q=0,C2=l,a」/」若取22,则得另一种形式的二阶龙格-库塔公式，也称中点公式:乂+i=y*+g,“hh,s(8.16)e=/5*+5,乂+万左) 314从公式建立过程中可看到，二阶龙格-库塔公式的局部截断误差仍为。（&3）,是二阶精度的计算公式。类似地，可建立高阶的龙格-库塔公式，同时可知四阶龙格-库塔公式的局部截断误差为OS3,是四阶精度的计算公式。欧拉法是低精度的方法，适合于方程的解或其导数有间断的情况以及精度要求不高的情况，当解需要高精度时，必须用高阶的龙格-库塔等方法。四阶龙格-库塔方法应用面较广，具有自动起步和便于改变步长的优点，但计算量比一般方法略大。为了保证方法的收敛性，有时需要步长取得较小，因此，不适于解病态方程。8.2.2四阶龙格-库塔公式下面列出常用的三阶、四阶龙格-库塔计算公式。三阶龙格-库塔公式h乂+1（总+4月+内）占匆）(8.17)方%+2妩2）乂+1=乂+（（用+3总）左=/1办）占=/区+（〃,稣+3方左）⑵卜"产）(8.18)八+i=Z.+J(2&+3芍+%)左=/(2*),11e=/(/+5〃,〃+万咕)33%=/（勺+彳氏乂+了岫）(8.19)四阶龙格-库塔公式 315h,,c,2,、%=%+工（左+2&+2&+匕）o左=«为=/（x*+，〃，%片）（8.20）%=/（々+?，乂+［桃）K=」（/+机乂+蚂）乂+1=乂+/+%+3质+々）O左一1/（',%）<用=/（X*+g〃,/+；，左）用=/（/+|〃/+，左+%）K=/（/+队匕+咕-%j+/）例8.3用四阶龙格・库塔公式（8.20）解初值问题:dy204x40.8—=ycosx 316计算结果列表&3中。表8.3计算结果>2X*y»双）1乂-双）110.21.247891.247920.0000320.41.637621.637780.0001630.62.296182.296960.0007840.83.533893.538020.004138.2.3步长的自适应欧拉方法和龙格-库塔方法在计算居+1时仅用到前一步乂的值，我们称这样的方法为单步法。在单步法计算中根据需要可以取等步长或变步长。对于y。）变化平缓的区段，步长可以取大一些；而对于y（x）变化剧烈的区段，则步长可以取小一些。以计算为必为例，已知外的数值，取"="。为初始步长，£是给定的控制精度，用欧拉或龙格-库塔方法的计算公式。记用,为由必和取步长〃算出的乃，记乂,"为由必和取步长〃/2算出的必，即乂"，泮⑶均是取+〃）的近似值，当［1#）-乂""|<£,逐步放大步长，22A,直到时为止,最后取〃/2为步长的值。而当1乂"-乂""|<£,逐步缩小步长，22A,直到旧£时为止，取〃为步长的值。下面给出自适应算法的描述：给定控制精度£和初始步长对于力=1,2,…，加，城-乂"邓£THEN 317WHILEENDWHILEh=2fhELSEWHILEBh=21hENDWHILEh=hENDIF8.3线性多步法常微分方程初值问题［式（8.1）］与积分'"）"丁"）+工‘^"）"等价，在&1.3节中用一些简单的数值积分公式建立了欧拉和梯形等数值方法，下面给出一些更一般的用数值积分工具求解常微分方程初值问题的方法。数值积分公式当积分节点中包含时称为隐式公式，否则称为显式公式。分别取为/为剖分点/+1'乂-，，有7（^+i）=+广y（x）dx若用积分节点…'/r构造插值多项式近似V（x）,在区间上计算数值积分J、"，则称为隐式公式。与通常数值积分不同的是，在线性多步法公式中，有两个控制量「和北尸控制积分区间，q控制插值节点。特别取P=°,有vD=双）+0y（x）dx 318若积分03用关于积分节点…的数值积分近似，就可得到显式的阿达姆斯公式。例8.4构造0=°，，=1的显式公式:双+1)=y(xj+0及⑶火/)+R(x)]dx这里广*4(x)dx7MX**Tdx=-hdx=--h2J%J%x„2心=广R5dx-p22>--)(x-j)dx=⑷J，j/NINy(%i)=y(%)+与/❷,居)-/(、_"(x.i))]+4+1即2截断1+i可得到y(、+i)近似值居+i的计算公式：h稣+i=〃+A[31/(“0-/(4_1,居_1)]2(8.22)上式称为二阶显式阿达姆斯公式。类似可得三阶显式阿达姆斯公式：h乂+广以+丘[23〃“*)_i6X*，k)+5RMm)](&23)O四阶显式阿达姆斯公式:h 319乂+i=外+力55」（5，乂）-59/（;—）+37」（x—Ji）-9」（x*-3,加）]（8.24）“系仆）「y（x）dxXX…X若积分k用关于积分节点，A】'*''*+1的数值积分近似，就可得到隐式的阿达姆斯公式。隐式格式需要通过迭代求解。它们分别为二阶隐式阿达姆斯公式（也称为梯形公式）：h%+i=居+三【/（、,乂）+/（/+1，%+D]2（8.25）彳+L（中⑴⑶三阶隐式阿达姆斯公式：hy»+1—蔚+yr[5/（\+卜居+i）+x）/（与一1'y»-i）l12（8.26）四阶隐式阿达姆斯公式：h乂+i=〃+汇W（x*+i，乂+i）+19/（x*,乂）15」（工1J1）+/（/-2J»-2）1（g27）和=啕"丁'© 320例8.5建立°=Lq=2显式公式。解:用+「匕1+。口式x)」(x*M)+4(x)/(xLi)+，2(X)/(/_2，*.2)]dx居+1=%-1+如%/(X*必)+//(入1,/_1)+。2/(*2办-2)],(X-、_1)0-/_，,7aQh=\-^-dx=-hJx(x*-*])(、-k2)3,小《(x-x_)(x-x__2).2,印=f--——dx=—-h%(/-「々)(*1-*2)32=广(x-Q(x-x)dx」〃JJ(X”2-X”)(4-2-X*T)3计算格式：hy^i-A-i+-Uf^,yx)_2/(%J*-)+/(xx_2,7,.2)]n=2,3,4,---,w-l,A=---其中：w.截断误差：7产*wws,4+1=九―h(x-%2)(X-%1)(x-xjdx=!亡/⑶©(x-x”2)(x-5-1)(x-x*)dx6JiJ 321(/⑶(")=/%))这是一个三步三阶的显式格式，可用龙格-库塔公式计算必'为。例8.6构造「=2,g=2显式公式。解：取〔'-2，々+J为积分区间,以5+1J（%）），&J（4）），（%，/（kJ）构造拉格朗日插值多项式。[%+1=片-2+跳凤/（勺+1，%+1）+/V（X*/）+为（3，%-1）]眦=广”（x-x.x-ki）dx=lha,滔《（x-5）（x-x*）9=\——dx^-hJ-（k-x*+i）（x*_「/）4h%+1=〃-2+;[3/（%,乂+1）+9」（X1，兀1）]计算格式：4截断误差：1X看+1=7匚+y⑴（⑤（x-%）（x-5）。-xQdx0J3=-），*©o为了避免迭代，可用预估-校正公式替换上面的计算公式。丸+1=%-2+1（7/（^，乂）-2/（x*_i，1yI）+/（々-2,乂-2）[h- 322%+1=〃-2+/3/（々+1又1）+9/（*1，居-1）]8.4常微分方程组的数值解法8.4.1一阶常微分方程组的数值解法将由溶个一阶方程组成的常微分方程初值问题: 323争匕,…,〃)aX孕=Zl(XJl，乃,)学■=力区乃必,)34x45)aX乃⑷=/匕3)=%.%3)=&(8,28)写成向量形式:rdY_ 324预估-校正公式:J）港（勺，Xe，zDj伊+i]=仔"+乂4）+（勺+卜乂+14+]）1J|zj2[|g（G%zjJ加iZ+i）四阶龙格・库塔公式：及「中二」（％+九"+媚”4+%婷））,4岑J苫（々+〃，M+的,Z*+碣2）例8.7两种果树寄生虫，其数量分别是"=""），u=u⑥，其中一种寄生虫以吃另一种寄生虫为生,两种寄生虫的增长函数如下列常微分方程所示，预测3年后这一对寄生虫的数量。 325—=0.09u(l--)-0.45uudt20duu«—=0.06u(l-—)-0.00山uw(0)=1.6尸(0)=1.2=0.09〃(1-合-0.45uug(ufu)=0.06u(l---0.00山u解：记IIJ在本题中/0〃，")=」(〃，u)，g(E，〃，u)=gQ，u)。用欧拉预估.校正公式(8.31)：.+J[%'■/(%，4+1)'[g(“*+l，Ux+l)i取力=1,计算结果见表8.4。t(年)W)必）11.61.220.287661.3310330.2913081.4777640.000217271.63947表8.4计算结果8.4.2高阶常微分方程数值方法卜.面以三阶常微分方程为例说明高阶常微分方程的数值计算步骤。 326ax<丁(。)=4°)(aix 327定性是把方法用到最为典型的稳定微分方程(8.35)上进行的。血” 328图8.4欧拉方法的绝对稳定区域例8.9对方程(8.35)讨论向后欧拉方法的稳定性。解：用向后欧拉方法计算方程(8.35)的公式：居+「居+兀机+1误差方程：%+。；+1=居+兄私+d)P*4=P*+刃计算相邻两步误差的比值：同”尚=1—1<1当我(㈤,1^1，当绝对稳定区域是左半平面时，则称在这个区域数值方法是上稳定的，或称无条件绝对稳定。 329例8.10对方程(8.35)讨论中心差分方法的稳定性。解：用中心差分方法计算方程（8.35）的公式:%+1=%-1+2物误差方程：0*+1=2-1+2兀力4（8.38）不失一般性，我们考虑由0,Ql对以后心的影响，差分方程（8.36）的特征方程为:记-2助4-1=0它的两个根为：1=助+Jl+（&）2备=助-Jl+（&）2由差分方程理论（略），°”的一般解可由下式表达：Q*=ag+J+（">+以a-J1+«&）2C（&39）其中°e可°。，巧决定。Xh-J1+（双）2=兀〃-1+0G1方）2由于Re（砌<0,故|助7|〉1以a-J1+（曲2）”因力增大而恶性发展，所以方法为不稳定的。可以证明，龙格-库塔方法以及显式或隐式的阿达姆斯方法均是绝对稳定的。 3308.6程序示例程序8.1用改进的欧拉公式来求解常微分方程初值问题:'y'(x)=/(")(a 331scaf("％lf”,&yO)；printf(**nlnputmvalue(divide(%f,%f)：”,a,b|；scanf&m)；if(m<=0){printf("Pleaseinputanumberlargerthen1. 332")；return1；)h=(b-a)/m：xn=a；y=yO；for(i=1；i<=m：i++){xnl=xn+h；ynlb=yn+hf(xn,yn)；*ynl=yn+h/2(f(xn,yn)+f(xnl,ynlb))；printfCx%d=%f,y%d=%f 333”,i,xnl,i,ynl)；xn=xnl；yn=ynl；}return0：}计算实例用改进的欧拉公式，求解常微分方程初值问题的解：'dx(04x40.4)7(0)=1L.程序输入输出Inputthevalueat0.000000:1 334Inputmvalue(divide(0.000000,0.400000)|：4xl=0.100000,yl=l.110500x2=0.200000,y2=1.248276x3=O.3OOOOO,y3=1.424760x4=0.400000,y4=1.658736程序8.2用四阶龙格-库塔方法求解常微分方程初值问题:'y'(x)=/(")(a 335inti；doublea,b,yO；doublexn,yn»ynl；doublekl,k2,k3»k4；doubleh；printf('、nlnputthebeginandendofx：")；scanf("％lf%lf”,&a,&b)；printf("Inputtheyvalueat%f：”,a)；scafC%lfM,&y0)；printf(44nlnputmvalue[divide(%f,%f)：”，a,b]；scanf("％d“，&m)；if(m<=0)(printf("Pleaseinputanumberlargerthen1. 336")；return1；)h=(b-a)/m；xn=a；yn=yO；for(i=1；i<=m；i++)(kl=f(xn,yn)；k2=f((xy+h/2),(yn+hkl/2))； 337k3=f((xn+h/2),(yn+hk2/2)；k4=f((xn+h),(yn+hk3)；ynl=yn+h/6(k1+2k2+2k3+k4)；xn+=h；printf(4 3381.用向前欧拉式解初值问题:'力=+2〈五一'0.14x40.5，取为=0.1双0)=12.用向后欧拉公式解初值问题：dy2.—=x+y— 339“(0)=0.193u(O)=0.083
此文档下载收益归作者所有
