线性方程组的数值解法公开课一等奖优质课大赛微课获奖课件.pptx

数值计算办法第二章 线性方程组数值解法第1页第1页2/114第二章 线性方程组数值解法2.0 引引 言言2.1 Gauss消去法消去法2.2 矩阵三角分解矩阵三角分解2.3 QR分解和奇异值分解分解和奇异值分解第2页第2页3/114 在在自自然然科科学学和和工工程程技技术术中中诸诸多多问问题题处处理理经经常常归归结结为为解解线线性代数方程组。性代数方程组。比如：电学中网络问题比如：电学中网络问题 用最小二乘法求试验数据曲线拟合问题用最小二乘法求试验数据曲线拟合问题 解非线性方程组问题解非线性方程组问题 用差分法或者有限元办法解常微分方程用差分法或者有限元办法解常微分方程 偏微分方程边值问题等偏微分方程边值问题等都造成求解线性代数方程组。都造成求解线性代数方程组。2.0 引引 言言第3页第3页4/114这些方程组系数矩阵大体分为两种这些方程组系数矩阵大体分为两种一个是低阶稠一个是低阶稠密矩阵密矩阵(比如，阶数大约为比如，阶数大约为150)另一个是大型稀疏矩阵另一个是大型稀疏矩阵(即矩阵阶数高且零元素较多即矩阵阶数高且零元素较多)2.0 引引 言言第4页第4页5/114设有线性方程组设有线性方程组Ax=b，其中，其中 为非奇异阵，为非奇异阵，关于线性方程组数值解法普通有两类：关于线性方程组数值解法普通有两类：直接法与迭代法。直接法与迭代法。2.0 引引 言言第5页第5页6/1141.直接法直接法 就就是是通通过过有有限限步步算算术术运运算算，可可求求得得方方程程组组准准确确解解办办法法(若计算过程中没有舍入误差若计算过程中没有舍入误差)。但但实实际际计计算算中中由由于于舍舍入入误误差差存存在在和和影影响响，这这种种办办法法也也只只能求得线性方程组近似解。能求得线性方程组近似解。本章将阐述这类算法中最基本高斯消去法及其一些变形。本章将阐述这类算法中最基本高斯消去法及其一些变形。这这类类办办法法是是解解低低阶阶稠稠密密矩矩阵阵方方程程组组有有效效办办法法，近近十十几几年年来来直直接接法法在在求求解解含含有有较较大大型型稀稀疏疏矩矩阵阵方方程程组组方方面面取取得得了了较较大进展。大进展。2.0 引引 言言第6页第6页7/1142.迭代法迭代法 就是用某种极限过程去逐步迫近线性方程组准确解办法。就是用某种极限过程去逐步迫近线性方程组准确解办法。迭迭代代法法含含有有需需要要计计算算机机存存贮贮单单元元较较少少、程程序序设设计计简简朴朴、原原始始系系数数矩矩阵阵在在计计算算过过程程中中始始终终不不变变等等长长处处，但但存存在在收收敛敛性及收敛速度问题。性及收敛速度问题。迭迭代代法法是是解解大大型型稀稀疏疏矩矩阵阵方方程程组组(尤尤其其是是由由微微分分方方程程离离散后得到大型方程组散后得到大型方程组)主要办法。主要办法。第第6 6章简介迭代法解线性方程组。章简介迭代法解线性方程组。2.0 引引 言言第7页第7页8/114 高斯高斯(Gauss)(Gauss)消去法是解线性方程组最惯用办法之一消去法是解线性方程组最惯用办法之一 它它基基本本思思想想是是通通过过逐逐步步消消元元（行行初初等等变变换换），把把方方程程组组化化为为系系数数矩矩阵阵为为三三角角形形矩矩阵阵同同解解方方程程组组，然然后后用用回回代代法法解解此三角形方程组（简朴形式）得原方程组解。此三角形方程组（简朴形式）得原方程组解。比如：比如：直接法基础直接法基础2.1 Gauss消去法消去法第8页第8页9/114 下面讨论普通解下面讨论普通解n阶方程组高斯消去法。阶方程组高斯消去法。1.消去过程消去过程 将原方程组记为将原方程组记为 A(1)(1)x=b(1)(1)其中其中A(1)(1)=(=(aij(1)(1)n n=(=(aij)n n ，b(1)(1)=b(1)(1)第一次消元。第一次消元。2.1 Gauss消去法消去法第9页第9页10/1141.消去过程消去过程(1)(1)第一次消元。第一次消元。其中其中2.1 Gauss消去法消去法注注：若若a11(1)=0，则则在在第第1列列中中至至少少有有一一个个元元素素不不为为0，可可互互换换行后再消元行后再消元第10页第10页11/114(2)(2)第第k次消元。次消元。2.1 Gauss消去法消去法第11页第11页12/1141.消去过程消去过程(2)(2)第第k次消元。次消元。注：为减少计算量，令注：为减少计算量，令 ，则则2.1 Gauss消去法消去法第12页第12页13/1141.消去过程消去过程(3)(3)当当k=n 1时得时得完毕第完毕第n 1次消元后得到与原方程组等价三角形方程组次消元后得到与原方程组等价三角形方程组A(n)x=b(n)注：当注：当det(A)0时，显然有时，显然有aii(i)0，(i=1，n)称称aii(i)为主元素。为主元素。2.1 Gauss消去法消去法第13页第13页14/1142.回代过程回代过程求解三角形方程组求解三角形方程组A(n)x=b(n)，得到求解公式得到求解公式注：求解过程称为回代过程。注：求解过程称为回代过程。2.1 Gauss消去法消去法第14页第14页15/1143.算法设计算法设计令令 i=k+1,k+2,n2.1 Gauss消去法消去法输输入入A，n，epsFor(k=1 to n-1)for(i=k+1 to n)A(i，k)=A(i，k)/A(k，k)for(j=k+1 to n+1)A(i,j)=A(i,j)-A(i,k)*A(k,j)A(n，n+1)=A(n，n+1)/A(n，n)for(k=n 1 to 1 step-1)w=0for(j=k+1 to n)w=w+A(k，j)*A(j，n+1)A(k，n+1)=A(k，n+1)wA(k，n+1)=A(k，n+1)/A(k，k)第15页第15页16/1143.算法设计算法设计function x=gauss(a)m=size(a);n=m(1);for k=1:n-1 for i=k+1:n a(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k);endendfor j=n:-1:1 a(j,n+1)=(a(j,n+1)-a(j,j+1:n)*a(j+1:n,n+1)/a(j,j)endx=a(1:n,n+1);2.1 Gauss消去法消去法输输入入A，n，epsFor(k=1 to n-1)for(i=k+1 to n)A(i，k)=A(i，k)/A(k，k)for(j=k+1 to n+1)A(i,j)=A(i,j)-A(i,k)*A(k,j)A(n，n+1)=A(n，n+1)/A(n，n)for(k=n 1 to 1 step-1)w=0for(j=k+1 to n)w=w+A(k，j)*A(j，n+1)A(k，n+1)=A(k，n+1)wA(k，n+1)=A(k，n+1)/A(k，k)第16页第16页17/1143.算法设计算法设计 i=k+1,k+2,n2.1 Gauss消去法消去法第17页第17页18/1143.算法设计算法设计 在在Excel中中1A1B1C1D1E12A2B2C2D2E23A3B3C3D3E34A4B4C4D4E35=A1=B1=C1=D1=E16=B2-B1*A2/A1=C2-C1*A2/A1=D2-D1*A2/A1=E2-E1*A2/A17=B3-B1*A3/A1=C3-C1*A3/A1=D3-D1*A3/A1=E3-E1*A3/A18=B4-B1*A4/A1=C4-C1*A4/A1=D4-D1*A4/A1=E4-E1*A4/A19=A5=B5=C5=D5=E510=B6=C6=D6=E611=C7-C6*B7/B6=D7-D6*B7/B6=E7-E6*B7/B612=C8-C6*B8/B6=D8-D6*B8/B6=E8-E6*B8/B613=A9=B9=C9=D9=E9=(E13-D13*F16-C13*F15-B13*F14)/A1314=B10=C10=D10=E10=(E14-D14*F16-C14*F15)/B1415=C11=D11=E11=(E15-D15*F16)/C1516=D12-D11*C12/C11=E12-E11*C12/C11=E16/D162.1 Gauss消去法消去法第18页第18页19/1141A1B1C1D1E12A2B2C2D2E23A3B3C3D3E34A4B4C4D4E35=A1=B1=C1=D1=E16=B2-B$1*$A2/$A$1=C2-C$1*$A2/$A$1=D2-D$1*$A2/$A$1=E2-E$1*$A2/$A$17=B3-B$1*$A3/$A$1=C3-C$1*$A3/$A$1=D3-D$1*$A3/$A$1=E3-E$1*$A3/$A$18=B4-B$1*$A4/$A$2=C4-C$1*$A4/$A$1=D4-D$1*$A4/$A$1=E4-E$1*$A4/$A$19=A5=B5=C5=D5=E510=B6=C6=D6=E611=C7-C$6*$B7/$B$6=D7-D$6*$B7/$B$6=E7-E$6*$B7/$B$612=C8-C$6*$B8/$B$6=D8-D$6*$B8/$B$6=E8-E$6*$B8/$B$613=A9=B9=C9=D9=E9=(E13-D13*F16-C13*F15-B13*F14)/A1314=B10=C10=D10=E10=(E14-D14*F16-C14*F15)/B1415=C11=D11=E11=(E15-D15*F16)/C1516=D12-D$11*$C12/$C$11=E12-E$11*$C12/$C$11=E16/D162.1 Gauss消去法消去法3.算法设计算法设计 在在Excel中中第19页第19页20/1144.Gauss消去法计算量消去法计算量以乘除法次数为主以乘除法次数为主(1)消元过程：消元过程：第第k步时步时(n k)+(n k)(n k+1)=(n k)(n k+2)共有共有2.1 Gauss消去法消去法第20页第20页21/1144.Gauss消去法计算量消去法计算量(1)消元过程：消元过程：(2)回代过程：回代过程：求求xi中中,乘乘ni次次,除除1 1次次,共共ni+1+1次次(i=1,n1)共有共有2.1 Gauss消去法消去法第21页第21页22/1144.Gauss消去法计算量消去法计算量(1)消元过程：消元过程：(2)回代过程：回代过程：(3)总次数为总次数为注注：当当n=20时时约约为为26702670次次，比比克克莱莱姆姆法法则则9.7 1020次次大大大减少。大减少。2.1 Gauss消去法消去法第22页第22页23/1145.阐明阐明(1)若消元过程中出现若消元过程中出现akk(k)=0=0，则无法继续，则无法继续(2)若若akk(k)0 0，但但较较小小，则则小小主主元元做做除除数数将将产产生生大大误误差差(3)每每次次消消元元都都选选择择绝绝对对值值最最大大者者作作主主元元，称称为为高高斯斯主主元元消去法消去法2.1 Gauss消去法消去法第23页第23页24/1145.阐明阐明(4)通常第通常第k步选择步选择第第k列列主主对对角角元元下下列列元元素素绝绝对对值值最最大大者者作作主主元元（该该行行与第与第k行对调），称为列主元消去法。行对调），称为列主元消去法。2.1 Gauss消去法消去法第24页第24页25/1145.阐明阐明【例例1】用舍入三位有效数字求解线性方程组】用舍入三位有效数字求解线性方程组（准确解是（准确解是x1=10.0，x2=1.00）解：解：1)不选主元不选主元GaussGauss消去法计算结果：消去法计算结果：x1=-10.0，x2=1.01，此解无效；此解无效；2)按按列列选选主主元元GaussGauss消消去去法法计计算算结结果果：x1=10.0，x2=1.00.2.1 Gauss消去法消去法第25页第25页26/1145.阐明阐明【例例2】求解线性方程组】求解线性方程组 和和（准确解是（准确解是x1=1.0，x2=1.0）2.1 Gauss消去法消去法第26页第26页27/1145.阐明阐明(5)行标度化行标度化 当当线线性性方方程程组组系系数数矩矩阵阵元元素素大大小小相相差差很很大大时时，则则也也许许引引进进因因丢丢失失有有效效数数字字而而产产生生误误差差，并并且且舍舍入入误误差差影影响响严严重重，即使用即使用GaussGauss主元消去法得到解也不可靠主元消去法得到解也不可靠 为为避避免免这这一一问问题题，可可将将系系数数矩矩阵阵行行元元素素按按百百分分比比增增减减以以使元素改变减小使元素改变减小2.1 Gauss消去法消去法第27页第27页28/1145.阐明阐明(5)行标度化行标度化 如如用用每每行行元元素素最最大大模模除除该该行行各各元元素素，使使它它们们模模都都小小于于1 1，这叫行标度化，其目的是要找到真正主元，这叫行标度化，其目的是要找到真正主元 消消元元过过程程仍仍是是对对原原方方程程组组进进行行，只只但但是是在在GaussGauss列列主主元元消去法算法中，按列选主元消去法算法中，按列选主元ck时，应修改为时，应修改为其中其中2.1 Gauss消去法消去法第28页第28页29/1146.列主元消去法算法列主元消去法算法输输入入A，n，epsFor(k=1 to n-1)选主元A(I0，k)拟定行号p=A(I0，k)若若|p|=epsText=1，退出循，退出循环环若若I0 kT换换行（行（I0 k）消元消元若若|A(n，n)|A(I0，k)|TI0=ifor(i=k+1 to n)A(i，k)=A(i，k)/A(k，k)for(j=k+1 to n+1)A(i,j)=A(i,j)-A(i,k)*A(k,j)A(n，n+1)=A(n，n+1)/A(n，n)for(k=n 1 to 1 step-1)w=0for(j=k+1 to n)w=w+A(k，j)*A(j，n+1)A(k，n+1)=A(k，n+1)wA(k，n+1)=A(k，n+1)/A(k，k)for(j=k to n)n=A(k，j),A(k，j)=A(I0，j),A(I0，j)=n2.1 Gauss消去法消去法第29页第29页30/1146.列主元消去法算法列主元消去法算法function x=gaussa(a)m=size(a);n=m(1);x=zeros(n,1);for k=1:n-1 c,i=max(abs(a(k:n,k);q=i+k-1;if q=k d=a(q,:);a(q,:)=a(k,:);a(k,:)=d end for i=k+1:n a(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k)endendfor j=n:-1:1 x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n)/a(j,j)end2.1 Gauss消去法消去法输输入入A，n，epsFor(k=1 to n-1)选主元A(I0，k)拟定行号p=A(I0，k)若若|p|=epsText=1，退出循，退出循环环若若I0 kT换换行（行（I0 k）消元消元若若|A(n，n)|行标行标k时时,lkm=0=0 j=k,k+1,.,nk=1,2,.,n2.2 矩阵三角分解矩阵三角分解第35页第35页36/1142.2.1 LU分解和分解和LDU分解分解2.LU分解分解设设A各阶顺序主子式不为各阶顺序主子式不为0，且，且A=LU，即，即(2)(2)主对角线下边：当行标主对角线下边：当行标m 列标列标k时时,umk=0 i=k+1,k+2,.,nk=1,2,.,n1欲求欲求lik和和ukj2.2 矩阵三角分解矩阵三角分解第36页第36页37/1142.2.1 LU分解和分解和LDU分解分解2.LU分解分解(1)(1)主对角线主对角线(含含)上边上边:当列标当列标m行标行标k时时,lkm=0 j=k,k+1,.,nk=1,2,.,n(2)(2)主对角线下边：当行标主对角线下边：当行标m 列标列标k时，时，umk=0 i=k+1,k+2,.,nk=1,2,.,n1设设k=1,=1,a1 1j=l1111u1 1j u1 1j=a1 1j j=1,2,.,n ai1 1=li1 1u1111 li1 1=ai1 1/u1111 i=2,3,.,n欲求欲求lik和和ukj2.2 矩阵三角分解矩阵三角分解第37页第37页38/1142.2.1 LU分解和分解和LDU分解分解2.LU分解分解普通地，普通地，最后：最后：u11u12.u1n第第1步步l21u22.u2n第第2步步l31l32.ln1ln2unn第第n步步2.2 矩阵三角分解矩阵三角分解第38页第38页39/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组下三角下三角 Ly=b：上三角上三角 Ux=y：2.2 矩阵三角分解矩阵三角分解第39页第39页40/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组【例【例3】用杜丽特尔法解方程组】用杜丽特尔法解方程组解：解：2.2 矩阵三角分解矩阵三角分解第40页第40页41/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组【例【例3】用杜丽特尔法解方程组】用杜丽特尔法解方程组解：解：2.2 矩阵三角分解矩阵三角分解第41页第41页42/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组【例【例3】用杜丽特尔法解方程组】用杜丽特尔法解方程组解：解：2.2 矩阵三角分解矩阵三角分解第42页第42页43/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组【例【例3】用杜丽特尔法解方程组】用杜丽特尔法解方程组解：解：2.2 矩阵三角分解矩阵三角分解第43页第43页44/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组【例【例3】用杜丽特尔法解方程组】用杜丽特尔法解方程组解：解：2.2 矩阵三角分解矩阵三角分解第44页第44页45/1142.2.1 LU分解和分解和LDU分解分解3.求解方程组求解方程组【例【例3】用杜丽特尔法解方程组】用杜丽特尔法解方程组解：解：2.2 矩阵三角分解矩阵三角分解第45页第45页46/1142.2.1 LU分解和分解和LDU分解分解4.紧凑格式紧凑格式2.2 矩阵三角分解矩阵三角分解第46页第46页47/1142.2.1 LU分解和分解和LDU分解分解5.代码代码function x=lua(a,b)m=size(a);n=m(1);x=zeros(n,1);y=zeros(n,1);for k=1:n a(k,k:n)=a(k,k:n)-a(k,1:k-1)*a(1:k-1,k:n)a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k)/a(k,k)endU=triu(a,0)L=eye(n)+tril(a,-1)for i=1:n y(i)=b(i)-L(i,1:i-1)*y(1:i-1)endfor i=n:-1:1 x(i)=(y(i)-U(i,i+1:n)*x(i+1:n)/U(i,j)end2.2 矩阵三角分解矩阵三角分解第47页第47页48/1142.2.1 LU分解和分解和LDU分解分解5.代码代码紧凑格式：紧凑格式：function x=luj(a)m=size(a);n=m(1);x=zeros(n,1)for k=1:n a(k,k:n+1)=a(k,k:n+1)-a(k,1:k-1)*a(1:k-1,k:n+1)a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k)/a(k,k)endU=triu(a,0)for j=n:-1:1 x(j)=(a(j,n+1)-U(j,j+1:n)*x(j+1:n)/U(j,j)end2.2 矩阵三角分解矩阵三角分解第48页第48页49/11457910266810918710892257659=A1=B1=C1=D1=E1=(E5-(B5*F6+C5*F7+D5*F8)/A5=A2/A$5=B2-$A6*B5=C2-$A6*C5=D2-$A6*D5=E2-$A6*E5=(E6-(C6*F7+D6*F8)/B6=A3/A$5=(B3-A7*B$5)/B$6=C3-($A7*C5+$B7*C6)=D3-($A7*D5+$B7*D6)=E3-($A7*E5+$B7*E6)=(E7-D7*F8)/C7=A4/A$5=(B4-A8*B$5)/B$6=(C4-A8*C$5-B8*C$6)/C$7=D4-($A8*D5+$B8*D6+$C8*D7)=E4-($A8*E5+$B8*E6+$C8*E7)=E8/D82.2 矩阵三角分解矩阵三角分解2.2.1 LU分解和分解和LDU分解分解6.Excel实现实现第49页第49页50/1142.2 矩阵三角分解矩阵三角分解2.2.1 LU分解和分解和LDU分解分解6.Excel实现实现第50页第50页51/1142.2.1 LU分解和分解和LDU分解分解7.LDU分解分解设设A=LU，令，令D=diag(u11,u22,.,unn)，则，则U1=D-1U仍是一个单位上三角阵仍是一个单位上三角阵故矩阵三角分解除了故矩阵三角分解除了杜丽特尔分解，尚有下列：杜丽特尔分解，尚有下列：(1)克克劳劳特特分分解解：A=LU，L为为下下三三角角阵阵，U为为单单位位上上三三角阵角阵(2)LDU分分解解：A=LDU，L为为单单位位下下三三角角阵阵，D为为对对角角阵阵，U为单位上三角阵为单位上三角阵上述分解均含有唯一性！上述分解均含有唯一性！2.2 矩阵三角分解矩阵三角分解第51页第51页52/1142.2.2 乔列斯基分解乔列斯基分解1.原理原理定理定理：若：若A对称对称，则，则A可唯一分解为可唯一分解为A=LDLT T其中其中L为单位下三角，为单位下三角，D为元素全是正数对角阵为元素全是正数对角阵证证实实：设设A=LDU，L为为单单位位下下三三角角阵阵，D为为对对角角阵阵，U为为单位上三角阵。单位上三角阵。由对称性由对称性LDU=UTDLT，由唯一性知，由唯一性知L=UT，因此，因此A=LDLT T2.2 矩阵三角分解矩阵三角分解第52页第52页53/1142.2.2 乔列斯基分解乔列斯基分解1.原理原理定理定理：若：若A对称对称，则，则A可唯一分解为可唯一分解为A=LDLT T其中其中L为单位下三角，为单位下三角，D为元素全是正数对角阵为元素全是正数对角阵证实：设证实：设A=LDLT，D=diag(d1,d2,.,dn)对于对于ei=(0,.,0,1,0,.,0)T,存在存在xi 0,使得使得LTxi=ei另外另外,xiTAxi=xiT(LDLT)xi =(LTxi)TD(LTxi)=eiTDei=di由于由于A正定，有正定，有di=xiTAxi 0 2.2 矩阵三角分解矩阵三角分解第53页第53页54/1142.2.2 乔列斯基分解乔列斯基分解1.原理原理定理定理：若：若A对称对称，则，则A可唯一分解为可唯一分解为A=LDLT T其中其中L为单位下三角，为单位下三角，D为元素全是正数对角阵为元素全是正数对角阵乔列斯基分解乔列斯基分解:若若A正定正定，则，则A可唯一分解为可唯一分解为A=GGT T其中其中G为下三角阵为下三角阵记记则有则有A=LDLT=LD1/2D1/2LT=(LD1/2)(LD1/2)T=GGT2.2 矩阵三角分解矩阵三角分解第54页第54页55/1142.2.2 乔列斯基分解乔列斯基分解2.算法算法设设A=GGT则有：则有：从而从而2.2 矩阵三角分解矩阵三角分解第55页第55页56/1142.2.2 乔列斯基分解乔列斯基分解2.算法算法设设A=GGT则有：则有：注：由上式看出注：由上式看出 ，因此，因此 即算法是稳定即算法是稳定2.2 矩阵三角分解矩阵三角分解第56页第56页57/1142.2.2 乔列斯基分解乔列斯基分解2.算法算法Ax=b GGTx=b Gy=b，GTx=y求解公式：求解公式：上述办法称为上述办法称为“平方根法平方根法”注：不能选主元作行互换注：不能选主元作行互换破坏对称性破坏对称性2.2 矩阵三角分解矩阵三角分解第57页第57页58/1142.2.2 乔列斯基分解乔列斯基分解【例【例4】用平方根法解方程组】用平方根法解方程组解解:2.2 矩阵三角分解矩阵三角分解第58页第58页59/1142.2.2 乔列斯基分解乔列斯基分解3.代码代码function x=choleskey(a,b)m=size(a);n=m(1);x=zeros(n,1);y=zeros(n,1);G=zeros(m)for j=1:n G(j,j)=sqrt(a(j,j)-G(j,1:j-1)*G(j,1:j-1)G(j+1:n,j)=(a(j+1:n,j)-G(j+1:n,1:j-1)*G(j,1:j-1)/G(j,j)endfor i=1:n y(i)=(b(i)-G(i,1:i-1)*y(1:i-1)/G(i,i)endfor i=n:-1:1 x(i)=(y(i)-G(i+1:n,i)*x(i+1:n)/G(i,i)end2.2 矩阵三角分解矩阵三角分解第59页第59页60/1142.2.2 乔列斯基分解乔列斯基分解4.ExcelExcel中计算中计算2.2 矩阵三角分解矩阵三角分解第60页第60页61/1142.2.3 追赶法追赶法在实际问题中，经常碰到下列形式方程组在实际问题中，经常碰到下列形式方程组这种方程组系数矩阵称为三对角矩阵。这种方程组系数矩阵称为三对角矩阵。下下列列针针对对该该方方程程组组特特点点提提供供一一个个简简便便有有效效算算法法 追追赶赶法法2.2 矩阵三角分解矩阵三角分解第61页第61页62/1142.2.3 追赶法追赶法作克劳特分解，作克劳特分解，A=LU，易知，易知L，U具下列形式：具下列形式：比较第比较第i行元素得行元素得2.2 矩阵三角分解矩阵三角分解第62页第62页63/1142.2.3 追赶法追赶法作克劳特分解，作克劳特分解，A=LU，易知，易知L，U具下列形式：具下列形式：计算过程：计算过程：l1 1u1 1l2 2u2 2.ln2.2 矩阵三角分解矩阵三角分解第63页第63页64/1142.2.3 追赶法追赶法有有了了ALU分分解解后后，线线性性方方程程组组Ax=d等等价价于于两两个个简简朴朴方方程程组：组：Ly=d，Ux=y Ly=f 即即2.2 矩阵三角分解矩阵三角分解第64页第64页65/1142.2.3 追赶法追赶法有有了了ALU分分解解后后，线线性性方方程程组组Ax=d等等价价于于两两个个简简朴朴方方程程组：组：Ly=d，Ux=y Ux=y即即2.2 矩阵三角分解矩阵三角分解第65页第65页66/1142.2.3 追赶法追赶法分解公式是分解公式是计算计算y公式是公式是计算过程：计算过程：l1 1u1 1y1 1l2 2u2 2y2 2.lnyn (追追)2.2 矩阵三角分解矩阵三角分解第66页第66页67/1142.2.3 追赶法追赶法分解公式是分解公式是计算计算y公式是公式是计算过程：计算过程：2.2 矩阵三角分解矩阵三角分解第67页第67页68/1142.2.3 追赶法追赶法分解公式是分解公式是计算计算x公式是公式是(赶赶)2.2 矩阵三角分解矩阵三角分解第68页第68页69/1142.2.3 追赶法追赶法2.2 矩阵三角分解矩阵三角分解第69页第69页70/1142.2.3 追赶法追赶法【例【例5】用追赶法解】用追赶法解解：追解：追2.2 矩阵三角分解矩阵三角分解第70页第70页71/1142.2.3 追赶法追赶法【例【例5】用追赶法解】用追赶法解解：赶解：赶2.2 矩阵三角分解矩阵三角分解第71页第71页72/1142.2.3 追赶法追赶法2.2 矩阵三角分解矩阵三角分解输入a,b,c,dl1=b1,y1=d1/l1i=1 to n 1ui=ci/lili+1=bi+1 ai+1uiyi+1=(di+1 ai+1yi)/li+1xn=yni=n 1 to 1xi=yi uixi+1第72页第72页73/1142.2.3 追赶法追赶法function x=chase(a,b,c,d)m=size(b);n=m(2);L=linspace(0,0,n);U=L;x=L;y=L;L(1)=b(1);y(1)=d(1)/L(1);for k=1:n-1 U(k)=c(k)/L(k)L(k+1)=b(k+1)-a(k)*U(k)y(k+1)=(d(k+1)-a(k)*y(k)/L(k+1)endx(n)=y(n)for i=n-1:-1:1 x(i)=y(i)-U(i)*x(i+1)end2.2 矩阵三角分解矩阵三角分解输入a,b,c,dl1=b1,y1=d1/l1i=1 to n 1ui=ci/lili+1=bi+1 ai+1uiyi+1=(di+1 ai+1yi)/li+1xn=yni=n 1 to 1xi=yi uixi+1第73页第73页74/1142.2.3 追赶法追赶法2.2 矩阵三角分解矩阵三角分解输入a,b,c,dd1=d1/b1i=1 to n 1ci=ci/bibi+1=bi+1 ai+1cidi+1=(di+1 ai+1di)/bi+1i=n 1 to 1di=di cidi+1输入a,b,c,dl1=b1,y1=d1/l1i=1 to n 1ui=ci/lili+1=bi+1 ai+1uiyi+1=(di+1 ai+1yi)/li+1xn=yni=n 1 to 1xi=yi uixi+1第74页第74页75/1142.2.3 追赶法追赶法function d=chase(a,b,c,d)m=size(b);n=m(2);d(1)=d(1)/b(1);for k=1:n-1 c(k)=c(k)/b(k);b(k+1)=b(k+1)-a(k)*c(k);d(k+1)=(d(k+1)-a(k)*d(k)/b(k+1);endfor i=n-1:-1:1 d(i)=d(i)-c(i)*d(i+1);end2.2 矩阵三角分解矩阵三角分解输入a,b,c,dd1=d1/b1i=1 to n 1ci=ci/bibi+1=bi+1 ai+1cidi+1=(di+1 ai+1di)/bi+1i=n 1 to 1di=di cidi+1第75页第75页76/1142.2.3 追赶法追赶法2.2 矩阵三角分解矩阵三角分解第76页第76页77/114定定理理：若若A为为对对角角占占优优(diagonally(diagonally dominant)dominant)三三对对角角阵阵，且满足且满足|b1|c1|0，|bi|ai|+|ci|，aici 0，i=2，3，n-1|bn|an|0.则追赶法可解以则追赶法可解以A为系数矩阵方程组为系数矩阵方程组注：注：1.1.假假如如A是是严严格格对对角角占占优优阵阵，则则不不要要求求三三对对角角线线上上所所有有元素非零。元素非零。2.2 矩阵三角分解矩阵三角分解第77页第77页78/114注：注：1.1.假假如如A是是严严格格对对角角占占优优阵阵，则则不不要要求求三三对对角角线线上上所所有有元素非零。元素非零。2.2.依据不等式依据不等式可可知知：分分解解过过程程中中，矩矩阵阵元元素素不不会会过过度度增增大大，算算法法确确保保稳稳定。定。3.3.运算量为运算量为 O(6 O(6n)。2.2 矩阵三角分解矩阵三角分解第78页第78页79/1142.3 QR分解和奇异值分解 矩矩阵阵分分解解是是将将矩矩阵阵分分解解为为数数个个含含有有特特殊殊性性质质矩矩阵阵因因子子乘乘积积.除除了了三三角角分分解解以以外外，尚尚有有本本节节要要简简介介QR分分解解（QR Factorization）和和 奇奇 异异 值值 分分 解解 法法（Singular Value Decompostion）.第79页第79页80/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵【定义【定义2.3.1】若矩阵】若矩阵QRn n,且满足且满足QQT=QTQ=I则称矩阵则称矩阵Q为正交矩阵为正交矩阵正交矩阵正交矩阵Q有下列性质：有下列性质：Q-1=QT;det(Q)=1;Qx长度与长度与x长度相等长度相等.下面简介几类特殊正交矩阵下面简介几类特殊正交矩阵.第80页第80页81/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵1.置换矩阵置换矩阵将将单单矩矩阵阵任任意意两两行行(列列)互互换换得得到到矩矩阵阵，称称为为置置换换矩矩阵阵.譬譬如如,将单位矩阵第将单位矩阵第i行和第行和第j列互换列互换,得到置换矩阵得到置换矩阵Pij：任意个置换矩阵乘积仍然是置换矩阵任意个置换矩阵乘积仍然是置换矩阵.第81页第81页82/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）对于某个角度对于某个角度，记，记s=sin,c=cos 那么，那么，是一个正交矩阵是一个正交矩阵.记记w=(x,y)T为为二二维维平平面面中中一一个个向向量量，用用极极坐坐标标表表示示为为w=(rsin,rcos)T.那么那么即即Gw表表示示将将向向量量w顺顺时时针针旋旋转转 角角所所得得到到向向量量，如如图图2-1所所表示表示.第82页第82页83/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）推广到推广到n n情形，形如情形，形如矩矩阵阵称称为为Givens矩矩阵阵或或Givens变变换换，或或称称(平平面面)旋旋转转矩矩阵阵（旋转变换），其中（旋转变换），其中 为旋转角度为旋转角度.显然，显然，G(i,j,)也是正交矩阵也是正交矩阵.第83页第83页84/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）若若xRn n，y=G(i,j,)x，则，则y分量为分量为假如要使假如要使yj=0只要选择只要选择 满足满足即可即可第84页第84页85/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）【例【例6】用】用Givens变换将海森伯格变换将海森伯格(Hdssenberg)型矩阵型矩阵化为上三角矩阵。化为上三角矩阵。解：首先消去解：首先消去A(2,1)，结构，结构Givens变换变换G(1,2,),其中其中第85页第85页86/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）【例【例6】用】用Givens变换将海森伯格变换将海森伯格(Hdssenberg)型矩阵型矩阵化为上三角矩阵。化为上三角矩阵。解：从而解：从而G(1,2,)=第86页第86页87/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）【例【例6】用】用Givens变换将海森伯格变换将海森伯格(Hdssenberg)型矩阵型矩阵化为上三角矩阵。化为上三角矩阵。解：解：A1=G(1,2,)A=第87页第87页88/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）【例【例6】用】用Givens变换将海森伯格变换将海森伯格(Hdssenberg)型矩阵型矩阵化为上三角矩阵。化为上三角矩阵。解解：另另一一方方面面消消去去A1(3,2),结结构构Givens变变换换G(2,3,)，其其中中第88页第88页89/1142.3 QR分解和奇异值分解2.3.1 正交矩阵正交矩阵2.旋转矩阵（旋转矩阵（Givens变换）变换）【例【例6】用】用Givens变换