线性代数问题求解_第1页
线性代数问题求解_第2页
线性代数问题求解_第3页
线性代数问题求解_第4页
线性代数问题求解_第5页
已阅读5页,还剩89页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

线性代数问题求解第1页,共94页,2023年,2月20日,星期二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))%size返回m、n两个值

第2页,共94页,2023年,2月20日,星期二随机元素矩阵若矩阵随机元素满足[0,1]区间上的均匀分布生成nm阶标准均匀分布伪随机数矩阵:

A=rand(n,m)

生成nn阶标准均匀分布伪随机数方阵:

A=rand(n)第3页,共94页,2023年,2月20日,星期二对角元素矩阵已知向量生成对角矩阵:

A=diag(V)

已知矩阵提取对角元素列向量:

V=diag(A)

生成主对角线上第k条对角线为V的矩阵:

A=diag(V,k)第4页,共94页,2023年,2月20日,星期二例:diag()函数的不同调用格式>>C=[123];V=diag(C)%生成对角矩阵V=100020003>>V1=diag(V)'%将列向量通过转置变换成行向量V1=123>>C=[123];V=diag(C,2)%主对角线上第k条对角线为C的矩阵V=0010000020000030000000000第5页,共94页,2023年,2月20日,星期二生成三对角矩阵:>>V=diag([1234])+diag([234],1)+diag([543],-1)V=1200523004340034第6页,共94页,2023年,2月20日,星期二Hilbert矩阵及逆Hilbert矩阵

生成n阶的Hilbert矩阵:

A=hilb(n)

求取逆Hilbert矩阵:

B=invhilb(n)第7页,共94页,2023年,2月20日,星期二Hankel(汉克)矩阵

其中:第一列的各个元素定义为C向量,最后一行各个元素定义为R。H为对称阵。

H1=hankel(C)

由Hankel

矩阵反对角线上元素相等得出一下三角阵均为零的Hankel

矩阵第8页,共94页,2023年,2月20日,星期二Vandermonde(范德蒙)矩阵

第9页,共94页,2023年,2月20日,星期二伴随矩阵其中:P(s)为首项系数为1的多项式。

第10页,共94页,2023年,2月20日,星期二

例:考虑一个多项式2*x^4+4*x^2+5*x+6,试写出该多项式的伴随矩阵。>>P=[20456];A=compan(P)A=0-2.0000-2.5000-3.00001.000000001.000000001.00000第11页,共94页,2023年,2月20日,星期二符号矩阵的输入数值矩阵A转换成符号矩阵:

B=sym(A)例:>>A=hilb(3)A=1.00000.50000.33330.50000.33330.25000.33330.25000.2000>>B=sym(A)B=[1,1/2,1/3][1/2,1/3,1/4][1/3,1/4,1/5]第12页,共94页,2023年,2月20日,星期二4.1.2矩阵基本概念与性质行列式格式:d=det(A)例:求行列式>>A=[162313;511108;97612;414151];det(A)ans=0第13页,共94页,2023年,2月20日,星期二例:>>tic,A=sym(hilb(20));det(A),tocans=1/2377454716768534509091644243427616440175419837753486493033185331234419759310644585187585766816573773440565759867265558971765638419710793303386582324149811241023554489166154717809635257797836800000000000000000000000000000000000elapsed_time=2.3140高阶的Hilbert矩阵是接近奇异的矩阵。第14页,共94页,2023年,2月20日,星期二矩阵的迹格式:t=trace(A)矩阵的秩格式:r=rank(A)%用默认的精度求数值秩

r=rank(A,)%给定精度下求数值秩矩阵的秩也表示该矩阵中行列式不等于0的子式的最大阶次。可证行秩和列秩(线性无关的)应相等。第15页,共94页,2023年,2月20日,星期二例>>A=[162313;511108;97612;414151];rank(A)ans=3该矩阵的秩为3,小于矩阵的阶次,故为非满秩矩阵。例>>H=hilb(20);rank(H)%数值方法ans=13det(A)=0则A为奇异矩阵>>H=sym(hilb(20));rank(H)%解析方法,原矩阵为非奇异矩阵ans=20第16页,共94页,2023年,2月20日,星期二矩阵范数第17页,共94页,2023年,2月20日,星期二矩阵的范数定义:格式:

N=norm(A)

%求解默认的2范数

N=norm(A,选项)

%选项可为1,2,inf等第18页,共94页,2023年,2月20日,星期二例:求一向量、矩阵的范数>>a=[162313];>>[norm(a),norm(a,2),norm(a,1),norm(a,Inf)]ans=2.092844953645635e+0012.092844953645635e+0013.400000000000000e+0011.600000000000000e+001>>A=[162313;511108;97612;414151];>>[norm(A),norm(A,2),norm(A,1),norm(A,Inf)]ans=34343434

符号运算工具箱未提供norm()函数,需先用double()函数转换成双精度数值矩阵,再调用norm()函数。第19页,共94页,2023年,2月20日,星期二特征多项式格式:C=poly(A)例:>>A=[162313;511108;97612;414151];>>poly(A)%直接求取ans=1.000000000000000e+000-3.399999999999999e+001-7.999999999999986e+0012.719999999999999e+003-2.819840539024018e-012>>A=sym(A);poly(A)%运用符号工具箱

ans=x^4-34*x^3-80*x^2+2720*x第20页,共94页,2023年,2月20日,星期二矩阵多项式的求解第21页,共94页,2023年,2月20日,星期二符号多项式与数值多项式的转换格式:

f=poly2sym(P)

或f=poly2sym(P,x)

格式:P=sym2poly(f)第22页,共94页,2023年,2月20日,星期二例:>>P=[123456];%先由系数按降幂顺序排列表示多项式>>f=poly2sym(P,'v')%以v为算子表示多项式f=v^5+2*v^4+3*v^3+4*v^2+5*v+6>>P=sym2poly(f)P=123456第23页,共94页,2023年,2月20日,星期二矩阵的逆矩阵格式:C=inv(A)例:>>formatlong;H=hilb(4);H1=inv(H)H1=1.0e+003*0.01600000000000-0.119999999999990.23999999999998-0.13999999999999-0.119999999999991.19999999999990-2.699999999999761.679999999999840.23999999999998-2.699999999999766.47999999999940-4.19999999999961-0.139999999999991.67999999999984-4.199999999999612.79999999999974第24页,共94页,2023年,2月20日,星期二检验:>>H*H1ans=1.000000000000010.00000000000023-0.000000000000450.000000000000230.000000000000011.00000000000011-0.000000000000110.000000000000110.0000000000000101.0000000000001100.000000000000000.00000000000011-0.000000000000111.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页,共94页,2023年,2月20日,星期二>>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:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=2.339949e-018.ans=53.23696008570294>>H2=invhilb(13);norm(H*H2-eye(size(H)))ans=11.37062973181391对接近于奇异矩阵,高阶一般不建议用inv(),可用符号工具箱。第26页,共94页,2023年,2月20日,星期二>>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页,共94页,2023年,2月20日,星期二例:奇异阵求逆>>A=[162313;511108;97612;414151];>>formatlong;B=inv(A)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.306145e-017.B=1.0e+014*0.938249922368852.81474976710656-2.81474976710656-0.938249922368852.814749767106568.44424930131968-8.44424930131968-2.81474976710656-2.81474976710656-8.444249301319688.444249301319682.81474976710656-0.93824992236885-2.814749767106562.814749767106560.93824992236885>>norm(A*B-eye(size(A)))%检验ans=1.64081513306419>>A=sym(A);inv(A)

%奇异矩阵不存在一个相应的逆矩阵,用符号工具箱的函数也不行???Errorusing==>sym/invError,(ininverse)singularmatrix第28页,共94页,2023年,2月20日,星期二同样适用于含有变量的矩阵求逆。例:>>symsa1a2a3a4;>>C=[a1a2;a3a4];>>inv(C)

ans=

[-a4/(-a1*a4+a2*a3),a2/(-a1*a4+a2*a3)][a3/(-a1*a4+a2*a3),-a1/(-a1*a4+a2*a3)]第29页,共94页,2023年,2月20日,星期二矩阵的相似变换与正交矩阵

其中:A为一方阵,B矩阵非奇异。%上式B’相似变换后,X矩阵的秩、迹、行列式与特征值等均不发生变化,其值与A矩阵完全一致。对于一类特殊的相似变换满足如下条件,称为正交基矩阵。第30页,共94页,2023年,2月20日,星期二例:>>A=[5,9,8,3;0,3,2,4;2,3,5,9;3,4,5,8];>>Q=orth(A)Q=-0.61970.7738-0.0262-0.1286-0.2548-0.15510.94900.1017-0.5198-0.5298-0.1563-0.6517-0.5300-0.3106-0.27250.7406>>norm(Q'*Q-eye(4))ans=4.6395e-016>>norm(Q*Q'-eye(4))ans=4.9270e-016第31页,共94页,2023年,2月20日,星期二>>C=Q'*A*QC=17.92516.4627-4.4714-2.0354-0.02821.71944.6816-5.07350.6800-0.93861.06740.6631-0.05490.36580.17760.2882>>det(A),det(C)ans=120ans=120.0000第32页,共94页,2023年,2月20日,星期二>>trace(A),trace(C)ans=21ans=21.0000>>rank(A),rank(C)ans=4ans=4第33页,共94页,2023年,2月20日,星期二>>eig(A),eig(C)%特征值求解ans=17.82051.1908+2.6499i1.1908-2.6499i0.7979ans=17.82051.1908+2.6499i1.1908-2.6499i0.7979

第34页,共94页,2023年,2月20日,星期二例:>>A=[16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1];>>Q=orth(A)%A为奇异矩阵,故得出的Q为长方形矩阵Q=-0.50000.67080.5000-0.5000-0.2236-0.5000-0.50000.2236-0.5000-0.5000-0.67080.5000>>norm(Q'*Q-eye(3))%Q*Q‘=~Ians=1.0140e-015第35页,共94页,2023年,2月20日,星期二4.2线性方程组直接解法

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

x=A\b第36页,共94页,2023年,2月20日,星期二例:解方程组>>A=[.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129;.3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927];>>b=[0.40430.15500.4240-0.2557]';>>x=A\b;x'ans=-0.1819-1.66302.2172-0.4467第37页,共94页,2023年,2月20日,星期二4.2.2线性方程组直接求解-判定求解第38页,共94页,2023年,2月20日,星期二第39页,共94页,2023年,2月20日,星期二例:>>A=[1234;4321;1324;4132];B=[51;42;33;24];>>C=[AB];rank(A),rank(C)ans=4ans=4>>x=inv(A)*B%A\B

x=-1.80002.40001.8667-1.26673.8667-3.2667-2.13332.7333第40页,共94页,2023年,2月20日,星期二检验>>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第41页,共94页,2023年,2月20日,星期二原方程组对应的齐次方程组的解求取A矩阵的化零矩阵:格式:Z=null(A)求取A矩阵的化零矩阵的规范形式:格式:Z=null(A,‘r’)第42页,共94页,2023年,2月20日,星期二例:判断可解性>>A=[1234;2211;2468;4422];B=[1;3;2;6];>>C=[AB];[rank(A),rank(C)]ans=22>>Z=null(A,'r')%解出规范化的化零空间Z=2.00003.0000-2.5000-3.5000%

1.0000001.0000第43页,共94页,2023年,2月20日,星期二>>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.4409e-015第44页,共94页,2023年,2月20日,星期二解析解>>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页,共94页,2023年,2月20日,星期二验证得出的解>>a1=randn(1);a2=rand(1);%取不同分布的随机数>>x=a1*Z(:,1)+a2*Z(:,2)+x0;norm(double(A*x-B))ans=0通解>>symsa1a2;>>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/131]第46页,共94页,2023年,2月20日,星期二

摩尔-彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。第47页,共94页,2023年,2月20日,星期二4.2.3线性方程组的直接求解分析LU分解

第48页,共94页,2023年,2月20日,星期二第49页,共94页,2023年,2月20日,星期二第50页,共94页,2023年,2月20日,星期二格式

[l,u,p]=lu(A)

L是一个单位下三角矩阵,u是一个上三角矩阵,

p是代表选主元的置换矩阵。故:Ax=y=>PAx=Py=>LUx=Py=>PA=LU

[l,u]=lu(A)其中l等于P-1L,u等于U,所以(P-1L)U=A第51页,共94页,2023年,2月20日,星期二例:对A进行LU分解>>A=[123;241;467];>>[l,u,p]=lu(A)l=1.0000000.50001.000000.25000.50001.0000u=4.00006.00007.000001.0000-2.5000002.5000p=001010100第52页,共94页,2023年,2月20日,星期二>>[l,u]=lu(A)%l=P-1Ll=0.25000.50001.00000.50001.000001.000000u=4.00006.00007.000001.0000-2.5000002.5000第53页,共94页,2023年,2月20日,星期二QR分解

将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。求得正交矩阵Q和上三角阵R,Q和R满足A=QR。

格式:

[Q,R]=qr(A)第54页,共94页,2023年,2月20日,星期二例:>>A=[123;456;789;101112];>>[Q,R]=qr(A)Q=-0.0776-0.83310.5456-0.0478-0.3105-0.4512-0.69190.4704-0.5433-0.0694-0.2531-0.7975-0.77620.31240.39940.3748R=-12.8841-14.5916-16.29920-1.0413-2.082600-0.0000000第55页,共94页,2023年,2月20日,星期二Cholesky(乔里斯基)分解若矩阵A为n阶对称正定阵,则存在唯一的对角元素为正的三角阵D,使得

第56页,共94页,2023年,2月20日,星期二格式:

D=chol(A)第57页,共94页,2023年,2月20日,星期二例:进行Cholesky分解。>>A=[1648;45-4;8-422];>>D=chol(A)D=41202-3003第58页,共94页,2023年,2月20日,星期二利用矩阵的LU、QR和cholesky分解求方程组的解

(1)LU分解:

A*X=b变成L*U*X=b所以X=U\(L\b)这样可以大大提高运算速度。例:求方程组的一个特解。解:>>A=[42-1;3-12;1130];>>B=[2108]';>>D=det(A)D=0第59页,共94页,2023年,2月20日,星期二>>[L,U]=lu(A)L=0.3636-0.50001.00000.27271.000001.000000U=11.00003.000000-1.81822.0000000.0000第60页,共94页,2023年,2月20日,星期二>>X=U\(L\B)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=2.018587e-017.X=1.0e+016*%结果中的警告是由于系数行列式为零产生的。

-0.4053%可以通过A*X验证其正确性。

1.48621.3511>>A*X%Matlab7.0显示没有解ans=088第61页,共94页,2023年,2月20日,星期二(2)Cholesky分解若A为对称正定矩阵,则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积,方程A*X=b变成R’*R*X=b所以X=R\(R’\b)(3)QR分解对于任何长方矩阵A,都可以进行QR分解,其中Q为正交矩阵,R为上三角矩阵的初等变换形式,即:A=QR方程A*X=b变形成QRX=b所以X=R\(Q\b)

这三种分解,在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。第62页,共94页,2023年,2月20日,星期二三个变换在线性方程组的迭代求解中,要用到系数矩阵A的上三角矩阵、对角阵和下三角矩阵。此三个变换在MATLAB中可由以下函数实现。上三角变换:格式triu(A,1)对角变换:格式diag(A)下三角变换:格式tril(A,-1)例:对此矩阵做三种变换。第63页,共94页,2023年,2月20日,星期二>>A=[12-2;111;221];%>>triu(A,1)ans=02-2001000>>tril(A,-1)ans=000100220>>b=diag(A);b'ans=111第64页,共94页,2023年,2月20日,星期二4.3迭代解法的几种形式

5.3.1Jacobi迭代法方程组Ax=bA可写成A=D-L-U

其中:D=diag[a11,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.第65页,共94页,2023年,2月20日,星期二functiony=jacobi(a,b,x0)%x0为初值D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;%n为迭代次数whilenorm(y-x0)>=1.0e-6x0=y;y=B*x0+f;n=n+1;endn第66页,共94页,2023年,2月20日,星期二例:用Jacobi方法求解,设x(0)=0,精度为10-6。>>a=[10-10;-110-2;0-210];>>b=[9;7;6];>>jacobi(a,b,[0;0;0])n=11ans=0.99580.95790.7916第67页,共94页,2023年,2月20日,星期二4.3.2Gauss-Seidel迭代法由原方程构造迭代方程

x(k+1)=Gx(k)+f其中:G=(D-L)-1U,f=(D-L)-1bD=diag[a11,a22,…,ann],-L、-U分别为A的严格下、上三角部分(不包括对角线元素).第68页,共94页,2023年,2月20日,星期二functiony=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;whilenorm(y-x0)>=1.0e-6x0=y;y=G*x0+f;n=n+1;endn第69页,共94页,2023年,2月20日,星期二例:对上例用Gauss-Seidel迭代法求解>>a=[10-10;-110-2;0-210];>>b=[9;7;6];>>seidel(a,b,[0;0;0])n=7ans=0.99580.95790.7916例:分别用Jacobi和G-S法迭代求解,看是否收敛。第70页,共94页,2023年,2月20日,星期二>>a=[12-2;111;221];b=[9;7;6];>>jacobi(a,b,[0;0;0])n=4ans=-27268>>seidel(a,b,[0;0;0])n=1011ans=1.0e+305*-InfInf-1.7556第71页,共94页,2023年,2月20日,星期二4.3.3SOR迭代法

在很多情况下,J法和G-S法收敛较慢,所以考虑对G-S法进行改进。于是引入一种新的迭代法-逐次超松弛迭代法(SuccesiseOver-Relaxation),记为SQR法。迭代公式为:

X(k+1)=(D-wL)-1((1-w)D+wU)x(k)

+w(D-wL)-1b

其中:w最佳值在[1,2)之间,不易计算得到,因此w通常有经验给出。第72页,共94页,2023年,2月20日,星期二functiony=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;whilenorm(y-x0)>=1.0e-6x0=y;y=M*x0+f;n=n+1;endn第73页,共94页,2023年,2月20日,星期二例:上例中,当w=1.103时,用SOR法求解原方程。>>a=[10-10;-110-2;0-210];>>b=[9;7;6];>>sor(a,b,1.103,[0;0;0])n=8ans=0.99580.95790.7916第74页,共94页,2023年,2月20日,星期二4.3.4两步迭代法

当线性方程系数矩阵为对称正定时,可用一种特殊的迭代法来解决,其迭代公式为:(D-L)x(k+1/2)=Ux(k)+b

(D-U)x(k+1)=Lx(k+1/2)+b=>x(k+1/2)=(D-L)-1Ux(k)+(D-L)-1bx(k+1)=(D-U)-1Lx(k+1/2)+(D-U)-1b第75页,共94页,2023年,2月20日,星期二functiony=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;whilenorm(y-x0)>=1.0e-6x0=y;y=G1*x0+f1;y=G2*y+f2;n=n+1;endn第76页,共94页,2023年,2月20日,星期二例:求解方程组>>a=[10-120;-111-13;2-1103;03-18];>>b=[6;25;-11;15];>>twostp(a,b,[0;0;0;0])n=7ans=1.07911.9824-1.40440.9560第77页,共94页,2023年,2月20日,星期二4.4线性方程组的符号解法

在MATLAB的SymbolicToolbox中提供了线性方程的符号求解函数,如

linsolve(A,b)

等同于X=sym(A)\sym(b).

solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')第78页,共94页,2023年,2月20日,星期二例:>>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]第79页,共94页,2023年,2月20日,星期二例:>>[x,y]=solve('x^2+x*y+y=3','x^2-4*x+3=0','x','y')x=[1][3]

y=[1][-3/2]

第80页,共94页,2023年,2月20日,星期二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第81页,共94页,2023年,2月20日,星期二例:>>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=0000010000010000010000010第82页,共94页,2023年,2月20日,星期二例: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第83页,共94页,2023年,2月20日,星期二(3,4)1(4,4)4(5,4)1(4,5)1(5,5)4>>full(a)ans=4100014100014100014100014第84页,共94页,2023年,2月20日,星期二格式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=4100014100014100014100014第85页,共94页,2023年,2月20日,星期二格式:spconvert(dd)

对于无规律的稀疏矩阵,可使用此命令由外部数据转化为稀疏矩阵。调用形式为:先用load函数加载以行表示对应位置和元素值的.dat文本文件,再用此命令转化为稀疏矩阵。例:无规律稀疏矩阵的建立。首先编制文本文件sp.dat如下:515.00358.00442.00550第86页,共94页,2023年,2月20日,星期二>>loadsp.dat>>spconvert(sp)ans=(5,1)5(4,4)2(3,5)8>>full(ans)ans=000000000000

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论