您的当前位置:首页正文

数值分析上机实验指导书

2023-02-19 来源:步旅网
“数值计算方法”上机实验指导书

实验一 误差分析

实验1.1(病态问题)

实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。通过本实验可获得一个初步体会。

数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。

问题提出:考虑一个高次的代数多项式

p(x)(x1)(x2)(x20)(xk)k120(1.1)

显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。现考虑该多项式的一个扰动

p(x)x190(1.2)

其中是一个非常小的数。这相当于是对(1.1)中x19的系数作一个小的扰动。我们希望比较(1.1)和(1.2)根的差别,从而分析方程(1.1)的解对扰动的敏感性。

实验内容:为了实现方便,我们先介绍两个MATLAB函数:“roots”和“poly”。

uroots(a)

其中若变量a存储n+1维的向量,则该函数的输出u为一个n维的向量。设a的元素依次为a1,a2,,an1,则输出u的各分量是多项式方程

a1xna2xn1anxan10

的全部根;而函数 bpoly( v的输出b是一个n+1维向量,它是以n维向量v的各分量为根的多项式的系数。可见“roots”和“poly”是两个互逆的运算函数。

ess0.000000001;vezeros(1,21);ve(2)ess;roots(poly(1:20)ve)

上述简单的MATLAB程序便得到(1.2)的全部根,程序中的“ess”即是(1.2)中的。 实验要求:

(1)选择充分小的ess,反复进行上述实验,记录结果的变化并分析它们。如果扰动

项的系数很小,我们自然感觉(1.1)和(1.2)的解应当相差很小。计算中你有什么出乎意料的发现?表明有些解关于如此的扰动敏感性如何?

(2)将方程(1.2)中的扰动项改成x18或其它形式,实验中又有怎样的现象出现? (3)(选作部分)请从理论上分析产生这一问题的根源。注意我们可以将方程(1.2)

写成展开的形式,

p(x,)x20x190(1.3)

同时将方程的解x看成是系数的函数,考察方程的某个解关于的扰动是否敏感,与研究它关于的导数的大小有何关系?为什么?你发现了什么现象,哪些根关于的变化更敏感?

思考题一:(上述实验的改进)

在上述实验中我们会发现用roots函数求解多项式方程的精度不高,为此你可以考虑用符号函数solve来提高解的精确度,这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考MATLAB的帮助。

思考题二:(二进制产生的误差)

用MATLAB计算0.1100。结果居然有误差!因为从十进制数角度分析,这一计算

i11000应该是准确的。实验反映了计算机内部的二进制本质。

思考题三:(一个简单公式中产生巨大舍入误差的例子) 可以用下列式子计算自然对数的底数

1(1)n ee1limnn1这个极限表明随着n的增加,计算e值的精度是不确定的。现编程计算f(n)(1)n与

nexp(1)值的差。n大到什么程度的时候误差最大?你能解释其中的原因吗?

相关MATLAB函数提示:

poly(a) 求给定的根向量a生成其对应的多项式系数(降序)向量 roots(p) 求解以向量p为系数的多项式(降序)的所有根 poly2sym(p) 将多项式向量p表示成为符号多项式(降序) sym(arg) 将数字、字符串或表达式arg转换为符号对象 syms arg1 arg2 argk 将字符arg1,arg2,argk定义为基本符号对象 solve('eq1') 求符号多项式方程eq1的符号解

实验二 插值法

实验2.1(多项式插值的振荡现象)

问题提出:考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多,插值多项式的次数就越高。我们自然关心插值多项式的次数增加时,Ln(x)是否也更加靠近被逼近的函数。龙格(Runge)给出一个例子是极著名并富有启发性的。设区间[-1,1]上函数

f(x)1 2125x实验内容:考虑区间[-1,1]的一个等距划分,分点为 xi12i,i0,1,2,,n n则拉格朗日插值多项式为 Ln(x)1lx( )2ixii0125n其中的li(x),i0,1,2,,n是n次拉格朗日插值基函数。 实验要求:

(1) 选择不断增大的分点数目n=2,3….,画出原函数f(x)及插值多项式函数Ln(x)在[-1,1]上的图像,比较并分析实验结果。

(2)选择其他的函数,例如定义在区间[-5,5]上的函数

h(x)x,g(x)arctanx 1x4重复上述的实验看其结果如何。 (3)区间[a,b]上切比雪夫点的定义为 xkbaba(2k1)cos,k1,2,,n1 222(n1)以x1,x2,xn1为插值节点构造上述各函数的拉格朗日插值多项式,比较其结果,试分析原因。

实验2.2(样条插值的收敛性)

问题提出:多项式插值是不收敛的,即插值的节点多,效果不一定就好。对样条函数插值又如何呢?理论上证明样条插值的收敛性是比较困难的,但通过本实验可以验证这一理论结果。

实验内容:请按一定的规则分别选择等距或者非等距的插值节点,并不断增加插值节

点的个数。考虑实验2.1中的函数或选择其他你有兴趣的函数。

实验要求:

(1)随节点个数增加,比较被逼近函数和样条插值函数误差的变化情况。分析所得结果并与拉格朗日多项式插值比较(可以用MATLAB的函数“spline”作此函数的三次样条插值,取n=10、20,分别画出插值函数及原函数的图形)。

(2)样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题:某汽车制造商用三次样条插值设计车门的曲线,其中一段的数据如下: xk yk yk’

要求:i.自己编程计算(用三弯矩、三转角方程均可) ii.主函数myspline(x,y,边界类型,边界值,xi ) 其中:x 节点 y 节点上的函数值 xi 未知节点 返回:S(xi)

iii.三对角方程组用追赶法求解(书P160)。

实验2.3:(一维插值的应用—画图) 画你自己的手的形状,在MATLAB中输入

figure('position',get(0,'screensize')) axes('position',[0 0 1 1]) [x,y]=ginput;

0 0.0 0.8 1 2 3 4 5 6 7 8 9 10 0.2 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 将你的手掌张开放在计算机屏幕上,然后使用计算机鼠标选取一系列点勾勒出手的轮廓,按回车键结束ginput过程,这样就获得了一系列你的手掌外形数据点(x,y)。也可以这样获得数据点(x,y),先把手放在一张白纸上,并用笔画出它的轮廓线,然后将纸贴在计算机屏幕上,透过纸能看到平面上的鼠标,并通过ginput记录下轮廓上的点。

将x和y坐标值看作是两个独立变量的函数,独立变量的取值为从1到记录的点的数目。利用MATLAB的插值函数进行插值,并画出你的手掌外形轮廓。

思考题:(二维插值)

1. 在一丘陵地带测量高程,x和y方向每隔100米测一个点,得高程数据如下。试用MATLAB的二维插值函数“interp2”进行插值,并由此找出最高点和该点的高程。 y x 100 200 300 400 100 636 697 624 478 200 698 712 630 478 300 680 674 598 412 400 662 626 552 334 2. 画出山区地貌图。

利用MATLAB的peaks函数生成某山区的一些地点及其高度三维数据(单位:m)。命令格式:[x,y,z]=peaks(n),生成的n阶矩阵x,y,z为测量的山区地点三维数据(n>=30)。根据peaks函数生成的数据,利用Matlab二维插值画出该山区的地貌图和等值线图(提示函数:interp2、meshgrid、plot3等)。

相关MATLAB函数提示: plot(x,y) 作出以数据(x(i),y(i))为节点的折线图,其中x,y为同长度的向量 subplot(m,n,k) 将图形窗口分为m*n个子图,并指向第k幅图 yi=interp1(x,y,xi) 根据数据(x,y)给出在xi的分段线性插值结果yi pp=spline(x,y) 返回样条插值的分段多项式(pp)形式结构 pp=csape(x,y,‘边界类型’,‘边界值’) 生成各种边界条件的三次样条插值 yi=ppval(pp,xi) pp样条在xi的函数值 ZI=interp2(x,y,z,xi,yi) x,xi为行向量,y,yi为列向量,z为矩阵的双线性二维插值 ZI=interp2(…,'spline') 使用二元三次样条插值 ZI=griddata(x,y,z,xi,yi) x,y,z均为向量(不必单调)表示数据,xi,yi为网格向量的三角形线性插值(不规则数据的二维插值) 6

实验三 函数逼近与曲线拟合

实验3.1(曲线逼近方法的比较)

问题提出:曲线的拟合和插值,是逼近函数的基本方法,每种方法具有各自的特点和特定的适用范围,实际工作中合理选择方法是重要的。

实验内容:考虑实验2.1中的著名问题。下面的MATLAB程序给出了该函数的二次和三次拟合多项式。

x=-1:0.2:1;

y=1/(1+25*x.*x); xx=-1:0.02:1; p2=polyfit(x,y,2); yy=polyval(p2,xx); plot(x,y,’o’,xx,yy); xlabel(‘x’); ylabel(‘y’); hold on;

p3=polyfit(x,y,3); yy=polyval(p3,xx); plot(x,y,’o’,xx,yy); hold off;

适当修改上述MATLAB程序,也可以拟合其他你有兴趣的函数。 实验要求:

(1)将拟合的结果与拉格朗日插值及样条插值的结果比较。 (2)归纳总结数值实验结果,试定性地说明函数逼近各种方法的适用范围,及实际应用中选择方法应注意的问题。

实验3.2:(最小二乘拟合的经验公式和模型)

1.(已知经验公式):某类疾病发病率为y‟和年龄段x(每五年为一段,例如0~5岁为第一段,6~10岁为第二段……)之间有形如yae的经验关系,观测得到的数据表如下 x 1 2 y 0.898 2.38 3 3.07 12 6.53 4 1.84 13 10.9 5 2.02 14 16.5 6 1.94 15 22.5 7 2.22 16 35.7 8 2.77 17 50.6 9 4.02 18 61.6 19 81.8 bxx y 10 4.76 11 5.46 实验要求: bxyae(1)用最小二乘法确定模型中的参数a和b(提示函数:lsqcurvefit,

lsqnonlin)。

bxyae(2)利用MATLAB画出离散数据及拟合函数图形。

7

(3)利用MATLAB画出离散点处的误差图,并计算相应的均方误差。 2.(最小二乘拟合模型未知) 某年美国轿车价格的调查资料如表,其中xi表示轿车的使用年数,yi表示相应的平均价格,实验要求:试分析用什么形式的曲线来拟合表中的数据,并预测使用4.5年后轿车的平均价格大致为多少? xi yi 1 2 3 4 5 6 7 8 9 10 2615 1943 1494 1087 765 538 484 290 226 204 实验3.3(研究最佳平方逼近多项式的收敛性质)

x实验内容和要求:取函数f(x)e,在[-1,1]上以勒让德多项式为基函数,

(x)f(x)pn(x)对于n0,1,,10构造最佳平方逼近多项式pn(x),令n,将n(x)~x的曲线画在一个图上。

令Enmaxf(x)pn(x),画出En~n的曲线。做出En~n之间的最小二乘

1x1曲线,能否提出关于收敛性的猜测。

思考题一:(病态)考虑将[0,1]30等分节点,用多项式y1xx5生成数据,再用polyfit求其3次、5次、10次、15次拟合多项式,并分析误差产生的原因。

思考题二:调研Matlab拟合工具箱的使用。

相关MATLAB函数提示: p=polyfit(x,y,k) 用k次多项式拟合向量数据(x,y),返回多项式的降幂系数,当k>=n-1时,实现多项式插值,这里n是向量维数 8

实验四 数值积分与数值微分

实验4.1 (高斯数值积分方法用于积分方程求解)

问题提出:线性的积分方程的数值求解,可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数,与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下,高斯数值积分方法有较高的代数精度,用它通常会得到较好的结果。

实验内容:求解第二类Fredholm积分方程

y(t)k(t,s)y(s)dsf(t),atb

ab首先将积分区间[a,b]等分成n份,在每个子区间上离散方程中的积分就得到线性代数方程组。

实验要求:分别使用如下方法,离散积分方程中的积分 1.复化梯形方法;2.复化辛甫森方法;3.复化高斯方法。求解如下的积分方程。

21tttey(s)dsee1)y(t),方程的准确解为; 0e1112)y(t)esty(s)ds1e2tet,方程的准确解为et;

02 比较各算法的计算量和误差以分析它们的优劣。

实验4.2(高维积分数值计算的蒙特卡罗方法)

问题提出:高维空间中的积分,如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式,可以用一元函数积分的数值方法来计算高维空间的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。

实验内容:对于一般的区域,计算其测度(只要理解为平面上的面积或空间中的体积)的一般方法是:先找一个规则的区域A包含,且A的测度是已知的。生成区域A中m个均匀分布的随机点pi,i1,2,,m,如果其中有n个落在区域中,则区域的测度m()为n/m。函数f(x)在区域上的积分可以近似为:区域的测度与函数f(x)在中n个随机点上平均值的乘积,即

f(x)dxm()1f(pk) npk 实验要求:假设冰琪淋的下部为一锥体而上面为一半球,考虑冰琪淋体积问

题:计算锥面z2x2y2上方和球面x2y2(z1)21内部区域的体积。如果使用球面坐标,该区域可以表示为如下的积分:

/42cos02002sinddd

用蒙特卡罗方法可以计算该积分。

另一方面,显然这样的冰琪淋可以装在如下立方体的盒子里

1x1,1y1,0z2

9

而该立方体的体积为8。只要生成这个盒子里均匀分布的随机点,落入冰琪淋锥点的个数与总点数之比再乘以8就是冰琪淋锥的体积。比较两种方法所得到的结果。

类似的办法可以计算复杂区域的测度(面积或体积)。试求由下列关系所界定区域的测度:

0x1,1y2,1z3(1)exy sin(z)y01x3,1y4(2)x3y329 yex20x1,0y1,0z1(3)x2sinyz xzex1

相关MATLAB函数提示: diff(x) 如果x是向量,返回向量x的差分;如果x是矩阵,则按各列作差分 diff(x,k) k阶差分 q=polyder(p) 求得由向量p表示的多项式导函数的向量表示q Fx=gradient(F,x) 返回向量F表示的一元函数沿x方向的导函数F'(x),其中x是与F同维数的向量 z=trapz(x,y) x表示积分区间的离散化向量;y是与x同维数的向量,表示被积函数;z返回积分的近似值 z=guad(fun,a,b,tol) 自适应步长Simpson积分法求得Fun在区间[a,b]上的定积分,Fun为M文件函数句柄,tol为积分精度 z=dblquad(fun,a,b,c,d,tol,method) 求得二元函数Fun(x,y)的重积分 z=triplequad(fun,a,b,c,d,e,f,tol,method) 求得三元函数Fun(x,y,z)的重积分 10

实验五 解线性方程组的直接方法

实验5.1 (主元的选取与算法的稳定性)

问题提出:Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss消去法作为数值算法的稳定性呢?Gauss消去法从理论算法到数值算法,其关键是主元的选择。主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。

实验内容:考虑线性方程组

Axb,ARnn,bRn

编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss消去过程。

实验要求:

61786115(1)取矩阵A,b,则方程有解x*(1,1,,1)T。

861158614取n=10计算矩阵的条件数。让程序自动选取主元(顺序消元),结果如何?

(2)现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。

(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。

(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。重复上述实验,观察记录并分析实验结果。

实验5.2(线性代数方程组的性态与条件数的估计) 问题提出:理论上,线性代数方程组Axb的摄动满足

xcon(dA)x1A1AAbAb 矩阵的条件数确实是对矩阵病态性的刻画,但在实际应用中直接计算它显然不现实,因为计算A1通常要比求解方程Axb还困难。

实验内容:MATLAB中提供有函数“condest”可以用来估计矩阵的条件数,

它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端,使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动A和b,使得A和b充

11

分小。

实验要求:

ˆbb,以1-范数,(1)假设方程Ax=b的解为x,求解方程(AA)x给出

xxˆxxx的计算结果。

(2)选择一系列维数递增的矩阵(可以是随机生成的),比较函数“condest”

所需机器时间的差别.考虑若干逆是已知的矩阵,借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。

(3)利用“condest”给出矩阵A条件数的估计,针对(1)中的结果给出

xx的理论估计,并将它与(1)给出的计算结果进行比较,分析所得结果。注意,如果给出了cond(A)和A的估计,马上就可以给出A1的估计。 (4)估计著名的Hilbert矩阵的条件数。

H(hi,j)nn,hi,j1,ij1i,j1,2,,n

思考题一:(Vadermonde矩阵)设

11 A11x0x1x2xn2x0x122x22xnnix00nix0nxin1x1i0nx2,bni,

x2i0nxnnixni0其中,xk10.1k,k0,1,,n,

(1)对n=2,5,8,计算A的条件数;随n增大,矩阵性态如何变化?

(2)对n=5,解方程组Ax=b;设A的最后一个元素有扰动10-4,再求解Ax=b (3)计算(2)扰动相对误差与解的相对偏差,分析它们与条件数的关系。 (4)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗?

12

相关MATLAB函数提示: zeros(m,n) 生成m行,n列的零矩阵 ones(m,n) 生成m行,n列的元素全为1的矩阵 eye(n) 生成n阶单位矩阵 rand(m,n) 生成m行,n列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量x的元素构成的对角矩阵 tril(A) 提取矩阵A的下三角部分生成下三角矩阵 triu(A) 提取矩阵A的上三角部分生成上三角矩阵 rank(A) 返回矩阵A的秩 det(A) 返回方阵A的行列式 inv(A) 返回可逆方阵A的逆矩阵 [V,D]=eig(A) 返回方阵A的特征值和特征向量 norm(A,p) 矩阵或向量的p范数 cond(A,p) 矩阵的条件数 [L,U,P]=lu(A) 选列主元LU分解 R=chol(X) 平方根分解 Hi=hilb(n) 生成n阶Hilbert矩阵 13

实验六 解线性方程组的迭代法

实验6.1(病态的线性方程组的求解)

问题提出:理论的分析表明,求解病态的线性方程组是困难的。实际情况是否如此,会出现怎样的现象呢?

实验内容:考虑方程组Hx=b的求解,其中系数矩阵H为Hilbert矩阵,

H(hi,j)nn,hi,j1,ij1i,j1,2,,n

这是一个著名的病态问题。通过首先给定解(例如取为各个分量均为1)再计算出右端b的办法给出确定的问题。

实验要求:

(1)选择问题的维数为6,分别用Gauss消去法、列主元Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何?

(2)逐步增大问题的维数(至少到100),仍然用上述的方法来解它们,计算的结果如何?计算的结果说明了什么?

(3)讨论病态问题求解的算法

实验6.2 书上P211计算实习题2,其中N=100,第二小问改为用Jacobi迭代、G-S迭代、红黑排序的G-S迭代求解,并比较他们之间的收敛速度;进一步,用BSOR迭代求解,试找出最优松弛因子。

14

实验七 非线性方程求根

实验7.1(迭代法、初始值与收敛性)

实验目的:初步认识非线性问题的迭代法与线性问题迭代法的差别,探讨迭代法及初始值与迭代收敛性的关系。

问题提出:迭代法是求解非线性方程的基本思想方法,与线性方程的情况一样,其构造方法可以有多种多样,但关键是怎样才能使迭代收敛且有较快的收敛速度。

实验内容:考虑一个简单的代数方程

x2x10

针对上述方程,可以构造多种迭代法,如

2xn1xn1(7.1)

(7.2)

xn111xnxn1xn1(7.3)

在实轴上取初始值x0,请分别用迭代(7.1)-(7.3)作实验,记录各算法的迭代过程。

实验要求:

(1)取定某个初始值,分别计算(7.1)-(7.3)迭代结果,它们的收敛性如何?重复选取不同的初始值,反复实验。请自选设计一种比较形象的记录方式(如利用MATLAB的图形功能),分析三种迭代法的收敛性与初值选取的关系。

(2)对三个迭代法中的某个,取不同的初始值进行迭代,结果如何?试分析迭代法对不同的初值是否有差异?

(3)线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异,有何结论和问题。

相关MATLAB函数提示: x=fzero(fun,x0) 返回一元函数fun的一个零点,其中fun为函数句柄,x0为标量时,返回在x0附近的零点;x0为向量[a,b]时,返回函数在[a,b]中的零点 [x,f,h]=fsolve(fun,x0) 返回一元或多元函数x0附近fun的一个零点,其中fun为函数句柄,x0为迭代初值;f返回fun在x的函数值,应该接近0; h返回值如果大于0,说明计算结果可靠,否则不可靠 15

实验八 常微分方程初值问题数值解法

实验8.1(Lorenz问题与混沌) 问题提出:考虑著名的Lorenz方程

dxdts(yx)dy(8.1) rxyxzdtdzxybzdt其中s,r,b为变化区域有一定限制的实参数。该方程形式简单,表面上看并无惊人之处,但由该方程提示出的许多现象,促使“混沌”成为数学研究的崭新领域,在实际应用中也产生了巨大的影响。

实验内容:先取定初值y0=(0,0,0),参数s=10,r=28,b=8/3,用MATLAB的数值求常微分方程函数ods45编程对(8.1)进行求解

实验要求:

(1)对目前取定的参数值s,r和b,选取不同的初值y0进行运算,观察计算的结果有什么特点?解的曲线是否有界?解的曲线是不是周期的或趋于某个固定点?

(2)在问题允许的范围内适当改变其中的参数值s,r,b,再选取不同的初始值y0进行运算,观察并记录计算的结果有什么特点?是否发现什么不同的现象?

思考题一:考虑如下常微分方程

()k0k1cosmt() x''rx'bsinax先固定其中参数:r=0.9,b=3,a=1,k0=0,m=1,初始值为原点,t的范围为0-500

分别取不同的k1,如1,1.5,2,2.5,3,3.5等,绘图观察相位图(x,x')的变化情况,能找到其中的规律吗?

试试改变其它参数呢?

16

相关MATLAB函数提示: [t,y]=ode45(odefun,tspan,y0) odefun表示f(t,y)的函数句柄,t是标量,y是标量或向量;tspan是二维向量[t0,tf],表示自变量初值t0和终值tf;y0表示初值向量,若无输出参数,则作出图形 17

附录

MATLAB简介

这里介绍MATLAB一些入门知识,包括MATLAB桌面和窗口,MATLAB命令格式、数据格式、数据文件和变量管理,MATLAB的数组和矩阵运算,MATLAB的字符串、元胞和结构等数据类型,MATLAB的程序设计方法,MATLAB作图方法在线帮助的使用和程序文件和目录的管理等。

表一 MATLAB的基本命令 主题词 含义 主题词 含义 format 设置数据显示格式 feval 函数求值 who 显示变量名 input 提示输入 whos 显示变量信息 disp 输出 clear 清除内存变量 tic 启动秒表 save 保存工作变量到文件 toc 时间读数(秒) load 从文件装载变量 help 帮助 linspace 区间等分 lookfor 查找 length 获取数组长度 type 列程序清单 size 矩阵大小 which 查找文件目录 max 最大值 double 双精度 min 最小值 str2num 字符串转化为数值 sum 求和 num2str 数值转化为字符串 find 条件检索 一、MATLAB桌面

启动MATLAB后,就进入MATLAB的桌面,图1为MATLAB6.1的默认(Default)桌面。第一行为菜单栏,第二行为工具栏,下面是三个最常用的窗口。右边最大的是命令窗口(Command Window),左上方前台为发行说明书窗口(Launch pad),后台为工作空间(Workspace),左下方为命令历史(Command History)后台为当前目录(Current Directory)。

1.窗口

(1)命令窗口

该窗口是进行MATLAB操作最主要的窗口。窗口中“>>”为命令输入提示符,其后输入运算命令,按回车键就可执行运算,并显示运算结果.。

18

图1

(2)发行说明书窗口

发行说明书窗口是MATLAB所特有的,用来说明用户所拥有的Mathworks公司产品的工具包、演示以及帮助信息。

(3)工作空间

在默认桌面,位于左上方窗口前台,列出内存中MATLAB工作空间的所有变量的变量名、尺寸、字节数。用鼠标选中变量,击右键可以打开、保存、删除、绘图等操作。

(4)当前目录

在默认桌面,位于左下方窗口后台,用鼠标点击可以切换到前台。该窗口列出当前目录的程序文件(.m)和数据文件(.mat)等。用鼠标选中文件,击右键可以进行打开、运行、删除等操作。

(5)命令历史(Command History)

该窗口列出在命令窗口执行过的MATLAB命令行的历史记录。用鼠标选中命令行,击右键可以进行复制、执行(Evaluate Selection)、删除等操作。

除上述窗口外,MATLAB常用窗口还有编程器窗口、图形窗口等。

二、数据和变量

1.表达式

在命令窗口作一些简单的计算,就如同使用一个功能强大的计算器,使用变量无须预先定义类型。

19

4例如,设球半径为r=2,求球的体积Vr3。

3>>r=2 %表达式将2赋予变量r r= %系统返回r的值

2

>>v=4/3*pi*r^3 %pi为内置常量,乘方用^表示 v=

33.5103

几个表达式可以写在一行,用分号(;)或逗号(,)分割,用分号(;)使该表达式运算结果不显示,而逗号(,)则显示结果。也可以将一个长表达式分在几行上写,用三点(…)续行。

若需要修改已执行过的命令行,可以在命令历史中找到该命令行复制,再粘贴至命令窗口修改。也可以直接使用键盘↑↓调出已执行过的命令行修改。

2.数据显示格式

MATLAB默认的数据显示格式为短格式(short):当结果为整数,就作为整数显示;当结果是实数,以小数点后四位的长度显示。若结果的有效数字超出一

定范围,以科学计数法显示(如3.2000e-006表示3.2106)。数据显示格式可使用命令Format改变。例如:

>> format long;v %长格式,16位 v =

33.51032163829112

>> format short;v %短格式 v =

33.5103

>> format rational;v %有理格式,近似分数 v =

6501/194 3.复数

MATLAB中复数可以如同实数一样,直接输入和计算。例如: >> a=1+2i;b=5-4*i;c=a/b c =

-0.0732 + 0.3415i 4.预定义变量

MATLAB有一些预定义变量(表1),启动时就已赋值,可以直接使用,如前我们使用的圆周率pi和虚数单位i.

表1 常用预定义变量 变量名 说 明 i或j 虚数单位1 pi 圆周率3.14159… eps 浮点数识别精度2^(-52)=2.22041016 realmin 最小正实数2.225110308 realmax 最大正实数1.797710308 inf 无穷大

20

NaN 没有意义的数 预定义变量在工作空间观察不到。如果预定义变量被用户重新赋值,则原来的功能暂不能使用。当这些用户变量被清除(clear)或MATLAB重新启动后,这些功能得以恢复。

5.用户变量

MATLAB变量名总以字母开头,以字母、数字或下划线组成,区分大小写,有效字符长度为63个。如A,a,a1,a_b都是合法的,且a与A表示不同变量。在Command Window中使用的变量一旦被赋值,就会携带这个值存在于工作空间,直到被清除或被赋予新的值。

ans是系统一个特别的变量名。若一个表达式运算结果没有赋予任何变量,系统自动用ans存放答案。例如:

>> A=5+4i;b=5-4*i;B=1;A*b %没有定义A*b的输出变量 ans = 41 %ans来接受计算结果,注意这是大写A与小写b的乘积,尽管我们可以使用工作空间来查询和清除变量,但使用下列命令方式更快捷:

>> whos %查询Workspace中的变量列表 Name Size Bytes Class

A 1x1 16 double array (complex) B 1x1 8 double array

a 1x1 16 double array (complex) ans 1x1 8 double array

b 1x1 16 double array (complex) c 1x1 16 double array (complex) Grand total is 6 elements using 80 bytes

>> A %查询变量A的值 A =

5.0000 + 4.0000i >> clear A %清除变量A >> A %再查询A的值,已经不存在了 ??? Undefined function or variable 'A'. >> clear %清除Workspace中所有变量 >> whos %Workspace中已没有任何变量了

三、数组和矩阵运算

MATLAB基本数据单元是无需指定维数的数组。数组运算是MATLAB最鲜明的特点,一方面可以使得计算程序简明易读,另一方面可以提高计算速度。

1.数组的输入

最常用的数组是双精度数值数组(double array)。一维数组相当于向量,二维数组相当于矩阵,一维数组可以视为二维数组的特例。二维数组的第一维称为“行”,第二维称为“列”。MATLAB数组无需预先定义维数。直接输入数组的

21

元素,用中括号([])表示一个数组,同行元素间用空格或逗号分隔,不同行间用分号或回车分隔,例如:

>> clear;a=[1,2,3;4,5,6;7,8,9] a =

1 2 3 4 5 6 7 8 9 或

>> a=[1 2 3 %这种方式特别适用于大型矩阵 4 5 6 7 8 9] a =

1 2 3 4 5 6 7 8 9

对于等差数列构造的一维数组,可用冒号运算生成,也可用函数linspace生成。

>> b=0:3:10 %初值:增量:终值 b =

0 3 6 9 >> b=0:10 %增量为1可省略 b =

0 1 2 3 4 5 6 7 8 9 10

> >b=10:-3:0 %递减 b =

10 7 4 1 >> b=linspace(0,10,4) %将区间[0,10]等分为4-1=3份 b =

0 3.3333 6.6667 10.0000 >> length(b) %查询b的长度 ans = 4 >> b(3) %查询b的第三个元素 ans =

6.6667 >> b([1,end]) %查询b的首和尾元素 ans =

0 10

二维数组元素双下标编址按通常方式,单下标编址按列排序。 >> size(a) %查询数组a的尺寸 ans =

3 3 >> a(3,2),a(6)

22

ans = 8 ans = 8

>> c=a([1 3],[2 3]) %提取a的第一、第三行和第二、第三列(分块矩阵)

c =

2 3 8 9 >> d=a(2,:) %提取a的第二行 d =

4 5 6 >> a(:) %将a所有元素按单下标顺序排为列向量ans = 1 4 7 2 5 8 3 6 9

一些特殊的二维数组可以用函数产生,例如: >> a=zeros(2,4) %生成2行4列零矩阵 a =

0 0 0 0 0 0 0 0 >> b=ones(1,4) %生成1行4列1矩阵 b =

1 1 1 1 >> c=[a;b] %拼接 c =

0 0 0 0 0 0 0 0 1 1 1 1 >> c(2,1)=100 %修改部分元素 c =

0 0 0 0 100 0 0 0 1 1 1 1 >> reshape(c,2,6) %按2行6列重排矩阵元素 ans =

0 1 0 0 1 0

23

100 0 1 0 0 1

注意:数组下标对应矩阵的行和列,编址一律从1开始,不能用0. 矩阵输入也可用“load”命令从外部数据文件导入 2.数组运算

数组运算是指数组对应元素之间的运算,也称点运算。矩阵的乘法、乘方和除法有特殊的数学含义,并不是数组对应元素的运算,所数组乘法、乘方和除法的运算符前特别加了一个点。特别要区分数组运算在乘法、乘方和除法上的意义和表示上与矩阵运算的不同。

表2 数组运算符 运算 符号 说明 数组加与减 A+B与A-B 对应元素之间加减 数乘数组 k*A或A*k k乘A的每个元素 数与数组加减 k+A或k-A k加(减)A的每个元素 数组乘数组 A.*B 点运算只有点乘、点乘方、点除三个,表示对应元素之间的运算;(.*)是一个整数组乘方 A.^k,k.^A 体,点(.)不能漏掉,(.)和(*)之间也数除以数组 k./A 数组除法 左除A.\\B,右除B./A 不能有空格 >> clear;A=[1 -1;0 2];B=[0 1;1 -1]; >> A.*B %注意不是A*B ans =

0 -1 0 -2 >> A.\\B,A./B

Warning: Divide by zero. ans =

0 -1.0000 Inf -0.5000 Warning: Divide by zero. ans =

Inf -1 0 -2 >> A.^2 ans =

1 1 0 4 >> 1./A

Warning: Divide by zero. ans =

1.0000 -1.0000 Inf 0.5000 3.矩阵运算

矩阵是一个二维数组,所以矩阵的加、减、数乘等运算与数组运算是一致的。但是有两点需要注意:

24

(1)对于乘法、乘方和除法等三种运算,矩阵运算与数组运算的运算符及含义不同:矩阵运算按线性变换定义,使用通常符号;数组运算按对应元素运算定义,使用点运算符;

(2)数与矩阵加减、矩阵除法在数学上是没有意义的,在MATLAB中为简便起见,定义了这两类运算,其含义见表3.

表3 矩阵运算符 运算 符号 说明 转置 A' 加与减 A+B与A-B 同数组运算 数乘矩阵 k*A或A*k 同数组运算 矩阵乘法 A*B 矩阵乘方 A^k 数与矩阵加减 k+A与k-A k+A等价于k*ones(size(A))+A 矩阵除法 左除A\\B,右除B/A 它们分别为矩阵方程AX=B和XA=B的解 >> A=[1 2;3 4];B=[4 3;2 1]; >> 100+A ans =

101 102 103 104 >> A*B,A.*B %注意矩阵运算和数组运算的区别 ans =

8 5 20 13 ans =

4 6 6 4

>> A\\B,B/A,A.\\B,B./A %注意矩阵运算和数组运算的区别 ans =

-6.0000 -5.0000 5.0000 4.0000 ans =

-3.5000 2.5000 -2.5000 1.5000 ans =

4.0000 1.5000 0.6667 0.2500 ans =

4.0000 1.5000 0.6667 0.2500 4.数学函数

数组的数学函数也是按每个元素的运算,使用通常的函数符号,常用数学函数见表4

25

函数 意义 意义 sin 正弦 向0取整 cos 余弦 模余 tan 正切 除法余数 cot 余切 绝对值(模) asin 反正弦 指数函数 acos 反余弦 自然对数 sqrt 开方 以10为底的对数 >> A=[4 -1;3 2]; >> B=exp(A) B =

54.5982 0.3679 20.0855 7.3891 >> C=fix(B) C =

54 0 20 7 >> D=sin(C) D =

-0.5588 0 0.9129 0.6570 >> E=log(D)

Warning: Log of zero. E =

-0.5820 + 3.1416i -Inf -0.0911 -0.4201 5.关系与逻辑运算

MATLAB的关系运算和逻辑运算符都是对于元素的操作,其结果是特殊的逻辑数组(logical array)表5,“真”用1表示,“假”用0表示,而逻辑运算中,所有非零元素作为1(真)处理。

表5 关系运算和逻辑运算 运算符 含义 运算符 含义 < 小于 & 与 <= 小于等于 | 或 > 大于 ~ 非 >= 大于等于 all = = 等于 any ~= 不等于 >> A=-2:4,B=4:-1:-2 A =

-2 -1 0 1 2 3 4 B =

4 3 2 1 0 -1 -2

26

表4 数学函数 函数 fix mod rem abs exp log log10 >> A>B ans =

0 0 0 0 1 1 1 >> A==B ans =

0 0 0 1 0 0 0

>> A&B %逻辑运算中,所有非零元素作为1(真)处理 ans =

1 1 0 1 0 1 1 >> A|B ans =

1 1 1 1 1 1 1 >> find(abs(A)>=2) %返回绝对值大于或等于2的元素的下标 ans =

1 5 6 7 >> any(abs(A)>5) %若存在绝对值大于5的元素,返回1 ans = 0

>> all(abs(A)>5) %若所有元素绝对值大于5,返回1 ans = 0

四、字符串、元胞和结构

除数值(double)以外,常用的数据类型还有字符(char)、元胞(cell)和结构(structure),由此进一步组成字符数组(char array)、元胞数组(cell array)和结构数组(structure array).

1.字符串

MATLAB字符串用单引号对来标识,其数据类型为字符数组。 >> a1='Hello everyone' a1 =

Hello everyone >> a2='各位好' %注意单引号,不是双引号 a2 = 各位好

>> a=[a1,'.',a2,'.'] %字符串拼接 a =

Hello everyone.各位好. >> size(a) ans =

1 19 %共19个字符 >> double(a) ans =

27

Columns 1 through 8

72 101 108 108 111 32 101 118 Columns 9 through 16

101 114 121 111 110 101 46 47351 Columns 17 through 19

52923 47811 46 %中文ASCII码很大 >> a=char(ans) a =

Hello everyone.各位好.

数字字符串与数值之间也可以用num2str和str2num转换。一个数组的元素要么都是数值,要么都是字符串,数值转换后可以与字符串出现在同一数组中。

>> a=12;b=sqrt(a);

>> ['a=',num2str(a),',a的开方=',num2str(b)] ans =

a=12,a的开方=3.4641

MATLAB命令可以定义成一个字符串,使用eval可以使该字符串所表达的MATLAB命令得到执行。

>> fun='x.^2.*sin(x)'; >> x=1;eval(fun) ans =

0.8415

>> x=1:3;eval(fun) ans =

0.8415 3.6372 1.2701 2.元胞和结构

不管是数值数组还是字符数组,其数据结构必须是整齐的。首先数值和字符不能混合,其次小数组拼接成大数组时,其尺寸(size)必须相符(agree)。例如:

>> A=['first';'second'] %错误 ??? Error using ==> vertcat

All rows in the bracketed expression must have the same number of columns.

将不同类型、不同尺寸的数组,加大括号({}),可构成一个元胞。几个元胞可以构成元胞数组。

>> Ac1={'first';1:3};Ac2={'second';[1 2;3 4]}; >> Ac=[Ac1,Ac2] Ac =

'first' 'second' [1x3 double] [2x2 double] >> size(Ac) ans =

2 2 >> Ac(2,1) %小括号,查询Ac的第二行第一列元素 ans =

[1x3 double]

28

>> Ac{2,1} %大括号,查询Ac的第二行第一列元素的具体内容

ans =

1 2 3

一个结构通过“域”来定义,比元胞更丰富、更灵活。几个结构可以合成一个结构数组,但这些结构的域名必须一致。

>> As1.f1='first';As1.f2='second';As2.f1=1:3;As2.f2=[1 2;3 4]; >> As=[As1;As2] As =

2x1 struct array with fields: f1 f2 >> size(As) %注意其size结果与元胞数组不同 ans =

2 1 >> As(2,1).f1 ans =

1 2 3 >> As.f1 ans = first ans =

1 2 3

元胞数组与结构数组可以用struct2cell和cell2struct函数进行适当的转换。 >> Bc=struct2cell(As) %注意结果与Ac的不同 Bc =

'first' [1x3 double] 'second' [2x2 double]

>> Bs=cell2struct(Ac,{'one','two'},1) %定义域名,并指定取域名的维 Bs =

2x1 struct array with fields: one two

看一看Workspace有哪些类型,并观察其字节数。

五、程序设计

1.控制流

前面我们用的命令都是顺序结构的,对于复杂的计算,需要循环和分支等复杂的程序结构。MATLAB控制流语法都以end结尾。常用控制流见表6.

29

类型 循环语句 循环语句 分支语句 分支语句 中断语句 中断语句 中断语句 中断语句 表6 MATLAB常用控制流 语法 解释 for 循环变量=数对于循环变量依次取数组中的值,循环执行命令组, 组直到循环变量遍历数组。数组最常用的形式是命令组 初值:增量:终值 end while 条件式, 当条件式满足,循环执行命令组直到条件不满 命令组 足。使用While语句要注意避免出现死循环 end if 条件1, 如果条件式1满足,则执行命令组1,且结束该 命令组1; 语句;否则检查条件式2,若满足则执行命令组elesif 条件2, 2,且结束该语句;„„;若所有条件式都不满 命令组2; 足,则执行命令组k,并结束该语句。最常用的„„; 格式是 else, if 条件式 命令组k; 命令组 end end switch 分支变量 若分支变量的取值1,则执行命令组1,且结束case 值1, 该语句;若分支变量的取值2,则执行命令组2, 命令组1; 且结束该语句;„„若分支变量不取所列出的case 值2, 值,则执行命令组k 命令组2; „„; otherwise 命令组k; end pause 暂停执行,直到击键盘,pause(n)为暂停n秒后再继续 break 中断执行,用在循环语句内表示跳出循环 return 中断执行该程序,回到主调函数或命令窗口 error(字符串) 提示错误并显示字符说明 100例1 计算s1 2nn1>> clear;s=0;

>> for n=1:100 s=s+1/n/n; end >> s s =

1.6350 2.M脚本文件

复杂程序结构在命令窗口调试保存都不方便,所以进行复杂的运算大都使用程序文件。从命令窗口用命令“edit”就进入MATLAB的程序编辑器窗口,用

30

以编写用户的M文件。M文件可分为两类:M脚本文件和M函数文件。

将多条MATLAB语句写在编辑器中,并以.m文件保存在适当的目录中(这个目录须为MATLAB的搜索目录),就得到一个M脚本。如我们将例1中的几条语句写在编辑器中,

保存为naega_1,然后在命令窗口执行: >>naega_1 s=

1.6350

执行M脚本文件也可以在程序编辑器的Debug菜单选Run。使用编辑器也可打开和修改M文件、观察变量值、调试程序等。

注意:M文件名一律以字母开头,以字母、数字或下划线组成,不要含有空格、减号等,并要防止它与系统的变量名、系统内部的M函数名冲突。例如1.m,ega-1.m,ega.1.m都是不合法的。另外,别忘了每次修改程序后都要存盘。

3.M函数文件

M脚本文件没有参数传递功能,当我们需要修改程序中的某些变量值,必须修改M文件。而M函数文件使得我们可以进行参数传递。

M函数文件以function开头,格式如下: function 输出变量=函数名(输入变量) 语句 例如,写函数文件: % M函数naega_1f.m function s=f(m) s=0;

for n=1:m s=s+1/n/n; end

保存为naega_1f.m,在命令窗口执行: >>clear;naega_1f(100),naega_1f(1000) ans=

1.6350 ans=

1.6439

注意:在MATLAB中,使用M函数是以该函数的磁盘文件主名调用,而不是文件中的函数名,但为了增强程序可读性,最好两者同名。

M函数不能像M脚本那样在编辑器窗口用Debug\\run执行,因为M函数必须给予输入参数值。M函数常常被M脚本或其他M函数调用。

4.函数句柄和内嵌函数

M函数除了直接用其函数名调用之外,也可以作为一个参数调用。调用时使用所谓函数句柄(handle)方式。MATLAB命令feval用于执行函数的参数方式。例如:

>>fname=@naega_1f;feval(fname,1000) ans=

1.6439

比较简单的函数表达式可以不用写成外部M函数,而是用更简捷的内嵌

31

(inline)函数方式。inline的使用格式如下: fun=inline(expr,arg1,agr2,…) expr为函数表达式字符串,arg1,arg2,….为自变量名字符串 例如: >> fname=inline('sum(1./(1:n).^2)','n') fname =

Inline function:

fname(n) = sum(1./(1:n).^2) >> feval(fname,1000) ans =

1.6439 5.其他 (1)注释

为了增强程序的可读性,程序中常常需要注释语句。M文件开头一般应有一段注释。注释用%开头,顶格书写,对本行后面字符起作用,说明文件的功能和使用方法。注释语句不参与运算,只起说明作用。使用Help可看到。注释符(%)也常用于程序调试。

(2)对话

Input在交互式执行程序中用于提示键盘输入,Disp用于屏幕显示。 例2 编写一个脚本文件,使对键盘提示输入的向量求得元素总和。 % M文件naega_2.m

%用途:本程序提示输入一个向量,并求得元素总和 %用法:输入向量用中括号,元素之间用逗号 clear A;

A=input('Enter a vector:'); d=sum(A);

disp(['The sum is',num2str(d)]); 然后在命令窗口执行: >>naega_2

Enter a vector:[1 2 3 4] The sum is 10 >>help naega_2

用途:本程序提示输入一个向量,并求得元素总和 用法:输入向量用中括号,元素之间用逗号 (4)子函数

M函数中允许使用子函数。M函数中第一个function为主函数,其他function为子函数。子函数只能被同一文件的主函数和其他子函数调用,不能被外部函数调用。

(5)全程变量与局部变量

M函数中所有变量为局部变量,而脚本文件中所有变量同命令窗口的命令一样都是全程变量。M函数变量值传递主要通过其输入输出变量,但也可以用global定义全程变量。它的意义与普通全程变量稍有区别,只对有定义的文件起作用。

(6)nargin与nargout

32

在M函数内,nargin表示该函数的输入变量个数,nargout表示该函数的输出变量个数。

(7)提高速度

MATLAB软件主要缺点是执行循环语句时速度慢。好的M程序文件应尽量使用数组运算和内部函数,少用循环语句,以提高运算速度。尽管MATLAB数组无须定义尺寸,但经常改变数组尺寸会影响速度,采取一些预分配方法可提高运算速度。另外,减少运行过程中不必要的结果显示也可提高速度。

(8)强行中断

使用快捷键Ctrl+C可以强行中断程序运行。

例3 编写一个M函数,对于任意输入的向量x,可以计算下列分段函数值构成的向量:

x2,x1f(x)1,1x1

32x,x1%M函数naega_3a.m function y=naega_3a(x) n=length(x); for i=1:n if x(i)>1 y(i)=x(i)^2; elseif x(i)>-1 y(i)=1; else y(i)=3+2*x(i); end end %M函数naega_3b.m function y=naega_3b(x) y=zeros(size(x)); k1=find(x>1);y(k1)=x(k1).^2; k2=find(x>-1&x<=1);y(k2)=1; k3=find(x<=-1);y(k3)=3+2*x(k3); 程序naega_3a.m是按通常思路所编程序,可读性好,但速度慢。程序naeag_3b.m是根据MATLAB数组运算设计的程序,程序简短且速度快,数组维数越高,效率改进越显著。下列程序中tic表示计时开始,toc为读数。

>>clear;tic;x=-2:0.0005:2;y=naega_3a(x);toc %时间因电脑不同而有区别 elapsed_time = 1.0400

>> clear;tic;x=-2:0.0005:2;y=naega_3b(x);toc elapsed_time = 0.4400

33

六、作图

表7 常用作图命令和函数 主题词 含 义 主题词 含 义 plot 基本二维图形 clabel 等高线高度标志 fplot 一元函数图像 grid 格栅 ezplot 画二维曲线的符号命令 hold 图形保持 plot3 空间曲线 axis 定制坐标轴 meshgrid 网格数据生成 view 改变视点 mesh 网面图 subplot 子图 surf 曲面图 figure 新图形窗口 contour 等高线图 clf 清除图形 contour3 三维等高线图 close 关闭图形窗口 title 标题 ylabel y轴说明 xlabel x轴说明 zlabel z轴说明 1.曲线图 plot(x,y) 作出以数据(x(i),y(i))为节点的折线图,其中x,y为同长度的向量 fplot('fun',[a,b]) 作出函数fun在区间[a,b]上的函数图,fun可以是M函数主名,也可是字符串 ezplot(fun,lims) 绘制字符串fun(可是显函数、隐函数或参数方程)指定的函数 plot3(x,y,z) 空间曲线图,其中x,y,z为同长度向量 >> plot([1 4 2 5],[1 3 -1 2]) 依次将(1,1),(4,3),(2,-1),(5,2)连接起来得一条折线(图2) 图2

图形显示在图形窗口。在图形窗口可以使用File菜单保存(Save)为M文件,导出(Export)为图形文件。也可利用图形窗口Edit菜单Copy figure作为图片复制到剪贴板,从而进一步粘贴到Word或其他应用程序中。

图形的线型、标记、颜色均可根据要求设定。常用的见表8

34

b g r m y k 颜色 蓝(默认) 绿 红 洋红 青 黑

表8 图形元素设定 线型 - 实线(默认) : 虚线 -. 点划线 -- 划线 0.2 . o x + * 标记 无标记(默认) 点 圈 叉 十字 星 例4 两个一元函数yx3x1和yxsin(5x)在区间-1% M文件naega_4.m

fplot('x^3-x-1',[-1,2]);

hold on; %在作下一幅图时保留已有图像 x=-1:0.2:2;

y=abs(x).^2.*sin(5*x); %注意数组运算 plot(x,y,':ro');

hold off; %释放hold on

2.曲面图 [x,y]=meshgrid(xa,ya) 当xa,ya分别为m维和n维行向量,得到x和y均为n行m列矩阵。meshgrid常用于生成X-Y平面上的网格数据 mesh(x,y,z) 绘制网面图,是最基本的曲面图形命令,其中x,y,z是同阶矩阵,表示曲面三维数据 surf(x,y,z) 绘制曲面图,与mesh用法类似 contour(x,y,z) 绘制等高线图,与mesh用法类似 contour3(x,y,z) 绘制三维等高线图,与mesh用法类似 例如: >> xa=6:8;ya=1:4; %生成x,y各自的节点 >> [x,y]=meshgrid(xa,ya); %生成X-Y面上网格

35

>> z=x.^2+y.^2; %计算X-Y面上各网格点的z轴高度 >> mesh(x,y,z) >> [x,y,z] ans =

6 7 8 1 1 1 37 50 65 6 7 8 2 2 2 40 53 68 6 7 8 3 3 3 45 58 73 6 7 8 4 4 4 52 65 80 这3组数据构成空间网面的12格点坐标

例5 二元函数图zxexp(x2y2) %M文件 naega_5.m

clear;close;

% close关闭当前图形窗口

xa=-2:0.2:2;ya=xa; [x,y]=meshgrid(xa,ya); z=x.*exp(-x.^2-y.^2); mesh(x,y,z);pause %网格图,pause暂停直到按任意键 surf(x,y,z);pause %网面图 contour(x,y,z);pause %等高线图

contour(x,y,z,[0.1 0.1]); %z=0.1的一条等高线

36

3.图形说明和定制 title('字符串') 图形标题说明 xlabel,ylabel,zlabel 用法类似于title,分别说明坐标轴x,y,z grid on/off 显示/不显示格栅 box on/off 使用/不使用坐标框 hold on/off 保留/释放现有图形 axis off/on 不显示/显示坐标轴 axis([a,b,c,d]) 定制二维坐标轴范围 a%M文件naega_6.m clear;close;

t=0:0.1:20;r=exp(-0.2*t);th=0.5*pi*t; x=r.*cos(th);y=r.*sin(th);z=sqrt(t); subplot(1,2,1) plot3(x,y,z); title('螺旋线'); subplot(1,2,2) plot3(x,y,z);

axis([-1 1 -1 1 0 4]) grid on;

37

七、在线帮助和文件管理

1.在线帮助 help 显示MATLAB主题目录 help 子目录名 显示子目录中所有MATLAB系统命令及函数 help 命令或函数 显示该命令或函数的说明部分 lookfor 关键字 显示与该关键字有关的命令和函数 type M文件名 显示M文件程序代码 which M文件名 显示指定的MATLAB文件的路径 demo 演示MATLAB功能 MATLAB提供了两种形式的帮助系统:纯文本帮助、HTML帮助。纯文本帮助是最常用、最经济的在线帮助,不需要另外安装,是边用边学计算机软件的最有效的方法。

>>help %显示MATLAB及其工具箱的主题目录,其中有graph3d >>help graph3d %显示3维图形主题目录内所有M命令和函数,其中有mesh

>>help mesh %显示M函数mesh的用法说明,即其M文件的注释部分

>>which mesh %显示M函数mesh所在的目录 >>type mesh %显示函数mesh的M文件程序代码

>>lookfor surface %显示MATLAB搜索路径中凡是第一行注释含surface的

M命令和函数,其中有函数mesh

若我们现在要解决一个线性规划问题,但不知道怎么用MATLAB求解,可

利用线性规划的关键字programming:

>>lookfor programming

可以找到有关programming的很多命令,其中有一个LINPROG是线性规划

(Linear programming).再用:

>>help linprog

可以得到使用linprog解线性规划问题用法的详细说明。进一步使用: >>type linprog

可以看到linprog的M文件程序代码。

2.文件和目录管理

MATLAB文件有M、Mat、Mex等。其中M文件是最重要的,MATLAB绝大多数内部命令和函数是M文件,用户自编的程序一般也是M文件。MATLAB只执行当前目录(Current Directory)和搜索路径(写在文件pathdef.m)中的M文件。

初学者在M文件的保存上经常出现下列几种错误: (1)文件修改后没有保存;

(2)文件保存的目录不在当前目录和MATLAB搜索路径中; (3)文件名使用了常数或内存中的变量,如1.m,pi.m等;

(4)文件名用了减号、空格等非法字符,如eg2-1.m,eg2.1.m等;

(5)文件名与MATLAB内建(build-in)函数和其他内部函数冲突,如

38

mesh.m,fitfun.m等。

39

因篇幅问题不能全部显示,请点此查看更多更全内容