《最新在科学计算中的应用4幻灯片》由会员分享，可在线阅读，更多相关《最新在科学计算中的应用4幻灯片（87页珍藏版）》请在金锄头文库上搜索。
1、在科学计算中的应用在科学计算中的应用4 44.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) 2Vandermonde(范德蒙)矩阵 9伴随矩阵其中：P(s)为首项系数为1的多项式。 10 例：考虑一个多项式2*x4+4*x2+5*x+6，试写出该多项式的伴随矩阵。 P=2 0 4 5 6;A=compan(P)A = 0 -2.0000
2、-2.5000 -3.0000 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 011符号矩阵的输入 数值矩阵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/5124.1.2 矩阵基本概念与性质行列式 格式 ：d=det(A)例：求行列式 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1;
3、 det(A)ans = 013例： tic, A=sym(hilb(20); det(A), toc ans = 1/2377454716768534509091644243427616440175419837753486493033185331234419759310644585187585766816573773440565759867265558971765638419710793303386582324149811241023554489166154717809635257797836800000000000000000000000000000000000elapsed_time =
4、 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) % 解析
5、方法，原矩阵为非奇异矩阵ans = 2016矩阵范数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; nor
6、m(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.819
7、840539024018e-012 A=sym(A); poly(A) 运用符号工具箱 ans = x4-34*x3-80*x2+2720*x20矩阵多项式的求解返回返回n次多项式在次多项式在x处的值处的值, x可以是一个矩阵矩阵或者一个向量,在这两种情况下，polyval计算在x中任意元素处的多项式a的估值. 21X=1 1;2 2p=poly(X)C1=poly(p,X)C2=polyvalm(p,X)C1= -2 -2 -2 -2C2= 0 0 0 022符号多项式与数值多项式的转换格式： f=poly2sym(P) 或 f=poly2sym(P,x) 格式： P=sym2poly(f)
8、23例： 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 624矩阵的逆矩阵格式： 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.199999999999
9、90 -2.69999999999976 1.67999999999984 0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961 -0.13999999999999 1.67999999999984 -4.19999999999961 2.7999999999997425检验： H*H1ans = 1.00000000000001 0.00000000000023 -0.00000000000045 0.00000000000023 0.00000000000001 1.00000000000011 -0.0
10、0000000000011 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-01426 H=hilb(10); H1=inv(H); norm(H*H1-e
11、ye(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(
12、size(H)ans = 11.37062973181391对接近于奇异矩阵，高阶一般不建议用inv( )，可用符号工具箱。27 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, -72
13、765000, 62092800, -2018016048510, -1940400, 18711000, -72765000, 133402500, -115259760, 37837800-38808, 1596672, -15717240, 62092800, -115259760, 100590336, -3329726412012, -504504, 5045040, -20180160, 37837800, -33297264, 11099088 H=sym(hilb(30); norm(double(H*inv(H)-eye(size(H)ans = 028例：奇异阵求逆 A=1
14、6 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.4442493013
15、1968 -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/invError, (in inverse) sin
16、gular matrix29同样适用于含有变量的矩阵求逆。例： 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)30矩阵的相似变换与正交矩阵 其中：A为一方阵，B矩阵非奇异。 相似变换后，X矩阵的秩、迹、行列式与特征值等均不发生变化，其值与A矩阵完全一致。 对于一类特殊的相似变换满足如下条件，称为正交基矩阵。31例： A=5,9,8,3; 0,3,2,4; 2,3,5,9; 3,4,5,8; Q=ort
17、h(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-01632 C=Q*A*QC = 17.9251 6.4627 -4.4714 -2.0354 -0.0282 1.7194 4.6816 -5.0735 0.6800 -0.9386 1.0674 0.6631
18、 -0.0549 0.3658 0.1776 0.2882 det(A),det(C)ans = 120ans = 120.000033 trace(A),trace(C)ans = 21ans = 21.0000 rank(A),rank(C)ans = 4ans = 434 eig(A), eig(C)ans = 17.8205 1.1908 + 2.6499i 1.1908 - 2.6499i 0.7979 ans = 17.8205 1.1908 + 2.6499i 1.1908 - 2.6499i 0.7979 35例： A=16,2,3,13; 5,11,10,8; 9,7,6,1
19、2; 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) % Q*Q=Ians = 1.0140e-015364.2 线性方程组直接解法4.2.1线性方程组直接求解矩阵除法关于线性方程组的直接解法，如Gauss消去法、选主元消去法、平方根法、追赶法等等，在MATLAB中，只需用“”或“”就解决问题。它内部实际包含着许许多多的自适应算法，如对超定方程
20、用最小二乘法，对欠定方程时它将给出范数最小的一个解，解三对角阵方程组时用追赶法等等。 格式： x=Ab37例：解方程组 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; xans = -0.1819 -1.6630 2.2172 -0.4467384.2.2线性方程组直接求解判定求解3940例： A=1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2; B
21、=5 1; 4 2; 3 3; 2 4; C=A B; rank(A), rank(C)ans = 4ans = 4 x=inv(A)*B % AB x = -1.8000 2.4000 1.8667 -1.2667 3.8667 -3.2667 -2.1333 2.733341检验 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 = 042原方程组对应的齐次方程组的解求取A矩阵的化
22、零矩阵： 格式： Z=null(A)求取A矩阵的化零矩阵的规范形式： 格式： Z=null(A, r )null是用来求齐次线性方程组的基础解系的，加上r则求出的是一组 最小正整数解，如果不加，则求出的是解空间的规范正交基。 43例：判断可解性 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.000044 x0=pinv(A)*B % 得出
23、一个特解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-01545解析解 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 46验证得出的解 a1=randn(1); a2=rand(1); % 取不同分布的随机数 x=a1*Z(:,1)
24、+a2*Z(:,2)+x0; norm(double(A*x-B)ans = 0通解 syms a1 a2; x=a1*Z(:,1)+a2*Z(:,2)+x0x = 2*a1+3*a2+125/131 -5/2*a1-7/2*a2+96/131 a1-10/131 a2-39/13147 摩尔彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。484.2.3 线性方程组的直接求解分析LU分解 495051格式 l,u,p=lu(A) L是一个单位下三角矩阵，u是一个上三角矩阵， p是代表选主元的置换矩阵。故：Ax=y = PAx=Py = LUx=Py = PA=LU l,u=lu(A)其中
25、l等于P-1 L，u等于U，所以(P-1 L)U=A52例：对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.0000u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000p = 0 0 1 0 1 0 1 0 053 l,u=lu(A) lP-1 Ll = 0.2500 0.5000 1.0000 0.5000 1.0000 0 1.0000 0 0u = 4.0000 6.0000 7.0000 0 1.0000
26、 -2.5000 0 0 2.500054QR分解 将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。 求得正交矩阵Q和上三角阵R，Q和R满足A=QR。 格式： Q,R = qr(A)55例： A = 1 2 3;4 5 6; 7 8 9; 10 11 12; Q,R = qr(A)Q = -0.0776 -0.8331 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.3748R = -12.8841 -14.5916 -16.2992 0
27、 -1.0413 -2.0826 0 0 -0.0000 0 0 056Cholesky(乔里斯基 )分解 若矩阵A为 n阶对称正定阵,则存在唯一的对角元素为正的三角阵D，使得 57格式： D=chol(A)58例：进行Cholesky分解。 A=16 4 8; 4 5 -4; 8 -4 22; D=chol(A)D = 4 1 2 0 2 -3 0 0 359利用矩阵的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
28、=2 10 8; D=det(A)D = 060 L,U=lu(A)L = 0.3636 -0.5000 1.0000 0.2727 1.0000 0 1.0000 0 0U = 11.0000 3.0000 0 0 -1.8182 2.0000 0 0 0.000061 X=U(LB)Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.018587e-017.X = 1.0e+016 * % 结果中的警告是由于系数行列式为零产生的。 -0.4053 % 可以通过A
29、*X验证其正确性。 1.4862 1.3511 A*X % Matlab7.0显示没有解ans = 0 8 862（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) 这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。63三个变换 在线性方程组的迭代求解中，要用到系数矩阵
30、A的上三角矩阵、对角阵和下三角矩阵。此三个变换在MATLAB中可由以下函数实现。上三角变换： 格式 triu(A,1)对角变换： 格式 diag(A)下三角变换： 格式 tril(A,-1)例：对此矩阵做三种变换。64 A=1 2 -2;1 1 1;2 2 1; triu(A,1)ans = 0 2 -2 0 0 1 0 0 0 tril(A,-1)ans = 0 0 0 1 0 0 2 2 0 b=diag(A); bans = 1 1 1654.4 线性方程组的符号解法 在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数，如 linsolve(A,b) 等同于
31、X = sym(A)sym(b). solve(eqn1,eqn2,.,eqnN,var1,var2,.,varN )66例： A=sym(10,-1,0;-1,10,-2;0,-2,10)； b=(9;7;6)； linsolve(A,b) ? 仅支持数值型的，用支持数值型的，用doubleans = 473/475 91/95 376/475 vpa(ans)ans = .99578947368421052631578947368421 .95789473684210526315789473684211 .7915789473684210526315789473684267例： x,y =
32、 solve(x2 + x*y + y = 3,x2 - 4*x + 3 = 0,x,y) x = 1 3 y = 1 -3/2 684.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) 469例： a2=sparse(2:n, 1:n-1,ones(1,n-1),n,n)a2 = (2,1) 1 (3,2) 1 (4,3)
33、1 (5,4) 1 full(a2)ans = 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 070例：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 71 (3,4) 1 (4,4) 4 (5,4) 1 (4,5) 1
34、(5,5) 4 full(a)ans = 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 472格式 A=spdiags(B,d,m,n) 生成一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 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 473格式： spconvert(dd) 对于无规律的稀疏矩阵，可使用此命令由外部数据转化为稀疏矩阵。调用形式
35、为：先用load函数加载以行表示对应位置和元素值的.dat文本文件，再用此命令转化为稀疏矩阵。例：无规律稀疏矩阵的建立。首先编制文本文件sp.dat如下：5 1 5.003 5 8.004 4 2.005 5 074 load sp.dat spconvert(sp)ans = (5,1) 5 (4,4) 2 (3,5) 8 full(ans)ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 0 2 0 5 0 0 0 075稀疏矩阵的计算： 同满矩阵比较，稀疏矩阵在算法上有很大的不同。具体表现在存储空间减少，计算时间减少。例：比较求解下面方程组n1000时两种方
36、法的差别。76 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; b=ones(1000,1); tic; x=ab; t1=toct1 = 0.4800 a=full(a); tic; x=ab; t2=toct2 = 1.3220774.6 矩阵的特征值问题 4.6.1一般矩阵的特征值与特征向量格式： d=eig (A) 只求解特征值。格式： V, D=eig (A) 求解特征值和特征向量。78例：直接求解： A=16 2 3 13; 5 11 10 8
37、; 9 7 6 12; 4 14 15 1; eig(A)ans = 34.0000 8.9443 -8.9443 0.000079精确解： eig(sym(A)ans = 0 34 4*5(1/2) -4*5(1/2)高精度数值解： vpa(ans,70)ans = 0 34. 8.944271909999158785636694674925104941762473438446102897083588981642084 -8.9442719099991587856366946749251049417624734384461028 9708358898164208480同时求出特征值与特征向量
38、：直接求解： 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 0.6708 -0.5000 0.3764 -0.8236 0.2236d = 34.0000 0 0 0 0 8.9443 0 0 0 0 -8.9443 0 0 0 0 0.0000 norm(A*v-v*d)81解析解： 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)
39、+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) 824.6.2 矩阵的广义特征向量问题 若B=I，则化成普通矩阵特征值问题。格式： d=eig (A，B) 求解广义特征值。格式： V, D=eig (A，B) 求解广义特征值和特征向量。83例：直接求解： 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
40、, 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 84D = 4.7564 0 0 0 0 0.0471 + 0.1750i 0 0 0 0 0.0471 - 0.1750i 0 0 0 0 -0.0037 检验： norm(A*V-B*V*D)ans = 1.3897e-014 符号运算工具箱中的eig( )函数不支持广义特征值的运算。85E N D86结束语结束语谢谢大家聆听！谢谢大家聆听！87
