MATLAB入门线性代数.ppt

《MATLAB入门线性代数.ppt》由会员分享，可在线阅读，更多相关《MATLAB入门线性代数.ppt（41页珍藏版）》请在三一办公上搜索。

1、第三章 线性代数,3.1 常用矩阵函数,norm：矩阵或向量范数,For vectors.norm(V,P)=sum(abs(V).P)(1/P).norm(V)=norm(V,2).norm(V,inf)=max(abs(V).norm(V,-inf)=min(abs(V).,例 x=1 2 3 4 5;x=1 20 300 40 5;x=10000 2000 300 40 5;norm(x,1),norm(x,2),norm(x,3),norm(x,inf),For matrices.norm(X)is the largest singular value of X,max(svd(X).

2、norm(X,2)is the same as norm(X).norm(X,1)is the 1-norm of X,the largest column sum,=max(sum(abs(X).norm(X,inf)is the infinity norm of X,the largest row sum,=max(sum(abs(X).norm(X,fro)is the Frobenius norm,sqrt(sum(diag(X*X).norm(X,P)is available for matrix X only if P is 1,2,inf or fro.,例 x=1 20 300

3、;400 50 6;70 800 9;norm(x,1),norm(x,2),norm(x,inf),norm(x,fro),dot(x,y)向量的内积 det(A)方阵的行列式；rank(A)矩阵的秩；trace(A)矩阵的迹；rref(A)初等变换化矩阵A为阶梯矩阵inv(A)矩阵的逆；即 A-1pinv(A)矩阵的广义逆A+null(A)零空间的基阵 roth(A)值空间的基阵 orth(A)将A标准正交化 cond(A,flag)矩阵的条件数,flag=2,1,inf,fro;,例：分别求x=1 3 7 8-2,y=3 9 3-3 9的长度与它们的夹角。x=1 3 7 8-2;y=3

4、9 3-3 9;xx=norm(x,2);yy=norm(y,2);theta=acos(dot(x,y)/(xx*yy);s=xx,yy,theta,例：给定一组线性无关的向量，将其标准正交化a=magic(5);b=orth(a),d=eig(A)：方阵的特征值；V,D=eig(A):A*V=V*D V,J=jordan(A):A*V=V*J c=condeig(A)：向量c中包含矩阵A关于各 特征值的条件数 V,D,c=condeig(A)：,例:A=1 0 0;1 2 0;1 2 3,d=eig(A),V,D=eig(A),C=condeig(A),V,D,C=condeig(A),例

5、：观察7阶随机矩阵特征值的分布a=rands(7,7)%产生7阶随机矩阵e=eig(a)title(特征值的分布);plot(real(e),imag(e),o)xlabel(实轴);ylabel(虚轴);,注：本例验证了如下定理：实方阵的特征值或为实 数或呈共轭对出现。,例：观察正交矩阵的特征值分布a=rands(7,7);b=orth(a);%构造一个正交矩阵theta=0:0.01:2*pi;e=eig(b);plot(real(e),imag(e),r*,cos(theta),sin(theta);axis equaltitle(正交矩阵特征值的分布);xlabel(实轴);ylabe

6、l(虚轴);,注：本例验证了正交矩阵的特征值分布在复平面的单 位圆上。,3.2 矩阵的运算,一、矩阵的转置、乘积，逆,A=1 0 0;1 2 0;1 2 3,A_trans=AH=1 2 3;2 1 0;1 2 3,K=1 2 3;2 1 0;2 3 1HK=H*KH_inv=inv(H),K_inv=K-1,二、矩阵的左除和右除,左除“”：求矩阵方程AX=B的解；（A、B的行要保持一致）解为 X=AB；当A为方阵且可逆时有X=AB=inv(A)*B；右除“/”：求矩阵方程XA=B的解（A、B的列要保持一致）解为 X=B/A，当A为方阵且可逆时有X=B/A=B*inv(A),【例】“求逆”法和

7、“左除”法解恰定方程的性能对比（1）构造一个条件数很大的高阶恰定方程randn(state,0);A=gallery(randsvd,100,2e13,2);x=ones(100,1);b=A*x;cond(A)ans=1.9990e+013,（2）用“求逆”法求解ticxi=inv(A)*b;ti=toceri=norm(x-xi)rei=norm(A*xi-b)/norm(b),ti=0.4400eri=0.0469rei=0.0047,（3）用“左除”法求解tic;xd=Ab;td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)td=0.0600e

8、rd=0.0078red=2.6829e-015,ti=0.4400eri=0.0469rei=0.0047,求矩阵方程：设A、B满足关系式：AB2B+A,求B。其中A=3 0 1;1 1 0;0 1 4。解：有(A-2I)BA程序：A=3 0 1;1 1 0;0 1 4;B=inv(A-2*eye(3)*A,BB=(A-2*eye(3)A观察结果：,三.矩阵函数的计算,给定n阶方阵A，求exp(A),sqrt(A),log(A)矩阵函数的Matlab函数分别为 expm(A),sqrtm(A),logm(A),例：求三个特殊矩阵函数a=magic(3);e=expm(a)%求矩阵函数ee=e

9、xp(a)%求矩阵中每个元素的函数值s=sqrtm(a)ss=sqrt(a)l=logm(a)ll=log(a),L,U,P=lu(A)：PA=LUr=chol(A):A=LLT,r=LT Q,R=qr(A),Q,R=qr(A,0):A=QRU,S,V=svd(A):A=USVTQ,R=schur(A):QTAQ=RP,H=hess(A):PAP-1=H,3.3 常用矩阵分解函数,实方阵的初等化简分解:A1=3 1 3;2 5 2;1 2 3,L,U=lu(A1)A2=1 1 2;1 2 3;1 2 1;1 1 6,Q,R=qr(A2),Q1,R1=qr(A2,0)A3=2 1 1;1 4-1

10、;1-1 3,r=chol(A3),奇异值分解:A=1 1;1 1;0 0,U,D,V=svd(A),U1,D1,V1=svd(A,0),实方阵的正交相似化简:Q,R=schur(A)A1=3 1 0;-4-1 0;4 8-2,Q1,R1=schur(A1)A2=9-31 49 30;1 0 0 0;1 1 0 0;0 0 1 0,Q2,R2=schur(A2)A3=2 1 1;1 4-1;1-1 3,Q3,R3=schur(A3),实方阵的上Hessenberg分解:P,H=hess(A)A1=3 1 0;-4-1 0;4 8-2,P,H=hess(A1)P*H*inv(P),3.4 稀疏矩

11、阵,一.稀疏矩阵的特点和存储,在编辑器内建立一个数据文本文件:st.m行 列 aij1 1 32 1 21 2 12 2 13 2 54 2 1 3 43 3 2,load st.msa=spconvert(st)sa=(1,1)3(2,1)2(1,2)1(2,2)1(3,2)5(4,2)1(2,3)4(3,3)2,二.稀疏矩阵的运算,例:a=full(sa)a=3 1 0 2 1 4 0 5 2 0 1 0pa=0 1 0 0;1 0 0 00 0 0 1;0 0 1 0pa=0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0,pa*saans=2 1 4 3 1 0 0 1 0

12、 0 5 2,1.稀疏矩阵的初等变换,行初等变换pv=2 1 4 3;s1=sa(pv,:),列初等变换pv=2 1 3;s2=sa(:,pv),s1=spa*sas1=(1,1)2(2,1)3(1,2)1(2,2)1(3,2)1(4,2)5(1,3)4(4,3)2,spa=sparse(pa)spa=(2,1)1(1,2)1(4,3)1(3,4)1,2.稀疏矩阵的分解函数,L,U,P=lu(sa)L,U=luinc(sa,0)r=chol(sa)r=cholinc(sa,0)Q,R=qr(sa)eigs(sa)svds(sa),例A=4 1 0 0 0;1 4 1 0 0;0 1 4 1 0

13、;0 0 1 4 1;0 0 0 1 4;sa=sparse(A);L,U,P=lu(sa),用下面例子说明稀疏矩阵的处理优点,xs.m%设n=3000,输入n=3000;b=1:n;a1=sparse(1:n,1:n,1,n,n);a2=sparse(2:n,1:n-1,1,n,n);a=4*a1+a2+a2;%输出用稀疏矩阵求解的时间t1tic;x=ab;t1=tocaa=full(a);%输出用满矩阵求解的时间t2tic;xx=aab;t2=toc%为检验x与xx是否相同分别输出其分量之和y=sum(x)yy=sum(xx)%结果y与yy相同，而t1与t2相差巨大,1.不用预优矩阵的共轭

14、斜量法 x=pcg(a,b,tol,kmax)2.用预优矩阵的共轭斜量法(1)x=pcg(a,b,tol,kmax,m)(2)r=chol(m)x=pcg(a,b,tol,kmax,r,r,x0)3.未给定预优矩阵的共轭斜量法 r=cholinc(sa,0)x=pcg(a,b,tol,kmax,r,r,x0),3.5 用预条件共轭斜量法求解线性方程组,其中a,m为对称正定矩阵,1.不用预优矩阵的双共轭梯度法 x=bicg(a,b,tol,kmax)2.用预优矩阵的双共轭梯度法 l,u=lu(m)x=bicg(a,b,tol,kmax,l,u,x0)3.未给定预优矩阵的双共轭梯度法 l,u=lu

15、inc(sa,0)x=bicg(a,b,tol,kmax,l,u,x0),3.6 求解大型稀疏非对称正定的线性方程组,一.双共轭梯度法,二.广义极小残差法 gmres,3.7 不可解问题,例：考察下面三个线性方程组的解(A1)-x+y=1(A2)-x+y=1(A3)x+2y=-2-2x+2y=2-x+y=0-x+y=1 2x-y=0,A1=-1 1;-2 2,b1=1;2,x1=A1b1A10=-1 1,b10=1,x10=A10b10,A2=-1 1;-1 1,b2=1;0,x2=A2b2,A3=1 2;-1 1;2-1,b3=-2;1;0,x3=A3b3,注：系数矩阵A必须是行满秩或列满秩

16、,3.8 病态问题,例：考察舍入误差对线性方程组解的影响 0.12065x+0.98775y=2.01045 0.12032x+0.98755y=2.00555,A=0.12065 0.98775;0.12032 0.98755,b=2.01045;2.00555,x=Ab,A=0.12065 0.98775;0.12032 0.98755,b1=2.01145;2.00555,x1=Ab1,A1=0.12165 0.98775;0.12032 0.98755,b=2.01144;2.00555,x2=A1b,3.9 关于多项式的MATLAB命令,一、多项式表达方式的约定多项式 用行向量表示

17、用比较习惯的方式显示多项式：pp=poly2str(p,x),【例】多项式 可表示为 p=2 1 4 5 pp=poly2str(p,x),二、多项式运算函数,r=roots(p)：求多项式的零点p=poly(r):以r为零点的多项式p=poly(A):A的特征多项式PA=polyval(p,S):按数组运算规则，计算多项式的值 其中S，PA为矩阵PM=polyvalm(p,S):按矩阵运算规则，计算多项式的值,其中S，PM为矩阵 p=conv(p1,p2)：多项式的乘积q,r=deconv(p1,p2)：多项式的除法，p1/p2 p1(x)=p2(x)q(x)+r(x),【例】由给定根向量求

18、多项式系数向量。R=-0.5,-0.3+0.4*i,-0.3-0.4*i;P=poly(R)PPR=poly2str(P,x)P=1.0000 1.1000 0.5500 0.1250PPR=x3+1.1 x2+0.55 x+0.125,【例】求多项式 的零点。r=roots(1-6 15-20 15-6 1)r=1.0042+0.0025i 1.0042-0.0025i 1.0000+0.0049i 1.0000-0.0049i 0.9958+0.0024i 0.9958-0.0024i,注：尽管利用MATLAB使得从系数转换到零点或从零点转换到系数都非常容易，但是使用时一定要注意计算的精度

19、。如果存在重根，这种转换可能会降低精度。对于数值计算，计算重根是最困难的问题之一。,【例】求3阶方阵A的特征多项式。A=11 12 13;14 15 16;17 18 19;PA=poly(A)PPA=poly2str(PA,x)PA=1.0000-45.0000-18.0000-0.0000PPA=x3-45 x2-18 x-2.8387e-015,【例】求 的“商”及“余”多项式。p1=conv(1,0,2,conv(1,4,1,1);p2=1 0 1 1;q,r=deconv(p1,p2);cq=商多项式为;cr=余多项式为;disp(cq,poly2str(q,x)disp(cr,po

20、ly2str(r,x)商多项式为 x+5余多项式为 5 x2+4 x+3,【例】两种多项式求值指令的差别。S=pascal(4)P=poly(S);PP=poly2str(P,s)PA=polyval(P,S)PM=polyvalm(P,S),S=1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20PP=s4-29 s3+72 s2-29 s+1,PM=1.0e-010*0.0016 0.0033 0.0090 0.0205 0.0045 0.0101 0.0286 0.0697 0.0095 0.0210 0.0653 0.1596 0.0163 0.0387 0.1226

21、0.3019,PA=1.0e+004*0.0016 0.0016 0.0016 0.0016 0.0016 0.0015-0.0140-0.0563 0.0016-0.0140-0.2549-1.2089 0.0016-0.0563-1.2089-4.3779,可以用命令polyval或 polyvalm计算多项式的值,p=2 1 4 5;xi=2.5;yi=polyval(p,xi)yi=52.5000 yi=polyvalm(p,xi)yi=52.5000,【例】求多项式的积分 poly_itg.m function py=poly_itg(p)n=length(p)py=p.*n:-1:1.(-1),0,【例】将Chebyshev多项式展开为幂级数形式 Cheby_pw.m function pn=Cheby_pw(n)pbb=1;if n=0,pn=pbb;break;end pb=1 0;if n=1,pn=pb;break;end for i=2:n pn=2*pb,0-0,0,pbb;pbb=pb;pb=pn;end,