第2讲线性方程组的数值方法

线性代数方程组的数值解法(1)Gauss消去法bxAnnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111(DemosinMatlab:airfoilin2D)线性代数方程组的数值解法直接法：Gauss消去法，SuperLU迭代法：定常迭代（Jacobi,GS,SOR,SSOR）Krylov子空间方法（CG,MINRES,GMRES,QMR,BiCGStab）TheLandscapeofAx=bSolversPivotingLUGMRES,QMR,…CholeskyConjugategradientDirectA=LUIterativey’=AyNon-symmetricSymmetricpositivedefiniteMoreRobustLessStorageMoreRobustMoreGeneral刘徽（约220-280）Gaussianelimination,whichfirstappearedinthetextNineChaptersontheMathematicalArtwrittenin200BC,wasusedbyGaussinhisworkwhichstudiedtheorbitoftheasteroidPallas.UsingobservationsofPallastakenbetween1803and1809,Gaussobtainedasystemofsixlinearequationsinsixunknowns.GaussgaveasystematicmethodforsolvingsuchequationswhichispreciselyGaussianeliminationonthecoefficientmatrix.(TheMacTutorHistoryofMathematics,)Gauss（1777-1855）263234323923zyxzyxzyx今有上禾三秉，中禾二秉，下禾一秉，实三十九斗；上禾二秉，中禾三秉，下禾一秉，实三十四斗；上禾一秉，中禾二秉，下禾三秉，实二十六斗。问上、中、下禾实一秉各几何？答曰：上禾一秉九斗四分斗之一。中禾一秉四斗四分斗之一。下禾一秉二斗四分斗之三。-------《九章算术》一个两千年前的例子263213413239123130803912338343135437,417,411xyz263234323923zyxzyxzyx138392338343135zyzyzyx533512313583923zzyzyx5335123135008039123Basicidea:AddmultiplesofeachrowtolaterrowstomakeAuppertriangular一个两千年前的例子(2)nnnnnnaaaaaaaaa212222111211bxAnnbbbxxx2121Solvinglinearequationsisnottrivial.Forsythe(1952)AaaaaaaaaaAnnnnnn)1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11)1()2()2(2)2(2)2(22)1(1)1(12)1(11)2(00nnnnnaaaaaaaA)3()3(3)3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11)3(00000nnnnnnaaaaaaaaaaaA),,3,2,(),,3,2()1(11)1()2()1(11)1(11njialaaniaaljiijijii),,4,3,(),,4,3()2(22)2()3()2(22)2(22njialaaniaaljiijijii),,1,(),,1()()1()1()()(nkjialaankiaalkkjikkijkijkkkkikik)()(1,)()(,1)(1,1)(,1)()(1,)()1(,1)1(1,1)1(,1)1(1,1)1(1)1(1,1)1(1)1(1,1)1(11)(000000knnkknknkknkkkkkkkkknkkkkkkknkkkkkkkkkknkkkkaaaaaaaaaaaaaaaaaaA)1()1(1,)1(,1)1(1,1)()(1,)()1(,1)1(1,1)1(,1)1(1,1)1(1)1(1,1)1(1)1(1,1)1(11)1(00000000knnkknknkkkkkknkkkkkkknkkkkkkkkkknkkkkaaaaaaaaaaaaaaaaA)()3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11nnnnnnaaaaaaaaaa)()3(3)2(2)1(1321nnnbbbbxxxxAfterk=1Afterk=2Afterk=3Afterk=n-1Gauss消去过程图示)1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11)1(nnnnnnaaaaaaaaaA)2()2(2)2(2)2(22)1(1)1(12)1(11)2(00nnnnnaaaaaaaA)3()3(3)3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11)3(00000nnnnnnaaaaaaaaaaaA111)1(11def)1(~AcraAAT222)2(221)1(11def11)1(11)1(1)2(~0~0)1(1111AcraraAraALATTarcTT22)2(221)1(11)2(2)3(~00)2(2222AraraALAarcTTTTnacelIIL22def22)2(22211TnacelIIL11def11)1(1111用矩阵变换表达消去过程UaaaaaaALLLAnnnnnnndef)()2(2)2(22)1(1)1(12)1(11)1(121)(LUULLLAndef111211利用Gauss变换矩阵的性质：TnnTTnelelelILLLL112211111211)()2()1(11jielelIelIelILLelIelILTjjTiiTjjTiijiTiiTiii用矩阵变换表达消去过程(2)单位下三角形1111)3(33)3(3)2(22)2(2)1(11)1(1)2(22)2(32)1(11)1(31)1(11)1(21aaaaaaaaaaaannn用Gauss消去法求解Ax=b---LU分解A=LU(cost=2/3n3flops)---求解Ly=b(cost=n2flops)---求解Ux=y(cost=n2flops)版本一fork=1ton-1fori=k+1tonm=A(i,k)/A(k,k)forj=ktonA(i,j)=A(i,j)-m*A(k,j)算法实现：GaussEliminationAlgorithmfork=1ton-1…对第k列,消去对角线以下元素…(通过每行加上第k行的倍数)fori=k+1ton…对第k行以下的每一行iforj=kton…第k行的倍数加到第i行A(i,j)=A(i,j)-(A(i,k)/A(k,k))*A(k,j)版本二:在内循环中去掉常量A(i,k)/A(k,k)的计算上一版本GaussEliminationAlgorithm(2)fork=1ton-1fori=k+1tonm=A(i,k)/A(k,k)forj=ktonA(i,j)=A(i,j)-m*A(k,j)fork=1ton-1fori=k+1tonm=A(i,k)/A(k,k)forj=k+1tonA(i,j)=A(i,j)-m*A(k,j)版本三:第k列对角线以下为0,无需计算上一版本GaussEliminationAlgorithm(3)fork=1ton-1fori=k+1tonm=A(i,k)/A(k,k)forj=k+1tonA(i,j)=A(i,j)-m*A(k,j)fork=1ton-1fori=k+1tonA(i,k)=A(i,k)/A(k,k)forj=k+1tonA(i,j)=A(i,j)-A(i,k)*A(k,j)版本四:将乘子m存储在对角线以下备用上一版本fork=1ton-1fori=k+1tonA(i,k)=A(i,k)/A(k,k)forj=k+1tonA(i,j)=A(i,j)-A(i,k)*A(k,j)fork=1ton-1fori=k+1tonA(i,k)=A(i,k)/A(k,k)fori=k+1tonforj=k+1tonA(i,j)=A(i,j)-A(i,k)*A(k,j)版本五:SplitloopGaussEliminationAlgorithm(4)上一版本版本六:用矩阵运算fork=1ton-1fori=k+1tonA(i,k)=A(i,k)/A(k,k)fori=k+1tonforj=k+1tonA(i,j)=A(i,j)-A(i,k)*A(k,j)GaussEliminationAlgorithm(5)fork=1ton-1A(k+1:n,k)=A(k+1:n,k)/A(k,k)…BLAS1(scaleavector)A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n)…BLAS2(rank-1update)Whatwehaven’ttoldyou选主元策略(当主元A(k)(k,k)为0或很小时)向后误差分析并行技术块算法SparseLU,BandLU最新进展(F.Gustavson&S.Toledo,RecursiveAlgorithm)还可用于矩阵求逆,求行列式,秩定理:主元A(k)(k,k)不为0的充要条件是顺序主子矩阵非奇异定理:分解的存在性和唯一性Matlab中的相应函数Linpack中对应的函数invlu\C,Fortran,Matlab代码sgea.fsgefa.ffunctionx=lsolve(A,b)%x=lsolve(A,b)returnsthesolutiontotheequationAx=b,%whereAisann-by-bmatrixandbisacolumnvectorof%lengthn(oramatrixwithsever