您所在位置：网站首页 > 学术论文 > 其它学术论文 > MATLAB及在控制系统课程中的应用PPT教学课件（完整版）

MATLAB及在控制系统课程中的应用PPT教学课件（完整版）.ppt

357页

卖家[上传人]：公****

文档编号：584888067

上传时间：2024-09-01

文档格式：PPT

文档大小：2.94MB

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

20金贝

下载

/ 357 举报版权申诉马上下载

文本预览

下载提示

常见问题

MATLAB及在控制系统课程中的应用2010年12月20日2024/9/11参考书目参考书目1 1、、MatlabMatlab及在电子信息课程中的应用（第二版）及在电子信息课程中的应用（第二版）陈怀琛陈怀琛电子工业出版社电子工业出版社 2004 2004年年1 1月月2 2、、控制系统仿真与计算机辅助设计控制系统仿真与计算机辅助设计薛定宇（东北大学）机械工业出版社 2005年1月3 3、、控制系统数字仿真与控制系统数字仿真与CAD (CAD (第二版第二版) ) 张晓华（哈尔滨工业大学）机械工业出版社 2006年5月4 4、、控制系统的数字仿真与计算机辅助设计控制系统的数字仿真与计算机辅助设计钱积新等化学工业出版社 2003年5月5、基于基于MATLABMATLAB的系统分析与设计的系统分析与设计——控制系统控制系统楼顺天等西安电子科技大学出版社6 6、、MATLAB6.X MATLAB6.X 教程教程7 7、、MATLAB与控制系统仿真实践与控制系统仿真实践定价：34 元作者：赵广元书号：978-7-81124-787-9 北京航空航天大学出版社2024/9/12目录•第1章 MATLAB 语言概述•第2章基本语法•第3章 MATLAB 的开发环境和工具•第4章 MATLAB 的其他函数库•第5章 MATLAB的的SIMULINK仿真仿真•第6章 MATLAB在自动控制原理中应用2024/9/13第1章 MATLAB语言概述•1.1 MATLAB语言的发展•1.2 MATLAB语言的特点•1.3 MATLAB的工作环境 1.3.1 命令窗 1.3.2 图形窗 1.3.3 文本编辑窗•1.4 演示程序•1.5 网络资源2024/9/141.1 MATLAB语言的发展•1.1.1 MATLAB 概述 MATLAB是集数值计算、符号运算及图形处理等强大功能于一体的科学计算语言，是一种交互式的以矩阵为基础的系统计算平台,它用于科学和工程的计算与可视化。

它的优点在于快速开发计算方法，而不在于计算速度 MATLAB已成为一门高校必修的课程，也是最为普遍的计算工具之一2024/9/151.1 MATLAB语言的发展(续)•1.1.2 Matlab的发展§ MATLAB名字由MATMATrix和 LABLABoratory 两词的前三个字母组合而成那是20世纪七十年代，时任美国新墨西哥大学计算机科学系主任的Cleve Moler出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK矩阵软件工具包库程序的的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB§1984年由Little、Moler、Steve Bangert合作成立MathWorks公司，并把MATLAB正式推向市场从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能2024/9/16§1997年仲春，年仲春，MATLAB5.0版问世，紧接着是版问世，紧接着是5.1、、5.2，以及和，以及和1999年春的年春的5.3版现今的版现今的MATLAB拥有更丰拥有更丰富的数据类型和结构、更友善的面向对象、更加快速精富的数据类型和结构、更友善的面向对象、更加快速精良的图形可视、更广博的数学和数据分析资源、更多的良的图形可视、更广博的数学和数据分析资源、更多的应用开发工具。

应用开发工具§2000年末又推出6.0版本无论在界面的设计上还是在内容上较以前版本都有很大的进展2024/9/171.1 MATLAB语言的发展(续)•1.1.3 Matlab的版本演化1.Matlab 1.02.Pc matlab->matlab 3863.Matlab3.5+simulink4.Matlab 4.0：simlink内嵌（1992）5.Matlab 5.0 :全面的面向对象6.Matlab 5.1~5.3 （1999）7.Matlab 6.0 （2000）8.Matlab 6.5：购并了MATRIXx9.Matlab 7.0: （2004）2024/9/181.2 MATLAB语言的特点•友好的工作平台和编程环境•简单易用的程序语言•强大的科学计算及数据处理能力•出色的图形处理功能•应用广泛的模块集和工具箱•实用的程序接口和发布平台•模块化的设计和系统级的仿真2024/9/191.3 MATLAB的工作环境1.3.1 命令窗(Command Window) 单行命令执行方式执行结果直接显示1.3.2 图形窗(Figure Window) 用图形方式表示计算结果1.3.3 文本编辑窗(File Editor) 多行命令组成语言组，可以文件方式存盘•下面就具体看一下MATLAB 的工作环境演示。

2024/9/1101.4 演示程序•在MATLAB的命令窗中键入 demo或demos2024/9/111%pend.mplot([-0.2,0.2],[0;0],'color','y','linestyle','-','linewidth',10);g=0.98;l=1;theta0=pi/6;x0=l*sin(theta0);y0=-l*cos(theta0);axis([-0.75,0.75,-1.25,0]);axis('off');head=line(x0,y0,'color','r','linestyle','.','erasemode','xor','markersize',40);body=line([0;x0],[0,y0],'color','b','linestyle','-','erasemode','xor');t=0;dt=0.01;while t<=50 t=t+dt; theta=theta0*cos(sqrt(g/l)*t); x=l*sin(theta);y=-l*cos(theta); set(head,'xdata',x,'ydata',y); set(body,'xdata',[0;x],'ydata',[0;y]); drawnow;end2024/9/112部分命令的演示例1.求的算术运算结果。

12+2*(7-4))/3^2 ans = 2 例2.简单矩阵的输入步骤A = [1,2,3; 4,5,6; 7,8,9] A = 1 2 3 4 5 6 7 8 92024/9/113例例3. 矩阵的分行输入矩阵的分行输入A=[1,2,34,5,67,8,9] A = 1 2 3 4 5 6 7 8 9 例例4. 指令的续行输入指令的续行输入S=1–1/2+1/3–1/4+1/5–1/6+1/7 …-1/8 S = 0.63452024/9/114例5. 复数表达，及计算 1）z1= 3 + 4i z1 = 3.0000 + 4.0000i （2）z2 = 1 + 2 * iz3=2*exp(i*pi/6)z=z1*z2/z3 z2 = 1.0000 + 2.0000iz3 = 1.7321 + 1.0000iz = 0.3349 + 5.5801i 2024/9/115例6. 复数矩阵的生成及运算•A=[1,3;2,4]-[5,8;6,9]*i•B=[1+5i,2+6i;3+8*i,4+9*i] •C=A*B •A =1.0000 - 5.0000i 3.0000 - 8.0000i• 2.0000 - 6.0000i 4.0000 - 9.0000i•B =• 1.0000 + 5.0000i 2.0000 + 6.0000i• 3.0000 + 8.0000i 4.0000 + 9.0000i•C =• 1.0e+002 *• 0.9900 1.1600 - 0.0900i• 1.1600 + 0.0900i 1.3700 2024/9/116例7 . 求上例复数矩阵C的实部、虚部、模和相角。

•C_real=real(C)•C_imag=imag(C)•C_magnitude=abs(C)•C_phase=angle(C)*180/pi •C_real = 99 116• 116 137•C_imag = 0 -9• 9 0•C_magnitude = 99.0000 116.3486• 116.3486 137.0000•C_phase = 0 -4.4365• 4.4365 0 2024/9/117•例8. 用MATLAB计算能得到 –2 吗？•（1）a=-8; r=a^(1/3) •r = 1.0000 + 1.7321i •（2）全部方根计算如下•m=[0,1,2];•R=abs(a)^(1/3);•Theta=(angle(a)+2*pi*m)/3;•rrr=R*exp(i*Theta) •rrr =• 1.0000 + 1.7321i -2.0000 + 0.0000i 1.0000 - 1.7321i 2024/9/118•（3）图形表示•t=0:pi/20:2*pi;x=R*sin(t);y=R*cos(t);•plot(x,y,'b:'),grid•hold on•plot(rrr(1),'.','MarkerSize',30,'Color','r')•plot(rrr([2,3]),'o','MarkerSize',15,'Color','b')•axis([-3,3,-3,3]),axis square•hold off 2024/9/119例9. 画出衰减振荡曲线及其它的包络线。

t的取值范围是 t=0:pi/50:4*pi;y0=exp(-t/3);y=exp(-t/3).*sin(3*t);plot(t,y,'-r',t,y0,':b',t,-y0,':b')2024/9/120•例10.画出所表示的三维曲面的取值范围是[-8,8].•clear;•x=-8:0.5:8; y=x';•X=ones(size(y))*x;•Y=y*ones(size(x));•R=sqrt(X.^2+Y.^2)+eps;•Z=sin(R)./R; •mesh(X,Y,Z);colormap(hot) xlabel('x'),ylabel('y'),zlabel('z') 2024/9/1211.5 网络资源1.USENET新闻组–MATLAB的新闻组是comp.soft-sys.MATLAB浏览器指向– –2.网络上的工具箱– –1.5 网络资源(续)3.BBS–哈尔滨工业大学bbs: telnet://–上海交通大学bbs: telnet://–清华大学bbs mathtools4.www服务–––Matlab 大观园:– http://www.matlab-2024/9/123第2章基本语法•2.1 变量及其赋值•2.2 矩阵的初等运算•2.3 元素群运算•2.4 逻辑判断及流程控制•2.5 基本绘图方法•2.6 M文件及程序调试2024/9/1242.1 变量及其赋值•2.1.1 标识符与数标识符是标识变量名、常量名、函数名、文件名的字符串的总称。

1、表示符第1个字符必须是字母 2、长度不超过31个 3、区分大小写 4、变量中不能含有标点符号 5、变量可直接参与计算 6、变量一般无需事先定义2024/9/1252.1.1 标识符与数(续)7、特殊变量2024/9/1262.1.1 标识符与数(续)8、数值显示格式•MATLAB中所有的量为双字长浮点数，显示按下面显示规则：1.在缺省情况下，当结果为整数，作为整数显示；当结果为实数，以小数后4位的精度近似显示2. 如果结果中的有效数字超出了这一范围，以科学计数法显示结果3.format命令改变显示格式，常用的的格式有•long (16位) bank(2个十进制位) hex(十六进制)•short(缺省) short e(5位加指数) +(符号)• long e(16位加指数) rat(有理数近似)2024/9/1272.1.2 矩阵及其元素的赋值•矩阵获取格式：变量=表达式（或数）1、直接输入：A=[1 2 3；4 5 6；7，8，9] *矩阵用中括号括起矩阵用中括号括起 *元素间用空格隔开，或用逗号隔开。

元素间用空格隔开，或用逗号隔开 *每行用分号；号表示回车每行用分号；号表示回车2、行向量 B=[1 2 3 4 5]3、列向量 C=[1；2；3；4；5]; 每行命令后面的分号；表示结果不显示每行命令后面的分号；表示结果不显示2024/9/1282.1.2 矩阵及其元素的赋值(续)4、元素可用表达式表示 D=[-1.3 sqrt(3) (1+2+3)/5+1]5、用语句生成行向量 E=from:step:to 即E=开始数：步长：结束数 E=1：2：10 得E=[1 3 5 7 9]6、矩阵连接 B=[a b] V=[a;b]2024/9/1292.1.2 矩阵及其元素的赋值(续)7、用函数创建如: zeros(m,n) ones(m,n) eye(m,n) zeros(3); zeros(3,3); zeros(2,3); zeros(3,2); •ones(3); ones(3,3); ones(2,3); ones(3,2); •eye(3); eye(3,3); eye(3,4); eye(4,3);2024/9/1302.1.2 矩阵及其元素的赋值(续)rand(m,n) %产生均匀分布随机数（产生均匀分布随机数（0，，1））rand(‘state’,0) %把均匀分布伪随机发生器置为把均匀分布伪随机发生器置为0状态状态randn(m,n) %产生正态分布随机数产生正态分布随机数 magic(m) %产生魔方数组（对高维不适用）产生魔方数组（对高维不适用） %即每行、每列及对角元素之和为即每行、每列及对角元素之和为(n^3+n)/2linspace(a,b,n) %在在a和和b之间均匀产生之间均匀产生n个点的值个点的值如：如：f=linspace(0,1,5) 则则 f=0 0.25 0.5 0.75 1.0logspace(a,b,n) %在在a和和b之间对数分布产生之间对数分布产生n个点的值个点的值如：如：f=logspace(0,1,5) 则则 f=1.0000 1.7783 3.1623 5.6234 10.00002024/9/1312.1.2 矩阵及其元素的赋值(续)•矩阵中的元素（用圆括号中数字来注明）1. A( i, j ) 表示第i 行，第j列元素。

2. A( i ) 表示第i个元素矩阵中元素的排序如右所示矩阵中元素的排序如右所示3. A( i, j)=常量，表示给A中元素赋值当下标超出原矩阵的尺寸，则自动扩展行列并补零当下标超出原矩阵的尺寸，则自动扩展行列并补零2024/9/1322.1.2 矩阵及其元素的赋值(续)4. A( : , j ) 表示A阵中第j 列所有元素5. A( i , : ) 表示A阵中第i 行所有元素6. A(2:3,4:6) 表示第2行到第3行，第4列到第6列的子矩阵7. A(3:7) 指A阵中第3个到第7个元素（列优先）（列优先）矩阵的序号编址：按列计数8. A(2)=[ ] 表示去除矩阵中元素此时矩阵变为行矩阵9. A( : ) 指A阵中所有元素组成列向量2024/9/1332.1.3 复数•复数的虚部部分用i 或j表示如：2+3i ,3-4j•复数可直接计算• 如如：：z=[2+3i;3-4j] 或 f=z+[2+j;3];•复数的实部和虚部可分别赋值但 i和j需先清除• 如：clear i j• f=[1,3;5,7]+[2,4;6,8]*j2024/9/1342.1.3 复数•B=Z’ 表示共轭转置。

•B=conj(Z)表示共轭如： Z=[1+2i,3-4j]则： B=Z’ 有 B= 1-2i 3+4j B=conj(Z) 有 B= [1-2i,3+4j]2024/9/1352.1.3 复数•B=conj(Z)’表示转置•B=Z.’ 表示非共轭复数转置如： Z=[1+2i,3-4j]则： B=conj(Z)’ 有 B= 1+2i 3-4j B= Z.’ 有 B= 1+2i 3-4j2024/9/1362.1.4 变量的查询，存储，提取•变量的查询 who 或 whos•变量的存储 save 文件名[.mat] 变量列表如：save sar a b c 变量中间用空格隔开，不能加逗号 •变量的提取 load 文件名•变量的清除 clear 变量列表•清除所有变量 clear all2024/9/1372.1.5 基本赋值矩阵•为了方便给大量元素赋值，MATLAB提供了一些基本矩阵。

见书中表2.1•如：A=zeros(m,n) 全0矩阵B=ones(m,n) 全1矩阵C=eye(m,n) 单位矩阵D= rand(m,n) 0—1之间随机数均匀分布randn(‘state’,0’); %把随机数发生器置0E=randn(m,n) 均值为0的，单位方差正态分布随机矩阵 F= magic(m) 魔方矩阵2024/9/138G= linspace(a,b,n) 线性分隔，a,b之间均匀产生n个数H= logspace(a,b,n) 对数分隔， a,b之间产生n个数K=diag(A); 取A中对角线元素得到列向量列向量P=diag(diag(A)) 产生对角阵 a=[1 2 3 4]; b=diag(a) 产生对角阵如 A=[1 2 3;4 5 6;7 8 9] B=diag(A) 则： B=[1；5；9]；2024/9/1392.2 矩阵的初等运算•2.2.1 矩阵的加减乘除矩阵的加减乘除1、、+，，-，，*，，/，， \ 2、点乘：、点乘：.* 右除：右除：./ 左除：左除：.\•C=A+B; C=A-B C=A*B 注意：矩阵注意：矩阵必须相匹配必须相匹配• X=A\B 表示表示AX=B X=A-1B 即 X=inv(A)*B• X=A/B 表示表示XB=A X=AB-1 即 X=A*inv(B)•[m,n]=size(A) 计算矩阵计算矩阵A的行列大小的行列大小•K=length(A) 计算矩阵计算矩阵A的行列大小中最大的数的行列大小中最大的数2024/9/1402.2.1 矩阵的加减乘除•点乘、点除•C=A.*B 对应元素间相乘。

2024/9/141•C=A./B 对应元素间相除•C=A.\B2024/9/1422.2.2 矩阵除法及线性方程组的解•方阵的行列式 B=det(A) 即B=|A|•方阵的求逆 B=inv(A) 即B= A-1 条件|A|≠0•方阵的伪逆矩阵 B=pinv(A) 条件|A|=0•方阵的伴随矩阵 B=inv(A)*det(A) 即B= A-1 |A|2024/9/1432.2.3 矩阵的乘方和幂次函数1、、^矩阵乘方矩阵乘方2、、.^元素对元素的乘方元素对元素的乘方•C=A^n 表示A阵自乘n次•C=A^(-n) 表示A阵的逆矩阵自乘n次•C=A.^n 表示A阵中每个元素自乘n次•C=A.^（-n）表示A阵中每个元素自乘n次后的逆阵2024/9/144•如 •C=A^2•C=A .^2•C=A ^(-2)====inv(A)^2•C=A .^(-2)2024/9/1452.2.4 矩阵结构形式的提取与变换•B=fliplr(A) %将A矩阵左右翻转 •B=flipud(A) %将A矩阵上下翻转 •B=reshape(A,m,n) %将A阵重组为mxn矩阵•B=rot90(A) %将A矩阵逆时针翻转90度•B=diag(A) %提取A矩阵的对角组成列向量•B=tril(A) %提取A矩阵的左下三角部分•B=triu(A) %提取A矩阵的右上三角部分2024/9/146•如：•B=fliplr(A)•B=flipud(A)•B=rot90(A)•B=tril(A)2024/9/1472.3 元素群运算•2.3.1 数组及其赋值1、t=[初值：步长：终值]；如t=0:0.1:1 tt=[10:-1:1]2、t=linspace(初值,终值,点数）如：tr=linspace(0, 2*pi, 9)3、t=logspace(初值,终值,点数）如：tp=logspace(0, 1, 11)2024/9/1482.3.2 元素群的四则运算•表示对矩阵中每个元素进行运算如如 X=[1 2 3]; Y=[4 5 6]•Z=X.*Y Z=[4 10 18]•Z=X.\Y Z=[4 2.5 2]•Z=X.^Y Z=[1 32 729] •Z=X.^N N=2 Z=[1 4 9]•Z=2.^[X Y] Z=[2 4 8 16 32 64]2024/9/1492.3.3 元素群的函数•等命令可以直接MATLAB中exp、sprt、sin、cos使用在矩阵上，这种运算只是定义在矩阵的单个元素上，即分别对矩阵的每个元素进行运算。

MATLAB中也提供了基本的三角函数 •注意其中的取整注意其中的取整函数名函数名含含义abs绝对值或者复数模sqrt平方根real实部imag虚部conj复数共轭round4舍5入到整数fix舍入到最接近0的整数floor舍入到最接近-∞的整数ceil舍入到最接近∞的整数2024/9/1502.3.3 元素群的函数函数名函数名含含义sign符号函数rem留数sin正弦cos余弦tan正切asin反正弦acos反余弦atan反正切atan2第四象限反正切函数名函数名含含义sinh双曲正弦cosh双曲余弦tanh双曲正切exp自然指数log自然对数log10以10为底的对数bessel贝赛尔函数gamma伽吗函数rat有理逼近2024/9/1512.4 逻辑判断及流程控制2.4.1 关系操作符•MATLAB常用的关系操作符有：<(小于)、<=(小于或等于)、>(大于)、>=(大于或等于)、 = =(等于)、 ~=(不等于)•MATLAB的关系操作符可以用来比较两个大小相同的数组，或者比较一个数组和一个标量在与标量比较时，结果和数组大小一样•»a=1:9;•b=a>4•b = 0 0 0 0 1 1 1 1 1•»c=a(a>4)•c = 5 6 7 8 92024/9/1522.4.1 关系操作符•矩阵查找和排序•子矩阵的查找使用find命令完成，它返回关系表达式为真的下标。

例如：•»a=10:20;•»find(a>15)•ans =• 7 8 9 10 11•矩阵的排序使用sort函数，它将矩阵按照升序排列2024/9/1532.4.2 逻辑运算•逻辑操作符定义了一种与或非的关系表达式MATLAB的逻辑操作符有&(与)、|(或)、~(非)、xor(异或）例如：•»c=~(a>4)•c =1 1 1 1 0 0 0 0 0•»c=(a>4)&(a<7)•c =0 0 0 0 1 1 0 0 0•C=xor(A,B)2024/9/1542.4.3 其他关系与逻辑函数•xor(x,y) 异或运算x或y非零(真)返回1，x和y都是零(假)或都是非零(真)返回0•any(x) 如果在一个向量x中，任何元素是非零，返回1；矩阵x中的每一列有非零元素，返回1•all(x) 如果在一个向量x中，所有元素非零，返回1；矩阵x中的每一列所有元素非零，返回12024/9/155%逻辑函数的运用示例randn('state',1),R=randn(3,6) %创建正态随机阵 L=abs(R)<0.5|abs(R)>1.5 %不等式条件运算，结果给出逻辑数组 R(L)=0%"逻辑1"对应的元素赋0值。

s=(find(R==0))'%利用find获得符合关系等式条件的元素"单下标" R(s)=111%利用"单下标"定位赋值 [ii,jj]=find(R==111);%利用find获得符合关系等式条件的元素"双下标"disp(ii'),disp(jj')2024/9/156【例】关系运算运用之一：求近似极限，修补图形缺口t=-2*pi:pi/10:2*pi;y=sin(t)./t;subplot(1,2,1),plot(t,y),axis([-7,7,-0.5,1.2]),xlabel('t'),ylabel('y'),title('残缺图形')tt=t+(t==0)*eps;yy=sin(tt)./tt;subplot(1,2,2),plot(tt,yy),axis([-7,7,-0.5,1.2])xlabel('t'),ylabel('yy'),title('正确图形') Warning: Divide by zero. 2024/9/1572024/9/158【例】逻辑操作应用之一：逐段解析函数的计算和表现本例演示削顶整流正弦半波的计算和图形绘制•t=linspace(0,3*pi,500);y=sin(t); •z1=((t2*pi)).*y;•w=(t>pi/3&t<2*pi/3)+(t>7*pi/3&t<8*pi/3);•w_n=~w;•z2=w*sin(pi/3)+w_n.*z1;•subplot(1,3,1),plot(t,y,':r'),ylabel('y')•subplot(1,3,2),plot(t,z1,':r'),axis([0 10 -1 1])•subplot(1,3,3),plot(t,z2,'-b'),axis([0 10 -1 1]) 2024/9/1592024/9/160【例】写出生成下图所示波形的MATLAB脚本文件M文件。

图中虚线为正弦波，要求它的负半波被置零，且在处被削顶 2024/9/161t=linspace(0,3*pi,500);y=sin(t); z1=((t2*pi)).*y;w=(t>pi/3&t<2*pi/3)+(t>7*pi/3&t<8*pi/3);w_n=~w;z2=w*sin(pi/3)+w_n.*z1;plot(t,y,'r-.',t,z2,'k-')xlabel('x');ylabel('z=f(t)');title('逐段解析函数');legend('y=sin(t)','z=f(t)');2024/9/1622.4.4 流程控制语句•脚本•对于简单问题，使用直接输入命令简单有效；对稍复杂和多次重复的应用，直接输入命令比较麻烦MATLAB提供了逻辑解决方案，它允许用户把多个命令放在一个简单的文本文件中，如同在MATLAB中键入命令一般，这种文件称为脚本文件，由于脚本文件以m为扩展名，它常称为M文件脚本文件为文本形式的，对跨平台处理十分有利•使用脚本文件，可以把命令保存在磁盘上，便于以后的访问；同时对使用大的数组也带来的方便；增加注释可以为脚本中的命令作文挡以免以后忘记。

2024/9/1632.4.4 流程控制语句•为了便于应用，MATLAB提供了一些流程控制的命令这些命令对脚本编写带来了一些方便，但是需要注意的是，尽量不要使用这些流程控制命令，尤其是循环控制命令•If 语句•很多情况下，命令的序列必须根据关系的检验有条件的执行，它由if-else-end结构提供它的结构如下：–if expression1– commands1–elseif expression2– commands2–elseif …– …–else– commands–end2024/9/1642.4.4 流程控制语句•在执行过程中，MATLAB依次检查各个表达式，只执行第一个表达式为真的命令串，接下来的关系表达式不检验，跳过其余的if-else-end结构，而且，最后的else命令可有可无2024/9/165求x=input('x=');if x>=10 y=x^2+3;elseif x>=0 y=x^3+4*x;else y=x^5+x;end yx=input('x=');if x>=10 y=x^2+3;else if x>=0 y=x^3+4*x; else y=x^5+x; endendy2024/9/1662.4.4 流程控制语句 for循环•for循环允许一组命令以固定的次数重复，它的一般形式是–for x=array– command–end•for 和end之间的命令串按数组array的每一列执行一次，直到n次后终止。

•如：for j=1:2:10 y=j+j.^2; end2024/9/1672.4.4 流程控制语句1.for循环不能使用内部重新赋值循环变量而终止；2.for循环内部接受任何有效的MATLAB数组；3.for循环可以嵌套；4.只要有矩阵形式可以解决的问题，不要使用for循环使用for循环的算法执行很慢，一个好的MATLAB算法不应当出现循环语句Tic/toc5.循环可以使用break跳出，但只跳出所在的循环，不跳出整个嵌套结构2024/9/1682.4.4 流程控制语句•while循环•与for循环以固定的次数求一组指令相反，while循环以不定的次数求一组语句的值While循环的一般形式为：–while expression– commonds–end•只要表达式expression里的所有元素为真，就执行命令串commands通常表达式求值给一个标量值，单数组值也同样有效2024/9/169求y=0;for x=1:100 y=y+x;endyN=input('N=');y=0;for i=1:N for j=1:N y=y+1/(i+j); endendyN=input('N=');y=0;i=1;while i<=N j=1; while j<=N y=y+1/(i+j); j=j+1; end i=i+1;endy2024/9/170【例】Fibonacci数组的元素满足Fibonacci 规则：，；且。

现要求该数组中第一个大于10000的元素•a(1)=1; a(2)=1;•i=2;•while a(i)<=10000• a(i+1)=a(i-1)+a(i); • i=i+1;•end;•i•a(i) •i = 21 •ans = 10946 2024/9/171用for循环指令来寻求Fibonacc数组中第一个大于10000的元素•n=100;a=ones(1,n);•for i=3:n• a(i)=a(i-1)+a(i-2);• if a(i)>=10000• a(i)• break; • end;•End•i •ans = 10946•i = 21 2024/9/1722.4.4 流程控制语句Switch 语句是一种均衡实现的多分支语句–Switch expression–Case 值1– commands1–Case 值2 – commands2 – …–Otherwise– commandsN–end2024/9/173学生的成绩管理，用来演示switch结构的应用。

clear;for i=1:10 a{i}=89+i;b{i}=79+i;c{i}=69+i;d{i}=59+i;end;c=[d,c];Name={' Jack','Marry','Peter',' Rose',' Tom'};Mark={72,83,56,94,100};Rank=cell(1,5);S=struct('Name',Name,'Marks',Mark,'Rank',Rank);2024/9/174for i=1:5 switch S(i).Marks case 100 S(i).Rank='满分满分'; case a S(i).Rank=' 优秀优秀'; case b S(i).Rank=' 良好良好'; case c S(i).Rank=' 及格及格'; otherwise S(i).Rank='不及格不及格'; endenddisp([‘学生姓名学生姓名 ’,‘ 得分得分 ’,‘ 等级等级’]); disp(' ')for i=1:5; disp([S(i).Name,blanks(6),num2str(S(i).Marks),blanks(6),S(i).Rank]);end; 结果：结果：学生姓名得分等级 Jack 72 及格 Marry 83 良好 Peter 56 不及格 Rose 94 优秀 Tom 100 满分 2024/9/1752.5 基本绘图方法2.5.1 直角坐标中的两维曲线•plot(y) 以y的下标作为x坐标，以y值作为y坐标。

•plot(x,y) 数组x和y的长度应匹配•每次绘制将清除以前的图形2024/9/1762.5.1 直角坐标中的两维曲线•图形的标注和图例1、title(‘text’) %给图形加上标题2、xlabel(‘text’) %给X轴加上说明3、ylabel(‘text’) %给Y轴加上说明4、zlabel(‘text’) %给Z轴加上说明5、text(x,y,’string’) %在图形指定位置加上说明6、gtext(‘string’) %利用鼠标在图形加上说明7、legend(‘string1’,’string2’,..) %给图形加图例8、legend off %关闭图例2024/9/177如：作y=sin(t)的二维图形•t=linspace(0,3*pi,200);•y=sin(t);•plot(t,y);•title('y=sin(t)');•xlabel('t/s');•ylabel('y=sin(t)');•text(3,0.4,'y=sin(t)');•legend('y=sin(t)');•gtext('y=sin(t)')2024/9/1782.5.2 线型、点型和颜色•plot(x,y,’r:’) 后面是颜色和线型标识符颜色标识符线型标识符线型yYellow 黄.点S正方形标记mMagenta 品红o圆圈D菱形标记cCyan 青xX号^朝上三角形r Red 红++号V朝下三角形gGreen 绿-实线>朝右三角形bBlue 蓝*星号<朝左三角形wWhite 白:虚号P五角星kBlack 黑-.点划线H六角星- -虚线none无符号标记2024/9/1792.5.3 多条曲线的绘制1、plot(x1,y1,x2,y2);2、plot(x1,y1,’r’,x2,y2);3、plot(x1,y1) hold 是乒乓切换是乒乓切换 hold on %图形保持 plot(x2,y2,’r’) hold off %解除保持4、plot(t,[y1;y2;y3]) %自动给颜色和线型。

5、plotyy(x1,y1,x2,y2) %可画2个不同纵坐标的图2024/9/180•t=0:0.1:3*pi;•y1=sin(t);•y2=cos(t);•plot(t,y1,'r-.',t,y2,'k');•xlabel('t/s');•ylabel('y1=sin(t), y2=cos(t)'); •title('y1=sin(t), y2=cos(t)');•text(3,0.4,'y1=sin(t)');•text(2,0,'y2=cos(t)');•legend('y1=sin(t)','y2=cos(t)'); 2024/9/1812.5.4 屏幕控制与其他2维绘图1. figure %打开图形窗口2. figure(n) %打开指定图形窗口3. close %关闭当前图形窗口4. close all %关闭所有图形窗口5. close(n) %关闭指定图形窗口6. clf %清除窗口内所有内容2024/9/1822.5.4 屏幕控制与其他2维绘图(续)8.subplot(m,n,p) %图形分为m x n个子图，并指定第p个。

排号从左到右，从上到下排号从左到右，从上到下9.stem(t,y) %绘脉冲图10.stairs(t,y) %绘阶梯图11. bar(t,y) %绘条形图12.errorbar(t,y) %绘误差条形图13.hist(y) %绘直方图14.fill(t,y,’r’) %绘填充图2024/9/183•如y=exp(-0.1t)*sin(t)•t=0:0.3:4*pi;•y=exp(-0.1*t).*sin(t);•figure(3)•plot(t,y,'k*');•figure(4)•subplot(2,2,1);stem(t,y,'k.');title('stem(t,y)');•subplot(2,2,2);stairs(t,y,'b');title('stairs(t,y)');•subplot(2,2,3);bar(t,y,'g');title('bar(t,y)');•subplot(2,2,4);fill(t,y,'r');title('fill(t,y,''r'')');2024/9/184•hist(y)•t=0:0.1:4*pi;•y=exp(-0.1*t).*sin(t);•y1=5.*y.*sin(t);•plotyy(t,y,t,y1);2024/9/1852.5.4 屏幕控制与其他2维绘图(续)15. pause %暂停16. grid on %增加网格17. grid off %取消网格18. grid %乒乓增加和取消网格19. loglog %双对数坐标log1020. semilogx %半对数坐标，x轴半对数21. semilogy %半对数坐标，y轴半对数22. polar(theta,rho) %极坐标图2024/9/1862.5.4 屏幕控制与其他2维绘图(续)23、虚数的绘图--------- Z为虚数 plot(Z) %实部为x坐标，虚部为y轴 plot(t,Z) %虚部丢失24.axis([xmin,xmax,ymin,ymax]) %定义坐标25.axis square %两轴坐标长度相等26.axis equal %两轴坐标刻度相同27.axis tight %坐标区域和图形吻合28.set(gca,’xtick’,[-1,3,7,11]) %在x轴指定处标记刻度2024/9/187clear,clft=0:2*pi/99:2*pi;x=1.15*cos(t);y=3.25*sin(t); %y为长轴，x为短轴subplot(2,3,1);plot(x,y),axis normal,grid on,title('Normal and Grid on')subplot(2,3,2);plot(x,y),axis equal,grid on,title('Equal')subplot(2,3,3);plot(x,y),axis square,grid on,title('Square')subplot(2,3,4);plot(x,y),axis image,box off,title('Image and Box off')subplot(2,3,5);plot(x,y),axis image fill,box offtitle('Image and Fill')subplot(2,3,6);plot(x,y),axis tight,box off,title('Tight') 2024/9/1882024/9/1892.5.5 三维曲线和曲面1.plot3(x,y,z,’r’); %画三维曲线Plot3(x1,y1,z1,’r’,x2,y2,z2,’b’)t=0:0.02*pi:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z,'b-',x,y,z,'bd');view([-82,58]);box onlegend('链链','宝石宝石');2024/9/1902.5.5 三维曲线和曲面（续）2. mesh(z) %画三维网格曲线,z为x,y的函数 mesh(x,y,z) %常用画三维网格曲线 mesh(x,y,z,’r’) %带颜色的三维图x=-8:0.5:8;y=x’;X=ones(size(y))*x;Y=y*ones(size(x));R=sqrt(X.*X+Y.*Y);Z=sin(R )./R;mesh(Z);2024/9/1912.5.5 三维曲线和曲面（续）3. surf(Z) %由多个小面组成表面视图 surf(x,y,z) surf(x,y,z,’r’)x=-8:0.5:8;y=x';X=ones(size(y))*x;Y=y*ones(size(x));R=abs(X)+abs(Y)+eps;Z=sin(R )./R;surf(Z);2024/9/1922.5.5 三维曲线和曲面（续）4. Meshgrid(x,y) %生成网格点坐标函数x=-4:4;y=x;[X,Y]=meshgrid(x,y);Z=X.^2+Y.^2;surf(X,Y,Z);colormap(hot) %Black-red-yellow-white颜色颜色hold onstem3(X,Y,Z,'bo');2024/9/1932.5.5 三维曲线和曲面（续）5. view(方位角，俯仰角) %改变视角6. shading flat %把曲面上的小格平滑掉 shading interp %更平滑7. rotate3d %旋转8. contour3(Z) %画等高线9. meshc, surfc %带等高线的三维作图10.colormap(hot) %hot,cool,gray,copper,pink,jet,prism11. colorbar %画彩色条12. hidden off %透视被叠压的图形 hidden on %消隐被叠压的图形2024/9/1942.5.5 三维曲线和曲面（续）clear,clf[X0,Y0,Z0]=sphere(30); %产生单位球面的三维坐标X=2*X0;Y=2*Y0;Z=2*Z0;%产生半径为2的球面的三维坐标surf(X0,Y0,Z0);%画单位球面shading interp%采用插补明暗处理hold on;mesh(X,Y,Z);colormap(hot);hold off %采用hot色图hidden off%产生透视效果axis equal,axis off%不显示坐标轴 2024/9/1952.5.5 三维曲线和曲面（续）13. Moviein , getframe, movie %动画axis equal M=moviein(16);for j=1:16 plot(fft(eye(j+16))); M(:,j)=getframe;endmovie(M,30);14. alpha(v) %透明度控制 v∈[0,1]2024/9/196•clear; clf; shg,•x=3*pi*(-1:0.05:1); y=x; • [X,Y]=meshgrid(x,y);•R=sqrt(X.^2+Y.^2)+eps; • Z=sin(R)./R;•h=surf(X,Y,Z); colormap(jet);•axis off•n=12;•mmm=moviein(n); %预设画面矩阵。

•for i=1:n•rotate(h,[0 0 1],25);%使图形绕z轴旋转25度/每次•mmm(:,i)=getframe;%捕获画面•end•movie(mmm,5,10) %以每秒10帧速度，重复播放5次2024/9/1972.6 M文件及程序调试•M文件可以分为两种：一种是主程序，一种是子程序即函数文件一个较复杂的程序往往是由这两种程序混合组成•2.6.1 主程序文件•主程序一般用clear，close all开头•程序主体•程序存盘的文件名 2024/9/1982.6.2 函数文件•把一个比较大的任务分解为多个比较小的任务，它们之间通过调用实现参数传递，小任务可以是函数•格式：•function [out1,out2,…]=函数名(in1,in2,…)1、函数调用•常见的函数调用形式为：–[out1,out2,…]=function(in1,in2,…)•一个函数可以嵌套，也可以调用其它的函数，甚至调用自己（也就是递归调用）•函数文件，函数名称和文件名一般相同2024/9/1992.6.2 函数文件• 与脚本文件不同与脚本文件不同，函数文件犹如一个，函数文件犹如一个“黑箱黑箱”，把一些数据送进并经加工处理，再把结果送，把一些数据送进并经加工处理，再把结果送出来。

出来§MATLAB提供的函数指令大部分都是由函数文提供的函数指令大部分都是由函数文件定义的件定义的§M函数文件的特点是：函数文件的特点是：• 从形式上看从形式上看，与脚本文件不同，与脚本文件不同，函数文件的，函数文件的笫一行总是以笫一行总是以 “function”引导的引导的“函数申明函数申明行行”2024/9/11002.6.2 函数文件•从运行上看从运行上看，与脚本文件运行不同，与脚本文件运行不同，每当函数，每当函数文件运行，文件运行， MATLAB就会专门为它开辟一个临就会专门为它开辟一个临时工作空间，称为时工作空间，称为函数工作空间函数工作空间（（ Function workspace））当执行文件最后一条指令时当执行文件最后一条指令时，，就结束该函数文件的运行，同时该临时函数空就结束该函数文件的运行，同时该临时函数空间及其所有的中间变量就立即被清除间及其所有的中间变量就立即被清除•MATLAB允许使用比允许使用比 “标称数目标称数目 ”较少的输较少的输入输出宗量，实现对函数的调用入输出宗量，实现对函数的调用 2024/9/11012.6.2 函数文件§由于从结构上看由于从结构上看，脚本文件只是比函数文件少，脚本文件只是比函数文件少一个一个“函数申明行函数申明行”，所以只须描述清楚函数，所以只须描述清楚函数文件的结构文件的结构。

•典型典型 M函数文件的结构如下函数文件的结构如下：：• 函数申明行：位于函数文件的首行，以关键字函数申明行：位于函数文件的首行，以关键字 functio 开头，函数名以及函数的输入输出宗量开头，函数名以及函数的输入输出宗量都在这一行被定义都在这一行被定义•笫一注释行：紧随函数申明行之后以笫一注释行：紧随函数申明行之后以%开头笫一开头笫一注释行该行供注释行该行供lookfor关键词查询和关键词查询和 help帮助使用帮助使用 2024/9/11022.6.2 函数文件•帮助文本区帮助文本区：笫一注释行及其之后的连续以：笫一注释行及其之后的连续以%开头开头的所有注释行构成整个帮助文本的所有注释行构成整个帮助文本•编写和修改记录：与帮助文本区相隔一个编写和修改记录：与帮助文本区相隔一个“空空”行，行，也以也以%开头，标志编写及修改该开头，标志编写及修改该M文件的作者和日期等文件的作者和日期等 • 函数体：为清晰起见，它与前面的注释以函数体：为清晰起见，它与前面的注释以“空空”行相行相隔2024/9/11032.6.2 函数文件•需要注意函数文件的放置位置，一般自己的函数文件放在当前目录；如果对一个专题有了足够多的函数，可以生成一个工具箱，放在一个固定的目录下，并在MATLAB中加入这个目录路径即可。

•使用函数可以加快计算速度MATLAB首次执行一个函数时，它将打开的文件编译为存储器内部形式，加速了执行速度普通的m文件不被编译，在每次编译时，文件将逐行解释执行•函数的前一部分注释为帮助行，在使用help命令是看到的为这些注释行2024/9/11042.6.2 函数文件2、参数传递•MATLAB函数的输入输出数目都可以变化，通过这个特性，可以实现一些自定义的功能函数的输入输出参数数目可以通过变量nargin和nargout获得函数调用中可以使用少于规定的输入输出参数数目，但是不能更多•在MATLAB中，参数具有自己的专有工作空间函数中的参数和命令行参数不在一个空间中，它们的唯一联系为函数的输入输出变量输入参数在函数中是可读的，但任何改动不会传递回上一级空间•使用global命令可以将变量说明为全局的，则在函数、命令行等都可以共享这些变量在实际应用中，应当尽量避免使用全局变量2024/9/11052.6.2 函数文件3、函数注意•函数可以按少于函数M 文件中所规定的输入和输出变量进行调用，但不能用多于函数M 文件中所规定的输入和输出变量数目如果输入和输出变量数目多于函数M 文件中function 语句一开始所规定的数目，则调用时自动返回一个错误。

•当调用一个函数时，所用的输入和输出的参量的数目，在函数内是规定好的函数工作空间变量nargin 包含输入参量个数；函数工作空间变量nargout 包含输出参量个数事实上，这些变量常用来设置缺省输入变量，并决定用户所希望的输出变量在M 文件函数里，变量nargout 可用来检验输出参量的个数，并按要求修正输出变量的创建2024/9/11062.6.2 函数文件•函数有它们自己的专用工作空间，它与MATLAB 的工作空间分开函数内变量与MATLAB 工作空间之间唯一的联系是函数的输入和输出变量如果函数任一输入变量值发生变化，其变化仅在函数内出现，不影响MATLAB 工作空间的变量函数内所创建的变量只驻留在函数的工作空间，而且只在函数执行期间临时存在，以后就消失因此，从一个调用到下一个调用，在函数工作空间变量存储信息是不可能的2024/9/11072.6.2 函数文件•当调用一个函数时，输入变量不会拷贝到函数的工作空间，但使它们的值在函数内可读然而，改变输入变量内的任何值，那么数组就拷贝到函数工作空间进而，按缺省，如果输出变量与输入变量相同，例如，函数x=fun(x, y, z) 中的x ，那么就将它拷贝到函数的工作空间。

因此，为了节约存储和增加速度，最好是从大数组中抽取元素，然后对它们作修正，而不是使整个数组拷贝到函数的工作空间2024/9/11082.6.2 函数文件•如果变量说明是全局的,函数可以与其它函数、MATLAB 工作空间和递归调用本身共享变量为了在函数内或MATLAB 工作空间中访问全局变量，在每一个所希望的工作空间，变量必须说明是全局的 •实际编程中，无论什么时候应尽量避免使用全局变量要是用了全局变量，建议全局变量名要长，它包含所有的大写字母，并有选择地以首次出现的M 文件的名字开头如果遵循建议，则在全局变量之间不必要的互作用减至最小2024/9/11092.6.2 函数文件• MATLAB 以搜寻脚本文件的同样方式搜寻函数M 文件例如，输入» cow ，MATLAB 首先认为cow 是一个变量如果它不是，那么MATLAB 认为它是一个内置函数如果还不是，MATLAB 检查当前cow.m 的目录或文件夹如果它不存在，MATLAB 就检查cow.m 在MATLAB 搜寻路径上的所有目录或文件夹 •从函数M 文件内可以调用脚本文件在这种情况下，脚本文件查看函数工作空间，不查看MATLAB 工作空间。

从函数M 文件内调用的脚本文件不必用调用函数编译到内存函数每调用一次，它们就被打开和解释因此，从函数M 文件内调用脚本文件减慢了函数的执行2024/9/11102.6.2 函数文件•当MATLAB 运行时，它缓存了存储在Toolbox 子目录和Toolbox 目录内的所有子目录中所有的M 文件的名字和位置这使MATLAB 很快地找到和执行函数M 文件被缓存的M 文件函数当作是只读的如果执行这些函数，以后又发生变化，MATLAB 将只执行以前编译到内存的函数，不管已改变的M 文件而且，在MATLAB 执行后，如果M 文件被加到Toolbox 目录中，那么它们将不出现在缓存里，因此不可利用所以，在M 文件函数的使用中，最好把它们存储在Toolbox 目录外，或许最好存储在MATLAB 目录下，直至它们被认为是完备的当它们是完备时，就将它们移到一个只读的Toolbox 目录或文件夹的子目录内最后，要确保MATLAB 搜索路径改变，以确认它们的存在2024/9/11112.6.2 函数文件•MATLAB 函数error 在命令命令窗口显示一个字符串，放弃函数执行，把控制权返回给键盘。

这个函数对提示函数使用不当很有用，如在以下文件片段中：–if length(val)>1–error(' VAL must be a scalar. ')–end•这里，如果变量val 不是一个标量，error 显示消息字符串，把控制权返回给命令命令窗口和键盘2024/9/1112M函数文件示例[circle.m]function sa = circle(r,s)%CIRCLE plot a circle of radii r in the line specified by s.% r指定半径的数值% s指定线色的字符串% sa圆面积%% circle(r)利用蓝实线画半径为 r 的圆周线.% circle(r,s)利用串 s 指定的线色画半径为 r 的圆周线.% sa=circle(r) 计算圆面积，并画半径为 r 的蓝色圆面.% sa=circle(r,s)计算圆面积，并画半径为 r 的 s 色圆面. % 编写于2006年4月7日，修改于2006年6月27日2024/9/1113if nargin>2 error('输入宗量太多');end;if nargin==1 s='b';end;clf;t=0:pi/100:2*pi;x=r*exp(i*t);if nargout==0 plot(x,s);else sa=pi*r*r; fill(real(x),imag(x),s)endaxis('square') 2024/9/1114第3章 MATLAB 的开发环境和工具 •3.1 MATLAB与其它软件的接口关系•3.2 MATLAB的文件管理系统•3.3 MATLAB 6.x的开发环境2024/9/11153.1 MATLAB与其它软件的接口关系•3.1.1 与磁盘操作系统的接口关系1、变量的存储和下载如：save aa a b c %将内存变量a, b, c内容以文件aa.mat的方式存储在磁盘中。

或save aa a b c –ascii %以ASCII码格式存储2024/9/11163.1.1 与磁盘操作系统的接口关系 load aa 表示将磁盘上存储的aa.mat 数据文件内容取回到工作空间即内存中此时内存中的变量与存储时的变量相同此时必须注意，原来内存中不能有与提取文件中的变量相同的变量，否则原来内存中的变量内容将被取代而丢失2024/9/11173.1.1 与磁盘操作系统的接口关系2、工作日志的记录 diary 命令可以把MATLAB工作过程中的全部屏幕文字和数据以文本方式记录下来，成为一个工作记录 diary on %默认文件名diary.txt diary bbb %文件名为bbb.txt diary off %结束记录2024/9/11183.1.1 与磁盘操作系统的接口关系3、日期和时间命令 t0=clock; %提取年月日时分秒数据并求差值 y=inv(rand(100,100)); etime(clock,t0) t=cputime; %以开机时间为基准 y=inv(rand(100,100)); cputime-t tic; %秒表置零，求经历的时间 y=inv(rand(100,100)); toc2024/9/11193.1.1 与磁盘操作系统的接口关系4、不退出MATLAB环境运行其他软件格式：！命令 2024/9/11203.1.2 与文字处理系统WINWORD 的关系•利用剪贴板进行交互•文字编辑器的使用•Notebook软件工具2024/9/11213.1.3 图形文件的转储•可以利用图形窗口中 figure copy来粘贴。

•可以利用图形窗口File菜单中的导出子菜单Export来选择需要存储的图形文件格式•可以利用图形窗口的工具对图形进行一些相关处理2024/9/11223.1.4 低层输入输出函数库•主要是实现文件内容的相互交换具体函数可以参考P55页的表3.3•如：[X，map]=imread(‘aa.bmp’,’bmp’)就是将图象数据读入X中，颜色读入map中每条命令的使用可以通过HELP查阅2024/9/11233.2 MATLAB的文件管理系统•MATLAB 自身的用户文件格式1、程序文件 .m2、数据文件.mat3、可执行文件.mex4、仿真模型文件.mdl5、仿真文件.s2024/9/11243.2 MATLAB的文件管理系统•who(whos)•cd(chdir,pwd)•dir(ls)•type•what•which•clc•edit•!•echo•load•clear/pack•save•diary2024/9/11253.3 MATLAB 6.X的开发环境•这里主要通过MATLAB 6.0的演示来获得开发环境的了解1、命令窗口 2、历史命令窗口3、资源目录本 4、当前路径浏览器5、工作空间浏览器 6、帮助浏览器7、数组编辑器 8、程序编辑器要掌握最基本的环境使用。

2024/9/1126第4章 MATLAB的其他函数库•4.1 数据分析函数库datafun•4.2 矩阵的分解与变换matfun•4.3 多项式函数库polyfun•4.4 函数功能和数值积分函数库funfun•4.5 字符串函数库strfun•4.6 稀疏矩阵函数库sparfun•4.7 图形界面函数库guitools•4.8 数据类型函数库datatypes 2024/9/11274.1 数据分析函数库datafun•corrcoef(x)求相关系数•cov(x)协方差矩阵•cplxpair(x)把向量分类为复共轭对•cross(x, y)向量的向量积•cumprod(x)列累计积•cumsum(x)列累计和•del2(A)五点离散拉氏算子•diff(x) 计算元素之间差•dot(x, y)向量的点积•gradient(Z, dx, dy)近似梯度•histogram(x)直方图和棒图•max(x), max(x, y)最大分量•mean(x)均值或列的平均值•median(x)列的中值•min(x), min(x, y)最小分量•prod(x)列元素的积•rand(x)均匀分布随机数•randn(x)正态分布随机数•sort(x)按升序排列•std(x)列的标准偏差•subspace(A, B)两个子空间之间的夹角•sum(x)各列的元素和2024/9/11284.1 数据分析函数库datafun•4.1.1 基本的数据分析1. max %求各列最大值2.min %求各列最小值 3.mean %求各列平均值4.std %求各列标准差5.median %求各列中间元素6.sum %求各列和7.trapz %梯形法求积分 8.diff %求差分9.sort %排序 2024/9/1129max(A)==[0.9501 0.8214 0.9218 0.9355]min(A)==[0.2311 0.0185 0.1763 0.4057]mean(A)=[0.6331 0.5006 0.6487 0.7124]std(A)==[0.2963 0.3197 0.2861 0.2783]median(A)==[0.6068 0.4565 0.7382 0.8936]sum(A)=[3.1654 2.5032 3.2437 3.5620]trapz(A)==[2.2447 1.8998 2.8478 2.9123]diff(A)==[-0.7190 -0.3056 0.1765 0.5298 0.3757 -0.4380 0.1299 -0.0186 -0.1209 0.8029 -0.1836 -0.5066 0.4053 -0.3767 -0.5619 0.4834]sort(A)==[0.2311 0.0185 0.1763 0.4057 0.4860 0.4447 0.6154 0.4103 0.6068 0.4565 0.7382 0.8936 0.8913 0.7621 0.7919 0.9169 0.9501 0.8214 0.9218 0.9355]2024/9/11304.1 数据分析函数库datafun•4.1.2 用于场论的数据分析函数1、gradient %求梯度2、del2 %拉普拉斯算子3、cross %求矢量积4、dot %求数量积2024/9/11314.1 数据分析函数库datafun•4.1.3 用于随机数据分析的函数1、rand(m,n)2、randn(m,n)3、hist(x) %画直方图4、hist(x,N) %N等份画直方图2024/9/1132x=rand(1,1000);plot(x)hist(x)2024/9/1133x=randn(1,1000);plot(x)hist(x); hist(x,50);2024/9/11344.1 数据分析函数库datafun•4.1.4 用于相关和傅立叶分析的函数1、corrcoef(x,y) %相关系数2、cov(x,y) %协方差 3、conv(x,y) %卷积和多项式相乘4、deconv(x,y) %卷积和多项式相除5、filter(b,a,x) %一维数据滤波6、X=fft(x,N) %快速傅立叶变换7、X=ifft(X) %逆傅立叶变换8、sound(u,s) %向量变为声音9、filter2, conv2,deconv, fft2,ifft2等2024/9/1135•x=rand(1,1000);•y=randn(1,1000); R= 1.0000 0.0094•R=corrcoef(x,y) 0.0094 1.0000 %对角线是自相关系数•Z=cov(x,y) Z = 0.0835 0.0026• 0.0026 0.8928%对角线是x和y的均方差•P=conv(x,y); %卷积 P长度为Mx+Ny-1如 x=[1 3 2 1]; y=[1 2 2];•[Q ,R]=deconv(x,y)则：Q = 1 1•R = 0 0 -2 -12024/9/1136%构造受噪声污染的信号clear； randn('state',0)； t=linspace(0,10,512);y=3*sin(5*t)-6*cos(9*t)+5*randn(size(t)); plot(t,y) %绘制幅频曲线Y=fft(y);Ts=t(2)-t(1)%时间信号的采样周期Ws=2*pi/Ts;%时间信号的采样频率Wn=Ws/2 %Nyquest频率w=linspace(0,Wn,length(t)/2);%半采样频率中相应的刻度Ya=abs(Y(1:length(t)/2));plot(w,Ya) %绘制局部放大的幅频曲线ii=find(w<=20);plot(w(ii),Ya(ii))grid,xlabel('Frequency Rad/s') 2024/9/1137%信号滤波clear,randn('state',1)ws=1000;%采样频率t=0:1/ws:0.4;x=sin(2*pi*10*t)+cos(2*pi*100*t) …+0.2*randn(size(t)); %生成带噪声的多频率信号wn=ws/2; %Nyquest频率[B,A]=butter(10,30/wn);%截止频率为30/wn的10阶ButterWorth低通滤波器y=filter(B,A,x); %进行（初值为0的）滤波处理plot(t,x,'b-',t,y,'r.','MarkerSize',10)legend('Input','Output',0) 2024/9/11384.2 矩阵的分解与变换matfun主要考虑以下函数det %求行列式 rank %求矩阵的秩inv %求逆矩阵 trace %求矩阵的迹cond %求条件数 [l,u]=lu(a) %三角分解 a=准下三角l*上三角u[q,r]=qr(a) %正交分解 a=正交方阵q*上三角阵r[u,s v]=svd(a) %奇异值分解 a=正交方阵u*对角阵s*正交方阵v[e,r]=eig(a); %特征值分解； e为特征向量，r为特征值2024/9/1139A=magic(3);B=det(A) ========= B=-360C=inv(A) ===== C = 0.1472 -0.1444 0.0639 -0.0611 0.0222 0.1056 -0.0194 0.1889 -0.1028R=rank(A) ===== R=3T=trace(A) ====== T=15CD=cond(A) ====== CD=4.3301[l,u]=lu(A)[q,r]=qr(A)[u,s,v]=svd(A)[e,r]=eig(A)2024/9/1140[l,u]=lu(A)准下三角必须交换二行才成为真下三角l = 1.0000 0 0 0.3750 0.5441 1.0000 0.5000 1.0000 0u = 8.0000 1.0000 6.0000 0 8.5000 -1.0000 0 0 5.29412024/9/1141[q,r]=qr(A)q = -0.8480 0.5223 0.0901 -0.3180 -0.3655 -0.8748 -0.4240 -0.7705 0.4760r = -9.4340 -6.2540 -8.1620 0 -8.2394 -0.9655 0 0 -4.63142024/9/1142[u,s,v]=svd(A) 其中|u|=1，|v|=-1u = -0.5774 0.7071 0.4082 -0.5774 0.0000 -0.8165 -0.5774 -0.7071 0.4082s = 15.0000 0 0 0 6.9282 0 0 0 3.4641v = -0.5774 0.4082 0.7071 -0.5774 -0.8165 -0.0000 -0.5774 0.4082 -0.70712024/9/1143[e,r]=eig(A)e = -0.5774 -0.8131 -0.3416 -0.5774 0.4714 -0.4714 -0.5774 0.3416 0.8131r = 15.0000 0 0 0 4.8990 0 0 0 -4.89902024/9/11444.3 多项式函数库polyfun1、P=roots(a)和a=poly(P) 根和求多项式 a=[a(1),a(2) ,a(3), …., a(n), a(n+1)]2、c=conv(a,b) 相当于 c(x)=a(x)*b(x)3、[q,r]=deconv(a,b) 多项式相除，q为商，r为余子式4、e=polyder(a)多项式求导5、f=polyval(a,xv) 多项式求值 x=xv2024/9/11456、P=polyfit(x,y,n)多项式拟合 n为阶次7、yi=interp1(x,y,xi,’method’) 插值函数有 linear,cubic,spline8、tc= interp2(x,y,z,wx,dy,’method’) 二维插值9、[r,p,k]=residue(b,a) 求线性微分方程的解。

2024/9/11461.求方程f(x)=x^3 + 1.1 x^2 + 0.55 x + 0.125=0的根A=[1 ,1.1 ,0.55, 0.125];p=roots(A)p = -0.5000 -0.3000 + 0.4000i -0.3000 - 0.4000i2024/9/11472.由给定根向量求多项式系数向量R=[-0.5,-0.3+0.4*i,-0.3-0.4*i];P=poly(R)PR=real(P) PPR=poly2str(PR,'x') P = 1.0000 1.1000 0.5500 0.1250PR = 1.0000 1.1000 0.5500 0.1250PPR = x^3 + 1.1 x^2 + 0.55 x + 0.125 2024/9/11483.求3阶方阵A的特征多项式•A=[11 12 13;14 15 16;17 18 19];PA=poly(A) PPA=poly2str(PA,'s') • PA = 1.0000 -45.0000 -18.0000 -0.0000PPA = s^3 - 45 s^2 - 18 s - 2.8387e-015 2024/9/11494.求的“商”及“余”多项式。

p1=conv([1,0,2], conv([1,4],[1,1]));p2=[1 0 1 1];[q,r]=deconv(p1,p2);cq=‘商多项式为商多项式为 ’;cr='余多项式为余多项式为 ';disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')]) 商多项式为 s + 5余多项式为 5 s^2 + 4 s + 32024/9/11505.求多项式f(x)=x^3 + 4 x^2 + 5 x +6的导数和f '(2)a=[1 4 5 6];e=polyder(a)fx=poly2str(e,'x')f=polyval(e,2)e = 3 8 5fx = 3 x^2 + 8 x + 5f= 332024/9/1151% 求频率特性a=[2,4,6,8],b=[3,6,9];clear jw=linspace(0,10);A=polyval(a,j*w);B=polyval(b,j*w);subplot(2,2,1);h=plot(w,abs(B./A));set(h,'linewidth',2)subplot(2,2,3);h=plot(w,angle(B./A));,set(h,'linewidth',2)w1=logspace(-1,1)F=polyval(b,j*w1)./polyval(a,j*w1);subplot(2,2,2);loglog(w1,abs(F),'linewidth',2)subplot(2,2,4);semilogx(w1,angle(F),'linewidth',2)2024/9/11526.对于给定数据对x0 , y0 ，求拟合三阶多项式，并图示拟合情况。

•x0=0:0.1:1;•y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];•n=3;•P=polyfit(x0,y0,n)•xx=0:0.01:1;•yy=polyval(P,xx);•plot(xx,yy,'-b',x0,y0,'.r','MarkerSize',20),xlabel('x') •P =• 56.6915 -87.1174 40.0070 -0.90432024/9/1153例：曲线拟合x=0:0.1:1;y=[-.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];subplot(2,3,1),plot(x,y,'o');title('原始数据点');a1=polyfit(x,y,1);xi=linspace(0,1);y1=polyval(a1,xi);subplot(2,3,2),plot(x,y,'o',xi,y1,'r'); title('线性拟合');a2=polyfit(x,y,2);y2=polyval(a2,xi);subplot(2,3,3),plot(x,y,'o',xi,y2,'r');title('二次拟合'); a3=polyfit(x,y,3);y3=polyval(a3,xi);subplot(2,3,4),plot(x,y,'o',xi,y3,'r'); title('三次拟合'); a9=polyfit(x,y,9);y9=polyval(a9,xi);subplot(2,3,5),plot(x,y,'o',xi,y9,'r'); title('九次拟合'); a10=polyfit(x,y,10);y10=polyval(a10,xi);subplot(2,3,6),plot(x,y,'o',xi,y10,'r') ;title('十次拟合'); 2024/9/11542024/9/11557.以上例所给数据，研究一维插值，并观察插值与拟合的区别。

•x0=0:0.1:1;•y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22]; •xi=0:0.02:1;•yi=interp1(x0,y0,xi,'cubic'); •plot(xi,yi,'-b',x0,y0,'.r','MarkerSize',20),xlabel('x') 2024/9/11569.部分分式展开a=[1,3,4,2,7,2];b=[3,2,5,4,6];[r,s,k]=residue(b,a) •r = 1.1274 + 1.1513i ； 1.1274 - 1.1513i；• -0.0232 - 0.0722i； -0.0232 + 0.0722i；• 0.7916 •s = -1.7680 + 1.2673i； -1.7680 - 1.2673i；• 0.4176 + 1.1130i； 0.4176 - 1.1130i；• -0.2991 •k = [] 2024/9/1157求线性微分方程在输入u(t)为单位脉冲及单位阶跃信号时的响应a=[1,5,4,7];b=[3,0.5,4];[r,p,k]=residue(b,a)t=0:0.2:10;yi=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t);subplot(1,2,1),h=plot(t,yi);set(h,'linewidth',2)a=[1,5,4,7,0];b=[3,0.5,4];[r,p,k]=residue(b,a);ys=0;for i=1:length(r) ys=ys+r(i)*exp(p(i)*t);endsubplot(1,2,2), plot(t,ys,'linewidth',2)2024/9/1158响应结果图2024/9/11594.4 函数功能和数值积分函数库funfun•%绘制函数曲线•fplot(‘函数名’，[初值x0，终值xf]）•fplot(FUN,LIMS,TOL,N,'LineSpec',P1,P2,...) function [y]=jszero(t,a,b) y=sin(t).^2.*exp(-a.*t)-b.*abs(t);fplot('jszero',[-10,10],[],[],[],0.1,0.5)2024/9/1160subplot(2,2,1), fplot(@humps,[0 1])•f = inline('abs(exp(-j*x*(0:9))*ones(10,1))');•subplot(2,2,2), fplot(f,[0 2*pi])subplot(2,2,3), fplot('[tan(x),sin(x),cos(x)]',2*pi*[-1 1 -1 1])•subplot(2,2,4), fplot('sin(1 ./ x)', [0.01 0.1],1e-3)2024/9/11612024/9/1162%求函数极值fmin(‘函数名’，初值x0，终值xf） X = fmin('F',x1,x2,OPTIONS,P1,P2,...) 如：m=fmin(‘humps’,[0,1.5]);%求函数零点fzero(‘函数名’,初猜值x0)X = fzero(@sin,3) X = fzero(inline('sin(3*x)'),2);2024/9/1163求的零点 P1=0.1;P2=0.5;%用字符串表示被处理函数用字符串表示被处理函数yC='sin(x).^2.*exp(-P1*x)-P2*abs(x)'; x=-10:0.01:10; %作图观察函数零点分布作图观察函数零点分布Y=eval(yC); %对字符串表达式进行计算并输出对字符串表达式进行计算并输出clf,plot(x,Y,‘r’); hold onplot(x,zeros(size(x)),'k');xlabel('t');ylabel('y(t)'),hold off 2024/9/1164%利用zoom和ginput 指令获得零点的初始近似值zoom on[tt,yy]=ginput(5);zoom off%计算tt(4)附近的零点[t4,y4]=fzero(yC,tt(4),[ ],P1,P2);for i=1:length(tt) [xup(i),yup(i)]=fzero(yC,tt(i),[ ],P1,P2)end2024/9/1165或[z,fz,exitflag]=fzero(fun,x0,options,p1,p2,..)Exitflag>0找到零点后退出Options为优化迭代所采用的参数选项function [y]=jszero(t,a,b) y=sin(t).^2.*exp(-a.*t)-b.*abs(t); [z,fz]=fzero('jszero',1.6,[],0.1,0.5)2024/9/1166quad() 采用递推自适应Simpson法quad(‘函数名’，初值x0，终值xf）求定积分quad(fun,a,b,TOL,TRACE,P1,P2,...) tol 控制绝对误差,trace 取非零将逐点画。

F = inline('1./(x.^3-2*x-5)'); Q = quad(F,0,2);Q = quad(@myfun,0,2); function y = myfun(x) y = 1./(x.^3-2*x-5);2024/9/1167quadl()采用递推自适应Lobatto法 F = inline('1./(x.^3-2*x-5)'); Q = quadl(F,0,2); Q = quadl(@myfun,0,2); function y = myfun(x) y = 1./(x.^3-2*x-5); 2024/9/1168ode45采用4，5阶龙格库塔法求解微分方程[x,y]=ode45(‘函数名’，自变量初值x0，自变量终值xf，因变量初值x0); [t,y]=ode45(@vdp1,[0 20],[2; 0]); plot(t,y(:,1));figure; plot(y(:,1),y(:,2));function [yp]=vdp1(t,y) r=2; yp(1)=y(2); yp(2)=-r*(y(1).^2-1)*y(2)-y(1); yp=yp'; % yp=[y(2);-r*(y(1).^2-1)*y(2)-y(1)];2024/9/1169[t,y]=ode45('vdp2',[0, 30],[1;0;0;1 ])plot(t,y(:,1))function [yp]=vdp2(t,y) yp(1)=y(2); yp(2)=y(3); yp(3)=y(4);; yp(4)=-y(4)-(y(1)^2+y(1))*y(3)-(y(1)+1)*y(2)-y(1); yp=yp'; %yp=[y(2);y(3);y(4);-y(4)-(y(1)^2+y(1))*y(3)-(y(1)+1)*y(2)-y(1)]2024/9/1170•y(1)=x;y(2)=dx/dt;•function ydot=DyDt(t,x) mu=2; ydot=[x(2);mu*(1-x(1)^2)*x(2)-x(1)];[tt,yy]=ode45(@DyDt,[0,30],[1;0]);plot(tt,yy(:,1));title(‘x(t)’);2024/9/11714.5 字符串函数库strfun•eval(string) 作为一个MATLAB 命令求字符串的值•eval(try,catch)•blanks(n) 返回一个n 个零或空格的字符串•deblank 去掉字符串中后拖的空格•feval 求由字符串给定的函数值•findstr 从一个字符串内找出字符串•isletter 字母存在时返回真值2024/9/1172•isspace 空格字符存在时返回真值•isstr 输入是一个字符串，返回真值•lasterr 返回上一个所产生MATLAB 错误的字符串•strcmp 字符串相同，返回真值•strrep 用一个字符串替换另一个字符串•strtok 在一个字符串里找出第一个标记2024/9/1173eval的使用P1=0.1;P2=0.5;y='sin(x).^2.*exp(-P1*x)-P2*abs(x)'; x=-10:0.01:10;Y=eval(y);plot(x,Y);grid on2024/9/11744.5 字符串函数库strfun•字符串转换•abs字符串到ASCII转换•dec2hex十进制数到十六进制字符串转换•fprintf把格式化的文本写到文件中或显示屏上•hex2dec十六进制字符串转换成十进制数•hex2num十六进制字符串转换成IEEE浮点数•int2str整数转换成字符串•lower字符串转换成小写2024/9/1175•num2str数字转换成字符串•setstr ASCII转换成字符串•sprintf用格式控制，数字转换成字符串•sscanf用格式控制，字符串转换成数字•str2mat字符串转换成一个文本矩阵•str2num字符串转换成数字•upper字符串转换成大写2024/9/1176串数组的属性和标识 •1、创建字符串数组、创建字符串数组---用单引号对用单引号对A=‘This is an example.’•2、字符串大小、字符串大小---字符个数字符个数[m,n]=size(A)•3、串数组的元素标识、串数组的元素标识A14=A(1:4) %提取子字符串ra=A(end:-1:1) %字符串倒排2024/9/1177•4、串数组的、串数组的ASCII码码ascii_a=double(A); %产生ASCII码char(ascii_a) %把ASCII码变回字符串•5、对字符串、对字符串ASCII码数组的操作码数组的操作W=find(A>=‘a’&A<=‘z’);ascii_a(W)=ascii_a(W)-32;char(ascii_a) %将小写变大写2024/9/1178•6、中文字符串数组、中文字符串数组A=‘这是一个算式’；A_s=size(A);A56=A([5 6]);ASCII_a=double(A);char(ASCII_a);2024/9/1179•7、创建带单引号的字符串、创建带单引号的字符串b=‘Example ’’3.1.2_1’’’;•8、由小串构成长串、由小串构成长串Ab=[a(1:7),’ ’,b,’.’];Ab=This is Example ‘3.1.2._1’.2024/9/1180复杂串数组的创建 1、多行串数组的直接创建•clear•S=['This string array '• 'has multiple rows.'] •但长度必须相等但长度必须相等2024/9/1181•2、由专门函数char , str2mat , strvcat创建多行串数组 S1=char('This string array','has two rows.') S1 =This string arrayhas two rows. 2024/9/1182•S2=str2mat('这这','字字符符','串串数数组组','由由4行行组组成成') •S2 =•这 •字符 •串数组 •由4行组成 2024/9/1183•S3=strvcat('这这','字字符符','串串数数组组',' ','由由4行行组组成成') •S3 =•这 •字符 •串数组 • •由4行组成 •size(S3) •ans =• 5 52024/9/1184转换函数产生数码字符串 •最常用的数组/字符串转换函数int2str , num2str , mat2str A=eye(2,4);A_str1=int2str(A) %非整数四舍五入非整数四舍五入 rand('state',0)B=rand(2,4); B3=num2str(B,3) %将非整数转换为串将非整数转换为串 2024/9/1185•把数值数组转化成输入形态的串数组把数值数组转化成输入形态的串数组•B_str=mat2str(B,4) •Expression=['exp(-',B_str,')'];•eval(Expression)2024/9/1186在MATLAB计算生成的图形上标出图名和最大值点坐标。

•Y=exp(-2t)*sin(3t)• 2024/9/1187clear;a=2;w=3;t=0:0.01:10;y=exp(-a*t).*sin(w*t);[ymax,imax]=max(y);ttext=['t=',num2str(t(imax))];ytext=['y=',num2str(ymax)];maxtext=char('maximum',ttext,ytext);tit=['y=exp(-',num2str(a),'t)*sin(',num2str(w),'t)'];plot(t,zeros(size(t)),'k');hold onplot(t,y,'b')plot(t(imax),ymax,'r.','MarkerSize',20)text(t(imax)+0.3,ymax+0.05,maxtext)title(tit),xlabel('t'),ylabel('y'),hold off2024/9/1188fprintf, sprintf, sscanf的用法示例•rand('state',0);a=rand(2,2);•s1=num2str(a)•s_s=sprintf('%.10e\n',a) •fprintf('%.5g\\',a) •0.95013\0.23114\0.60684\0.48598\ •s_sscan=sscanf(s_s,'%f',[3,2]) •s_sscan = %浮点格式把串转换成3X2数值数组• 0.9501 0.4860• 0.2311 0• 0.6068 0 2024/9/1189符号计算符号计算•f=sym(arg) % 将arg转化为符号对象。

•syms arg2 arg3 arg4如y=sym(‘2*sin(x)*cos(x)’)yy=simple(y); %化为最简形式如syms xy=2*sin(x)*cos(x)yy=simple(y); 2024/9/1190合并同类项z=sym(‘3*x^2*y+2*x*(y-2)+2*x^2-x*(y-2)’)z = 3*x^2*y+2*x*(y-2)+2*x^2-x*(y-2)p=collect(z)p = (3*y+2)*x^2+x*(y-2)2024/9/1191求积分求积分 •int(P) 对表达式P进行不定积分• int(P,v) 以v为积分变量对P进行不定积分• int(P,v,a,b) 以v为积分变量,以a为下限,b为上限对P进行定积分2024/9/1192•syms a1 a2 a3 a4•A=[a1 a2;a3 a4];•DA=det(A);IA=inv(A);EA=eig(A);•求定积分 int(f,x,t1,t2)•syms A t tao w•yf=int(A*exp(-i*w*t),t,-tao/2,tao/2);•yf=simple(yf)2024/9/1193置换操作•[RS]=subs(S,old, new);•[RS]=subs(S, new);•sym a x;•f=a*sin(x)+5; f1=subs(f,’sin(x)’,sym(‘y’));•f2=subs(f,{a,x},{2,sym(pi/3)})•f3=int(f,a);•f4=int(f,x,0,2);2024/9/1194求•syms x y z•f2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2)•vf2=vpa(f2) %指定精度•vpa(x,n) %指定n位相对精度2024/9/1195求导数求导数•diff(S,v) 求表达式S对变量v的一阶导数。

•diff(S,v,n) 求表达式S对变量v的n阶导数如求sin(x)+ex的三阶导数，键入命令 diff('sin(x)+x*exp(x)',3) 得 ans = -cos(x)+3*exp(x)+x*exp(x) 2024/9/1196•sym a t x; f=[a,t^3;t*cos(x);log(x)];•df=diff(f); dfdt2=diff(f,t,2); dfdxdt=diff(diff(f,x),t);2024/9/1197求和求和•symsum(S) 对通项S求和，其中k为变量且从0变到k-1•symsum(S,v) 对通项S求和，指定其中v为变量且v从0变到v-1•symsum(S,a,b) 对通项S求和，其中k为变量且从a变到b•symsum(S,v,a,b) 对通项S求和，指定其中v为变量且v从a变到b2024/9/1198例：键入k=sym(‘k’); symsum(k) 得 • ans = 1/2*k^2-1/2*k 又例如：键入 symsum(k^2,0,10)得 • ans = 385 又例如：键入•symsum('x'^k/sym('k!'),k,0,inf)得 • ans = exp(x) • 这最后的一个例子是无穷项求和。

2024/9/1199求极限求极限•limit(P) 表达式P中自变量趋于零时的极限•limit(P,a) 表达式P中自变量趋于a时的极限•limit(P,x,a,'left') 表达式P中自变量x趋于a时的左极限•limit(P,x,a,'right') 表达式P中自变量x趋于a时的右极限2024/9/1200例如：键入 P=sym('sin(x)/x'); limit(P) 得• ans = 1 键入 P=sym('1/x'); limit(P,x,0,'right') 得• ans = inf 键入 P=sym('(sin(x+h)-sin(x))/h');h=sym('h');•limit(P,h,0) 得 ans = cos(x) 键入 v=sym('[(1+a/x)^x,exp(-x)]');• limit(v,x,inf,'left') 得• ans = [ exp(a), 0] 2024/9/1201线性方程组的求解线性方程组的求解线性方程组的形式为A*X=B；其中A至少行满秩•X=linsolve(A,B) 输出方程的特解X。

•例如：键入•A=sym('[cos(t),sin(t);sin(t),cos(t)]');•B=sym('[1;1]');•c=linsolve(A,B) •c =•[ 1/(sin(t)+cos(t))]•[ 1/(sin(t)+cos(t))] 2024/9/1202代数方程的求解• solve(P,v) 对方程P中的指定变量v求解v可省略• solve(P1,P2,…，Pn,v1,v2,…，vn) •对方程P1，P2，…Pn中的指定变量v1, v2…vn求解2024/9/1203•例：可输入 solve('p+sin(x)=r') 得：•ans =-asin(p-r) •例输入：P1='x^2+x*y+y=3';P2='x^2-4*x+3=0'; • [x,y]=solve(P1,P2) 得：• x = [ 1; 3] y = [1; -3/2] •可输入：P1='a+u^2+v^2=0';P2='u-v=1';•[u,v]=solve(P1,P2,'u','v') 得：•u = [ 1/2+1/2*(-1-2*a)^(1/2); 1/2-1/2*(-1-2*a)^(1/2)]•v =[ -1/2+1/2*(-1-2*a)^(1/2) ;-1/2-1/2*(-1-2*a)^(1/2)]2024/9/1204•对于有些无法求出解析解的非线性方程组，MATLAB只给出一个数值解。

这一点可以从表示解的数字不被方括号括住而确定例如：键入：•[x,y]=solve('sin(x+y)-exp(x)*y=0','x^2-y=2') 得：• x = -6.0173272500593065641097297117905• y = 34.208227234306296508646214438330 • 由于这两个数字没有被[ ]括住，所以它们是数值解2024/9/1205•另外，可利用solve来解线性方程组的通解例如：键入• P1='2*x1+7*x2+3*x3+x4=6';• P2='3*x1+5*x2+2*x3+2*x4=4';• P3='9*x1+4*x2+x3+7*x4=2';• u=solve(P1,P2,P3,'x1','x2','x3','x4') • Warning: 3 equations in 4 variables.• u = x1: [1x1 sym] x2: [1x1 sym] x3: [1x1 sym]•x4: [1x1 sym] 2024/9/1206•可以看到：屏幕提示“有3个方程4个变量”，意为解不唯一有时会提示解不唯一）且输出的是解的结构形式。

为进一步得到解，可输入：• u.x1,u.x2,u.x3,u.x4, 得：• ans = x1• ans = -5*x1-4*x4• ans = 11*x1+9*x4+2• ans = x4 •这样就得到了原方程组的通解2024/9/1207 解符号微分方程解符号微分方程 dsolve('eq1','eq2',…) 其中eq表示相互独立的常微分方程、初始条件或指定的自变量默认的自变量为t如果输入的初始条件少于方程的个数，则在输出结果中出现常数c1,c2,等字符关于微分方程的表达式有如下的约定：字母y表式函数，Dy表示y对t的一阶导数；• Dny表示y对t的n阶导数2024/9/1208•例如：求 • •的解。

可键入：•[x,y]=dsolve('Dx=y','Dy=-x') 得x = cos(t)*C1+sin(t)*C2 y = -sin(t)*C1+cos(t)*C2 2024/9/1209例如求 , , f(0)=0 , g(0)=1 的解可输入指令：P='Df=3*f+4*g,Dg=-4*f+3*g';v='f(0)=0,g(0)=1';[f,g]=dsolve(P,v) f = exp(3*t)*sin(4*t)g = exp(3*t)*cos(4*t)•注意：微分方程表达式中字母D必须大写2024/9/1210求解微分方程 y=dsolve('D3y=-y','y(0)=1,Dy(0)=0,D2y(0)=0','x') 得: y =(1/3+2/3*exp(1/2*x)*cos(1/2*3^(1/2)*x)*exp(x))/exp(x) 2024/9/1211最后看一个解非线性微分方程的例子： dsolve('(Dy)^2+y^2=1','y(0)=0','x') ans = [ sin(x)] [ -sin(x)] 对于无法求出解析解的非线性微分方程，屏幕将提示出错信息。

2024/9/1212泰勒（taylor）级数的展开格式1：R=taylor(f)%默认x为变量的n-1阶麦克劳林f的多项式格式2：R=taylor(f,n,v)%对符号变量为v的n-1阶麦克劳林f的多项式2024/9/1213syms xf=exp(-x);ff=taylor(f)结果：ff = 1-x+1/2*x^2-1/6*x^3+1/24*x^4-1/120*x^52024/9/1214syms xf=log(x);ff=taylor(f,6,1)结果：ff = x-1-1/2*(x-1)^2+1/3*(x-1)^3-1/4*(x-1)^4 +1/5*(x-1)^52024/9/1215syms xf=sin(x);ff=taylor(f,6,pi/2)结果： ff = 1-1/2*(x-1/2*pi)^2+1/24*(x-1/2*pi)^42024/9/1216syms x tf=x^t;ff=taylor(f,3,t) %展开到t的二次方项结果：ff = 1+log(x)*t+1/2*log(x)^2*t^22024/9/12174.6 稀疏矩阵函数库sparfun2024/9/12184.7 图形界面函数库guitools例：对于传递函数为的归一化二阶系统，制作一个能绘制该系统单位阶跃响应的图形用户界面。

本例演示：（A）图形界面的大致生成过程；（B）静态文本和编辑框的生成；（C）坐标网格控制键的形成；（D）如何使用该界面2024/9/1219（（1）产生图形窗和轴位框）产生图形窗和轴位框： clf resetH=axes('unit','normalized','position',[0,0,1,1],'visible','off');set(gcf,'currentaxes',H);str='\fontname{隶书}归一化二阶系统的阶跃响应曲线';text(0.12,0.93,str,'fontsize',13);h_fig=get(H,'parent');set(h_fig,'unit','normalized','position',[0.1,0.2,0.7,0.4]);h_axes=axes('parent',h_fig,... 'unit','normalized','position',[0.1,0.15,0.55,0.7],... 'xlim',[0 15],'ylim',[0 1.8],'fontsize',8);2024/9/12202024/9/1221（（2）在坐标框右侧生成作解释用的）在坐标框右侧生成作解释用的““静态静态文本文本””和可接收输入的和可接收输入的““编辑框编辑框”” h_text=uicontrol(h_fig,'style','text','unit','normalized','position',[0.67,0.73,0.25,0.14]，'horizontal','left','string',{'输入阻尼比系数','zeta ='});h_edit=uicontrol(h_fig,'style','edit','unit','normalized','position',[0.67,0.59,0.25,0.14], 'horizontal','left', 'callback',[ 'z=str2num(get(gcbo,''string'')); 't=0:0.1:15;','for k=1:length(z);','y(:,k)=step(1,[1 2*z(k) 1],t);', 'plot(t,y(:,k));', 'if (length(z)>1) ,hold on,end,', 'end;','hold off,']); 2024/9/12222024/9/1223（（3）形成坐标网格控制按键）形成坐标网格控制按键 h_push1=uicontrol(h_fig,'style','push',... 'unit','normalized','position',[0.67,0.37,0.12,0.15],... 'string','grid on','callback','grid on');h_push2=uicontrol(h_fig,'style','push',... 'unit','normalized','position',[0.67,0.15,0.12,0.15],... 'string','grid off','callback','grid off');2024/9/12242024/9/1225（（4）输入阻尼比系数，可得单位阶跃响应曲线）输入阻尼比系数，可得单位阶跃响应曲线 2024/9/1226利用控制对象面板设计1在Matlab中选中File----New----GUI2选中Blank GUI(Default)3选中Axes作坐标，并存盘4选Static Text并在图中选择，双击写字符串5选择Edit编辑文本框，Tag取名fdisp,同时编写fdisp.m 并在Callback调用fdisp.m6增加Grid on 和Grid off按钮2024/9/1227%fdisp.m脚本文件clear all% z=str2num(get(findobj(gcbf,'tag','fdisp'),'string'));z=str2num(get(gcbo,'string'))t=0:0.1:15;for k=1:length(z) y(:,k)=step(1,[1 2*z(k) 1],t); plot(t,y(:,k)); if (length(z)>1) hold on endendhold off2024/9/12282024/9/1229实例分析1先创建一个图形窗口，键入figure，并保存为ex01.fig2在Matlab菜单中选中File----New----GUI建立一个新的控制对象面板3从控制面板菜单中选File-New-Open--ex014选中Axes作坐标5添加编辑文本框，Tag取名ex_edit,同时编写ex_edit.m 并在Callback调用ex_edit.m2024/9/12306增加Load, Apply, Grid on ,Grid off ,Close按钮，并增加弹出菜单（Popup Mean）;增加两个静态文本Input data和Color choose7加入popmenu,string为spring,summer, autumn, winter, callback和tag为ex_color,8选mean editor—new menu---label和tag输ex_view,再增加view_2D,view_3D,Callback为view(2),view(3)2024/9/1231ex_edit.m ex_load.m%ex_edit.mfunction ex_editct=get(findobj(gcbf,'tag','ex_edit‘),'string');save mydata ct%ex_load.mfunction ex_loadload mydataset(gcbf,'userdata',ct);2024/9/1232ex_apply%ex_applyfunction ex_applyct=get(gcbf,'userdata');eval(ct);plot(x,y)%Closeclose(gcbf)2024/9/1233%ex_color%ex_color function ex_colorpopstr={'spring','summer','autumn','winter'};vpop=get(findobj(gcf,'tag','ex_color'), 'value');colormap(eval(popstr{vpop}));2024/9/12342024/9/12354.8 数据类型函数库datatypes元胞数组元胞数组C_str=char('这是这是','元胞数组创建算例元胞数组创建算例 1');R=reshape(1:9,3,3);Cn=[1+2i];S_sym=sym('sin(-3*t)*exp(-t)') A(1,1)={C_str};A(1,2)={R};A(2,1)={Cn};A(2,2)={S_sym};A={C_str,R,Cn,S_sym}A2024/9/1236结果：A= [2x10 char] [3x3 double] [1.0000+ 2.0000i] [1x1 sym]显示内容：A{1,1}A{1,2}A{1,4}2024/9/1237元胞数组在存放和操作字符串上的应用a='MATLAB 5 ';b='introduces new data types:';c1='◆◆Multidimensional array';c2='◆◆User-definable data structure';c3='◆◆Cell arrays';c4='◆◆Character array';c=char(c1,c2,c3,c4);C={a;b;c};adisp([C{1:2}]); disp(' ')disp(C{3}) MATLAB 5 introduces new data types: ◆Multidimensional array ◆User-definable data structure◆Cell arrays ◆Character array 2024/9/1238构架数组本例通过温室数据（包括温室名、容积、温度、湿度等）演示：单构架的创建和显示。

green_house.name='一号房一号房';green_house.volume='2000立方米立方米';green_house.parameter.temperature=[31.2 30.4 31.6 28.7 29.7 31.1 30.9 29.6];green_house.parameter.humidity=[62.1 59.5 57.7 61.5 62.0 61.9 59.2 57.5];2024/9/1239green_house green_house = name: '一号房' volume: '2000立方米' parameter: [1x1 struct] green_house.parameter ans = temperature: [2x4 double] humidity: [2x4 double] 2024/9/1240green_house.parameter.temperature ans = 31.2000 30.4000 31.6000 28.7000 29.7000 31.1000 30.9000 29.6000 2024/9/1241利用构造函数创建构架数组 a=cell(2,3);green_house_1=struct('name',a,'volume',a,'parameter',a(1,2)) green_house_1 = 2x3 struct array with fields: name volume parameter 2024/9/1242利用构造函数创建构架数组green_house_2=struct('name',a,'volume',[],'parameter',[]) green_house_2 = 2x3 struct array with fields: name volume parameter 2024/9/1243利用构造函数创建构架数组green_hopuse_3(2,3)=struct('name',[],'volume',[],'parameter',[]) green_hopuse_3 = 2x3 struct array with fields: name volume parameter 2024/9/1244利用构造函数创建构架数组a1={'六号房六号房'};a2={'3200立方米立方米'};green_house_4(2,3)=struct('name',a1,'volume',a2,'parameter',[]);T6=[31.2,30.4,31.6,28.7;29.7,31.1,30.9,29.6];green_house_4(2,3).parameter.temperature=T6;green_house_4 ans = 2x3 struct array with fields: name volume parameter 2024/9/1245第五章第五章 MATLAB的的SIMULINK仿真仿真§SIMULINK是一个进行动态系统建模、仿真和综合分是一个进行动态系统建模、仿真和综合分析的集成软件包。

它可以处理的系统包括：线性、非线析的集成软件包它可以处理的系统包括：线性、非线性系统；离散、连续及混合系统；单任务、多任务离散性系统；离散、连续及混合系统；单任务、多任务离散事件系统事件系统•在在 SIMULINK 提供的图形用户界面提供的图形用户界面GUI上，只要进行上，只要进行鼠标的简单拖拉操作就可构造出复杂的仿真模型它外鼠标的简单拖拉操作就可构造出复杂的仿真模型它外表以方块图形式呈现，且采用分层结构表以方块图形式呈现，且采用分层结构•从建模角度讲，这既适于自上而下（从建模角度讲，这既适于自上而下（Top-down））的设的设计流程（概念、功能、系统、子系统、直至器件），又计流程（概念、功能、系统、子系统、直至器件），又适于自下而上（适于自下而上（Bottum-up））逆程设计逆程设计2024/9/1246MATLAB的的SIMULINK仿真仿真•从分析研究角度讲，这种 SIMULINK 模型不仅能让用户知道具体环节的动态细节，而且能让用户清晰地了解各器件、各子系统、各系统间的信息交换，掌握各部分之间的交互影响§在 SIMULINK 环境中，用户将观察到现实世界中非线性因素和各种随机因素对系统行为的影响。

§在 SIMULINK 环境中，用户可以在仿真进程中改变感兴趣的参数，实时地观察系统行为的变化2024/9/1247MATLAB的的SIMULINK仿真仿真§在在MATLAB 6.0 版版中，可直接在中，可直接在 SIMULINK 环境中运环境中运作的工具包很多，已覆盖通信、控制、信号处理、作的工具包很多，已覆盖通信、控制、信号处理、DSP、、电力系统等诸多领域，所涉内容专业性极强电力系统等诸多领域，所涉内容专业性极强§本讲由浅入深地讲述本讲由浅入深地讲述 SIMULINK 对各种数学、工程问对各种数学、工程问题的建模、仿真和分析的基本方法，采用题的建模、仿真和分析的基本方法，采用“算例算例”作为作为主体，配以适量的归纳性表述主体，配以适量的归纳性表述2024/9/1248模型的创建和模型文件模型的创建和模型文件SIMULINK 模型是什么？模型是什么？SIMULINK 模型有以下几层含义：模型有以下几层含义：•在视觉上表现为直观的方框图；在视觉上表现为直观的方框图；•在文件上则是扩展名为在文件上则是扩展名为 mdl 的的ASCII代码；代码；•在数学上表现为一组微分方程或差分方程；在数学上表现为一组微分方程或差分方程；• 在行为上则模拟了实际系统的动态特性在行为上则模拟了实际系统的动态特性。

2024/9/1249§SIMULINK SIMULINK 模型通常包含三种模型通常包含三种 “ “组件组件””：：•信源（信源（ SourcesSources）：）：可以是常数、时钟、白噪声、正弦可以是常数、时钟、白噪声、正弦波、阶梯波、扫频信号、脉冲生成器、随机数产生器等波、阶梯波、扫频信号、脉冲生成器、随机数产生器等信号源；信号源；• 系统（系统（ SystemSystem）：）：即指被研究系统的即指被研究系统的 SIMULINK SIMULINK 方框方框图；图；• 信宿信宿（（ S Sinkink）：）：可以是示波器、图形记录仪等可以是示波器、图形记录仪等§对于具体的对于具体的 SIMULINK SIMULINK 模型而，不一定完全地包含这三模型而，不一定完全地包含这三大组件例如：研究初始条件对系统影响就不必包含信大组件例如：研究初始条件对系统影响就不必包含信源组件2024/9/1250SIMULINK 模型的创建模型的创建•创建模型文件；创建模型文件；•选择对象；选择对象；• 模块的操作；模块的操作；• 连线的操作；连线的操作；• 对模型的注释；对模型的注释；• 创建子系统；创建子系统；• 仿真的配置仿真的配置；；§ 保存模型；保存模型；•仿真和结果分析。

仿真和结果分析2024/9/1251仿真运行仿真运行•使用菜单进行仿真使用菜单进行仿真§设置仿真参数和选择求解器设置仿真参数和选择求解器•通过选择菜单通过选择菜单 Simulation 下的下的 Parameters 菜单项，用菜单项，用来设置仿真参数和选择求解器其中有三个页面管理这来设置仿真参数和选择求解器其中有三个页面管理这些仿真参数些仿真参数•在在 Solver 页面，设置开始和停止时间，选择求解器和页面，设置开始和停止时间，选择求解器和指定求解器（指定求解器（solver））的参数，另外还可以选择一些输的参数，另外还可以选择一些输出选项•在在 Solver options 中，中， SIMULINK 模型的仿真涉及到模型的仿真涉及到一组常微分方程（一组常微分方程（ODEs））的数值积分如果模型是连的数值积分如果模型是连续系统，使用续系统，使用ode45方法方法；；如果模型不是连续系统，使如果模型不是连续系统，使用用discrete方法2024/9/1252仿真运行仿真运行•在在 Workspace I/O Workspace I/O 页面，管理对页面，管理对 MATLAB MATLAB 工作空间的工作空间的输入和输出。

输入和输出 •在在 Diagnostics Diagnostics 页面，可以选择在仿真期间显示的警页面，可以选择在仿真期间显示的警告信息的层次告信息的层次•通过命令行运行仿真通过命令行运行仿真§通过命令行运行仿真与通过菜单运行仿真相比通过命令行运行仿真与通过菜单运行仿真相比，有如，有如下的下的一些优点：一些优点：•可以不理睬模块中的初始条件（参数可以不理睬模块中的初始条件（参数 x0 x0 ）；）；•可以定义任何外部输入（用参数可以定义任何外部输入（用参数 utut ）；）；2024/9/1253仿真运行仿真运行•可以由一个可以由一个M 文件来启动一个仿真，并且允许模块中文件来启动一个仿真，并且允许模块中的参数发生改变的参数发生改变 •用来进行仿真的命令有四个：用来进行仿真的命令有四个：• 使用使用 set_param 命令：开始、停止或者继续仿真或者命令：开始、停止或者继续仿真或者更新模块的方框图更新模块的方框图 get_param 命令来检查一个仿真命令来检查一个仿真的状态• 使用使用 sim 命令：启动仿真命令；命令：启动仿真命令；• 使用使用 simset 命令：用来向命令：用来向 sim 命令产生或者编辑仿真命令产生或者编辑仿真参数和积分法属性的命令；参数和积分法属性的命令；• 使用使用 simget 命令：可以得到选项结构体属性和参数。

命令：可以得到选项结构体属性和参数2024/9/1254SIMULINK模型窗的组成 L 工具条工具条：最左边9个图标实现标准的Windows操作其余图标含义如下：打开库浏览器模型浏览器单双窗外形切换展现当前系统的父系统打开调试器仿真的启动或继续暂停（在仿真执行过程中出现）结束仿真显示库连接观察封装子系统2024/9/1255如何调用MATLAB工作空间中的信号矩阵作为模型输入本例所需的输入为（（1）编写一个产生信号矩阵的）编写一个产生信号矩阵的M函数文件函数文件function TU=source82_1(T0,N0,K)t=linspace(0,K*T0,K*N0+1);N=length(t);u1=t(1:(N0+1)).^2;u2=(t((N0+2):(2*N0+1))-2*T0).^2;u3(1:(N-(2*N0+2)+1))=0;u=[u1,u2,u3];TU=[t',u'];2024/9/1256（（2）构造简单的接收信号用的实验模型）构造简单的接收信号用的实验模型 2024/9/1257（（3）模块的参数设置）模块的参数设置双击S82_1，在对话框中的Data中填写TU（（4））在在指指令令窗窗中中，，运运行行以以下下指指令令，，在在MATLAB工作空间中产生工作空间中产生TU信号矩阵信号矩阵。

TU=source82_1(1,100,4);（（5））选选中中模模型型窗窗菜菜单单【【Simulation：：Start】】，，示波器呈现图右图信号示波器呈现图右图信号2024/9/1258复位积分器的功用示例 2024/9/1259从实际抽象出初始状态为0的二阶微分方程，是单位阶跃函数本例演示如何用积分器直接构搭求解该微分方程的模型 2024/9/1260利用存放在MATLAB工作空间中的仿真数据所绘制的曲线Clftt=ScopeData.time;%为书写简单，把构架域的时间数据另赋给为书写简单，把构架域的时间数据另赋给ttxx=ScopeData.signals.values;%目的同上目的同上[xm,km]=max(xx);plot(tt,xx,'r','LineWidth',4);hold onplot(tt(km),xm,'b.','MarkerSize',36);hold offstrmax=char('最大值最大值',['t = ',num2str(tt(km))],['x = ',num2str(xm)]);text(6.5,xm,strmax),xlabel('t'),ylabel('x')2024/9/1261利用存放在MATLAB工作空间中的仿真数据所绘制的曲线2024/9/1262直接利用传递函数模块求解方程 2024/9/12632024/9/1264假设上式中的输入函数是单位脉冲函数，利用状态方程研究该系统的位移变化。

2024/9/1265利用使能子系统实现半波整流的仿真模型 2024/9/1266利用触发子系统获得零阶保持的采样信号本例演示：触发子系统工作原理；在MATLAB指令窗中运行SIMULINK模型 2024/9/1267假设有某过程的离散状态方程式中u(k)是输入该过程的采样周期为0.1秒控制器应用采样周期为0.25秒的比例控制器；显示系统的更新周期为0.5秒UD1,UD2单位延迟模块的采样周期设为0.1sZOH1的采样周期设为0.25s ，ZOH2的采样周期设为0.5s 2024/9/1268Sources 库库Band-Limited White Noise（（限带白噪声限带白噪声））Chirp Signal（（扫频信号扫频信号））Clock（（时钟时钟））Constant（（常量）常量）Digital Clock（（数字时钟数字时钟）） Discrete Pulse Generator（（离散脉冲生成器离散脉冲生成器））From Workspace（（从工作空间读取数据从工作空间读取数据））From File（（从文件读数据从文件读数据））Pulse Generator（（脉冲生成器脉冲生成器））Ramp（（倾斜）倾斜）Random Number（（随机数产生器随机数产生器））Repeating Sequence（（重复序列）重复序列）2024/9/1269Sources 库（续）库（续）Signal Generator（（信号发生信号发生器器））Sine Wave（（正弦波正弦波））Step（（阶跃阶跃））Uniform Random Number（（均匀分布随机数）均匀分布随机数）2024/9/1270Sinks Sinks 库库库库Display（（显示显示））Scope（（示波器示波器））Stop Simulation（（停止仿真）停止仿真）To File（（写入文件写入文件））To Workspace（（写到工作空间写到工作空间））XY Graph（（显示平面图形显示平面图形））2024/9/1271Discrete 库库Discrete Filter（（离散滤波器离散滤波器））Discrete State Space（（离散状态空间离散状态空间））Discrete-Time Integrator（（离散时间积分器离散时间积分器））Discrete Transfer Fcn（（离散传递函数离散传递函数））Discrete Zero-Pole（（数字零极点函数数字零极点函数））First-Order Hold（（一阶保持）一阶保持）Zero-Order Hold（（零阶保持零阶保持））Unit Delay（（单位延迟单位延迟））2024/9/1272Continuous 库库Derivative（（导数导数））Integrator（（积分器积分器））Memory（（记忆）记忆）State Space（（状态空间状态空间）） Transfer Fcn（（传递函数传递函数）） Transport Delay（（传递延迟传递延迟）） Variable Transport Delay（（可变传输延迟）可变传输延迟） Zero-Pole（（零零-极点极点）） 2024/9/1273Math 库库Abs（（绝对值绝对值））Algebraic ConstraintCombinatorial Logic（（组合组合逻辑逻辑））Complex to Magnitude-Angle Complex to Real-Image Dot Product（（点乘点乘））Gain（（增益增益）） Logical Operator（（逻辑逻辑运算运算）） Magnitude--Angle to Complex（）（）Math Function（（数学函数数学函数））Matrix Gain（（矩阵增益矩阵增益））MinMax（（最大最小值最大最小值））2024/9/1274Math 库（续）库（续）Product（（乘积乘积））Real--Image to ComplexRelational Operator（（关系关系运算运算））Rounding Function（（圆整函数圆整函数））Sign（（符号符号）） Slider Gain（（滑块增益滑块增益）） Sum（（和和）） Trigonometric Function（（三角函数三角函数）） 2024/9/1275Nonlinear 库库Backlash 模块模块Coulomb and Viscous Friction（（库仑和粘性摩擦库仑和粘性摩擦））Dead Zone（（死区）死区）Manual Switch（（手动开关手动开关）） Multiport Switch（（多路转换开关多路转换开关）） Quantizer（（量化）量化） Rate Limiter（（限速器限速器））Relay（（继电器继电器）） Saturation（（饱和）饱和）Switch（（选择开关选择开关）） 2024/9/1276Signals & Systems Signals & Systems 库库Bus SelectorBus Selector（（总线选择器总线选择器））Configurable SubsystemConfigurable Subsystem（（可配置子系统可配置子系统））Data Store MemoryData Store Memory（（数据存储器数据存储器））Data Store ReadData Store Read（（读数据存储读数据存储）） Data Store WriteData Store Write（（写数据存储写数据存储）） Data Type ConversionData Type Conversion（（数据类型转换数据类型转换）） DemuxDemux（（解混）解混） EnableEnable（（激活）激活） FromFrom（（导入）导入）GotoGoto（（传出传出））GotoGoto Tag Visibility Tag Visibility（（传出标记符的可见性传出标记符的可见性））GroundGround（（接地接地））Hit CrossingHit Crossing（（捕获穿越点捕获穿越点））2024/9/1277Signals & Systems 库（续）库（续）IC（（初始状态初始状态））Inport（（输入端口输入端口））Merge（（合并合并））Model Info （（模型信息）模型信息）Mux （（混合混合））Outport（（输出端口输出端口）） Probe （（探测器探测器））Selector（（选择器选择器）） Subsystem（（子系统）子系统）Terminator（（终结器终结器））Trigger（（触发器触发器））Width（（宽度宽度））Function-Call Generator（（函数调用发生器函数调用发生器））2024/9/1278Functions & Tables 库库Fcn（（函数表达式函数表达式））Look-Up Table（（查找表查找表））Look-Up Table （（2-D）（）（二维查找表二维查找表））MATLAB Fcn （（MATLAB 函数函数））S-Function（（ S 函数函数）） 2024/9/1279S-函数 •S-函数是系统函数(System Function)的简称，是指采用非图形化的方式(即计算机语言，区别于Simulink的系统模块)描述的一个功能块。

用户可以采用MATLAB代码，C，C++，FORTRAM或Ada等语言编写S-函数S-函数由一种特定的语法构成，用来描述并实现连续系统、离散系统以及复合系统等动态系统；S-函数能够接受来自Simulink求解器的相关信息，并对求解器发出的命令作出适当的响应，这种交互作用非常类似于Simulink系统模块与求解器的交互作用一个结构体系完整的S-函数包含了描述动态系统所需的全部能力，所有其他的使用情况都是这个结构体系的特例往往S-函数模块是整个Simulink动态系统的核心2024/9/1280S-函数•根据S-函数代码使用的编程语言，S-函数可以分成 M文件 S-函数（即用MATLAB语言编写的S-函数）、C语言S-函数、C++语言S-函数、Ada语言S-函数以及Fortran语言S-函数等通过S-函数创建的模块具有与Simulink模型库中的模块相同的特征，它可以与Simulink求解器进行交互，支持连续状态和离散状态模型 2024/9/1281S-函数•S-函数作为与其他语言相结合的接口，可以使用这个语言所提供的强大能力例如，Matlab语言编写的S-函数可以充分利用MATLAB所提供的丰富资源，方便地调用各种工具箱函数和图形函数；使用C语言编写的S-函数可以实现对操作系统的访问，如实现与其他进程的通信和同步等。

2024/9/1282S-函数•用用户户可可能能会会有有如如下下的的疑疑问问：：SimulinkSimulink已已经经提提供供了了大大量量的的内内置置的的系系统统模模块块，，并并且且允允许许用用户户自自定定义义模模块块，，那那么么为为何何还还要要使使用用S-S-函函数数呢呢？？诚诚然然，，对对于于大大多多数数动动态态系系统统仿仿真真分分析析语语言言，，使使用用SimulinkSimulink提提供供的的模模块块即即可可实实现现，，而而无无需需使使用用S-S-函函数数但但是是，，当当需需要要开开发发一一个个新新的的通通用用的的模模块块作作为为一一个个独独立立的的功功能能单单元元时时，，使使用用S-S-函函数数实实现现则则是是一一种种相相当当简简便便的的方方法法另另外外，，由由于于S-S-函函数数可可以以使使用用多多种种语语言言编编写写，，因因此此可可以以将将已已有有的的代代码码结结合合进进来来，，而而不不需需要要在在SimulinkSimulink中中重重新新实实现现算算法法，，从从而而在在某某种种程程度度上上实实现现了了代代码码移移植植此此外外，，在在S-S-函函数数中中使使用用文文本本方方式式输输入入公公式式、、方方程程，，非非常常适适合合复复杂杂动动态态系系统统的的数数学学描描述述，，并并且且在在仿仿真真过过程程中中可可以对仿真进行更精确的控制。

以对仿真进行更精确的控制2024/9/1283•simulinksimulink的仿真过程（以便理解的仿真过程（以便理解s s函数），函数），simulinksimulink的仿真有两个阶段：一个为初始的仿真有两个阶段：一个为初始化，这个阶段主要是设置一些参数，像系化，这个阶段主要是设置一些参数，像系统的输入输出个数、状态初值、采样时间统的输入输出个数、状态初值、采样时间等；第二个阶段就是运行阶段，这个阶段等；第二个阶段就是运行阶段，这个阶段里要进行计算输出、更新离散状态、计算里要进行计算输出、更新离散状态、计算连续状态等等，这个阶段需要反复运行，连续状态等等，这个阶段需要反复运行，直至结束直至结束. .2024/9/1284•在在matlabmatlab的的workspaceworkspace里打里打edit edit sfuntmplsfuntmpl( (这是这是matlabmatlab自己提供的自己提供的s s函数模板函数模板) )，我们看它来具体分析，我们看它来具体分析s s函数的函数的结构它的第一行是这样的：它的第一行是这样的：•function [sys,x0,str,ts]=function [sys,x0,str,ts]=sfuntmpl(t,x,u,flagsfuntmpl(t,x,u,flag) )输入与输出变量的含义：输入与输出变量的含义：t t是采样时间，是采样时间，x x是状态变量，是状态变量，u u是输入（是做成是输入（是做成simulinksimulink模块的输入）模块的输入）,flag,flag是仿真过是仿真过程中的状态标志（以它来判断当前是初始化还是运行等）程中的状态标志（以它来判断当前是初始化还是运行等）；；syssys输出根据输出根据flagflag的不同而不同（下面将结合的不同而不同（下面将结合flagflag来来讲讲syssys的含义），的含义），x0x0是状态变量的初始值，是状态变量的初始值，strstr是保留参是保留参数（数（mathworksmathworks公司还没想好该怎么用它，一般在初始公司还没想好该怎么用它，一般在初始化中将它置空就可以了化中将它置空就可以了, ,strstr=[])=[])，，tsts是一个是一个1 1××2 2的向量，的向量，ts(1)ts(1)是采样周期，是采样周期，ts(2)ts(2)是偏移量。

是偏移量2024/9/1285•switch flagswitch flag, , % %判断判断flagflag，看当前处于哪个状，看当前处于哪个状态态case 0case 0, , [sys,x0,str,ts]=[sys,x0,str,ts]=mdlInitializeSizesmdlInitializeSizes; ;flag=0flag=0表示处于初始化状态，此时用函数表示处于初始化状态，此时用函数mdlInitializeSizesmdlInitializeSizes进行初始化，此函数在进行初始化，此函数在 sfuntmpl.msfuntmpl.m的的149149行行我们找到他，在初始化状态下，我们找到他，在初始化状态下，syssys是一个结构体，用它来是一个结构体，用它来设置模块的一些参数，各个参数详细说明如下设置模块的一些参数，各个参数详细说明如下 size = size = simsizessimsizes;%;%用于设置模块参数的结构体用用于设置模块参数的结构体用simsizessimsizes来生成来生成 sizes.NumContStatessizes.NumContStates = 0;% = 0;%模块连续状态变量的模块连续状态变量的个数个数 sizes.NumDiscStatessizes.NumDiscStates = 0;% = 0;%模块离散状态变量的模块离散状态变量的个数个数 sizes.NumOutputssizes.NumOutputs = 0;% = 0;%模块输出变量的模块输出变量的个数个数 sizes.NumInputssizes.NumInputs = 0;% = 0;%模块输入变量模块输入变量的个数的个数 sizes.DirFeedthroughsizes.DirFeedthrough = 1;% = 1;%模块是否存在直接贯模块是否存在直接贯通通 sizes.NumSampleTimessizes.NumSampleTimes = 1;% = 1;%模块的采样时间个数，至少是模块的采样时间个数，至少是一个一个 sys = sys = simsizes(sizessimsizes(sizes);); % %设置完后赋给设置完后赋给syssys输出输出2024/9/1286•case 1, sys=mdlDerivatives(t,x,u);flag=1表示此时要计算连续状态的微分， •case 2, sys=mdlUpdate(t,x,u);flag=2表示此时要计算下一个离散状态，表示此时要计算下一个离散状态， 2024/9/1287•case 3,case 3, sys= sys=mdlOutputs(t,x,umdlOutputs(t,x,u););flag=3flag=3表示此时要计算输出，表示此时要计算输出，•case 4,case 4, sys=sys=mdlGetTimeOfNextVarHit(t,x,umdlGetTimeOfNextVarHit(t,x,u););flag=4flag=4表示此时要计算下一次采样的时间，只在离散采表示此时要计算下一次采样的时间，只在离散采样系统中有用样系统中有用( (即上文的即上文的mdlInitmdlInit ializeSizesializeSizes中提到的中提到的tsts设置设置ts(1)ts(1)不为不为0)0)连续系统中只需在连续系统中只需在mdlGetTimeOfNextVarHitmdlGetTimeOfNextVarHit函数中写上函数中写上sys=[];sys=[];这个函数主要用于变步长的设置，具体实现大这个函数主要用于变步长的设置，具体实现大家可以用家可以用edit edit vsfuncvsfunc看看vsfunc.mvsfunc.m这个例子这个例子2024/9/1288•case 9, sys=mdlTerminate(t,x,u);flag=9表示此时系统要结束，一般来说写上在mdlTerminate函数中写上sys=[]就可，如果你在结束时还要设置什么，就在此函数中写2024/9/1289•s s函数还可以带用户参数，下面给个例子，和函数还可以带用户参数，下面给个例子，和simulinksimulink下的下的gaingain模块功能一样模块功能一样function [sys,x0,str,ts] = function [sys,x0,str,ts] = sfungain(t,x,u,flag,gainsfungain(t,x,u,flag,gain) )switch flag,switch flag,case 0,case 0, sizes = sizes = simsizessimsizes; ; sizes.NumContStatessizes.NumContStates = 0; = 0; sizes.NumDiscStatessizes.NumDiscStates = 0; = 0; sizes.NumOutputssizes.NumOutputs = 1; = 1; sizes.NumInputssizes.NumInputs = 1; = 1; sizes.DirFeedthroughsizes.DirFeedthrough = 1; = 1; sizes.NumSampleTimessizes.NumSampleTimes = 1; = 1; sys = sys = simsizes(sizessimsizes(sizes);); x0=[]; x0=[]; strstr=[];=[]; tsts=[0,0];=[0,0];case 3,case 3, sys=gain*u; sys=gain*u;case {1,2,4,9},case {1,2,4,9}, sys = []; sys = [];end end 2024/9/1290第六章第六章 Matlab在控制系统中的应用在控制系统中的应用1、控制系统的基本理论2、控制系统工具箱函数3、控制系统分析与设计2024/9/1291控制系统模型传递函数模型零极点增益模型状态空间模型2024/9/1292模型的转换1[num,den]=ss2tf(a,b,c,d) [num,den]=ss2tf(a,b,c,d) [a,b,c,d]=tf2ss(num,den)[a,b,c,d]=tf2ss(num,den)[z,p,k]=tf2zp(num,den)[z,p,k]=tf2zp(num,den)[num,den]=zp2tf(z,p,k)[num,den]=zp2tf(z,p,k)[a,b,c,d]=zp2ss(z,p,k)[a,b,c,d]=zp2ss(z,p,k)[z,p,k]=ss2zp(a,b,c,d)[z,p,k]=ss2zp(a,b,c,d)2024/9/1293例求传函和状态空间模型 k=6;z=[-3];p=[-1,-2,-5]; pzmap(p,z)[num,den]=zp2tf(z,p,k)tf(num,den)[A,B,C,D]=zp2ss(z,p,k)[zz,pp,kk]=ss2zp(A,B,C,D)2024/9/1294s1=tf([3 4 5],[1 3 5 7 9],'InputName','U','OutputName','Y')s2=tf([3 4 5],[1 3 5 7 9],0.1,'InputName','U','OutputName','Y')Transfer function from input "U" to output "Y":Transfer function from input "U" to output "Y": 3 s^2 + 4 s + 5 3 s^2 + 4 s + 5----------------------------------------------------------s^4 + 3 s^3 + 5 s^2 + 7 s + 9s^4 + 3 s^3 + 5 s^2 + 7 s + 9Transfer function from input "U" to output "Y":Transfer function from input "U" to output "Y": 3 z^2 + 4 z + 5 3 z^2 + 4 z + 5----------------------------------------------------------z^4 + 3 z^3 + 5 z^2 + 7 z + 9z^4 + 3 z^3 + 5 z^2 + 7 z + 9Sampling time: 0.1Sampling time: 0.12024/9/1295模型的转换2连续时间系统连续时间系统===== =离散时间系统离散时间系统连续到离散[ad,bd]=c2dc2d(a,b,Ts)[ad,bd,cd,dd]=c2dtc2dt(a,b,c,d,Ts,lambda) 带输入延时[ad,bd,cd,dd]=c2dmc2dm(a,b,c,d,Ts,’method’)[num,den]=c2dmc2dm(num,den,Ts,’method’)方法：‘zoh’默认零阶保持器 ‘foh’默认一阶保持器(无逆变) ‘tustin’利用双线性逼近导数‘prewrap’利用频率预变的双线性来逼近‘matched’利用匹配零-极点方法将SISO系统变换2024/9/1296模型的转换3连续时间系统连续时间系统===== =离散时间系统离散时间系统离散到连续[ac,bc]=d2cd2c(ad,bd, Ts)[ac,bc,cc,dc]=d2cd2c(ad,bd,cd,dd,Ts)[ac,bc,cc,dc]=d2cmd2cm(ad,bd,cd,dd,Ts,’method’)[num,den]=d2cmd2cm(num,den,Ts,’method’)2024/9/1297[a,b,c,d]=ord2(1,0.2)step(a,b,c,d);hold on[aa,bb,cc,dd]=c2dm(a,b,c,d,0.5,'tustin')dstep(aa,bb,cc,dd)c2dm(a,b,c,d,0.5,'tustin')2024/9/1298用零阶保持器和双线性变换求离散传递函数format compactf=[-4 1]; g=[1 2 10]; ts=0.2sc=tf(f,g)disp('零阶保持器')sd1=c2d(sc,ts)disp('双线性变换')sd2=c2d(sc,ts,'t')零阶保持器 Transfer function: -0.5994 z + 0.6313----------------------z^2 - 1.351 z + 0.6703Sampling time: 0.2双线性变换 Transfer function:-0.3 z^2 + 0.01538 z + 0.3154----------------------------- z^2 - 1.385 z + 0.6923Sampling time: 0.22024/9/1299% 连续和离散系统的多种输出响应曲线clear,clf[a,b,c,d]=rmodel(4);s1=ss(a,b,c,d);Ts=0.2;sd1=c2d(s1,Ts,'t')t=0:Ts:25;u=sin(0.5*t);for i=1:2 if i==1 s=s1; else s=sd1; end figure(i) subplot(2,2,1) impulse(s,5);grid subplot(2,2,2) lsim(s,u,t); subplot(2,2,3) x0=[1,-1,0,2]; initial(s,x0,5),grid subplot(2,2,4) pzmap(s),gridend2024/9/1300典型系统的生成s=rss(n) 或s=rss(n，p出，m入)随机生成n阶稳定的连续状态空间模型[num,den]=rmodel(n)随机生成n阶稳定的连续线性模型系数s=drss(n)随机生成n阶稳定的离散状态空间模型[num,den]=drmodel(n)随机生成n阶稳定的离散线性模型系数[num,den]=ord2(wn,z)生成固有频率为wn，阻尼系数为z的二阶系统系数2024/9/1301系统的建模1系统的并联parallel[a,b,c,d]=parallel(a1,b1,c1,d1,a2,b2,c2,d2);[num,den]=parallel(num1,den1,num2,den2)[a,b,c,d]=parallel(a1,b1,c1,d1,a2,b2,c2,d2,in1,in2,out1,out2)将输入in1和in2合并为一个输入将输出out1和out2合并为一个输出2024/9/1302系统的建模2系统的串联series[a,b,c,d]=series(a1,b1,c1,d1,a2,b2,c2,d2);[num,den]=series(num1,den1,num2,den2)[a,b,c,d]=series(a1,b1,c1,d1,a2,b2,c2,d2,out1,in2)将系统1的输出与系统2的输入进行串联2024/9/1303系统的建模3系统的反馈feedback[a,b,c,d]=feedback(a1,b1,c1,d1,a2,b2,c2,d2)[a,b,c,d]=feedback(a1,b1,c1,d1,a2,b2,c2,d2,signsign);[num,den]=feedback(num1,den1,num2,den2)[num,den]=feedback(num1,den1,num2,den2,sign)sign缺省时，默认为负，即负反馈[a,b,c,d]=feedback(a1,b1,c1,d1,a2,b2,c2,d2,inp1,out1inp1,out1);将系统1的指定输出通过系统2反馈到系统1的输入。

2024/9/1304系统的建模4系统的闭环cloop[ac,bc,cc,dc]=cloop(a,b,c,d,sign)[numc,denc]=cloop(num,den,sign)[ac,bc,cc,dc]=cloop(a,b,c,d,out,in)将系统的输出反馈到输入，形成单位反馈2024/9/1305a1=[0,1;-1,-2]; b1=[0;1];c1=[1,3]; d1=1;a2=[0,1;-1,-3]; b2=[0;1];c2=[1,4]; d2=[0];[a,b,c,d]=series(a1,b1,c1,d1,a2,b2,c2,d2)[a,b,c,d]=parallel(a1,b1,c1,d1,a2,b2,c2,d2)[a,b,c,d]=feedback(a1,b1,c1,d1,a2,b2,c2,d2)[a,b,c,d]=feedback(a1,b1,c1,d1,a2,b2,c2,d2,1)2024/9/1306sys1=tf([2 1],[3 2 1]);sys2=tf(1,[3 2]);H=feedback(sys1,sys2);p=pole(H);ni=find(real(p)>0);n=length(ni);if n>0 disp(‘系统不稳定’)else disp(‘系统稳定’)end 2024/9/1307模型特性1可控性矩阵可控性矩阵 [co]=ctrb(a,b) cco=length(a)-rank(co) 如果cco==0表示可控可观性矩阵可观性矩阵 [ob]=obsv(a,c) oob=length(a)-rank(ob) 如果oob==0表示可观2024/9/1308模型特性2[aa,bb,cc,T,k]=ctrbf(a,b,c) 可控性分解[aa,bb,cc,T,k] obsvf(a,b,c) 可观性分解sysT=ss2ss(sys,T) 相似变换[csys,T]=canon(sys,’type’) 规范分解type=model 约当矩阵形式type=companion 伴随矩阵形式2024/9/1309系统可控和可观判断a=[-3,1;1,-3]; b=[1,1;1,1]c=[1,1;1,-1] ; d=0cam=ctrb(a,b); cco=length(a)-rank(cam)if cco==0 disp('系统可控')else disp('系统不可控')endoam=obsv(a,c); oob=length(a)-rank(oam)if oob==0 disp('系统可观')else disp('系统不可观')end2024/9/1310系统可控和可观分解% % 化为可控阶梯形化为可控阶梯形clear,A=[-2,2,-1;0,-2,0;1,-4,3];B=[0;0;1];C=[1,-1,1];D=0; %系数矩阵赋值s1=ss(A,B,C,D);% 构成LTI模型[Abar,Bbar,Cbar,T,k]=ctrbf(s1.a,s1.b,s1.c) % 化为可控阶梯型rA=rank(A)% 状态方程系数矩阵的秩rc=sum(k)% 判断可控矩阵的秩[Abaro,Bbaro,Cbaro,To,ko]=obsvf(s1.a,s1.b,s1.c);% 化为可观阶梯型ro=sum(ko)% 判断可观矩阵的秩2024/9/1311时域响应1求连续系统的单位阶跃响应求连续系统的单位阶跃响应[y,x,t]=step(a,b,c,d)[y,x,t]=step(a,b,c,d,iu)[y,x,t]=step(a,b,c,d,iu,t)[y,x,t]=step(num,den)[y,x,t]=step(num,den,t)不带输出变量引用时绘制出响应曲线t为用户指定时间绘制曲线，iu为指定输入的输出2024/9/1312时域响应2求离散系统的单位阶跃响应求离散系统的单位阶跃响应[y,x]=dstep(a,b,c,d)[y,x]=dstep(a,b,c,d,iu)[y,x]=dstep(a,b,c,d,iu,n)[y,x]=dstep(num,den)[y,x]=dstep(num,den,n)不带输出变量引用时绘制出响应曲线n为用户指定时间绘制曲线，iu为指定输入的输出2024/9/1313时域响应3求连续系统的单位冲激响应求连续系统的单位冲激响应[y,x,t]=impulse(a,b,c,d)[y,x,t]=impulse(a,b,c,d,iu)[y,x,t]=impulse(a,b,c,d,iu,t)[y,x,t]=impulse(num,den)[y,x,t]=impulse(num,den,t)不带输出变量引用时绘制出响应曲线t为用户指定时间绘制曲线，iu为指定输入的输出2024/9/1314时域响应4求离散系统的单位冲激响应求离散系统的单位冲激响应[y,x]=dimpulse(a,b,c,d)[y,x]=dimpulse(a,b,c,d,iu)[y,x]=dimpulse(a,b,c,d,iu,n)[y,x]=dimpulse(num,den)[y,x]=dimpulse(num,den,n)不带输出变量引用时绘制出响应曲线n为用户指定时间绘制曲线，iu为指定输入的输出2024/9/13151、求时单位阶跃响应wn=6;kosi=[0.1:0.1:1,2.0]figure(1)hold onfor ks=kosi num=wn.^2; den=[1,2*ks*wn,wn.^2] step(num,den) %impulse(num,den)endtitle('Step Response'); hold off2024/9/13162、求时单位阶跃响应kosi =0.7;wn =[2:2:12]figure(2)hold onfor wn=wn num=wn.^2; den=[1,2*kosi*wn,wn.^2] step(num,den) %impulse(num,den)endtitle('Step Response'); hold off2024/9/1317脉冲响应和阶跃响应：脉冲响应和阶跃响应：a = [-0.5572 -0.7814;0.7814 0];b = [1 -1;0 2];c = [1.9691 6.4493];sys = ss(a,b,c,0);impulse(sys)step(sys)2024/9/1318时域响应5求连续系统的零输入响应求连续系统的零输入响应[y,x,t]=initial(a,b,c,d,x0)[y,x,t]=initial(a,b,c,d,x0,t)不带输出变量引用时绘制出响应曲线t为用户指定时间绘制曲线x0为指定初始状态2024/9/1319时域响应6求离散系统的零输入响应求离散系统的零输入响应[y,x]=dinitial(a,b,c,d,x0)[y,x]=dinitial(a,b,c,d,x0,n)不带输出变量引用时绘制出响应曲线t为用户指定时间绘制曲线x0为指定初始状态2024/9/1320[a,b,c,d]=ord2(0.7,0.35)[y1]=step(a,b,c,d)x0=[1;1][y2]=initial(a,b,c,d,x0)figure(1)plot(y1);hold onplot(y2,'k')plot(y1+y2,'r')2024/9/1321时域响应7对任意输入的连续系统进行仿真对任意输入的连续系统进行仿真[ [y,x]=y,x]=lsim(a,b,c,d,u,tlsim(a,b,c,d,u,t) )[ [y,x]=lsim(a,b,c,d,u,t,x0)y,x]=lsim(a,b,c,d,u,t,x0)[ [y,x]=y,x]=lsim(num,den,u,tlsim(num,den,u,t) )num=[2,5,1]; den=[1 2 3];num=[2,5,1]; den=[1 2 3];t=[0:0.1:10]; period=4; t=[0:0.1:10]; period=4; % %输入为方波输入为方波T=4T=4u=(u=(rem(t,periodrem(t,period)>=period./2);)>=period./2);lsim(num,den,u,tlsim(num,den,u,t) )2024/9/1322时域响应8对任意输入的离散系统进行仿真对任意输入的离散系统进行仿真[ [y,x]=y,x]=dlsim(a,b,c,d,udlsim(a,b,c,d,u) )[ [y,x]=dlsim(a,b,c,d,u,x0)y,x]=dlsim(a,b,c,d,u,x0)[ [y,x]=y,x]=dlsim(num,den,udlsim(num,den,u) )2024/9/1323求系统对100点随机噪声的响应num=[2 ,-3.4,1.5]den=[1 -1.6 0.8];u=rand(100,1);dlsim(num,den,u)title('Noise Response')2024/9/1324频率响应1求连续系统的求连续系统的BodeBode频率响应频率响应[mag,phase,w]=bode(a,b,c,d)[mag,phase,w]=bode(a,b,c,d,iu)[mag,phase,w]=bode(a,b,c,d,iu,w)[mag,phase,w]=bode(num,den)[mag,phase,w]=bode(num,den,w)函数中w为指定绘图频段,无返值直接绘制图2024/9/1325[a,b,c,d]=ord2(1,0.2);bode(a,b,c,d);grid on[mag,phase,w]=bode(a,b,c,d)title('Bode Plot')figure(2); subplot(2,1,1);semilogx(w,20*log10(mag)); gridsubplot(2,1,2);semilogx(w,phase);grid2024/9/1326频率响应2求离散系统的求离散系统的BodeBode频率响应频率响应[mag,phase,w]=dbode(a,b,c,d,Ts)[mag,phase,w]=dbode(a,b,c,d,Ts,iu)[mag,phase,w]=dbode(a,b,c,d,Ts,iu,w)[mag,phase,w]=dbode(num,den,Ts)[mag,phase,w]=dbode(num,den,Ts,w)函数中w为指定绘图频段,无返值直接绘制图2024/9/1327num=[2 ,-3.4,1.5]den=[1 -1.6 0.8];dbode(num,den,0.3);gridtitle('Discrte Bode Plot')[mag,phase,w]=dbode(num,den,0.3)2024/9/1328频率响应3求连续系统的求连续系统的NyquistNyquist频率频率曲线（极坐标图）曲线（极坐标图）[re,im,w]=nyquist(a,b,c,d)[re,im,w]=nyquist(a,b,c,d,iu)[re,im,w]=nyquist(a,b,c,d,iu,w)[re,im,w]=nyquist(num,den)[re,im,w]=nyquist(num,den,w)函数中w为指定绘图频段,无返值直接绘制2024/9/1329num=[2 5 1];den=[1 2 3]nyquist(num,den)title('Nyquist Plot')2024/9/1330频率响应4求求离散离散系统的系统的NyquistNyquist频率频率曲线（极坐标图）曲线（极坐标图）[re,im,w]=dnyquist(a,b,c,d,Ts)[re,im,w]=dnyquist(a,b,c,d,Ts,iu)[re,im,w]=dnyquist(a,b,c,d,Ts,iu,w)[re,im,w]=dnyquist(num,den,Ts)[re,im,w]=dnyquist(num,den,Ts,w)函数中w为指定绘图频段,无返值直接绘制2024/9/1331num=[2 ,-3.4,1.5]den=[1 -1.6 0.8];dnyquist(num,den,0.3);title('Discrte Nyquist Plot')2024/9/1332频率响应5求连续系统的求连续系统的增益和相角裕度增益和相角裕度[ [gm,pm,wc,wggm,pm,wc,wg]=]=margin(mag,phase,wmargin(mag,phase,w) )[ [gm,pm,wc,wggm,pm,wc,wg]=margin(num,den)]=margin(num,den)[ [gm,pm,wc,wggm,pm,wc,wg]=margin(a,b,c,d)]=margin(a,b,c,d)函数中w为指定绘图频段,无返值直接绘制2024/9/1333num=[5,5]den=[1 2 3 0];bode(num,den);grid[mag,phase,w]=bode(num,den)[gm,pm,wcp,wcg]=margin(mag,phase,w)[gm,pm,wcp,wcg]=1.2823e+004 38.5930 253.1864 2.36602024/9/1334% 四阶连续和离散系统开闭环bode图clear,Ts=0.1s=zpk(-6,[0,-1,-10,-10],200) % 生成四阶连续系统ssd=c2d(s,Ts) % 变换为四阶离散系统sdsb=feedback(s,1)% 把连续系统闭合，生成闭环系统sbsbd=feedback(sd,1)% 把离散系统闭合，生成闭环系统sbdfigure(1),gridbode(s,'--',sb,'.-')% 在图1中同时绘出连续系统的开闭环频率特性figure(2),gridbode(sd,'--',sbd,'.-')% 在图2中同时绘出离散系统的开闭环频率特性damp(sb)% 闭环sb的根的固有频率和阻尼系数damp(sbd)% 判断闭环sbd的根的固有频率和阻尼系数[Gm,Pm,wcg,wcp]=margin(s)% 判断闭环sb的稳定裕度[Gmd,Pmd,wcgd,wcpd]=margin(sd) % 判断闭环sb的稳定裕度2024/9/1335根轨迹1绘制系统的零极点图绘制系统的零极点图[p,z]=pzmap(a,b,c,d)[p,z]=pzmap(num,den) pzmap(p,z) % %绘制给定零极点图绘制给定零极点图2024/9/1336num=[2 5 1]den=[1 2 3]pzmap(num,den)title('Pole-Zero Map')[p,z]=pzmap(num,den)p = z = -1.0000 + 1.4142i -2.2808 -1.0000 - 1.4142i -0.21922024/9/1337根轨迹2求系统的根轨迹求系统的根轨迹[r,k]=rlocus(num,den)[r,k]=rlocus(num,den,k)[r,k]=rlocus(a,b,c,d)[r,k]=rlocus(a,b,c,d,k)无返回变量直接绘制曲线，函数中的k为指定k值的根轨迹2024/9/1338根轨迹3sgrid %绘制阻尼系数和自然频率栅格sgrid(z,wn) %绘制指定指定阻尼系数和自然频率栅格sgrid(‘new’) %重画阻尼系数和自然频率栅格zgrid %绘制离散系统阻尼系数和自然频率栅格zgrid(z,wn) %绘制指定指定阻尼系数和自然频率栅格zgrid(‘new’) %重画阻尼系数和自然频率栅格2024/9/1339根轨迹4计算给定一组根的根轨迹增益计算给定一组根的根轨迹增益[k,ploes]=rlocfind(a,b,c,d)[k,ploes]=rlocfind(a,b,c,d,p)[k,ploes]=rlocfind(num,den)[k,ploes]=rlocfind(num,den,p)%可指定要得到增益的根矢量p2024/9/1340num=[2 5 1] % num=[1] num=[1 2]den=[1 2 3] %den=[1 16 36 80 0] den=conv([1 4 3],[1 4 3])rlocus(num,den)title('Root Locus')[r,k]=rlocus(num,den)[k,poles]=rlocfind(num,den)sgrid2024/9/1341num=[1]den=[1 5 8 6 0]rlocus(num,den)title('Root locus')s=tf(num,den)sd=c2d(s,0.5,'t')figure(2)rlocus(sd)zgrid 1------------------------------s^4 + 5 s^3 + 8 s^2 + 6 sTransfer function: 1--------------------- =s(s+3)(s^2+2s+2) 2024/9/1342利用李亚普诺夫方程确定线性时不变系统的稳定性a=[-1 -2;1 -4]q=[q 0;0 1];p=lyap(a,q)detp=det(p)2024/9/1343极点配置和观测器设计k=place(a,b,p)k=acker(a,b,p)采用全反馈u=-kx；p为希望配置的闭环极点g=place(a',c',p)g=acker(a',c',p)采用全反馈u=-Gx；p为希望配置的观测器闭环极点2024/9/1344设计u=-kx 使闭环系统的极点为a=[0 1 0;0 0 1;-6 –11 –6];b=[0;0;10];p=[-2+2*sqrt(3)*i -2-2*sqrt(3)*i -10]k=place(a,b,p)2024/9/1345控制系统的设计方法控制系统的设计方法这里，所谓的设计设计控制系统，就是在已知被控对象的数学模型（一般为传递函数）的基础上，引入校正环校正环节节，使得系统的闭环暂态和稳态性能暂态和稳态性能满足要求。

工业上，一般称校正环节为控制器一般，设计系统要包括：控制系统结构设计，控制器控制算法设计、控制器参数设计等对于用传递函数描述的控制系统，常用的经典设计方法是根轨迹与频域法2024/9/1346例被控对象传函为，要求：（1）速度误差系数为10 （2）相角裕量为45度解：设则绘出KcG(jw)的伯德图可得：稳定裕量大约为32度，须补13度画伯德图：num=2000;den=[1,30,200,0];bode(num,den)grid on计算角度裕量：[m,p,w]=bode(num,den)找m=0对应p为-147度2024/9/1347考虑校正装置影响剪切频率的位置，再加上5度的裕量则fi=18度则校正网络在ωm处增益为10lgα，原伯德图上增益为-10lgα=-2.77dB可由数据中找到-2.77dB所对应的角频率大约为9.2rad/s，则T=0.08s因此可得校正装置为校正后，系统相角裕度为41度，接近设计要求。

2024/9/1348校正前后系统的动态仿真结果如下图所示可见：校正改善系统动态性能，系统响应速度加快但不利于稳态性能，降低抗干扰能力是高通滤波器，而噪声为高频）•其它常用校正方法有：•滞后校正：改善稳态性能，但引入相位滞后，对动态性能不利•超前-滞后校正：•反馈校正2024/9/1349控制系统的优化设计控制系统的优化设计所谓优化设计，是指在所有可能的设计方案中寻找具有最优目标的设计方法包含两方面内容：一是控制系统参数最优化问题，即在系统构成确定系统构成确定情况下选择适当的参数，使系统的某种性能某种性能达到最佳如PID控制中的参数整定另一个是系统控制器结构的最优化问题，即在被控对象被控对象确定确定情况下选择适当的控制结构（或控制规律），以使系统的某种性能达到最佳2024/9/1350控制系统的优化设计控制系统的优化设计1 1、几个概念、几个概念（1）设计变量某些有待选择的量值，如系统参数其初值可任意设定，不影响优化结果，只影响优化设计效率，如计算时间2）约束条件优化设计的结果不唯一，有些结果可能超出了某些设计的限制条件，应舍弃，继续寻求满足条件的优化结果这些限制条被称为约束条件。

3）目标函数所谓“最优”都是相对的，是一定条件下的最优这个条件，就是目标函数目标函数的选取是整个优化设计过程中的最重要的决策之一4）综合目标函数用来比较满足不同约束条件下各个最优方案之间的好坏如目标是快速性，限制条件有稳定、超调超调的地位比稳定要小，即权重小2024/9/13512 2、优化设计原理、优化设计原理----单纯形法单纯形法•优化设计就是要寻求一组最优的设计变量使目标函数最小或最大，即数学中求（多元）函数的极值问题，工程上称为“寻优问题”但是，由于工程问题往往难以列写出函数关系，由导数求极值法不适用•单纯形：单纯形：变量空间内最简单的规则形体如二维平面内的正三角形，三维空间的正四面体•寻优原理寻优原理：取一个正三角形，其三个顶点相当于三个设计方案，可以计算相应的目标函数值称最大的为“坏点”，最小的为“好点”寻优过程为：（1）寻优规则抛弃坏点，求其对称点，构成新正三角形坏点重复，抛弃次坏点2）终点判别三角形围绕同一点转圈时，即好点重复时，若满足T.=1.65N，为终点其中N为维数，T为翻转次数3）精度调整减小三角形边长，可提高精度单纯形顶点坐标计算：单纯形顶点坐标计算：新坐标点=2*留下n点坐标之和/n-抛弃点的坐标2024/9/13523 3 3 3、目标函数的选取、目标函数的选取、目标函数的选取、目标函数的选取对不同的优化设计问题，目标函数的选取方式是不一样的。

对不同的优化设计问题，目标函数的选取方式是不一样的针对如图控制系统，常有以下几种：针对如图控制系统，常有以下几种：（（1 1））IAE (Integral absolute error)IAE (Integral absolute error)（（2 2））ISE (Integral square error)ISE (Integral square error)（（3 3））ITAEITAE（（4 4））ITSEITSE（（5 5））ISTAEISTAE（（6 6））ISTSEISTSEPIDG(s)目标函数2024/9/13534、例：已知对象传函为，采用PI调节器，希望对单位阶跃信号具有快速响应特性，试设计调节器参数Kp，Ki值解：由于ITSE准则可以反映系统的快速性（ITAE也可以），选取ITSE准则下的目标函为输出建立simulink下仿真模型，命名为e59.mdl用命令：[t,x,y]=sim(‘e59.mdl’, 仿真时间,初值); 2024/9/1354•优化目标函数：优化目标函数：optm.moptm.m•function function ssss= =optm(xoptm(x) )•global global kpkp; global ; global kiki; ; global i;global i;•kpkp=x(1);ki=x(2);=x(1);ki=x(2);•i=i+1;i=i+1;•[ [tt,xx,yytt,xx,yy]=sim(‘e59.mdl’,]=sim(‘e59.mdl’,40,[]);40,[]);•yylongyylong= =length(yylength(yy););•ssss= =yy(yylongyy(yylong););•主程序：主程序：•global global kpkp; ;• global global kiki; ; •global i;global i;•i=1;i=1;•result=fmins(‘optm’,[2 result=fmins(‘optm’,[2 1]);1]);2024/9/13552024/9/1356MajpjMVcyzj21HLfrvy96dv02lPPfYgxUS7IYmZkyEmZ0kGeYZS3bpLCkYH1lt4EK7CxmUX3ijoYSOer7ZuaVWYgz4EpZrUirVpMzzvNtf1XZw5oswSXOtFaejnOcmfE1lZgnN1RSXg8wLCG8CVQ3XPJMvodPFWcpiYJgZazNSEPNIaklYSu7qSd1UpaxmZDlpN9zW7kljfsLCLi26Yv109ffbnDH8LbUN1G6ACURQ39eG12KHL9tXsZ1jzgoCK8g1kuNOh5eFvcmVT5ZYVQt9zk3rp3qLnf02FovEXxVRxjCcFRNppiJljNiOuk6fONnyX7fyGg7sXZ49BmCN5oy9VesHpKzdjTKwjrkCEQCFDehVmGax3lrOEbw63VscA3YSijtUKoCyiLzAlVRp7l4QgPNHxvJFFDyjUVN3oHlMah0XBd4uTbkfPIhHtw0evPmYOrdhEDoPwvYhzlGplU1AU9mpyiCXH8gpPCBRYjq77VcnbXumNE1yGfyTsbSj89J63kRTKDkKUg3mdS5sJ4X5cQ8dK7oW9IkScssECQdz2O9UTlpRjAFPChjhLdzopQzwxQf8ozdzOhogwAooXpUF83BX4C3jRgjDJiiXEUDMaNz4vQ4n164vspddHvOIVuBBdMA4xp1YhiHk0vOJ8TL1BxogzVlMpmod6ianYGmksQq6NWCEd56hZF4wfaNyZcrGfNxnPiG6ZAxSkfmhJAKtNMajpjMVcyzj21HLfrvy96dv02lPPfYgxUS7IYmZkyEmZ0kGeYZS3bpLCkYH1lt4EK7CxmUX3ijoYSOer7ZuaVWYgz4EpZrUirVpMzzvNtf1XZw5oswSXOtFaejnOcmfE1lZgnN1RSXg8wLCG8CVQ3XPJMvodPFWcpiYJgZazNSEPNIaklYSu7qSd1UpaxmZDlpN9zW7kljfsLCLi26Yv109ffbnDH8LbUN1G6ACURQ39eG12KHL9tXsZ1jzgoCK8g1kuNOh5eFvcmVT5ZYVQt9zk3rp3qLnf02FovEXxVRxjCcFRNppiJljNiOuk6fONnyX7fyGg7sXZ49BmCN5oy9VesHpKzdjTKwjrkCEQCFDehVmGax3lrOEbw63VscA3YSijtUKoCyiLzAlVRp7l4QgPNHxvJFFDyjUVN3oHlMah0XBd4uTbkfPIhHtw0evPmYOrdhEDoPwvYhzlGplU1AU9mpyiCXH8gpPCBRYjq77VcnbXumNE1yGfyTsbSj89J63kRTKDkKUg3mdS5sJ4X5cQ8dK7oW9IkScssECQdz2O9UTlpRjAFPChjhLdzopQzwxQf8ozdzOhogwAooXpUF83BX4C3jRgjDJiiXEUDMaNz4vQ4n164vspddHvOIVuBBdMA4xp1YhiHk0vOJ8TL1BxogzVlMpmod6ianYGmksQq6NWCEd56hZF4wfaNyZcrGfNxnPiG6ZAxSkfmhJAKtNmCqbRmppeXp8inz4eq3HkWCMSORyMMX522xpHG6basNr6KQfbZsFbHjzyNlJrruLolKFcC84dqfijBO5Dy2NaBcNEBPgQrT12PgpcKx2or2YChN5DPjs80zzdtdAdTKuW4uVv9bbZu3K2SZ2aEhTlIC1UqrIWibkzwHh6p8gLv26zr01mJybfOzFc4T7kQH1IpPwOzMDnAKPLsLrznXGjFNIA9bSWWms6ibKZwQIKrMzalwbFrQJvOP1rPH8rx2KkyYqrtQk5VRwM1HSXmCqbRmppeXp8inz4eq3HkWCMSORyMMX522xpHG6basNr6KQfbZsFbHjzyNlJrruLolKFcC84dqfijBO5Dy2NaBcNEBPgQrT12PgpcKx2or2YChN5DPjs80zzdtdAdTKuW4uVv9bbZu3K2SZ2aEhTlIC1UqrIWibkzwHh6p8gLv26zr01mJybfOzFc4T7kQH1IpPwOzMDnAKPLsLrznXGjFNIA9bSWWms6ibKZwQIKrMzalwbFrQJvOP1rPH8rx2KkyYqrtQk5VRwM1HSX2024/9/1357。

点击阅读更多内容