教程-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=122334>>Null_M=[

第式>>>>C=[1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+1]5.4000+2>>R=[123;456],M=[111213;141516]R=123456MCN=1.00002.00003.00004.00005.00006.0000symsyms，先定义一>>sym_matrix=sym（'[abc；Jack，HelpMe!，NOWAY!]，'）sym_matrix= HelpMe! >>sym_digits=sym（'[123；abc；sin（x）cos（y）tan（z）]'）sym_digits= 1-4>>symsabc>>M1=>>M2=sym（'>>M3=>>syms_matrix=[abc；M1，M2，M3；int2str（[235]）]syms_matrix= [ClassicalJazzBlues] 5]供了一个将数值型转化成符号型令，即sym。>>Digit_Matrix=[1/3sqrt（2）3.4234；exp（0.23）log（29）23^（->>Syms_Matrix= 对于大型矩阵，一般创建M文件，以便于修改：例1-6 用M文件创建大矩阵，文件名为example.mexm=[285596在窗口输入 %显示exm5 %表示exm56列函数格式说明n=1n=2时分别构造[A1；A2]和[A1，A2]，都是二维数组，而n=3时可以构造>>A1=[1,2,3;4,5,6;7,8,9];A2=A1';A3=A1-A4(:,:,1)=123456789A4(:,:,2)147258369A4(:,:,3) - - - 1-8>>A1=[1,2,3;4,5,6;7,8,9];A2=A1';A3=A1-A5(:,:,1)= A5(:,:,2) A5(:,:,3) - - - 命令全零阵函数zeros格式B= B Bzeros([m B=zeros(d1,d2,d3…)%生成d1×d2×d3×…全零阵或数组B=zeros([d1d2d3…])%生成d1×d2×d3×…全零阵或数组B=zeros(size(A)) %生成与矩阵A相同大小的全零阵命令单位阵函数eye格式Y %生成n×n单位 %生成与矩阵A相同大小的单位阵命令全1阵函数格式Y= %生成n×n全1Y Yones([m Y=ones(d1,d2,d3…) %生成d1×d2×d3×…全1阵或数组Y=ones([d1d2d3…])%生成d1×d2×d3×…全1阵或数组Y= 命令均匀分布随机矩阵函数rand格式Y= %生成n×n随机矩阵，其元素在（0，1）Y Yrand([m Y= Yrand([mn Y %生成与矩阵A相同大小的 s %产生包括均匀发生器当前状态的35个元素', ', ', rand('state',sum(100*clock)) 例1-9 产生一个3×4随机矩阵R >>x=命令正态分布随机矩阵函数randn格式Y=randn(n) %生成n×n正态分布随机矩阵Y=randn(m,n) Y=randn([m Yrandn(m,n,p,…)%m×n×pY=randn([mnp…]) Y=randn(size(A)) %生成与矩阵A相同大小的正态分布随机矩阵 s %产生包括正态发生器当前状态的2个元素的srandn('state srandn('state srandn('state s=randn('state',sum(100*clock)) 例1-11 产生均值为0.6，方差为0.1的4阶矩阵>>mu=0.6;x=命令产生随机排列函数randperm格式p= %产生1～n之间整数的随机排ans= 命令产生线性等分向量函数linspace格式y= %在(a,b)上产生100个线性等分y=linspace(a,b,n) %在(ab)上产生n个线性等分点命令产生对数等分向量函数格式y= %在（10ay=logspace(a,pi)命令计算矩阵中元素个数

n=numel(a) 命令产生以输入元素为对角线元素的矩阵函数格式out=blkdiag(a,b,c,d,…) 例1-13>>out=blkdiag(1,2,3,4)out=1000020000300004命令友矩阵函数compan格式A= %u为多项式系统向量，A为友矩阵，A的第1行元素 求多项式(x1)(x2)(x3)x37x6的友矩阵和>>u=[10-7 A07-100010 ans-命令hadamard矩阵函数hadamard格式H=hadamard(n) %返回n阶hadamard矩阵例1-15h11111-1-11--1--1命令Hankel方阵函数hankel格式H= %第1列元素为c，反三角以下元素为0H=hankel(c,r) >>c r h=12382389389H= 产生一个3阶Hilbert矩>>formatH=1命令逆Hilbert矩阵函数invhilb格式H=invhilb(n) %产生n阶逆Hilbert矩阵命令Magic(魔方)矩阵函数格式M=magic(n) %产生n阶魔方矩阵例1-18M=816357492PascalA=%产生nPascalA=A=%返回Pascal(n,1)的转置和交A=11111234A=1001-01-1A=111--0100命令托普利兹矩阵函数格式TT=%用向量r生成一个对称的托普利>>c=[1234>>r=[1.52.53.54.5T=121321432154321命令Wilkinson特征值测试阵函数wilkinson格式W=wilkinson(n) %返回n阶Wilkinson特征值测试阵例1-21>>W1001100110>>01W31000001210000011100000101000001110000012100000131-22>>A=[1,1,1;1,2,3;1,3,>>B=[8,1,6;3,5,7;4,9,>>A-Ｂ=A-9274758- --- --- >>X=3412211100100]；8533则显示：a=468244C=C=相同；若为矩阵，则AB有相同的维数。例>>X=[-10>>Y=[-2-1则显示：Z44在中，用函数cross实现。函数Ccross(A,B)%A、BABC=A×B，A、B3A、B3×nAB对应列的叉积，A、B3×n矩C=cross(A,B,dim) 相同的维数，size(A,dim)和size(B,dim)必须是3。 计算垂直于向量(1,2,3)和(4,5,6)的向量>>a=[12>>b=[45- -可得垂直于向量(1,2,3)和(4,5,6)的向量为±(-3,6,3) 计算向量a=(1,2,3)、b=(4,5,6)和c=(-3,6,-3)的混合积a(b>>a=[123];b=[456];c=[-36->>x=dot(a,cross(b,结果显示：x函数conv格式w= 说明长度为munv的卷积(Convolution)kwk)∑u(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)= 展开多项式(s22s2)(s4)(s解：w %将w表示成多项Ps^4+7s^3+16s^2+18s+8函数deconv格式[q,r]= %多项式v除以多项式u，返回商多项式q和余多项式r (x32x23x4)(10x220x30)，则其卷积>>u= >>v= c= q=r 0函数格式C=kron %A为m×n矩阵，B为p×q矩阵，则C为mp×nq矩阵说明ABCAB

a22B

ABBMMOML1 A

B

AB >>A=[12;34];B=[123;456;78C=1232464568789集合运函数intersect格式c= %返回向量a、b的公共部分，即c=a∩bc [c,ia,ib] >>A=[1234;1246;6714]A=123412466714>>B=[1238;1146;6714]B=123811466714C= >>A=[19620];B=[1234610>>[c,ia,ib]=c16ia134ib157数ismember格式k= k= >>S=[024681012141618>>a=[12345k=

>>A=[1234;1246;671>>B=[1238;1146;671k=00%1数setdiff格式c= c= [c,i >>A=[179620];B=[1234610c= >>A=[1234;1246;671>>B=[1238;1146;671c=12341246数setxor格式c= c= [c,ia,ib] %ia、ib表示caA)、b(B)>>A=[123>>B=[245C= >>A=[1234;1246;6714]A=123412466714>>B=[1238;1146;6714]B=123811466714C=11 6ia12ib21函数union格式c= %返回a、b的并集，即c=a∪bc= [c,ia,ib]=union(…) 1-37>>A=[123>>B=[245>>c

>>A=[1234;1246]A= >>B=[1238;1146]B=12381146c=11246ia12ib21格式b=unique b=unique [b,i,j]=unique(…) %i、j体现b中元素在原向量（矩阵）中的位置1-39>>A=[112244646]A= c=ij

>>A=[1224;1146;1146]A=122411461146c= i31j211除法运x=b/a是方程x*a=b的解。例：a=[123;426;749]b=[4;1;-a\b=inv(a)*bb/a=矩阵乘小于0的整数时，A^P表示A-1P次方。dd当A为方阵，p为非整数时，则A^

V1VA的特征 为特征值对角矩阵。如果有重根，以上指令不成立标量的矩阵乘方PA，标量的矩阵乘方定义为PA V1式中V，Dpdn自特征值分解AV=AD

pa11

p标量的数组乘方P.^A，标量的数组乘方定义为P.^A 数组乘方p Lpam函数expm格式Y= %使用泰勒级数计算eA %使用特征值和特征向量计算命令矩阵的对数

Y=logm(X) %计算矩阵X的对数，它是expm(X)的反函数。[Y,esterr]=logm(X) >>A=[110;002;00-Y=000A=0000-F=是任意基本函数，如sincos等等，例如： 为结果所产生的相对误X=sqrtm(A)%AA1/2X*X=AXA的特征值XAX为复矩阵；若A为奇异矩阵，则X不存在。 [X,alpha,condest 命令矩阵A的多项式函数polyvalm格式polyvalm(P, 矩阵转函数格式ddet(XX1-42>>A=[123;456;78A D=0逆与伪命令函数格式 %求方阵X的逆矩阵。若X为奇异阵或近似奇异阵将给出警告信息1231-A221的逆矩343>>A=[123;221;34>>Y=inv(A)或Y=A^(-Y----1B3

321 3

0>>B=[1,2,3,1,0,0;2,2,1,0,1,0;3,4,3,0,0, >>X=C %取矩阵C中的A^(-1)部0000-00--00-----X>>A=[21-1;212;1-1>>formatD=

00--0命令伪逆函数格式B= Bpinv(A %tol种程度上代表矩阵的逆，若A为非奇异矩阵，则pinv(A)=inv(A)。 A=A(:,1:4) %54列元素构成矩阵A。A185746 %计算AX--------函数格式b=trace %返回矩阵A的迹，即A的对角线元和命令向量的范数函数norm∑|xk格式n= %X为向量，求欧几里德范数，即||∑|xknnorm(X,inf)%求-范数，即||X||maxabsXnnorm(X,1)%1-范数，即||X||1∑|xk|k)%p∑|xk|pn=norm(X,p)%求p-范数，即||Xp∑|xk|p命令矩阵的范数函数norm格式n=norm(A)%A为矩阵，求欧几里德范数||A||2A的最大奇异值。nnorm(A,1)%A的列范数||A||1A1-范数的最大值。n=norm(A,2)%A的欧几里德范数||A||2norm(A)相同。nnorm(A,inf)%求行范数||A||A1-即：max(sum(abs(A')))∑∑|aij n=norm(A,'fro' %求矩阵A的Frobenius∑∑|aij 即sqrt(sum(diag(A'*A)))，不能用矩阵p-范数的定义来求。命令范数的估计值函数格式nrm= 差小106nrm [nrm,count 条件命令矩阵的条件数函数cond格式c= %求X的2-范数的条件数即X的最大奇异值和最小奇异值的商c 义为：cond(A) |A||||A1||命令1-函数格式c=condest [c,v]=condest( %v为向量，满足||Av||||A||||v||c=norm(A,1)*norm(v,1)/c[c,v]condestA,t)%cv，同时显示出关于计算的步骤信息。如果t=1，则计算的每步都显示出来；如果t=-1，则给出商命令矩阵可逆的条件数估值函数rcond格式c=rcond(A) %对于差条件矩阵A来说，给出一个接近于0的数；对于好条件矩阵A，则给出一个接近于1的数。命令特征值的条件数函数condeig格式c= %DA的特征值对角阵，VA函数格式k=rank %求矩阵A的krank 函数diagXdiag(v,k)%vXkk=0时，vXk>0时，vkk<0时，v为下方第k条对角线。Xdiag(v)%v0Xvdiag(X,k)%Xkv。k=0：抽取主对角线元素；k>0k条对角线元素；k<0k条v=diag(X) v。1-46>>v=[12x=0000100002000030>>A=[123;456;789]A=123456789v=26函数tril 格式L= L=tril(X,k) 为主对角线以上；k<0为主对角线以下。函数 格式U= U=triu(X,k) 为主对角线以上；k<0为主对角线以下。 A111111111111 >>L=tril(A,1) L1100111011111111>>U=triu(A,-1) U=1111111101110011（1-48>A=[123456;678901]A=123456678901B=111111111111B1格式B=reshape(A,m,n) %返回以矩阵A的元素构成的m×n矩阵BB=reshape(A,m,n,p,…) %将矩阵A变维为m×n×p×…Breshape(A,[mn B >>b= 格式B=rot90 B=rot90(A,k) 1-50>>A=[123;456;789]A=123456789>>Y1=rot90(A),Y2=rot90(A,-Y1 Y2 格式B= 格式B= >>A=[123;456]A=12345632165321654= 123函数flipdim格式B= >>A=[123;456]A=12345645612456123= 654函数repmat格式B= %将矩阵Am×n块，即B由m×n块A平铺而成Brepmat(A,[m Brepmat(A,[mn >>A=[12;56]A= B=121212125656565612121212565656561212121256565656若关系满足，则将结果矩阵中该位置元素置为1，否则置0。1-><=~>>A=[1234;5678];B=[0214;077>>C1=A==B,C2=A>=B,C3=A~=BC1=01010010C211111011C310101101函数floor格式 函数ceil格式 函数round格式round %将A中元素按最近的整数取整，即四舍五入取整函数fix格式fix(A) %将A中元素按离0近的方向取整例1-55A--200-200-2-011=31103-122=200-2-122=20002-011格式[n,d]=rat n=--d

C=rem(A, (A./x)。允许模x为小数。 格式A&Band(A说明AB010格式A|Bor(A说明AB001格式~Anot格式xorA,B)说明AB0001。1-57>>A=[0234;1350],B=[1053;150023402341350= 1505BC1=00111100C211111111C310000001C4110000116.x4.2版中为符号矩阵设计的复杂函数形式，把符号矩阵的四则运算\）等或：符号矩阵的和（symadd、差（symsub、乘(symmul)1-58Asym([1x1/(x11/(x2),1/(x>>Bsym([x,1;x2,0])>>C=B-x-1/x1- （det（inv（rank（^）和指数（expexpm）等都与数值矩阵相同函数格式 %将A转化为符号矩阵A= B=[2/3, >>B(2,3) ans=4947709893870346*2^(->>B(2,3)='log(7)' B= 2/3,sqrt(2), 7/5,100/23,函数factor格式 说明S为符号矩阵或符号表达式，常用于多项式的因式分解。1-61x9-1分解因式在命令窗口键入symsx则显示：ans(x- 问“入”取何值时，齐次方程组2x1(3)x2x30有非0解 symsA=[1-k-24;23-k1;111-Dans=-k*(k-2)*(-函数expand格式 %符号表达式s的展开函说明：s为符号矩阵或表达式。常用在多项式的因式分解中，也常用于三角函数，指数例1-63将(x+1)3、sin(x+y)展开 symsxpq函数Collect格式 %函数simple或simplify 格式simple(s) %s是矩阵或表达式[R,how]=simple(s) %R为返回的最简形，how为简化过程中使用的主要方法。说明Simple(s)sPretty格式 %使表达式s更加精

1111abcdd1111abcdddsymsabcA=[1111;abcd;a^2b^2c^2d^2;a^4b^4c^4d^4]; %让表达式d2符合人们的书写习惯d1b*c^2*d^4-b*d^2*c^4-b^2*c*d^4+b^2*d*c^4+b^4*c*d^2-b^4*d*c^2-^2*d^4-a*b^2*c^4-a*b^4*d^2+a*b^4*c^2+a^2*c*d^4-a^2*d*c^4-a^2*b*d^4+a^2*b*c^4+a^2*b^4*d-d2=(-d+c)*(b-d)*(b-c)*(-d+a)*(a-c)*(a-(-d+c)(b-d)(b-c)(-d+a)(a-c)(a-函数格式n=numel(a) 例1-65>>A=[1234;567n=8函数格式R=chol(X) %如果X为n阶对称正定矩阵，则存在一个实的非奇异上三角阵R，满足R'*R=X；若X非正定，则产生错误信息。[R,p]=chol(X) 若X非正定，则p为正整数，R是有序的上三角阵。>>X%4阶pascal矩11 12 R=111101230013p00010LU上三角矩阵U的乘积，即A=LU。函数格式[L,U]=lu(X) %U为上三角阵，L为下三角阵或其变换形式，满足LU=X。[L,U,P]= LU=PX>>A=[123;456;78L=000U000000000000P001100010 QRU函数qr格式[Q,R]=qr(A) %求得正交矩阵Q和上三角阵R，Q和R满足A=QR。[Q,R,E]=qr(A) %求得正交矩阵Q和上三角阵R，E为单位矩阵的变换形式，[Q,R]=qr(A,0) %产生矩阵A的“经济大小”分解[Q,R,E]=qr(A,0) %E的作用是使得R的对角线元素降序，且Q*R=A(:,E)。R=qr(A) %稀疏矩阵A的分解，只产生一个上三角阵R，满足R'*R=[C,R]=qr(A,b) qr(A,b)，x=R\c。R [C,R]=qr(A,b,0) 1-68>>A=[123;456;789;1011-------------=- - -0--00000R函数格式[Q,R]=qrdelete(Q,R,j) 例1-69>>A=[-149-50-154;537180546;-27-9-Q=---R-000 %将A3列去掉后进行qr分解------=557.9418000R函数格式[Q,R]= %在矩阵A中第j列插入向量x后的新矩阵进行qr>>A=[-149-50-154;537180546;-27-9->>x=[3510Q=- - -R--0-00Schur函数格式T=schur(A) %产生schur矩阵T即T的主对角线元素为特征值的三角阵。T=schur(A,flag) %若A有复特征根，则flag='complex'，否则flag='real'。[U,T]=schur(A,…) %返回正交矩阵U和schur矩阵T，满足A=U*T*U'。>>H=[-149-50-154;537180546;-27-9-25U= - - T 2.0000- 函数格式[U,T]=rsf2csf(U,T) 例1-72>>A=[1113;1211;1131;-211u=--------t----0--0000--0.2756-0.2133--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 Td=d=%计算A的特征值对角阵D和特VAV=VD成立[V,Deig(A,'nobalance')%A该指令可能更精确。'nobalance'起误差调节[V,Deig(A,B)%计算广义特征值向量阵V和广义特征值阵D[V,Deig(A,B,flag)%flagDV，flag的可能值为：'chol表示对BCholesky分解算法，这里AHermitian矩阵，B为正定阵。'qzQZ算法，这里A、B为非对称或非Hermitian矩阵。AxxAxBx函数格式s=svd [U,S,VsvdX)%XSUV，[U,S,V]=svdX,0)%U>>A=[12;34;56;7U=-----------S000000V---U=-------S00V---函数格式[U,V,X,C,S]= II为单位矩阵)；AB的列数必须相同，行数可[U,V,X,C,S]=gsvd(A,B,0) sigma=gsvd(A,B) A %产生4阶魔方B8 ---= - --------=0------- 000000000000VXCS0 00 特征值QZ函数格式[AA,BB,Q,Z,V]= 阵。且满足：Q*A*Z=AAQ*B*ZBB。[AA,BB,Q,Z,V]= %flag决定的分解结果，flag'copx'（默认'ra'：如果矩H的第一子对角线下元素0H为海森伯格(Hessenberg)矩阵。如果矩阵是对称矩阵，则它的海森伯格形式是对角三角阵。可以通过相似变换将矩阵变换函数格式H= %P为酉矩阵，满足：APHP'P'Peye(size(A))>>A=[-149-50-154;537180546;-27-9-P= - H- 42.2037--537.6783152.5511- 若系数矩阵的秩r<n，则可能有无穷解；=求线性方程组的唯一解或特解（第一类问题5x1 x15x2 求方程 x25x3 0的解x35x4 x4 >>A=[50000156001B=[1000 R_A5X-- %由系数矩阵和常数列构成增广矩阵 %将C化成行最简R00000000-00000000-0000x1x23x3x4 求方程组3x1x23x34x44的一个特解x15x29x38x4>>A=[11-3-1;3-1-34;15-9->>B=[14 X=[00-0.53330.6000]’（一个特解近似值）。若用rref求解，则比较>>A=[11-3-1;3-1-34;15-9-B=[14 R=0-0---00000由此得解向量 – LU换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则 变成所以X=U\(L\b) 命令[L，U]=lu(A) 求方程

4x12x2x323xx2x10 11x13x2

A 11 >>A=[42-1;3-12;113>>B=[210D0--000=00-00UWarning:MatrixisclosetosingularorbadlyResultsmaybeinaccurate.RCOND=2.018587e-InD:\ \pujun\lx0720.matline4X=1.0e+016-即：ARR 其中R为上三角阵。方程A*X=b

RR*XXR\(R\QR方程 变形 所 上例中[Q 格式z= %z的列向量为方程组的正交规范基，满足ZZI x12x22x3x4 求解方程组的通解：2x1x22x32x4x1x24x33x4>>A=[1221;21-2-2;1-1-4->>format B2--1001B=0--00000symsk1 X[2*k1+5/3*k2][-2*k1-4/3*k2] %[2k1+5/3k2 [-2k1- 第二步：求AX=b的一个特解x12x23x3x4 求解方程组3x1x25x33x42x1x22x32x4解：在中建立M文件如下A=[1-23-1;3-15-3;212-b=[12formatratifR_A==R_B&R_A==n elseifR_A==R_B&R_A<n %求AX=0的基础解系elseX='equitionnosolve' R_AX=

3equitionno说明

x1x23x3x4 求解方程组的通解：3x1x23x34x4x15x29x38x4 A=[11-3-1;3-1-34;15-9-b=[14formatratifR_A==R_B&R_A==nelseX='Equationhasnosolves'R_A2R_B2Warning:Rankdeficient,rank=2tol InD:\\pujun\lx0723.matline11X=00-C-1001

3/3/

3/ 7/ 所以原方程组的通解为 8

3/A=[11-3-1;3-1-34;15-9-b=[14B=[A C10-01---00000

3/3/，10

2

3/7/01

非齐次方程组的特解为：5/1/* 所以，原方程组的通解为：X=k11+k22+*00线性方LQ函数xsymmlq(A,b)%AX=bX。An阶对称方阵，b为n元列向量。AafunA*X的函数。如果示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。 %指定误差tol，默认值是1e-6 %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元列 %指定误差tol，默认值是1e-6 %maxit指定最大迭代次数 %M为用于对称正定矩阵的预处理因子 %x0为初始估计值，默认值为0。[x,flag]bicg(A,b,…)%flag的取值为：0表示在指定迭代次数之内按要求精度收敛；1表示在指定迭代次数内不收敛；2表示M为坏条件的预处理因子；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表示每次迭代的残差：norm(b- 调用6.0数据文件west0479>>load %将数据取为系数矩阵Ab=sum %将A的各行求和，构成一列向量 %用“\AX=b>>norm(b-A*X)/norm(b) ans=1.2454e- %用bicg函数x=（0，由于太长，不显示出来）flag= relres= %相对残差relresnorm(b-A*x)/norm(b)norm(b)/norm(b)1。iter %表明解法不当，使得初始估计值0向量比后来所有迭代值0resvec=（略 >> xlabel('iteration %x轴为迭代次数ylabel('relative %y1-1双共轭梯度法相对误差函数格式x[x,flag]=bicgstab(A,b,…)[x,flag,relres]=bicgstab(A,b,…)1-84>>loadx1的值全为0，flag=2表示预处理因子为坏条件的预处理因子。>>[L2,U2]=luinc(A,1e- %稀疏矩阵的不完全LU>>[x2,flag2,relres2,iter2,resvec2]=bicgstab(A,b,1e-%指定最大迭代次数为10次，预处理因子M=L*U>>semilogy(0:0.5:iter2,resvec2/norm(b),'- %每次迭代的相对残差图形，见图1-2 结果 flag2 relativerelres2relative iter2 resvec2 10复共轭梯度平方法解方程 10函数格式x

iteration[x,flag]=cgs(A,b,…)[x,flag,relres]=cgs(A,b,…)

1-2稳定双共轭梯度方法的相对误差LSQR函数格式x[x,flag]=lsqr(A,b,…)[x,flag,relres]=lsqr(A,b,…)1-85>>n=>>on=- >>b=>>tol1e- >>maxit >>M1=spdiags([on/(-2)on],->>M2spdiags([4*on flag relres3.5241e- iter 函数格式x[x,flag]=gmres(A,b,…)[x,flag,relres]=gmres(A,b,…)[x,flag,relres,iter]=gmres(A,b,…)式与命令bicg一样。>>load>>A=>>b=>>[x,flag]=[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)[x8,flag8,relres8,iter8,resvec8]=gmres(A,b,8,tol,3,L2,U2)函数格式x[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一样。>>n=100;on=>>b=>>tol=1e->>maxit=>>M1=flag0relres=iter=

4.6537e-resvec=（略）函数格式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,…)函数格式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,…)1-88>>load>>A=>>b=>>[x,flag]=>>semilogy(0:iter2,resvec2/norm(b),'- flag2= relres2= iter2= resvec2= 1.0e+005*1-3准最小残差法相对误差A的特征值，非零向量xA”的特征向量。 求矩阵A 的特征值和特征向 >>A=[-211;020;-41V- - - - D-00020002 - - 求矩阵A 0的特征值和特征向量 >>A=[-110;-430;10V0-0--D200010001说明当特征值为1(二重根)时，对应特征向量都是k -函数格式[T,B]= %求相似变换矩阵T和平衡矩阵B，满足BT1ATB %求平衡矩阵函数格式[V,D]=cdf2rdf >>A=[123;045;0-5v= -0.0191-0.4002i-0.0191+00-0+0d0004.0000+0 4.0000-V=--00-00D0000-命令格式B=orth(A) %将矩阵A正交规范化，B的列与A的列具有相同的空间，B的列向量是正交向量，且满足：B'*B=eye(rank(A))。4 将矩阵A0

31正交规范 >>A=[400;031;01P000-0Q000000 求一个正交变换X=PY，把二次0 1 11 1 0f2x1x22x1x32x0 1 11 1 0AA=[011-1;10-11;1-101;-111symsy1y2y3y4 f=[y1y2y3-----00=1000010000-00001DX[.79*y1+.21*y2+.50*y3-.29*y4][.21*y1+.79*y2-.50*y3+.29*y4][.56*y1-.56*y2-.50*y3+.29*y4] fy1^2+y2^2- f=y2y23y2 函数格式k= k 求向量组1(1 3)，2 13)，3(1 3)4 3)5 4)的秩，并判断其线性相关性。>>A=[1-223;-24-13;-1203;0623;2-63k3<交换两行rirj（i、第j两行交换i行的K倍kri行的Kj行上去rjkr函数rrefrrefmovie格式R=rref(A) %用高斯—约当消元法和行主元法求A的行最简行矩阵R[R,jb]= [R,jb] %给出每一步化简 >>-2>>a2=[-4- >>a3=[-2062-6a2A1--02-426-2-02333334R=jb

- - >>ans1-0-462-2333即 a4为向量组的一个基函数格式S=sparse(A) %将矩阵A转化为稀疏矩阵形式，即由A的非零元素和下标构成稀疏矩阵S。若A本身为稀疏矩阵，则返回A本身。S S=sparse(i,j,s) 量或与i，j长度相同的向量，表示在(i,j)位置上的元素。S= si，mmax(i)nmax(j)S= S，nzmaxiS=123456789S= 函数格式A=full(S) %S为稀疏矩阵，A为满矩阵。例1-97S=4567>>8A4000005000006000007000008函数格式k= [i,j [i,j,v] >>I=[234j=[234v567函数格式 %D是只有3列或4列的矩load函数把外部数据（.mat文件或.dat文件）装载于内存空间中的变T；T数组的行维为nnznnz+1，列维为3（对实数而言）或列维为4（对复数而言）T数组的每一行（以[i,j,Sre,Sim]形式）指定一个稀疏矩阵元素。>>D254;346;3123254346367S=3647>>D=[1234;2540;3469;3674];D=1234254034693674S= 3.0000+ 6.0000+ 7.0000+函数spdiags格式[B,d]= B A=spdiags(B,d,A) A=spdiags(B,d,m,n) 在由d指定的对角线位置上。>>A= 000000000000>>[B,d]=spdiags(A)B=00d- %表示B1列元素在A中主对角线下3条对角线 %表示B2A >>B=[12356791011131415>>d=[-201>>ans=270050数speye格式S= S %生成n×n的单位函数sprand格式R= R R=sprand(m,n,density,rc) 数sprandn格式R= R= R= %生成一个近似的条件数为1/rc、大小为m×n函数sprandsym格式R= %生成稀疏对称随机矩阵，其下三角和对角线与S具有相R= %生成n×nR %生成近似条件数为1/rc的稀疏对称R 表示矩阵由一正定对角矩阵经随机Jacobi旋阵为外积的换位和，其条件数近似等于1/rc；kind=3表示生成一个与矩阵S结构相同的稀疏随机矩阵，条件数近似为1/rc，density被函数nnz格式n= 数nonzeros格式s= >>A=[ >>s=[2176511101615函数nzmax格式n=nzmax(S) 函数格式S= %产生一个m×n阶只有nzmax个非零元素的稀疏矩阵函数spfun格式f= %用S中非零元素对函数'function'求值，如果'function'例1-103 4阶稀疏矩阵对角矩阵S f=f

%即指数e的非零元素函数格式R= %将稀疏矩阵S中的非零元素全换为函数格式 %参数与上面相同000>>load>>函数colamd格式p= %返回稀疏矩阵S的列 nz=1-4稀疏矩阵函数colmmd格式p=colmmd(S)%S的列的最小度排序向量p，按p排列后的矩阵为S(:,p)。 比较稀疏矩阵S与排序后的矩阵S(:,>>load>>>>>>>>>>>>1-5稀疏矩阵的排序函数colperm格式j=colperm(S)%S的列变换的向量。列按非0元素升序排列。有时这是LU分解前有用的变换：LU(S(:,j))。如果S是一个对称矩阵，对行和列进行排序，有利于Cholesky分解：函数dmperm格式p=dmperm A(p,:)是具有非0对角线元素的[p,q,r]=dmperm(A) %A为方阵p为行排列向量q为列排列向量，使得A(p,q)是上三角块形式，r为索引向量。[p,q,r,s] >>A=[110130;410044;022024;00630]A=000000000p= q r >>ans000000000数randperm格式p=randperm %对正整数1，2，3，…，n的随机排列，可以用来创建随p= 函数symamd格式p= 函数symrcm格式r=symrcm(S) %返回S的对称逆Cuthill-McKee排序r，使S的非0元素集中函数symmmd格式p=symmmd(S) %返回S的对称最小度排列向量p，S为对称正定矩阵。例1-108>>B=>>r=>>p=>>R=>>S=0 nz=2400 nz=

0 nz=2400 nz=1-6稀疏对称最小度命令矩阵A1-范数估计值函数condest格式c= [c,v]= %A1-cv||Av||(||A||||v||c，即命令2-函数格式nrm=normest(S) nrm=normest(S,tol) %tol为指定的相对误差，而不用默认误差10-6。[nrm,count]=normest(…) %count为给出的计算范数迭代的次数LU分解函数luinc格式[L,U]= 三角矩阵；'0'是一种分[L,U,P [L,U]= 掉的列元素。ugiag1表示droptol值代替上三角0。thresh为[L,U [L,U,P]=luinc(X,options)[L,U,P]=>>S=[110130;410044;022024;00630]S=000000000S= 1111LUpCholesky分解函数cholinc格式R= droptolR 果michol=1，则从对角线上抽取出被去掉的元素。rdiagdroptol值代替上三角分解因子中的对角线上的R %'0'是一种分解[R,p]= %不产生任何出错信息，如果R存在，则p=0；如果R d=d=%求稀疏矩阵的广义特征值问题。满足AV=BVD，其中d d deigs(A,k,sigma)%sigma取值：'lm'表示最大数量的特征值；'sm'最小数量特征值；对非对称和复数问题：'lr'表示最大实部；'sr表示最小实部；'li'表示最大虚部；'si'表示最小虚部。d %d d d d=d=d=d=d=d= %D6个最大特征值对角阵，V的列向量为对应特征向量[V,D,flag]= 2章数值计算与数据分析三角函数与双曲函函数sin、功能正弦函数与双曲正弦函数格式Y=sin(X) 度分量的正弦值Y，所有分量的角度单位为弧度。Y Z=x+iy，函数的定义为：sin(x+iy)=sin(x)*cos(y)sin(z) sin(z) eizsin(z) sin(z)2

ezezx=-pi:0.01:pi;x=-5:0.01:5;2-1正弦函数与双曲正弦函数函数asin、格式Yasin(X)%X（可以是向量、矩阵）中每一个元素的反正弦函数值YX中有的分量处于[-1,1]Y=asin(X)对应的分量处于[-π/2,π/2]之间，若X中有分量在区间[-1,1]之外，则Y=Yasinh(X)%X说明反正弦函数与反双曲正弦函数的定义为：asinziln(iz1z2asinhzln(z

1z2x=-1:.01:1;x=-5:.01:5;函数cos、

2-2反正弦函数与反双曲正弦函数功能余弦函数与双曲余弦函数格式Ycos(X)%X（可以是向量、矩阵，元素可以是复数）中每一个角cos(pi/2)并不是精确的零，而是与浮点精度有关的无穷小量eps，因为pi仅仅是精确值π浮点近似的表示值而已。Y Xzx+iy，则函数定义为：cos(x+iycos(x)*cos(ycosz

eizeiz

，coshz

ezezx=-pi:0.01:pi;x=-5:0.01:5;2-3余弦函数与双曲余弦函数函数acos、格式Y= %返回参 X（可以是向量、矩阵）中每一个元素的反余弦函YX中有的分量处于[-1,1]之间Y=acos(X)对应的分Yasinh(X)%X说 反余弦函数与反双曲余弦函数定义为：acosziln(izi1z2)acoshzln(z

z2x=-1:.01:1;x=-5:.01:5;2-4反余弦函数与反双曲余弦函数函数tan、功能正切函数与双曲正切函数格式Y=tan(X) %计算参量X（可以是向量、矩阵，元素可以是复数）中每一个角tan(pi/2)并不是精确的零，而是与浮点精度有关的无穷小量eps，因为pi仅仅是精确值π浮点近似的表示值而已。Y xpi/2)+0.01:0.01:(pi/2)-0.01;%稍微缩小定义x=-5:0.01:5;2-5正切函数与双曲正切函数函数atan、格式Y=atan(X)%X（可以是向量、矩阵）中每一个元素的反正切函数Y。若X中有的分量Y=atan(X)对应的分量处于[-Yatanh(X)%XYatanzilniz，atanhz1ln1 i 1x=-20:0.01:20;x=-0.99:0.01:0.99;2-6反正切函数与反双曲正切函数函数cot、功能余切函数与双曲余切函数格式Y=cot(X) 度分量的余切值Y，所有角度分量的单位为弧度。Y x1pi+0.01:0.01:-0.01;%去掉奇点x0x2=0.01:0.01:pi-0.01;%做法同上2-7余切函数与双曲余切函数函数acot、格式Y=acot(X) Y=acoth(X)%返回参量X中每一个元素的反双曲余切函数值Yx12*pi:pi/30:-0.1;x20.1:pi/30:2*pi;%去掉奇异点x0x1=-30:0.1:-1.1;x2=2-8反余切函数与反双曲余切函数函数sec、功能正割函数与双曲正割函数格式Y=scX)%X（）度分量的正割函数值Y，所有角度分量的单位为弧度。我们要的是，ec(p2)并不是无穷大，而是与浮点精度有关的无穷小量eppi仅仅是精确值πYsech(X)%Xx1pi/2+0.01:0.01:pi/2-0.01;%去掉奇异点xpi/2x2=pi/2+0.01:0.01:(3*pi/2)-0.01;x=-2*pi:0.01:2*pi;2-9正割函数与双曲正割函数函数asec、格式Y=asec(X) Y x1=-5:0.01:-1;x2=x=0.01:0.001:1;2-10反正割函数与反双曲正割函数csc、功能余割函数与双曲余割函数格式Y=csc(X) 角度分量的余割函数值Y，所有角度分量的单位为弧度。Y x1pi+0.01:0.01:-0.01;x20.01:0.01:pi-0.01;%去掉奇异点函数acsc、

2-11余割函数与双曲余割函数格式Y=asec(X) 值YY x110:0.01:-1.01;x21.01:0.01:10;%去掉奇异点xx1=-20:0.01:-1;x2=2-12反余割函数与反双曲余割函数功能四象限的反正切函数P=atan2(Y,X)%XYXY元素的实数部分对应的、元素对元素的四象限的反正切函数阵列P，其中X和Y的虚数部分将忽略。阵列P中的元素分布在闭区间[-pi,pi]上。特定的象限将取决于sign(Y)与sign(X)。yy-xr=z=r*exp(i*theta)feather(z);holdonaxisequal;holdthetaz1.0000+2-13四象限的反正切函数函数功能格式Bfix(A)%对AA同维的数组。A，则返回一复数，其分量的实数与虚数部分分别取>>A=[-1.9,-0.2, ,5.6,7.0,>>B=BColumns1through- Columns5through 2.0000+函数功能朝最近的方向取整。格式Y=roun(X)%对XX的数组。对于复数参量X，则返回一复数，其分量的实数与虚数>>A=[-1.9,-0.2, ,5.6,7.0,>>Y=YColumns1through- Columns5through 2.0000+函数功能朝负无穷大方向取整格式B=floor(A) %对A的每一个元素朝负无穷大的方向取整数部分，返回与A同维的数组。对于复数参 A，则返回一复数，其分量的实数与

>>A=[-1.9,-0.2, ,5.6,7.0,>>F=FColumns1through- - Columns5through 2.0000+函数功能求作除法后的剩余数Rrem(X,Y)%Xfix(X./Y).*YX、YX、Y为是不可意料的。fix(X./Y)为商数X./Y朝零方向取的整数部分。若X与Yrem(X,Y)返回的结果与mod(X,Y)相同，Xrem(-X,Ymod(-X,Y)-Y。该命令返回的结果在区间[0，sign(X)*abs(Y)]，若Y中有量，则相应地NaN。

>>X=[122334>>Y=[372>>R=R函数

功能朝正无穷大方向取整格式B=floor(A)%A的每一个元素朝正无穷大的方向取整数部分，返回与A同维的数组。对于复数参量A，则返回一复数，其分量的实数与虚>>A=[-1.9,-0.2, ,5.6,7.0,>>B=BColumns1through- Columns5through 3.0000+函数功能以e格式Y=exp(X) %对参量X的每一分量，求以e为底数的指数函数Y。X中的分量可以为复数。对于复数分量如，z=x+i*y，则相应地计算：e^ze^x*(cos(y)i*sin(y))>>A=[-1.9,-0.2, ,5.6,7.0,>>Y=Y1.0e+003Columns1through Columns5through -0.0099-函数功能求矩阵的以e格式Y=expm(X) %计算以e为底数、x的每一个元素为指数的指数函数值。若矩阵x有小于等于零的特征值，则返回复数的结果。Taylorexpm2.m中。但这种一般计在文件expm3.m中，先是将矩阵对角线化，再把函数计算出相应的的特征向量，>>Y=Y功能自然对数，即以e格式Y=log(X)%XX中的元素可以log对复数的计算如下：log(z)=log(abs(z))+i*atan2(y,x)2-21下面的语句可以得到无理数π>>Pi=abs(log(-Pi函数功能常用对数，即以10格式Y=log10(X) %计算X中的每一个元素的常用对数，若X中出现复数，则可>>L1log10(realmax)%由此可得特殊变量realmax的近>>L2log10(eps)%由此可得特殊变量eps的近似>>M=>>L3=L10

-函数格式B=sort(A)%AA中的元素。A可以是字符串的、实数的、复数的单元数组。对于A中完AA为则再按它们在区间[-π,π]的幅角从小到大排列；若A中有元素为若A为二维矩阵则按列的方向进行排列若A为数组sort(A)B= %沿着矩阵A（向量的、矩阵的或的）中指定维[B,INDEXsort(A,…)%输出参量BINDEXsize(A)A中列向量的A中有重复出现的相同的>>A=[-1.9,-0.2, ,5.6,7.0,>>[B1,INDEX]=>>M=>>B2=B1Columns1through- - 2.4000+Columns5through INDEX B2 函数

功能数值的绝对值与复数的幅值格式Y=abs(X) %返回参量X的每一个分量的绝对值；若X为复数的，则返回每一分量的幅值：abs(X)=sqrt(real(X).^2+imag(X).^2)。>>A=[-1.9,-0.2, ,5.6,7.0,>>Y=Y 函数功能格式ZC= conj(Z)=real(Z)-函数功能格式Y=imag(Z) %返回输入参量Z的每一个分量的虚数部分。例2-25ans3函数功能格式Y=real(Z) %返回输入参量Z的每一个分量的实数部分。例2-26ans2函数功能格式P=angle(Z) %返回输入参量Z的每一复数元素的、单位为弧度的相角，其值说明angle(zimaglog(zatan2>>Z=[1-i,2+i,3-i,>>P=P--------函数功能用实数与虚数部分创建复数格式c= （同为向量、矩阵、或阵列。该命令比下列形式的复数输入更有用：ai*baj*b因为ij可能被用做其他的变量(不等于sqrt(-1))，或者ab不是双精度的。c=complex(a) 2-28>>a=>>b=>>c=c1.0000+2.0000+3.0000+4.0000+函数功能模数（带符号的除法余数用法Mmod(X,Y)%X、Y应为整数，此时返回余数X-Y.*floor(X./Y)，等于rem(X,Y。总之，对于整数x,ymod(-x,y)=>>M1=>>M2=>>M3=M1M2=M3

函数n<15函数C= %参量n,k为非负整数，返回n!/((n-k)!k!)，即一次从nkC=nchoosek(v,k) kn-kk!)k>>C=C2468246248268468函数功能生成元素均匀分布于(0,1)用法Y=rand(n) Y=rand(m,n)、Y=rand([m 于区间(0,1)上矩阵YY=rand(m,n,p,…)、Y=rand([mn Y=rand(size(A)) %生成一与阵列A同型的随机均匀阵列Y s= 2-命命令 设置状态为设设第k个状态(k为整数在每次使用时的状态都不同（因为clock每次都不同>>R1=>>a=10;b=>>R2ab-a*rand(5)%生成元素均匀分布于(10,50)R1R2函数功能生成元素服从正态分布（N(0,1)）格式Y=randn(n) %返回n*n阶的方阵Y，其元素服从正态分布N(0,1)。若n不是Y=randn(m,n)、Y=randn([m (0,1)YY=randn(m,n,p,…)、Y=randn([mn Y=randn(size(A)) %生成一与阵列A同型的随机正态阵列Y s=randn('state') 前状态。该改变的当前状态，见表2-2。2- 设置状态为设设第k个状态(k为整数设在每次使用时的状态都不同（因为clock每次都不同>>R1=>>R2=0.6+sqrt(0.1)*R2yyf(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的yφ()y=φ()x=()在x点的值。寻找这样的函数()，办法是很多的。()可以是一个代数多项式，或是三角多项式，也φx)（任意阶导数连续的近效果。许多应用中，常要用一个解析函（一、二元函插值命命令 功能一维数据插值（表格查找f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图为图2-14Y：原始数据点Yi：插值2-14数据点与插值点关系示意格式yi= %返回插值向量yi，每一元素对应于参量xi，同时若Y为一矩阵，则按Y的每列计算。yi是阶数为yi yi=interp1(x,Y,xi,method) 值：’nearest’：最近邻点插值，直接完成计 ，数。命令spline用它们执行三次样条函数插值；’pchip’：分段三次Hermite插值。对于该方法，命令interp1调用函数pchip，用xy执行分段三次内插值。该方法保留单调性与’cubic’：与’pchip’操作相 yi 殊的外插值法extrapyi= >>x=0:10;y=>>xx=0:.25:10;yy=插值图形为图2-15。2-32>>year=>>product=[75.99591.972105.711123.203131.669150.697179.323203.212249.633256.344267.893>>p1995=>>x=>>y=p1995图2-15一元函数插值图 图2-16离散数据的一维插值命令 功能二维数据内插值（表格查找格式ZI= %返回矩阵ZI，其元素包含对应于参量XI与YI（以是向量、或同型矩阵）的元素，即Zi(i,j←[Xi(i,j),yi(i,j)]Xi与Yi，此时，输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩X、YZ确定的二维函数Z=f(X,Y)。参XY必须是单调的，且相同的meshgridXiYi中有在X与Y范围之外的点，则相应地nanNumberZI= %缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再ZI= 等价于interp2(z,1)。ZI 法；’nearest’：最插值；>>[X,Y]=meshgrid(->>Z=>>ZZ=>>axis([-33-33-520]);shading>>hold2-34>>years=>>service=>>wage=[150.697199.592179.323195.072203.212179.092226.505153.706249.633120.281>>w=w

2-17二维插值命令 功能三维数据插值（查表格式VI=interp3(X,Y,Z,V,XI,YI,ZI) 在点（XI,YI,ZI）的值。参量XI,YI,ZI是同型阵列或向量。若向量参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令(X,Y,Z)之外的点，则相应地返回特殊变量值NaN。VI= %缺省地X=1:NY=1:MZ=1:P，其中[M,N,P]=size(V)，再按上面的情形计算。VI= interp3(V)等价于interp3(V,1)VI >>[x,y,z,v]=>>vv=2-18三维插值命令 功