算法与编程手册-matlab6.0数学

1章矩阵及其基本运，即“矩阵它是以矩阵为基本运算单元。因此，本书从最基本的运算单元出发，介绍令及其用法。的强大功能之一体现在能直接处理向量或矩阵。当然首要任务是输入待处理（）（）一方括号（[]）内；当矩阵是（三维以上，且方括号内的元素是维数较低的矩阵时，>>Time=[1112123456789Time1112123456789>>X_Data=[2.323.43；4.37X_Data2.434.37>>vect_a=[12345]vect_a=1234>>Matrix_B=[12 234；34Matrix_B=23234345>>Null_M=[第式1->>>>C=[1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+1]5.4000+21->>R=[123;456],M=[111213;141516]R=123456M>>CN=R+i*MCN=1.00002.00003.00004.00005.00006.0000在中输入符号向量或者矩阵的方法和输入数值类型的向量或者矩阵在形式上symsyms，先定义一sym1->>sym_matrix=sym（'[abc；Jack，HelpMe!，NOWAY!]，'）sym_matrix= HelpMe! >>sym_digits=sym（'[123；absym_digits syms1->>symsabc>>M1=>>M2=sym（'>>M3=>>syms_matrix=[abc；M1，M2，M3；int2str（[235]）]syms_matrix= [ClassicalJazzBlues] 3 5]数值型和符号型在中是不相同的，它们之间不能直接进行转化。供了一个将数值型转化成符号型令，即sym。1->>Digit_Matrix=[1/3sqrt（2）3.4234；exp（0.23）log（29）23^（->>Syms_Matrix=Digit_MatrixSyms_Matrix 例1- exm=[285596在窗口输入 %显示exm5 %表示exm56列函数格式说明n=1n=2时分别构造[A1；A2]和[A1，A2]n=3时可以构造1-A4(:,:,1)=123456789A4(:,:,2)147258369A4(:,:,3) - - - 1->>A5(:,:,1)=A1,A5(:,:,2)=A2,A5(:,:,3)=A3A5(:,:,1)=123456789A5(:,:,2)147258369A5(:,:,3) - - - B=%生成n×nB=zeros(m,n)B=zeros([mn])B=zeros(d1,d2,d3…)B=zeros([d1d2d3…])B=zeros(size(A))%m×n%m×n%生成d1×d2×d3× %n×n 1%m×nY=%n×n1Y=ones(m,n)Y=ones([mn])Y=ones(d1,d2,d3…)Y=ones([d1d2d3…])Y=ones(size(A))%m×n1%m×n1%A1格式Y= Y Yrand([m %m×nY=rand(m,n,p,…) Y=rand([mn Y s rand('state rand('state rand('state %jjrand('statesum 例1- 产生一个3×4随机矩>>R=rand(3,4)R=例1- 产生一个在区间[10,20]内均匀分布的4阶随机矩>>>>x=a+(b-a)*rand(4)x=命令函数格式Y=randn(n) %生成n×n正态分布随机矩阵Y=randn(m,n) Y=randn([m Yrandn(m,n,p,…)%m×n×p×Y=randn([mnp…]) Y=randn(size(A)) %生成与矩阵A相同大小的正态分布随机矩阵 s= srandn('state srandn('state srandn('state %jjsrandn('state 例1- 产生均值为0.6，方差为0.1的4阶矩>>mu=0.6;>>xp=1->>randperm(6)ans= 命令函数格式y= %在(a,b)上产生100个线性等分y= %在(ab)n命令函数格式y= %在（10ay=logspace(a,b,n)y=logspace(a,pi)命令

10b）之间产生50个对数等分向n= %A命令函数格式out= 1->>out=blkdiag(1,2,3,4)out=1000020000300004A=%u为多项式系统向量，A为友矩阵，A的第-u2:n)/u(1)u2:n)u2n个元素，A为特例1- 求多项式(x1)(x2)(x3)x37x6的友矩阵和>>u=[10-7A07-100 >> %Aans-命令hadamard函数格式H= %返回n阶hadamard矩1->>h=hadamard(4)h=11111-111-1--1命令Hankel函数格式H= %第1列元素为c，反三角以下元素为0H=hankel(c,r) %1crc的最后一个元素与rc1->>c r >>h=hankel(c,r)h=12382389389命令Hilbert函数格式H= %返回n阶Hilbert矩阵，其元素为H(i,j)=1/(i+j-1)例1- 产生一个3阶Hilbert矩>>format>>H=hilb(3)H=%以有理形式输1命令逆Hilbert函数格式H= %产生n阶逆Hilbert矩命令Magic(魔方)函数格式M=1->>M=magic(3)M=%816357492命令Pascal函数格式A=A=A=1-A=11111234A=1001-01-1A=111-0100T=1T 1->>c=[1234>>r=[1.52.53.54.5>>T=toeplitz(c,r)T=121321432154321命令Wilkinson函数格式W= %返回n阶Wilkinson特征值测试1->>W100110011001>>W31000001210000011100000101000001110000012100000131->>A=[1,1,1;1,2,3;1,3,>>B=[8,1,6;3,5,7;4,9,9274758-0---41->>X=3412211100100]；8533则显示：a=468244A.*B表示AB函数格式C= %若A、B为向量，则返回向量A与B的点积，A与B长相同；若为矩阵，则ABC %dim维数中给出AB >>X=[-10>>Y=[-2-1>>Z=dot(X,则显示：Z44在中，用函数cross实现。函数Ccross(A,B)%A、B为向量，则返回ABC=A×B，A、B3A、B3×nAB对应列的叉积，A、B3×n矩C=cross(A,B,dim) %dimA与B的叉积。AB必须具有相同的维数，size(A,dim)size(B,dim)3。例1- 计算垂直于向量(1,2,3)和(4,5,6)的向量>>a=[12>>b=[45 -可得垂直于向量(1,2,3)和(4,5,6)的向量为±(-3,6,3)例1- 计算向量a=(1,2,3)、b=(4,5,6)和c=(-3,6,-3)的混合积a>>a=[123];b=[456];c=[-36->>x=dot(a,cross(b,结果显示：x函数格式w= 说明munv的卷积(Convolution)kwku(jv(k1j式中：w向量序列的长度为(m+n-1)m=nw(1)=w(2)=w(3)=…w(n)=u(1)*v(n)+u(2)*v(n-1)+……w(2*n-1)=例1- 展开多项式(s22s2)(s4)(s解：>>w >> %将wPs^4+7s^3+16s^2+18s+8函数deconv格式[q,r]= %多项式v除以多项式u，返回商多项式q和余多项式r(x32x23x4)(10x220x例1- >>u= >>v= >>c=c >>[q,r]=deconv(c,u)q=r8

函数格式C=kron %A为m×n矩阵，B为p×q矩阵，则C为mp×nq矩阵说明ABCAB

ABB mp×nq矩阵，但一般地ABBA

am2B

A 44

B

6 6

AB>>A=[12;34];B=[123;456;78>>C=kron(A,B)C=1232464568789 函数格式c= %返回向量a、b的公共部分，即c=a∩bc= [c,ia,ib] ibb1->>A=[1234;1246;6714]A=123412466714>>B=[1238;1146;6714]B=12381146 >>C=intersect(A,B,'rows')C=

>>A=[19620];B=[1234610>>[c,ia,ib]=intersect(A,B)c= ia ib 函数格式k= %当a中元素属于S时，k取1，否则，k取0k= 1->>S=[024681012141618>>a=[12345>>k=ismember(a,S)k=

%1表示相同元素的位>>A=[1234;1246;671>>B=[1238;1146;671>>k=ismember(A,B,'rows')k=00 %1表示元素相同的函数格式c=setdiff(a,b) %返回属于a但不属于b的不同元素的集合，C=a-b。c=setdiff(A,B,'rows') %返回属于A但不属于B的不[c,i %c与前面一致，i表示cA1->>A=[179620];B=[1234610>>c=setdiff(A,B)c= 1->>A=[1234;1246;671>>B=[1238;1146;671>>c=setdiff(A,B,'rows')c=12341246函数格式c= %返回集合a、b交集的c= %A、B交集的非，A、B[c,ia,ib] %ia、ib表示caA)、b(B)1->>A=[123>>B=[245>>C=setxor(A,B)C= 1->>A=[1234;1246;6714]A=123412466714>>B=[1238;1146;6714]B=123811466714>>[C,ia,ib]=setxor(A,B,'rows')C=1 ia12ib21函数格式c= %返回a、b的并集，即c=a∪bc= [c,ia,ib] %ia、ib分别表示c中行向量在原矩阵(向量)1->>A=[123>>B=[245>>则结果c

>>A=[1234;124A12341246>>B=[1238;1146]B=12381146>>[c,ia,ib]=union(A,B,'rows')c=1 ia12ib21格式b=unique %取集合a的不重复元素构成的向b=unique [b,i,junique %i、jb中元素在原向量（矩阵）1->>A=[112244646]A= >>[c,i,j]=unique(A)c=ij1-

>>A=[1224;1146;1146]A=122411461146>>[c,i,j]=unique(A,'rows')c= i31j211（/x=b/ax*a=b例：a=[123;426;74b=[4;1;则显示-a为非奇异矩阵，则a\bb/aaba\b=inv(a)*bb/a=A./BABA为方阵，P0的整数时，A^P表示AP次方，即AP次；P为0的整数时，A^P表示A-1P次方。当A为方阵，p为非整数时，则A^P

V1VA量，

PA，标量的矩阵乘方定义AV=AD。

ppA

pdnn

pp

a1npp标量的数组乘方P.^A，标量的数组乘方定义为P.^A 数组乘方A.^P：表示AP

pamnY=%使用特征值和特征向量计算格式Y=logm(X) %计算矩阵X的对数，它是expm(X)的反函数。[Y,esterr]=logm(X) %esterr为相对残差的估计值：norm(expm(Y)-X)/norm(X)1->>A=[110;002;00->>Y=expm(A)Y=000A=0000F=sincos’[F,esterr]=命令函数X=sqrtm(A)%A的平方根A1/2X*X=AXA的特征值XAXAX不存在。[X,resnorm resnorm[X,alpha,condest 命令矩阵A函数格式polyvalm(P, 运算规则：若矩阵AAAA对应元素的共轭复数构成。函数格式d= %返回方阵X的多项式的1->>A=[123;456;789]A=123456789>>D=det(A)D=0命令函数格式 例1- 求A

2214>>A=[123;221;34Y----B

232143

0 00>>B=[1,2,3,1,0,0;2,2,1,0,1,0;3,4,3,0,0, >>X=C %取矩阵CA^(-1)部C00-00- -X--1->>A=[21-1;212;1-1>>format>>D=inv(A)D=00-0B=%AB=pinv(A,%tol说明AX=I和XA=IA的伪逆能在某种程度上代表矩阵的逆，若Apinv(A)=inv(A)。1- %5阶魔方阵A=A(:,1:4) %5阶魔方阵的前4列元素构成矩阵A。A185746>> %计算AX----函数格式b=trace %返回矩阵A的迹，即A的对角线元和命令函数|xk格式n= %X为向量，求欧几里德范数，即|||xkn= %求-范数，即||X||max(absXn= %1-范数，即||X||1|xk|knnorm(X,- %求向量X||X||min(absXp|xk|pn=norm(X, p|xk|p命令函数格式n=norm(A) %A为矩阵，求欧几里德范数||A||2，等于A的最大奇异值。n=norm(A,1) %求A的列范数||A||1，等于A的列向量的1-范数的最大值。n=norm(A,2) %求A的欧几里德范数||A||2，和norm(A)相同。n= %求行范数||A||A1-|aij n=norm(A,'fro' %求矩阵A的Frobenius|aij 命令函数格式nrm= 差小于106nrm %tol[nrm,count 命令函数格式c= %求X的2-范数的条件数即X的最大奇异值和最小奇异值的商c 说明线性方程组AX=b1cond

|A||||A1命令1-函数格式c=condest [c,v]=condest( %v||Av||||A||||v||，即c[c,v]condestA,t)%cvt=1，则计算的每步都显示出来；如果t=-1，则给出商命令函数格式crcond(A)%对于差条件矩阵A0的数；对于好条件A1的数。命令函数格式c= %返回矩阵A的特征值的条件 %D为A的特征值对角阵，VA函数格式k=rank %求矩阵A的krank %tol函数格式Xdiag(v,k)%v的元素作为矩阵Xkk=0时，vXk>0时，vkk<0时，vk条对角线。Xdiag(v)%v0Xvdiag(X,k)%Xkv。k=0：抽取主对角线元素；k>0k条对角线元素；k<0k条vdiag(X)%v1->>v=[12x=0000100002000030>>A=[123;456;789]A=123456789>>v=diag(A,1)v=26函数 格式L= %抽取X的主对角线的下三角部分构成矩阵L= %Xk条对角线的下三角部分；k=0函数 格式U= U=triu(X,k) %Xk条对角线的上三角部分；k=0为主对角线；k>0为主对角线以上；k<0为主对角线以下。1- %产生4阶全1A1111111111111111>>L=tril(A,1) L1100111011111111>>U=triu(A,-1) U=1111111101110011reshape（11->A=[123456;678901]A=123456678901>>B=ones(3,4)B=111111111111>>B1（2）Reshape格式B=reshape(A,m,n) %返回以矩阵A的元素构成的m×n矩阵BB=reshape(A,m,n,p,…) %将矩阵A变维为m×n×p×…Breshape(A,[mn B %siz决定变维的大小，元素个数与A例1- 矩阵变>>>>b=reshape(a,2,6)b=

格式B=rot90 %将矩阵A逆时针方向旋转Brot90 %将矩阵A逆时针方向旋转(k×90°)，k1->>A=[123;456;789]A=123456789Y1 369258147Y27 8 9 函数格式B= %将矩阵A左右翻函数格式B= %将矩阵A上下翻1->>A=[123;456]A=123456>>B1=fliplr(A),B2=flipud(A)B1=321654B2456123函数格式B= 1->>A=[123;456]A=123456>>B1=flipdim(A,1),B2=flipdim(A,2)B1=456123B2321654函数格式B= %将矩阵Am×n块，即B由m×n块A平铺而成Brepmat(A,[m %Brepmat(A,[mn %Bm×n×p×A %Aaam×n1->>A=[12;56]A= >>B=repmat(A,3,4)B=12121212565656561212121256565656121212125656565610。1-含含>大于关<大于关=等于关~1->>A=[1234;5678];B=[0214;077>>C1=A==B,C2=A>=B,C3=A~=BC1= C2001011111011C310101101对于小数构成的矩阵A函数格式 %将A中元素按-∞方向取整，即取不足整数函数格式 %将A中元素按+∞方向取整，即取过剩整数函数格式round %将A中元素按最近的整数取整，即四舍五入取整0函数格式fix %将A中元素按离0近的方向取1-A=--2002002-011=31103122=2002-122=20002011函数格式[n,d]=rat 例1- -d函数格式C=rem(A,x)%Axx=0rem(A,0)=NaNx≠0则整fix(A./x)表示C=A-x.*fix(A./x)。允许模x为小数。设矩阵A和B都是m×n矩阵或其中之一为标量， 格式A&B或and(A,说明AB010格式A|Bor(A,说明AB001格式~Anot说明若A010格式xor说明A 1505>>A=[02302341350= 1505B>>C1=A&B,C2=A|B,C3=~A,C4=xor(A,B)C1=00111100C21111C3111110000001C4110000116.x4.2版中为符号矩阵设计的复杂函数形式，把符号矩阵的四则运算（＋（－（×（symadd（symsub1-58Asym([1x1/(x11/(x21/(x>>Bsym([x,1;x2,0])则显示x-1/x1-x+2- （'（det（inv（rank（^）和指数（exp和expm）等都与数值矩阵相同函数格式 %将A转化为符号矩阵1-A=>>B=sym(A)B=[ 例1- >>B(2,3) ans=>>B(2,3)='log(7)' B= 2/3,sqrt(2), 7/5,100/23,函数格式 说明S例1-61 将x9-1分解因式在命令窗口键入syms则显示：ans(1)x12x24x3例1- 问“入”取何值时，齐次方程组2x(3)xx0有非0解 x1x2(1)x3解：在编辑器中建立M文件symsA=[1-k-24;23-k1;111-Dans=k=0、k=2或k=30函数expand格式 说明：s为符号矩阵或表达式。常用在多项式的因式分解中，也常用于三角函数，指数例1-63 在编辑器中建立M文件symsxpq函数格式 %将s中的变量v的同幂项系数合 %s是矩阵或表达式，此命令对由命令findsym函数返回的默认变函数simple或 格式simple %s是矩阵或表达( 说明Simple(s)将表达式s格式 %使表达式s更加精1111abcd1111abcddd

在编辑器中建立M文件symsabcA=[1111;abcd;a^2b^2c^2d^2;a^4b^4c^4d^4]; %化简表达式 %让表达式d2符合人们的书写习惯d1d2=函数格式n= %计算矩阵A中元素的个1->>A=[1234;567>>n=numel(A)n=8Cholesky函数格式R=chol(X) 阵R，满足R'*R=X；若X非正定，则产生错误信息。[R,p]=chol(X) Xp为正整数，R是有序的上三角阵。1->>X%产生4阶pascal矩11 12 >>[R,p]=chol(X)R=1111012300130001p0LU矩阵的三角分解又称LUL和一个上三角矩阵U的乘积，即A=LU。函数格式[L,U]=lu(X) %U为上三角阵，L为下三角阵或其变换形式，满足LU=X。[L,U,P]=lu(X) %U为上三角阵，L为下三角阵，P为单位矩阵的行变换矩阵，1->>A=[123;456;780000= 000U000000= 000001UP1000 QR10函数格式[Q,R]=qr(A) %求得正交矩阵Q和上三角阵R，Q和R满足A=QR。[Q,R,E]=qr(A) %求得正交矩阵Q和上三角阵R，E为单位矩阵的变换形式，R的对角线元素按大小降序排列，满足AE=QR。[Q,R]=qr(A,0) %产生矩阵A的“经济大小”分解[Q,R,E]=qr(A,0) %E的作用是使得R的对角线元素降序，且Q*R=A(:,E)。R= %稀疏矩阵A的分解，只产生一个上三角阵R，满足R'*R[C,R]=qr(A,b) qr(A,b)，x=R\c。R %针对稀疏矩阵A[C,R 1->>A=[123;456;789;1011>>[Q,R]=------- - -0--00000R函数格式[Q,R]= %返回将矩阵A的第j列移去后的新矩阵的qr分1->>A=[-149-50-154;537180546;-27-9->>[Q,R]=qr(A)Q=-R557.9418187.0321000 %将A3列去掉后进行qr分解Q--- R557.9418 函数格式[Q,R]= %在矩阵A中第j列插入向量x后的新矩阵进行qrj大于AAx1->>A=[-149-50-154;537180546;-27-9->>x=[3510>>[Q,R]=qrinsert(Q,R,4,x)Q=-R-000Schur函数格式T=schur(A) %产生schur矩阵即T的主对角线元素为特征值的三角阵。T=schur(A,flag) %若A有复特征根，则flag='complex'，否则flag='real'。[U,T]=schur(A,…) %返回正交矩阵U和schur矩阵T，满足A=U*T*U'。1->>H=[-149-50-154;537180546;-27-9-25>>[U,T]=schur(H)U= - T 2.0000- 实Schur分解转化成复Schur函数格式[U,T]=rsf2csf 1->>A=[1113;1211;1131;-211>>[u,t]=schur(A)u=---t--0--0000--0.2756---0.2756-0.2133+---0.1012+-0.1046+-0.1842+-0.1867---0.2635-0.3134--0.9697+-0.5212+2.0051i-01.9202+ 0.1117+001.9202- 0.8002+00 T函数格式d= d=eig(A,B) %A、B为方阵，求广义特征值d，以向量形式存放d。[V,D]=eig(A) %计算A的特征值对角阵D和特征向量V使AV=VD成立。 %当矩阵A中有与截断误差数量级相差不远的值时[V,D]=eig(A,B)%计算广义特征值向量阵V和广义特征值阵D，满足[V,D]eig(A,B,flag)%由flagD和特征向量V，flag的A为对称Hermitian矩阵，B为正定阵。'qz'表示使用QZ算法，这里A、BHermitian矩阵。说明AxxAxBx函数格式s=svd %返回矩阵X的奇异值向[U,S,V]=svd(X) %XSUV，且满足U*S*V'Am×nUm×m阵，V[U,S,V]=svd( %得到一个“有效大小”的分解，只计算出矩阵U列，矩阵Sn×n1->>A=[12;34;56;7---000000 SV>>[U,S,V]=svd(A,0)U=S00V函数格式[U,V,X,C,S]= %返回酉矩阵U和V、一个普通方阵X、非负对角II为单位矩阵)；AB的列数必须相同，行数可[U,V,X,C,S]=gsvd(A,B,0) sigma=gsvd(A,B) 1- %产生34列矩阵，元素由1，2，…，12构成A >> %产生4B8 U= - V------0----000000000000000000000CS特征值问题的QZ函数格式[AA,BB,Q,Z,V]= %A、B为方阵，产生上三角阵AA和BB，阵。且满足：Q*A*Z=AA和Q*B*Z=BB。[AA,BB,Q,Z,V]= %flag决定的分解结果，flag阵是对称矩阵，则它的海森伯格形式是对角三角阵。可以通过相似变换将矩阵变函数格式H= %P为酉矩阵，满足：APHP'且P'Peye(size(A))1->>A=[-149-50-154;537180546;-27-9-P=000-0H- 42.2037--537.6783152.5511- r=n（n为方程组中未知变量的个数，则有唯一解；r<n，则可能有无穷解；=对应齐次方程组的通解+型稀疏矩阵——迭代法。5x1 x15x2 例1- 求方程组

x5x3x35x4x4

0>>A=[50000156001B=[1000 R_A5%满X或用函数rref>> %由系数矩阵和常数列构成增广矩阵>> %将C化成行最简R00000000000000000000x1x23x3x4例1- 求方程组3x1x23x34x44的一个特解x15x29x38x4>>A=[11-3-1;3-1-34;15-9->>B=[14 X=[00-0.53330.6000]’（一个特解近似值）。>>A=[11-3-1;3-1-34;15-9-B=[14>> >>R=rref(C)R=0- 0-- -00000 – LU、QR和choleskyLULUGauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则 所以 命令[L，U]=lu4x12x2x3例1- 求方程组3xx2x10的一个特解 A

220

11x13x2>>A=[42-1;3-12;113>>B=[210D0L - U - Warning:MatrixisclosetosingularorbadlyResultsmaybeinaccurate.RCOND=2.018587e->InD:\\pujun\lx0720matline4X=1.0e+016说明结果中的警告是由于系数行列式为零产生的。可以通过A*XCholeskyAR

R方程A*X=b

RR*XXR\(R\QR对于任何长方矩阵AQRQ为正交矩阵，R为上三角矩阵的初方程 变形 所 上例中[Q,说明这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘在中，函数null用来求解零空间，即满足A·X=0的解空间，实际上是求出格式z= %z的列向量为方程组的正交规范基，满足ZZI zAX=0x12x22x3x4例1- 求解方程组的通解：2x1x22x32x4x1x24x33x4>>A=[1221;21-2-2;1-1-4->>format B2-1001>>B=rref(A)B=0--00000symsk1 X[2*k1+5/3*k2][-2*k1-4/3*k2] %下面是其简化形[2k1+5/3k2 [-2k1- AX=b是否有解，若有解则进行第二步第二步：求AX=b的一个特解第三步：求AX=0第四步：AX=b的通解AX=0的通解+AX=bx12x23x3x4例1- 求解方程组3x1x25x33x42x1x22x32x4解：在中建立M文件如下A=[1-23-1;3-15-3;212-b=[12B=[Ab];formatratif elseifR_A==R_B&R_A<n %求AX=0的基础解系elseX='equitionnosolve' R_AR_B=X=

3equitionno说明

x1x23x3x4例1- 求解方程组的通解：3x1x23x34x4x15x29x38x4解法一：在编辑器中建立M文件如下A=[11-3-1;3-1-34;15-9-b=[14B=[Ab];formatratifR_A==R_B&R_A==nelseifR_A==R_B&R_A<nelseX='Equationhasnosolves'R_A2R_B2Warning:Rankdeficient,rank=2tol >InD:\\pujun\lx0723matline11X=00-C -

3/3/

3/4 7/4 所以原方程组的通解为X=k1 +k2

+ 1 8解法二：用rref

0 0

3/5A=[11-3-1;3-1-34;15-9-b=[14B=[A C10-01-0000310-01-00003/03/ 10

2

5/41/*

所以，原方程组的通解为：X=k11+k22+*00LQ函数格式xsymmlq(A,b)%求线性方程组AX=bX。An阶对称方阵，b为n元列向量。AafunA*X的函数。如果norm(b-A*x)/norm(b)和计算终止的迭代次数。 %maxit指定最大迭代次数 %M为用于对称正定矩阵的预处理因子 %x0为初始估计值，默认值为0。[x,flag]= 精度收敛；1表示在指定迭代次数内不收敛；2表示相同；4表示标量参数太小或太大；5表示预处理因[x,flag,relres]=symmlq(A,b,…) %relres表示相对误差norm(b-A*x)/norm(b)[x,flag,relres,iter]=symmlq(A,b,…) %iter表示计算x的迭代次数[x,flag,relres,iter,resvec]=symmlq(A,b,…) %resvec表示每次迭代的残差：[x,flag,relres,iter,resvec,resveccg]=symmlq(A,b,…) %resveccg表示每次迭代共函数x=bicg(A,b)%AX=bX。An阶方阵，bn元列向量。A可以是由afunA*X的函数。如果收敛，将 %maxit指定最大迭代次数 %M为用于对称正定矩阵的预处理因子 %x0为初始估计值，默认值为0。[x,flag]bicg(A,b,…)%flag的取值为：0表示在指定迭代次数之内按要求精度件的预处理因子；3表示两次连续迭代完全相同；4表示[x,flag,relres]=bicg(A,b,…) %relres表示相对误差norm(b-A*x)/norm(b)[x,flag,relres,iter]=bicg(A,b,…) %iter表示计算x的迭代次数[x,flag,relres,iter,resvec]=bicg(A,b,…) %resvec表示每次迭代的残差：例1- >>load %将数据取为系数矩阵Ab=sum %将A的各行求和，构成一列向量>> %用“\”求AX=b>>norm(b-A*X)/norm(b) ans=x,flag,relres,iter,resvec] %用bicg函数求解x（全为0，由于太长，不显示出来）flag= relres= %相对残差relresnorm(b-A*x)/norm(b)norm(b)/norm(b1。iter %表明解法不当，使得初始估计值0向量比后来所有迭代值0resvec=（略 >> xlabel('iteration %x轴为迭代次数ylabel('relative %yrelativerelative1-1双共轭梯度法相对误差函数格式x[x,flag]=bicgstab(A,b,…)[x,flag,relres]=bicgstab(A,b,…)[x,flag,relres,iter]=bicgstab(A,b,…)bicg一样。1-84>>load>>[L1,U1]=luinc(A,1e->>[x1,flag1]=bicgstab(A,b,1e-ALUL和UM=L*U，其结果是x10，flag=2表示预处理因子为坏条件的预处理因子。若改 %稀疏矩阵的完全LU分解 指定最大迭代次数为10次，预处理因子M=L*U。

relativerelative

1-2稳定双共轭梯度方法的相对误差x2=（1，略）flag2= relres2= iter2 resvec2 %每次迭代的残函数格式x[x,flag]=cgs(A,b,…)[x,flag,relres]=[x,flag,relres,iter]=cgs(A,b,…)[x,flag,relres,iter,resvec]=cgs(A,b,…)共轭梯度的LSQR函数格式x[x,flag]=lsqr(A,b,…)[x,flag,relres]=lsqr(A,b,…)[x,flag,relres,iter]=1->>n=>>on=>>Aspdiags([-2*on4*onon],- >>b=>>tol1e- >>maxit >>M1=spdiags([on/(-2)on],->>M2spdiags([4*on flag

x的值全为 relres3.5241e- iter 函数格式x[x,flag]=gmres(A,b,…)[x,flag,relres]=gmres(A,b,…)[x,flag,relres,iter]=gmres(A,b,…)[x,flag,relres,iter,resvec]=gmres(A,b,…)restart（An）可以给出外，调用方式和返回的结果形bicg一样。1->>load>>A=>>b=>>[x,flag]=结果显示flag=1，表示在默认精度和默认迭代次数下不收敛。>>[L1,U1]=luinc(A,1e->>[x1,flag1]=gmres(A,b,5,1e-[L2,U2]=luinc(A,1e-tol=1e-[x4,flag4,relres4,iter4,resvec4]=gmres(A,b,4,tol,5,L2,U2)[x6,flag6,relres6,iter6,resvec6]=gmres(A,b,6,tol,3,L2,U2)函数格式xminres(A,b,tolmaxit,M)[x,flag]=minres(A,b,…)[x,flag,relres]=minres(A,b,…)[x,flag,relres,iter]=minres(A,b,…)[x,flag,relres,iter,resvec]=minres(A,b,…)[x,flag,relres,iter,resvec,resveccg]=minres(A,b,…)bicg一样。1->>n=100;on=>>A=spdiags([-2*on4*on-2*on],->>b=>>tol=1e->>maxit=>>M1=flag0relresiter

resvec=（略）resveccg（略）函数格式x[x,flag]=pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)[x,flag,relres]=pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)[x,flag,relres,iter]=pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)bicg一样，这里A函数格式x[x,flag]=qmr(A,b,…)[x,flag,relres]=qmr(A,b,…)[x,flag,relres,iter]=qmr(A,b,…)[x,flag,relres,iter,resvec]=qmr(A,b,…)bicg一样，这里A1->>load>>A=>>b=>>[x,flag]=结果显示flag=1，表示在缺省情况下不收敛>>[L1,U1]=luinc(A,1e->>[x1,flag1]=qmr(A,b,1e-flag=2，表示是坏条件的预处理因子，不收敛。>>[L2,U2]=luinc(A,1e- x=（全为1，略）flag2= relres2= iter2= resvec2= 1.0e+005*

1-3准最小残差法相对误差Ann维列向量x使得关系式Axx成立，则称为方AxA”的特征向量。1.3.51.3.6例1- 求矩阵A

20 >>A=[-211;020;-4100=- 020002D即：特征值-1对应特征向量(- -特征值2对应特征向量(- -0.9701)T和(- - 例1- 求矩阵A

0>>A=[-110;-430;10000-200010001D说明当特征值为1(二重根)时，对应特征向量都是k -函数格式[T,B]= %求相似变换矩阵T和平衡矩阵B，满足BT1ATB 函数格式[V,D]=cdf2rdf 1->>A=[123;045;0-5v= -0.0191-0.4002i-0.0191+00-0+0d0004.0000+0 4.0000-V=-00D000000-命令格式B=orth(A) %将矩阵A正交规范化，B的列与A的列具有相同的空间，B的列向量是正交向量，且满足：B'*B=eye(rank(A))。例1- 将矩阵A

31 >>A=[400;031;01P000-0Q000000例1- 求一个正交变换X=PY，把二次f2x1x22x1x32x1x42x2x32x2x42x3x4化成标准形。 11A

1 1 10 0在编辑器中建立M文件如下A=[011-1;10-11;1-101;-111symsy1y2y3y4 %vpa表示可变精度计算，这里取2位精f=[y1y2y3-----00= 01000000001DX[.79*y1+.21*y2+.50*y3-.29*y4][.21*y1+.79*y2-.50*y3+.29*y4][.56*y1-.56*y2-.50*y3+.29*y4] f f=y2y23y2 AA中最高阶非零子式的阶数；向量组的秩通常由该向量组构成的矩函数格式k= %返回矩阵A的行（或列）向量中线性无关个k %tol例1-

1

3)，2

13)，3(1

4

3)，5

>>A=[1-223;-24-13;-1203;0623;2-63k33<

ri

（ij两行交换i行的K

kri行的Kj

rjkr将矩阵化成行最简形令是rref或rrefmovie。函数rref格式R=rref(A) %用高斯—约当消元法和行主元法求A的行最简行矩阵R[R,jb]=rref(A) A的列向量基；jb[R,jb] %tol 例1- >>-24-206-6a2A1--02-426-2-02333334>>[R,jb]=rref(A)R=jb

- - >>ans1-0-462-2333即 a4为向量组的一个基函数格式S=sparse(A) %将矩阵A转化为稀疏矩阵形式，即由A的非零元素和下标构成稀疏矩阵S。若A本身为稀疏矩阵，则返回A本身。S= %m×n0S=sparse(i,j,s)%i，js定义的稀疏矩阵S，其i，j长度相同的向量，表示在(i,j)位置上的元素。S= si，mmax(i)nmax(j)S= %m×nnzmaxS，nzmaxi1->>S=sparse(1:10,1:10,1:10)S=123456789>>S=sparse(1:10,1:10,5)S= 函数格式 %S为稀疏矩阵，A为满矩阵1->>S=sparse(1:5,1:5,4:8)S=45678>>A4000005000006000007000008函数格式k= [i,j %Xi[i,j,v] %Xij例1- 上例>>I=[134j=[134v=[4 函数格式 %D是只有3列或4列的矩说明：load函数把外部数据（.mat文件或.dat文件）装载于内存空间中的变量T；T数组的行维为nnznnz+1，列维为3（对实数而言）或列维为4（对复数而言；T数组的每一行（以[i,j,Sre,Sim]形式）指定一个稀疏矩阵元素。1->>D254;346;3123254346367>>S=spconvert(D)S=3647>>D=[1234;2540;3469;3674];D=1234254034693674>>S=spconvert(D)S= 3.0000+ 6.0000+ 7.0000+函数格式[B,d]= BdB %从AdBA=spdiags(B,d,A) %B中的列替换AdA=spdiags(B,d,m,n) %m×nAB中的列元素放d指定的对角线位置上。1->>A= 00000000000000>>[B,d]=spdiags(A)B=00d- %表示B的第1列元素在A中主对角线下方第3 %表示B的第2列在A %表示B的第3列在A的主对角线上方第21->>B=[12356791011131415>>d=[-201>>>>full(A)ans=270050函数格