第2章 矩阵与MATLAB.doc

第2章 矩阵与MATLAB2.1 概述MATLAB是一门语言，它的运算指令和语法基于一系列基本的矩阵运算以及它们的扩展运算。2.2 矩阵和向量矩阵a的大小可用如下命令获得:size(a)-size可显示出两个值:第一个值为行数(m),第二个值为列数(n)。m,n=size(a)-m,n表示矩阵的行赋给m，列数赋给n。2.2.1 方阵当m=n时,a称为方阵。2.2.2 对角阵当aij=0,ij,且m=n时,得对角阵:(nn)如果a=1 3 8,那么,a=1 3 8；z=diag(a)产生如下对角阵：z= 1 0 0 0 3 0 0 0 8如果b是一个方阵，那么，diag(b)-返回b的对角阵元素。单位矩阵I的MATLAB表达式为：eye(n)-其中n是方阵的大小。c=eye(3)产生：c= 1 0 0 0 1 0 0 0 12.2.3 列矩阵和行矩阵(向量)当aij=ai1(即只有1列时),a称为列矩阵(列向量),记做:当aij=a1j时,即矩阵只有1行,则矩阵称为行矩阵或行向量,记做:a=a11 a12 a1n=a1 a2 an在MATLAB中,这是向量的默认定义。2.2.4 矩阵和向量的转置矩阵和向量的转置用()表示。向量a的长度，即元素的个数：L=length(a) 或 L=size(a)矩阵的大小只能用命令size获得。使用size的好处是不需要事先知道a是向量还是矩阵。但对于向量来说，length命令更方便。2.3 向量的创建在MATLAB中给出了几种方法用来对向量和矩阵中的元素进行赋值。第一种方法是使用冒号来指定数值范围和相邻值的步长；第二种方法是指定数值的范围和期望的数值个数。在第一种方法中，步长最为重要。在第二种方法中，数值的个数最为重要。第一种方法： x=s:d:f其中，s=初值，d=增量或减量值，f=结束值或终值。可产生如下的行向量x:x =s+d,s+2d,s+nd其中，s+nd =f 。当d省略时，MATLAB默认d=1。即：x=s:f 产生向量： x =s+1,s+2,s+n第二种方法：（指定n为从s到f的等间隔值）s=linspace(s,f,n)其中，增量或减量值d由MATLAB通过下式计算得出：如果没有指定n的值，MATLAB则默认n=100。如果需要用对数等间隔，则：x=logspace(s:f:n)产生向量：x =10s,10s+d,10s+2d,10f如果没有指定n的值，MATLAB则默认n=50。例1：若想创建-2,1,3,5,7,9,10的向量x。x=-2,1:2:9,10或 x=-2,1,3,5,7,9,10若 z=x-1 z=-3,0,2,4,6,8,9 z=x；z(2)=z(2)/2；（括号里2代表某列的值） z=-2,0.5,3,5,7,9,10 z=x；z(2)=z(2)/2；z(3:4)=z(3:4)*3-1；（第三个和第四个数乘三减一） z=-2,0.5,8,14,7,9,10 例2：如果， y=-1,6,15,-7,31,2,-4,-5 x=y(3:5) x=15,-7,31 x=y(1),y(2),y(7),y(8)或x=y(1,2,7,8) 或 index=1 ,2,7,8；x=y(index) x=-1,6,-4,-5例3：如果y=-1,6,15,-7,31,2,-4,-5；z=10,20,30,40,50,60,70,80用sort函数把向量y按升序排列，然后根据向量y中的新元素顺序重新排列向量z。函数sort的调用格式为： ynew,index=sort(y)其中，ynew是重新排列后的向量；index是一个由元素位置所构成的向量，它表示排列后向量y中元素的原始位置号。y=-1,6,15,-7,31,2,-4,-5；z=10,20,30,40,50,60,70,80；ynew,index=sort(y),znew=z(index)执行结果为：ynew -7 -5 -4 -1 2 6 15 31index 4 8 7 1 6 2 3 5znew 40 80 70 10 60 20 30 50可以通过find命令来扩展这种功能。find命令可找到满足用户指定条件或表达式的向量（或矩阵）的元素的位置（不是元素值）。例4：使用find命令在y中选出非正元素，组成一个新的向量s。y=-1,6,15,-7,31,2,-4,-5；index=find(y=0)；s=y(index)s -1 -7 -4 -5MATLAB中隐含向量和矩阵说明的一个最大好处就是可帮助用户实现对于一组织的一系列操作。例5：确定sin(x)的值，其中x在-x的区间内以/5的步长变化。x=-pi:pi/5:pi；y=sin(x)y 0.0000 -0.5878 -0.9511 -0.9511 -0.5878 0.0000 0.5878 0.9511 0.5878 0.00002.4 矩阵的创建矩阵可以通过以下几种方法创建。（下列矩阵是什么形式的复制不了？）第一种方法：第二种方法：（直接创建）第三种方法：（更为形象）= 11 12 13;21 22 23; 31 32 33;41 42 43第四种方法：(右表的下标是怎么写上的？)= 11 12 1321 22 23 31 32 33 41 42 43两个非常有用的函数可以生成矩阵的元素。1、创建一个元素为1的(rc)阶矩阵： one=ones(r,c)2、创建一个元素为0的(rc)阶矩阵： zero=zeros(r,c)例6：试创建如下矩阵（不会）（35） = 3:2:11；linspace(20,21,5)；ones(1,5) (1,1) 3 , (3,4) 1 , (:,2) 5 20.25 1(2,:) 20 20.25 20.5 20.75 21(1:3,3:5) 7 9 11；20.5 20.75 21；1 1 1例7: 创建一个大小与矩阵相同且元素都等于4的矩阵：= 3:2:11；linspace(20,21,5)；ones(1,5)；z = 4*ones(size(a)下面介绍一个特殊函数：magic(n)-魔方矩阵，矩阵中任意行或任意列中元素之和，以及两条对角线上的元素之和都相等。例8：z=magic(4)结果显示：16 2 3 135 11 10 89 7 6 124 14 15 1例9：z=magic(4)； z(2,:)=z(2,:)/2； z(:,4)=z(:,4)+z(:,2)例10：z=magic(4)； z=z-diag(diag(z) % 将矩阵z的对角线元素置为0例11：z=magic(4)； z=z-diag(diag(z)+5*eye(4) % 将矩阵z的对角线元素置为5例12：z=magic(4)； z=z-diag(diag(z)+diag(11 23 54 61)将矩阵z的对角线元素分别置为11、23、54、61。MATLAB提供了两个创建矩阵的函数，它们是通过复制标量、列向量、行向量和块矩阵来创建矩阵的。1、函数repmat：调用格式为： s=repmat(x,r,c)其中，x 可以是标量、向量或矩阵；r 是x的行数；c是x的列数。函数repmat可创建任意长度的列向量和行向量，而且每一个元素都具有相同的值。例13：w = repmat(45.72,1,6) w = 45.72, 45.72, 45.72, 45.72, 45.72, 45.72 或 w(1,1:6)=45.72（？再看看怎么回事） w = repmat(45.72,3,3) 或 w(1:3,1:3)=45.72（？再看看怎么回事）例14：考虑向量 v=repmat(s,3,1) v=repmat(s,3,2) v=repmat(,1,3) v=repmat(,2,3)（转置矩阵怎么写？）2、函数meshgrid:调用格式为： u,v=meshgrid(s,t)其中，s、t为两个行向量，u、v为调用meshgrid后生成的新矩阵。例15：s=1,2,3,4 , t=5,6,7u,v=meshgrid(s,t)而w=meshgrid(s,t)仅返回一个矩阵，如： w=meshgrid(s,t) w=uMATLAB还提供了两个矩阵操作函数fliplr(a)和flipud(a)。它们可实现对矩阵行和列的旋转操作。例16：考虑矩阵（fliplr(a)表示行不变，但是列左右变）（flipud(a)表示每行的左右顺序不变，但是行之间上下颠倒）a=1,2,3,4,5 6,7,8,9,10 fliplr(a)= 5,4,3,2,1 10,9,8,7,6flipud(a)=6,7,8,9,10 1,2,3,4,5flipud(fliplr(a)=10,9,8,7,6 5,4,3,2,1c=a,fliplr(a)=1,2,3,4, 5, 5,4,3,2,1 6,7,8,9,10,10,9,8,7,6如果想去掉其中一行，可以通过表达式使其中一行为空，其中方括号中没有空格用表达式：c(length(a),:)= 1,2,3,4, 5, 4,3,2,1 6,7,8,9, 10,9,8,7,6下面通过MATLAB三种不同操作的结果来进一步说明上述用法。例17：首先创建两个（25）阶矩阵a和b。a=1,2,3,4,5 6,7,8,9,10 b=11,12,13,14,1516,17,18,19,20考虑它们在三种MATLAB操作中的用法：相加：c=a+bc= 1+11, 2+12, 3+13, 4+14, 5+15 6+16, 7+17, 8+18, 9+19,10+20 （25）列扩展：c=a,bc= 1, 2, 3, 4, 5, 11, 12, 13, 14, 15 6, 7, 8, 9,10, 16, 17, 18, 19, 20 （210）行扩展：c=a;bc= 1, 2, 3, 4, 5 6, 7, 8, 9,10 11,12,13,14,15 16,17,18,19,20 （45）更进一步，如果 x=1,2,3 y=4,5,6那么，z=x,y 或z=x；y z= 1, 4 2, 5 3, 6 （23） z=x；y 1,2,3,4,5,6 （）2.5 点运算点的运算也就是单个元素之间的运算。为了与矩阵作为整体的运算符号相区别，要在运算符“*”、“/”、“”、“”前加一个“.”，以表示进行点运算。矩阵的加减法本来就是对元素进行的，故不再有点的运算符。参与点运算的两个矩阵必须是同阶的（标量除外）。考虑下面的（34）阶矩阵：a= 1 2 3 4 m=13 14 15 16 5 6 7 8 17 18 19 20 9 10 11 12 21 22 23 24直接写出MATLAB点运算：Zm=x.*m= 1*13 2*14 3*15 4*16 5*17 6*18 7*19 8*20 9*21 10*22 11*23 12*24Zd=x./m= 1/13 2/14 3/15 4/16 5/17 6/18 7/19 8/20 9/21 10/22 11/23 12/24对于点乘来说，如果x=xO（xO为标量），则点乘和乘法运算是一样的。xO.*m = xO*13 xO*14 xO*15 xO*16 xO 13 14 15 16 xO*17 xO*18 xO*19 xO*20 = 17 18 19 20 xO*21 xO*22 xO*23 xO*24 21 22 23 24=xO*m类似地，当m=mO=标量常数时，则 x.*m=x*mO对于点除运算来说，当m=mO且mO为标量常数时：x./mO = 1/mO 2/mO 3/mO 4/mO 1/mO 1 2 3 4 5/mO 6/mO 7/mO 8/mO = 5 6 7 8 9/mO 10/mO 11/mO 12/mO 9 10 11 12=x/mO此时点的运算是不必要的。当x=xO为标量常数时：xO./m =xO/13 xO/14 xO/15 xO/16 xO 1/13 2/14 3/15 4/16 xO/17 xO/18 xO/19 xO/20 = 5/17 6/18 7/19 8/20 xO/21 xO/22 xO/23 xO/24 9/21 10/22 11/23 12/24=xO./m这里的点运算是必需的。很显然，当m=mO为标量常数或x=xO为标量常数时，指数运算总要用到点运算。即： z=xO.m z=x.mO例17：为说明指数的点运算，考虑计算2j，其中j=1,2,8。x=1:8；y=2.x 2 4 8 16 32 64 128 256上面的程序可简写为：y=2.1:8 或y=2.(1:8)反之，如求取x2的值：y=(1:8).2 或y=1:8.2 1 4 9 16 25 36 49 64如果令f(x)代替任意一个函数，如sin、cos、cosh、贝塞尔函数等，则当矩阵y为一个（32）阶矩阵时：z= f(x)=f(y11) f(y12) f(y21) f(y22) f(y31) f(y32)在矩阵阶数相同的情况下，可实现与点运算的联合运算。例如，如果a、b、c、d和f都是（32）阶矩阵，则表达式： 可写为： z=(tan(a)-f.*(b./c).d).2z中各元素的数值实际上是通过以下的表达式计算得出：点的另一个应用是函数meshgrid。利用下述语句： u,v=meshgrid(s,t)通过向量s=s1 s2 s3 s4 和 t=t1 t2 t3 可生成两个（34）阶矩阵：若让u和v中的相关元素相乘，则点乘： z=u.*v z中的元素是向量s和t中元素的所有可能乘积值。当进行加、减、除和指数运算时，结果相似。只不过是将乘法符号（*）用其他运算符所替代而已。检验两个函数：(这两个函数经常用于点操作之中)1、sum当sum(a)函数的参数为向量时，sum函数将计算向量中的所有元素之和，并返回此值。当参数为矩阵时，函数将矩阵中的元素按列相加，并返回一个长度等于原始矩阵列数的行向量，即：2、cumsum设向量v有n个元素vj，cumsum函数对其操作生成另一个长度为n的向量，其元素为：如果w是由元素wjk组成的一个(mn)阶矩阵，则cumsum(w)为如下矩阵：例18：计算 m=1:4； z=sum(m.m) 或z=sum(1:4).(1:4)（为什么前后数值结果不一样？）(因为后边的先求和后求其1 ，2，3，4次方)例19: 通过计算下式对本节的内容做一小结，式中N=305，且x在0，2的区间内按5等分取值。并将计算值与其理论值比较。（没看懂）nn=1:2:305；xx=linspace(0,2,5)； 1.0021 1.0000x,n=meshgrid(xx,nn)； 0.8889 0.9968s=4*pi*sum(n.*(-1).(n-1)/2)./(n*pi).2+4*x.2)；0.6501 0.6481se=sech(xx)； 0.4272 0.4251compare=sw 0.2679 0.2658或x,n=meshgrid(linspace(0,2,5),1:2:3050；s=4*pi*sum(n.*(-1).(n-1)/2)./(n*pi).2+4*x.2)；compare=ssech(linspace(0,2,5)2.6 矩阵的数学运算在本节将介绍矩阵的几种运算：加、减、乘、转置、行列式、方程式及其根（特征值）。这些运算经常用来求解各种工程问题。2.6.1 加法和减法如果有两个阶矩阵a和b，则MATLAB中矩阵加、减法的表达式为: c=ab cij=aijbij2.6.2 乘法 如有阶矩阵a和阶矩阵b,则MATLAB中矩阵乘积的表达式为: c=a*b 矩阵乘积的扩展运算为：f=c*d=a*b*d 。通常情况下abba 。当m=n时，且c=ab，则：c=(ab)=ba。如果a=I，且m=n ,则：Ib=bI=b 。用矩阵乘法来计算级数考虑下列级数：x,y具有一系列的值，设x=x1，x2，xm ；而y=y1，y2，yn 。令f为阶矩阵：且g为阶矩阵：对于x和y的任意一种数值组合，矩阵的乘积运算都可实现级数的求和。注意，在级数的求和运算中，除了要用到矩阵乘法外，还要用到矩阵的点运算和meshgrid命令。可以看出，对于x和y的任何元素组合，矩阵乘积运算为级数求和提供了一种非常简便的方法。考虑矩阵乘积的三种特殊情况：、行向量和列向量乘积。、列向量和行向量乘积。、行向量和矩阵的乘积。在这三种情况的讨论中，提供了一种解决工程问题的方法，这种方法充分利用MATLAB的优点，如简洁的表达式和矩阵的多种运算方法。第种情况：行向量和列向量的乘积。令a为行向量，b为列向量：a=a1 a2 ak b=b1 b2 bk则a与b相乘的结果为标量：d=ab=ajbj j=1,2,k在MATLAB中上述两向量矩阵相乘的表达式定义为：d=a*b 或d=dot(a,b)如果c是一个阶矩阵，x是一个阶列向量，则：即结果也为标量值。上式的MATLAB表达式为：f=这种组合运算在求解特征值问题时极为有用。第种情况：列向量和行向量的乘积。令b为阶列向量，a为阶行向量，则乘积为：f=ba 上式的MATLAB表达式为： f=b*a例20：探讨极坐标到笛卡尔坐标的变换。即 设半径值向量和角度值向量，则相应得笛卡尔坐标值为：这样就把极坐标映射到笛卡尔坐标中去。例21：环形薄膜的形状考虑固态环形薄膜的形状模型，沿外边界受压，外边界r=1:其中J1(x)为第一类1阶贝塞尔函数，为薄膜上任意一点的极坐标值。贝塞尔函数由下式确定。 besselj(n,x)其中n为阶数，x为参数。薄膜中心处为坐标原点，薄膜的固有频率参数为3.8316。该模型外形可由MATLAB命令mesh(x,y,z)绘制，其中(x,y)为表面函数z(x,y)上一个点的平面坐标。以和为增量绘制表面图，程序如下：r=0:0.05:1； % (21x1)phi=0:pi/20:2*pi； % (1x41)x=r*cos(phi)； % (21x41)y=r*sin(phi)； % (21x41)z=besselj(1,3.8316*r)*cos(phi)；% (21x41)mesh(x,y,z)例:22 拉普拉斯方程求解在边界条件下, 拉普拉斯方程的解为:其中表面的程序如下：n=(1:25)*pi; % (1x25)eta=0:0.025:1; % (1x41)xi=0:0.05:7; % (1x15)tempc=meshgrid(1-cos(n)./n.3,xi); % (15x25)tempe=exp(-n*xi); % (15x25)tempec=tempc.*tempe; % (15x25)temps=sin(n*eta); % (25x41)z=4*tempec*temps; % (15x41)mesh(eta,xi,z)第种3情况：行向量和矩阵的乘积。令b为阶矩阵且a为阶行向量。则乘积为： g=ab MATLAB表达式为： g=a*b例23：傅里叶级数求和脉冲宽度为d，周期为T的矩形脉冲的傅里叶级数由下式给出：其中。对的150个值求和并绘制当d/T=0.25，-1/21/2时的图形。图形的绘制可通过函数plot(x，y)实现，其中x=，y=。程序为：n=1:150; % (1x150)tau=linspace(-0.5,0.5,100); % (1x100)sn=sin(pi*n/4)./(pi*n/4); % (1x150)cntau=cos(2*pi*n*tau); % (150x100)f=0.25*(1+2*sn*cntau); % (1x100)plot(tau,f)例24:正态累积分布函数一个近似正态（高斯）累积概率分布函数，在0Xx区间内，概率由下式给出：其中，并且：b1=0.319381520，b2=-0.356563782，b3=1.781477937，b4=-1.821255978，b5=1.330274429由1-P()获得区域，其中。计算并绘制在区域内增量的累积分布的程序为：b=0.319381530 -0.356563782 1.781477937 -1.821255978 1.330274429 m=1:length(b); %(1x5) x=0:0.2:3; %(1x16) mm=meshgrid(m,x); %(16x5) bb=meshgrid(b,x); %(16x5) temp=(1./(1+0.231641*x); %(16x1) zxx=meshgrid(temp,m); %(16x1) bzm=sum(bb.*(zxx.mm); %sum(16x5)-sum(5x16)-(1x16) px=1-bzm.*exp(-0.5*x.2)/sqrt(pi*2); %(1x16) plot(x,px,k,-fliplr(x),fliplr(1-px),k)例25:级数的双重求和考虑下列级数。该级数描述了一个正方形薄膜相对于其原始位置进行规则位移的无量纲相应。其中，b为薄膜的边长，为薄膜的波动速度，t为时间。在时计算级数的程序为：tau=0;n=0:10;m=n;xi=0:0.1:1;eta=xi;nden=meshgrid(2*n+1).3,xi);mden=meshgrid(2*m+1).3,eta);snx=sin(pi*(2*n+1)*xi);sme=sin(pi*(2*m+1)*eta);xterm=snx./nden;yterm=sme./mden;n2,m2=meshgrid(2*n+1,2*m+1,1);w=xterm*cos(pi*tau*sqrt(n2.2+m2.2)*yterm;mesh(xi,eta,w)2.6.3 矩阵的逆方阵a的逆运算: a-1a=aa-1=I上式表明a为非奇异阵，也就是说，它的行列式不等于0（0）。获得矩阵a逆的MATLAB表达式为： inv(a) 或 a-12.6.4 行列式行列式的MATLAB表达式为: det(a)一些工程应用中经常有如下这类问题: 其中a和b是阶矩阵, 为方程的根(或特征值)。线性方程组的解可由下式获得： lambda=eig(a,b)例26:多项式变换有如下形式的多项式:ax2+by2+cz2+2dxy+2exz+2gyz其中a ,b ,c ,d ,e ,g为实数。多项式可变换为实对角阵形式：其中为另一座标系，座标系原点也在（0,0,0）,为行列式：的根，A为实对称阵：考虑多项式：4x2+3y2-z2-12xy+4xz-8yz则：为求取根rj ,可输入： r=eig(4,-6,2；-6,3,-4；2,-4,-1 r=-4.0000,-1.0000,11.0000或 r=-sort(-eig(4,-6,2；-6,3,-4；2,-4,-1) 按降序排列。也可用flipud和fliplr。2.6.5 方程组求解 方程组的矩阵形式可写成: ax=b其中a为阶矩阵，x和b为阶向量。方程的解渴在方程两边同乘以a-1求出，即： a-1ax=a-1b x=a-1b因为a-1a=I且Ix=x，故推荐采用如下表达式求解此方程组。 x=ab其中反斜线运算符代表矩阵的除法，在MATLAB中称之为矩阵的左除。与以下公式相比，左除具有更好的数值稳定性： x=a-1*b 或 x=inv(a)*b例如： MATLAB程序为：a=8,1,6；3,5,7；4,9,2；b=7.5；4；12x=1.2931 0.8972 -0.6236x=ab程序程序可简化为: x=8,1,6；3,5,7；4,9,2 7.5；4；12例2.7 相关方程组的解以正方形平板,四边固定,表面载荷均匀。则求解其静态变形时，首先应先从从下述无限方程组中通过舍去一些位获得常数Em值：其中， 且 。如果仅取方程组的前4个，则：可用如下程序求这4个方程的解：m=1:2:7;i=m;alp=m*pi/2;ai=(tanh(alp)+alp./(cosh(alp).2)./i; %(1x4)ci=4.*(alp./cosh(alp).2-tanh(alp)./(pi3)*i.4); %(1x4)ii,mm=meshgrid(i,m); %(1x4)bim=(8/pi)*ii./(1+(ii.2)./(mm.2).2).*mm.3); %(1x4)format long eee=(diag(ai)+bim)ciformat short练习 2.1 建立两个向量a和b向量，向量a元素为2n-1，向量元素为2n+1，n =0，1，7。(a)求a与b的和?(b)a和b的区