matlab在科学计算中的应用4(姜志鹏).ppt

第四章线性代数问题求解,矩阵线性方程组的直接解法线性方程组的迭代法线性方程组的符号解法稀疏矩阵技术特征值与特征向量,1,4.1矩阵4.1.1特殊矩阵的输入,数值矩阵的输入零矩阵、幺矩阵及单位矩阵生成nn方阵：A=zeros(n),B=ones(n),C=eye(n)生成mn矩阵：A=zeros(m,n),B=ones(m,n),C=eye(m,n)生成和矩阵B同样位数的矩阵：A=zeros(size(B),2,随机元素矩阵若矩阵随机元素满足0,1区间上的均匀分布生成nm阶标准均匀分布伪随机数矩阵：A=rand(n,m)生成nn阶标准均匀分布伪随机数方阵：A=rand(n),3,对角元素矩阵已知向量生成对角矩阵：A=diag(V)已知矩阵提取对角元素列向量：Vdiag(A)生成主对角线上第k条对角线为V的矩阵：A=diag(V,k),4,例：diag()函数的不同调用格式C=123;V=diag(C)%生成对角矩阵V=100020003V1=diag(V)%将列向量通过转置变换成行向量V1=123C=123;V=diag(C,2)%主对角线上第k条对角线为C的矩阵V=0010000020000030000000000,5,生成三对角矩阵:V=diag(1234)+diag(234,1)+diag(543,-1)V=1200523004340034,6,Hilbert矩阵及逆Hilbert矩阵生成n阶的Hilbert矩阵：A=hilb(n)求取逆Hilbert矩阵：B=invhilb(n),7,Hankel(汉克)矩阵其中：第一列的各个元素定义为C向量，最后一行各个元素定义为R。H为对称阵。H1=hankel(C)由Hankel矩阵反对角线上元素相等得出一下三角阵均为零的Hankel矩阵,8,Vandermonde(范德蒙)矩阵,9,伴随矩阵其中：P(s)为首项系数为1的多项式。,10,例：考虑一个多项式2*x4+4*x2+5*x+6，试写出该多项式的伴随矩阵。P=20456;A=compan(P)A=0-2.0000-2.5000-3.00001.000000001.000000001.00000,11,符号矩阵的输入数值矩阵A转换成符号矩阵：B=sym(A)例：A=hilb(3)A=1.00000.50000.33330.50000.33330.25000.33330.25000.2000B=sym(A)B=1,1/2,1/31/2,1/3,1/41/3,1/4,1/5,12,4.1.2矩阵基本概念与性质,行列式格式：d=det(A)例：求行列式A=162313;511108;97612;414151;det(A)ans=0,13,例：tic,A=sym(hilb(20);det(A),tocans=1/2377454716768534509091644243427616440175419837753486493033185331234419759310644585187585766816573773440565759867265558971765638419710793303386582324149811241023554489166154717809635257797836800000000000000000000000000000000000elapsed_time=2.3140高阶的Hilbert矩阵是接近奇异的矩阵。,14,矩阵的迹格式：t=trace(A)矩阵的秩格式：r=rank(A)用默认的精度求数值秩r=rank(A，)给定精度下求数值秩矩阵的秩也表示该矩阵中行列式不等于0的子式的最大阶次。可证行秩和列秩（线性无关的）应相等。,15,例A=162313;511108;97612;414151;rank(A)ans=3该矩阵的秩为3，小于矩阵的阶次，故为非满秩矩阵。例用数值方法和解析方法分别求2020的Hilbert矩阵的秩H=hilb(20);rank(H)数值方法ans=13H=sym(hilb(20);rank(H)%解析方法，原矩阵为非奇异矩阵ans=20,16,矩阵范数,17,矩阵的范数定义：格式：N=norm(A)求解默认的2范数N=norm(A，选项)选项可为1,2,inf等,18,例：求一向量、矩阵的范数a=162313;norm(a),norm(a,2),norm(a,1),norm(a,Inf)ans=2.092844953645635e+0012.092844953645635e+0013.400000000000000e+0011.600000000000000e+001A=162313;511108;97612;414151;norm(A),norm(A,2),norm(A,1),norm(A,Inf)ans=34343434符号运算工具箱未提供norm()函数，需先用double()函数转换成双精度数值矩阵，再调用norm()函数。,19,特征多项式格式：C=poly(A)例：A=162313;511108;97612;414151;poly(A)直接求取ans=1.000000000000000e+000-3.399999999999999e+001-7.999999999999986e+0012.719999999999999e+003-2.819840539024018e-012A=sym(A);poly(A)运用符号工具箱ans=x4-34*x3-80*x2+2720*x,20,矩阵多项式的求解,21,符号多项式与数值多项式的转换格式：f=poly2sym(P)或f=poly2sym(P,x)格式：P=sym2poly(f),22,例：P=123456;%先由系数按降幂顺序排列表示多项式f=poly2sym(P,v)%以v为算子表示多项式f=v5+2*v4+3*v3+4*v2+5*v+6P=sym2poly(f)P=123456,23,矩阵的逆矩阵格式：C=inv(A)例：求Hilbert矩阵的逆矩阵formatlong;H=hilb(4);H1=inv(H)H1=1.0e+003*0.01600000000000-0.119999999999990.23999999999998-0.13999999999999-0.119999999999991.19999999999990-2.699999999999761.679999999999840.23999999999998-2.699999999999766.47999999999940-4.19999999999961-0.139999999999991.67999999999984-4.199999999999612.79999999999974,24,检验：H*H1ans=1.000000000000010.00000000000023-0.000000000000450.000000000000230.000000000000011.00000000000011-0.000000000000110.000000000000110.0000000000000101.0000000000001100.000000000000000.00000000000011-0.000000000000111.00000000000011计算误差范数：norm(H*inv(H)-eye(size(H)ans=6.235798190375727e-013H2=invhilb(4);norm(H*H2-eye(size(H)ans=5.684341886080802e-014,25,H=hilb(10);H1=inv(H);norm(H*H1-eye(size(H)ans=0.00264500826202H2=invhilb(10);norm(H*H2-eye(size(H)ans=1.612897415528547e-005H=hilb(13);H1=inv(H);norm(H*H1-eye(size(H)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=2.339949e-018.ans=53.23696008570294H2=invhilb(13);norm(H*H2-eye(size(H)ans=11.37062973181391对接近于奇异矩阵，高阶一般不建议用inv()，可用符号工具箱。,26,H=sym(hilb(7);inv(H)ans=49,-1176,8820,-29400,48510,-38808,12012-1176,37632,-317520,1128960,-1940400,1596672,-5045048820,-317520,2857680,-10584000,18711000,-15717240,5045040-29400,1128960,-10584000,40320000,-72765000,62092800,-2018016048510,-1940400,18711000,-72765000,133402500,-115259760,37837800-38808,1596672,-15717240,62092800,-115259760,100590336,-3329726412012,-504504,5045040,-20180160,37837800,-33297264,11099088H=sym(hilb(30);norm(double(H*inv(H)-eye(size(H)ans=0,27,例：奇异阵求逆A=162313;511108;97612;414151;formatlong;B=inv(A)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.306145e-017.B=1.0e+014*0.938249922368852.81474976710656-2.81474976710656-0.938249922368852.814749767106568.44424930131968-8.44424930131968-2.81474976710656-2.81474976710656-8.444249301319688.444249301319682.81474976710656-0.93824992236885-2.814749767106562.814749767106560.93824992236885norm(A*B-eye(size(A)检验ans=1.64081513306419A=sym(A);inv(A)奇异矩阵不存在一个相应的逆矩阵，用符号工具箱的函数也不行?Errorusing=sym/invError,(ininverse)singularmatrix,28,同样适用于含有变量的矩阵求逆。例：symsa1a2a3a4;C=a1a2;a3a4；inv(C)ans=-a4/(-a1*a4+a2*a3),a2/(-a1*a4+a2*a3)a3/(-a1*a4+a2*a3),-a1/(-a1*a4+a2*a3),29,矩阵的相似变换与正交矩阵其中：A为一方阵，B矩阵非奇异。相似变换后，X矩阵的秩、迹、行列式与特征值等均不发生变化，其值与A矩阵完全一致。对于一类特殊的相似变换满足如下条件，称为正交基矩阵。,30,例：A=5,9,8,3;0,3,2,4;2,3,5,9;3,4,5,8;Q=orth(A)Q=-0.61970.7738-0.0262-0.1286-0.2548-0.15510.94900.1017-0.5198-0.5298-0.1563-0.6517-0.5300-0.3106-0.27250.7406norm(Q*Q-eye(4)ans=4.6395e-016norm(Q*Q-eye(4)ans=4.9270e-016,31,C=Q*A*QC=17.92516.4627-4.4714-2.0354-0.02821.71944.6816-5.07350.6800-0.93861.06740.6631-0.05490.36580.17760.2882det(A),det(C)ans=120ans=120.0000,32,trace(A),trace(C)ans=21ans=21.0000rank(A),rank(C)ans=4ans=4,33,eig(A),eig(C)ans=17.82051.1908+2.6499i1.1908-2.6499i0.7979ans=17.82051.1908+2.6499i1.1908-2.6499i0.7979,34,例：A=16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1;Q=orth(A)A为奇异矩阵，故得出的Q为长方形矩阵Q=-0.50000.67080.5000-0.5000-0.2236-0.5000-0.50000.2236-0.5000-0.5000-0.67080.5000norm(Q*Q-eye(3)%Q*Q=Ians=1.0140e-015,35,4.2线性方程组直接解法4.2.1线性方程组直接求解矩阵除法,关于线性方程组的直接解法，如Gauss消去法、选主元消去法、平方根法、追赶法等等，在MATLAB中，只需用“”或“”就解决问题。它内部实际包含着许许多多的自适应算法，如对超定方程用最小二乘法，对欠定方程时它将给出范数最小的一个解，解三对角阵方程组时用追赶法等等。格式：x=Ab,36,例：解方程组A=.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129;.3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927;b=0.40430.15500.4240-0.2557;x=Ab;xans=-0.1819-1.66302.2172-0.4467,37,4.2.2线性方程组直接求解判定求解,38,39,例：A=1234;4321;1324;4132;B=51;42;33;24;C=AB;rank(A),rank(C)ans=4ans=4x=inv(A)*B%ABx=-1.80002.40001.8667-1.26673.8667-3.2667-2.13332.7333,40,检验norm(A*x-B)ans=7.4738e-015精确解x1=inv(sym(A)*Bx1=-9/5,12/528/15,-19/1558/15,-49/15-32/15,41/15检验norm(double(A*x1-B)ans=0,41,原方程组对应的齐次方程组的解求取A矩阵的化零矩阵：格式：Z=null(A)求取A矩阵的化零矩阵的规范形式：格式：Z=null(A,r),42,null是用来求齐次线性方程组的基础解系的，加上r则求出的是一组最小正整数解，如果不加，则求出的是解空间的规范正交基。,例：判断可解性A=1234;2211;2468;4422;B=1;3;2;6;C=AB;rank(A),rank(C)ans=22Z=null(A,r)%解出规范化的化零空间Z=2.00003.0000-2.5000-3.50001.0000001.0000,43,x0=pinv(A)*B%得出一个特解x0=0.95420.7328%全部解-0.0763-0.2977验证得出的解a1=randn(1);a2=rand(1);%取不同分布的随机数x=a1*Z(:,1)+a2*Z(:,2)+x0;norm(A*x-B)ans=4.4409e-015,44,解析解Z=null(sym(A)Z=2,3-5/2,-7/21,00,1x0=sym(pinv(A)*B)x0=125/13196/131-10/131-39/131,45,验证得出的解a1=randn(1);a2=rand(1);%取不同分布的随机数x=a1*Z(:,1)+a2*Z(:,2)+x0;norm(double(A*x-B)ans=0通解symsa1a2;x=a1*Z(:,1)+a2*Z(:,2)+x0 x=2*a1+3*a2+125/131-5/2*a1-7/2*a2+96/131a1-10/131a2-39/131,46,摩尔彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。,47,4.2.3线性方程组的直接求解分析,LU分解,48,49,50,格式l,u,p=lu(A)L是一个单位下三角矩阵，u是一个上三角矩阵，p是代表选主元的置换矩阵。故：Ax=y=PAx=Py=LUx=Py=PA=LUl,u=lu(A)其中l等于P-1L，u等于U，所以(P-1L)U=A,51,例：对A进行LU分解A=123;241;467;l,u,p=lu(A)l=1.0000000.50001.000000.25000.50001.0000u=4.00006.00007.000001.0000-2.5000002.5000p=001010100,52,l,u=lu(A)lP-1Ll=0.25000.50001.00000.50001.000001.000000u=4.00006.00007.000001.0000-2.5000002.5000,53,QR分解将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。求得正交矩阵Q和上三角阵R，Q和R满足A=QR。,54,格式：Q,R=qr(A),例：A=123;456;789;101112;Q,R=qr(A)Q=-0.0776-0.83310.5456-0.0478-0.3105-0.4512-0.69190.4704-0.5433-0.0694-0.2531-0.7975-0.77620.31240.39940.3748R=-12.8841-14.5916-16.29920-1.0413-2.082600-0.0000000,55,Cholesky(乔里斯基)分解若矩阵A为n阶对称正定阵,则存在唯一的对角元素为正的三角阵D，使得,56,57,格式：D=chol(A),58,例：进行Cholesky分解。A=1648;45-4;8-422;D=chol(A)D=41202-3003,利用矩阵的LU、QR和cholesky分解求方程组的解,（1）LU分解：A*X=b变成L*U*X=b所以X=U(Lb)这样可以大大提高运算速度。例：求方程组的一个特解。解：A=42-1;3-12;1130;B=2108;D=det(A)D=0,59,L,U=lu(A)L=0.3636-0.50001.00000.27271.000001.000000U=11.00003.000000-1.81822.0000000.0000,60,X=U(LB)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=2.018587e-017.X=1.0e+016*%结果中的警告是由于系数行列式为零产生的。-0.4053%可以通过A*X验证其正确性。1.48621.3511A*X%Matlab7.0显示没有解ans=088,61,（2）Cholesky分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，方程A*X=b变成R*R*X=b所以X=R(Rb)（3）QR分解对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR方程A*X=b变形成QRX=b所以X=R(Qb)这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。,62,三个变换在线性方程组的迭代求解中，要用到系数矩阵A的上三角矩阵、对角阵和下三角矩阵。此三个变换在MATLAB中可由以下函数实现。上三角变换：格式triu(A,1)对角变换：格式diag(A)下三角变换：格式tril(A,-1)例：对此矩阵做三种变换。,63,A=12-2;111;221;triu(A,1)ans=02-2001000tril(A,-1)ans=000100220b=diag(A);bans=111,64,4.4线性方程组的符号解法,在MATLAB的SymbolicToolbox中提供了线性方程的符号求解函数，如linsolve(A,b)等同于X=sym(A)sym(b).solve(eqn1,eqn2,.,eqnN,var1,var2,.,varN),65,例：A=sym(10,-1,0;-1,10,-2;0,-2,10)；b=(9;7;6)；linsolve(A,b)ans=473/47591/95376/475vpa(ans)ans=.99578947368421052631578947368421.95789473684210526315789473684211.79157894736842105263157894736842,66,例：x,y=solve(x2+x*y+y=3,x2-4*x+3=0,x,y)x=13y=1-3/2,67,4.3迭代解法的几种形式4.3.1Jacobi迭代法,方程组Ax=bA可写成A=D-L-U其中:D=diaga11,a22,ann,-L、-U分别为A的严格下、上三角部分（不包括对角线元素）.由Ax=bx=Bx+f由此可构造迭代法：x(k+1)=Bx(k)+f其中：B=D-1(L+U)=I-D-1A,f=D-1b.,functiony=jacobi(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);B=D(L+U);f=Db;y=B*x0+f;n=1;whilenorm(y-x0)=1.0e-6x0=y;y=B*x0+f;n=n+1;endn,例：用Jacobi方法求解，设x(0)=0,精度为10-6。a=10-10;-110-2;0-210;b=9;7;6;jacobi(a,b,0;0;0)n=11ans=0.99580.95790.7916,4.3.2Gauss-Seidel迭代法,由原方程构造迭代方程x(k+1)=Gx(k)+f其中：G=(D-L)-1U,f=(D-L)-1bD=diaga11,a22,ann,-L、-U分别为A的严格下、上三角部分（不包括对角线元素）.,functiony=seidel(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G=(D-L)U;f=(D-L)b;y=G*x0+f;n=1;whilenorm(y-x0)=1.0e-6x0=y;y=G*x0+f;n=n+1;endn,例：对上例用Gauss-Seidel迭代法求解a=10-10;-110-2;0-210;b=9;7;6;seidel(a,b,0;0;0)n=7ans=0.99580.95790.7916例：分别用Jacobi和G-S法迭代求解,看是否收敛。,a=12-2;111;221;b=9;7;6;jacobi(a,b,0;0;0)n=4ans=-27268seidel(a,b,0;0;0)n=1011ans=1.0e+305*-InfInf-1.7556,4.3.3SOR迭代法,在很多情况下，J法和G-S法收敛较慢，所以考虑对G-S法进行改进。于是引入一种新的迭代法逐次超松弛迭代法（SuccesiseOver-Relaxation),记为SQR法。迭代公式为：X(k+1)=(D-wL)-1(1-w)D+wU)x(k)+w(D-wL)-1b其中：w最佳值在1,2）之间，不易计算得到，因此w通常有经验给出。,functiony=sor(a,b,w,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);M=(D-w*L)(1-w)*D+w*U);f=(D-w*L)b*w;y=M*x0+f;n=1;whilenorm(y-x0)=1.0e-6x0=y;y=M*x0+f;n=n+1;endn,例：上例中，当w=1.103时，用SOR法求解原方程。a=10-10;-110-2;0-210;b=9;7;6;sor(a,b,1.103,0;0;0)n=8ans=0.99580.95790.7916,4.3.4两步迭代法,当线性方程系数矩阵为对称正定时，可用一种特殊的迭代法来解决，其迭代公式为：（D-L）x(k+1/2)=Ux(k)+b（D-U）x(k+1)=Lx(k+1/2)+b=x(k+1/2)=(D-L)-1Ux(k)+(D-L)-1bx(k+1)=(D-U)-1Lx(k+1/2)+(D-U)-1b,functiony=twostp(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G1=(D-L)U;f1=(D-L)b;G2=(D-U)L;f1=(D-U)b;y=G1*x0+f1;y=G2*y+f2;n=1;whilenorm(y-x0)=1.0e-6x0=y;y=G1*x0+f1;y=G2*y+f2;n=n+1;endn,例：求解方程组a=10-120;-111-13;2-1103;03-18;b=6;25;-11;15;twostp(a,b,0;0;0;0)n=7ans=1.07911.9824-1.40440.9560,4.5稀疏矩阵技术,稀疏矩阵的建立：格式S=sparse(i,j,s,m,n)生成一mxn阶的稀疏矩阵，以向量i和j为坐标的位置上对应元素值为s。例：n=5;a1=sparse(1:n,1:n,4*ones(1,n),n,n)a1=(1,1)4(2,2)4(3,3)4(4,4)4(5,5)4,81,例：a2=sparse(2:n,1:n-1,ones(1,n-1),n,n)a2=(2,1)1(3,2)1(4,3)1(5,4)1full(a2)ans=0000010000010000010000010,82,例：n=5,建立主对角线上元素为4，两条次对角线为1的三对角阵。n=5;a1=sparse(1:n,1:n,4*ones(1,n),n,n);a2=sparse(2:n,1:n-1,ones(1,n-1),n,n);a=a1+a2+a2a=(1,1)4(2,1)1(1,2)1(2,2)4(3,2)1(2,3)1(3,3)4(4,3)1,83,(3,4)1(4,4)4(5,4)1(4,5)1(5,5)4full(a)ans=4100014100014100014100014,84,格式A=spdiags(B,d,m,n)生成一mxn阶的稀疏矩阵，使得B的列放在由d指定的位置。例：n=5b=spdiags(ones(n,1),4*ones(n,1),ones(n,1),-1,0,1,n,n);full(b)ans=4100014100014100014100014,85,格式：spconvert(dd)对于无规律的稀疏矩阵，可使用此命令由外部数据转化为稀疏矩阵。调用形式为：先用load函数加载以行表示对应位置和元素值的.dat文本文件，再用此命令转化为稀疏矩阵。例：无规律稀疏矩阵的建立。首先编制文本文件sp.dat如下：515.00358.00442.00550,86,loadsp.datspconvert(sp)ans=(5,1)5(4,4)2(3,5)8full(ans)ans=0000000000000080002050000,87,88,稀疏矩阵的计算：同满矩阵比较，稀疏矩阵在算法上有很大的不同。具体表现在存储空间减少，计算时间减少。例：比较求解下面方程组n1000时两种方法的差别。,n=1000;a1=sparse(1:n,1:n,4*ones(1,n),n,n);a2=sparse(2:n,1:n-1,ones(1,