第五章 解线性代数方程组.doc

求解线性方程组的直接解法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)/djjendend存储量,运算量同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-i ci-1,i=2,3,n2解Lg=f求g:g1=f1,gi=fi-ifi-1, i=2,3,n3解Ux=g求x:xn=gn/n,x i=(gi-cixi+1)/i,i=n-1,n-2,1 编程时，A可用三个一维数组，f用一个一维数组.L,U存入A。g，x存入f。还有一种计算格式,消去时用主元素除主行元素,即分解A为下三角矩阵和单位上三角矩阵之积,相当于对AT作LU分解.括号中是单位上三角矩阵的上对角元素.计算步骤:1分解A=LU:1=b1,1=c1/1,i=bi-aii-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,x i=gi-ixi+1,i=n-1,n-2,1三对角矩阵是带形矩阵的特例.所谓带形矩阵是那些主对角线附近几条对角线以外元素皆零的矩阵,即aij0,仅当-m1j-i0使mxxMx可根据范数的连续性来证明它.由定理1可得。定理2.，其中为向量的任一种范数。此时称x(k)收敛于x,记作x(k) x(k),或。 矩阵的范数定义2.设,满足1. 正定性:X0,X=0 iff X=02. 齐次性:cX=cX, 3. 三角不等式:X+YX+Y4. 相容性: XYXY则称Cnn中定义了矩阵范数,X为矩阵X的范数.注意:矩阵X可视为n2维向量,故有前三条性质.因此定理1,2中向量的等价性和向量序列收敛的概念与性质等也适合于矩阵.第四条,是考虑到矩阵乘法关系而设.AxAx所谓由向量范数导出的矩阵范数与该向量范数就是相容的.定理3. 设A是nn矩阵,是n维向量范数则A=maxAx/x=1= maxAx/x，x0是一种矩阵范数,称为由该向量范数导出的矩阵范数或算子范数,它们具有相容性或者说是相容的。单位矩阵的算子范数为1。可以证明任一种矩阵范数总有与之相容的向量范数.例如定义:x=X,X=(xxx)常用的三种向量范数导出的矩阵范数是1-范数:A1= maxAx1/x1=1= 2-范数:A2=maxAx2/x2=1=,1是ATA的最大特征值. -范数:A=maxAx/x=1=此外还有Frobenius范数:.它与向量2-范数相容. 矩阵譜半径定义3.设A是nn矩阵,i是其特征值,i=1,2,n.称 为A的譜半径.譜半径是矩阵的函数,但非矩阵范数.对任一矩阵范数有如下关系:(A)A因为任一特征对,x,Ax=x,令X=(xxx),可得AX=X.两边取范数,由矩阵范数的相容性和齐次性就导出结果.定理3. 矩阵序列I,A,A2,Ak,收敛于零的充分必要条件是(A)1.5.5误差分析 病态现象例3给出一个方程组顺序消去法解的误差很大,主元素法解的误差很小.该方程组数据有微小变化时解的变化也小.但有些方程组不是这样的,数据有微小变化时解的变化大.换句话说后一种方程组对数据变化敏感,前一种方程组对数据变化不敏感,这两种方程组(和相应的矩阵)分别称为病态的和良态的.例5. 病态方程组例6. 病态矩阵H4取五位有效数字,其逆误差在前面第二、三位上: 扰动分析与矩阵条件数现在考虑系数、右端项有扰动时解的变化,也就是数据有误差时解的误差.设Ax=b,右端项有扰动A(x+x)=b+b,A可逆解皆存在惟一,其差x=A-1b,xA-1b,x/x(A-1A)b/b再考虑系数有扰动(A+A)(x+x)=b.首先,当A可逆,A-1A1时A+A可逆.因为此时(A-1A)A-1AA-1A1,I+A-1A可逆,从而A+A=A(I+A-1A)可逆.原方程与扰动方程解皆存在惟一,二方程相减有Ax= -A(x+x), x= -A-1A(x+x)两边取范数可得xA-1A(x+x)从而有 类似的方法不难导出一般情况下,即系数、右端项都有扰动时的估计:注意到估计式表明:A-1A不大,对解的影响也不大; A-1A越大,扰动对解的影响也越大.这就是说该向量是方程组敏感性以及病态或良态的度量，称为矩阵的条件数,记为Cond(A)=A-1A.它有如下性质:1. Cond(A)12. Cond(cA)=Cond(A),c03. Cond(A)2=A-12A2=称为谱条件数。1,n分别是AHA的最大和最小特征值.故正交矩阵,酉矩阵的谱条件数为1.在例1中有Cond (A)1=2.00012104.例2中Cond (H4)1=28000.另外,计算机计算解可归结为数据有一定扰动的准确解,因而可据以事先估计计算解的误差(向后误差分析). 事后误差分析计算解的误差还可根据下列定理用计算解的剩余量估计.定理4.设x和x*分别为非奇异方程组Ax=b(0)的准确解和近似解,r为x*的剩余量r = b-Ax*则因为b=AxAx,x*-x=A-1rA-1r.由此可见对病态方程组剩余量小时误差还可能很大.例7.解方程组0.780x1+0.563x2=0.2170. 913x1+0.659 x2=0.254解 x=(1,-1)T，x*=(0.341,-0.087)T，r=(-0.000 001,0)T，x*-x =(-0.659,0.913)T 二 实验部分本章实验内容：实验题目：Gauss消元法,追赶法,范数。实验内容：编制用Gauss消元法求解线性方程组Ax=f的程序。 编制用追赶法求解线性方程组Ax=f的程序。 编制向量和矩阵的范数程序。实验目的：了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程组的算法,并能计算行列式的值。 掌握追赶法，能利用追赶法求解线性方程组。理解向量和矩阵范数定义,性质并掌握其计算方法.编程要求：利用Gauss消元法,追赶法解线性方程组。分析误差。计算算法:Gauss消元法:1. 消元过程设,对计算回代过程追赶法：1.分解Ax=f: ( )2.解Lg=f,求g: ()3解Ux=g,求x: （）范数：常用向量范数有：（令x=( x1,x2,xn)）1-范数： x1=x1+x2+xn2-范数： x2=(x12+x22+xn2)1/2-范数： x=max(x1,x2,xn)常用的三种向量范数导出的矩阵范数是：1-范数:A1= maxAx1/x1=1=2-范数:A2=maxAx2/x2=1=,1是ATA的最大特征值.-范数:A=maxAx/x=1=实验例题：用Gauss消元法解方程组 实验例题: 用追赶法解三对角方程组Ax=b,其中 实验例题：设 ， 计算A的行列范数.程序:Gauss消元法function x=Gauss(A,b) %A是线性方程组的系数矩阵,b为自由项.n=length(A)for k=1:n-1 m(k+1:n,k)=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m(k+1:n,k)*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k);endx=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endx=x;disp(sprintf(k x(k);for i=0:n disp(sprintf(%d %f ,i,x(i+1);end数值结果：x=Gauss(A,b)n =3k x(k)0 1.111111 1 0.777778 2 2.555556程序:追赶法function x,y,beta=zhuiganfa(a,b,c,f)%a,b,c是三对角阵的对角线上的元素,f是自由项.n=length(b);beta(1)=c(1)/b(1);for i=2:n beta(i)=c(i)/(b(i)-a(i)*beta(i-1);endy(1)=f(1)/b(1);for i=2:n y(i)=(f(i)-a(i)*y(i-1)/(b(i)-a(i)*beta(i-1);endx(n)=y(n);for i=n-1:-1:1 x(i)=y(i)-beta(i)*x(i+1);enddisp(sprintf(k x(k) y(k) beta(k);for i=0:n disp(sprintf(%d %f ,i,x(i+1),y(i+1),beta(i+1);end数值结果:a=0 -1 -1 -1 -1;b=2 2 2 2 2;c=-1 -1 -1 -1 0;f=1 0 0 0 0;x,y,beta=zhuiganfa(a,b,c,f)k x(k) y(k) beta(k)0 0.833333 5.000000e-001 -0.500000 1 0.666667 3.333333e-001 -0.666667 2 0.500000 2.500000e-001 -0.750000 3 0.333333 2.000000e-001 -0.800000 4 0.166667 1.666667e-001 0.000000程序： 1.列范数:function fan=lie(A)%A为已知矩阵n=length(A)for j=1:n x(j)=0 for i=1:n x(j)=x(j)+abs(A(i,j); endendfan=max(x)disp(sprintf(n x(n);for i=0:ndisp(sprintf( %d %f,i,x(i+1);end数值结果: fan=lie(A) fan =0.8000n x(n)0 0.7000001 0.8000002.行范数:function fan=hang(A)%A为已知矩阵n=length(A)for i=1:n x(i)=0 for j=1:n x(i)=x(i)+abs(A(i,j); endendfan=max(x)disp(sprintf(n x(n);for i=0:n disp(sprintf( %d %f,i,x(i+1);end数值结果: fan=hang(A) fan =1.1000 n x(n)0 1.1000001 0.400000总结:从代数上看,直接分解法和Gauss消去法本质上一样,但如果我们采用”双精度累加”计算,那么直接三角分解法的精度要比Gauss消去法为高. 求线性方程组的直接法,其算式繁杂,给人以枯燥沉闷的感觉.为了改善教学效果,本章着重介绍了三对角方程组的追赶法.三对角方程组以及其拓广形式的带状方程组有着广泛的实际应用.追赶法是解三对角线方程组(对角