Matlab计算与仿真技术

现在是1页\一共有81页\编辑于星期一Matlab计算与仿真技术现在是2页\一共有81页\编辑于星期一物理与光电工程学院大连理工大学Matlab计算与仿真技术现在是3页\一共有81页\编辑于星期一什么是MATLAB?MATLAB代表MATrixLABoratory它的首创者是美国新墨西哥大学计算机系的系主任CleveMoler博士，他在教授线性代数课程发现其他语言很不方便，篇构思开发了MATLAB。最初采用FORTRAN语言编写，20世纪80年代后出现了MATLAB的第二版，全部采用C语言编写.1984年Moler博士和一批数学家及软件专家创建了MathWorks公司，专门开发MATLAB。现在是4页\一共有81页\编辑于星期一MATLAB的发展历史1984年，MATLAB第1版(DOS版)

1992年，MATLAB4.0版

1994年，MATLAB4.2版

1997年，MATLAB5.0版

1999年，MATLAB5.3版

2000年，MATLAB6.0版

2001年，MATLAB6.1版

2002年，MATLAB6.5版

2004年，MATLAB7.0版

2008年，MATLAB7.6版（支持4线程计算）

2009年，MATLAB7.8版（支持8线程计算）是一种解释性语言，一种演草纸式的科学计算语言现在是5页\一共有81页\编辑于星期一MATLAB能干什么？MATLAB可以进行：数学计算、算法开发、数据采集建模、仿真、原型数据分析、开发和可视化科学和工程图形应用程序的开发，包括图形用户界面的创建。MATLAB广泛应用于：数值计算、图形处理、符号运算、数学建模、系统辨识、小波分析、实时控制、动态仿真等领域。现在是6页\一共有81页\编辑于星期一掌握MATLAB……MATLAB的构成：MATLAB开发环境：进行应用研究开发的交互式平台MATLAB数学与运算函数库：用于科学计算的函数MATLAB语言：进行应用开发的编程工具图形化开发：二维、三维图形开发的工具应用程序接口(API)：用于与其他预言混编面向专门领域的工具箱：小波工具箱、神经网络工具箱、信号处理工具箱、图像处理工具箱、模糊逻辑工具箱、优化工具箱、鲁棒控制工具箱等几十个不同应用的工具箱。现在是7页\一共有81页\编辑于星期一MATLAB界面：命令窗口、图形窗口、编辑窗口、帮助窗口。现在是8页\一共有81页\编辑于星期一常用MATLAB系统命令help在线帮助who显示当前变量whos显示当前变量的详细信息clear清空工作间的变量和函数pack整理工作间的内存load把文件调入变量到工作间save把变量存入文件中echo命令回显what显示指定的matlab文件dir显示目录内容现在是9页\一共有81页\编辑于星期一基本数学运算在MATLAB下进行基本数学运算，只需将运算式直接打在提示号>>后面，并按Enter键。MATLAB将计算的结果以ans显示。【例】求的算术运算结果。（1）用键盘在MATLAB指令窗中输入以下内容>>(12+2*(7-4))/3^2

（2）在上述表达式输入完成后，按【Enter】键，该就指令被执行。（3）在指令执行后，MATLAB指令窗中将显示以下结果。ans=

2

我们也可给运算式的结果设定一个变量x：x=(5*2+1.3-0.8)*10^2/25x=42变量x的值可以在下个语句中调用：y=2*x+1y=85

现在是10页\一共有81页\编辑于星期一例：计算在命令窗口中输入sqrt（18）+（7*sin（pi/6））/（5+cos（pi/6））注意：（1）“》”为命令提示符；（2）字母要区分大小写，所有命令和函数名都必须用小写。现在是11页\一共有81页\编辑于星期一变量命名规则1.变量名的大小写敏感。2.变量的第一个字符必须为英文字母，而且不能超过31个字符。3.变量名可以包含下连字符、数字，但不能为空格符、标点。注：数值型变量和符号型变量在MATLAB中是不相同的，它们之间不能直接进行转化。MATLAB提供了一个将数值型转化成符号型的命令，即sym（syms）。现在是12页\一共有81页\编辑于星期一MATLAB提供基本的算术运算有：加(+)、减(-)、乘(*)、除(/)、幂次方(^)，范例为：5+3,5-3,5*3,5/3,5^3MATLAB书写表达式的规则与“手写算式”差不多相同，但要求所有表达式都是以纯文本形式输入。如果一个指令过长可以在结尾加上...（代表此行指令与下一行连续），例如：>>

1*2+3*4+5*6+7*8+9*10+11*12+...13*14+15*16ans=744现在是13页\一共有81页\编辑于星期一MATLAB常用数学函数三角函数和双曲函数名称含义名称含义名称含义sin正弦csc余割atanh反双曲正切cos余弦asec反正割acoth反双曲余切tan正切acsc反余割sech双曲正割cot余切sinh双曲正弦csch双曲余割asin反正弦cosh双曲余弦asech反双曲正割acos反余弦tanh双曲正切acsch反双曲余割atan反正切coth双曲余切atan2四象限反正切acot反余切asinh反双曲正弦sec正割acosh反双曲余弦现在是14页\一共有81页\编辑于星期一MATLAB常用数学函数指数函数名称含义名称含义名称含义expE为底的指数log1010为底的对数pow22的幂log自然对数log22为底的对数sqrt平方根复数函数名称含义名称含义名称含义abs绝对值conj复数共轭real复数实部angle相角imag复数虚部现在是15页\一共有81页\编辑于星期一其他函数名称含义名称含义min最小值max最大值mean平均值median中位数std标准差diff相邻元素的差sort排序length个数norm欧氏（Euclidean）长度sum总和prod总乘积dot内积cumsum累计元素总和cumprod累计元素总乘积cross外积现在是16页\一共有81页\编辑于星期一矩阵的MATLAB表示直接输入

>>A=[12,3;456;7,89]冒号操作符

>>a=0:1:10>>a=linspace(0,1,10)>>a=logspace(1,2,10)

现在是17页\一共有81页\编辑于星期一MATLAB下矩阵的运算矩阵的代数运算转置B=A’加减乘A+BA-BA*B左除A\B即AX=B的解X=A-1B右除A/B即XB=A的解X=AB-1翻转fliplrflipudrot90乘方A^B点运算A.*BA./BA.\BA.^BA.’函数现在是18页\一共有81页\编辑于星期一矩阵的逻辑运算运算符&(与)|(或)~(非)结果是一个0-1矩阵。当逻辑表达式的值为真时，赋值1，否则为0

逻辑函数all：当某列的元素都为真时，返回值为1，否则返回0。最终运算结果为一个0-1行向量any：当向量中至少有一个元素为真时，返回值为1，否则返回0。最终运算结果为一个0-1行向量find：用于查找向量中的真元素的下标，返回由所有真元素下标构成的列向量。现在是19页\一共有81页\编辑于星期一矩阵的比较关系在MATLAB里共有六个关系运算符

<小于 >大于

<=小于等于>=大于等于

==等于 ~=不等于关系运算符将生成一个0-1矩阵，当运算数相应元素为真时，对应位置上生成1，否则为0。

现在是20页\一共有81页\编辑于星期一矩阵元素的数据变换floor(A)向－∞取整ceil(A)向＋∞取整round(A)四舍五入取整fix(A)按离0近的方向取整现在是21页\一共有81页\编辑于星期一下标操作矩阵下标（按列优先排列）用两个下标来表示。第i行j列的元素用A(i,j)来表示用一个下标来表示。对于矩阵按列操作，可用单下标引用A(i)。在下标的表达式中使用冒号表示矩阵的一部分。下标引用从序号1开始。现在是22页\一共有81页\编辑于星期一下标操作A=123456789>>A(4)+A(8)现在是23页\一共有81页\编辑于星期一行列删除>>B=[123;456;789];>>B(2,:)=[]>>B=123789现在是24页\一共有81页\编辑于星期一矩阵的尺寸为了判断矩阵的尺寸，Matlab提供了以下函数：[m,n]=size(A) %返回矩阵的行数m和列数nM=size(A) %在矩阵M中返回矩阵的行数、列数

len=length(A) %返回矩阵A行数和列数中的最大值

现在是25页\一共有81页\编辑于星期一特殊矩阵单位矩阵

eye(n)

%生成n维单位阵

eye(m,n)

%生成m×n的单位阵

eye([m,n])%生成m×n的单位阵

eye(size(A))%生成与A矩阵同样大小的单位矩阵

oneszeros

现在是26页\一共有81页\编辑于星期一特殊矩阵全1矩阵

ones(n)ones(m,n)ones([m,n]) ones(size(A))全零矩阵

zeros(n)zeros(m,n)zeros([m,n]) zeros(size(A))随机矩阵

rand(n)rand(m,n)rand([m,n])rand(size(A))生成的随机矩阵的各个元素值在0和1之间。对角阵diag(n)对角阵魔术矩阵

magic(n)魔术矩阵的所有行元素和、列元素和、对角线元素和均相等现在是27页\一共有81页\编辑于星期一阵列与矩阵MATLAB的运算事实上是以阵列(array)及矩阵(matrix)方式在做运算.阵列强调元素对元素的运算，而矩阵则采用线性代数的运算方式.

宣告一变量为阵列或是矩阵时，须用中括号[]

将元素置于其中。阵列为一维元素所构成，而矩阵为多维元素所组成.

例如:»x=[12345678];%

一维1x8阵列

【例】简单矩阵的输入步骤。（1）在键盘上输入下列内容:(以;区隔各列的元素)

A=[1,2,3;4,5,6;7,8,9]

（2）按【Enter】键，指令被执行。（3）在指令执行后，MATLAB指令窗中将显示以下结果：

A=123456789现在是28页\一共有81页\编辑于星期一阵列的运算符号

除了加减符号外其余的阵列运算符号均须多加.符号。阵列运算功能(注意:一定要多加.符号)+加-减.*乘./左除.\右除.^次方.’

转置矩阵的几种基本变换操作（1）

通过在矩阵变量后加’的方法来表示转置运算

>>a=[10,2,12;34,2,4;98,34,6];>>a'ans=10349822341246现在是29页\一共有81页\编辑于星期一（2）

矩阵求逆>>inv(a)ans=-0.01160.0372-0.00150.0176-0.10470.03450.0901-0.0135-0.0045（3）

矩阵的特征值>>[u,v]=eig(a)u=-0.29600.3635-0.3600-0.2925-0.41280.7886-0.9093-0.83520.4985v=48.8395000-19.8451000-10.9943现在是30页\一共有81页\编辑于星期一（4）

取出上三角和下三角

>>triu(a)ans=10212024006>>tril(a)ans=1000342098346（5）LU分解

>>[l,u]=lu(a)l=0.10200.15001.00000.34691.000001.000000u=98.000034.00006.00000-9.79591.91840011.1000现在是31页\一共有81页\编辑于星期一MatLab图形绘制功能基本绘图函数命令含义plot建立向量或矩阵各队队向量的图形loglogx、y轴都取对数标度建立图形semilogxx轴用于对数标度，y轴线性标度绘制图形semilogyy轴用于对数标度，x轴线性标度绘制图形title给图形加标题xlabel给x轴加标记ylabel给y轴加标记text在图形指定的位置上加文本字符串gtext在鼠标的位置上加文本字符串grid打开网格线现在是32页\一共有81页\编辑于星期一x=0:0.001:10;%0到10的1000个点的x座标

y=sin(x);%对应的y座标

plot(x,y);%绘图（可以设定显示范围）现在是33页\一共有81页\编辑于星期一画椭圆

a=0:pi/50:2*pi; X=cos(a)*3;

%参数方程Y=sin(a)*2;plot(X,Y);xlabel('x'),ylabel('y');title('椭圆')现在是34页\一共有81页\编辑于星期一线条和颜色控制符:.点y黄。小圆圈m棕色xx标记c青色+加号r红色-实线g绿色*星号b蓝色：虚线w白色-.点划线h黑色现在是35页\一共有81页\编辑于星期一条形图当资料点数量不多时，长条图是很适合的表示方式：closeall;%关闭所有的图形视窗x=1:10;y=rand(size(x));bar(x,y)现在是36页\一共有81页\编辑于星期一误差图x=linspace(0,2*pi,30);y=sin(x);e=std(y)*ones(size(x));%std:标准差errorbar(x,y,e)现在是37页\一共有81页\编辑于星期一极坐标图形theta=linspace(0,2*pi);r=cos(4*theta);polar(theta,r);现在是38页\一共有81页\编辑于星期一x=linspace(0,10,50);y=sin(x).*exp(-x/3);stairs(x,y);%阶梯图x=linspace(0,10,50);y=sin(x).*exp(-x/3);stem(x,y);%针状图现在是39页\一共有81页\编辑于星期一符号表达式画图在许多的场合，将表达式可视化是有利的。MATLAB提供了函数ezplot来完成该任务。>>y='16*x^2+64*x+96'y= 16*x^2+64*x+96>>ezplot(y)现在是40页\一共有81页\编辑于星期一三维曲线和曲面

1．空间曲线绘制函数plot3()格式一：plot3(x,y,z)功能：若x、y、z为相同长度的向量，则根据向量x、y、z绘制空间三维曲线。如果x、y、z为同阶矩阵，则绘制对应列的多条曲线。格式二：plot3(x,y,z,’s’)功能：按字符串s设置的线型、颜色、标记符号绘制三维空间曲线。字符串s的设置与plot命令的设置相同。

2．三维网线图函数mesh、meshc和meshzMATAL提供的函数mesh用来绘制三维网线图。格式一：mesh(X,Y,Z)

功能：根据矩阵X、Y和Z绘制彩色的空间三维网线图。X、Y和Z中对应的元素为三维 空间上的点，点与点之间用线连接。其中网线的颜色随着网点高度的改变而改 变。格式二：mesh(x,y,Z)

功能：n维向量x、m维向量y和m×n矩阵Z绘制网线图，节点的坐标为 (x(j),y(i),Z(i,j))，网线的颜色随着网点高度的改变而改变。现在是41页\一共有81页\编辑于星期一三维曲面图函数surf及surfc

格式一：surf(X,Y,Z,C)功能：绘制由四个矩阵所指定的带色参数的网状表面图。视角是由view所指定。轴的刻度决定于X，Y及Z的范围，或当前对轴的设定；颜色范围由C指定。格式二：surf(X，Y，Z)功能：将C设为与Z相等，则颜色与网的高度成正比。现在是42页\一共有81页\编辑于星期一三维螺旋线t=0:pi/50:10*pi;plot3(sin(t),cos(t),t)grid现在是43页\一共有81页\编辑于星期一空间中的曲线t=linspace(0,20*pi,501);plot3(t.*sin(t),t.*cos(t),t);%注意用点乘.*t=linspace(0,10*pi,501);plot3(t.*sin(t),t.*cos(t),t,t.*sin(t),t.*cos(t),-t);现在是44页\一共有81页\编辑于星期一立体网状图x=linspace(-2,2,25);%在x轴上取25点y=linspace(-2,2,25);%在y轴上取25点[xx,yy]=meshgrid(x,y);%xx和yy都是21x21的矩阵zz=xx.*exp(-xx.^2-yy.^2);%计算函数值，zz也是21x21的矩阵mesh(xx,yy,zz);%画出立体网状图

现在是45页\一共有81页\编辑于星期一网状图填补x=linspace(-2,2,25);%在x轴上取25点y=linspace(-2,2,25);%在y轴上取25点[xx,yy]=meshgrid(x,y);%xx和yy都是21x21的矩阵zz=xx.*exp(-xx.^2-yy.^2);%计算函数值，zz也是21x21的 矩阵surf(xx,yy,zz);%画出立体曲面图现在是46页\一共有81页\编辑于星期一peaks函数peaks函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点，其方程式为：[x,y,z]=peaks;waterfall(x,y,z);现在是47页\一共有81页\编辑于星期一等高线图meshc同时画出网状图与等高线：[x,y,z]=peaks;meshc(x,y,z);surfc同时画出曲面图与等高线：[x,y,z]=peaks;surfc(x,y,z);现在是48页\一共有81页\编辑于星期一剔透玲珑球[X0,Y0,Z0]=sphere(30);%产生单位球面的三维坐标X=2*X0;Y=2*Y0;Z=2*Z0;%产生半径为2的球面的三维坐标surf(X0,Y0,Z0);%画单位球面shadinginterp%采用插补明暗处理holdon;mesh(X,Y,Z);holdoff%画外球面hiddenoff%产生透视效果axisoff%不显示坐标轴现在是49页\一共有81页\编辑于星期一彗星状轨迹图

【*例】简单二维示例。（请指令窗中运行以下指令）n=10;t=n*pi*(0:0.0005:1);x=sin(t);y=cos(t);plot(x,y,'g');holdon;comet(x,y,0.02);holdoff现在是50页\一共有81页\编辑于星期一卫星返回地球的运动轨线示意R0=1; %以地球半径为一个单位a=12*R0;b=9*R0;T0=2*pi; %T0是轨道周期T=5*T0;dt=pi/100;t=[0:dt:T]';f=sqrt(a^2-b^2); %地球与另一焦点的距离th=12.5*pi/180; %卫星轨道与x-y平面的倾角E=exp(-t/20); %轨道收缩率

x=E.*(a*cos(t)-f);y=E.*(b*cos(th)*sin(t));z=E.*(b*sin(th)*sin(t));plot3(x,y,z,'g') %画全程轨线[X,Y,Z]=sphere(30);X=R0*X;Y=R0*Y;Z=R0*Z;%获得单位球坐标gridon,holdon,surf(X,Y,Z),shadinginterp %画地球x1=-18*R0;x2=6*R0;y1=-12*R0;y2=12*R0;z1=-6*R0;z2=6*R0;axis([x1x2y1y2z1z2]) %确定坐标范围view([11737]),comet3(x,y,z,0.02),holdoff %设视角、画运动轨线

现在是51页\一共有81页\编辑于星期一求解线性方程组MATLAB语句实现：

a=[2,3,-1；8,2,3；45,3,9]；

b=[2；4；23]；

x=inv(a)*b现在是52页\一共有81页\编辑于星期一求多项式方程的根roots(p)solve(p)例.求方程的所有根

p=[1-49-10]r=roots(p)

或s1=sym('x^3-4*x^2+9*x-10');solve(s1)现在是53页\一共有81页\编辑于星期一求方程组的根和超越方程解[x,y]=solve('x+y=1','x-11*y=5','x','y')求超越方程的根，其中p，r为常数ff=sym('p*sin(x)=r')solve(ff,'x')现在是54页\一共有81页\编辑于星期一求极限symsx;y1=(1+4*x)^(1/x);y2=(exp(x)-1)/x;limit(y1,x,0)limit(y2,x,0)现在是55页\一共有81页\编辑于星期一求导数与偏导数

symsx;f=3*x^3+5*x+1;diff(f,1)

symsxy;%对x偏导数

z=x^2*sin(2*y);B=diff(z,x)现在是56页\一共有81页\编辑于星期一求解不定积分symsxb;f=cos(b*x);f1=int(f,x)f1=1/b*sin(b*x)ezplot(f1)现在是57页\一共有81页\编辑于星期一求定积分多重积分

symsx;y=(x^2+sin(x))/(1+x^2);int(y,x,-1,1)

symsxy;int(int(x^2+y^2+1,y,x,x+1),x,0,1)现在是58页\一共有81页\编辑于星期一求微分方程通解和特解dsolve('Dy=1+y^2','x')dsolve('Dy=1+y^2','y(0)=1','x')dsolve('x^2*D2y+x*Dy+(x^2-1/2)*y=0','y(pi/2)=2','Dy(pi/2)=-2/pi','x')现在是59页\一共有81页\编辑于星期一MATLAB语言编程介绍MATLAB中各种命令可以完成许多单一的任务，对于某些较为复杂的问题，仅靠现有的命令或函数来解决，往往是难以达到目的。为此，要运用MATLAB编程语言编制程序，形成M-文件。程序是使计算机完成各项运算的命令集，运行一个编制好的程序，计算机会从第一条命令行开始，一行接一行地执行相应的命令，直到终止。程序一般包含：数据输入，数据处理和结果输出三部分，其中数据处理是核心。程序编写调试完成后，需要存盘，形成永久性文件，可以随时对它进行调用或修改。文件名以字母开头，但不能用专用变量名，如pi等。MATLAB中每一个命令都是一个M-文件。现在是60页\一共有81页\编辑于星期一MATLAB程序类型脚本文件函数文件functiony=mean(x)%MEANAverageormeanValue%Forvectors,Mean(x)isthemeanvalueofX%Formatrices,Mean(x)isarowvectorcontaining%themeanvalueofeachcolumn[m,n]=size(x)ifm==1m=n;endy=sum(x)/m;函数定义行H1行函数帮助文件函数体现在是61页\一共有81页\编辑于星期一M-文件中输入、输出命令

直接赋值输入：提示对话输入(input命令)x=2;x=input('请输入参数x=');a=[1,2;3,4];a=input('请输入矩阵a=');s='anystring!';s=input('Pleaseinputs=');

直接输出：格式控制输出(fprintf命令)xfprintf('x=%.0f,y=%.5f\n',pi,pi);disp([a,b]);fprintf('x=%5g,y=%10.5g\n',N,pi);现在是62页\一共有81页\编辑于星期一M-文件中循环控制命令(for命令)格式：fori=n1:(step):n2 commands;end作用：重复执行命令集commands.现在是63页\一共有81页\编辑于星期一例2：求奇数和：s=1+3+5+…+(2k-1)n=input('请输入n='); s=0;fori=1:2:n s=s+i;fprintf('i=%.0f,s=%.0f\n',i,s);end例1：求和：s=1+2+…+nn=input('请输入n='); s=0;fori=1:n s=s+i;fprintf('i=%.0f,s=%.0f\n',i,s);end现在是64页\一共有81页\编辑于星期一例3：求阶乘：p=1×2×3×…×n=n!n=input('请输入n='); p=1;fori=1:n p=p*i;fprintf('i=%.0f,p=%.0f\n',i,p);end例4：求e：e=1+1+1/2!+1/3!+…+1/n!n=input('请输入n='); p=1;e=1;fori=1:n p=p*i;p1=1/p;e=e+p1;fprintf('i=%.0f,p=%.0f,e=%.8f\n',i,p,e);end现在是65页\一共有81页\编辑于星期一例5：求pi：pi/4=1-1/3+1/5-1/7+1/9-1/11+……n=input('请输入n='); pi=0;k=-1;fori=1:2:n k=(-1)*k;pi=pi+k*1/i;fprintf('i=%.0f,pi=%.8f\n',i,4*pi);end现在是66页\一共有81页\编辑于星期一M-文件中条件循环命令(while命令)格式：while(conditionistrue) commands;end作用：重复执行命令集commands.现在是67页\一共有81页\编辑于星期一例6：求和：s=1+2+…+n<3000 s=0;k=1;whiles<3000 s=s+k;fprintf('k=%.0f,s=%.0f\n',k,s);k=k+1;end例7：求e：e=1+1+1/2!+1/3!+…(1/n!<1.0e-8) p=1;e=1;r=1;i=1;whiler>=1.0e-8 p=p*i;r=1/p;e=e+r;fprintf('i=%.0f,r=%.9f,e=%.9f\n',i,r,e);i=i+1;end现在是68页\一共有81页\编辑于星期一M-文件中选择控制命令(if命令)单项选择控制格式：if(conditionistrue) commands;end作用：若条件成立，则执行命令集

commands.否则，不执行。现在是69页\一共有81页\编辑于星期一例8：求n个实数中最大的数M.a=input('请输入数组a[n]=');[m,n]=size(a);M=a(1); fori=2:n ifM<a(i)M=a(i);endfprintf('M=%.5f,a(%.0f)=%.5f\n',M,i,a(i));end现在是70页\一共有81页\编辑于星期一M-文件中选择控制命令(if命令)多项选择控制格式：if(conditionistrue) commands;elseif(conditionistrue)commands;elsecommands;end作用：若条件成立，则执行命令集

commands.否则，不执行。现在是71页\一共有81页\编辑于星期一例9：建立符号函数sign(x)x=input('x=');ifx>0sn=1;elseifx==0sn=0;elsesn=-1;endfprintf('x=%.5f,sn=%.0f\n',x,sn);现在是72页\一共有81页\编辑于星期一Matlab编程(续)程序设计技术循环向量化计算模式目的:向量计算替代循环模式实例:i初始化

ticticfort=0:0.01:100t=0.0:0.01:100

i=i+1;y(i)=sin(t)y=sin(t)e