《matlab在科学计算中的应用4(姜志鹏).ppt》由会员分享，可在线阅读，更多相关《matlab在科学计算中的应用4(姜志鹏).ppt（98页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第四章第四章 线性代数性代数问题求解求解矩矩阵线性方程性方程组的直接解法的直接解法线性方程性方程组的迭代法的迭代法线性方程性方程组的符号解法的符号解法稀疏矩稀疏矩阵技技术特征特征值与特征向量与特征向量14.1 矩矩阵特殊矩特殊矩阵的的输入入数数值矩矩阵的的输入入零矩零矩阵、幺矩、幺矩阵及及单位矩位矩阵 生成生成n n n n方方阵：A=zeros(n),B=ones(n),C=eye(n)A=zeros(n),B=ones(n),C=eye(n)生成生成m m n n矩矩阵：A=zeros(m,n),B=ones(m,n),A=zeros(m,n),B=ones(m,n),C=eye(m,n)
2、C=eye(m,n)生成和矩生成和矩阵B B同同样位数的矩位数的矩阵：A=zeros(size(B)A=zeros(size(B)2随机元素矩随机元素矩阵若矩若矩阵随机元素随机元素满足足0,10,1区区间上的均匀分布上的均匀分布 生成生成n n m m阶标准均匀分布准均匀分布伪随机数矩随机数矩阵：A=rand(n,m)A=rand(n,m)生成生成n n n n阶标准均匀分布准均匀分布伪随机数方随机数方阵：A=rand(n)A=rand(n)3对角元素矩角元素矩阵 已知向量生成已知向量生成对角矩角矩阵：A=diag(V)A=diag(V)已知矩已知矩阵提取提取对角元素列向量：角元素列向量：V
3、Vdiag(A)diag(A)生成主生成主对角角线上第上第k k条条对角角线为V V的矩的矩阵：A=diag(V,k)A=diag(V,k)4例：例：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 05生成三生成三对角矩角矩阵
4、: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 46HilbertHilbert矩矩阵及逆及逆HilbertHilbert矩矩阵 生成生成n n阶的的HilbertHilbert矩矩阵：A=hilb(n)A=hilb(n)求取逆求取逆HilbertHilbert矩矩阵：B=invhilb(n)B=invhilb(n)7Hankel(Hankel(汉克克 )矩矩阵 其中：第一列的各个元素定其中：第一列的各个元素定义为C C向量，最后一行各向量，最后一行各个元素定个元素定义为R R。H H为对
5、称称阵。H1=hankel(C)由由 Hankel 矩矩阵反反对角角线上元素相等得出一下三角上元素相等得出一下三角阵均均为零的零的Hankel 矩矩阵8VandermondeVandermonde(范德蒙范德蒙)矩矩阵 9伴随矩伴随矩阵其中：其中：P(s)P(s)为首首项系数系数为1 1的多的多项式。式。10 例：考例：考虑一个多一个多项式式2*x4+4*x2+5*x+6，试写出写出该多多项式的伴随矩式的伴随矩阵。P=2 0 4 5 6;A=compan(P)A=0 -2.0000 -2.5000 -3.0000 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 011
6、符号矩符号矩阵的的输入入 数数值矩矩阵A A转换成符号矩成符号矩阵：B=sym(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 4.1.2 矩矩阵基本概念与性基本概念与性质行列式行列式 格式格式 ：d=det(A)d=det(A)例：求行列式例：求行列式 A=16 2 3 13;5 11 10 8;9 7 6 12;4 14 15 1;det(A)ans=013例：例
7、：tic,A=sym(hilb(20);det(A),toc ans=elapsed_time=2.3140高高阶的的Hilbert矩矩阵是接近奇异的矩是接近奇异的矩阵。14矩矩阵的迹的迹 格式：格式：t=trace(A)t=trace(A)矩矩阵的秩的秩格式：格式：r=rank(A)r=rank(A)用默用默认的精度求数的精度求数值秩秩 r=rank(Ar=rank(A，)给定精度下求数定精度下求数值秩秩 矩矩阵的秩也表示的秩也表示该矩矩阵中行列式不等于中行列式不等于0 0的子式的最大的子式的最大阶次。可次。可证行秩和列秩（行秩和列秩（线性无关的）性无关的）应相等。相等。15例例 A=16
8、2 3 13;5 11 10 8;9 7 6 12;4 14 15 1;rank(A)ans=3该矩矩阵的秩的秩为3，小于矩，小于矩阵的的阶次，故次，故为非非满秩矩秩矩阵。例例 用数用数值值方法和解析方法分方法和解析方法分别别求求2020的的Hilbert矩矩阵阵的秩的秩 H=hilb(20);rank(H)数数值方法方法ans=13 H=sym(hilb(20);rank(H)%解析方法，原矩解析方法，原矩阵为非奇异矩非奇异矩阵ans=2016矩矩阵范数范数17矩矩阵的范数定的范数定义：格式：格式：N=norm(A)求解默求解默认的的2范数范数 N=norm(A，选项)选项可可为1,2,in
9、f等等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()函数
10、，需先用函数，需先用double()函数函数转换成双精度数成双精度数值矩矩阵，再，再调用用norm()函数。函数。19特征多特征多项式式格式：格式：C=poly(A)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 A=sym(A);poly(A)运用符号工具箱运用符号工具箱 ans=x4-34*x3-80*x2+2720*x20矩矩阵多多项式的求解式的求解21符号多符号多项式与数式与数值多多项式的式的转换
11、格式：格式：f=poly2sym(P)f=poly2sym(P)或或 f=poly2sym(P,x)f=poly2sym(P,x)格式：格式：P=sym2poly(f)P=sym2poly(f)22例：例：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 623矩矩阵的逆矩的逆矩阵格式：格式：C=inv(A)C=inv(A)例：求例：求HilbertHilbert矩矩阵的逆矩的逆矩阵
12、 format long;H=hilb(4);H1=inv(H)H1=1.0e+003*24检验：H*H1ans=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.0000000000001
13、1计算算误差范数：差范数：norm(H*inv(H)-eye(size(H)ans=6.235798190375727e-013 H2=invhilb(4);norm(H*H2-eye(size(H)ans=5.684341886080802e-01425 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=H=hilb(13);H1=inv(H);norm(H*H1-eye(size(H)Warning:Matrix is close to
14、 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,-19
15、40400,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,3
16、7837800,-33297264,11099088 H=sym(hilb(30);norm(double(H*inv(H)-eye(size(H)ans=027例：奇异例：奇异阵求逆求逆 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*norm(A*B-eye(size(A)检验ans=1.6408151
17、3306419 A=sym(A);inv(A)奇异矩奇异矩阵不存在一个相不存在一个相应的逆矩的逆矩阵，用符号工具箱，用符号工具箱的函数也不行的函数也不行?Error using=sym/invError,(in inverse)singular matrix28同同样适用于含有适用于含有变量的矩量的矩阵求逆。求逆。例：例：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矩矩阵的相似的相似变换与正交矩与正交矩阵 其中
18、：其中：A A为一方一方阵，B B矩矩阵非奇异。非奇异。相似相似变换后，后，X X矩矩阵的秩、迹、行列式与特征的秩、迹、行列式与特征值等等均不均不发生生变化，其化，其值与与A 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
19、.3106 -0.2725 0.7406 norm(Q*Q-eye(4)ans=4.6395e-016 norm(Q*Q-eye(4)ans=4.9270e-01631 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 -0.0549 0.3658 0.1776 0.2882 det(A),det(C)ans=120ans=120.000032 trace(A),trace(C)ans=21ans=21.0000 rank(A),rank(C)an
20、s=4ans=433 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 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.5000 0.6708 0.5000 -0.5000 -0.2236 -0.5000 -0.5000 0.2236 -0.5000 -0.5000 -0.6708 0.5000
21、 norm(Q*Q-eye(3)%Q*Q=Ians=1.0140e-015354.2 线性方程性方程组直接解法直接解法线性方程性方程组直接求解矩直接求解矩阵除法除法关于关于线性方程性方程组的直接解法，如的直接解法，如Gauss消去法、消去法、选主元消去法、平方根法、追赶法等等，在主元消去法、平方根法、追赶法等等，在MATLAB中，只需用中，只需用“”或或“”就解决就解决问题。它内部它内部实际包含着包含着许许多多的自适多多的自适应算法，算法，如如对超定方程用最小二乘法，超定方程用最小二乘法，对欠定方程欠定方程时它将它将给出范数最小的一个解，解三出范数最小的一个解，解三对角角阵方方程程组时用追赶法
22、等等。用追赶法等等。格式：格式：x=Ab36例：解方程例：解方程组 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.446737线性方程性方程组直接求解判定求解直接求解判定求解3839例：例：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
23、),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.733340检验 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=041原方程原方程组对应的的齐次方程次方程组的解的解求取求取A矩矩阵的化零矩的化零矩阵：格式：格式：Z=null(A)求取求取A矩矩阵的化零矩的化零矩
24、阵的的规范形式：范形式：格式：格式：Z=null(A,r)42null是用来求齐次线性方程组是用来求齐次线性方程组的基础解系的，加上的基础解系的，加上r则求则求出的是一组出的是一组 最小正整数解，最小正整数解，如果不加，则求出的是解空如果不加，则求出的是解空间的规范正交基间的规范正交基。例：例：判断可解性判断可解性 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
25、 0 0 1.000043 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-01544解析解解析解 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 45验证得出的解得出的解 a1=randn
26、(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/13146 摩摩尔彭彭罗斯广斯广义逆求解出的方程最小二乘解逆求解出的方程最小二乘解不不满足原始代数方程。足原始代数方程。474.2.3 线性方程性方程组的直接求解分析的直接求解分析LU分解分解 484950格式格式 l,u,p=lu(A)
27、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=A51例：例：对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
28、0 052 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 -2.5000 0 0 2.500053QR分解分解 将将矩矩阵A分分解解成成一一个个正正交交矩矩阵与与一一个个上上三三角角矩矩阵的乘的乘积。求求得得正正交交矩矩阵Q和和上上三三角角阵R，Q和和R满足足A=QR。54格式：格式：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.8331 0.5456 -0.0478
29、-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 -1.0413 -2.0826 0 0 -0.0000 0 0 055Cholesky(乔里斯基里斯基 )分解分解 若矩若矩阵A为 n阶对称正定称正定阵,则存在唯存在唯一的一的对角元素角元素为正的三角正的三角阵D，使得，使得 5657格式：D=chol(A)58例：进行例：进行Cholesky分解。分解。A=16 4 8;4 5-4;8-4 22;D=c
30、hol(A)D=4 1 2 0 2 -3 0 0 3利用矩利用矩阵的的LU、QR和和cholesky分解求方程分解求方程组的解的解（1）LU分解：分解：A*X=b A*X=b 变成成 L*U*X=bL*U*X=b所以所以 X=U(Lb)X=U(Lb)这样可以大大提高运算速度。可以大大提高运算速度。例：例：求方程求方程组 的一个特解。的一个特解。解：解：A=4 2-1;3-1 2;11 3 0;A=4 2-1;3-1 2;11 3 0;B=2 10 8;B=2 10 8;D=det(A)D=det(A)D=D=0 059 L,U=lu(A)L=0.3636 -0.5000 1.0000 0.27
31、27 1.0000 0 1.0000 0 0U=11.0000 3.0000 0 0 -1.8182 2.0000 0 0 0.000060 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*X验证其正确性。其正确性。1.4862 1.3511 A*X%Matlab7.0显示没有解示没有解ans=0 8 8
32、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)这三三种种分分解解，在在求求解解大大型型方方程程组时很很有有用用。其其优点点是运算速度快、可以是运算速度快、可以节省磁省磁盘空空
33、间、节省内存。省内存。62三个三个变换 在在线性方程性方程组的迭代求解中，要用到系数的迭代求解中，要用到系数矩矩阵A的上三角矩的上三角矩阵、对角角阵和下三角矩和下三角矩阵。此三个此三个变换在在MATLAB中可由以下函数中可由以下函数实现。上三角上三角变换：格式格式 triu(A,1)对角角变换：格式格式 diag(A)下三角下三角变换：格式格式 tril(A,-1)例：例：对此矩此矩阵做三种做三种变换。63 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);
34、bans=1 1 1644.4 线性方程性方程组的符号解法的符号解法 在在MATLAB的的Symbolic Toolbox中提供了中提供了线性方程的符号求解函数，如性方程的符号求解函数，如 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/475 91/95 376/475 vpa(ans)ans=66例：例：x,y=solve(x2+x*y+y=3,x2-4
35、*x+3=0,x,y)x=1 3 y=1 -3/2 674.3 迭代解法的几种形式迭代解法的几种形式4.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=
36、-tril(a,-1);B=D(L+U);f=Db;y=B*x0+f;n=1;while norm(y-x0)=1.0e-6 x0=y;y=B*x0+f;n=n+1;endn例：用例：用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=11ans=0.9958 0.9579 0.79164.3.2 Gauss-Seidel迭代法迭代法由原方程构造迭代方程由原方程构造迭代方程 x(k+1)=G x(k)+f其中：其中：G=(D-L)-1 U,f=(D-L)-1 b D=diag
37、a11,a22,ann,-L、-U分分别为A的的严格下、上三角部格下、上三角部分（不包括分（不包括对角角线元素）元素）.function y=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;while norm(y-x0)=1.0e-6 x0=y;y=G*x0+f;n=n+1;endn例：例：对上例用上例用Gauss-Seidel迭代法求解迭代法求解 a=10-1 0;-1 10-2;0-2 10;b=9;7;6;seidel(a,b,0;0;0)n=7ans=0.9
38、958 0.9579 0.7916例：分例：分别用用Jacobi和和G-S法迭代求解法迭代求解,看是否收看是否收敛。a=1 2-2;1 1 1;2 2 1;b=9;7;6;jacobi(a,b,0;0;0)n=4ans=-27 26 8 seidel(a,b,0;0;0)n=1011ans=1.0e+305*-Inf Inf -1.75564.3.3 SOR迭代法迭代法 在很多情况下，在很多情况下，J法和法和G-S法收法收敛较慢，慢，所以考所以考虑对G-S法法进行改行改进。于是引入一。于是引入一种新的迭代法逐次超松弛迭代法种新的迭代法逐次超松弛迭代法（Succesise Over-Relaxa
39、tion),记为SQR法。法。迭代公式迭代公式为：X(k+1)=(D-wL)-1(1-w)D+wU)x(k)+w(D-wL)-1 b 其中：其中：w最佳最佳值在在1,2）之）之间，不易，不易计算得到，因此算得到，因此 w通常有通常有经验给出。出。function y=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-6 x0=y;y=M*x0+f;n=n+1;endn例：上例中，当例：上例中，当w
40、=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=8ans=0.9958 0.9579 0.79164.3.4 两步迭代法两步迭代法 当当线性方程系数矩性方程系数矩阵为对称正定称正定时，可用一种特殊的迭代法来解决，其迭代可用一种特殊的迭代法来解决，其迭代公式公式为：（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 b x(k+1)=(D-U)-1 Lx(k+1/2)+(D-U
41、)-1 bfunction 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 norm(y-x0)=1.0e-6 x0=y;y=G1*x0+f1;y=G2*y+f2;n=n+1;endn例：求解方程例：求解方程组 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=7ans=1.0791 1.98
42、24 -1.4044 0.95604.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)481例：例：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 0 1 0
43、 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 082例：例：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)4 full(a)ans=4 1 0 0 0 1 4 1 0 0
44、0 1 4 1 0 0 0 1 4 1 0 0 0 1 484格式格式 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 485格式：格式：spconvert(dd)对于无于无规律的稀疏矩律的稀疏矩阵，可使用此命令，可使用此命令由外部数据由外部数据转化化为稀疏矩稀疏矩阵。调
45、用形式用形式为：先用：先用load函数加函数加载以行表示以行表示对应位置和元素位置和元素值的的.dat文本文件，再用此命文本文件，再用此命令令转化化为稀疏矩稀疏矩阵。例：无例：无规律稀疏矩律稀疏矩阵的建立。的建立。首先首先编制文本文件制文本文件sp.dat如下：如下：5 1 5.003 5 8.004 4 2.005 5 086 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 08788稀疏矩阵的计算稀疏矩阵的计算：同满矩
46、阵比较，稀疏矩阵在算法上有很大的不同满矩阵比较，稀疏矩阵在算法上有很大的不同。具体表现在存储空间减少，计算时间减少。同。具体表现在存储空间减少，计算时间减少。例：比较求解下面方程组例：比较求解下面方程组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;b=ones(1000,1);tic;x=ab;t1=toct1=0.4800 a=full(a);tic;x=ab;t2=toct2=1.3220894.6 4.6 矩矩
47、阵的特征的特征值问题 一般矩一般矩阵的特征的特征值与特征向量与特征向量格式：格式：d=eig(A)只求解特征只求解特征值。格式：格式：V,D=eig(A)求解特征求解特征值和特征向量。和特征向量。90例：例：直接求解：直接求解：A=16 2 3 13;5 11 10 8;9 7 6 12;4 14 15 1;eig(A)ans=34.0000 8.9443 -8.9443 0.000091精确解：精确解：eig(sym(A)ans=0 34 4*5(1/2)-4*5(1/2)高精度数高精度数值解：解：vpa(ans,70)ans=0 34.28 9708358898164208492同同时求出
48、特征求出特征值与特征向量：与特征向量：直接求解：直接求解：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)93解析解：解析解：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
49、,-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)944.6.2 4.6.2 矩矩阵的广的广义特征向量特征向量问题 若若B=I，则化成普通矩化成普通矩阵特征特征值问题。格式：格式：d=eig(A，B)求解广求解广义特征特征值。格式：格式：V,D=eig(A，B)求解广求解广义特征特征值和特征向量。和特征向量。95例：例：直接求解：直接求解：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
50、;-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 96D=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