解线性代数方程组.pdf

《解线性代数方程组.pdf》由会员分享，可在线阅读，更多相关《解线性代数方程组.pdf（11页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、求解线性方程组的直接解法求解线性方程组的直接解法5.35.3特殊矩阵的三角分解特殊矩阵的三角分解实对称矩阵的实对称矩阵的 LDLLDLT T分解分解设 A A 是实对称阵,且 A A 的所有顺序主子式均不为零，则 LDR 分解中R R=L LT,故可用以作 LDLT分解.这就是说,当 A A 的对角元素非零时，我们可以作 LU分解,也就得到LDLT分解,L L相同,是单位上三角阵,U U的对角元素构成 D D.不过没有利用对称性,存储量运算量都未能节省预计是一半。试用 n=3 的计算表格说明如何实现节省。d1=u11u12=a12u13=a13=a11l21=u12/d1l31=u13/d1d

2、2=u22=a22-l21u12u23=a23-l21u13l32=u23/d2u33=a33-l31u13-l32u23这样,可用上半部元素逐列计算 D,LD,LT。也可用下半部元素逐行计算 L,DL,D。引进輔助量 t1,t2代替 u1j,u2j,并利用对称性得到：d1=a11t1=a21l21=t1/d1t1=a31l31=t1/d1d2=a22-t1l21t2=a32-t1l21l32=t2/d2d3=a33-t1l31-t2l32据此不难写出 LDLT分解 A A=LDLLDLT的计算公式和程序(逐行计算 L,DL,D).d1=a11forfori=2:nforfor j=1:i-1

3、tj=aij-lj1t1-lj2t2-lj,j-1tj-1lij=tj/djendenddi=aii-li1t1-li2t2-li,i-1ti-1endend存储约 n(n+1)/2 单元,乘加运算各约 n3/6.利用 LDLT分解解 AxAx=b b 分四步:1分解 A A=LDLLDLT2解LgLg=b b求 g g3解DyDy=g g求 y y4解L LTx x=y y求 x x实对称正定矩阵的实对称正定矩阵的 LLLLT T分解分解A A 实对称正定时顺序主子式皆正,可作 LDLT,D 的对角元素皆正,有正1/11的平方根。因此有 LLT分解 A A=LLLLT,L L 下三角阵,对角

4、元素皆正,是 LDLT中的 LDLD1/2.乃可用上半部元素逐列计算 L LT.l11=a111/2l21=a12/l11l31=a13/l11l11=a111/2l21=a21/l11l31=a31/l11l11=a111/2forfori=2:nforfor j=1:i-1lij=(aij-li1lj1-li2lj2-li,j-1lj,j-1)/djjendend221/2lii(aiili21li2li,i1)l22=(a22-l212)1/2l32=(a23-l21l31)/l22l33=a33-l312-l322也可用下半部元素逐行计算 L L.计算表格和算法安排如下:l33=(a3

5、3-l312-l322)1/2l22=(a22-l212)1/2l32=(a32-l31l21)/l22endend存储量,运算量同 LDLT分解,但要 n 次求平方根.利用 LLT分解解 AxAx=b b 分三步:1分解 A A=LLLLT2解LgLg=b b求 g g3解L LTx=g g求 x x三对角方程组的追赶法三对角方程组的追赶法消去法或 LU 分解用于三对角方程组有特殊形式,即称追赶法追赶法.设 AxAx=f f:b1x1+c1x2=f1aixi-1+bixi+cixi+1=fii=2,3,n-1anxn-1+bnxn=fnA A 是三对角阵,则 L L,U U 同样结构.L L

6、 的对角元素为2,3,，n,U U 的对角元素为1,2,n,上对角元素同 A A.1分解 A A=LULU:1=b1,i=ai/i-1,i=bi-i ci-1,i=2,3,n2解LgLg=f f 求 g g:g1=f1,gi=fi-ifi-1,i=2,3,n3解U Ux=g g 求 x x:xn=gn/n,x i=(gi-cixi+1)/i,i=n-1,n-2,1编程时，A A 可用三个一维数组，f f 用一个一维数组.L L,U U 存入 A A。g g，x x存入 f f。还有一种计算格式,消去时用主元素除主行元素,即分解 A A 为下三角矩阵和单位上三角矩阵之积,相当于对 A AT作 L

7、U 分解.2/11b1c1a2b2c2 cn1anbnf11(1)f2a22(2)(n1)fnanng1g2gn括号中是单位上三角矩阵的上对角元素.计算步骤:1分解 A A=LULU:1=b1,1=c1/1,i=bi-aii-1,i=ci/i,i=2,3,n2解LgLg=f f求 g g:g1=f1/1,gi=(fi-aigi-1)/i,i=2,3,n3解U Ux=g g求 x x:xn=gn,x i=gi-ixi+1,i=n-1,n-2,1三对角矩阵是带形矩阵的特例.所谓带形矩阵是那些主对角线附近几条对角线以外元素皆零的矩阵,即 aij0,仅当-m1j-i0 使mx xx xMx x3/11

8、可根据范数的连续性来证明它.由定理 1 可得。定理定理 2 2.lim x(k)x lim x(k)x 0，其中为向量的任一种范数。此时kk称x x(k)收敛收敛于 x x,记作 x x(k)x x(k),或lim x(k)x。k矩阵的范数矩阵的范数定义定义 2 2.设:X Cnn X R,满足1.1.正定性:X X0,X X=0 iff X X=0 02.2.齐次性:cX X=cX X,cC3.3.三角不等式:X X+Y YX X+Y Y4.4.相容性:XYXYX XY Y则称 Cnn中定义了矩阵范数,X X为矩阵 X X 的范数.注意注意:矩阵 X X 可视为 n2维向量,故有前三条性质.

9、因此定理 1,2 中向量的等价性和向量序列收敛的概念与性质等也适合于矩阵.第四条,是考虑到矩阵乘法关系而设.AxAxA Ax x所谓由向量范数导出的矩阵范数与该向量范数就是相容的.定理定理 3 3.设 A A 是 nn 矩阵,是 n 维向量范数则A A=maxAxAx/x x=1=maxAxAx/x x，x x0 0是一种矩阵范数,称为由该向量范数导出的矩阵范数或算子范数,它们具有相容性或者说是相容的。单位矩阵的算子范数为单位矩阵的算子范数为 1 1。可以证明任一种矩阵范数总有与之相容的向量范数.例如定义:x x=X X,X X=(xxxxx x)常用的三种向量范数导出的矩阵范数是1-1-范数

10、范数:A A1=maxAxAx1/x x1=1=maxajj1 jni1n2-2-范数范数:A A2=maxAxAx2/x x2=1=1,1是A ATA A的最大特征值.-范数范数:A A=maxAxAx/x x=1=maxaij1inj1n此外还有 FrobeniusFrobenius 范数范数:AF(aij)i,j1n212.它与向量 2-范数相容.矩阵譜半径矩阵譜半径定义定义 3 3.设 A A 是 nn 矩阵,i是其特征值,i=1,2,n.称4/11i(A)m ax1in为 A A 的譜半径.譜半径是矩阵的函数,但非矩阵范数.对任一矩阵范数有如下关系:(A A)A A因为任一特征对,x

11、 x,AxAx=x x,令 X X=(xxxxx x),可得 AXAX=X X.两边取范数,由矩阵范数的相容性和齐次性就导出结果.2k k定理定理 3.3.矩阵序列 I I,A A,A A,A A,收敛于零的充分必要条件是(A A)1.5.55.5误差分析误差分析病态现象病态现象例 3 给出一个方程组顺序消去法解的误差很大,主元素法解的误差很小.该方程组数据有微小变化时解的变化也小.但有些方程组不是这样的,数据有微小变化时解的变化大.换句话说后一种方程组对数据变化敏感,前一种方程组对数据变化不敏感,这两种方程组(和相应的矩阵)分别称为病态的和良态的.例例 5.5.病态方程组11221121,1

12、 1.0001201 1.00012.00011例例 6.6.病态矩阵1/4120240140 161201200270016801/51,H4 2402700648042001/61/7140168042002800H4取五位有效数字,其逆误差在前面第二、三位上:16.248122.72246.49122.721229.9-2771.3 246.49-2771.36650.1-144.201726.1-4310.0-144.201726.1-4310.02871.1 11/2H41/31/41/21/31/41/51/31/41/51/6扰动分析与矩阵条件数扰动分析与矩阵条件数现在考虑系数、

13、右端项有扰动时解的变化,也就是数据有误差时解的误差.设 AxAx=b,右端项有扰动 A A(x x+x x)=b b+b b,A A 可逆解皆存在惟一,其差x x=A A-1b,x xA A-1b b,-1x x/x x(A AA A)b b/b b再考虑系数有扰动(A A+A A)(x x+x x)=b.b.首先,当 A A 可逆,A A-1A A1 时A A+A A 可逆.因为此时(A A-1A A)A A-1A AA A-1A A1,I I+A A-1A A 可逆,从而 A A+A=AA=A(I I+A A-1A A)可逆.原方程与扰动方程解皆存在惟一,二方程相减有A Ax x=-A A

14、(x x+x x),x x=-A A-1A A(x x+x x)两边取范数可得x xA A-1A A(x x+x x)从而有5/11xxA A11 A1AAA类似的方法不难导出一般情况下,即系数、右端项都有扰动时的估计:xxA A11 A1AA xA AxA注意到估计式表明:A A-1A A不大,对解的影响也不大;A A-1A A越大,扰动对解的影响也越大.这就是说该向量是方程组敏感性以及病态或良态的度量，称为矩阵的条件数条件数,记为 Cond(A A)=A A-1A A.它有如下性质:1.1.Cond(A A)12.2.Cond(cA A)=Cond(A A),c03.3.Cond(A A)

15、2=A A-12A A2=12称为谱条件数。1,n分别是 A AHA A 的最大和最小特征值.故正交矩阵,酉矩阵的谱条件数为 1.在例 1 中有 Cond(A A)1=2.00012104.例2中 Cond(H H4)1=28000.另外,计算机计算解可归结为数据有一定扰动的准确解,因而可据以事先估计计算解的误差(向后误差分析).事后误差分析事后误差分析计算解的误差还可根据下列定理用计算解的剩余量估计.定理定理 4 4.设 x x 和 x x*分别为非奇异方程组 AxAx=b b(0 0)的准确解和近似解,r r 为 x x*的剩余量 r r=b b-AxAx*则x x因为b b=AxAxA

16、Ax x,x x*-x x=A A-1r rA A-1r r.x Cond(A)rb由此可见对病态方程组剩余量小时误差还可能很大.例例 7.7.解方程组0.780 x x1+0.563x x2=0.2170.913x x1+0.659 x x2=0.254解解x x=(1,-1)T，x x*=(0.341,-0.087)T，r r=(-0.000 001,0)T，x x*-x x=(-0.659,0.913)T二二实验部分实验部分本章实验内容：本章实验内容：实验题目：实验题目：GaussGauss 消元法消元法,追赶法追赶法,范数。范数。实验内容：实验内容：编制用编制用 GaussGauss

17、消元法求解线性方程组消元法求解线性方程组 Ax=fAx=f 的程序。的程序。编制用追赶法求解线性方程组编制用追赶法求解线性方程组Ax=fAx=f 的程序。的程序。编制向量和矩阵的范数程序。编制向量和矩阵的范数程序。实验目的：实验目的：了解了解 GaussGauss 消元法原理及实现条件消元法原理及实现条件,熟练掌握熟练掌握 GaussGauss6/11消元法解方程组的算法消元法解方程组的算法,并能计算行列式的值。并能计算行列式的值。掌握追赶法，能利用追赶法求解线性方程组。掌握追赶法，能利用追赶法求解线性方程组。理解向量和矩阵范数定义理解向量和矩阵范数定义,性质并掌握其计算方性质并掌握其计算方法

18、法.编程要求：编程要求：利用利用 GaussGauss 消元法消元法,追赶法解线性方程组。追赶法解线性方程组。分析误差。分析误差。计算算法计算算法:GaussGauss 消元法消元法:1.消元过程消元过程(k)设设akk 0,对对k 1,2,n 1计算计算(k)(k)mik aik/akk(k1)(k)(k)aij aijmikakj(k1)bi(k)mikbk(k)bii,j k 1,k 2,n回代过程回代过程(n)(n)xn bn/annn(i)(i)(i)x (bax)/aiiijjiiji1i n 1,2,1追赶法：追赶法：1.1.分解分解 Ax=f:Ax=f:1 c1/b1,i ci

19、/(biaii1)(i 2,3,n 1)2.2.解解 Lg=f,Lg=f,求求 g:g:g1 f1/b1,gi(fiaigi1)/(biaii1)(i 2,3,n)3 3解解 Ux=g,Ux=g,求求 x:x:xn yn,xi yiixi1（i n 1,n 2,2,1）范数：范数：常用向量范数有：常用向量范数有：（令（令x x=(x1,x2,xn)T）1-1-范数范数：x x1=x1+x2+xn2-2-范数范数：x x2=(x1+x2+xn)2221/2-范数范数：x x=max(x1,x2,xn)常用的三种向量范数导出的矩阵范数是：常用的三种向量范数导出的矩阵范数是：7/111-1-范数范数

20、:A A1=maxAxAx1/x x1=1=maxajj1 jni1n2-2-范数范数:A A2=maxAxAx2/x x2=1=1,1是是A ATA A的的最大特征值最大特征值.-范数范数:A A=maxAxAx/x x=1=maxaij1inj1n实验例题：用实验例题：用GaussGauss消元法解方程组消元法解方程组 211 x14 1 23x251 31x36实验例题实验例题:用追赶法解三对角方程组用追赶法解三对角方程组Ax=bAx=b,其中其中00 2101121000 A 01210,b 0 001210001200实验例题：设实验例题：设0.60.5A 0.10.3，计算计算 A

21、 A 的行列范数的行列范数.程序程序:Gauss:Gauss 消元法消元法function x=Gauss(A,b)function x=Gauss(A,b)%A%A是线性方程组的系数矩阵是线性方程组的系数矩阵,b b为自由项为自由项.n=length(A)n=length(A)for k=1:n-1for k=1:n-1m(k+1:n,k)=A(k+1:n,k)/A(k,k);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);A(k+1:n,k+1:n)=A(k+1:n,k+1:n

22、)-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);b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k);endendx=zeros(n,1);x=zeros(n,1);x(n)=b(n)/A(n,n);x(n)=b(n)/A(n,n);for k=n-1:-1:1for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endendx=x;x=x;8/11disp(sprintf(kdisp(s

23、printf(kx(k);x(k);for i=0:nfor i=0:ndisp(sprintf(%ddisp(sprintf(%d%f,i,x(i+1);%f,i,x(i+1);endend数值结果：数值结果：x=Gauss(A,b)x=Gauss(A,b)n=3n=3k kx(k)x(k)0 01.1111111.1111111 10.7777780.7777782 22.5555562.555556程序程序:追赶法追赶法function x,y,beta=zhuiganfa(a,b,c,f)function x,y,beta=zhuiganfa(a,b,c,f)%a,b,c%a,b,c是

24、三对角阵的对角线上的元素是三对角阵的对角线上的元素,f f是自由项是自由项.n=length(b);n=length(b);beta(1)=c(1)/b(1);beta(1)=c(1)/b(1);for i=2:nfor i=2:nbeta(i)=c(i)/(b(i)-a(i)*beta(i-1);beta(i)=c(i)/(b(i)-a(i)*beta(i-1);endendy(1)=f(1)/b(1);y(1)=f(1)/b(1);for i=2:nfor i=2:ny(i)=(f(i)-a(i)*y(i-1)/(b(i)-a(i)*beta(i-1);y(i)=(f(i)-a(i)*y(

25、i-1)/(b(i)-a(i)*beta(i-1);endendx(n)=y(n);x(n)=y(n);for i=n-1:-1:1for i=n-1:-1:1x(i)=y(i)-beta(i)*x(i+1);x(i)=y(i)-beta(i)*x(i+1);endenddisp(sprintf(kdisp(sprintf(kx(k)x(k)y(k)y(k)beta(k);beta(k);for i=0:nfor i=0:ndisp(sprintf(%ddisp(sprintf(%d%f,i,x(i+1),y(i+1),beta(i+1);%f,i,x(i+1),y(i+1),beta(i+1

26、);endend数值结果数值结果:a=0-1-1-1-1;a=0-1-1-1-1;b=2 2 2 2 2;b=2 2 2 2 2;c=-1-1-1-1 0;c=-1-1-1-1 0;f=1 0 0 0 0;f=1 0 0 0 0;x,y,beta=zhuiganfa(a,b,c,f)x,y,beta=zhuiganfa(a,b,c,f)k kx(k)x(k)y(k)y(k)beta(k)beta(k)0 00.833333 5.000000e-0010.833333 5.000000e-001-0.500000-0.5000001 10.666667 3.333333e-0010.666667

27、 3.333333e-001-0.666667-0.6666672 20.500000 2.500000e-0010.500000 2.500000e-001-0.750000-0.7500003 30.333333 2.000000e-0010.333333 2.000000e-001-0.800000-0.8000004 40.166667 1.666667e-0010.166667 1.666667e-0010.0000000.000000程序：程序：1.1.列范数列范数:functionfunctionfan=lie(A)fan=lie(A)%A%A为已知矩阵为已知矩阵9/11n=le

28、ngth(A)n=length(A)for j=1:nfor j=1:nx(j)=0 x(j)=0for i=1:nfor i=1:nx(j)=x(j)+abs(A(i,j);x(j)=x(j)+abs(A(i,j);endendendendfan=max(x)fan=max(x)disp(sprintf(ndisp(sprintf(nx(n);x(n);for i=0:nfor i=0:ndisp(sprintf(%ddisp(sprintf(%d%f,i,x(i+1);%f,i,x(i+1);endend数值结果数值结果:fan=lie(A)fan=lie(A)fan=0.8000fan=

29、0.8000n nx(n)x(n)0 00.7000000.7000001 10.8000000.8000002.2.行范数行范数:function fan=hang(A)function fan=hang(A)%A%A为已知矩阵为已知矩阵n=length(A)n=length(A)for i=1:nfor i=1:nx(i)=0 x(i)=0for j=1:nfor j=1:nx(i)=x(i)+abs(A(i,j);x(i)=x(i)+abs(A(i,j);endendendendfan=max(x)fan=max(x)disp(sprintf(ndisp(sprintf(nx(n);x(

30、n);for i=0:nfor i=0:ndisp(sprintf(%ddisp(sprintf(%d%f,i,x(i+1);%f,i,x(i+1);endend数值结果数值结果:fan=hang(A)fan=hang(A)fan=1.1000fan=1.1000n nx(n)x(n)0 01.1000001.1000001 10.4000000.400000总结总结:从代数上看从代数上看,直接分解法和直接分解法和GaussGauss消去法本质上一样消去法本质上一样,但如果我们但如果我们采用采用”双精度累加双精度累加”计算计算aibi,那么直接三角分解法的精度要比那么直接三角分解法的精度要比G

31、aussGauss消去法为高消去法为高.求线性方程组的直接法求线性方程组的直接法,其算式繁杂其算式繁杂,给人以枯燥沉闷的感觉给人以枯燥沉闷的感觉.为了改善教学效果为了改善教学效果,本章着重介绍了三对角方程组的追赶法本章着重介绍了三对角方程组的追赶法.三对三对10/11角方程组以及其拓广形式的带状方程组有着广泛的实际应用角方程组以及其拓广形式的带状方程组有着广泛的实际应用.追赶法是解三对角线方程组追赶法是解三对角线方程组(对角元占优势对角元占优势)的有效方法的有效方法,它它具有计算量少具有计算量少,方法简单方法简单,算法稳定等优点算法稳定等优点,具有鲜明的对称美具有鲜明的对称美.复习思考题复习思考题 1.用消去法解线性方程组为什么最好选主元？怎样的方程组可以不用选主元？2.用高斯约当消去法求矩阵的逆，其理论根据是什么？3.求矩阵A的LU分解的紧凑格式有何规律？写出LU分解法解Ax=b的算法。4.乔累斯基分解法适用于哪类线性方程组？有何优点？5.追赶法有何优点？哪类方程组可用追赶法求解？6.向量范数与矩阵范数是如何定义的？各有何性质？常用范数有哪些？7.什么叫病态方程组？什么叫矩阵的条件数？条件数刻画了矩阵的什么性质？能否用选主元的方法求得病态方程组的高精度的解？11/11