第四章 线性代数问题求解.ppt

上传人:bonesoil321 文档编号:374244 上传时间:2018-10-05 格式:PPT 页数:94 大小:918KB
下载 相关 举报
第四章 线性代数问题求解.ppt_第1页
第1页 / 共94页
第四章 线性代数问题求解.ppt_第2页
第2页 / 共94页
第四章 线性代数问题求解.ppt_第3页
第3页 / 共94页
第四章 线性代数问题求解.ppt_第4页
第4页 / 共94页
第四章 线性代数问题求解.ppt_第5页
第5页 / 共94页
亲,该文档总共94页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述

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),随机元素矩阵若矩阵随机元素满足0,1区间上的均匀分布生成nm阶标准均匀分布为随机数矩阵:A=rand(n,m)生成nn阶标准均匀分布为随机数方阵:A=rand(n),对角

2、元素矩阵已知向量生成对角矩阵:A=diag(V)已知矩阵提取对角元素列向量:Vdiag(A)生成主对角线上第k条对角线为V的矩阵:A=diag(V,k),例:diag( )函数的不同调用格式 C=1 2 3; V=diag(C) % 生成对角矩阵 V =1 0 00 2 00 0 3 V1=diag(V) % 将列向量通过转置变换成行向量 V1 =1 2 3 C=1 2 3; V=diag(C,2) % 主对角线上第 k条对角线为C的矩阵 V =0 0 1 0 00 0 0 2 00 0 0 0 30 0 0 0 00 0 0 0 0,生成三对角矩阵: V=diag(1 2 3 4)+diag

3、(2 3 4,1)+diag(5 4 3,-1)V =1 2 0 05 2 3 00 4 3 40 0 3 4,Hilbert矩阵及逆Hilbert矩阵生成n阶的Hilbert矩阵: A=hilb(n)求取逆Hilbert矩阵:B=invhilb(n),Hankel(汉克 ) 矩阵其中:第一列的各个元素定义为C向量,最后一行各个元素定义为R。H为对称阵。H1=hankel(C)由 Hankel 矩阵反对角线上元素相等得出一下三角阵均为零的Hankel 矩阵,Vandermonde(范德蒙)矩阵,伴随矩阵其中:P(s)为首项系数为一的多向式。,符号矩阵的输入数值矩阵A转换成符号矩阵:B=sym(

4、A) 例: A=hilb(3) A =1.0000 0.5000 0.33330.5000 0.3333 0.25000.3333 0.2500 0.2000 B=sym(A) B = 1, 1/2, 1/3 1/2, 1/3, 1/4 1/3, 1/4, 1/5,4.1.2 矩阵基本概念与性质,行列式格式 :d=det(A)例:求行列式 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; det(A) ans =0,例: tic, A=sym(hilb(20); det(A), toc ans = 1/23774547167685345090916442

5、43427616440175419837753486493033185331234419759310644585187585766816573773440565759867265558971765638419710793303386582324149811241023554489166154717809635257797836800000000000000000000000000000000000 elapsed_time =2.3140 高阶的Hilbert矩阵是接近奇异的矩阵。,矩阵的迹格式: t=trace(A)矩阵的秩格式:r=rank(A) 用默认的精度求数值秩r=rank(A, )

6、 给定精度下求数值秩矩阵的秩也表示该矩阵中行列式不等于0的子式的最大阶次。可证行秩和列秩(线性无关的)应相等。,例 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; rank(A) ans =3 该矩阵的秩为3,小于矩阵的阶次,故为非满秩矩阵。例 H=hilb(20); rank(H) 数值方法 ans =13 H=sym(hilb(20); rank(H) % 解析方法,原矩阵为非奇异矩阵 ans = 20,矩阵范数,矩阵的范数定义:格式:N=norm(A) 求解默认的2范数N=norm(A,选项) 选项可为1,2,inf等,例:求一向量、矩阵的范数

7、a=16 2 3 13; norm(a), norm(a,2), norm(a,1), norm(a,Inf) ans =2.092844953645635e+001 2.092844953645635e+001 3.400000000000000e+001 1.600000000000000e+001 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; norm(A), norm(A,2), norm(A,1), norm(A,Inf) ans =34 34 34 34符号运算工具箱未提供norm( )函数,需先用double( )函数转换成双精度数值矩

8、阵,再调用norm( )函数。,特征多项式格式: C=poly(A)例: A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; poly(A) 直接求取 ans = 1.000000000000000e+000 -3.399999999999999e+001 -7.999999999999986e+001 2.719999999999999e+003-2.819840539024018e-012 A=sym(A); poly(A) 运用符号工具箱ans =x4-34*x3-80*x2+2720*x,矩阵多项式的求解,符号多项式与数值多项式的转换格式:f=po

9、ly2sym(P) 或 f=poly2sym(P,x) 格式: P=sym2poly(f),例: P=1 2 3 4 5 6; % 先由系数按降幂顺序排列表示多项式 f=poly2sym(P,v) % 以 v 为算子表示多项式 f = v5+2*v4+3*v3+4*v2+5*v+6 P=sym2poly(f) P =1 2 3 4 5 6,矩阵的逆矩阵格式: C=inv(A)例: format long; H=hilb(4); H1=inv(H) H1 =1.0e+003 *0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999

10、999999999-0.11999999999999 1.19999999999990 -2.69999999999976 1.679999999999840.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961-0.13999999999999 1.67999999999984 -4.19999999999961 2.79999999999974,检验: H*H1 ans =1.00000000000001 0.00000000000023 -0.00000000000045 0.000000000000230

11、.00000000000001 1.00000000000011 -0.00000000000011 0.000000000000110.00000000000001 0 1.00000000000011 00.00000000000000 0.00000000000011 -0.00000000000011 1.00000000000011 计算误差范数: norm(H*inv(H)-eye(size(H) ans =6.235798190375727e-013 H2=invhilb(4); norm(H*H2-eye(size(H) ans =5.684341886080802e-014,

12、 H=hilb(10); H1=inv(H); norm(H*H1-eye(size(H) ans =0.00264500826202 H2=invhilb(10); norm(H*H2-eye(size(H) ans =1.612897415528547e-005 H=hilb(13); H1=inv(H); norm(H*H1-eye(size(H) Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 2.339949e-018. ans =53.2369600857

13、0294 H2=invhilb(13); norm(H*H2-eye(size(H) ans =11.37062973181391 对接近于奇异矩阵,高阶一般不建议用inv( ),可用符号工具箱。, H=sym(hilb(7); inv(H)ans = 49, -1176, 8820, -29400, 48510, -38808, 12012 -1176, 37632, -317520, 1128960, -1940400, 1596672, -504504 8820, -317520, 2857680, -10584000, 18711000, -15717240, 5045040 -294

14、00, 1128960, -10584000, 40320000, -72765000, 62092800, -20180160 48510, -1940400, 18711000, -72765000, 133402500, -115259760, 37837800 -38808, 1596672, -15717240, 62092800, -115259760, 100590336, -33297264 12012, -504504, 5045040, -20180160, 37837800, -33297264, 11099088 H=sym(hilb(30); norm(double(

15、H*inv(H)-eye(size(H) ans =0,例:奇异阵求逆 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; format long; B = inv(A) Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 1.306145e-017. B =1.0e+014 *0.93824992236885 2.81474976710656 -2.81474976710656 -0.938249922368852.81474976

16、710656 8.44424930131968 -8.44424930131968 -2.81474976710656-2.81474976710656 -8.44424930131968 8.44424930131968 2.81474976710656-0.93824992236885 -2.81474976710656 2.81474976710656 0.93824992236885 norm(A*B-eye(size(A) 检验 ans =1.64081513306419 A=sym(A); inv(A) 奇异矩阵不存在一个相应的逆矩阵,用符号工具箱的函数也不行 ? Error us

17、ing = sym/inv Error, (in inverse) singular matrix,同样适用于含有变量的矩阵求逆。 例: syms a1 a2 a3 a4; C=a1 a2;a3 a4; inv(C)ans = -a4/(-a1*a4+a2*a3), a2/(-a1*a4+a2*a3) a3/(-a1*a4+a2*a3), -a1/(-a1*a4+a2*a3),矩阵的相似变换与正交矩阵其中:A为一方阵,B矩阵非奇异。相似变换后,X矩阵的秩、迹、行列式与特征值等均不发生变化,其值与A矩阵完全一致。对于一类特殊的相似变换满足如下条件,称为正交基矩阵。,例: A=5,9,8,3; 0

18、,3,2,4; 2,3,5,9; 3,4,5,8; Q=orth(A) Q =-0.6197 0.7738 -0.0262 -0.1286-0.2548 -0.1551 0.9490 0.1017-0.5198 -0.5298 -0.1563 -0.6517-0.5300 -0.3106 -0.2725 0.7406 norm(Q*Q-eye(4) ans =4.6395e-016 norm(Q*Q-eye(4) ans =4.9270e-016,例: A=16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1; Q=orth(A) A为奇异矩阵,故得出的Q为长方形

19、矩阵 Q =-0.5000 0.6708 0.5000-0.5000 -0.2236 -0.5000-0.5000 0.2236 -0.5000-0.5000 -0.6708 0.5000 norm(Q*Q-eye(3) ans =1.0140e-015,4.2 线性方程组直接解法 4.2.1线性方程组直接求解矩阵除法,关于线性方程组的直接解法,如Gauss消去法、选主元消去法、平方根法、追赶法等等,在MATLAB中,只需用“”或“”就解决问题。它内部实际包含着许许多多的自适应算法,如对超定方程用最小二乘法,对欠定方程时它将给出范数最小的一个解,解三对角阵方程组时用追赶法等等。格式: x=Ab

20、,例:解方程组 A=.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129; .3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927; b=0.4043 0.1550 0.4240 -0.2557; x=Ab; x ans =-0.1819 -1.6630 2.2172 -0.4467,4.2.2线性方程组直接求解判定求解,例: A=1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2; B=5 1; 4 2; 3 3; 2 4; C=A B; rank(A), rank(C) ans =4 ans

21、 =4 x=inv(A)*B x =-1.8000 2.40001.8667 -1.26673.8667 -3.2667-2.1333 2.7333,检验 norm(A*x-B) ans =7.4738e-015 精确解 x1=inv(sym(A)*Bx1 = -9/5, 12/5 28/15, -19/15 58/15, -49/15 -32/15, 41/15 检验 norm(double(A*x1-B) ans =0,原方程组对应的齐次方程组的解求取A矩阵的化零矩阵:格式: Z=null(A)求取A矩阵的化零矩阵的规范形式:格式: Z=null(A, r ),例:判断可解性 A=1 2

22、3 4; 2 2 1 1; 2 4 6 8; 4 4 2 2; B=1;3;2;6; C=A B; rank(A), rank(C) ans =2 2 Z=null(A,r) % 解出规范化的化零空间 Z =2.0000 3.0000-2.5000 -3.5000 1.0000 00 1.0000, 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.4

23、409e-015,解析解 Z=null(sym(A) Z = 2, 3 -5/2, -7/2 1, 0 0, 1 x0=sym(pinv(A)*B) x0 = 125/131 96/131 -10/131 -39/131,验证得出的解 a1=randn(1); a2=rand(1); % 取不同分布的随机数 x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(double(A*x-B) ans =0 通解 syms a1 a2; x=a1*Z(:,1)+a2*Z(:,2)+x0 x = 2*a1+3*a2+125/131 -5/2*a1-7/2*a2+96/131 a1-10/131

24、 a2-39/131,摩尔彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。,4.2.3 线性方程组的直接求解分析,LU分解,格式l,u,p=lu(A) L是一个单位下三角矩阵,u是一个上三角矩阵,p是代表选主元的置换矩阵。 故:Ax=y = PAx=Py = LUx=Py = PA=LUl,u=lu(A) 其中l等于P-1 L,u等于U,所以(P-1 L)U=A,例:对A进行LU分解 A=1 2 3; 2 4 1; 4 6 7; l,u,p=lu(A) l =1.0000 0 00.5000 1.0000 00.2500 0.5000 1.0000 u =4.0000 6.0000 7.

25、00000 1.0000 -2.50000 0 2.5000 p =0 0 10 1 01 0 0, l,u=lu(A) lP-1 Ll =0.2500 0.5000 1.00000.5000 1.0000 01.0000 0 0u =4.0000 6.0000 7.00000 1.0000 -2.50000 0 2.5000,QR分解 将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。求得正交矩阵Q和上三角阵R,Q和R满足A=QR。,格式:Q,R = qr(A),例: A = 1 2 3;4 5 6; 7 8 9; 10 11 12; Q,R = qr(A) Q =-0.0776 -0.83

26、31 0.5456 -0.0478-0.3105 -0.4512 -0.6919 0.4704-0.5433 -0.0694 -0.2531 -0.7975-0.7762 0.3124 0.3994 0.3748 R =-12.8841 -14.5916 -16.29920 -1.0413 -2.08260 0 -0.00000 0 0,Cholesky(乔里斯基 )分解若矩阵A为 n阶对称正定阵,则存在唯一的对角元素为正的三角阵D,使得,格式:D=chol(A),例:进行Cholesky分解。 A=16 4 8; 4 5 -4; 8 -4 22; D=chol(A) D =4 1 20 2

27、-30 0 3,利用矩阵的LU、QR和cholesky分解求方程组的解,(1)LU分解:A*X=b 变成 L*U*X=b 所以 X=U(Lb) 这样可以大大提高运算速度。例:求方程组 的一个特解。解: A=4 2 -1;3 -1 2;11 3 0; B=2 10 8; D=det(A) D =0, L,U=lu(A) L =0.3636 -0.5000 1.00000.2727 1.0000 01.0000 0 0 U =11.0000 3.0000 00 -1.8182 2.00000 0 0.0000, X=U(LB) Warning: Matrix is close to singula

28、r or badly scaled.Results may be inaccurate. RCOND = 2.018587e-017. X =1.0e+016 * % 结果中的警告是由于系数行列式为零产生的。 -0.4053 % 可以通过A*X验证其正确性。 1.48621.3511 A*X ans =088,(2)Cholesky分解若A为对称正定矩阵,则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积, 方程 A*X=b 变成 R*R*X=b 所以 X=R(Rb) (3)QR分解对于任何长方矩阵A,都可以进行QR分解,其中Q为正交矩阵,R为上三角矩阵的初等变换形式,即:A=QR

29、 方程 A*X=b 变形成 QRX=b 所以 X=R(Qb)这三种分解,在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。,三个变换在线性方程组的迭代求解中,要用到系数矩阵A的上三角矩阵、对角阵和下三角矩阵。此三个变换在MATLAB中可由以下函数实现。 上三角变换:格式 triu(A,1) 对角变换:格式 diag(A) 下三角变换:格式 tril(A,-1)例:对此矩阵做三种变换。, A=1 2 -2;1 1 1;2 2 1; triu(A,1) ans =0 2 -20 0 10 0 0 tril(A,-1) ans =0 0 01 0 02 2 0 b=diag(

30、A); b ans =1 1 1,4.3 迭代解法的几种形式 5.3.1 Jacobi迭代法,方程组 Ax=b A可写成 A=D-L-U其中:D=diaga11,a22,ann, -L、-U分别为A的严格下、上三角部分(不包括对角线元素).由 Ax=b x=Bx+f由此可构造迭代法:x(k+1)=Bx(k)+f其中:B=D-1(L+U)=I-D-1A, f=D-1b.,function y=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; while norm(y

31、-x0)=1.0e-6x0=y;y=B*x0+f;n=n+1; end n,例:用Jacobi方法求解, 设x(0)=0,精度为10-6。 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; jacobi(a,b,0;0;0) n =11 ans =0.99580.95790.7916,4.3.2 Gauss-Seidel迭代法,由原方程构造迭代方程x(k+1)=G x(k)+f其中:G=(D-L)-1 U, f=(D-L)-1 bD=diaga11,a22,ann,-L、-U分别为A的严格下、上三角部分(不包括对角线元素).,function y=seidel(

32、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; while norm(y-x0)=1.0e-6x0=y;y=G*x0+f;n=n+1; end n,例:对上例用Gauss-Seidel迭代法求解 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; seidel(a,b,0;0;0) n =7 ans =0.99580.95790.7916 例:分别用Jacobi和G-S 法迭代求解,看是否收敛。, a=1 2 -2; 1 1 1; 2 2 1

33、; b=9; 7; 6; jacobi(a,b,0;0;0) n =4 ans =-27268 seidel(a,b,0;0;0) n =1011 ans =1.0e+305 *-InfInf-1.7556,4.3.3 SOR迭代法,在很多情况下,J法和G-S法收敛较慢,所以考虑对G-S法进行改进。于是引入一种新的迭代法逐次超松弛迭代法(Succesise Over-Relaxation),记为SQR法。迭代公式为:X(k+1)= (D-wL)-1(1-w)D+wU)x(k) + w(D-wL)-1 b其中:w最佳值在1, 2)之间,不易计算得到,因此 w通常有经验给出。,function y

34、=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; while norm(y-x0)=1.0e-6x0=y;y=M*x0+f;n=n+1; end n,例:上例中,当w=1.103时,用SOR法求解原方程。 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; sor(a,b,1.103,0;0;0) n =8 ans =0.99580.95790.7916,4.3.4 两步迭代法,当线性方程系数矩阵

35、为对称正定时,可用一种特殊的迭代法来解决,其迭代公式为:(D-L)x(k+1/2) =U x(k) +b(D-U)x(k+1)=Lx(k+1/2) +b =x(k+1/2) =(D-L)-1 U x(k) + (D-L)-1 bx(k+1)= (D-U)-1 Lx(k+1/2) + (D-U)-1 b,function y=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; while n

36、orm(y-x0)=1.0e-6x0=y;y=G1*x0+f1; y=G2*y+f2;n=n+1; end n,例:求解方程组 a=10 -1 2 0; -1 11 -1 3; 2 -1 10 3; 0 3 -1 8; b=6; 25; -11; 15; twostp(a, b, 0; 0; 0; 0) n =7 ans =1.07911.9824-1.40440.9560,4.4 线性方程组的符号解法,在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数,如linsolve(A,b)等同于 X = sym(A)sym(b).solve(eqn1,eqn2,.,eqn

37、N,var1,var2,.,varN ),例: A=sym(10,-1,0;-1,10,-2;0,-2,10); b=(9; 7; 6); linsolve(A,b) ans = 473/475 91/95 376/475 vpa(ans) ans = .99578947368421052631578947368421 .95789473684210526315789473684211 .79157894736842105263157894736842,例: x,y = solve(x2 + x*y + y = 3,x2 - 4*x + 3 = 0,x,y) x = 1 3y = 1 -3/2

38、,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,例: 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) 1 full(a2) ans =0 0 0 0 01 0 0 0 00 1 0 0 00 0 1 0 00 0 0 1 0,

39、例: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+a2 a =(1,1) 4(2,1) 1(1,2) 1(2,2) 4(3,2) 1(2,3) 1(3,3) 4(4,3) 1,(3,4) 1(4,4) 4(5,4) 1(4,5) 1(5,5) 4 full(a) ans =4 1 0 0 01 4 1 0 00 1 4 1 00 0 1 4 10 0 0 1 4,格式 A=spdiags(B,d,m,n)生

40、成一mxn阶的稀疏矩阵,使得B的列放在由d指定的位置。 例: n=5 b=spdiags(ones(n,1),4*ones(n,1),ones(n,1), -1,0,1,n,n); full(b) ans =4 1 0 0 01 4 1 0 00 1 4 1 00 0 1 4 10 0 0 1 4,格式: spconvert(dd)对于无规律的稀疏矩阵,可使用此命令由外部数据转化为稀疏矩阵。 调用形式为:先用load函数加载以行表示对应位置和元素值的.dat文本文件,再用此命令转化为稀疏矩阵。 例:无规律稀疏矩阵的建立。 首先编制文本文件sp.dat如下: 5 1 5.00 3 5 8.00

41、4 4 2.00 5 5 0, load sp.dat spconvert(sp) ans =(5,1) 5(4,4) 2(3,5) 8 full(ans) ans =0 0 0 0 00 0 0 0 00 0 0 0 80 0 0 2 05 0 0 0 0,稀疏矩阵的计算:同满矩阵比较,稀疏矩阵在算法上有很大的不同。具体表现在存储空间减少,计算时间减少。例:比较求解下面方程组n1000时两种方法的差别。, n=1000; 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+a2;

42、b=ones(1000,1); tic; x=ab; t1=toc t1 =0.4800 a=full(a); tic; x=ab; t2=toc t2 =1.3220,4.6 矩阵的特征值问题 4.6.1一般矩阵的特征值与特征向量,格式: d=eig (A)只求解特征值。格式: V, D=eig (A)求解特征值和特征向量。,例:直接求解: A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; eig(A) ans =34.00008.9443-8.94430.0000,精确解: eig(sym(A) ans = 0 34 4*5(1/2) -4*5(1/

43、2) 高精度数值解: vpa(ans,70) ans = 0 34. 8.944271909999158785636694674925104941762473438446102897083588981642084 -8.94427190999915878563669467492510494176247343844610 28 97083588981642084,同时求出特征值与特征向量: 直接求解: v, d = eig(A) v =-0.5000 -0.8236 0.3764 -0.2236-0.5000 0.4236 0.0236 -0.6708-0.5000 0.0236 0.4236

44、0.6708-0.5000 0.3764 -0.8236 0.2236 d =34.0000 0 0 00 8.9443 0 00 0 -8.9443 00 0 0 0.0000,解析解: v,d=eig(sym(A) v = -1, 1, -8*5(1/2)-17, 8*5(1/2)-17 -3, 1, 4*5(1/2)+9, -4*5(1/2)+9 3, 1, 1, 1 1, 1, 4*5(1/2)+7, -4*5(1/2)+7 d = 0, 0, 0, 0 0, 34, 0, 0 0, 0, 4*5(1/2), 0 0, 0, 0, -4*5(1/2),4.6.2 矩阵的广义特征向量问题

45、,若B=I,则化成普通矩阵特征值问题。格式: d=eig (A,B)求解广义特征值。格式: V, D=eig (A,B)求解广义特征值和特征向量。,例:直接求解: A=5,7,6,5; 7,10,8,7; 6,8,10,9; 5,7,9,10; B=2,6,-1,-2; 5,-1,2,3; -3,-4,1,10; 5,-2,-3,8; V, D = eig(A, B) V =0.3697 -0.3741 + 0.6259i -0.3741 - 0.6259i 1.0000 0.9948 -0.0674 - 0.2531i -0.0674 + 0.2531i -0.6090 0.7979 0.9239 + 0.0264i 0.9239 - 0.0264i -0.2316 1.0000 -0.6599 - 0.3263i -0.6599 + 0.3263i 0.1319,

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 教学课件 > 大学教育

copyright@ 2008-2019 麦多课文库(www.mydoc123.com)网站版权所有
备案/许可证编号:苏ICP备17064731号-1