1 / 94

第四章 线性代数问题求解

第四章 线性代数问题求解. 矩阵 线性方程组的直接解法 线性方程组的迭代法 线性方程组的符号解法 稀疏矩阵技术 特征值与特征向量. 4.1 矩阵 4.1.1 特殊矩阵的输入. 数值矩阵的输入 零矩阵、幺矩阵及单位矩阵 生成 n  n 方阵： A=zeros(n), B=ones(n), C=eye(n) 生成 m  n 矩阵： A=zeros(m,n), B=ones(m,n), C=eye(m,n) 生成和矩阵 B 同样位数的矩阵： A=zeros(size(B)). 随机元素矩阵 若矩阵随机元素满足 [0,1] 区间上的均匀分布

gitano
Download Presentation

第四章 线性代数问题求解

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 第四章 线性代数问题求解 • 矩阵 • 线性方程组的直接解法 • 线性方程组的迭代法 • 线性方程组的符号解法 • 稀疏矩阵技术 • 特征值与特征向量

  2. 4.1 矩阵4.1.1特殊矩阵的输入 • 数值矩阵的输入 • 零矩阵、幺矩阵及单位矩阵 生成nn方阵： A=zeros(n), B=ones(n), C=eye(n) 生成mn矩阵： A=zeros(m,n), B=ones(m,n), C=eye(m,n) 生成和矩阵B同样位数的矩阵： A=zeros(size(B))

  3. 随机元素矩阵 若矩阵随机元素满足[0,1]区间上的均匀分布 生成nm阶标准均匀分布为随机数矩阵： A=rand(n,m) 生成nn阶标准均匀分布为随机数方阵： A=rand(n)

  4. 对角元素矩阵 已知向量生成对角矩阵： A=diag(V) 已知矩阵提取对角元素列向量： V＝diag(A) 生成主对角线上第k条对角线为V的矩阵： A=diag(V,k)

  5. 例：diag( )函数的不同调用格式 >> C=[1 2 3]; V=diag(C) % 生成对角矩阵 V = 1 0 0 0 2 0 0 0 3 >> V1=diag(V)' % 将列向量通过转置变换成行向量 V1 = 1 2 3 >> C=[1 2 3]; V=diag(C,2) % 主对角线上第 k条对角线为C的矩阵 V = 0 0 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0

  6. 生成三对角矩阵: >> V=diag([1 2 3 4])+diag([2 3 4],1)+diag([5 4 3],-1) V = 1 2 0 0 5 2 3 0 0 4 3 4 0 0 3 4

  7. Hilbert矩阵及逆Hilbert矩阵 生成n阶的Hilbert矩阵： A=hilb(n) 求取逆Hilbert矩阵： B=invhilb(n)

  8. Hankel(汉克 ) 矩阵 其中：第一列的各个元素定义为C向量，最后一行各个元素定义为R。H为对称阵。 H1=hankel(C) 由 Hankel矩阵反对角线上元素相等得出一下三角阵均为零的Hankel矩阵

  9. Vandermonde(范德蒙)矩阵

  10. 伴随矩阵 其中：P(s)为首项系数为一的多向式。

  11. 符号矩阵的输入 数值矩阵A转换成符号矩阵： B=sym(A) 例： >> A=hilb(3) A = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.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]

  12. 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

  13. 例： >> tic, A=sym(hilb(20)); det(A), toc ans = 1/2377454716768534509091644243427616440175419837753486493033185331234419759310644585187585766816573773440565759867265558971765638419710793303386582324149811241023554489166154717809635257797836800000000000000000000000000000000000 elapsed_time = 2.3140 高阶的Hilbert矩阵是接近奇异的矩阵。

  14. 矩阵的迹 格式： t=trace(A) • 矩阵的秩 格式：r=rank(A) ％用默认的精度求数值秩 r=rank(A， ) ％给定精度下求数值秩 矩阵的秩也表示该矩阵中行列式不等于0的子式的最大阶次。可证行秩和列秩（线性无关的）应相等。

  15. 例 >> 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

  16. 矩阵范数

  17. 矩阵的范数定义： 格式： N=norm(A)％求解默认的2范数 N=norm(A，选项)％选项可为1,2,inf等

  18. 例：求一向量、矩阵的范数 >> 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( )函数转换成双精度数值矩阵，再调用norm( )函数。

  19. 特征多项式 格式： 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 = x^4-34*x^3-80*x^2+2720*x

  20. 矩阵多项式的求解

  21. 符号多项式与数值多项式的转换 格式： f=poly2sym(P)或 f=poly2sym(P,x) 格式： P=sym2poly(f)

  22. 例： >> P=[1 2 3 4 5 6]; % 先由系数按降幂顺序排列表示多项式 >> f=poly2sym(P,'v') % 以 v 为算子表示多项式 f = v^5+2*v^4+3*v^3+4*v^2+5*v+6 >> P=sym2poly(f) P = 1 2 3 4 5 6

  23. 矩阵的逆矩阵 格式： C=inv(A) 例： >> format long; H=hilb(4); H1=inv(H) H1 = 1.0e+003 * 0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999 -0.11999999999999 1.19999999999990 -2.69999999999976 1.67999999999984 0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961 -0.13999999999999 1.67999999999984 -4.19999999999961 2.79999999999974

  24. 检验： >> H*H1 ans = 1.00000000000001 0.00000000000023 -0.00000000000045 0.00000000000023 0.00000000000001 1.00000000000011 -0.00000000000011 0.00000000000011 0.00000000000001 0 1.00000000000011 0 0.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

  25. >> 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.23696008570294 >> H2=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, -504504] [8820, -317520, 2857680, -10584000, 18711000, -15717240, 5045040] [-29400, 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(H*inv(H)-eye(size(H)))) ans = 0

  27. 例：奇异阵求逆 >> 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.93824992236885 2.81474976710656 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 using ==> sym/inv Error, (in inverse) singular matrix

  28. 同样适用于含有变量的矩阵求逆。 例： >> 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)]

  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.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

  31. 例： >> A=[16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1]; >> Q=orth(A) ％A为奇异矩阵，故得出的Q为长方形矩阵 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

  32. 4.2 线性方程组直接解法4.2.1线性方程组直接求解－矩阵除法 • 关于线性方程组的直接解法，如Gauss消去法、选主元消去法、平方根法、追赶法等等，在MATLAB中，只需用“／”或“\”就解决问题。它内部实际包含着许许多多的自适应算法，如对超定方程用最小二乘法，对欠定方程时它将给出范数最小的一个解，解三对角阵方程组时用追赶法等等。 格式： x=A\b

  33. 例：解方程组 >> 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=A\b; x' ans = -0.1819 -1.6630 2.2172 -0.4467

  34. 4.2.2线性方程组直接求解－判定求解

  35. 例： >> 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 = 4 >> x=inv(A)*B x = -1.8000 2.4000 1.8667 -1.2667 3.8667 -3.2667 -2.1333 2.7333

  36. 检验 >> norm(A*x-B) ans = 7.4738e-015 精确解 >> x1=inv(sym(A))*B x1 = [ -9/5, 12/5] [ 28/15, -19/15] [ 58/15, -49/15] [ -32/15, 41/15] 检验 >> norm(double(A*x1-B)) ans = 0

  37. 原方程组对应的齐次方程组的解 求取A矩阵的化零矩阵： 格式： Z=null(A) 求取A矩阵的化零矩阵的规范形式： 格式： Z=null(A, ‘ r ’)

  38. 例： 判断可解性 >> A=[1 2 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 0 0 1.0000

  39. >> x0=pinv(A)*B % 得出一个特解 x0 = 0.9542 0.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

  40. 解析解 >> 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]

  41. 验证得出的解 >> 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] [ a2-39/131]

  42. 摩尔－彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。摩尔－彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。

  43. 4.2.3 线性方程组的直接求解分析 • LU分解

  44. 格式 [l,u,p]=lu(A) L是一个单位下三角矩阵，u是一个上三角矩阵， p是代表选主元的置换矩阵。 故：Ax=y => PAx=Py => LUx=Py => PA=LU [l,u]=lu(A) 其中l等于P-1 L，u等于U，所以(P-1 L)U=A

  45. 例：对A进行LU分解 >> A=[1 2 3; 2 4 1; 4 6 7]; >> [l,u,p]=lu(A) l = 1.0000 0 0 0.5000 1.0000 0 0.2500 0.5000 1.0000 u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000 p = 0 0 1 0 1 0 1 0 0

  46. >> [l,u]=lu(A) ％ l＝P-1 L l = 0.2500 0.5000 1.0000 0.5000 1.0000 0 1.0000 0 0 u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000

  47. QR分解 将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。 求得正交矩阵Q和上三角阵R，Q和R满足A=QR。 格式： [Q,R] = qr(A)

More Related