1、第第2 2章章 解线性方程组直接法解线性方程组直接法 在工程技术、自然科学和社会科学中,在工程技术、自然科学和社会科学中,经常碰到许多问题最终都可归结为解线性经常碰到许多问题最终都可归结为解线性方程组,如用最小二乘法求试验数据曲线方程组,如用最小二乘法求试验数据曲线拟合问题,工程中三次样条函数插值问题,拟合问题,工程中三次样条函数插值问题,经济运行中投入产出问题以及大地测量、经济运行中投入产出问题以及大地测量、机械与建筑结构设计计算问题等等,都归机械与建筑结构设计计算问题等等,都归结为求解线性方程组结为求解线性方程组(或非线性方程组或非线性方程组)数数学问题。所以线性方程组求解对于实际问学问题
2、。所以线性方程组求解对于实际问题是极其主要。题是极其主要。第1页 本章讨论n元线性方程组(2.1)直接解法。方程组(2.1)矩阵形式为 A Ax=b=b其中 第2页若矩阵A A非奇异,即det(A A)0,则方程组(2.1)有唯一解。所谓直接解法是指,若不考虑计算过程中舍入误差,经过有限次算术运算就能求出线性方程组准确解方法。但因为实际计算中舍入误差存在,用直接解法普通也只能求出方程组近似解。Cramer法则是一个不实用直接法,下面介绍几个实用直接法。1 Gauss(1 Gauss(高斯高斯)消去法消去法 Gauss Gauss消去法是一个规则化加减消元法,其基本思想是经过逐次消元计算,把普通
3、线性方程组求解转化为等价上三角形方程组求解。1.1 1.1 次序次序GaussGauss消去法消去法 为了清楚起见,先看一个简单例子.考虑线性方程组第3页消去后两个方程中x1得再消去最终一个方程x2得消元结束,经过回代得解:第4页上述求解消元过程可用矩阵表示为:(A,bA,b)=这是GaussGauss消去法计算形式,新增广矩阵对应线性方程组就是上三角形方程组,可进行回代求解。现在介绍求解线性方程组(2.1)次序GaussGauss消去法:记则,线性方程组(2.1)增广矩阵为第5页 第一步.设 ,依次用乘矩阵第1行加到第i行,得到矩阵:第6页其中 第二步.设 ,依次用乘矩阵第2行加到第i行,得
4、到矩阵:其中第7页如此继续消元下去,第n-1步结束后得到矩阵:这就完成了消元程。对应方程组变成:对此方程组进行回代,就可求出方程组解。第8页 次序Gauss消去法求解n元线性方程组乘除运算量是:+(n-1)2-1 +22-1 n2-1 -1+1+2+n n=20时,次序Gauss消去法只需3060次乘除法运算.次序Gauss消去法通常也简称为GaussGauss消去法消去法.次序Gauss消去法中 称为主元素主元素.主元素都不为零矩阵A A各阶次序主子式都不为零.第9页 1.2 1.2 列主元列主元GaussGauss消去法消去法 (用十进制四位浮点计算):(用Cramer法则可得准确解x1*
5、=1.00010,x2*=0.99990)解解 用次序Gauss消去法,消元得 回代得解:x2=1.00,x1=0.00 若将方程组改写成:例例1 解线性方程组第10页 用次序Gauss消去法,消元得 回代得解:x2=1.00,x1=1.00 为了提升计算数值稳定性,在消元过程中采取选择主元方法.常采取是列主元消去法列主元消去法和和全主元消去法全主元消去法.给定线性方程组AxAx=b b,记A A(1)=A A,b b(1)=b b,列主元列主元GaussGauss消去法消去法详细过程以下:首先在增广矩阵B B(1)=(A A(1),b b(1)第一列元素中,取 然后进行第一步消元得增广矩阵B
6、 B(2)=(A A(2),b b(2).再在矩阵B B(2)=(A A(2),b b(2)第二列元素中,取 第11页然后进行第二步消元得增广矩阵B(3)=(A(3),b(3).按此方法继续进行下去,经过n-1步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n).则方程组A(n)x=b(n)是与原方程组等价上三角形方程组,可进行回代求解.易证,只要|A|0,列主元Gauss消去法就可顺利进行.采取十进制四位浮点计算,分别用次序Gauss消去法和列主元Gauss消去法求解线性方程组:例例2第12页 方程组含有四位有效数字准确解为 x1*=17.46,x2*=-45.76,x3*=5.5
7、46 解解 1.1.用次序Gauss消去法求解,消元过程为回代得:x3=5.546,x2=100.0,x1=-104.0第13页 2.2.用列主元Gauss消去法求解,消元过程为第14页回代得:x3=5.545,x2=-45.77,x1=17.46 可见,列主元Gauss消去法是在每一步消元前,在主元所在一列选取绝对值最大元素作为主元素.而全主元Gauss消去法是在每一步消元前,在全部元素中选取绝对值最大元素作为主元素.但因为运算量大增,实际应用中并不经常使用.第15页2 2 矩阵三角分解法矩阵三角分解法 2.1 Gauss2.1 Gauss消去法矩阵运算消去法矩阵运算 对矩阵若 令 记第16
8、页则有 A(2)=记 令若则有第17页 A(3)=如此进行下去,第n-1步得到:A(n)=其中第18页也就是:A(n)=Ln-1A(n-1)其中 =Ln-1Ln-2A(n-2)=Ln-1Ln-2L2L1A(1)第19页所以有:A=A(1)=L1-1L2-1Ln-1-1A(n)=LU而且有其中L=L1-1L2-1Ln-1-1,U=A(n).L L为单位下三角矩阵;U U是上三角矩阵.第20页 式 A=LU称为矩阵A三角分解三角分解.2.2 Doolittle2.2 Doolittle(直接三角)分解法(直接三角)分解法 设n阶方阵A各阶次序主子式不为零,则存在唯一单位下三角矩阵L和上三角矩阵U使
9、A=LU.证实证实 只证唯一性,设有两种分解 A=LU则有 =E所以得于是 Ax=b LUx=b 令 Ux=y 得定理定理2.1第21页 下面介绍矩阵三角分解Doolittle分解方法,则得对k=2,3,n,计算 设 akj=lk1u1j+lk2u2j+lkk-1uk-1j+ukj aik=li1u1k+li2u2k+likukk第22页第23页对k=2,3,n,计算第24页由可得这就是求解方程组AxAx=b bDoolittle三角分解方法。第25页 利用三角分解方法解线性方程组 解解 因为所以例例3第26页先解 ,得再解 ,得 解线性方程组Ax=bDoolittle三角分解法计算量约为(1
10、/3)n3,与Gauss消去法基本相同.其优点在于求一系列同系数线性方程组Ax=bk,(k=1,2,m)时,可大大节约运算量.比如,求上例中矩阵A A逆矩阵.可分别取常向量 b b1=(1,0,0)T,b b2=(0,1,0)T,b b3=(0,0,1)T 第27页由所以 第28页 为了提升数值稳定性,可考虑列主元三角分解法,设已完成A=LUk-1步分解计算,矩阵分解成设 令 rkri 相当于取 为第k步分解主元素.但要注意方程组常数项也要对应变换.第29页 比如,用列主元三角分解解例3中方程组.则有第30页 设A为对称正定矩阵,则有唯一分解A=LU,且ukk0.则有 A=LDM又因为 (LD
11、M)T=MTDLT=LDM 所以 M=LT =LDLT则有 2.3 平平 方方 根根 法法第31页 分解A=GGT称为对称正定矩阵Cholesky分解分解.Ax=b 转换为 Gy=b,GTx=y 若记G G=(gij),则有:对k=1,2,n 实际计算时,可采取紧凑格式 _平方根法平方根法.第32页 解三角方程 Gy=b,GTx=y 可得 解解 例例4 解线性方程组第33页 平方根法是求对称正定系数线性方程组三角分解法,对称正定矩阵Cholesky分解计算量和存贮量均约为普通矩阵LU分解二分之一.且Cholesky分解含有数值稳定性.第34页 追赶法是求三对角线性方程组三角分解法.即方程 三对
12、角矩阵A各阶次序主子式都不为零一个充分条件是:|a1|c1|0;|an|dn|0;|ai|ci|+|di|,cidi 0,i=2,3,n-1.在此条件下,A=LDM=TM,称之为矩阵ACrout分解分解.对三对角矩阵A进行Crout分解,有2.4 2.4 追追 赶赶 法法第35页其中 解三角方程 Ty=b,Mx=y 可得称之为解三对角方程组追赶法追赶法.第36页 解解 例例5 解线性方程组第37页 当满足条件|a1|c1|0;|an|dn|0;|ai|ci|+|di|,cidi 0,i=2,3,n-1.时,追赶法是数值稳定,追赶法含有计算程序简单,存贮量少,计算量小优点.第38页3 3 向量和
13、矩阵范数向量和矩阵范数 3.1 3.1 向量范数向量范数 定义定义2.12.1 设是向量空间Rn上实值函数,且满足条件:(1)非负性非负性:对任何向量xRn,x0,且x=0当 且仅当x=0 0 (2)齐次性齐次性:对任何向量x Rn 和实数,x=|x (3)三角不等式三角不等式:对任何向量x,yRn x+yx+y则称为空间Rn上范数,x为向量x范数.第39页 记x=(x1,x2,xn)T,惯用向量范数有:向量向量1-1-范数范数:x1=|x1|+|x2|+|xn|向量向量2-2-范数范数:x2=向量向量-范数范数:x=例例6 6 设向量x=(2,-4,3,1)T,求向量范数xp,p=1,2,.
14、解解 由定义x1=10,x2=,x=4.即使不一样范数值可能不一样,但它们间存在等价关系.定理定理2.2 (范数等价性范数等价性)对于 Rn 上任何两种范数和,存在正常数m,M,使得 m x x M x,xRn第40页 惯用三种向量范数满足以下等价关系 x x1 nx,xRn 定义定义2.2 设向量序列 k=1,2,向量 假如 则称向量序列x(k)收敛于向量x*,记作 易见,第41页3.2 3.2 矩阵范数矩阵范数 定义定义2.32.3 设是以n阶方阵为变量实值函数,且满足条件:(1)非负性非负性:A0,且A=0当且仅当A=0 0 (2)齐次性齐次性:A=|A,R (3)三角不等式三角不等式:
15、A+BA+B (4)三角不等式三角不等式:ABAB则称A为矩阵A范数.记A=(aij),惯用矩阵范数有:矩阵矩阵1-1-范数范数:A1 ,也称矩阵列范数列范数.矩阵矩阵2-2-范数范数:A2 ,也称为谱范数谱范数.第42页 矩阵矩阵-范数范数:A ,也称为行范数行范数.矩阵矩阵F F-范数范数:AF 例例7 设矩阵求矩阵A范数Ap,p=1,2,F.解解 A1=4,A=5,AF 第43页 设是一个向量范数,则定义称之为由向量范数派生矩阵算子范数矩阵算子范数.矩阵算子范数满足 AxAx,xRn把满足上式矩阵范数称为与向量范数相容矩阵范数向量范数相容矩阵范数.对于p=1,2,矩阵范数Ap是由向量范数
16、xp派生矩阵算子范数,所以Ap是与xp相容矩阵范数.但AF不是一个算子范数,却与x2是相容.设是一个算子范数,则第44页 矩阵范数与矩阵特征值之间也有亲密联络.设是矩阵A特征值,x是对应特征向量,则有 Ax=x利用向量和矩阵范数相容性,则得|x=x=AxAx 于是|A 设n阶矩阵An个特征值为1,2,n,则称 为矩阵A谱半径谱半径.对矩阵任何一个相容范数都有 (A)A 另外,0,存在一个相容范数,使 A(A)+第45页 任何两种矩阵范数也含有等价性 m A A M A,ARnn 矩阵序列收敛性也定义为4 4 线性方程组固有性态与误差分析线性方程组固有性态与误差分析4.1 4.1 线性方程组固有
17、性态线性方程组固有性态考虑线性方程组其准确解为 x*=(-9800b1+9900b2 ,9900b1-10000b2)T第46页若把线性方程组变为解为x=(-9800b1+9900b2-19700 ,9900b1-10000b2+19900)T可见 x-x*=(-19700 ,19900)T 解误差可能放大到常数项误差近2万倍。若把线性方程组变为则线性方程组可能无解.这种因为原始数据微小改变而造成解严重失真线性方程组称为病态方程组病态方程组,对应系数矩阵称为病态矩阵病态矩阵.第47页 设线性方程组 Ax=b系数矩阵是准确,常数项有误差b,此时记解为x+x,则 A(x+x)=b+b于是 Ax=b
18、所以 x=A-1b A-1 b又因为 b=Ax A x所以 x b AA-1bx即第48页 再设b是准确,A有误差A,此时记解为x+x,则 (A+A)(x+x)=b则有 Ax+A(x+x)=0所以 x=-A-1A(x+x)于是 x A-1Ax+x也就是记 Cond(A)=AA-1,称为方程组Ax=b或矩阵A条件数.第49页经常使用条件数有 Condp(A)=ApA-1p p=1,2,。当A为对称矩阵时,有 Cond2(A)=|1|/|n|其中1,n分别是A按绝对值最大和最小特征值。比如,对前面方程组系数矩阵A有 Cond1(A)=Cond(A)=39601,Cond2(A)39206 因为计算
19、条件数运算量较大,实际计算中若碰到下述情况之一,方程组就有可能是病态:第50页 (1)矩阵元素间数量级差很大,且无一定规律;(2)矩阵行列式值相对来说很小;(3)列主元消去法求解过程中出现量级很小主元素;(4)数值求解过程中,计算解x剩下向量r=b-Ax已经很小,但x仍不符合要求.4.2 4.2 预条件和迭代改进预条件和迭代改进1.1.线性方程组预条件处理线性方程组预条件处理 对病态方程组Ax=b,考虑线性方程组其中第51页 称之为预条件方程组预条件方程组,显然与原方程组等价.可逆矩阵C称为预条件矩阵预条件矩阵.矩阵C应满足条件 (1)条件数 (2)方程组Cz=d轻易求解。比Cond(A)显著
20、小.对于普通矩阵A没有十分有效方法去选择预条件矩阵.当A是对称正定矩阵时,可取C .2.2.线性方程组解迭代改进线性方程组解迭代改进 设已求得方程组Ax=b近似解x(1),计算剩下向量 r(1)=b-Ax(1)再求解余量方程组Ax=r(1),得到解 则x(1)迭代改进解为:第52页练习题练习题第第50页页 习题习题22-2(1),2-6(1)(用二位浮点计算)(用二位浮点计算)第53页练习题练习题第第50页页 习题习题22-3(1),2-4,2-5,2-8第54页练习题练习题第第50页页 习题习题22-10,2-11,2-12,2-14,2-16,2-17第55页课课 堂堂 练练 习习 1.对矩阵A=进行LU分解.2.设矩阵A=,求(A)和Cond(A).解 1.因为 ,故 A=2.因为|A-E|=2-3-10=(+2)(-5)所以:1=-2,2=5,于是(A)=5.A-1=,Cond(A)=|A|A-1|=76/10=4.2 第56页课间休息休息第57页