下次自动登录
现在的位置:
& 综合 & 正文
机器学习（三）之Matlab实现的函数总结
在学习了线性回归与逻辑回归后，对其进行Matlab实现并总结如下：
（一） 运用Matlab内函数进行回归
%linear regression with matlab inner function
x=[1;2;3;4;5;6;7];
y=[2.1;5;5.8;8.2;10.5;11;15];
temp=ones(7,1);
X=[temp x];
b=regress(y,X)
z=b(1)+X*b(2);
plot(x,y,'or',x,z,'g');
xlabel('x');
ylabel('y');
title('Linear Regression with Matlab inner function');
这是用Matlab内部的函数进行实现的，实现图如下：
这里需要说明的是：这个例子实现的是一元一次的回归，当对多特征数据进行回归时，运用regress函数，同样可以实现线性回归。
另外，函数内部中的polyfit函数可实一元多项式的回归，而rstool可以实现多元二项式的回归，具体可参见Matlab help文档。
（二）运用梯度下降法进行Linear Regression
代码如下：
function [J,Gradient]=costfunction(theta)
x=[1;2;3;4;5;6;7];
y=[2.1;5;5.8;8.2;10.5;11;15];
[m,n]=size(x);
h=theta(1)+theta(2)*x;
J=sum((h-y).^2)/(2*m);
Gradient(1)=sum(h-y)/m;
Gradient(2)=sum((h-y).*x)/m;
function [optTheta,functionVal,exitFlag]=Gradient_descent( )
options = optimset('GradObj','on','MaxIter',100);
initialTheta = zeros(2,1);
[optTheta,functionVal,exitFlag] = fminunc(@costfunction,initialTheta,options);
效果如下：
与（1）相比，结果是一样的。
这里需要说明的是optimset函数与fminunc函数，应该认真看下help说明。（运用这两个函数，函数内部自动取值learning rate）
（三）Logistic Regression
由于逻辑回归的实现与线性回归的实现在运用梯度下降法实现时，处cost function不同外其他皆同，故代码省略。
如果Logistic Regression的实现过程中用Regularizaion，只需要对cost function进项稍加改动即可。
（四）小结
1）这两个实例都是对简单数据的线性回归。
对Linear Regression而言：当处理多特征的数据时可以采用（a）regress（或profit等）函数进行回归，（b）运用梯度下降法实现的过程中增加cost function的元
数，另外再增加各个变量的梯度即可。
对Logistic Regression而言：其本质是二值分类问题{0,1}，这样就遇到一个问题。（a）当实际的训练数据集是 two classes时，不管是是多特征的训练数据还是特征比较
少的数据，都可以参考Linear Regression的实现方法进行。（b）当实际的训练数据时mult-classes时，这时可以建立多个hypothesis
function，当predict新的数据时，取使得hypothesis function最大的那个即可得到。
2）再次重申下：
a）当训练样本的特征长度不一时，在对其进行回归前应进行Feature Scaling或者Mean Normalization，以更快的收敛；
b）在进行Logistic Regression时为避免ovrfit，通常都会（1）减少数据中无关的特征（2）Regularization
（五）Matlab函数总结
常见的三维绘图函数：
mesh（可画出立体网状图）
surf（立体曲面图）
peaks（为了方便测试立体绘图， MATLAB提供了一个peaks函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点，其方程式为：
meshz（可将曲面加上围裙）
waterfall（可在x方向或y方向产水流效果）、
meshc（同时画出网状图与等高线）、
surfc（同时画出曲面图与等高线）
contour3（画出曲面在三度空间
中的等高线）、
contour（画出曲面等高线在XY平面的投影）、
plot3（可画出三度空间中的曲线）
Reference：
【上篇】【下篇】[转载]arcgis&grid数据与matlab相互转换
目标：将arcgis
grid数据在matlab下导入，处理完后，重新输出，以ArcGIS9软件和Matlab&R2008a软件为工具。
　　　　首先利用ArcGIS9.0中Conversion Tool中的Raster to
Float工具，将栅格数据grid转化为float浮点数据，转化结果中包含两个文件，一个为flt扩展名，代表浮点数据，另一个为hdr扩展名，包含了栅格的行数和列数及栅格格子数值大小，还有左下角的x,y坐标等信息。用行数和列数参数，在matlab中使用fread函数，以float方式就可将数据读入,变成具体的数据矩阵变量。
在Matlab环境中编写数据处理函数，进行相应的计算处理后，利用matlab中的fprintf函数，将计算结果转为ASCII&数据，在输出时，应将原先在读取Arcgis的栅格数据时，获取的hdr文件行列数信息，及左下角的x,y坐标等参数，先写入文件头部，再输出计算出的栅格数值，这样，可以保证输出的ASCII&数据具有空间坐标，（这一点很重要，不然，无坐标信息，就相当于无价值了）然后重新在ArcGIS9.0环境下利用Conversion
Tool中的ASCII To Raster to Float工具导入，成为新的栅格信息图层。
具体代码如下：详细的说明之类，就不写了，有兴趣，Q
hdr文件示例：
原文件格式：
cellsize 25
NODATA_value
读取hdr头的函数：
读取头函数
function&[xmin,xmax,ymin,ymax,dx,nx,ny]=read_AGaschdr(hdr_file)
[names,values] =
textread(hdr_file,'%s%f',5);
nx=values(1);
ny=values(2);
xmin=values(3);
ymin=values(4);
dx=values(5);
xmax=xmin+(nx-1)*
ymax=ymin+(ny-1)*
读取实际数值函数：
function&[dsm]=read_AGAsc(flt_file,nx,ny)
fid=fopen(flt_file,'r');
[dsm]=fread(fid,[nx,ny],'float');
imshow(dsm,[min(dsm(:))
max(dsm(:))]);
实际运行命令：
[xmin,xmax,ymin,ymax,dx,nx,ny]=read_AGAscHead('c:tempnew1.hdr');
dsm=read_AGAsc('c:tempnew1.flt',nx,ny);
将计算结果输出为ascii文件的方法：
writeGrid2Arc(fileName,nz,mz,xllcorner,yllcorner,cellsize,Z)
%用于将matlab中的运算数据存到arcgis中的ascii文件中
%filename:用于存放的文件名，nz:行大小，mz:列大小，xllcorner:x左下角，
%yllcorneer:y左下角，cellsize：网格大小，
%Z要写入的数值。
%这样写出来的数据具有空间坐标。
%使用示例：
%writeGrid2Arc('D:workDirnew6.txt',00599,30,dsm);
fid=fopen(fileName,'wt');
dc=3; �fault number of
decimals to output
isnumeric(dc)
&dc=['%.',sprintf('%d',dc),'f'];
elseif isnumeric(dc)
&dc=['%.',sprintf('%d',dc),'d'];
%write header
fprintf(fid,'%st','ncols');
fprintf(fid,'
fprintf(fid,'%st','nrows');
fprintf(fid,'
fprintf(fid,'%st','xllcorner');
fprintf(fid,[dc,'n'],xllcorner);
fprintf(fid,'%st','yllcorner');
fprintf(fid,[dc,'n'],yllcorner);
fprintf(fid,'%st','cellsize');
fprintf(fid,[dc,'n'],cellsize);
fprintf(fid,'%st','NODATA_value');
fprintf(fid,[dc,'n'],-9999);
for i=1:mz
&for j=1:nz
&fprintf(fid,[dc,'n'],Z(i,j));
&fprintf(fid,[dc,'t'],Z(i,j));
&%update waitbar
&waitbar(i/mz,h,['Writing file:
',[fname,ext],...
&sprintf(' &%d%%
complete...',round(i/mz*100))])
fclose(fid);
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。MatLab 讲义2002 年 9 月版MATLAB 讲义第一章 MATLAB 系统概述1.1 MATLAB 系统概述MATLAB（MATrix LABoratory）矩阵实验室的缩写，全部用 C 语言编写。 特点： （1）以复数矩阵作为基本编程单元，矩阵运算如同其它高级语言中的语言变量操作一样方便，而且 矩阵无需定义即可采用。 （2）语句书写简单。 （3）语句功能强大。 （4）有丰富的图形功能。如 plot,plot3 语句等。 （5）提供了许多面向应用问题求解的工具箱函数。目前，有 20 多个工具箱函数，如信号处理、图 像处理、控制系统、系统识别、最优化、神经网络的模糊系统等。 （6）易扩充。1.2 MATLAB 系统组成（1）MATLAB 语言 MATLAB 语言是高级的矩阵、矢量语言，具有控制流向语句、函数、数据结构、输入输出等功能。同 时 MATLAB 又具有面向对象编程特色。MATLAB 语言包括运算符和特殊字符、编程语言结构、字符串、文件 输入/输出、时间和日期、数据类型和结构等部分。 （2）开发环境 MATLAB 开发环境有一系列的工具和功能体，其中大部分具有图形用户界面，包括 MATLAB 桌面、命令 窗口、命令历史窗口、帮助游览器、工作空间、文件和搜索路径等。 （3）图形处理 图形处理包括二维、三维数据可视化，图像处理、模拟、图形表示等图形命令。还包括低级的图形 命令，供用户自由制作、控制图形特性之用。 （4）数学函数库 有求和、正弦、余弦等基本函数到矩阵求逆、求矩阵特征值和特征矢量等。 MATLAB 数学函数库可分为基本矩阵和操作、基本数学函数、特殊化数学函数、线性矩阵函数、数学 分析和付里叶变换、多项式和二重函数等。 （5）MATLAB 应用程序接口（API） MATLAB 程序可以和 C/C++语言及 FORTRAN 程序结合起来，可将以前编写的 C/C++、FORTRAN 语言程序 移植到 MATLAB 中。1.3 MATLAB 的应用范围包括：MATLAB 的典型应用包括： ? 数学计算 ? 算法开发 ? 建模、仿真和演算 ? 数据分析和可视化 ? 科学与工程绘图 ? 应用开发（包括建立图形用户界面） 以矩阵为基本对象第二章 Matlab 基础2.1 MATLAB 快速入门（1）搜索路径 搜索路径也被看作是 MATLAB 的路径，其包含的文件被认为在路径上。搜索路径设置存放在文件 pathdef.m 中，称为当前目录，当要在 MATLAB 中打开一个文件时，就以当前目录为开始点。1MatLab 讲义2002 年 9 月版当输入一变量 value 时，MATLAB 的搜索路径次序： value 是否为变量 value 是否为内部函数 当前目录中是否存在 value.m 文件 搜索路径上是否存在 value.m 文件 path 函数可以控制 MATLAB 的目录搜索路径，主要使用的格式： path 显示当前的搜索路径 p=path 把当前的搜索路径存到字符变量 P 中 path('newpath') 设置路径为'newpath' path(path,'newpath') 向当前路径添加一个新目录 addpath 函数向 MATLAB 的搜索目录中添加一个新目录。 addpath 路径名 path(path,’路径名’)：增加搜索路径 rmpath 函数从 MATLAB 的搜索路径删除一个目录。 rmpath 路径名：删除路径 还可以利用菜单：File-&setpath（路径浏览器） what：显示出搜索路径上的文件名 what 路径名：路径名中的文件名 type value：显示变量内容 edit 文件名：对 m 文件进行编辑 （2）工作空间（Workspace） 工作空间是一个重要而且比较抽象的概念， 它是指运行 MATLAB 程序或命令所生成和存储在内存中的 所有变量和 MATLAB 提供的常量构成的集合。通过使用函数、运行 M 文件和装载保存的工作空间，可以向 工作空间增加变量。 ? save 保存整个工作空间或一部分变量，使用方式： save workspace as 文件名 或 save 文件名 [变量名] ? load 恢复工作空间，使用方式： load workspace load 文件名 ? 工作空间浏览器：File-&Show Workspace ? 还有一组命令来管理这些变量。 who,whos：显示出工作空间中的变量列表。 clear [变量名]：清除变量 （3）MATLAB 命令窗口 ? 输入命令和输出结果。 如输入：help [函数名] a=62.2 矩阵、变量、运算和表达式（1）矩阵的输入 A．直接输入： 注意： （1）行元素间用空格或逗号（， ）隔开； （2）行与行之间用分号（； ）或回车； （3）整个元素列表用[]括起。 直接输入的矩阵为一全局变量，一直保存在内存中。 例： a=[1 2 3;4 5 6] a= 1 2 3 4 5 6 a=[1,2,3;4,5,6;7,8,9] ? a=[1 2 3; 4 5 6; 7 8 9]2MatLab 讲义2002 年 9 月版矩阵元素：可以灵活地描述矩阵元素， ? 矩阵元素 a[i,j] 按列存放 通过下标单独对元素赋值 例：a(1,1)=1,a(3,2)=a(1,1) 得到 a = 1 a = 1 0 0 0 0 1 即自动形成一个 3 行 2 列矩阵，对未赋值的元素充值 0。 ? 矩阵的元素可以用任意形式的表达式 例：算术表达式 x=[-1,sqrt(5),(2+7)^4] x = 1.0e+003 * -0.2 6.5610 ? 大矩阵可以用小矩阵作为元素 例：a=[1 2;3 4] b=[a a+5;a-5 zeros(size(a))] 例：A=[1,2,3;4,5,6] A = 1 2 3 4 5 6 B=[A;7,8,9] B = 1 2 3 4 5 6 7 8 9 ? 可以从矩阵中抽取某些元素构成新矩阵 C=A(1:2,:) C = 1 2 3 4 5 6 例：a=[3,4,5;6,7,8] b=[+2,4*5,6] c=[sin(0.5*pi),sqrt(4),0] d=[a;b;c] ? 复数的表示 MATLAB 支持复数的运算，复数的虚部用 i 或 j 表示。 例：a=1+2i 或 a=1+2j 二者表示的结果一样。 复数可以直接运算， 例：a=3+4i; b=5+6j a+b 输出：ans= 8.0i 复数运算的一些常用函数： ①abs 返回复数的模 ②angle 返回复数的相角 ③conj 返回共轭复数 ④imag 返回复数的实部3MatLab 讲义2002 年 9 月版⑤real 返回复数的虚部 B.用语句或函数产生： a=randn(5,5) 产生正态分布 5*5 的随机矩阵。 C.用 M-文件或外部数据文件产生： M-文件是一个以.m 为后缀的文本文件，文件内容为一系列 MATLAB 命令，在 MATLAB 环境下键入该文 件名（不包括后缀） ，文件中的全部命令会依次逐个执行；M-文件名（不包括后缀）相当于一个宏命令. 例如：一个名为 magik.m 的文件包含了如下的内容， （假设 magik.m 在当前目录下） A = [ 16.0 3.0 2.0 13.0 5.0 10.0 11.0 8.0 9.0 6.0 7.0 12.0 4.0 15.0 14.0 1.0 ] 在 Matlab 环境下执行如下命令： magik A A = 16 3 2 13 5 10 11 8 9 6 7 12 4 15 14 1 D.用矩阵编辑器创建和修改矩阵： 使用 File-&Show workspace （2）矩阵运算 运算符 +，-，*，/(右除)，\(左除) 和^(幂)。 右除：C=A/B 即 C 满足 CB=A，当 B 可逆时，A/B=AB-1 左除：C=A\B 即 C 满足 AC=B，当 A 可逆时，A\B=A-1B 幂 A^n = A*…*A; A 必须是方阵。 例：矩阵的加减法： a=[1:3;4:6;7:9] b=a; c=a+b; c=a-b 注：矩阵相加减必须有相同的维数。 例：矩阵的点乘运算，^运算时矩阵必须为方阵，且只能与数字运算。 d=a*b 必须符合 m*n 与 n*l 的结构。 d=a.*b 矩阵的点乘运算 例：\(左除)：A\B=inv(A)*B,其中 inv(A)表示 A 逆阵，例如求解 AX=B。 A=[1 0 0；0 4 0；0 0 9]； B=[1 2 3；0 1 0；0 1 1]； X=A\B /(右除): A/B=A*inv(B),例如求解 XA=B。 X=B/A (3)变量与表达式 ? Matlab 的赋值语句有两种形式: 其一为：&变量&=表达式； 其二为：表达式，将表达式的值赋于一个自动定义的变量 ans。 注：A：如果以；结尾，则不显示计算结果，否则显示计算结果。 B：除保留字外，变量可以用字母开头，后跟 19 个字母或数字。变量名区分大小写，变量使用 时不需要先定义，也不必定义变量的类型。 ? 可以用 who 或 whos 来显示已定义的变量 例如： who Your variables are: A B C a ans4MatLab 讲义2002 年 9 月版whos Name Size Bytes Class A 2x3 48 double array B 3x3 72 double array C 2x3 48 double array a 3x2 48 double array ans 1x1 8 double array Grand total is 28 elements using 224 bytes ? 一些常用的变量 pi 3. //π值 i sqrt(-1 ) //虚数单位 j same as i eps floating-point relative precision, 2. realmin smallest floating-point number, 2. realmax largest floating-point number, 1. inf infinity (任意一个非零数除以 0) nan Not-a-number （0/0 或 inf-inf） 如： r=1/0 r=inf 1/r ans=0//容量变量 //最小浮点数 //最大浮点数 //正无穷大 //非数(4)矩阵的其他简单运算: A’: 矩阵转置 inv(A)：A-1 sum(A)：得到一个行向量，其元素为 A 的每一列的和 a=[1 2 3;4 5 6] sum(a) sum(a’) diag(A)：得到一个列向量，其元素为 A 的对角元 sum(diag(a)) 冒号(:)运算符： a:b:c：生成一个由等差数列构成的行向量 X，X(i+1)-X(i)=b 例：0:pi/4:pi ans = 0 0.8 2.6 如果省略 b，则等差数列的公差为 1 a=0:0.05:1 x=linspace(0,1,75) a=1:4;b=1:2:7;c=[b,a] 2 0.2 等比数列：logspace(0,2,11) 创建起点为 10，终点为 10 ，11 个元素，公比为 10 矩阵的变换：rot90: 矩阵逆时针旋转 n*90 度。 fliplr: 矩阵左右翻转。 flipud: 矩阵上下翻转。 稀疏矩阵的存储： sparse(A):用于把完全矩阵压缩为稀疏矩阵。 A=[0，1，0，0；0，3，0，4；5，0，0，0；0，0，0，7] sparse(A) ans= (3,1) 5 (1,2) 1 (2,2) 3 (2,4) 45MatLab 讲义2002 年 9 月版(4,4) 7 sparse(i,j,u):函数直接造成稀疏矩阵，i,j 为向量分别对应行号和列号，u 也为向量，存 储非元素的值. i=[1,2,2,3,4] j=[2,2,4,1,4] u=[1,3,4,5,7] A=sparse(i,j,u) full 函数把稀疏矩阵还原为完全矩阵。 (5)数组及其运算： 数组可以看作是行向量，实质为阵列运算。是元素对元素的运算，用句号（.）来区别。 数组和矩阵之间的区别在于运算规则不同，矩阵运算由线性代数规则来定义。 运算符：+,-和.*, ./, .\,.^ A.*B：A 与 B 对应的元素相乘 A.\B：B 的元素除以 A 的相应元素 A./B: A 的元素除以 B 的相应元素 A.^B:A 的元素为底，B 的相应元素为幂的数组 如：a=[1:3;4:6;7:9] b=a; c=a+b; c=a-b 查看下列运算的结果： a*b a.*b a/b a./b a\b a.\b a^b(指数和底数均为矩阵，无法求解) a.^b a’ a.’2.3 基本数学函数abs（绝对值或复数模） sqrt （平方根） real（复数的实部） imag （复数的虚部） conj（复数的共轭） round （舍入为最接近的整数） //round(-0.5)=-1 round(0.4)=0 fix （向 0 方向舍入为整数） //fix(0.99)=0 fix(1.01)=1 floor （向负无穷大舍入为整数） //floor(-0.5)=-1 floor(0.5)=0 ceil （向正无穷大舍入为整数） //ceil(-0.5)=0 ceil(0.6)=1 sign （符号函数） rem（x,y） （取余数函数） //得到 x/y 的余数，rem(11,4)=3 sin cos tan asin atan //三角函数都是面向阵列中的元素操作，角度单位均为弦度。 atan2(y,x) //-pi &= atan2(y,x) &= pi sinh （双曲正弦） cosh tanh exp （以 e 为底的指数） log （自然对数） log10 （常用对数） bessel （贝塞尔函数） gamma （伽码函数） rat （无理数的分式有理逼近）[N, D] = rat(x,tol), 要求：abs(x-N/D)&=tol*abs(x), tol 的缺 省值为 tol = 1.e-6*norm(X(:),1) 其中：norm( X ( : ),1) = max(sum(abs((X))))2.4 关于矩阵的一些有用的工具产生矩阵的工具: ? ‘[]’表示空矩阵 空矩阵不包含任何元素，它的维数为 0*0；空矩阵可以在运算中传递。 例：A=[1，2，3；4，5，6；7，8，9] A（2， ： ）=[]6MatLab 讲义2002 年 9 月版???? ??输出： A= 1 2 3 7 8 9 空矩阵有矩阵缩维的作用。 eye:单位矩阵：对角线元素为 1。 A=eye(3,3) A= 1 0 0 0 1 0 0 0 1 zeros：所有的元素都是 0 Z = zeros(2,4) Z = 0 0 0 0 0 0 0 0 ones：所有的元素都是 1 Z= 5*ones(3,3) Z = 5 5 5 5 5 5 5 5 5 rand： 矩阵元素是用均匀分布在[0.0,1.0]内产生的随机数 rand(n)：产生一个 n 阶（0-1）的随机矩阵 Z=rand(3) Z = 0.9 0.8 0.3 0.5 0.5347 rand(m,n): 产生一个 m*n 阶（0-1）的随机矩阵 Z=rand(1,6) Z = 0.8 0.2 0.3969 产生[-a,a]之间均匀分布的随机数的公式：R=a-2*a*rand(m,n) randn：矩阵元素是用期望为 0,方差为 1 的正态分布产生的随机数 Z=randn(3) Z =0.1927-0.7 1.6 -1.6 0.9 0.3273 产生均值为 b,方差为 q2 的正态分布随机数的公式：S=q*randn(m,n)+b ? 其它特殊矩阵： 例：a=[1 2 3;4 5 6;7 8 9] 1．diag(a) %% 生成矩阵 a 的对角矩阵 ans = 1 5 9 2．compan:伴随矩阵 例： u = [1 0 -7 6]7MatLab 讲义2002 年 9 月版A = compan(u) A = 0 7 -6 1 0 0 0 1 0 eig(compan(u)) ans = -3.0 1.0000 即方程的根。 3．magic:魔方矩阵（矩阵的元素由 1~10 组成，行、列和对角线之和相等） 例：magic(5) ans = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 4．tril(a) %% 生成矩阵 a 的下三角矩阵 ans = 1 0 0 4 5 0 7 8 9 5．triu(a) %%生成矩阵 a 的上三角矩阵 ans = 1 2 3 0 5 6 0 0 9 注：1.下标引用： 在 Matlab 中，矩阵元素的引用可用两个下标来表示,如矩阵 A 中，第 i 行第 j 列的元素用 A(i,j)引 用，也可以用一个下标来表示，用单个下标表示元素并不只限于矢量。对于矩阵，由于 Matlab 的运算基 本上都是对列操作的，矩阵可以认为是按列优先排列的一个长的列矢量，从而可用单下标引用。 例如 2*2 的矩阵，A(1)表示第一列的第一个元素，A(2) 表示第一列的第二个元素，A(3)表示第二列 的第一个元素，A(4)表示第二列的第二个元素，当矩阵的下标超出矩阵的实际元素的下标时，将给出错 误信息。但是，当某个值被赋给矩阵的一个新的元素时，Matlab 会自动增加矩阵的维数大小。例如 4*4 的矩阵 A，执行命令 A(4,5)=17 后，矩阵 A 将变成 4*5 的矩阵。 在下标的表达式里使用冒号表示矩阵的一部分。例如矩阵 A(1:k,j)表示矩阵 A 的第 j 列的前 k 个元 素。例如 sum(a(1:4,4))，表示计算第四列的前四个元素元素的和。 A(:,j)表示矩阵 A 的第 j 列的所有元素。由于有了冒号运算符，例如求矩阵 A 的第 j 列元素之和 ， 表示式为 sum(A(:,j)。 2．如何建立多维矩阵： 在科学研究与工程运算中，要建立多维矩阵，方法之一是扩展二维矩阵。 例：先建立一个简单的二维矩阵： a=[5 7 8;0 1 9;4 3 6] 为 a 矩阵扩展第三维 a(:,:,2)=[1 0 4;3 5 6;9 8 7] 方法之二，使用 cat 函数。它将现有的矩阵按照指定的维数建立新的多维矩阵。 B=cat(dim,A1,A2,...) 例:b=cat(3,[2 8;0 5],[1 3;7 9])8MatLab 讲义2002 年 9 月版2.5 复杂的矩阵运算一些矩阵运算： size(A)： 得到矩阵的阶，得到矩阵的 n 行 m 列. rank(A)： 得到矩阵的秩 det(A)： 行列式 eig(A)： 特征值和特征向量 e=eig(A) 得到特征值 [v,d]=eig(A): A*v=v*d svd(A)： 矩阵的奇异值分解[u,s,v]=svd(A), u,v 正交阵，s 对角阵，A=usv’ cond(A)：得到矩阵的条件数 （最大奇异值和最小奇异值的比值） lu(A): LU 分解（A 必须是一个方阵） [l,u]=lu(A)： u 上三角阵，l：一个下三角阵和一个置换矩阵的积 [l,u,p]=lu(A)：u 上三角阵，l：下三角阵，p：置换矩阵 qr(A)： QR 分解 [q,r]=qr(A)：r：上三角阵，q：正交阵2.6 矩阵输出格式格式命令（数字） X=[4/3 1.2345e-6]; format short：以定点数形式显示，小数后面保留 4 位有效数字 X X = 1.0 format short e：以浮点数形式显示,小数后面保留 4 位有效数字 X X = 1.. format short g：以定点数或浮点数的最佳形式显示，小数后面保留 4 位有效数字 X X = 1.5e-006 format long：以定点数形式显示，小数后面保留 14 位有效数字 X X = 1.33 0.00 format long e 以浮点数形式显示，小数后面保留 14 位有效数字 X X = 1.333e+000 1.000e-006 format long g 以定点数或浮点数的最佳形式显示，小数后面保留 14 位有效数字 X X = 1.33 1. format bank 以货币方式显示 X X = 1.33 0.00 format rat 用有理数近似表示 X X = 4/3 1/8100459MatLab 讲义2002 年 9 月版format hex 用 16 进制表示 X X = 3ff5 3eb4b6231abfd271 format compact 以紧凑形式显示 （format loose 以松散形式显示） X X = 3ff5 3eb4b6231abfd271 长命令行 用“…”来表示连接第三章 图形在分析问题的结果的时候，图形往往是一种很有用的工具，它可以帮助我们直观地了解结果的某些 性态。Matlab 提供了很强的图形功能，能够在根据向量或矩阵给定数据来生成图形。 Matlab 图形有数据可视化和图形处理两大功能，在数据的可视化部分，Matlab 可使用户计算所得的 数据根据其不同情况转化成相应的图形。用户可以选择直角坐标、极坐标等不同的坐标系；它可以表现 出平面曲线、空间曲线，绘制直方图、向量图、柱状图及空间网状面图、空间表面图等。 图形函数有 4 种： 通用图形函数、二维图形函数、三维图形函数、特殊图形函数. 建立一个图形的步骤： 步 骤 典型代码 1）准备数据 x=0:0.2:12 y1=bessel(1,x); y2=bessel(2,x); y3=bessel(3,x); 2)选择窗口，决定绘图位置 figure(1) subplot(2,2,1) 3)调用基本绘图函数 h=plot(x,y1,x,y2,x,y3); 4）设置绘图线条样式和标记 set(h,'LineWidth',2,{'LineStyle'},... {'--';':';'-.'} set(h,{'Color'},{'r';'g';'b'}) 5）设置坐标轴范围、刻度和栅格线 axis[0 12 -0.5 1]) grid on 6）标记图形坐标轴、图形图例以及其它文字 xlable('Time') ylable('Amplitude') legend(h,'First','Second','Third') title('Bessel Functions)' [y,ix]=min(y1); text(x(ix),y,'First Min\rightarrow'... 'HoizontalAlignment','right') 7）打印图形 print -dps23.1 图形窗口图形窗口（Figure Window）是所有 Matlab 的图形输出的专用窗口。figure:用来打开一个绘图窗口， 以供后续绘图命令输出图形。 创建图形窗口的命令：figure 两种格式：figure figure(n) 图形窗口的名称是按照该窗口创建的时间顺序依次命名的：Figure No.1,Figure No.2.....Figure No.n,命令 figure 将创建一个名为 figure No.n+1 的新的空白图形窗口，而 figure(n)命令则若 figure10MatLab 讲义2002 年 9 月版No.n 窗口已存在，则成为当前窗口。 h=figure(n)可得到图形窗口的句柄。 两个查阅图形参数和参数值的命令： get(n)：get(n)将返回关于图形窗口 Figure No.n 的所有图像参数的名称和当前值。 set(n)：set(n)将返回关于图形窗口 Figure No.n 的所有图像参数的名称和其可能取的值。3.2 创建二维图形（Plot 函数）Plot 是最简单而且使用最广泛的一个线型绘图函数，是 Matlab 的内部函数。 作用：可以生成线段、曲线和参数方程曲线图形。 几种基本命令形式： (1)向量式：plot(v) 参数 v 可以是向量、实数阵和复数阵。v 是长度为 n 的数值向量，生成(i,v(i))的一条折线，坐标轴 的范围由 MATLAB 系统根据向量的长度元素的大小自动生成，当向量的元素充分多时，即可以得到一条外 观光滑的曲线。 例：x=0:pi/20:2*plot(x) 例：x=plot(x) (2)参数式：plot(x,y) x,y 都是长度为 n 的向量，它的作用是在坐标系中生成顺序连接顶点{x(i),y(i)}的折线，绘得的连 线图以 x 为横坐标，y 为纵坐标。这种调用可以用来生成参数方程的图形。 例：y=sin(x);plot(x,y) 例：x1=0:1:10* x2=0:0.01:10* y1=sin(x1).*x1; y2=sin(x2).*x2; plot(x1,y1) hold on plot(x2,y2) （3）矩阵式：plot(y) y 是一个(m*n)和矩阵，为矩阵的每一列画出一条线。同时以矩阵的行向量为基准对 x 轴进行分度和 标注，标准时，采用向量 1:m，这里 m 中矩阵的行数。 例：z= %%产生一个 49*49 的矩阵 plot(z) （4）混合式：plot(x,y) 如果 X,Y 均为向量，则长度必须相等，亦即参数式；如果 X 是向量，而 Y 是一个矩阵，X 的长度与矩 阵 Y 的行数或列数相等，则它的作用是将向量 X 与矩阵 Y 的每列或每行的向量相对应作折（曲）线，当 Y 是方阵时，则将向量 X 与矩阵 Y 的列向量相对应作图；如果 X 是矩阵，Y 是向量，Y 的长度等于 X 的行数 或列数，则将 X 的每列或每行的向量与 Y 相对应作图。同样，当 X 是方阵时，则将 X 的各列与 Y 相对应 作图；如果 X 和 Y 都是矩阵，且维数相同，那么按列与列的对应方式来作图。 例：y=1:length(peaks);plot(peaks,y) （5）复向量式：plot(z) z 为复向量 ，会忽略向量的虚部，相当于 plot(real(z),imag(z)) 例:随机矩阵特征值分布 plot(eig(real(20,20),’o’,’Markersize’,6) （6）综合调用：plot(x1,y1,x2,y2,...) 矩阵中有多个矩阵调用对，其中的每一对按前四种方式之一进行调用，不同的矩阵对之间，其维数 可以不同。 例：t=0:pi/10:2*11MatLab 讲义2002 年 9 月版y1=sin(t); plot(t,y1); 可以使用 plot 同时显示由若干个 x-y 对表示的图形： y2=sin(t-0.5); y3=sin(t-0.25); plot(t,y1,t,y2,t,y3); 从图上我们可以看出，为了体现不同 x-y 对表示的图形，Matlab 自动地使用了不同的颜色；此外我 们还可以使用如下格式来设置图形的性态。3.3 线型、顶点标记和颜色在调用 plot 函数时 MATLAB 自动安排作图的线型的线段的颜色， 包括线段顶点的标记。 事实上， MATLAB 的 plot 函数可以设置和管理曲线的线段类型、顶点标记和颜色（在 MATLAB 中，它们通称为线型 Line Style） 。 MATLAB 定义的线段类型、顶点标记和线段颜色如下表： 颜色（color） 线型顶点标记(style) 顶点标记(marker) 类型 符号 类型 符号 类型 符号 黄色 y 实线 实点标记 . 洋红色 m 点线 : 圆圈标记 O 蛋青色 c 点虚线 -. 乘号标记 X X 红色 r 虚线 -加号标记+ + 绿色 g 星号标记* * 蓝色 b 方块标记□ s 白色 w 钻石形标记◇ d 黑色 k 向下的三角形标记 V 向上的三角形标记 ^ 向左的三角形标记 & 向右的三角形标记 & 五角星标记☆ p 六边形标记 h 注：1.连接节点的线型如果空缺则表示点与点之间没有直线相连； 2．如果不指定作图的颜色，自动循环使用 y,m,c,r,g,b,w 7 种颜色画线 plot 的最典型的调用方式是三元组参数： plot(x,y,’color-style-marker’) 例：plot(t,y,'b+', t, y2,'k-X', t,y3,'r--*') 例： [x,y,z]= //多输出函数，函数可产生多个输出值，输出值之间用逗号分开。 contour(x,y,z,20,’k’) hold on pcolor(x,y,z) shading interp 其中：peaks：根据 Gauss 分布（正态分布）得到 3 个 49 阶矩阵 contour：用相同的颜色画 20 条等高线， pcolor：将(x,y)点上的颜色设置成 z shading interp：设置渲染方式， 例：x=0:0.25:5; y1=x.^0.1; y2=x.^0.5; y3=x.^0.8; y4=x; y5=x.^1.5; y6=x.^2; y7=cos(x); y8=sin(x); hold on plot(x,y1,’yo’,x,y2,’mx’,x,y3,’c+’,x,y4,’rs’); plot(x,y5,’gh’,x,y6,’bd’,x,y7,’w&’,x,y8,’kp’);12MatLab 讲义2002 年 9 月版3.4 subplot 函数在一个已存在图形的绘图窗口中增加一个图形(在这之前讲的都是覆盖方式)hold on，加到图形上， hold off 则取代当前窗口。 允许在同一绘图窗口中显示或在同一纸上打印多个图形。 subplot(m,n,p)：将一个绘图窗口分割成 m*n 个子窗口，每一小块有自己的坐标轴，p 为活动窗口， 并且在 p 个窗口绘制当前图形，如果坐标系已存在，则 subplot(p)设置当前坐标系。 例如： t=0:pi/10:2* [x,y,z]=cylinder(4*cos(t)); %% (生成一个单位圆柱体，采用[0，1]之间的等距节点，每一个节点的旋转半径是 4*cos(t), 相当 于 4*cos(t)母线 CYLINDER(R,N), CYLINDER(R)(此时, n=20)) subplot(2,2,1); mesh(x) %%(x=1:m, y=1:n, z=x, c=x) subplot(2,2,2); mesh(y); subplot(2,2,3); mesh(z); subplot(2,2,4); mesh(x,y,z) subplot(111)是个特殊的情况，与 subplot(1,1,1)不同，它使下一条绘图指令在窗口中执行 elf 和 reset 指令（即删除当前图形的所有子对象） ，然后在默认位置创建一个坐标系。 绘制复数图形：plot(z) ? plot(real(z), imag(z)) 例如： t=0:pi/10:2* plot(exp(i*t),'-o')3.5 坐标轴控制 axis 函数axis 函数可控制坐标轴的刻度，使两个图形在对比时，有相同的比例因子。也可用于复杂图形的局 部透视。 axis 命令的格式： axis([xmin xmax ymin ymax]) 设置 X，Y 轴的极限范围 axis square： 保持 X，Y 轴 axis equal： 设置 X，Y 轴等长 axis auto： 设置 X，Y 轴为自动格式 axis on 显示 X，Y 轴 axis off 关闭 X，Y 轴 grid off 图形网格不显示 grid on 图形网格显示 例：复杂函数的局部透视，利用 axis 函数调整 x 轴刻度，可以比较清楚地看到某一局部区域。 x=0:1/3000:1;y=cos(tan(pi*x)); figure(1)，subplot(2,1,1),plot(x,y) title('复杂函数') subplot(2,1,2),plot(x,y) axis([0.4 0.6 -1 1]); title('复杂函数的局部透视') 例：利用 axis square 使输出图形的保持 x,y 轴，这样在图形窗口中绘制圆或椭圆都可用圆来表示。 t=0:pi/20:2* figure(1),subplot(2,2,1),plot(sin(t),cos(t)) title('圆形轨迹') subplot(2,2,2),plot(sin(t),2*cos(t)),title('椭圆形轨迹') subplot(2,2,3),plot(sin(t),cos(t)),axis square title('调节后的圆形轨迹') subplot(2,2,4),plot(sin(t),2*cos(t)),axis square title('调节后的椭圆形轨迹') 坐标标题及图中的文本使用函数：xlabel、ylabel、zlabel、text、title。13MatLab 讲义2002 年 9 月版说明： （1）xlabel、ylabel、zlabel 为 x,y,z 轴的标记。 （2）text 在当前坐标系中建立文本对象。 格式：text(x,y,’string’) 在图形的(x,y)点上放置指定的字符串。 text(...,’PropertyName’,’PropertyValue’,...)为标题文本指定特性。 特征为：∧表示为上标、_为表示下标、 \leftarrow 表示为← \rightarrow 表示为→ \it{斜体文字} arrow 表示箭头 bottom 表示上下箭头。 text(x,y,z,’string’) 在三维图形上放置文字 text(x,y,’string’,’color’,’k’) 放置文字的同时设置颜色。 例： t=-pi:pi/100: y=sin(t); plot(t,y) axis([-pi pi -1 1]) xlabel('-\pi \leq \itt \leq \pi') ylabel('sin(t)') title('Graph of the sine function') text(1,-1/3,'\it{Note the odd symmetry,}','color','r') （3）gtext 与 text 的不同是可利用鼠标放置文本。 例： t=-pi:pi/20: y1=sin(t);y2=2*cos(t); plotyy(t,y1,t,y2),grid on title('双 Y 轴正弦曲线') gtext('sin(t)') text(pi/2,0,'\leftarrow 2cos(t)') （4）plotyy：双 Y 轴的图形，同一张图上表示两条曲线时，可拥有各自的 Y 轴。 （5）title 给当前坐标系加上标题。 例：在同一张图上绘制出双 y 轴的 y1=sin(x)和 y2=2cos(x)函数。 x=-pi:pi/20: y1=sin(x);y2=2*cos(x); plotyy(x,y1,x,y2); grid on title('双 y 轴正余弦曲线') text(0,0,'\leftarrow sin(t)') text(pi/2,0,'\leftarrow 2cos(t)')3.6 专门用于绘制一元函数曲线的命令 fplot在 plot 命令中，系统是将从外部输入或间接确认的数值矩阵转化为连线图的。而在实际应用中绘制 函数的二维曲线时，一般并不清楚函数的具体情况，因而在确定自变量 x 的取值间隔时，往往一律用平 均间隔，这样往往不大准确。取好了，可以表现出函数图像的大概情况，取差了，则会因某处 x 元素的 间隔太大而根本错绘了曲线，不能反映出函数的变化情况从而使绘图失败。 fplot 不是直接从外部接收数据， 而是通过其内部自适应算法而产生的。 即在函数值变化比较平稳处， 它所取的数值点就会自动相对稀疏一些，在函数值变化剧烈处，所取有数值点就会自动密集一些。所以 对于曲线起伏剧烈的函数，用 fplot 命令将比用一般等间距取点的 plot 命令绘得的曲线光滑准确一些。 fplot 命令的具体使用格式为： fplot(‘function_name’,limits,tol,’linespec’,p1,p2,...)14MatLab 讲义2002 年 9 月版[x,y]=fplot(‘function_name’,limits,tol,’linespec’,p1,p2,...) 其中：‘function_name’：待绘制函数曲线的函数名称 limits ：limits=[xmin,xmax],为 x 的取值空间 tol ：为 fplot 命令在进行运算中的相对误差。tol 越小，所绘得的曲线就越接近 实际曲线的情况，但系统要为此占用很大的资源。 linepsec : 线型设置 p1,p2... : 函数传递的参数 x,y ：输出数据点坐标 例：fplot 与 plot 命令的比较。 function y=funfplot(x) y=sin(1./tan(pi.*x));%%建立一个函数 funfplot [X,Y]=fplot(‘funfplot’,[-0.1,0.1],2e-4); n=size(X); x=-0.1:0.2/(n(1)+1):0.1; y=funfplot(x); plot(x,y),hold on plot(X,Y)3.7 专门用于绘制一元符号函数曲线的命令 ezplot如果符号函数中只含有一个变量，那么 Matlab 将其视为一个普通函数，画出其在某个区域内的图像 用 ezplot 命令。 格式：ezplot(sym_function,limits) 其中：sym_function 符号函数或代表它的符号变量。 limits limits=[x1,x2],为 x 的取值空间。其默认值为[-2pi,2pi] ezplot(f) %%f=f(x) -2π&&x&2π ezplot(f,[min,max]) %% min&x&max ezplot(f,[xmin,xmax,ymin,ymax]) %%xmin&x&xmax,ymin&y&ymax ezplot(x,y) ezplot(x,y,[tmin,tmax]) ezplot(...,figure) 例： 1.ezplot(‘x^2-y^4’) 2.ezplot(‘(x^2)^(sin(x)^2)’)3.8 二维特殊图形MATLAB 5.3 提供了很多绘制二维图形的指令，如下图： 函数名称 功能 area 填充的函数折（曲）线图 bar 直方图 barh 垂直的直方图 bar3 三维直方图 bar3h 垂直的三维直方图 comet 慧星轨迹状的图形 errorbar 误差棒图 feather 沿 X 轴分布的复数向量图 fill 平面多边形填色 hist 向量的统计直方图 pareto 带有标准的直方图 pie 饼图 pie3 三维饼图15MatLab 讲义2002 年 9 月版plotmatrix 矩阵折（曲）线图 ribbon 带状图 scatter 点图（与 plot 的绘制结果相似，但是只有数据点） stem 火柴杆图 stairs 阶梯图 例 1：用 area 函数绘制面积图。area 函数将矢量数据或矩阵列矢量数据显示为曲线，并且将曲线与 x 轴之间的区间填充指定颜色。 用面积图描述某公司
年销售与利润情况。 sales=[51.6 82.4 90.8 59.1 47.0]; x=90:94; profits=[19.3 34. 61.4 50.5 29.40]; area(x,sales,'FaceColor',[.5 .9 .6],'EdgeColor','b','Linewidth',2) hold on area(x,profits,'facecolor',[.9 .85 .7],'edgecolor','y','linewidth',2) hold off set(gca,'Xtick',[90:94]);set(gca,'Layer','top'); gtext('\leftarrow sales');gtext('profits'),gtext('expenses') xlabel('years','fontsize',10), ylabel('Sales in 1,000''s','fontsize',10) 面积图例： Y = [ 1, 5, 3; 3, 2, 7; 1, 5, 3; 2, 6, 1]; area(Y)：grid on colormap summer set(gca,'Layer','top') title 'Stacked Area Plot' 例 2：用 bar 函数绘制向量 y 的直方图 x=0:pi/10:2* y=cos(x); bar(y); %%取 Y 元素的下标作为 X 坐标值 hold on bar(x,y,’r’); 将 bar 换成 barh 成垂直的直方图、bar3 成三维直方图、bar3h 成垂直的三维直方图。 例 3：用随机函数 rand 产生一个矩阵，从而得到复杂的条形图。 y=round(rand(5,3)*10); figure(1) subplot(2,2,1),bar(y,’group’),title(‘Group’) subplot(2,2,2),bar(y,’stack’),title(‘Stack’) subplot(2,2,3),barh(y,’stask’),title(‘Stack’) subplot(2,2,4),bar(y,1.5),title(‘width=1.5’) 例 4：用 errorbar(x,y,e)绘制误差棒图 x=1:10; y=sin(x); e=std(y)*ones(size(x)); errorbar(x,y,e) 注：x,y,e 必须是同维； 该指令常用于数理统计离散数据有理化。 例 5： hist 函数可在二维平面上绘制出柱状图， 用来表示数据值的分布情况。 其常用格式 hist(y,x)。 并且可以绘制笛卡尔坐标系和极坐标系两种方式。 分布统计直方图 x=-2.9:0.1:2.9;16MatLab 讲义2002 年 9 月版y=randn(20000,1); figure(1),hist(y,x) title(‘柱状图表示数据分布’) 用 rose 绘制 12 小时的风向图 wdir=[45 90 90 45 380 335 360 270 335 270 335 335]; wdir=wdir*pi/180; rose(wdir) 例 6：pie 函数可绘制出饼图，其常用格式为 pie(x)，可按向量 x 中的值的大小绘制出饼图。 pie(x,explode)可利用 explode 指定分离出的切片。 x=[1.1 2.8 0.5 2.5 2]; explode=[0 1 0 0 0]; figure(1) colormap hsv pie(x,explode) title(‘饼图’) 注：当 x 各元素和大于 1 时，函数对数据进行归一化处理；当各元素和小于 1 时，函数不对矢量进 行归一处理，而只是绘制部分饼图。如：pie([0.1,0.2,0.6]) 例 7：对于离散数据，采用 stem 函数。 alpha=0.02;beta=0.5;t=0:4:200; y=exp(-alpha*t).*sin(beta*t); stem(t,y) xlabel('Time in \musecs');ylabel('Magnitude') stem3 函数在三维空间描述离散序列。 th=(0:127)/128*2*x=cos(th);y=sin(th); f=abs(fft(ones(10,1),128)); stem3(x,y,f','d','fill') view([-65 30]); stairs 函数描述离散序列。 alpha=0.01; beta=0.5; t=0:10; f=exp(-alpha*t).*sin(beta*t); stairs(t,f);hold on plot(t,f,'--*');hold off3.9 网格和曲面MATLAB 具有强大的三维图形处理功能，包括三维数据显示、空间、曲面、分块、填充以及曲面光顺 着色、视点变换、旋转、隐藏等功能和操作。 （1）plot3(x,y,z,Linespec) plot3 是一个可用来画单变量的三维函数，函数格式除了包括第三维的信息之外，其它与二维函数 plot 相同。 例：绘制一个三维螺旋线。 t=0:pi/50:8* figure(1) plot3(sin(t),cos(t),t) grid off,axis square title('三维螺旋线') 例：增加维数的 plot3 命令可以使多个二维图形沿一个轴排列起来，而不是直接将二维图形叠加到另一 个的上面。下例就是二维图形在三维空间的排列。 x=linspace(0,3*pi); z1=sin(x); z2=sin(2*x);17MatLab 讲义2002 年 9 月版z3=sin(3*x); y1=zeros(size(x)); y3=zeros(size(x)); y2=y3/2; plot3(x,y1,z1,x,y2,z2,x,y3,z3); grid,xlabel(‘x-axis’),ylabel(‘y-axis’),zlabel(‘z-axis’) title(‘sin(x),sin(2x),sin(3x)’) (2)平面网格点的生成 在数学上，函数 z=f(x,y)的图形是三维空间的曲面，在 MATLAB 中，总是假设函数 z=f(x,y)是定义 在一个矩形的区域 D=[x0,xm]*[y0,yn]上的。为了绘制在区域 D 上的三维曲面，MATLAB 的方法是首先将 [x0,xm]在 x 方向上分成 m 份，将[y0,yn]在 y 方向分成 n 份，由各分划点分别作平行于坐标轴的直线， 将区域 D 分成 m*n 个小矩阵块，计算出在网格点的函数值。对于每个小矩形，在空间中决定出 4 个顶点 (xi,yi,f(xi,yi))， 连接 4 个顶点得到一个空间中的四边形片。 所有这些四边形片一起构成函数 z=f(x,y) 定义在区域 D 上的空间网格曲面。 为方便起见，MATLAB 中用 meshgrid 函数产生 x,y 轴向的网格数据，meshgrid 函数的一般格式为： [X,Y]=meshgrid(x,y) 其中，x,y 分别指向 x 轴向和 y 轴向的数据点，当 x 为 n 维向量，y 为 m 维向量时，X，Y 均为 m*n 矩阵，X(i,j)和 Y(i,j)共同指向平面上的一点，当 x,y 轴取同一向量时，可简写成 [X,Y]=meshgrid(x) 例：数学函数 ，定义在区域[-8,8]*[-8,8]上。在生成网格点后， 计算网格点上的函数值。 x=-8:0.5:8; y=x; [X,Y]=meshgrid(x,y); R=sqrt(X.^2+Y.^2)+ Z=sin(R)./R （3）网格曲面(mesh 函数) mesh 函数绘制三维空间上的网格曲面，如上例最后加上 mesh(Z)，生成一曲面。 mesh 函数的其它格式为： mesh(X,Y,Z): X,Y,Z 为同维数的矩阵。 mesh(X,Y,Z,C)：C 称为颜色矩阵。 网格曲面的顶点对应于空间的顶点(X(i,j),Y(i,j),Z(i,j)), 而网格曲面的网格线的颜色由 C 值根据当前的色谱来着色。这种调用形式还可以用来生成参数曲面片。 mesh(x,y,Z,C):其中，x 和 y 是向量，Z 和 C 是同维数的矩阵，且向量 x 的长度等于矩阵 Z 的列 数，而向量 y 的长度等于矩阵 Z 的行数。在这种情况下，网格曲面的曲格顶点是(x(j),y(i),Z(i,j))， 网格线的长度由矩阵 C 决定。 mesh(Z,C)：Z 和 C 都是 m*n 矩阵，该形式与 mesh(x,y,Z,C)等价。 与 mesh 相关的另外两个函数是 meshc 和 meshz，它们的调用形式与 mesh 相同。meshc 函数可同时绘 制出轮廓图（即等高线图） 。meshz 的作用除了生成与 mesh 相同的网和曲面之外，还在曲面下面加上一个 长方体的台柱，使图形更加美观。 例：曲面图 [x,y]=meshgrid(-3:.125:3); z=peaks(x,y);c=ones(size(z)); figure(1) mesh(x,y,z,c),grid on title('多峰函数的网格曲线') 例: [x,y]=meshgrid(-8:0.5:8); r=sqrt(x.^2 + y.^2) + /*加 eps 是为了防止出现 0/0 的情形*/ z=sin(r)./r; mesh(x,y,z); 例：曲面与等高线18MatLab 讲义2002 年 9 月版x=-8:0.5:8; y=x; [X,Y]=meshgrid(x,y); R=sqrt(X.^2+Y.^2)+ Z=sin(R)./R subplot(2,2,1);meshc(Z) subplot(2,2,2);meshz(Z)3.10 实曲面的绘制实曲面就是对网格曲面的网格块（四边形片）区域进行着色的结果。MATLAB 的函数 surf 函数可提供 这种功能。 它的调用方式与 mesh 相同， 但 mesh 仅对网格线着色， 网格线用黑色标出 （可以修改） 。 而 surf 是对网格片着色，绘制的是二维曲面。 例：用 surf 函数表现三维数据。 k=4; n=2^k-1; theta=pi*(-n:2:n)/n; phi=(pi/2)*(-n:2:n)'/n; X=cos(phi)*cos(theta); Y=cos(phi)*sin(theta); Z=sin(phi)*ones(size(theta)); colormap([0 0 0;1 1 1]);C=hadamard(2^k); surf(X,Y,Z,C); axis square 可以用 shading 函数来改变着色方式。 如：shading faceted (默认的着色方式，网格线为黑色)； shading flat (网格线分块着色)； shading interp (光顺性着色)。 网格块的内部像素的颜色由该片的四个顶点的颜色做双线性插值得出。 例：[x,y]=meshgrid(-3:.125:3); z=peaks(x,y);c=ones(size(z)); figure(1) shading faceted %% shading flat / interp surf(x,y,z,c),grid on title('多峰函数的网格曲面') surfc 函数绘制二维曲面同时还绘制出轮廓图（即等高线图） 。surfl 函数能生成具有光照效应的表 面，使图形更加美观。同时还可以指定三维坐标系的光源点坐标[x1,y1,z1]，或光源方向的球面坐标值， 即经度和纬度（仰角）向量[azimuth,elevation]。 例：用 surfl 指令创建一个图形，其中光源的方向为经度=-10 度，纬度=50 度。 surfl(peaks(200),[-10,50]); colormap(gray); shading flat 还有两个与 meshc 和 meshz 对应的函数 surfc 和 surfz，这里不再详述。3.11 等高线图形MATLAB 支持二维和三维等高线图形。 函数 contour 和 contour3 被用来实现这种功能。 它们将输入的 矩阵变量 M 看作是定义在[1,m]*[1,n]上的函数，生成若干条常数值的等高线段。用户也可以指定等高线 的条数、坐标系的比例以及某值上的等高线。 等高线的线型、顶点标记、颜色类型类似于 plot 函数定义。contour3 适用于三维等高线的生成，调 用方式与 contour 相同。 例：二维和三维等高线的生成。 contour(peaks,30) %生成二维等高线 % contour3(peaks,20) %生成三维等高线3.12 改变视角 view函数 view 改变所有类型的二维和三维图形的图形视角。与 z=0 平面所成的方向角叫仰角，与 x=0 平19MatLab 讲义2002 年 9 月版面的夹角叫做方位角，默认的三维视角方向仰角为 30 度，方位角-37.5 为度。而默认的二维视角仰角为 90 度，方位角为 0 度。 函数 view(az,el)和 view([az,el])将视角改变到所指定的方位角 az 和仰角 rl。 例：仰角和方位角的改变。 x=linspace(0,3*pi).’ z=[sin(x) sin(2*x) sin(2*x)]; y=[zeros(size(x)) ones(size(x))/2 ones(size(x))]; %生成矩阵 z,y subplot(2,2,1),plot3(x,y,z) grid on ,xlabel(“X-axis”),ylabel(“Y-axis”),zlabel(“z-axis”) title(‘默认视角：方位角=37.5，仰角=30’) view(-37.5,30) %将视角改变到所指定的方位角-37.5 和仰角 30 subplot(2,2,2),plot3(x,y,z) grid on ,xlabel(“X-axis”),ylabel(“Y-axis”),zlabel(“z-axis”) title(‘方位角旋转至 52.5 度’) view(-37.5+90,30) %将视角改变到所指定的方位角 52.5 和仰角 30 subplot(2,2,3),plot3(x,y,z) grid on ,xlabel(“X-axis”),ylabel(“Y-axis”),zlabel(“z-axis”) title(‘仰角变为 60 度’) view(-37.5,60) %将视角改变到所指定的方位角-37.5 和仰角 60 subplot(2,2,4),plot3(x,y,z) grid on ,xlabel(“X-axis”),ylabel(“Y-axis”),zlabel(“z-axis”) title(‘方位角为 0 度， 仰角为 90 度’) view(0,90) %将视角改变到所指定的方位角 0 和仰角 903.13 透视效应 hiddenMATLAB 在绘制图形时，在默认方式下，前面的图形会挡住后面的图形，即消去后面的隐藏线，可用 命令 hidden 改变这种模式。 命令格式：hidden off hidden on 例：透视图形. subplot(1,2,1), mesh(peaks(20)+7); hidden on subplot(1,2,2), mesh(peaks(20)+7); hidden off3.15 曲面的裁剪方法因为曲面图不能做成透明，但在一些情况下可以很方便地移走一部分表面以便看到表面下面的部分。 在 MATLAB 中，可以将要观察部分的数据置为特定的 NaN 来实现，由于 NaN 没有任何值，所有的作图函数 都忽略 NaN 的数据点。即将曲面裁剪掉一部分。 例：利用 NaN 进行曲面裁剪 [x,y,z]=peaks(30); x1=x(1,:);y1=y(:,1); i=find(y1&0.8 & y1&1.2); j=find(x1&-.6 & x1&.5); z(i,j)=nan*z(i,j); surf(x,y,z)3.16 四维表现和切片图 slice对于一般的定义在 x,y,z 坐标系上的四维可视化，可以用指令 slice 来实现。 为了实现三元函数 r=f(x,y,z)的可视化表现，MATLAB 提供了一个绘制三维物体切片图指令及与之配 合的三维网格坐标生成指令。 [X,Y,Z]=meshgrid(x,y,z) % 三维网格坐标的生成 slice(X,Y,Z,V,xi,yi,zi,n) % 绘制三维物体切片图20MatLab 讲义2002 年 9 月版其中： x,y,z :决定“网格”的位置，分别是(1*n)、(1*m)和(1*p)向量 X,Y,Z :三维网格坐标，它们都是(n*m)*p 维数组 V :在网线节点上的三元函数值数组，维数也为(n*m)*p xi,i,zi :分别决定垂直于 x,y,z 轴切面的位置向量。它们的维数可以不同。当取 0 维空阵时，表 示没有切面存在。 切片上的函数值的大小可以用色轴上对应的颜色表示。 V 数组中的最大有限值和最小有限值定义了色 轴的范围。又由于切片的位置可以任意设置，因此 slice 通过三维坐标点上的色彩变化把图形的表现能 力扩展到四维 -(x^2+y^2+z^2) 例：函数 R=xe 的四维表现图. x=-2:0.1:2; y=-2:0.25:2; z=-2:0.25:2; n=length(x); [X,Y,Z]=meshgrid(x,y,z); V=X.*exp(-X.^2-Y.^2-Z.^2); xi=[-0.7,0.7]; yi=0.5; zi=-0.5; slice(X,Y,Z,V,xi,yi,zi); xlabel(‘x’);ylabel(‘y’);zlabel(‘z’);hold on colorbar(‘horiz’) %表明颜色与数组的关系 view([30 45])3.17 颜色板（1）数字图象是指将一幅二维的图象表示成一个数值矩阵，矩阵的元素解释为像素的颜色值,为了 显示由矩阵表示的数字图像。将矩阵的每个元素对应到当前色谱的某个行标号，并取该行的颜色值作为 图像相应点的颜色。 （2）颜色板 map 是一个 m*3 的矩阵，其每行的值在 0.0-1.0 之间，分别表示红、绿、蓝三种颜色， 颜色板的每一行定义了一种颜色(对应于 RGB 三元色)。 （3）两维数组可以用 image 的形式显示出来，数组中的每一个元素（i,j）的值表示该点的颜色或 亮度，image 用来显求图象,包括装入调色板，生成二维图象。image(X)表示求生成矩阵 X 所表求的图像。 （4）colormap(map)将颜色板设置成 map。 colormap 颜色板矩阵 简单颜色 常用函数查找表函数 红 绿 兰 颜色 函数名称 说明 0 0 0 黑 bone 蓝色调灰色 1 1 1 白 cool 青和品红浓淡色 1 0 0 红 copper 线性变化纯铜色调 0 1 0 绿 flag 红白蓝黑交错色 0 1 1 青兰 gray 线性灰度 hot 黑红黄白色 hsv 带饱和值的颜色查找表 jet 颜色查找表的一种变体 pink 淡粉红色 prism 光谱色 colormap(‘default’) 将颜色板设置成缺省值. Colorbar 在当前的图形上显求一个水平的或的直的颜色标尺. cmap=colormap 得到当前使用的颜色板矩阵. colormap(hsv(128))产生 128 种颜色的 hsv 颜色板。hsv 的颜色从红、黄、绿、青、兰、洋红再回到 红，循环变化。 在 MAYLAB 中，图像由一个数据矩阵和一个颜色查找表矩阵组成，根据矩阵数据与像素颜色关系的不 同，可以将图像分为三种，索引图像、亮度图像和真彩图像。索引图像是将数据矩阵的元素视为颜色查 找表的索引；亮度图像的数据矩阵的元素为一个亮度值，该值线性对应于颜色查找表的索引；真彩图像 的数据矩阵是一个 m*n*3 的三维矩阵，第三维表示每一个像素点的 RGB 三色，而不用颜色查找表。 demo 目录中有一个名为 durer.mat 的文件，其中包含了两个矩阵 X, map。21MatLab 讲义2002 年 9 月版例：人体脊骨的图像, load spine image(X) colormap(bone);title(‘人体脊骨图’) 例： （ex23） %color(flag) cmap=l=length(cmap); x=[1:l];y=x'*ones(size(x)); bar(x(1:2),y(1:2,:)) colormap('summer') 例：hot(8)从黑色到暗红、洋红、黄色、白色平滑过渡。3.18 图形打印和读写print [Cdevice] [-option] &filename&：打印当前绘图窗口中的图形 其中：option 表示操作方式： 操作方式：-dps 单色 postscript 格式的文件 -dpsc 彩色 postscript 格式的文件 MATLAB 能够读写多种格式的图形文件，包括 TIFF、JPEG、BMP、PCX、XWD 和 HDF 格式。函数 imread 能够读取其中任意一种格式的图像文件。 例：RGB=imread('ngc6543a.jpg'); figure('position',[100 100 size(RGB,2), size(RGB,1)]); image(RGB); set(gca,'position',[0 0 1 1 ]); 例：load clown imwrite(X,map,'clown.bmp');3.19 动画制作用 MATLAB 产生动画序列有二种方法：第一种方法是先保存多幅不同的图片，然后连续回播；第二种 方法是连续不断地擦除并重画屏幕对象，在重画时屏幕对象不断变化，实现动画效果。这二种方法各有 优缺点。前者适于来不及快速重画的场合，它只是回放预先准备的画面。后者用到了画、擦、重画技术， 适于表现精度不高能够快速重画的场合。 例 1：使用 moviein 函数建立动画图形。 axis equal m=moviein(8); set(gca,'nextplot','replacechildren') for j=1:8 plot(fft(eye(j+8))) m(:,j)= end set(gca,'nextplot','replacechildren') movie(m,10) 例 2: h = uicontrol('style','slider','position',... [10 50 20 300],'Min',1,'Max',16,'Value',1) for k = 1:16 plot(fft(eye(k+16))) axis equal set(h,'Value',k)22MatLab 讲义2002 年 9 月版M(k) = getframe(gcf); end 例：使用擦除模式绘制动画图形。 a=[-8/3 0 0; 0 -10 10; 0 28 -1]; y=[35 -10 -7]'; h=0.01; p=plot3(y(1),y(2),y(3),'.','erasemode','none','markersize',5); axis([0 50 -25 25 -25 25]); hold on for i=1:4000 a(1,3)=y(2); a(3,1)=-y(2); ydot=a*y; y=y+h* set(p,'Xdata',y(1),'Ydata',y(2),'Zdata',y(3)) drawnow i=i+1; end第四章 控制流4.1 关系运算：MATLIB 有六种关系运算符：& &= & &= == ~= 结果 1 表示 true, 0 表示 false。 其中& &= & &=只用于操作数的实部比较，== ~=用于比较实部和虚部。 A.当两个变量是标量时，a 和 b 的关系成立，结果为 1，否则结果为 0。 B.当比较 a 和 b 是两个维数相同的数组时，按相同位置比较，结果是一个维数和 a 相同的数组，其 元素由 1 和 0 组成。 C.当比较的一个是数组 a， 一个是标量 b 时， 则把标量 b 和数组 a 的每一个元素按标量关系逐个比较， 结果是一个维数和 a 相同的数组，其元素由 1 和 0 组成。 如：a=rand(5);b=a c=a==b; c=a&b c=a&1;c=a&0;c=a&0.5 x=(1:10);t=x&54.2 逻辑运算：|（或）,&（与）,~（非）和逻辑函数 or、and、not、xor 及 any、all、findMATLAB 提供三个逻辑操作符&、|、~，同时又存在三个相应的 M 文件：all、or、not，这二组的作用 是相同的，只是使用格式稍有差异。xor 是第四个逻辑运算函数，完成异或操作。 在逻辑操作中，所有输入元素的非零值都当作 1 处理。逻辑运算的优先级最低。 例：x=[25 C5; 0 0.001] ~x ans= 0 0 1 0 y=[1 0;1 0] z1=x&y；z2=and(x,y)；z3=xor(x,y)； z1= 1 0 0 0 z2= 1 0 0 0 z3= 0 123MatLab 讲义2002 年 9 月版1 1 MATLAB 还提供了许多测试用的逻辑函数，巧妙地使用这些函数，可得到意想不到的结果。 any: 若作用于一向量 x，则当 x 中至少有一个非零元素时返回 1，否则返回 0，若 x 是矩阵，则对每 一列执行 any,返回一个元素为 0 或 1 的向量。主要作用测试矩阵中是否有非零元素。 如：a=rand(15) any(a)=(1 1 1 1 1) all:测试向量，若向量的所有元素非零，返回 1，否则返回 0，若对矩阵：则对每一列执行 all 操作。 主要作用测定矩阵中是否有零元素。 all(a) t=~(x&5) 例：ex24.m x=linspace(0,5,100) //产生 0~5 之间均匀分布的 100 个数据 y=cos(x) z=(y&0).*y //将 cos 函数的负数置为 0 z=z+0.3*(y&0) //将 cos 函数的正数值增加 30% z=(x&=4).*z //将超过 4 时的 z 值置为 0。 plot(x,z) find:找出矩阵中非零元素及其下标。 例：a=zeros(5,20); a(3,7)=0.5;a(4,15)=-0.4; [i,j,v]=find(a) i= 3 4 j= 7 15 v= 0.0 运算符的优先级： 在一个表达式内可能有几个运算符。MATLAB 为确定运算的次序设定了运算符的优先级。在同一 优先级时，程序先左后右执行，在优先级不同时，先高级后低级执行。 最高 （） ~（取反） . .^ + .* ./ .\ */ \+ : & &= & &= ++ ~=最低 & |4.3 if,else,elseif,end 语句if 语句的格式： if &逻辑表达式& 语句集 elseif &逻辑表达式& 语句集 else 语句集 end 仅由 if 和 end 组成的语句，可根据逻辑表达式的值选择是否执行。 例：if rem(a,2)=0 disp(’a is even’)24MatLab 讲义2002 年 9 月版b=a/2 end 例：if x disp(‘矩阵 x 全为 0’) else disp(‘矩阵 x 不全为 0’) end 例：if and(a==1,b&5) ...... end if 是 MATLAB 中最常用的条件执行语句，它与 end 语句一起构成各种格式，4.4．for-end 循环for 循环,指定次数的重复循环执行语句。 格式：for &循环参数& = &初态&:&步长&:&终态& 循 环 体 end 当步长=1 时，步长可以省略。 例：n=4; for i=1:n for j=1:n if i==j a(i,j)=2; elseif min([i,j])==1 a(i,j)=1; else a(i,j)=0; end end end disp(‘The matrix A is’);disp(a) 执行时： The matrix A is 2 1 1 1 1 2 0 0 1 0 2 0 1 0 0 2 例： （test1.m） for i=1:m for j=1:n a(i,j)=1/(i+j-1); end end 例：for i=9.8:3:-9 i end 执行时：i=9.800 i=6.800 3.800 i=0.800 i=-2.200 i=-5.200 例：s=’abcdefghijk’ for i=s i end %i 分别等于 s 中和每一个字符。25i=-8.200MatLab 讲义2002 年 9 月版例：c={‘aaa’,’bbb’,’ccc’} for i=c i end i=’aaa’ i=’bbb’ i=’ccc’ 例:在循环中,可以改变循环变量的值,但它不会改变循环的次数. r=1 for k=1:19 k=k+1; r=r*k; end 例：可利用数组(阵列)任意指定循环变是的值.(ex25.m) varx=[7 3 10 5]; vary=zeros(size(v)); k=0; for x=varx k=k+1 vary(k)=x.^2; end disp([v,w]) 例： （ex31.m）求 100-200 之间的素数。 for m=101:2:200 k=fix(sqrt(m)); for i=2:k+1 if rem(m,i)==0 end end if i&=k+1 disp(m) end end4.5 while-end 循环不定次数重复的循环执行语句。 格式：while &逻辑表达式& 循环体 end 例：求 1000 以下的菲波纳契数(ex26.m) f=[1 1]; i=1; while f(i)+f(i+1) & 1000 f(i+2) = f(i) + f(i+1); i=i+1; end 例：a=150; while a&0.1 a=a/2; end26MatLab 讲义2002 年 9 月版4.6 switch-end 语句：情况切换语句，case 语句中可以采用多个值。 语句格式： switch (变量或表达式) case v1 语句 case v2 语句 ....... otherwise 语句 end 例： switch (rem(n,4)==0) + (rem(n,2)==0) case 0 ...... case 1 ...... case 2 ...... otherwise error('This is impossible') end 例：switch var2 case {-2,-1} disp(‘var2 is negative one or two’) case {1,2,3} ...... otherwise disp(‘var2 is other value’) end4.7 break 语句：一个跳出循环的命令,导致最内层的 while,for,if 语句终止。 例： （test3.m） a = 0; fa = - b = 3; fb = while b-a & eps*b x = (a+b)/2; fx = x^3-2*x-5; if fx == 0 break elseif sign(fx) == sign(fa) a = fa = else b = fb = end end x27MatLab 讲义2002 年 9 月版4.8 try..catch 语句用于捕获程序运行中出现的错误，语法如下： try 语句,...语句,catch 语句,...语句 end 在正常情况下， 程序只执行 try 和 catch 之间的语句， 在执行上述语句发生错误时， 程序将执行 catch 和 end 之间的语句，在这里可以利用 lasterr 获得错误信息进行处理，如再发生错误，除非嵌套了另一 个 try...catch，否则程序将中断执行。4.9 return 语句return 语句终止当前的命令序列，把控制返回到调用函数或键盘。4.10 输入/输出语句? 获得用户输入 input：显示一段信息，并读取用户输入。 语法：x=input(提示信息,&输入类型&)，缺省为 double 型。 例：a=input(‘Please input a member’) a=input(‘Please input a char’,’s’) ? 暂停执行，等待用户按键 pause:程序暂停执行。 语法：pause(n) pause 其中：n 为暂停的秒数。 ? 建立完整的图形界面：菜单输入函数 menu。 例： （ex27.m） s=menu('color selection','red','greed','blue','yellow','black') switch(s) case 1 scolor='red'; case 2 scolor='green'; case 3 scolor='blue'; case 4 scolor='yellow'; case 5 scolor='blue'; otherwise disp('error'); end scolor 运行结果：28MatLab 讲义2002 年 9 月版? 输出函数 print,printopt4.11 自定义函数：用 M 函数自定义函数，函数 M 文件必须以函数的名称来作为文件名，M 文件的格式如下： &因变量&=&函数名&(&自变量&) 自变量和因变量都可以是矩阵或几个矩阵。 M 文件的第一行包括 function,该文件就是函数文件，函数和命令文件的区别是：命令文件的变量用 完后保存在内存，而函数文件文件内定义的量仅在函数文件内部起作用。函数文件有多个输出变量时， 用[]括起，有多个输出变量时，用（）括起。 例：(quroot.m) function x=quroot(a) if abs(a(1)) & eps d=a(2)^2-4*a(1)*a(3); if d &= 0 e=sqrt(d); f=2*a(1); x(1)=(-a(2)+e)/f; x(2)=(-a(2)-e)/f; end elseif abs(a(2)) & eps x=-a(3)/a(2); else error('coefficient error!') end 函数调用为： a=[1 2 1] x=quroot(a) x= -1 -1 /*子函数调用*/ function [mean, stdev] = stat(x) n=length(x); mean = avg(x,n); stdev = sqrt(sum((x-avg(x,n)).^2)/n); function mean = avg(x,n)29MatLab 讲义2002 年 9 月版mean = sum(x)/n; 在函数文件中可以包含多个函数，其中第一个函数称为主函数，其函数名与文件名相同，它可以由其 它文件引用，其它函数称为子函数，它只能由函数中的主函数和其它子函数引用。 函数中有两个永久变量nargin,nargout 表示引用函数时给出的输入变量数和输出变量数。 function c=testrg1(a,b) if(nargin==1) c=a.^2; else (nargin==2) c=a*b; end 输入一个变量时，计算这一变量的平方数 输入二个变量时，计算这两个变量的积。 函数的递归调用,在调用一个函数的过程中又出现直接或间接地调用该函数本身。 例：求 n!（ex33.m） function y=ex33(n) %%求 n！ if n&0 error('n is smaller them 0') return end if n==0|n==1 y=1 else y=n*ex33(n-1); end4.12 全局变量局部变量和全局变量： 函数工作空间中，变量有三类： （1）由调用函数传递输入和输出数据的变量； （2）在函数内部临时产生的变量（局部变量） ； （3）由调用函数空间，基本工作空间或其它函数工作空间提供的全局变量。 global：让不同的函数访问同一个变量 只要在函数中：global 变量名;%% 该变量为全局变量 function h = falling(t) global GRAVITY h = 1/2*GRAVITY*t.^2; 2 2 例： （ex28.m）对于函数 z=α(x-1) +β(y+1) 。编写相应的函数文件，其中α和β采用全局变量进 行参数传递。 global alpha beta x=[0:0.02:2]; y=[-2:0.02:0]; [X,Y]=meshgrid(x,y); subplot(2,2,1);alpha=1;beta=1; z=fun1(X,Y);mesh(z); subplot(2,2,2);alpha=2;beta=1; z=fun1(X,Y);mesh(z) subplot(2,2,3);alpha=1;beta=2;30MatLab 讲义2002 年 9 月版z=fun1(X,Y);mesh(z) subplot(2,2,4);alpha=0.8;beta=0.5; z=fun1(X,Y);mesh(z) function z=fun1(x,y) global alpha beta z=alpha*(X-1).^2+beta*(Y+1).^24.13 命令/函数的双重性在使用 Matlab 命令时，我们一般使用如下的格式：命令 参数 例如 load August17.dat help mesh 有时我们也可以使用如下的方式： load(‘August17.dat’) help(‘mesh’) 这就是 Matlab 中的命令/函数双重性，即 命令 参数 也可以写成函数形式 命令（’参数’） 从下面的例子里我们可以看出这种写法的好处： 例：有一批数据文件，文件名分别为：August1.dat, August2.dat,......august31.dat.希望分别 对这些文件做相应的处理： for i= 1:31 s=[‘August’ int2str(i) ‘.dat’] load(s) %process the contents of the ith file end 注：int2str为把数值变换为字符串。4.14 执行字符串命令? eval 函数 eval 是一个非常有用的函数，相当于一个命令解释器，能够解释一个字符串表达式，即计算以字符 串表示 matlab 表达式。 形式为： 1.eval(‘表达式’) 计算出表达式的值 2.[a1,a2,....an]=eval(‘表达式’) 当表达式可产生多个结果时， 分别存入 a1,a2,....an 等变量 中。 3.eval(字符串，选择字符串) 例： a=’4*5-3’ eval(a) ans= 17 例：生成 12 个随机矩阵，M1-M12, 阶数从 1-12 for n = 1:12 eval(['M' num2str(n) ' = rand(n)']) end ? feval 函数 可以把用户的输入和函数结合起来执行。 基本语法如下： feval(函数名称，数值或变量列表)； [y1,y2,...]=feval(function,x1,...,xn) 如果函数名称存在则将用给定的参数执行函数。 例：a=magic(4) [v,d]=feval('eig',a)31MatLab 讲义2002 年 9 月版feval 函数可以用于使一个函数接收设定为要调用的函数名的字符串变量，执行灵活的功能。 例：为多种函数画出输入输出图。 function plotf(fun,x) y=feval(fun,x); plot(x,y)4.15 时间和日期MATLAB 提供的时间和日期函数有四类： 函数简解： ? 返回当前时间的函数： now：返回当前日期与时间的数字序列。 date：返回当前日期字符串。 clock：返回当前日期与矢量时间。 ? 时间数据转换函数: datenum：将字符串日期转换成数字序列日期。 datestr：将数字序列日期转换成字符串形式。 datevec：将数字序列日期转换成矢量形式。 ? 计时函数: cputime：返回 MATLAB 使用的 CPU 时间（秒） 。 tic、toc：计时秒表。 etime：计算时间间隔。 ? 其它常用函数: calendar：日历。 weekday：查询某日为星期几。 eomday：月的最后一天。 ? 时间日期格式： 字符串日期：02-0ct-1996 数字序列日期：729300 矢量日期（系统内部格式） ： 0 0 0 格式转换函数：datenum(日期格式变数字)，datestr（日期格式变字符串）,datevec（日期格式变 矢量） 。4.12 向量化MATLAB 的实际使用中，有些循环可直接转换成向量操作，这样可大大提高程序执行速度，这种技术 称为循环的向量化，在编写程序时，应尽量避免采用循环，将它转变为向量进行处理。以便最大限度地 提高用 Matlab 编写的程序的效率。 如何将循环变量变换成向量操作是一个比较复杂的问题。简单问题很容易转换为向量操作，但复杂 的问题不易直观转换，视实际问题而定。 例：比较如下两个算法： （1） （t1.m） tstart= x=0.01; for k=1:1000 y(k) = log10(x); x = x + 0.01; end tend= interval = tend C tstart （t2.m） tstart =32MatLab 讲义2002 年 9 月版x=0.01:0.01:10; y=log10(x); tend = interval = tend C tstart t2.m的执行速度远远低于t1.m的执行速度。 例：产生一个由A矩阵重复产生的矩阵B function B=repmat1(A,M,N) if nagin&2 reeor(‘Requires at least 2 inputs.’) elseif nagin == 2 N=M end [m,n]=size(A);mind(1:m)’; mind(mind(;ones(1:M)); nind(nind(:,ones(1,N));disp(mind),disp(nind) B=A(mind,nind); 函数输入为：A[1 23;4 56]; b=repmat1(a,3,2)4.13 字符数组和字符串操作1．MATLAB的字符数组 在MATLAB中，字符串都被视作字符数组。 字符串的建立：s=’字符串’，即用’’将输入的字符串括起来。注意：不是” “,这点与其它高级语言 不同。而要建立一个字符串矩阵，则可以这样建立： sa=[‘string1’,’string2’......] 与数组不同，字符串矩阵的每一行字符元素的个数可以不同，但是每一行的所有字符串中的字符的 总个数必须相同，如果不满足这个条件，即使每行中字符串的个数相同，也会出错。事实上，MATLAB将 一个行内的所有字符串都合并起来，构成一个字符串，单个字符串之间不加空格，这正是每行中输入的 字符串的个数可以不同的原因。 char()函数自动将最短的字符串添加空格，使各行长度相等。 例: a=['aaa' 'bbb';'cc' 'dd'] 系统报错,必须为 a=['aaa' 'bbb';'cc' 'dddd'] 利用这个特点，可以用[]将任意字符串连接起来。 例：sa=[‘aa ‘bb ‘cc’;’dd’ ‘ee’ ‘f’ ‘k’] sa= aabbcc ddeefk 例：[sa(1,:) sa(2,:)] sa= aabbccddeefk 注：字符串的每个字符（包括空格）都是矩阵的一个元素。 处理字符矩阵时当作数据矩阵来处理 2．字符串的比较 函数strcmp(str1,str2)：比较两个字符串是否相同，若相等则返回1值，若不相等则返回0值。 函数strncmp(str1,str2,n):比较两个字符串的前n个字符是否相同，若相等则返回1值，若不相等则 返回0值。 3．字符串的归类 函数isletter(S):判断串中元素是否为字母。 函数isspace(S):判断串中元素是否为空格。 这二个函数的返回结果为与字符串长度相等的矢量，1表示结果为真，0表示结果为假。 4．字符串的查找与替换 函数findstr(‘str1’,’str2’)：在字符串str1中查找子串str2，返回str2在str1中的起始位置。33MatLab 讲义2002 年 9 月版函数strrep(‘str1’,’str2’,’str3’)：字符串替换 5．字符串运算 字符串及字符串矩阵可以进行加、减、乘、除四则运算和其它数字运算（利用ASCII码） 。 例：’e’+’f’=203 6.一般通用字符串操作 通用字符串操作包括字符串与ABCII间的转换、字符与数据间的相互转换，字符串大小写间的转换、 字符串中空格的删除等。 函数名 函数功能 string(A) 将一个整数数组转换成字符串矩阵 abs(S) 将字符串转换成ASCII码 double(S) 将字符串转换成相应的ASCII码（双精度） isstr(S) 确认是否为字符串 deblank(C) 删除字符串结尾处的空格 disp(n) 输入n个空格 str2mat(s1,s2,s3,...) 将字符串变为字符串矩阵 upper(S) 将字符串进行大写转换 lower(S) 将字符串进行小写转换 a=eval(‘’) 将字符串作为命令执行 strmatch(‘substr’,S) 匹配字符串操作 strtok(‘string’,d) 得到指定的子串 int2str(A) 将整数转换为字符串 num2str(A) 将浮点数转换为字符串 str2num(S) 将字符串转换为浮点数 dec2hex(A) 把十进制整数转换为十六进制字符串 hex2dec(S) 把十六进制字符串转换为十进制整数 hex2num(S) 把十六进制字符串转换为浮点数 dec2bin(A) 把十进制数转换为二进制字符串 例：打印出ASCII码从32-126之间所有可显示的字符。 char(32:126) char(48:126)4.14 类与对象由于引进了类和对象，MATLAB允许用户增加新的数据类型和操作，变量的类型描述变量的数据结构， 指示数据的操作。对象是类的实现。 MATLAB内建5种数据类。 数据类 说明 double 双精度浮点数数组或矩阵 sparse 二维实数（复数）稀疏矩阵 char 字符数组 struct 结构数组 cell 单元数组 要定义一个MATLAB类，需经过下面几方面工作： 1．建立类路径： 类路径下包含描述类的所有方法的M文件。类路径的名称由字符“@”加类名称构成。例如，定义类 ploynom的M文件都放置在一个名为@polynom的路径下面,我们用此类来描述多项式，通过调用函数构造函 3 数，由多项式的系数来确定多项式。如：p=polynom([1 0 -2 -5]),表示多项式p(x)=x -2x-5。 2.定义数据结构： 所有对象均存储在MATLAB结构中，结构的每一个字段以及字段的详细操作仅对操作类的方法是可见 的。如上例，将各次幂的系数以降幂的顺序作为一个行矢量来描述多项式。类polynom的对象p就成了一 个仅有一个字段p.c的结构。而这个字段仅有@polynom路径下的函数可以使用。 3．定义构造函数：34MatLab 讲义2002 年 9 月版类的方法路径下必须包含构造函数，构造函数的名称与类名称相同，构造函数的主要作用是为对象初 始化数据结构，并指定类标记。 下面是类polynom的构造函数@polynom/polynom.m: function p=polynom(a) if nargin==0 p.c=[]; p=class(p,'polynom'); elseif isa(a,'polynom') p=a; else p.c=a(:); p.class(p,'polynom'); end 4.定义类型转换函数： 转换函数的调用形式为：b=class_name(a) 这里a是一个非class_name类的对象。在MATLAB中最为重要的两个转换函数是double和char。这里， 我们为polynom类建立这样两个转换函数。建立M文件@polynom/double.m，实现从polynom类到double类 的转换。格式如下： function c=double(p) c=p.c 要实现类polynom到类char的转换，相对来说要麻烦得多。 function s=char(p) c=p.c; if all(c==0) s='0'; else d=length(p.c)-1; s=[]; for a=c; if a~=0; if ~isempty(s) if a&0 s=[s '+']; else s=[s '-']; a=-a; end end if a~=1|d==0 s=[s num2str(a)]; if d&0 s=[s '*']; end end end if d&=2 s=[s 'x^' int2str(d)]; elseif d==1 s=[s 'x'];35MatLab 讲义2002 年 9 月版end d=d-1; end end 此时，调用char(p)将的如下结果： char(p) ans=x^3-2*x-55.定义显示输出 定义显示输出的文件名为display。很多类的显示输出定义非常简单，仅仅是输出变量名，并将变量的值 转换为字符（调用函数char）输出。现为类polynom定义一个简单的显示输出文件@polynom/display.m. function display(p) disp(' '); disp([inputname(1),' =']); disp(' '); disp([' ' char(p)]) disp(' ');4.15 重载与继承重载包括三个方面的内容：算法的重载，操作符的重载和函数的重载。重载，实际是建立一系列特定 的M文件，以指定MATLAB在操作自定义类对象时的具体操作。如polynom对象加法的重载例子。有polynom 对象p,q，当遇到p+q时，将自动在@polynom路径下寻找plus.m文件，如果文件存在，执行文件。要重载 这个文件，实质就是建立这个文件。 function r=plus(p,q) p=polynom(p); q=polynom(q); k=length(q,c)-length(p.c); r=polynom([zeros(1,k) p.c]+[zeros(1,-k) q.c]);第五章 MATLAB 基本应用领域5.1 线性代数主要进行线性方程的求解、矩阵求逆、LU 分解、QR 分解、矩阵求幂、矩阵指数函数、求特征值及奇 异值分解。 1.MATLIB 中的矩阵 矩阵定义为二维实或复阵列。矩阵的加、减、乘、除、转置运算是最基本的运算，它们应符合维数 一致的限制。 2.线性代数方程求解 一般线性方程可表示成： AX=B XA=B 当矩阵 A 为方阵时，它的解为 X=A\B （X=inv(A)B）或 X=B/A。 当矩阵 A 为非奇异时，线性方程的解唯一；当矩阵 A 为奇异时，线性方程的解要么不存在，要么不 唯一。 a=[1 -2 -3 -4;2 1 -1 1;-1 0 -1 2;3 -3 4 -5]; b=[1 1 1 1]' x=a\b 当矩阵 A 为(m*n)维矩阵，且 m&n 时，方程 AX=B 中，方程个数多于变量个数，因此应采用最小二乘36MatLab 讲义2002 年 9 月版法来求解。 例：(ex32.m)对一组测量数据： t=[0 0.3 0.8 1.1 1.6 2.3]’ y=[0.82 0.72 0.63 0.6 0.55 0.5]’ 拟用延迟指数函数来拟合这组数据： y(t)≈c1+c2e-t 得 6 个方程，而未知数仅有 c1、c2 两个。应采用最小二乘法来求解。 A=[ones(size(t)) exp(-t)] c=A\y T=[0:.1:2.5]’; Y=[ones(size(T))exp(-T)]*c; plot(T,Y,’-‘,t,y,’o’) 3.矩阵求逆 det(A):函数可求矩阵 A 的行列式值。 inv(A): 函数可求矩阵 A 的逆矩阵值。 pinv(A):用于计算非方阵的伪逆。 4.LU、QR 分解 通过高斯对消或 LU 分解法，可将任意方阵表示成一个下三角阵与一个上三角阵的乘积：A=LU。 a=[1:3;4:6;4 2 6 ] [l,u]=lu(a) 得到 det(a)=det(l)*det(u) inv(a)=inv(u)*inv(l) 5.矩阵求幂和矩阵指数 2 3 矩阵求幂如 A ,B 可以用 A^2,B^3 sqrtm(A)：计算 A 的开方 A expm(A)：计算矩阵 A 的指数 e 6.特征值与特征向量 矩阵 A 的特征值λ和特征矢量ν，满足 Aν=λν。 eig(A)求矩阵 A 的特征值λ。 7.一元多项式的运算 (1)一元多项式的表示和创建 任意一个多项式都可以用一个行向量来表示,它的系数是按降序排列： 例:x^3-6*x^2+11*x-6 表示为 p=[1 C6 11 C6] poly2sym(p,’x’) 则 ans=x^3-6*x^2+11*x-6 poly2sym(p,’x’)函数：多项式的符号表示。 poly()函数：若 A 是 n*n 矩阵，则 p=poly(A)，p 为 A 的特征多项式，是 n+1 维向量,第一个元素的 值一定是 1。而特征多项式的根就是 A 的特征值。 u=[r1,r2,r3,......rm] poly(u)=(x-r1)(x-r2).....(x-rm) 例：A=[6 C11 6;1 0 0;0 1 0];p=poly(A); p=[1.0 C6.0 11.0 C6.0] roots(p) %%求多项式的根与 eig(A)相同。 例：求 1，2，-3 为根的多项式 u=[1 2 C3] p=poly(u) poly2sym(p) %% x^3-7*x+6 (b)多项式的基本运算 ? 多项式的加减运算：37MatLab 讲义2002 年 9 月版MATLIB 中没有多项式加减的函数，可以自编函数。 例： （qq.m） function yy=qq(x,y) nx=length(x); % x=reshape(x,1,nx) ny=length(y); %y=reshape(y,1,ny); n=max(nx,ny) cc=zeros(1,n); if nx&ny cc(1,(nx-ny+1):nx)=y;yy=x+ elseif nx&ny cc(1,(ny-nx+1):ny)=x;yy=y+ else yy=x+y; end 函数的调用为：p1=[3 4 6],p2=[5 2 C4 7] 则：求多项式的和为：p3=qq(p1,p2) 求多项式的差为：p3=qq(p1,-p2) ? 多项式的乘法运算 w=conv(u,v) 即求多项式 u 和 v 的乘积，即求向量 u 和 v 的卷积，公式如下：若 m=length(u);n=length(v);则 w 的长度为 m+n+1。 且 w(k)=u(1)v(k)+u(2)v(k+1)+....u(k)v(2k-1) 如:p1=[3 4 6];p2=[5 2 C4 7]; p3=conv(p1,p2); poly2sym(p3,’x’) %显示多项式15*x^5+26*x^4+26*x^3+17*x^2+4*x+42 两个以上的多项式的乘法需要重复使用conv. ? 多项式的除法运算 [q,r]=deconv(v,u)表示多项式 v 除以多项式 u 得到的商 q 和余数 r，q,r 均为多项式，进行解卷积 的运算。公式如下：例：[q,r]=deconv(p3,p1); poly2sym(q,’x’) %%q=5*x^3+2*x^2-4*x+7 [q,r]=deconv(p3,p2); poly2sym(q,’x’) %%q=3*x^2+4*x+6 ? 求多项式的根 roots(c)返回多项式的根组成的向量，也是多项式的友元阵的特征值向量，系数为实数的多项式的 根，若根为复数，则复数是成对出现的。 如 roots(p); 求 p=x^3-6x^2+15x-4 的根 p=[1 C6 15 C4],x=roots(p) x= 2.8494 + 2.2726i38MatLab 讲义2002 年 9 月版2.8494 - 2.1 ? 多项式的数组运算 按数组运算规则计算多项式的值 polyval。 y=polyval(p,x)计算多项式在x处的值，x可以是矩阵或向量。 x=[1:3;4:6];p=[1 2 3 4]; polyval(p,x) C．多项式的因式分解 n 阶多项式有 n 个根 r1,r2,...rn， 分解因式为(x-r1)(x-r2)...(x-rn)有可能有复根，复根则合并，因此在实数范围内的分解后的因式 的最高次数不会超过 2，分解因式后返回 n*3 的矩阵，其中 n 为分解后的因子个数，每一行就代表一个因 子。 例：多项式的因式分解的自编函数(expuly.m) function y=expuly(p) s=p(1);r=roots(p);y=[]; while ~isempty(r) c=r(1);r(1)=[]; if imag(c)&eps cc=conv([1 -c],[1 -conj(c)]); y=[y;cc]; r(find(abs(r-conj(c))&eps))=[]; else cc=real(c);cc=[0 1 -c];y=[y;cc]; end end if abs(s-1)&eps cc=[0 0 s];y=[y;cc]; end D.最大公约式和最小公倍式 最大公约式和最小公倍式是整系数多项式的两个比较重要的概念，MATLAB 本身并没有提供来求多项 式的最大公因式和最小公倍式的函数。自编函数 gcd.m，利用循环相除法求最大公约式。最小公倍式利用 两个多项式的乘积除以二者的最大公因式即可以得到最小公约式。 例：最大公约式和最小公倍式的自编函数（gcd.m） function y=gcd(f,g) while ~isempty(f) if abs(f(1))&2*eps f(1)=[]; end end while ~isempty(g) if abs(g(1))&2*eps g(1)=[]; else39MatLab 讲义2002 年 9 月版 end end if isempty(f)|isempty(g) error('Arguments should not be null'); end f=f/f(1);g=g/g(1); while 1 [q,r]=deconv(f,g); while ~isempty(r) if abs(r(1))&2*eps r(1)=[]; end end if isempty(r) y=g/g(1); return else f=g; g=r; end if length(r)==1 y=1; end end E．其它： ▲poly 函数 可以把向量转化为多项式，向量的元素为多项式的根。 u=[r1 r2 r3 ... rm] poly(u)为多项式 (x-r1)(x-r2)....(x-rm) 例：求 1 2 C3 为根的多项式 则 u=[1 2 C3]；p=poly(u) p=[1 0 C7 6] ▲residue 函数 [r,l,k]=residue 用于多项式的部分分式展开。 其中 p，q 分别为分子,分母多项式,r 为留数,l 为极点，k 为余项。 例：num=10*[1 4 5 6 7]; %分子多项式 den=poly([-2;-1;-0.5]); %分子多项式 [res,poles,ki]=residue(num,den) res = -6.040MatLab 讲义2002 年 9 月版64.1667 poles = -2.0 -0.5000 ki = 10 5 上面的结果说明了这个问题： 4 3 2 10(s +4s +5s +6s+7） -6.6667 = (s+2)(s+1)(s+0.5) s+2+-60.0000 s+1+64.16667 s+0.5+10s+ 5▲polyfit(x,y,n)函数，用于 n 次多项式来拟合该函数，n 的次数越高，拟合就越好。 例：用三次多项式来拟合指数函数。 x=0:0.1:5;y=exp(x);n=3;p=polyfit(x,y,n); f=polyval(p,x);plot(x,y,’b’,x,f,’r―‘) 例：x=1:5;y=[5.5 43.1 128 290.7 498.4] p=polyfit(x,y,3) x1=1:0.1:5; y1=polyval(p,x1) plot(x,y,’o’,x1,y1);grid on ▲多项式函数 函数名称 功能简介 conv(a,b) 乘法 [q,r]=deconv(a,b) 除法 poly(r) 用根构造多项式 polyder(a) 对多项式或有理多项式求导 polyfit(x,y,n) 多项式数据拟合 polyval(p,x) 计算 x 点处多项式值 [r,p,k]=residue(a,b) 部分分式展开式 [a,b]=residue(r,p,k) 部分分式组合 roots(a) 求多项式的根5.2 泛函分析--非线性数值方法函数MATLAB 提供了一些有关积分、常微分方程求解等方面的许多泛函（可对函数进行操作的函数） ，如找 出函数在区间上的极小值、求函数的零值点、计算函数积分等。 优 化 与 求 根 fzero 单变量函数的零点（即方程之根