数值计算方法：03线性代数方程组的数值解法

第三章第三章引言自然科学和工程技术：电学中的网络问题船体数学放样中建立三次样条函数问题用最小二乘法求实验数据的曲线拟合问题解非线性方程组问题用差分法或者有限元法解常微分方程偏微分方程边值问题等都导致求解线性方程组线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。nnnnnnnnbbbbxxxXaaaaaaaaaA2121212222111211nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111（3.1）线性方程组数值解法的分类线性方程组数值解法的分类 直接法直接法（适用于中等规模的（适用于中等规模的n n阶线性方程组）阶线性方程组）GaussGauss消去法及其变形消去法及其变形 矩阵的三角分解法矩阵的三角分解法 迭代法迭代法（适用于高阶线性方程组）（适用于高阶线性方程组）JacobiJacobi迭代法迭代法 Gauss-SeidelGauss-Seidel迭代法迭代法 逐次超松弛法逐次超松弛法 共轭斜量法共轭斜量法克莱姆法则代替所得。列用的第是，其中，法则：biAAADniDDxGrameriiiii)det(0A)det(D,.,2,1克莱姆法则在理论上有着重大意义，但在实际应用中存在很大的困难，在线性代数中，为解决这一困难给出了高斯消元法高斯消元法。1 1三角形方程组的解法三角形方程组的解法-回代法回代法nnnnnnbxaxaxabxaxabxa221122221211111 （3.2）nnnnnnnnbxabxaxabxaxaxa2222211212111（3.3）2 2顺序高斯消去法顺序高斯消去法1,3322111,223232221211,11313212111nnnnnnnnnnnnnnaxaxaxaxaaxaxaxaxaaxaxaxaxa 基本思想：基本思想：通过消元将上述方程组通过消元将上述方程组 化为三角形方程组进行求解。化为三角形方程组进行求解。nkinkjaaaankkjaaakkjkikkijkijkkkkkjkkj,11,11,1,)()1()1()()1()1()(1,1 1)()(1,)(1,nkxaaxaxnkjjkkjknkknnnn9高斯消元法v例1.用消元法解方程组)3(1022)2(543)1(614321321321xxxxxxxxx10例题v第一步：-2xr1+r3得)4(114)2(54)1(63232321xxxxxxx11例题v第二步：1 x r2+r4v回代得：x=1,2,3T)5(62)2(54)1(6332321xxxxxx12高斯顺序消去法算法框图13高斯消去法的计算量 次除法。即做kn),.,1i(/:)()(nkaalkkkkkikik步第:)1)()1()(故总的消元计算量为次乘法需由knknAAkk11)52)(1(61)1)()(nknnnknknkn)1(21)()(nnbXAnn回代时乘除运算量为解)13(312nnnN即总计算量为顺序顺序GaussGauss消去法可执行的前提消去法可执行的前提定理定理 1 1 给定线性方程组给定线性方程组 ，如果，如果n n阶方阵阶方阵 的所有顺序主子式都不为零，即的所有顺序主子式都不为零，即 则按顺序则按顺序GaussGauss消去法所形成的各主元素消去法所形成的各主元素 均不为零，从而均不为零，从而Gauss 消去法可顺利执行。消去法可顺利执行。AxbA()(1,2,)kkkakn注：注：当线性方程组的系数矩阵为对称正定或严格对角当线性方程组的系数矩阵为对称正定或严格对角占优阵时，按占优阵时，按GaussGauss消去法计算是稳定的。消去法计算是稳定的。,0,0222112112111aaaaDaD01111kkkkkaaaaD15高斯顺序消去法条件nkakkk,.,2,1,0)(高斯顺序消去法要求0.)det()()2(22)1(11nnnaaaA有：.也就是此算法的缺点.,00:)()()(即数值不稳定做除数易产生解的失真用此时很小，但若即使kkkkkkkkkaaa16vGauss列主元消元法v从第一列中选出绝对值最大的元素1111maxiniaannnnniiniinbaaabaaabaaabaaa.21212112221111211交换3、列主元、列主元Gauss消去法计算步骤：消去法计算步骤：1、输入矩阵阶数输入矩阵阶数n，增广矩阵增广矩阵 A(n,n+1);2、对于对于nk,2,1(1)按列选主元：选取按列选主元：选取 l 使使 0maxikniklkaa(2)如果如果 ，交换，交换 A(n,n+1)的第的第k行与第行与第l 行元素行元素kl(3)消元计算消元计算 ：nkiaamkkikik,1 1,1 ,1,nkjnkiamaakjikijij3、回代计算回代计算1,1,11,nnixaaxnijjijnii18框图4 4无回代过程的主元消去法无回代过程的主元消去法算法：算法：第一步：选主元，在第一列中选绝对值最大的元素，设第第一步：选主元，在第一列中选绝对值最大的元素，设第k行为主元行，行为主元行，将主元行换至第一行，将第一个方程中将主元行换至第一行，将第一个方程中x1的系数变为的系数变为1，并从，并从 其余其余个方程中消去个方程中消去x1。第二步：在第二列后第二步：在第二列后n 1个元素中选主元，将第二个方程中个元素中选主元，将第二个方程中x2的的 系数变为系数变为1，并从其它，并从其它个方程中消去个方程中消去x2。第第k步：在第步：在第k列后列后n k个元素中选主元，换行，将第个元素中选主元，换行，将第k个方程个方程xk的系数的系数 变为变为1，从其它，从其它个方程中消去变量个方程中消去变量xk，nkkinkkjaaaankkjaaakkjkikkijkijkkkkkjkkj,1,1,11,1,1,1,)()1()1()()1()1()(消元公式为：消元公式为：对对k=1,2,按上述步骤进行到第按上述步骤进行到第n步后，方程组变为：步后，方程组变为：)(1,)(1,22)(1,11nnnnnnnnaxaxax即为所求的解即为所求的解注：注：无回代的无回代的GaussGauss消元法实际上就是将消元法实际上就是将 方程组的系数矩阵化为行最简形矩阵。方程组的系数矩阵化为行最简形矩阵。22算法5无回代消去法的应用无回代消去法的应用（1）解线性方程组系解线性方程组系设要解的线性方程组系为：设要解的线性方程组系为：AX=b1,AX=b2,AX=bmmjnibbbbxxXaaaaAijnijijnnnnn,.2,1；,2,1，)(,)(,2)(,1j11111上述方程组系可以写为上述方程组系可以写为AX=B=(b1,bm)因此因此X=A-1B 即为线性方程组系的解。即为线性方程组系的解。在计算机上只需要增加几组右端常数项的存贮单元，在计算机上只需要增加几组右端常数项的存贮单元，其结构和解一个方程组时一样。其结构和解一个方程组时一样。行行n21nnnnnnaaaaaaaaa212222111211)(1,)1(1,)(1,2)1(1,2)(1,1)1(1,1bmnnnnmnnmnnbbbbb系数系数右端右端（2）求逆矩阵求逆矩阵设设A=(aij)n n是非奇矩阵是非奇矩阵，A 0，且令且令nnijxAX)(1由于由于 AA-1=AX=I因此，求因此，求A-1的问题相当于解下列线性方程组的问题相当于解下列线性方程组1010,0012112111nnnnnxxxAxxxA相当于相当于（1）中中m=n,B=I 的情形。的情形。（3）求行列式的值求行列式的值用高斯消去法将用高斯消去法将 A化成化成)()1(11)()2(22)1(11nnnnnnaaaaaA3 3 矩阵的矩阵的三角分解法三角分解法 高斯消元法的矩阵形式高斯消元法的矩阵形式 每一步消去过程相当于左乘初等变换矩阵每一步消去过程相当于左乘初等变换矩阵L Lk k(2)(1)(2)(1)11112111113111 ,1101 2 3001L()()iin i,nbbAL ALllaall记：其中(3)(2)(3)(2)222222223222 ,10101 3 4001L()()iini,nbbAL ALlaall记：(3)(1)(3)(1)2121 ,bbAL L AL L-1ii1,1,110 10111L=01010101 i iiiininiLllll 列 i列11211121(n)()nn(n)()nn AL LL AbbL LL1111121()(n)(n)nLLUAL LL AAA 的的 LU 分解分解(LU factorization)定理定理2：（矩阵的三角分解）设（矩阵的三角分解）设A为为n n实矩阵，如果实矩阵，如果解解AX=b用高斯消去法能够完成（限制不进行行的交用高斯消去法能够完成（限制不进行行的交换，即换，即 ），则矩阵），则矩阵A可分解可分解为单位下三角矩阵为单位下三角矩阵L与上三角知阵与上三角知阵U的乘积。的乘积。A=LU且这种分解是唯一的。且这种分解是唯一的。nkakkk,2,1,0)(L 为单位下三角阵而为单位下三角阵而 U 为为一般一般上三上三角阵的分解称为角阵的分解称为Doolittle 分解分解 （2）L 为一般下三角阵而为一般下三角阵而 U 为为单位单位上三上三角阵的分解称为角阵的分解称为Crout 分解分解。Doolittle分解法分解法 ：通过比较法直接导出通过比较法直接导出L 和和 U 的计算公式。的计算公式。思思路路 nnnnnnnnuuullaaaa.1.11.1111211111 ),min(1jikjkkiul jiaLU 分解求解线性方程组分解求解线性方程组 LY b,UXY AXb112221313212(1)1111121n22222nnn111 1 UXnnnnn nnnybybLYbybxyxyYxylllllluuuuuu 直接三角分解法解直接三角分解法解AX=b的计算公式的计算公式niauii,2,111niualii,2,11111对于对于r=2,3,n计算计算（2）计算计算U的第的第r行元素行元素),1,(11nrriulaurkkirkriri（3）计算计算L的第的第r 列元素列元素(r n),1()(11nriuulalrrrkkrikirir(1),3,2(111niylbybykikikii)1,2,1(1niuxuyxuyxiiknikikiinnnn(4)(5)123 12314 2521831520100123 210014.3510024 (14,18,20)(14,10,72),(14,10,72)(1,2,3).TTTTxxxALULyyUxx 例：用直角三角分解法解解：用分解计算公式得求解4 4 平方根法平方根法1 1矩阵的矩阵的LDRLDR分解分解定理定理3：如果如果n阶矩阵阶矩阵A的所有顺序主子式均不等于零，的所有顺序主子式均不等于零，则矩阵则矩阵A存在唯一的分解式存在唯一的分解式A=LDR其中其中L和和R分别是分别是n阶单位下三角阵和单位上三角阵，阶单位下三角阵和单位上三角阵，D是是n阶对角元素阶对角元素的不为零的对角阵，上述分解也称为的不为零的对角阵，上述分解也称为A的的LDR分解分解。2平方根法平方根法 如果如果A为对称正定矩阵为对称正定矩阵，则存在一个实的非奇异下三则存在一个实的非奇异下三角矩阵角矩阵，使使A=LLT ，且当限定的对角元素为正时且当限定的对角元素为正时，这种分解是唯一的这种分解是唯一的。定理定理4：（对称正定矩阵的三角分解）（对称正定矩阵的三角分解）将将对称对称 正定阵正定阵 A 做做 LU 分解分解U=uij=u11uij/uii111u22unn记为记为UD A 对称对称TUL 即即TLDLA 记记 D1/2=11u22unnu2/1LDL 则则 仍是下三角阵，且有仍是下三角阵，且有TLLA 用平方根法解线性代数方程组的算法用平方根法解线性代数方程组的算法（1）对矩阵对矩阵A进行进行Cholesky分解分解，即即A=LLT，由矩阵乘法由矩阵乘法：对于对于 i=1,2,n 计算计算21112ikikiiiilal1,2,111ijlllaljjjkikikijij（2）求解下三角形方程组求解下三角形方程组 iikikikiilylby11（3）求解求解LTX=y)1,1,(1nnilxlyxiiknikkiii3改进平方根法改进平方根法 11111112231131222111323121nnnnnllllldddllllAikdlskkikik其中其中11112121222111nnnnnlldssdsd改进平方根法解对称正定方程组的算法改进平方根法解对称正定方程组的算法 111111112,3,jijijikkjkijijjjiiiiiikikkdainsas llsddas l对 于 计 算令令LTX=y，先解下三角形方程组先解下三角形方程组 LDY=b 得得),2,1(11nidyldbyiikikikkkii解上三角形方程组解上三角形方程组 LTX=Y 得得 )1,2,1,(1nnixlyxnikkikii5 5 向量和矩阵的范数向量和矩阵的范数 1 1向量的范数向量的范数定义定义1 1：设设X R n，表示定义在表示定义在Rn上的一个实值函数上的一个实值函数，称之为称之为X的范数的范数，它具有下列性质它具有下列性质：XaaX（3）三角不等式三角不等式：即对任意两个向量即对任意两个向量X、Y R n，恒有恒有 YXYX(1)(1)非负性非负性：即对一切即对一切X R n，X 0,0(2)(2)齐次性齐次性：即对任何实数即对任何实数a R，X R n，设设X=(x1,x2,xn)T，则有则有nxxxX211（1）222212nTxxxXXX（2）inixX1max（3）三个常用的范数：三个常用的范数：范数等价范数等价:设设A A 和和B B是是R R上任意两种范数，若存在上任意两种范数，若存在 常数常数 C C1 1、C C2 2 0 0 使得使得 ,则称则称 A A 和和B B 等价等价。定理定理5：定义在定义在Rn上的向量范数上的向量范数 是变量是变量X分量的分量的 一致连续函数一致连续函数。X()Xf X定理定理6 6：在在Rn上定义的任一向量范数上定义的任一向量范数 都与范数都与范数 等价等价，即存在正数即存在正数 M 与与 m(Mm)对一切对一切X Rn，不等式不等式X1X11XMXXm成立成立。推论推论：Rn上定义的任何两个范数都是等价的。上定义的任何两个范数都是等价的。111XXXnXnXX1XnXX2对常用范数，容易验证下列不等式：对常用范数，容易验证下列不等式：定义定义2：设给定设给定Rn中的向量序列中的向量序列 ，即即kX01,X,kXX其中其中TknkkkxxxX)()(2)(1,若对任何若对任何i(i=1,2,n)都有都有*)(limikikxx则向量则向量 TnxxX),(*1*limXXkk称为向量序列称为向量序列 的极限的极限，或者说向量序列或者说向量序列 依坐标收敛于向量依坐标收敛于向量 ，记为记为kXkX*X定理定理7：向量序列向量序列Xk依坐标收敛于依坐标收敛于X*的充要条件是的充要条件是0lim*XXkk向量序列依范数收敛与依坐标收敛是等价的。向量序列依范数收敛与依坐标收敛是等价的。2 2矩阵的范数矩阵的范数定义定义3：设设A为为n 阶方阵阶方阵，Rn中已定义了向量范数中已定义了向量范数 ，则称则称 为矩阵为矩阵A A的范数或模的范数或模，记为记为 。即。即AXx1supAAXAx1sup矩阵范数的基本性质矩阵范数的基本性质：（1）当）当A=0时，时，0，当，当A 0时，时，0AA（2）对任意实数对任意实数k 和任意和任意A，有，有AkkA（3）对任意两个对任意两个n阶矩阵阶矩阵A、B有有BABA（5）对任意两个对任意两个n阶矩阵阶矩阵A、B，有有BAAB（4）对任意向量）对任意向量X Rn，和任意矩阵和任意矩阵A，有有XAAX例例5:5:设设A A(a aijij)M)M.定义定义2,11|nijijAan证明证明：这样定义的非负实数不是相容的矩阵范数这样定义的非负实数不是相容的矩阵范数.证明：设1111,1111AB2222AB|1,|1,|2ABAB从而|ABAB定理定理8：设设n 阶方阵阶方阵A=(aij)n n，则，则（）与）与 相容的矩阵范数是相容的矩阵范数是1xniijjaA11max（）与）与 相容的矩阵范数是相容的矩阵范数是2x12A其中其中 1为矩阵为矩阵ATA的最大特征值。的最大特征值。（）与）与 相容的矩阵范数是相容的矩阵范数是xnjijiaA1max上述三种范数分别称为矩阵的上述三种范数分别称为矩阵的1-1-范数、范数、2-2-范数和范数和-范数。范数。可以证明可以证明,对方阵对方阵 和和 ，有，有nnRAnx R22|FAxAx ninjijFaA112|(向量向量|2 2的直接推广的直接推广)FrobeniusFrobenius范数范数:3矩阵矩阵的范数与特征值之间的关系的范数与特征值之间的关系定理定理9：矩阵矩阵A 的任一特征值的绝对值不超过的任一特征值的绝对值不超过A的范数，即的范数，即 定义定义4：矩阵矩阵A 的诸特征值的最大绝对值称为的诸特征值的最大绝对值称为A的谱半径，的谱半径，1()maxii nA 记为：记为：1maxii nA 21max(ii nA 谱范数）并且如果并且如果A A为对称矩阵，则为对称矩阵，则 注注:R Rn nn n中的任意两个矩阵范数也是等价的。中的任意两个矩阵范数也是等价的。定义定义5 5：设设|为为R Rn nn n上的矩阵范数，上的矩阵范数，A,BA,BR Rn nn n称称|A-B|A-B|为为A A与与B B之间的距离之间的距离。定义定义6 6：设给定设给定R Rn nn n中的矩阵序列中的矩阵序列 ，若，若lim0kkAA则称矩阵序列则称矩阵序列 收敛于矩阵收敛于矩阵A A，记为，记为limkkAAkAkA定理定理1010 设设BRBRn nn n，则由，则由B B的各幂次得到的的各幂次得到的 矩阵序列矩阵序列B Bk k,k=0,1,2k=0,1,2)收敛于零矩阵收敛于零矩阵 （）的充要条件）的充要条件 为为 。()1Blim0kkB 求解求解 时，时，A 和和 的误差对解的误差对解 有何影响？有何影响？bxA bx 设设 A 精确，精确，有误差有误差 ，得到的解为，得到的解为 ，即，即bb xx bbxxA )(bAx 1|1bAx 绝对误差放大因子绝对误差放大因子|xAxAb 又又|1bAx|1bbAAxx 相对误差放大因子相对误差放大因子6 线性方程组的性态和解的误差分析线性方程组的性态和解的误差分析6 Error Analysis for .bxA 设设 精确，精确，A有误差有误差 ，得到的解为，得到的解为 ，即，即bA xx bxxAA )(bxxAxxA )()()(1xxAAx|11AAAAAAxxx bxAAxAA )()(xAxAA )(xAxAAIA )(1xAAAAIx 111)(Wait a minute Who said that(I+A 1 A)is invertible?(只要只要 A充分小，使得充分小，使得)1|11 AAAA|1|1|1111AAAAAAAAAAAAxx 是关键是关键的误差放大因子，称为的误差放大因子，称为A的的条件数条件数，记为，记为cond(A),越越 则则 A 越病态，越病态，难得准确解。难得准确解。|1 AA大大定义定义7：设设A 为为n 阶非奇矩阵，称数阶非奇矩阵，称数 为矩阵为矩阵A的条件数，的条件数，AA1条件数的性质：条件数的性质：）cond(A)1）cond(kA)=cond(A),k 为非零常数为非零常数）若）若 ，则则1A1)(cond AA记为记为cond(A)。注注:condcond(A A)与与 所取的范数有关所取的范数有关常用条件数有：常用条件数有：cond(A)2)(/)(minmaxAAAATT特别地，若特别地，若 A 对称，则对称，则|min|max)(2 Acondcond(A)1=A1 11Acond(A)=A 1A例：例：Hilbert 阵阵 12111131211211nnnnnHcond(H2)=27cond(H3)748cond(H6)=2.9 106cond(Hn)as n 注：注：现在用现在用MatlabMatlab数学软件可以很方便求矩阵的状态数数学软件可以很方便求矩阵的状态数!定义定义2:2:设线性方程组的系数矩阵是非奇异的设线性方程组的系数矩阵是非奇异的,如果如果condcond(A A)越大越大,就称这个方程组就称这个方程组越病态越病态.反之反之,condcond(A A)越小越小,就称这个方就称这个方程组程组越良态越良态.一般判断矩阵是否病态，并不计算一般判断矩阵是否病态，并不计算A A 1 1，而由经验得出。，而由经验得出。行列式很大或很小（如某些行、列近似相关）；行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；元素间相差大数量级，且无规则；主元消去过程中出现小主元；主元消去过程中出现小主元；特征值相差大数量级。特征值相差大数量级。近似解的误差估计及改善：近似解的误差估计及改善：设设 的近似解为的近似解为，则一般有，则一般有bxA*x0*xAbrbrxxx|*|cond(A)误差上限误差上限 改善方法改善方法(1)：Step 1:近似解近似解 bxA;1xStep 2:;11xAbr Step 3:;111drdA Step 4:;112dxx 若若 可被精确解出，则有可被精确解出，则有 就是精确解了。就是精确解了。1dbAxAbAxx11112)(2x经验表明经验表明：若：若 A 不是非常病态（例如：不是非常病态（例如：），），则如此迭代可达到机器精度；但若则如此迭代可达到机器精度；但若 A 病态，则此算法也病态，则此算法也不能改进。不能改进。1)(Acond 改善方法改善方法(2)(2)对方程组进行预处理，即适当选择非奇异对角阵对方程组进行预处理，即适当选择非奇异对角阵D D，C,C,使求解使求解 Ax=b Ax=b 的问题转化为求解等价方程组的问题转化为求解等价方程组 DACCDACC-1-1x=Db,x=Db,且使且使DACDAC 的条件数得到改善。（的条件数得到改善。（P88P88，例，例3.103.10）用双精度进行计算，以便改善和减轻病态矩阵的影响用双精度进行计算，以便改善和减轻病态矩阵的影响。)()()()()(1111CcondAcondDcondCCAADDDACDACDACcond