计算方法教案77316

1、第1章误差分析与数值计算3§1.1 引言3§1.2 绝对误差与相对误差、有效数字9§1.3 近似数的简单算术运算12§1.4数值计算中误差分析的一些原则13第2章非线性方程（组）的近似解法15§2.1 引言15§2.2根的隔离16§2.3 对分法17§2.4 迭代法19§2.6 弦截法21§2.6 弦截法22§1.7 用牛顿法解方程组23本章小结25第3章线性方程组的解法26§3.1 引言26§3.2 高斯消去法28§3.3 矩阵的LU分解31§

2、3.4 对称矩阵的LDLT分解32§3.5 线性方程组解的可靠性33§3.6 简单迭代法34本章小结43第4章矩阵特征值与特征向量的计算44§4.1 引言44§4.2幂法和反幂法45§4.3 雅可比方法46§4.4 QR方法*51本章小结52第5章插值与拟合53§5.1 引言53§5.2 插值多项式的存在和唯一性54§5.3 拉格朗日插值多项式55§5.4 均差插值公式57§5.5 差分等距结点插值公式59§5.6 爱尔米特插值公式61§5.7 分段低次插值62&#

3、167;5.8 三次样条函数63§5.9 曲线拟合的最小二乘法67本章小结70第6章数值积分和数值微分71§6.1 引言71§6.2 牛顿一科特斯型积分公式72§6.3 复合求积公式74§6.4 龙贝格求积公式77§6.5 高斯求积公式78§6.6 二重积分的数值积分法80§6.7 数值微分81本章小结83第7章常微分方程的数值解法84§7.1 引言84§7.2 欧拉法和改进的欧拉法85§7.3 龙格-库塔方法86§7.4 线性多步法89§7.5 算法的稳定性与收敛

4、性91§7.6 微分方程组和高阶微分方程解法92本章小结94第1章 误差分析与数值计算§1.1引言1、课程任务和目的：在第七届国际软件工程学术会议上，“计算方法”被列入应用方法学的研究领域，强调了计算方法的研究应用与软件方法学的研究密切结合。这就说明了计算方法与软件之间的联系以及在应用软件研制中的地位与作用，计算方法是研究各种数学问题求解的数值计算方法。在计算机成为数值计算的主要工具的今天，则要求研究适合于计算机使用的数值计算方法。计算方法就是研究用计算机解决数学问题的数值方法及其理论，它的内容包括函数的数值逼近、数值微分与数值积分、非线性方程值解、线性方程组数值解、常微和

5、偏微数值解等，即都是以数学问题为研究对象的。因此，计算方法是数学的一个分支，只是它不象纯数学那样只研究数学本身的理论，是把理论与计算紧密结合，着重研究数学问题的数值方法及其理论，计算方法是计算机应用和软件研制开发的重要组成部分，通过本课程的学习和上机实习，使学生掌握利用计算机进行科学计算的基本理论和基本方法，并且学会将基本理论和基本方法应用于软件开发以及软件研制。2、本课程基本要求（1）掌握方法的基本原理和思想。（2）掌握方法处理的技巧及与计算机的结合。（3）掌握误差分析，收敛性及稳定性的基本理论。（4）学会进行可靠的理论分析，对近似计算要确保精度要求，要进行误差分析。（5）通过例子，学习使用

6、各种计算方法解决实际计算问题。（6）通过上机实践，能编写算法和实现算法。（7）掌握数值计算中一些最基本、最常用的计算方法和算法。3、本课程与各课程的关系：由于本课内容包括了微积分、代数、常微分方程的数值方法，学生必须掌握这几门课的基本内容才能学好这一课程，同时，学习此课程还必须具备计算机系统的初步知识，掌握一门常用的高级语言，如： BASIC、PASCAL、C语言等，并须具备一定的编程能力。4、本课程的特点：（1）面向计算机，要根据计算机特点提供实际可行的有效算法。即算法只能包括加、减、乘、除运算和逻辑运算，是计算机能直接处理的。（2）有可靠的理论分析，能任意逼近并达到精度要求，对近似算法要保

7、证收敛性和数值稳定性，还要对误差进行分析，而且都是建立在相应数学理论基础上的。（3）有好的计算复杂性。时间复杂性好是指节省时间；空间复杂性好是指节省存储量。这也是建立算法时要研究的问题，因为它关系到算法能否在计算机上完成。（4）要有数值实验。即任何一种算法除了从理论上要满足上述三点外，还要通过数值实验证明是行之有效的。计算方法最基本的立足点是容许误差，在误差容许的范围内对某一数学问题进行近似计算，得到能满足要求的近似结果。现实世界中误差是普遍存在的，由于世界上没有绝对精确的量具（绝对精确的量具是没有刻度的），因此人类通过量具采集的数据都是近似值，另一方面，我们的生产、实验工具都不是绝对精确的，

8、这就使得人类在生产和科学实验中必需容许误差。计算机的应用可以分为二个方面，即数值计算和非数值计算。利用计算机进行数值计算的过程如下图所示：在上图中，计算方法的任务是：由建立的数学模型给出可编程并由计算机能完成的计算方法，然后编程和上机求解。由于计算方法是编程后可由计算机求解的近似计算方法，如何确保近似解的精度显得尤为重要，必须深入讨论有关误差的基本概念和基本理论，为近似计算的精度分析打下基础。1、误差的来源（种类）误差的来源主要有以下四种（1）模型误差：建立数学模型时的误差。例如：在求重量的数学模型G=m*g 中，重量G不是仅与质量和重力加速度有关，它还与温度、测量地点的海拔、地层结构等众多因

9、素有关，为了使模型较为简单和实用，采用抓住主要矛盾的方法，去掉了大量对重量影响不大的次要因素，建立了上述重量的近似模型，由此产生了模型误差。（2）观测误差：采集数据时的误差。采集数据时，通常是依靠仪器和量具，由于没有绝对精确的仪器和量具，因此采集的数据有误差，此误差称为观测误差。（3）舍入误差：由于计算机字长有限而产生的误差。硬件再发展，计算机的字长总是有限的，在计算过程中，当数据的长度超过了计算机的字长时，计算机就会进行四舍五入，由此产生的误差称为舍入误差。（4）截断误差：无限形式的有限化而产生的误差。在计算中有时会运用无限形式的计算公式，例如台劳公式：显然此公式无法进行计算，因此必需根据实

10、际需要，从某一项起将后面的各项截断，即由此产生的误差称为截断误差。§1.2绝对误差与相对误差、有效数字为描述方便，首先约定x*是精确值x的近似值。引入误差的概念，其目的是为了衡量近似值x*的好坏。（1）绝对误差： x*-x 由于精确值x通常无法确定，因此绝对误差无法计算，由此引入绝对误差限的概念。绝对误差限：绝对误差的一个上界。即：若 | x*-x | £ e，则称e为x*的绝对误差限。绝对误差限的性质是：A.不唯一这是因为| x*-x | 的上界是不唯一的。B.可确定只要我们对x*的实际背景有一定的了解，就不难确定| x*-x |的上界。例如，x*表示身高，则| x*-x

11、 |的上界可为3米。当x*是你求出的，那么为了说明你的工作认真，你一定会将| x*-x | 的上界估计得尽量小，因此在这种意义上绝对误差限可用来衡量x*的好坏。由于绝对误差限没有考虑问题的规模，因此有时它也不能衡量x*的好坏。例如：x是地球与太阳的距离，y是分子中二个原子间的距离，若| x*-x | £1公里，| y*-y | £1厘米，则并不能说y*比x*精确。由此引入相对误差和相对误差限的概念。（2）相对误差：(x*-x ) / x* 相对误差限：相对误差绝对值的一个上界。3、有效数字这里我们必须搞清楚什么是有效数字以及如何确定x*有几位有效数字。（1）有效数字的定义若

12、|x*-x|<x*的某一位的半个单位，则称x* 精确到这一位，并从这一位开始，一直到前面第一个不为零的数都是x*的有效数字。此定义实际上定义了什么叫精确到某一位和什么叫有效数字。例如：若x*精确到小数点后第3位，即指| x*-x | £ 0.5 ´10-3。（2）有效数字的判定方法方法一：四舍五入此方法首先确定x*是由x的哪一位四舍五入产生的，然后从这一位的前一位开始一直到前面第一个不为零的数都是x*的有效数字。例1 若x=0.872596, x*=0.87，求x*的有效位数。解： x*是由x的小数点后第三位四舍五入产生的，所以x*有二位有效数字。注意，方法一判定有效

13、数字很简单，但有时会失效。例如，若x=0.272987 x*=0.273102，此时无法用方法一确定x*的有效位数，原因是x*不是由x四舍五入产生的，在这种情况下，必须用有效数字的定义来确定x*的有效位数。即方法二：用定义此方法首先计算| x*-x |，再判断它小于等于x*的哪一位的半个单位，然后从近一位开始，一直到第一个不为零的数都是有效数字。例2 若x=0.62073，x*=0.6207，确定x*的有效位数。解：因为| x*-x | £ 0.0003 £ 0.5´10- 4，x*精确到小数点后第4位，所以x*有四位有效数字。例3 若x=0.080199，x*=

14、0.802，确定x*的有效位数。解：因为| x*-x |=0.00001£0.5´10- 5，所以0.5´10- 3，推出x*有三位有效数字。例4 若x=6.28936，x*=7.3132，确定x*的有效位数。解：| x*-x |=0.02357£0.5´10- 1，所以x*有二位有效数字。§1.3近似数的简单算术运算近似数的加法设有n的近似数xk*>0(k=1,2,n),其准确值为xk>0,(k=1,2,n)的绝对误差限E(x*)(1)和的绝对误差等于各项绝对误差之和。(2) 和的绝对误差限不超过各项绝对误差限之和。类似

15、的可以得到：和的相对误差限介于各项相对误差限的最小者与最大者之间。近似数的乘法结论：乘积的相对误差限不超过各项相对误差限之和。近似数的除法结论：商的相对误差限不超过被除数与除数相对误差限之和。近似数的幂和根（见教材p9）近似数的对数（见教材p9）近似数的减法结论：差的绝对误差等于各项绝对误差之和。注意两个几乎相等的近似数相减会使结果的有效数字损失，影响整个计算工作的准确性，应尽量避免。（当|x|很小时，x1与x2很接近时）§1.4数值计算中误差分析的一些原则为保证计算结果的高精度，在进行数值计算时应遵循下述几个原则。（1）在进行除法时，要避免除数的绝对值<<被除数的绝对值

16、。为什么要“避免”？若不“避免”，则除出的结果很大，由于计算机字长有限，它装不下，因此会进行四舍五入，一个很大的数进行四舍五入时舍去的部分也会很大，这会使舍入误差变大。怎样“避免”？因为用户只关心最后的计算结果，当中间计算过程中出现了除数的绝对值<<被除数的绝对值时，就应该换一种计算方法，以避免这种情况的发生，以后我们将会针对具体的计算问题来讨论“避免”的方法。（2）在进行减法时，要避免二个相近的数相减。为什么要“避免”？若不“避免”，就可能失去大量的有效数字，例如：若a=30001和b=30000都有五位有效数字，因为a-b=1，所以结果至多有1位有效数字。怎么“避免”？“避免”

17、的思路与第1个原则中“避免”的思路相同，须针对具体计算问题来讨论。（3）要防止“大数吃小数”什么是“大数吃小数”？我们用一个例子为说明。计算8756294874，其中n=10 20，0< ai<10-6。此题是一个很大的数与很多很小的数相加，若采用将大数依次与a1,a2,¼,an相加，由于计算机字长有限，因此在与ai相加时会进行四舍五入将ai舍去，这样，最后的结果仍是大数，这就是大数将a1,a2,¼,an吃掉了。为什么要“避免”？尽管每个小数都很小，但它们很多，可能它们的和比大数还大，而最后计算工结果为大数，显然误差可能很大。怎样“避免”？有的同学提出先将小数相

18、加，然后再与大数相加，这个思路是对的，但有一个漏洞，因为小数相加到一定程度也会变成大数，它也开始吃小数了。可以采取分部相加的方法解决。第2章 非线性方程（组）的近似解法§2.1 引言方程f(x)=0的解称为方程的根。也叫做函数f(x)的零点。方程求根大致包括三个问题（1）方程有没有根？如果有根，有几个根？（2）哪里有根？求有根的区间，区间内的任意一点作为根的近似值。（3）根的精确化，已知一个根的近似值后设法逐步把根精确化，直到足够精确为止。本课程主要研究问题（2）和（3）。§2.2 根的隔离求方程f(x)=0的解的近似值时，首先要确定若干个区间，使每个区间内只有的一个根，这

19、个步骤称为根的隔离。对一般的方程，根的隔离有两种方法（1）试值法。求出f(x)在若干点上的函数值，观察函数值符号变化的情况，从而确定隔根区间。（2）作图法。画出y=f(x)的草图，观察曲线y=f(x)与x轴交点的大致位置，从而确定隔根区间。例1.2.1讨论方程f(x)=2x3-4x2+4x+2= 0 的根的位置。f1=inline('2*x3-4*x2+4*x+2'),fplot(f1,-1,1), f2=inline('log(x)-1/x'), fplot(f2,1,2)-1012-10-50510例1.2.2将方程xlog(x)= 1 的根进行隔离。11.

20、52-1-0.500.5§2.3 对分法设有方程f(x)=0在（a b）内有且仅有一个根，这时有f(a) f(b)<0可用对分法求的近似值，方法如下（1）准备：计算区间（a b）两个端点的函数值f(a)， f(b)（2）对分：取c=(a+b)/2为（a b）的中点，计算f(c)（3）判断：如果f(c)=0，则 c为f(x)=0的根，否则检验：若f(c)f(a)<0，则方程的根位于a c内，用c代替b，若f(c)f(b)<0，则方程的根位于c b内，用c代替a。（4）检验：若|b-a|<e (e为精度要求)此时计算结束，取=c，否则转（2）。例2.3.1用对分法

21、求方程f(x)=x3+2x-5= 0 在 1 2 内的根，e=10-3。f=inline('x3+2*x-5')ep=1e-3;a=1;b=2;while abs(b-a)>epc=(a+b)/2;if f(c)*f(a)>0a=c;elseb=c;enda,b,abs(b-a),pauseendfplot(f,1 2)有根区间 1.0000 2.0000 1.0000 1.5000 1.2500 1.5000 1.2500 1.3750 1.3125 1.3750 1.3125 1.3438 1.3281 1.3438 1.3281 1.3359 1.3281 1

22、.3320 1.3281 1.3301 1.3281 1.3291方程的解 x= 1.3286§2.4 迭代法设有方程f(x)=0在a b上有且仅有一个根，可用迭代法求的近似值，方法如下（1）将方程f(x)=0写成迭代形式x=j(x)（2）在a b上任取一个初始值x0。（3）计算x1=j(x0)（4）若| x1-x0|<e (e为精度要求)，此时计算结束=x1，否则令x0=x1转（3）。定理(收敛定理) 设函数j(x)在a b上连续，在(a b)内可导，且满足(1)当xa b时,j(x)a b。(2)当xa b时，|j¢(x)|m<1,m为一个常数。则以下结论成

23、立：(1) 在a b上j(x)有且仅有一个根。(2) 对任意x0a b,由迭代公式xn+1=j(xn)，n=0,1,2,产生的数列|xn|a b,且有xn®(n®)。(3) 成立误差估计式讨论：(1) 此定理说明只要|xn -xn-1|足够小，就可以保证xn充分接近方程的根。所以在实际计算时，用条件|xn -xn-1|e控制迭代过程的结束。(2)定理指当对任意x0a b作为初始条件，迭代过程都收敛，这种形式的收敛定理称为大范围收敛定理，是一个充分条件。当此条件不成立时，往往可以取初始值x0接近于方程的根，使迭代过程收敛。(见P21定理2.4.2)。(3) 从定理2.4.1中

24、可以看出，收敛速度与m值有很大关系，m越小，收敛速度越快。当m接近于1时，收敛速度很慢。根据|j¢(x)|m<1，所以采用迭代法求解方程的根时，应适当的构造j(x)使|j¢(x)|在区间a b内尽量的小。(见下面的例子)例2.4.1用迭代法解方程x= 10x-2 ，x0=1分别采用迭代格式 x= 10x-2 和x=log(x+2)，观察两个计算过程的区别。e=1e-3迭代过程: 1.0000 0.4771 0.3939 0.3791 0.3764 0.3759 迭代6次x= 0.3759f=inline('log10(x+2)')x=1x=f(x)例2

25、.4.2用迭代法求方程f(x)=x3+2x-5= 0的根，x0=1 。迭代过程:1.0000 1.4422 1.28371.3449 1.3220 1.3306 1.3274 1.3286 1.3281f=inline('(5-2*x)(1/3)')x=1x=f(x)迭代9次x= 1.3281§2.5 牛顿迭代法牛顿法是解方程f(x)=0的重要方法，它也是一种迭代法。设有方程f(x)=0在a b上有且仅有一个根，可用牛顿法xn+1=xn- f(xn)/ f¢(xn)求的近似值，具体计算步骤如下（1）求出f(x)的导数f¢(x)，在a b上任取一个初

26、始值x0。（2）计算x1= x0- f(x0)/ f¢( x0)（3）若| x1-x0|<e (e为精度要求)，此时计算结束=x1，否则令x0=x1转（2）。例2.5.1用牛顿法解方程f(x)=x3-2x2-4x-7= 0在 34 内的根 x0=4。迭代过程: 4.0000 3.6786 3.6329 3.6320 迭代4次x= 3.6320f=inline('x3-2*x2-4*x-7')fd= inline('3*x2-4*x-4')x=4x=x-f(x)/fd(x)§2.6 弦截法弦截法也是一种是解方程f(x)=0的迭代法，它的特

27、点是不需要计算f(x)的导数f¢(x)，且收敛速度也相当快，是工程计算中常用的算法之一。设有方程f(x)=0在a b上有且仅有一个根，可用弦截法求的近似值，方法如下（1）求函数f(x) 在区间a b的两个端点的函数值 f(x0) ,f(x1) ，其中a= x0 ,b= x1（2）计算x2= x1- f(x1) x1- x0/ f(x1)- f(x0)（3）若| x2-x1|<e (e为精度要求)，此时计算结束=x2，否则令x0=x1，x1=x2转（2）。f=inline('x3-2*x2-4*x-7')x0=3,x1=4x2=x1-f(x1)*(x1-x0)/(

28、f(x1)-f(x0)x0=x1,x1=x2例2.6.1用牛顿法解方程f(x)=x3-2x2-4x-7 = 0在 34 内的根。§1.7 用牛顿法解方程组设有非线性方程组u (x,y)=0，v (x,y)=0，在 (xn, yn)按台劳级数展开，取展开式的前两项得到其中(xn, yn)是根的第n次近似值，如果Jn¹0方程组的第n+1次近似值(xn+1, yn+1)可用以下公式计算例1.7.1用牛顿法解方程组u (x,y)=x3+y3-4= 0 v (x,y)=x4+y2-3= 0 (¶u/¶x= 3x2,¶u/¶y= 3y2,

29、2;v/¶x=4x3,¶v/¶y=2y)x0=1;y0=1.4;J0=det(3*x02,3*y02;4*x03,2*y0);for k=1:10 u0=x03+y03-4; v0=x04+y02-3; Jx=det(3*y02,u0;2*y0,v0); Jy=det(u0,3*x02;v0,4*x03); x1=x0+Jx/J0; y1=y0+Jy/J0; x0=x1;y0=y1; x0 y0,pauseend迭代初值x0=1, y0=1.4例用牛顿法解方程组u (x,y)=2x3-y2-1= 0v (x,y)=xy3-y-4= 0 (¶u/¶

30、;x= 6x2,¶u/¶y= -2x,¶v/¶x=6y3,¶v/¶y=2xy2-1)迭代初值x0=1.2, y0=1.7。例用牛顿法解方程组u (x,y)=x-cos(y)= 0，v (x,y)=y-sin(x)= 0 (¶u/¶x= 1,¶u/¶y= sin(y),¶v/¶x= -cos(x),¶v/¶y=1)迭代初值x0=0, y0=0。本章小结为了比较各种迭代方法的收敛速度，我们引入收敛阶的概念。设迭代过程xn+1=j(xn) 收敛于方程x=j(x)

31、的根，令e n= xn-，e n称为迭代误差，如果存在实数P³1和非零常数K，使得，则称该迭代过程为P阶收敛的。P=1称为线性收敛，P>1称为超线性收敛，P=2称为平方收敛，显然P越大，迭代过程收敛的越快。可以证明当是方程f(x)=0的单根时，牛顿法是平方收敛的。当是方程f(x)=0的重根时，牛顿法仅为线性收敛。对分法：收敛速度与公比为1/2的等比级数相同。收敛速度最慢，但简单有效，不存在发散问题。弦截法：弦截法的收敛阶P=1.618。收敛速度次之，不需要计算f(x)的导函数，有发散问题。牛顿法：牛顿法是平方收敛，收敛速度最快，但要计算f(x)的导函数，有发散问题。第3章 线性

32、方程组的解法§3.1 引言在科学实验和工程设计中，经常用到解线性方程组的问题。本章讨论用计算机求解线性方程组的两类主要方法：直接法和迭代法。解线性方程组的一般表达式 根据矩阵的运算的性质有简记为Ax=b 其中 方程组Ax=b有唯一解的充分必要条件是det(A)¹0。理论上，可以用克莱姆法则xk=k/,（k=1,2.n）其中=det(A)，是A中第k列换成向量所得到的矩阵的行列式。计算需要的乘法次数为(n-1)(n+1)!,当n较大时，计算量很大。这种方法效率低，在实际应用中很少使用。解线性方程组的方法可以分为两类：一类是直接法，其特点是只包含有限次的四则运算，在每次运算都无

33、舍入误差的情况下，所得到的是方程组的准确解。由于实际计算中总是有舍入误差，所以实际得到的也是近似解。令一类是迭代法，它首先选择一组初始值，再运用同样的计算步骤，重复计算，得到近似解。由于这类方法中出现了极限过程，必须研究迭代过程的收敛性。本章主要介绍：直接法中的高斯消去法和主元高斯消去法。迭代法中的简单迭代法和塞德尔单迭代法。§3.2 高斯消去法顺序高斯消去法以n=4为例说明高斯消去法的计算过程，设有线性方程组ÞÞÞ经过3次消元步骤，得到以上形式。从最后一个方程中解出x4,依次回代得到方程组的全部解。计算过程见教材(P37)算法3.2,1主元消去法例 解

34、线性方程组 解：按4位小数计算，得到x1=0.6667,x2=0.0000, 准确解(x1=1/3, x2=2/3)，误差很大。如果将方程组的顺序对换，得到按4位小数计算，得到x1=0.3334, x2=0.6667, 准确解(x1=1/3, x2=2/3)，误差很小。在顺序主元消去法中，如果遇到方程组Ax=b中A的主对角线元素为零时，计算过程不能进行。如果A的主对角线元素的绝对值很小，也会产生较大的舍人误差，使得最后得到的解与准确解有较大的误差。工程中使用较多的是列主元高斯消去法，计算过程见教材(P42)算法3.4ì 6x1+3x2+2x3=6例用顺序消去法解方程组í10

35、x1+5x2+6x3=0î 8x1+5x2+3x3=0方程组的增广矩阵A|b 6 3 2 6 10 5 6 0 8 5 3 0消元 6.0000 3.0000 2.0000 6.0000 0 0 2.6667 -10.0000 0 1.0000 0.3333 -8.0000方程组系数矩阵主对角线元素为零，消元过程无法进行!例用主元高斯消去法解例中的方程组。方程组的增广矩阵A|b 6 3 2 6 10 5 6 0 8 5 3 0选主元 10 5 6 0 6 3 2 6 8 5 3 0消元 10.0000 5.0000 6.0000 0 0 0 -1.6000 6.0000 0 1.00

36、00 -1.8000 0选主元 10.0000 5.0000 6.0000 0 0 1.0000 -1.8000 0 0 0 -1.6000 6.0000消元 10.0000 5.0000 6.0000 0 0 1.0000 -1.8000 0 0 0 -1.6000 6.0000回代得到方程组的解 5.6250 -6.7500 -3.7500§3.3 矩阵的LU分解定理设矩阵的各阶顺序主子式不等于Dk0 (k=1,2,n)，则A有唯一的LU分解，A=LU,其中L为单位下三角矩阵，U为上三角矩阵。对于线性方程组Ax=b如果已经得到了的LU分解，方程组的解可以通过以下方法得到LUx=b

37、,令Ux=y,则有Ly=b,所以，求解Ax=b如的问题等价于求解两的系数矩阵为三角方阵的方程组。求解方程组Ax=b的计算变得非常简单。例 解线性方程组Ax=b，PA=LU, LUx=Pb,Ly=Pb,Ux=yA=magic(3), b=1;2;3, L,U,P=lu(A), y=L(P*b)，x=Uy§3.4 对称矩阵的LDLT分解设矩阵的各阶顺序主子式不等于Dk0 (k=1,2,n)，则A有唯一的LU分解，A=LU,其中L为单位下三角矩阵，U为主对角线元素不为0的上三角矩阵。从而A=LDV,若为对称矩阵A=AT，则有A=AT=(LDV)T=VT(DLT)，其中为VT为单位下三角矩阵

38、，(DLT)为主对角元素不为0的上三角矩阵，由LU分解的唯一性得到，L=VT,于是A=LDLT。定理设对称矩阵的各阶顺序主子式不等于Dk0 (k=1,2,n)，则唯一确定一个单位下三角矩阵L和主对角元素不为0的对角矩阵D，使A=LDLT。例3.4.1 A=5 10 5 -5;10 24 -2 -6;5 -2 44 -11;-5 -6 -11 23, L,P=chol(A), L'*L如果线性方程组Ax=b系数矩阵是对称正定的，A=LTL, LTL x=b,令L x=y, LTy=b例b=5 10 -1 -19',y=L'b,x=Ly§3.5线性方程组解的可靠性

39、 向量范数和误差向量定义设| .|是定义在Rn上的实值函数，如果对于Rn中的任意向量x,y及任意非负实数k，满足(1)（非负性）对任意向量x都有|x|>0且|x|=0的充分必要条件是x=0。(2)（正齐性）对任意实数k和任意向量x都有|kx|=k|x|。(3)（三角不等式）对任意向量|x+y|x|+|y|。称| .|为向量范数。定理对Rn上的任意向量x=(x1,x2xn)T,记| x |1，| x |2，|x |都是向量范数。定义3.5.2设| .|是定义在Rnxn上的实值函数，如果对于Rnxn中的任意矩阵A,B及任意非负实数k，满足(1)（非负性）对任意矩阵A都有|A|>0且|A

40、|=0的充分必要条件是A=0。(2)（正齐性）对任意实数k和任意矩阵A都有|kA|=k|A|。(3)（三角不等式）对任意矩阵A,B都有 |A+B|A|+|B| , |AB|A| |B|称| .|为矩阵范数。定理3.5.2对Rnxn上的任意矩阵A记| A |1，| A |2，|A |都是矩阵范数。A=-5 -3 1;4 0 1;4 1 -2, x=2 -5 1'norm(x,1), norm(x,2), norm(x,inf)norm(A,1), norm(A,2), norm(A,inf)norm(A*x,1), norm(A*x,2), norm(A*x,inf)定理3.5.3对Rn

41、上给定一种向量范数，对于Rnxn中的任意矩阵A，则是一种矩阵范数，并且它与给定的向量范数相容。例计算各种范数。定义设| .|是定义在Rnxn上的范数，对于Rnxn中的任意矩阵A称cond(A)=| A | | A-1 | 为矩阵A的条件数。例3.5.3A=1 1.0001;1 1 ,b=2;2 , b1=2.0001;2 ,Ab, Ab1当一个线性方程组Ax=b的系数矩阵A的条件数很大时，称方程组是病态的，判断线性方程组Ax=b是否为病态的方法是：1、 当|det(A)|很小或A的一些行或列近似相关时，方程组可能是病态的；2、 使用列主元消去法，出现主元的绝对值很小时，方程组可能是病态的；3、

42、 当A的元素在数量级上有很大差别，且无一定规律时，方程组可能是病态的。对于病态的方程组可以采取下面的方法：1、 采用高精度的算术运算，例如采用双精度数运算。2、 当A的元素在数量级上差别很大时，采用行平衡或列平衡的方法可以降低A的条件数。A=1 1; 1 1e4 , b=2;1e4 , cond(A,inf), s=max(A) , d=diag(1./s) ,cond(d*A,inf)clear,n=4; %取n=6,8,10 等进行比较b=ones(n,1)x=hilb(n)bx0=invhilb(n)*bmax(abs(x-x0)cond(hilb(n)rats(hilb(4) %字符串

43、rats(invhilb(4) %字符串§3.6简单迭代法设有方程组Ax=b，变为迭代形式，x=Mx+f，或x(k+1)=Mx(k)+f，任取初始值x(0) 程迭代得到x(0)，x(1)，x(2) ，¼，x(k) ，¼若极限 存在，则x*就是原方程组的解。以n=4为例­­­­x(k+1)M x(k)f写成分量形式定理若存在某种范数|. |，使|M|<1,则简单迭代法收敛。定理 迭代公式x(k+1)=Mx(k)+f，对任意初始值x(0)和f都收敛的充分必要条件是矩阵M的各个特征值的模都小于1。在简单迭代法的基础上作改进x

44、(k+1)=M1x(k+1)+ M2x(k)+f，以n=4为例­­­­­­x(k+1) M1x(k+1)M2x(k)f§3.7雅可比迭代法与高斯-塞德尔迭代法雅可比迭代法设方程组Ax=b，的系数矩阵主对角线元素不为零，A可以分解为A=D+L+U,其中D，L，U分别取A的主对角线，下三角和上三角部分的元素。由(D+L+U)x=b，Dx=-(L+U)x+b, x=-D-1 (L+U)x+D-1b,由此可以得到迭代形式x(k+1)=-D-1(L+U)x(k)+D-1b(k=0,1,2) (3-67)雅可比迭代法的迭代矩阵 MJ=-

45、D-1(L+U)。写成分量形式(3-68)高斯-塞德尔迭代法在计算(3-68)式时，用已经计算出的xj进行后续的计算，可以得到(3-69)称式(3-69)为高斯-塞德尔迭代法。式(3-69)的矩阵表达形式为x(k+1)=-D-1Lx(k+1)-D-1U)x(k)+D-1b(k=0,1,2)x(k+1)=-(D+L)-1Ux(k)+(D+L)-1b(k=0,1,2) (3-70)高斯-塞德尔迭代法的迭代矩阵 MGS=-(D+L)-1U。雅可比迭代法和高斯-塞德尔迭代法的收敛性定理3.7.1雅可比迭代法和高斯-塞德尔迭代法收敛的充分必要条件是它们的迭代矩阵的各个特征值的模都小于1。定理若存在某种范

46、数|. |，使|MJ|<1,则雅可比迭代法收敛。若存在某种范数|. |，使|MGS|<1,则高斯-塞德尔迭代法收敛。例3.7.1分别用雅可比迭代法和塞德尔迭代法解方程组误差e<10-6雅可比迭代法迭代公式x(k+1)=-D-1(L+U)x(k)+D-1b高斯-塞德尔迭代法迭代公式 x(k+1)=-(D+L)-1Ux(k)+(D+L)-1bA= 10 -1 -2;-1 10 -2;-1 -1 5b=72 83 42'D=diag(diag(A);U=triu(A,1);L=tril(A,-1);M=-inv(D)*(L+U);b=inv(D)*b;x=b;for k=1

47、:10,x=M*x+b,endA= 10 -1 -2;-1 10 -2;-1 -1 5b=72 83 42'D=diag(diag(A);U=triu(A,1);L=tril(A,-1);M=-inv(D+L)*U;b=inv(D+L)*b;x=b;for k=1:5,x=M*x+b,end A= 2 -1 1;1 1 1;1 1 -2D=diag(diag(A);U=triu(A,1);L=tril(A,-1);MJ=-inv(D)*(L+U);MS=-inv(D+L)*U;abs(eig(MJ)abs(eig(MS) 例3.7.2考察用雅可比迭代法和塞德尔迭代法解方程组的收敛性。雅

48、可比迭代法迭代公式x(k+1)=-D-1(L+U)x(k)+D-1b高斯-塞德尔迭代法迭代公式 x(k+1)=-(D+L)-1Ux(k)+(D+L)-1b例3.7.3用雅可比迭代法解方程组A= 1 -2 2;-1 1 -1;-2 -2 1b=1 0 -2'D=diag(diag(A);U=triu(A,1);L=tril(A,-1);M=-inv(D)*(L+U);b=inv(D)*b;x=b;x=M*x+b 误差e<10-3雅可比迭代法迭代公式x(k+1)=-D-1(L+U)x(k)+D-1b本章小结本章讨论了解线性方程组的直接解法和迭代解法。直接解法比较适用与系数矩阵稠密（既

49、零元素较少）的中、小型线性方程组，但对系数矩阵是带状或近似带状的大型线性方程组也适用。直接解法中的列主元高斯消去法具有精度较高和省时的优点，是计算机中常用的算法。迭代解法中主要介绍了雅可比迭代法、高斯-塞德尔迭代法和松弛法。迭代法具有计算公式简单、程序设计容易、占用计算机内存较少的优点。适用于解大型稀疏矩阵（既零元素较多）线性方程组。高斯-塞德尔迭代法是在雅可比迭代法的基础上改进得到，在很多情况下可以加快收敛速度，但它的收敛域与雅可比迭代法不同，因此不能互相取代。松弛法可以加速迭代过程的收敛速度，但要适当选择松弛因子（0<w<2）。在选择迭代法时，要特别注意检验方法的收敛性问题。第

50、4章 矩阵特征值与特征向量的计算§4.1 引言求矩阵的特征值和特征向量，是代数计算中的重要问题。在自然科学和工程中的许多问题，例如电磁振荡、桥梁的振动，机械振动等都可以归结为求矩阵的特征值和特征向量问题。矩阵A的特征值和特征向量是指，如果数l和非零列向量x满足关系式Ax=lx，则数l称为A的特征值非零列向量x称为A的与特征值 l 对应的特征向量。计算n阶矩阵A的特征值，就是求特征方程|A-lI|=0的根li (i=1,2,¼,n)。齐次线性方程组 (A-liI)x=0的非零解xi，是li对应的特征向量。本章讨论一些在计算机上计算矩阵的特征值和特征向量的较为稳定的数值算法。&

51、#167;4.2幂法和反幂法1、幂法：计算n阶矩阵A的模最大的特征值（主特征值）及对应的特征向量。任取n维列向量x( 0)，用迭代公式 x( k+1)=A x(k) 计算得到x(0)，x(1)，x(2)，¼设x(0) =a1v1+ a2v2+¼ +anvn ，因为Avj=livj所以x(1) =Ax(0) =a1l1 v1+ a2l2 v2+¼ +anl n vnx(2) =Ax(1) =a1l12 v1+ a2l22v2+¼ +anl n2 vn一般地有x(k+1) =Ax(k) =a1l1 k v1+ a2l2 k v2+¼ +anl n

52、k vn=l1 ka1 v1+ a2(l2/l1) k v2+¼ +an(ln/l1)k vn 当k充分大时x(k+1)»a1l1 k+1 v1»l1x(k)向量x(k+1)与 x(k)向近似地只差一个倍数，这个倍数就是模最大的特征值l1。1、反幂法：Ax=lx，A-1Ax= A-1lx，A-1x=l-1x，即A的特征值的倒数l-1是A的逆矩阵A-1的特征值。用幂法求A-1的模最大的特征值，它的倒数就是A的模最小的特征值。§4.3 雅可比方法对于2阶方阵 令 对于n阶方阵（以n=3为例） 令 作变换矩阵 则有 令 作变换矩阵 则有 令 作变换矩阵 则有

53、一般地说，令，可以将A中的元素的aij和aji变为0。在实际计算中采用以下公式例4.3.1 用雅可比求对称矩阵的特征值和特征向量。（见教材P86）§4.4 QR方法*QR方法是求一般矩阵A的全部特征值和特征向量的一种迭代方法。其基本思路是利用矩阵A的QR分解，其中，R是一个上三角矩阵，Q是正交矩阵。通过迭代格式将A化为相似的上三角矩阵（或分块上三角矩阵），从而求出A的全部特征值和特征向量。例4.4.1用QR方法求矩阵 的特征值。特征向量a= 2 1 3;1 2 3;0 1 2 ,q,r=qr(a) ,a=r*q-0.8165 -0.6215 -0.6760-0.4082 -0.621

54、5 -0.67600.4082 0.4770 -0.2935特征值 本章小结本章介绍了求矩阵的特征值和对应的特征向量的几种方法。幂法可以求出矩阵的主特征值和对应的特征向量，优点是算法简单，但当 |l1/l2| »1时，收敛速度很慢。反幂法可以求出矩阵的模最小的特征值和对应的特征向量。雅可比方法是利用一系列正交相似变换（即平面旋转变换）把实对称矩阵A化为对角阵（近似），从而求出实对称矩阵全部特征值。QR方法是用镜向反射阵将矩阵A作QR分解，是一种求矩阵的全部特征值的有效方法。第5章 插值与拟合§5.1 引言已知表格函数 y=f (x)xix0x1x2¼xn-1xnf

55、(xi)y0y1y2¼yn-1yn构造一个公式p (x) 近似地表示f (x)，解决这个问题的方法有两类：一类是插值法，另一类是拟合法，又称为逼近法。已知函数 y=f (x) 在互异点x0, x1, x2, ¼, xn-1, xn上的函数值y0, y1, y2, ¼, yn-1, yn，构造一个函数p(x)使得p(xi)= yi这样的问题称为插值问题。y=f (x) 称为被插值函数，x0xn 称为插值区间，p (x) 称为插值函数，x0, x1, x2, ¼, xn-1, xn称为插值点，在插值区间内部用p (x) 代替f (x)称为内插，在插值区间外部

56、用p (x) 代替f (x)称为外推，R(x)=f(x)-p(x) 称为插值函数p(x)的误差。§5.2 插值多项式的存在和唯一性定理：给出n+1个插值点及函数值xix0x1x2¼xn-1xnf(xi)y0y1y2¼yn-1yn求一个n次多项式pn(x)=a0+a1x+a2x2+¼+ a nx n(x0,y0)(x1,y1)(x2,y2)(x3,y3)(x4,y4)满足插值条件pn (xi)= yi(i=0,1,2, ¼,n) 的n次插值多项式pn (x)是唯一的。§5.3 拉格朗日插值多项式1、给出2个插值点 (x0 , y0)，(x1 , y1)可以得到一次多项式 2、给出3个插值点 (x0 , y0)，(x1 , y1) ，(x2 , y2)可以得到二次多项式不难验证p2(x) 满足插值条件p2 (x0)= y0p2 (x1)= y1 p2 (x2)= y23、给出n+1个插值点 (x0 , y0)，(x1 , y1),¼,(xn ,