第五章解线性代数方程组

求解线性方程组的直接解法5.3特殊矩阵的三角分解①实对称矩阵的LDLT分解设A是实对称阵,且A的所有顺序主子式均不为零，则LDR分解中R=LT,故可用以作LDLT分解.这就是说,当A的对角元素非零时，我们可以作LU分解,也就得到LDLT分解,L相同,是单位上三角阵,U的对角元素构成D.不过没有利用对称性,存储量运算量都未能节省—预计是一半。试用n=3的计算表格说明如何实现节省。d1=u11=a11u12=a12l21=u12/d1u13=a13l31=u13/d1d2=u22=a22-l21u12u23=a23-l21u13l32=u23/d2u33=a33-l31u13-l32u23这样,可用上半部元素逐列计算D,LT。也可用下半部元素逐行计算L,D。引进輔助量t1,t2代替u1j,u2j,并利用对称性得到：d1=a11t1=a21l21=t1/d1d2=a22-t1l21t1=a31l31=t1/d1t2=a32-t1l21l32=t2/d2d3=a33-t1l31-t2l32据此不难写出LDLT分解A=LDLT的计算公式和程序(逐行计算L,D).d1=a11fori=2:nforj=1:i-1tj=aij-lj1t1-lj2t2-…-lj,j-1tj-1lij=tj/djenddi=aii-li1t1-li2t2-…-li,i-1ti-1end存储约n(n+1)/2单元,乘加运算各约n3/6.利用LDLT分解解Ax=b分四步:1．分解A=LDLT2．解Lg=b求g3．解Dy=g求y4．解LTx=y求x②实对称正定矩阵的LLT分解A实对称正定时顺序主子式皆正,可作LDLT,D的对角元素皆正,有正的平方根。因此有LLT分解A=LLT,L下三角阵,对角元素皆正,是LDLT中的LD1/2.乃可用上半部元素逐列计算LT.l11=a111/2l21=a12/l11l31=a13/l11l22=(a22-l212)1/2l32=(a23-l21l31)/l22l33=a33-l312-l322也可用下半部元素逐行计算L.计算表格和算法安排如下:l11=a111/2l21=a21/l11l22=(a22-l212)1/2l31=a31/l11l32=(a32-l31l21)/l22l33=(a33-l312-l322)1/2l11=a111/2fori=2:nforj=1:i-1lij=(aij-li1lj1-li2lj2-…-li,j-1lj,j-1)/djjend2/121,2221)(iiiiiiiilllalend存储量,运算量同LDLT分解,但要n次求平方根.利用LLT分解解Ax=b分三步:1．分解A=LLT2．解Lg=b求g3．解LTx=g求x③三对角方程组的追赶法消去法或LU分解用于三对角方程组有特殊形式,即称追赶法.设Ax=f:b1x1+c1x2=f1aixi-1+bixi+cixi+1=fii=2,3,n-1anxn-1+bnxn=fnA是三对角阵,则L,U同样结构.L的对角元素为α2,α3,…，αn,U的对角元素为β1,β2,…,βn,上对角元素同A.1．分解A=LU:β1=b1,αi=ai/βi-1,βi=bi-αici-1,i=2,3,…,n2．解Lg=f求g:g1=f1,gi=fi-αifi-1,i=2,3,…,n3．解Ux=g求x:xn=gn/βn,xi=(gi-cixi+1)/βi,i=n-1,n-2,…,1编程时，A可用三个一维数组，f用一个一维数组.L,U存入A。g，x存入f。还有一种计算格式,消去时用主元素除主行元素,即分解A为下三角矩阵和单位上三角矩阵之积,相当于对AT作LU分解.nnnnnnnngggaafffbaccbacb2112221121122211)()()(括号中是单位上三角矩阵的上对角元素.计算步骤:1．分解A=LU:β1=b1,γ1=c1/β1,βi=bi-aiγi-1,γi=ci/βi,i=2,3,…,n2．解Lg=f求g:g1=f1/β1,gi=(fi-aigi-1)/βi,i=2,3,…,n3．解Ux=g求x:xn=gn,xi=gi-γixi+1,i=n-1,n-2,…,1三对角矩阵是带形矩阵的特例.所谓带形矩阵是那些主对角线附近几条对角线以外元素皆零的矩阵,即aij≠0,仅当-m1j-im2.带形矩阵的LU分解也保持结构.5.4向量和矩阵的范数引入实数的绝对值和复数的模(也称绝对值)来表示实数和复数的”大小”,从而带来许多用处.例如,数列收敛的概念就是通过绝对值来表示的.范数这个概念就是这些表示”大小”的数值普遍化.它在研究数值计算方法的收敛性和稳定性中有着重要的应用.①向量的范数定义1.如果向量)(nnCRx或的某个实值函数xxN)(，满足条件：1.正定性:║x║≥0,║x║=0iffx=02.齐次性:║cx║=│c│║x║,Cc3.三角不等式:①║x+y║≤║x║+║y║②|yx|yx则称Cn中定义了向量范数║x║为向量x的范数。可见向量范数是向量的一种具有特殊性质的实值函数。常用向量范数有：（令x=(x1,x2,…,xn)T）1-范数：║x║1=│x1│+│x2│+…+│xn│2-范数：║x║2=(│x1│2+│x2│2+…+│xn│2)1/2∞-范数：║x║=max(│x1│,│x2│,…,│xn│)易得║x║≤║x║2≤║x║1≤n1/2║x║2≤n║x║P-范数：).,1[,)(11PxxniPPiP其中定理1.Cn中任意两种向量范数║x║,║x║是等价的,即m,M0使m║x║≤║x║≤M║x║可根据范数的连续性来证明它.由定理1可得。定理2.0limlim)()(xxxxkkkk，其中为向量的任一种范数。此时称{x(k)}收敛于x,记作x(k)→x(k→∞),或xxkk)(lim。②矩阵的范数定义2.设RXCXnn:,满足1.正定性:║X║≥0,║X║=0iffX=02.齐次性:║cX║=│c│║X║,Cc3.三角不等式:║X+Y║≤║X║+║Y║4.相容性:║XY║≤║X║║Y║则称Cnn中定义了矩阵范数,║X║为矩阵X的范数.注意:矩阵X可视为n2维向量,故有前三条性质.因此定理1,2中向量的等价性和向量序列收敛的概念与性质等也适合于矩阵.第四条,是考虑到矩阵乘法关系而设.║Ax║≤║A║║x║所谓由向量范数导出的矩阵范数与该向量范数就是相容的.定理3.设A是n×n矩阵,║║是n维向量范数则║A║=max{║Ax║/║x║=1}=max{║Ax║/║x║，x≠0}是一种矩阵范数,称为由该向量范数导出的矩阵范数或算子范数,它们具有相容性或者说是相容的。单位矩阵的算子范数为1。可以证明任一种矩阵范数总有与之相容的向量范数.例如定义:║x║=║X║,X=(xx…x)常用的三种向量范数导出的矩阵范数是1-范数:║A║1=max{║Ax║1/║x║1=1}=nijjnja11max2-范数:║A║2=max{║Ax║2/║x║2=1}=1,λ1是ATA的最大特征值.∞-范数:║A║=max{║Ax║/║x║=1}=njijnia11max此外还有Frobenius范数:njiijFaA1,212)(.它与向量2-范数相容.③矩阵譜半径定义3.设A是n×n矩阵,λi是其特征值,i=1,2,…,n.称iniA1max)(为A的譜半径.譜半径是矩阵的函数,但非矩阵范数.对任一矩阵范数有如下关系:ρ(A)≤║A║因为任一特征对λ,x,Ax=λx,令X=(xx…x),可得AX=λX.两边取范数,由矩阵范数的相容性和齐次性就导出结果.定理3.矩阵序列I,A,A2,…Ak,…收敛于零的充分必要条件是ρ(A)1.5.5误差分析①病态现象例3给出一个方程组顺序消去法解的误差很大,主元素法解的误差很小.该方程组数据有微小变化时解的变化也小.但有些方程组不是这样的,数据有微小变化时解的变化大.换句话说后一种方程组对数据变化敏感,前一种方程组对数据变化不敏感,这两种方程组(和相应的矩阵)分别称为病态的和良态的.例5.病态方程组110001.220001.1111,02220001.1111例6.病态矩阵28004200168014042006480270024016802700120012014024012016,7/16/15/14/16/15/14/13/15/14/13/12/14/13/12/11144HHH4取五位有效数字,其逆误差在前面第二、三位上:2871.14310.0-1726.1144.20-4310.0-6650.12771.3-246.491726.12771.3-1229.972.122144.20-246.4972.122248.16②扰动分析与矩阵条件数现在考虑系数、右端项有扰动时解的变化,也就是数据有误差时解的误差.设Ax=b,右端项有扰动A(x+x)=b+b,A可逆解皆存在惟一,其差x=A-1b,║x║≤║A-1║║b║,║x║/║x║≤(║A-1║║A║)║b║/║b║再考虑系数有扰动(A+A)(x+x)=b.首先,当A可逆,║A-1║║A║1时A+A可逆.因为此时ρ(A-1A)≤║A-1A║≤║A-1║║A║1,I+A-1A可逆,从而A+A=A(I+A-1A)可逆.原方程与扰动方程解皆存在惟一,二方程相减有Ax=-A(x+x),x=-A-1A(x+x)两边取范数可得║x║≤║A-1║║A║(║x║+║x║)从而有AAAAAAxx111类似的方法不难导出一般情况下,即系数、右端项都有扰动时的估计:AAxxAAAAAAxx111注意到估计式表明:║A-1║║A║不大,对解的影响也不大;║A-1║║A║越大,扰动对解的影响也越大.这就是说该向量是方程组敏感性以及病态或良态的度量，称为矩阵的条件数,记为Cond(A)=║A-1║║A║.它有如下性质:1.Cond(A)≥12.Cond(cA)=Cond(A),c≠03.Cond(A)2=║A-1║2║A║2=21称为谱条件数。λ1,λn分别是AHA的最大和最小特征值.故正交矩阵,酉矩阵的谱条件数为1.在例1中有Cond(A)1=2.00012×104.例2中Cond(H4)1=28000.另外,计算机计算解可归结为数据有一定扰动的准确解,因而可据以事先估计计算解的误差(向后误差分析).③事后误差分析计算解的误差还可根据下列定理用计算解的剩余量估计.定理4.设x和x*分别为非奇异方程组Ax=b(≠0)的准确解和近似解,r为x*的剩余量r=b-Ax*则brAxxx)(Cond因为║b║=║Ax║≤║A║║x║,║x*-x║=║A-1r║≤║A-1║║r║.由此可见对病态方程组剩余量小时误差还可能很大.例7.解方程组0.780x1+0.563x2=0.2170.913x1+0.659x2=0.254解x=(1,-1)T，x*=(0.341,-0.087)T，r=(-0.000001,0)T，x*-x=(-0.659,0.913)T二实验部分本章实验内容：实验题目：Gauss消元法,追赶法,范数。实验内容：①编制用Gauss消元法求解线性方程组Ax=f的程序。②编制用追赶法求解线性方程组Ax=f的程序。③编制向量和矩阵的范数程序。实验目的：①了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程组的算法,并能计算行列式的值。②掌握追赶法，能利用追赶法求解线性方程组。③理