Matlab上机学习指导.docx

1、Matlab上机学习指导五 程序设计MATLAB作为一种高级语言，它不仅可以如前几节所介绍的那样，以一种人机交互式的命令行的方式工作，还可以像BASIC、FORTRAN、C等其他高级计算机语言一样进行控制流的程序设计，即编制一种以 .m为扩展名的MATLAB程序,简称M文件,所谓M文件就是由MATLAB语言编写的可在MATLAB语言环境下运行程序源代码文件。由于商用的MATLAB软件是用C语言编写而成。因此，M文件的语法与C语言十分相似。对广大参加建模竞赛且学过C语言的同学来说，M文件的编写是相当容易的。8.4.1 M文件可以分为脚本文件(Script)和函数文件(Function)两种。1

2、脚本文件：（1）编写文档:点击MATLAB指令窗口上面最左端的图标，即新建文件，就可打开MATLAB文件编辑器，像word一样。用户即可在空白窗口中编写程序。例如输入下面的程序： x=linspace(0,2*pi,20); y=sin(x); plot(x,y,r+) title(正弦曲线) (2）点击文件编辑器上面工具条中的保存，命名（例如将上面的程序命名为picture），然后保存。这样的文件就是M-文件：picture.m (3)运行：i）在命令窗口中输入文件名(如上面的picture),然后执行。ii）或直接在文件编辑器上面的工具条中找到debug里面的run（即运行），点击即可（或

3、直接按F5）。iii）如果发现错误，在写好的程序里直接改正，然后再保存，再运行，注意一定要先保存，后运行。2 函数文件：我们经常用到的像sin、cos、exp这样的一些函数都是MATLAB软件自身所带的函数，因此直接应用即可，但有时我们为了解决一些问题需要自己编写函数。自己编写函数有两个基本要求i）必须在MATLAB文件编辑器中编写，也是M-文件。ii）函数名和文件名必须相同。 例1: 编写函数 , 计算 （1）打开MATLAB文件编辑器，即点图标，输入 function y= fun1(x) % 表示y是x的函数，x是自变量，fun1是我们自己命名的函数名y=(x3 - 2*x2 + x -

4、 6.3)/(x2 + 0.05*x - 3.14);然后保存。注：在自己编写的函数前都要写上function（关键字），表示这是自己定义的函数。fun1表示函数名，那么最后保存文件也要命名为fun1。 （2）这样在命令窗口中就可以像应用sin、cos那样来使用函数fun1，如：在命令窗口中输入 fun1(1)*fun1(2)+fun1(3)*fun1(3) ans = -12.6023例2：也可以多输入、多输出，如输入向量x，y，要求编写函数输出x、y各自的长度、x与y的数量积、向量积。(1)建立M文件如下function L1,L2,nei,wai=fun2(x,y); % 定义函数名为f

5、un2L1=sqrt(sum(x.2), % x的长度L2=sqrt(sum(y.2), % y的长度nei=dot(x,y), % x与y的数量积wai=cross(x,y), % x与y的向量积（2）保存命名为fun2.m（3）在命令窗口中输入如下语句 x=4 3 1;y=-2 1 6;fun2(x,y)运行结果为：L1 =5.0990L2 =6.4031nei = 1wai =17 -26 108.4.2 程序结构和任何计算机语言一样，matlab有三种程序结构：顺序结构、选择结构（或分支结构）、循环结构1顺序结构，即按照程序语句的顺序逐条运行命令，前面函数作图都是顺序结构，这里不再叙述

6、2选择结构例3：if-end语句，例：cost=10;number=12; % cost和number都是变量名if number8 % 如果number8执行后面的语句，否则跳出sums=number*0.95*cost;end,sums 例4：if-else-end语句，例：cost=10;number=5; % 改变number的初值，看结果有何不同if number8 sums=number*0.95*cost;else sums=number*0.5*cost;end,sums例5：定义函数，计算， 建立函数文件fff.mfunction y=fff(x)if x=2 & x8 &

7、x y=fff(0.1),fff(1),fff(9),fff(22),fff(2*pi)结果y =1.1000 2.0000 31.0000 -1.0088 18.84963循环结构for循环语句（这里的for语句与C语言中的for语句不同，要更简单一些）例6：一个简单的for循环示例for i=1:10; % i依次取1,2,10,.x(i)=2*i; % 对每个i值，重复执行该指令end; % 表示循环结束，每一个for要对应一个endx % 要求显示运行后数组x的值。输入后观察结果，体会for语句的作用。注：在MATLAB里（在C语言中也一样），“”的作用表示把等号右边的值送给左边的变量

8、，这和数学中相等的含义不同。while循环语句例7: Fibonacci 数列：1，1，2，3，5，8， 即：，（ 1，2，3）现要求该数列中第一个大于10000 的元素。a(1)=1;a(2)=1;i=2;while a(i)=10000 % 当a(i)=10000a(i),break; % 表示跳出循环end;end, i 9 练习：1 用循环语句求1到100的整数之和。2 定义函数然后计算3 （提示：建立两个m文件，一个用来定义函数f(x)，一个用来写程序，用while语句，当时，执行算法，当时，停止运行，输出结果）六 插值与拟合作为Matlab的简单应用，本节我们讲插值与拟合。插值与拟

9、合是来源于实际、又广泛应用于实际的两种重要方法。随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。1 一维插值在实际中遇到的函数，有许多我们并不知道其表达式，而只知道它在某些点上的函数值，如，但其它点上的值就不知道了。现在为了获得这些其它点上的函数值，我们设法构造一个比较简单的函数（例如多项式）来近似代替，并要求经过那些已知点，即与在这些点上的函数值相同。这样我们想要知道在某点的函数值就可以通过计算在该点的值来近似。这种方法就叫插值法，最简单的方法就是把各相邻已知点用线段连起来，叫做线性插值，工程上常用的还有二次插值、三次插值、三次样条插值等等，

10、具体方法可参看任何一本数值方法的教材，使用这些插值方法不需编制函数程序，MATLAB自身提供了内部函数interp1供我们直接调用，格式如下：yi=interp1(x,y,xi,method) 对一组点(x,y) 进行插值，计算插值点xi的函数值。x为节点向量值，y为对应的节点函数值。 method用来指定插值的算法。默认为线性算法。其值常用的可以是如下的字符串： nearest 线性最近项插值。 linear 线性插值。 spline 三次样条插值。 cubic 三次插值。所有的插值方法要求x是单调的（从小到大），但不要求等距。例1：考虑下列问题，12小时内，一小时测量一次室外温度。数据如下

11、：时间：1，2，3， 4， 5， 6， 7， 8， 9，10，11，12温度：5，8，9，15，25，29，31，30，22，25，27，24现在根据以上数据估计3.2，4.7等时刻的温度hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24;t=interp1(hours, temps, 3.2,4.7) % 一阶线性插值，如果只估计一个点的值，则无须加方括号运行结果t =10.2000 22.0000输入如下程序，画出插值曲线hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:1

12、2; %每隔0.1一个点，共111个点t=interp1(hours, temps, h) ; %给出这111个点的插值的结果plot(hours, temps, + , h, t)工程上常用所谓三次样条插值，效果比较好，只需在上面插值命令后加上spline即可t=interp1(hours, temps, 3.2，4.7, spline);t=interp1(hours, temps, h，spline) ;图像如下Matlab也能够完成二维插值运算，函数为interp2，用法与interp1基本相同，只是输入和输出参数均为矩阵，对应于平面上的数值点。2 曲线拟合在科学实验的统计方法研究中,

13、往往要从一组实验数据中寻找出自变量x 和因变量y之间的函数关系。由于观测数据往往不够准确，因此并不要求经过所有的点（前述插值法要求插值函数必须经过所有已知点）,而只要求在给定点上误差按照某种标准达到最小，通常采用欧氏范数作为误差量度的标准。这就是所谓的最小二乘法。我们找到的这条曲线叫做所给数据点的拟合曲线，一般先在平面直角坐标系中描出所给数据点的图形，称为散点图，其次观察散点图，根据经验确定一条曲线，尽可能好的反映x与y的变化关系（通常用多项式），设出其一般形式（如含待定系数的多项式），然后根据已知数据，具体求出这条曲线。在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行（多项式

14、拟合）。例2：两组数据如下：x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1;y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;先做出x与y的散点图（用命令plot(x,y,o)）估计变量x与y之间大概是什么样的函数关系，比如估计是一次函数，则用下面的命令：n=1;p1=polyfit(x,y,n) % n表示用n阶多项式拟合，n=1为线性拟合，输出p1为所求多项式的系数poly2sym(p1) % 前面的拟合命令只给出多项式的系数，用此命令则将前面的p1转化为我们熟悉的多项式vpa(poly2sym(p1),

15、5) % 上面命令给出的结果是有理数形式，此命令将结果转化为数值形式，5表示显示5位有效数字x1=0:0.01:1; % 由此以后三句是画出拟合曲线的图像y1=polyval(p1,x1); % 此句是在x1这些点处求出多项式的值，送给y1plot(x,y,o,x1,y1)运行结果p1 =10.3185 1.4400ans =5808773505743561/562949953421312*x+31679/22000ans =10.318*x+1.4400n=2（用二次多项式拟合）时的图形看来要比用一次函数拟合效果要好一些n=10时当n比较大时，如上面n=10，数据点之间出现大的波动（虽然曲线

16、经过所有点）。当企图进行高阶曲线拟合时，这种波动现象经常发生，并不利于我们认识两组数据之间的规律，因此并不是阶数越高越好，实际问题当中，适当选择一个即可。例3：非线性拟合，这里介绍函数lsqcurvefit来拟合，拟合下面两组数据，t为自变量t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.011)作散点图2)观察此图，我们想用指数函数来拟合，建立函数文件fun01.m用来计算函数值function c=fun01(x,t)c=x(1)*exp(x(2)*t);3)建立脚本文件，输入程序x0=10, 0.5;t=0.25 0.5

17、1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;x,norm,res=lsqcurvefit(fun01,x0,t,c)其中x是所求系数，norm是误差的平方和，res是误差向量，x0是迭代初始值运行结果x =20.2413 -0.2420 %即r=20.2413，k=-0.2420norm =1.0659res =-0.1568 -0.2152 0.5311 -0.0198 -0.4143 0.4744 0.2394 -0.5006 -0.08894)画图：tt=0:0.2:8;cc=x(1).*exp(x(2).*tt);plot(tt,cc,t,c,o)