数值计算方法第三章讲义课件

1、第3章 线性方程组的数值解法三、 矩阵的三角分解及其应用四、 线性代数方程组的性态与误差分析五、 迭代法一、 线性方程组二、 高斯(Gauss)消元法六、 共轭梯度法求解上一页 下一页 返回 一、线性方程组 实际问题中的线性方程组分类：按系数矩阵中零元素的个数：稠密线性方程组稀疏线性方程组按未知量的个数：高阶线性方程组低阶线性方程组(如1000)(80%)按系数矩阵的形状对称正定方程组三角形方程组三对角占优方程组 实际中，存在大量的解线性方程组的问题。很多数值方法到最后也会涉及到线性方程组的求解问题：如样条插值的M和m关系式，曲线拟合的法方程，方程组的Newton迭代等问题。解线性方程组的两类

2、数值方法：直接法: 经过有限次四则运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法。（一般有限步内得不到精确解）特点：准确，可靠，理论上得到的解是精确的。适合于中小规模的方程组。特点：速度快，但有误差。适合于大规模的方程组。Cramer法则是一种直接法，但无法实用 直接法概述直接法是将原方程组化为一个或若干个三角形方程组的方法对于线性方程组根据Cramer(克莱姆)法则,若则都是三角形方程组上述方法称为直接三角形分解法不论是Gauss消元法还是直接三角形分解法,最终都归结为解三角形方程组一、三角形线性代数方程组的直接法若记下三

3、角形线性方程组上三角形线性方程组 Gauss消去法二、直接法 其解为其解为:按回代过程求解对方程组，作如下的变换，解不变：交换两个方程的次序一个方程的两边同时乘以一个非0的数一个方程的两边同时乘以一个非0数，加到另一个方程二、Gauss消元法因此，对增广矩阵(A，b )作如下的变换，解不变：交换矩阵的两行某一行乘以一个非0的数某一行乘以一个非0数，加到另一行消元法就是对增广矩阵作上述行的变换，变为三角形线性方程组，而后求根。用矩阵表示：=首先将A化为上三角阵，再回代求解。例3.1.1解化为同解方程组化为上三角方程组按回代过程求解上述求解三元线性代数方程组的方法即为Gauss消元法。(1) 定义

4、行乘数相当于第i个方程-第一个方程行乘数新的第i方程同解！第一个方程不动！ 且(1) 定义行乘数相当于第i个方程-第二个方程行乘数新的第i方程同解！第一、二个方程不动！ 系数矩阵与常数项： Gauss消元法的运算量乘法次数：除法次数：全部回代过程需作乘除法的总次数为于是Gauss消元法的乘除法运算总的次数为数级类似地，可得Gauss消元法的加减法运算总的次数为由于计算机完成一次乘除运算所耗时间要远远多于完成一次加减运算的时间，而且按统计规律，在一个算法中，加减运算和乘除运算次数大体相当，故在衡量一个算法的运算量时只需统计乘除的运算次数。Gauss消元法乘除法约为2700次而如果用Cramer法

5、则的乘除法运算次数约为用行列式定义回代求解 从Gauss消元法的过程自然想到：如果在消元过程中把对角线上方未知量的系数也化为零，使方程组化成每个方程只含一个未知量的方程组，则无须回代就可得到解。这种消元法称为GaussJordan消元法。但经计算，它需要的乘除法次数为 显然，其运算量超过Gauss消元法。所以用GaussJordan消元法求解线性代数方程组是不可取的。不过，用它来求逆矩阵却有方便之处。三、选主元Gauss消元法例3.1.3用Gauss消元法解线性方程组(用3位十进制浮点数计算）解本方程组的精度较高的解为用Gauss消元法求解(用3位十进制浮点数计算)回代后得到主元9999与精确

6、解相比,该结果相当糟糕！究其原因,在求行乘数时用了很小的数0.0001作除数! 上述例题表明：Gauss消元法是不稳定的算法，其中小主元是不稳定的根源。因而在Gauss消元法中要避免小主元的出现。这就需要采用“选主元”的技术。所谓选主元，简单地说就是选取绝对值最大的元素作主元。 全选主元Gauss消元法在此范围内选主元需进行元素比较 列选主元Gauss消元法在此范围内选主元不必选主元的情况： 当系数矩阵A对称正定或严格对角占优或不可约对角占优时，可不必选主元。现用列选主元消元法重解例3.1.3：将方程组的1,2行交换，即得回代后得到这是一个相当不错的结果！0.9999例3.1.4解线性方程组(

7、用8位十进制尾数的浮点数计算)解这个方程组和例3.1.3一样，若用Gauss消元法计算会有小数作除数的现象，若采用换行的技巧（即列选主元消元法） ，则可避免。绝对值最大不需换行经过回代后可得事实上,方程组的准确解为 列选主元Gauss消元法的计算步骤? 列选主元Gauss消元法的流程图开始输出无解信息消元换行停机回代求解本算法主要由四个模块组成 一些特殊情况, 主元就在对角线上,不需选主元. 元素满足如下条件的矩阵 即对角线上每一元素的绝对值均大于同行其他各元素绝对值之和,这样的矩阵称为按行严格对角占优矩阵,简称严格对角占优矩阵.例：性质：严格对角占优矩阵必定非奇异.上一页 下一页 返回 补充

8、一、 LU分解法 /* LU Factorization. */就 n=3的情况分析顺序消去法的消元过程.上一页 下一页 返回 矩阵的三角(LU)分解法三、直接法 可见, 消元过程相当于下述矩阵乘法运算：由分块乘法可得：直接计算可得由(*)式可得上一页 下一页 返回 Step 1:记M(1) =，则Step n 1: 其中 M (k)= 以上分析推广到n阶线性方程组可得上一页 下一页 返回 记为L单位下三角阵/* unitary lower-triangular matrix */记 U =A 的 LU 分解/* LU factorization */也称 A 的Doolittle分解上一页

9、下一页 返回 道立特分解法 /* Doolittle Factorization */： LU 分解的紧凑格式 /* compact form */根据矩阵乘法法则,先比较等式两边第1行和第1列元素有：上一页 下一页 返回 设已定出U 的第1行到第k-1行的元素 L 的第1列到第k-1列的元素上一页 下一页 返回 (1),(2)两式就是 LU分解的一般计算公式, 其结果与高斯消去法所得结果完全一样. 避免了中间过程的计算,故称为A的直接分解公式.上一页 下一页 返回 按上图逐框求出矩阵A的LU分解,这种计算方法称为紧凑格式法。上一页 下一页 返回 定理 若矩阵A非奇异, 则A能分解为LU 的充

10、分必要条件是A的所有顺序主子式 均不为0.定理 若非奇异矩阵A有LU 分解,则此分解是唯一的.上一页 下一页 返回 矩阵的Doolittle分解法的Matlab程序function 1,u=lu_Doolittle1(A)%A为可逆方阵，l为返回下三角矩阵，u为上三角阵n=length(A);u=zeros(n);l=eye(n);u(1,:)=A(1,:);l(2:n,1)=A(2:n,1)/u(1,1);for k=2:n u(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n); l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1

11、:k-1,k)/u(k,k);end例：利用系数矩阵的LU分解, 求解方程组解：LU分解的紧凑格式为上一页 下一页 返回 在Matlab命令窗口执行A=1 1 1 1;1 2 1 3;4 3 2 1;6 11 3 7;l,u=lu_Doolittle1(A)则求解原方程组等价于求解下面两个方程组Ly=b及Ux=y上一页 下一页 返回 在Matlab命令窗口执行A=1 1 1 1;1 2 1 3;4 3 2 1;6 11 3 7;b=6 12 14 45;y,x=LU_s(A,b)A=LU若U为单位上三角矩阵，则称为Crout分解.Crout分解相应的求解公式可类似得到。上一页 下一页 返回 二

12、、平方根法 对称 正定矩阵的分解法将对称 正定阵 A 做 LU 分解记为定理 设n阶对称正定矩阵A，则存在唯一的单位下三角阵L及对角阵D 使得 。 称为矩阵A 的乔累斯基分解上一页 下一页 返回 定理 设矩阵A对称正定，则存在唯一的对角元为正的下三角阵 L，使得 。 称为对称正定矩阵A 的乔累斯基分解 利用乔累斯基（Cholesky）分解式来求解Ax=b的方法也称Cholesky方法或平方根法三、三对角方程组追赶法若A满足Gauss消去法可行的条件，则可用LU分解法求解其中：上一页 下一页 返回 解方程组Ax=d分为两步，即求解Ly=d和Ux=y，计算公式如下：上述方法为求解三对角方程组的追赶

13、法，也称Thomas算法.上一页 下一页 返回 追赶法公式简单，计算量和存储量都小，整个求解过程只需要5n-4次乘除运算。 在分析方程组的解的误差及下章中迭代法的收敛性时，常产生一个问题,即如何判断向量 x 的“大小”，对矩阵也有类似的问题. 本节介绍n维向量范数和nn矩阵的范数. 向量范数是三维欧氏空间中向量长度概念的推广，在数值分析中起着重要作用.四、线性代数方程组的性态与误差分析 首先将向量长度概念推广到Rn(或Cn)中.称为向量x, y的数量积(内积). 将非负实数称为向量x的欧氏范数. 定义3.1 设x=(x1,x2,xn)T, y= (y1,y2,yn)TRn , 将实数 下面定理

14、可在线性代数书中找到. 定理 设x, yRn (或Cn), 则 1. (x, x)=0, 当且仅当x=0时成立; 我们还可以用其他办法来度量Rn中向量的“大小”. 例如, 对于x=(x1,x2)TRn, 可以用一个x的函数N(x)=max|xi|来度量 x 的“大小”, 而且这种度量 x 的“大小”的方法计算起来比欧氏范数方便. 在许多应用中, 对度量x的“大小”的函数N(x)都要求是正定的、齐次的且满足三角不等式. 下面我们给出向量范数的一般定义. (1) |x|0(|x|=0当且仅当x=0) (正定性), (2) |x|=| |x|, 对任何R(或C)(齐次性), (3) |x+y| |x

15、|+|y| (三角不等式).则称f(x)=|x|是Rn(或Cn)上的一个向量范数(或模).定义3.2(向量的范数) 如果向量xRn(或Cn)的某个实值函数f(x)=|x|，满足条件: 由(3)可推出不等式. (4) | |x|-|y| | |x-y|. 下面给出几种常用的向量范数, 设x=(x1,x2,xn)T.容易证明前三种范数是的p-范数特殊情况, 其中向量的-范数(最大范数)向量的1-范数向量的p-范数(0p0, 使得对一切xRn有 证明 只要就|x|s=|x|证明上式成立即可, 即证明存在常数c1,c20, 使得 考虑n元函数记S=x|x|=1, xRn, 则S是一个有界闭集. 由于f

16、(x)为S上的连续函数, 所以f(x)于S上达到最大值和最小值, 即存在x, xS使得设xRn且x0, 则x/|x|S, 从而有即 可以证明常用范数满足下述关系以及 注意 定理3.2不能推广到无穷维空间. 由定理3.2可得结论: 如果在一种范数意义下向量序列收敛时, 则在任何一种范数意义下该向量序列均收敛. 定理3.3 , 其中|为向量的任一种范数. 证明 显然, 而对于Rn上任一种范数|, 由定理15, 存在常数c1,c20, 使于是又有 下面我们将向量范数概念推广到矩阵上去. 可以将Rnn中的矩阵A=(aij)nn当作n2维向量, 则由向量的2-范数可以得到Rnn中矩阵的一种范数称为A的F

17、robenius范数. |A|F显然满足正定性、齐次性及三角不等式. 下面我们给出矩阵范数的一般定义. 定义3.4(矩阵的范数) 如果矩阵ARnn的某个实值函数N(A)=|A|，满足条件: (1) |A| 0 (|A|=0当且仅当A=0) (正定性)； (2) |cA|=|c| |A|, c为实数 (齐次性)； (3) 对任意A,B有|A+B| |A|+|B| (三角不等式) (4) 对任意A,B有|AB| |A| |B| (相容性)；则称 N(A) 是Rnn上的一个矩阵范数(或模). 上面我们定义的F(A)=|A|F就是Rnn上的一个矩阵范数. 上述条件(即不等式)就称为矩阵范数与向量范数的

18、相容性条件. 由于在大多数与估计有关的问题中, 矩阵和向量会同时参与讨论, 所以希望引进一种矩阵的范数, 它是和向量范数相联系而且和向量范数相容的, 即对任何向量xRn及矩阵ARnn都成立|Ax|A| |x|. 为此我们再引进一种矩阵的范数. 定义3.5(矩阵的算子范数) 设xRn, ARnn, 给出一种向量范数|x|v(如v=1,2或), 相应地定义一个矩阵的非负函数可验证|A|v满足定义4(见下面定理), 所以|A|v是Rnn上矩阵的一个范数, 称为A的算子范数. 定理3.4 设|x|v是Rn上的一个向量范数, 则|A|v是Rnn上矩阵的范数, 且满足相容条件 证明 由(5.6)相容条件(

19、5.7)是显然的. 现只验证定义4中条件(4). 由(5.7), 有当x0, 有故显然这种矩阵的范数|A|v依赖于向量范数|x|v的具体含义. 也就是说, 当给出一种具体的向量范数|x|v时, 相应地就得到了一种矩阵范数|A|v. 定理3.5 设xRn, A=(aij)Rnn, 则有常用的算子范数（称为A的行范数）,（称为A的列范数）,（称为A的2-范数）.其中max(ATA)表示ATA的最大特征值, 由于矩阵2-范数与ATA 的特征值有关，所以又称为A的谱范数. 定义3.6 设ARnn的特征值为i (i=1,2,n) , 称为A的谱半径.例子 设, 求A的谱半径.解得A的特征值故得A的谱半径

20、为 由定理3.5看出, 计算一个矩阵的|A|, |x|1还是比较容易的, 而矩阵的2-范数在|x|2计算上不方便, 但是矩阵的2-范数具有许多好的性质, 它在理论上是非常有用的. 我们指出, 对于复矩阵(即ACnn)定理3.5中1,2显然也成立, 对于3应改为例7 设则求相应的各种范数.解因为令得ATA的特征值故可见定理3.7(特征值上界) 设ARnn, 则(A)|A|,即A的谱半径不超过A的任何一种算子范数(对|A|F亦成立). 证明 设是A的任一特征值, x为相应的特征向量,则Ax=x, 由(5.7)得注意到|x|0, 即得|A|. 定理3.6 , 其中|为矩阵的任一种范数.定理3.8定理

21、3.9 如果|A|1时,则(6.3)是“病态”的(即A是“病态”矩阵, 或者说A是坏条件的, 相对于方程组), 当A的条件数相对的小, 则(6.3)是“良态”的(或者说A是好条件的). 注意, 方程组病态性质是方程组本身的特性. A的条件数越大, 方程组的病态程度越严重, 也就越难用一般的计算方法求得比较准确的解.注: cond (A) 的具体大小与 | | 的取法有关，但相对大小一致。 cond (A) 取决于A，与解题方法无关。常用条件数有：cond 1(A)cond (A)cond2 (A)条件数的性质： A可逆，则 cond p (A) 1； A可逆， R 则 cond ( A) =

22、cond (A) ； A正交，则 cond 2 (A) =1； A可逆，R正交，则 cond 2 (RA) = cond 2 (AR) = cond (A)2 。上一页 下一页 返回 通常使用的条件数, 有 (1) cond(A) =|A-1|A| ; (2) A的谱条件数;当A为对称矩阵时其中1, n为A的绝对值最大和绝对值最小的特征值. 条件数的性质: 1. 对任何非奇异矩阵, 都有cond(A)v1. 事实上 cond(cA)v=cond(A)v. 2. 设A为非奇异矩阵且c0(常数), 则 3. 如果A为正交矩阵, 则cond(A)2=1; 如果A为非奇异矩阵, R为正交矩阵, 则 c

23、ond(RA)2=cond(AR)2=cond(A)2.精确解为例计算cond 2(A) 。A1 = 解：考察 A 的特征根39206 1 测试病态程度：给一个扰动，其相对误差为此时精确解为2.0102 200%上一页 下一页 返回 例：Hilbert 阵cond (H2) =27cond (H3) 748cond (H6) =2.9 106cond (Hn) as n 注：一般判断矩阵是否病态，并不计算A1，而由经验得出。 行列式的值很大或很小（如某些行、列近似相关）； 元素间的数量级相差大，且无规则； 主元消去过程中出现小主元； 特征值相差大数量级。上一页 下一页 返回 直接法: 经过有限

24、次运算后可求得方程组精确解的方法(不计舍入误差!) 直接法比较适用于中小型方程组。直接法得到的解理论上是准确的，但它们的计算量都是n3数量级，存储量为n2数量级，这在n比较小的时候还比较合适（n400）。对于现在的很多实际问题，往往要我们求解n很大的高阶方程组，而且这些矩阵往往是稀疏的。 对于这类矩阵，用直接法时在运算中很难保持稀疏性，程序设计复杂，且会耗费大量的时间和存储单元。因此我们有必要引入一类新的方法：迭代法。四、迭代法 迭代法能保持矩阵的稀疏性，具有计算简单，存储量小，编制程序容易等优点，并在许多情况下收敛较快。故迭代法能有效地解一些高阶方程组。 迭代法的基本思想是构造一串收敛到解的

25、序列，即建立一种从已有近似解计算新的近似解的规则。由不同的计算规则得到不同的迭代法，本节介绍的几种迭代法都各有其限定的应用范围，这是因为对于一个给定的方程组，某些迭代法收敛得很快，而有些迭代法可能不收敛，或收敛得很慢，以至于无实用价值。 迭代法: 从解的某个近似值出发,通过构造一个无穷序列去逼近精确解的方法。(一般有限步内得不到精确解) 迭代法的一般形式及其收敛性四、迭代法几种古典迭代算法（一） 雅可比(Jacobi)迭代法五、迭代法建立迭代格式：记 Jacobi迭代的矩阵形式易知，Jacobi 迭代有 Jacobi迭代法的计算过程Jacobi迭代法分量形式的Matlab程序function

26、x,k=jacobif(A,b,x0,ep,Nmax)n=length(A);k=0;if nargin5 Nmax=500;endif nargin5 ep=1e-5;endif narginep&kNmax,k=k+1;x0=x; for i=1:n y(i)=b(i); for j=I y(i)=y(i)-A(i,j)*x0(j); end; end if abs(A(i,i)1e-10|k=Nmax warning(A(i,i)太小); return end y(i)=y(i)/A(i,i); end x=y;end（二） 逐次代换法（GaussSeidel迭代） GaussSeidel迭代的矩阵形式 Gauss-Seidel