matlab在科學(xué)計(jì)算中的應(yīng)用第4章 線性代數(shù)問(wèn)題求解

第四章線性代數(shù)問(wèn)題求解矩陣線性方程組的直接解法線性方程組的迭代法線性方程組的符號(hào)解法稀疏矩陣技術(shù)特征值與特征向量14.1矩陣

特殊矩陣的輸入數(shù)值矩陣的輸入零矩陣、幺矩陣及單位矩陣生成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同樣位數(shù)的矩陣：A=zeros(size(B))2隨機(jī)元素矩陣假設(shè)矩陣隨機(jī)元素滿足[0,1]區(qū)間上的均勻分布生成nm階標(biāo)準(zhǔn)均勻分布偽隨機(jī)數(shù)矩陣：A=rand(n,m)生成nn階標(biāo)準(zhǔn)均勻分布偽隨機(jī)數(shù)方陣：A=rand(n)3對(duì)角元素矩陣向量生成對(duì)角矩陣：A=diag(V)矩陣提取對(duì)角元素列向量：V＝diag(A)生成主對(duì)角線上第k條對(duì)角線為V的矩陣：A=diag(V,k)4例：diag()函數(shù)的不同調(diào)用格式>>C=[123];V=diag(C)%生成對(duì)角矩陣V=100020003>>V1=diag(V)'%將列向量通過(guò)轉(zhuǎn)置變換成行向量V1=123>>C=[123];V=diag(C,2)%主對(duì)角線上第k條對(duì)角線為C的矩陣V=00100000200000300000000005生成三對(duì)角矩陣:>>V=diag([1234])+diag([234],1)+diag([543],-1)V=12005230043400346Hilbert矩陣及逆Hilbert矩陣

生成n階的Hilbert矩陣：

A=hilb(n)

求取逆Hilbert矩陣：B=invhilb(n)7Hankel(漢克)矩陣

其中：第一列的各個(gè)元素定義為C向量，最后一行各個(gè)元素定義為R。H為對(duì)稱陣。

H1=hankel(C)由Hankel矩陣反對(duì)角線上元素相等得出一下三角陣均為零的Hankel矩陣8Vandermonde(范德蒙)矩陣

9伴隨矩陣其中：P(s)為首項(xiàng)系數(shù)為1的多項(xiàng)式。

10例：考慮一個(gè)多項(xiàng)式2*x^4+4*x^2+5*x+6，試寫出該多項(xiàng)式的伴隨矩陣。>>P=[20456];A=compan(P)A=0-2.0000-2.5000-3.00001.000000001.000000001.0000011符號(hào)矩陣的輸入數(shù)值矩陣A轉(zhuǎn)換成符號(hào)矩陣：

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]124.1.2矩陣根本概念與性質(zhì)行列式格式：d=det(A)例：求行列式>>A=[162313;511108;97612;414151];det(A)ans=013例：>>tic,A=sym(hilb(20));det(A),tocans=elapsed_time=2.3140高階的Hilbert矩陣是接近奇異的矩陣。14矩陣的跡格式：t=trace(A)矩陣的秩格式：r=rank(A)％用默認(rèn)的精度求數(shù)值秩r=rank(A，)％給定精度下求數(shù)值秩矩陣的秩也表示該矩陣中行列式不等于0的子式的最大階次。可證行秩和列秩〔線性無(wú)關(guān)的〕應(yīng)相等。15例>>A=[162313;511108;97612;414151];rank(A)ans=3該矩陣的秩為3，小于矩陣的階次，故為非滿秩矩陣。例>>H=hilb(20);rank(H)％數(shù)值方法ans=13>>H=sym(hilb(20));rank(H)%解析方法，原矩陣為非奇異矩陣ans=2016矩陣范數(shù)17矩陣的范數(shù)定義：格式：

N=norm(A)％求解默認(rèn)的2范數(shù)N=norm(A，選項(xiàng))％選項(xiàng)可為1,2,inf等18例：求一向量、矩陣的范數(shù)>>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符號(hào)運(yùn)算工具箱未提供norm()函數(shù)，需先用double()函數(shù)轉(zhuǎn)換成雙精度數(shù)值矩陣，再調(diào)用norm()函數(shù)。19特征多項(xiàng)式格式：C=poly(A)例：>>A=[162313;511108;97612;414151];>>poly(A)％直接求取ans=1.000000000000000e+000-3.399999999999999e+001>>A=sym(A);poly(A)％運(yùn)用符號(hào)工具箱ans=x^4-34*x^3-80*x^2+2720*x20矩陣多項(xiàng)式的求解21符號(hào)多項(xiàng)式與數(shù)值多項(xiàng)式的轉(zhuǎn)換格式：

f=poly2sym(P)或f=poly2sym(P,x)

格式：P=sym2poly(f)22例：>>P=[123456];%先由系數(shù)按降冪順序排列表示多項(xiàng)式>>f=poly2sym(P,'v')%以v為算子表示多項(xiàng)式f=v^5+2*v^4+3*v^3+4*v^2+5*v+6>>P=sym2poly(f)P=12345623矩陣的逆矩陣格式：C=inv(A)例：>>formatlong;H=hilb(4);H1=inv(H)H1=1.0e+003*24檢驗(yàn)：>>H*H1ans=1.000000000000010.00000000000023-0.000000000000450.000000000000230.000000000000011.00000000000011-0.000000000000110.000000000000110.0000000000000101.0000000000001100.000000000000000.00000000000011-0.000000000000111.00000000000011計(jì)算誤差范數(shù)：>>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:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=2.339949e-018.ans=

53.23696008570294>>H2=invhilb(13);norm(H*H2-eye(size(H)))ans=

11.37062973181391對(duì)接近于奇異矩陣，高階一般不建議用inv()，可用符號(hào)工具箱。26>>H=sym(hilb(7));inv(H)ans=[49,-1176,8820,-29400,48510,-38808,12021][-1176,37632,-317520,1128960,-1940400,1596672,-504504][8820,-317520,2857680,-10584000,18711000,-15717240,5045040][-29400,1128960,-10584000,40320000,-72765000,62092800,-20210160][48510,-1940400,18711000,-72765000,133402500,-115259760,37837800][-38808,1596672,-15717240,62092800,-115259760,100590336,-33297264][12021,-504504,5045040,-20210160,37837800,-33297264,11099088]>>H=sym(hilb(30));norm(double(H*inv(H)-eye(size(H))))ans=027例：奇異陣求逆>>A=[162313;511108;97612;414151];>>formatlong;B=inv(A)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.306145e-017.B=1.0e+014*>>norm(A*B-eye(size(A)))％檢驗(yàn)ans=1.64081513306419>>A=sym(A);inv(A)％奇異矩陣不存在一個(gè)相應(yīng)的逆矩陣，用符號(hào)工具箱的函數(shù)也不行???Errorusing==>sym/invError,(ininverse)singularmatrix28同樣適用于含有變量的矩陣求逆。例：>>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矩陣的相似變換與正交矩陣

其中：A為一方陣，B矩陣非奇異。相似變換后，X矩陣的秩、跡、行列式與特征值等均不發(fā)生變化，其值與A矩陣完全一致。對(duì)于一類特殊的相似變換滿足如下條件，稱為正交基矩陣。30例：>>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-01631>>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.000032>>trace(A),trace(C)ans=21ans=21.0000>>rank(A),rank(C)ans=4ans=433>>eig(A),eig(C)ans=17.82051.1908+2.6499i1.1908-2.6499i0.7979ans=17.82051.1908+2.6499i1.1908-2.6499i0.7979

34例：>>A=[16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1];>>Q=orth(A)％A為奇異矩陣，故得出的Q為長(zhǎng)方形矩陣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-015354.2線性方程組直接解法

線性方程組直接求解－矩陣除法關(guān)于線性方程組的直接解法，如Gauss消去法、選主元消去法、平方根法、追趕法等等，在MATLAB中，只需用“／〞或“\〞就解決問(wèn)題。它內(nèi)部實(shí)際包含著許許多多的自適應(yīng)算法，如對(duì)超定方程用最小二乘法，對(duì)欠定方程時(shí)它將給出范數(shù)最小的一個(gè)解，解三對(duì)角陣方程組時(shí)用追趕法等等。格式：x=A\b36例：解方程組>>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.446737線性方程組直接求解－判定求解3839例：>>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.733340檢驗(yàn)>>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]檢驗(yàn)>>norm(double(A*x1-B))ans=041原方程組對(duì)應(yīng)的齊次方程組的解求取A矩陣的化零矩陣：格式：Z=null(A)求取A矩陣的化零矩陣的標(biāo)準(zhǔn)形式：格式：Z=null(A,‘r’)null是用來(lái)求齊次線性方程組的根底解系的，加上'r'那么求出的是一組最小正整數(shù)解，如果不加，那么求出的是解空間的標(biāo)準(zhǔn)正交基。42例：判斷可解性>>A=[1234;2211;2468;4422];B=[1;3;2;6];>>C=[AB];[rank(A),rank(C)]ans=22>>Z=null(A,'r')%解出標(biāo)準(zhǔn)化的化零空間Z=2.00003.0000-2.5000-3.5000％1.0000001.000043>>x0=pinv(A)*B%得出一個(gè)特解x0=0.95420.7328%全部解-0.0763-0.2977驗(yàn)證得出的解>>a1=randn(1);a2=rand(1);%取不同分布的隨機(jī)數(shù)>>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驗(yàn)證得出的解>>a1=randn(1);a2=rand(1);%取不同分布的隨機(jī)數(shù)>>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摩爾－彭羅斯廣義逆求解出的方程最小二乘解不滿足原始代數(shù)方程。474.2.3線性方程組的直接求解分析LU分解

484950格式

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

L是一個(gè)單位下三角矩陣，u是一個(gè)上三角矩陣，p是代表選主元的置換矩陣。故：Ax=y=>PAx=Py=>LUx=Py=>PA=LU

[l,u]=lu(A)其中l(wèi)等于P-1L，u等于U，所以(P-1L)U=A51例：對(duì)A進(jìn)行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=00101010052>>[l,u]=lu(A)％l＝P-1Ll=0.25000.50001.00000.50001.000001.000000u=4.00006.00007.000001.0000-2.5000002.500053QR分解

將矩陣A分解成一個(gè)正交矩陣與一個(gè)上三角矩陣的乘積。求得正交矩陣Q和上三角陣R，Q和R滿足A=QR。

格式：

[Q,R]=qr(A)54例：>>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.000000055Cholesky(喬里斯基)分解假設(shè)矩陣A為n階對(duì)稱正定陣,那么存在唯一的對(duì)角元素為正的三角陣D，使得

56格式：

D=chol(A)57例：進(jìn)行Cholesky分解。>>A=[1648;45-4;8-422];>>D=chol(A)D=41202-300358利用矩陣的LU、QR和cholesky分解求方程組的解

〔1〕LU分解：A*X=b變成L*U*X=b所以X=U\(L\b)這樣可以大大提高運(yùn)算速度。例：求方程組的一個(gè)特解。解：>>A=[42-1;3-12;1130];>>B=[2108]';>>D=det(A)D=059>>[L,U]=lu(A)L=0.3636-0.50001.00000.27271.000001.000000U=11.00003.000000-1.81822.0000000.000060>>X=U\(L\B)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=2.018587e-017.X=1.0e+016*%結(jié)果中的警告是由于系數(shù)行列式為零產(chǎn)生的。-0.4053%可以通過(guò)A*X驗(yàn)證其正確性。

1.48621.3511>>A*X%Matlab7.0顯示沒(méi)有解ans=08861〔2〕Cholesky分解假設(shè)A為對(duì)稱正定矩陣，那么Cholesky分解可將矩陣A分解成上三角矩陣和其轉(zhuǎn)置的乘積，方程A*X=b變成R’*R*X=b所以X=R\(R’\b)〔3〕QR分解對(duì)于任何長(zhǎng)方矩陣A，都可以進(jìn)行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：A=QR方程A*X=b變形成QRX=b所以X=R\(Q\b)

這三種分解，在求解大型方程組時(shí)很有用。其優(yōu)點(diǎn)是運(yùn)算速度快、可以節(jié)省磁盤空間、節(jié)省內(nèi)存。62三個(gè)變換在線性方程組的迭代求解中，要用到系數(shù)矩陣A的上三角矩陣、對(duì)角陣和下三角矩陣。此三個(gè)變換在MATLAB中可由以下函數(shù)實(shí)現(xiàn)。上三角變換：格式triu(A,1)對(duì)角變換：格式diag(A)下三角變換：格式tril(A,-1)例：對(duì)此矩陣做三種變換。63>>A=[12-2;111;221];％>>triu(A,1)ans=02-2001000>>tril(A,-1)ans=000100220>>b=diag(A);b'ans=111644.4線性方程組的符號(hào)解法在MATLAB的SymbolicToolbox中提供了線性方程的符號(hào)求解函數(shù)，如

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('x^2+x*y+y=3','x^2-4*x+3=0','x','y')

x=[1][3]

y=[1][-3/2]

674.5稀疏矩陣技術(shù)稀疏矩陣的建立：格式S=sparse(i,j,s,m,n)生成一mxn階的稀疏矩陣，以向量i和j為坐標(biāo)的位置上對(duì)應(yīng)元素值為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)468例：>>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=000001000001000001000001069例：n=5,建立主對(duì)角線上元素為4，兩條次對(duì)角線為1的三對(duì)角陣。>>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)170(3,4)1(4,4)4(5,4)1(4,5)1(5,5)4>>full(a)ans=410001410001410001410001471格式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=410001410001410001410001472格式：spconvert(dd)對(duì)于無(wú)規(guī)律的稀疏矩陣，可使用此命令由外部數(shù)據(jù)轉(zhuǎn)化為稀疏矩陣。調(diào)用形式為：先用load函數(shù)加載以行表示對(duì)應(yīng)位置和元素值的.dat文本文件，再用此命令轉(zhuǎn)化為稀疏矩陣。例：無(wú)規(guī)律稀疏矩陣的建立。首先編制文本文件sp.dat如下：515.00358.00442.0055073>>loadsp.dat>>spconvert(sp)ans=(5,1)5(4,4)2(3,5)8>>full(ans)ans=000000000000008000205000074稀疏矩陣的計(jì)算：同滿矩陣比較，稀疏矩陣在算法上有很大的不同。具體表現(xiàn)在存儲(chǔ)空間減少，計(jì)算時(shí)間減少。例：比較求解下面方程組n＝1000時(shí)兩種方法的差異。75>>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=a\b;t1=toct1=0.4800>>a=full(a);>>tic;x=a\b;t2=toct2=1.3220764.6矩陣的特征值問(wèn)題

一般矩陣的特征值與特征向量格式：d=eig(A)只求解特征值。格式：[V,D]=eig(A)求解特征值和特征向量。77例：直接求解：>>A=[162313;511108;97612;414151];>>eig(A)ans=34.00008.9443-8.94430.000078精確解：>>eig(sym(A))ans=[0][34][4*5^(1/2)][-4*5^(1/2)]高精度數(shù)值解：>>vpa(ans,70)ans=[0][34.]2897083588981642084]79同時(shí)求出特征值與特征向量：直接求解：>>[v,d]=eig(A)v=-0.5000-0.82360.3764-0.2236-0.50000.42360.0236-0.6708-0.50000.02360.42360.6708-0.50000.3764-0.82360.2236d=

34.00000000