Matlab上机学习指导

五 程序设计

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

8.4.1 M文件可以分为脚本文件(Script)和函数文件(Function)两种。

1 脚本文件：（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（即运行），

点击即可（或直接按F5）。

iii）如果发现错误，在写好的程序里直接改正，然后再保存，再运行，注意

一定要先保存，后运行。

2 函数文件：我们经常用到的像sin、cos、exp这样的一些函数都是MATLAB软件自身所带的函数，因此直接应用即可，但有时我们为了解决一些问题需要自己编写函数。自己编写函数有两个基本要求

i）必须在MATLAB文件编辑器中编写，也是M-文件。 ii）函数名和文件名必须相同。

32例1: 编写函数 f?x?2x?x?6.3, 计算f(1)f(2)?f2(3)

2x?0.05x?3.14 （1）打开MATLAB文件编辑器，即点图标，输入

function y= fun1(x) % 表示y是x的函数，x是自变量，fun1是我们

自己命名的函数名

y=(x^3 - 2*x^2 + x - 6.3)/(x^2 + 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); % 定义函数名为fun2 L1=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.0990 L2 =

6.4031 nei = 1 wai =

17 -26 10 8.4.2 程序结构

和任何计算机语言一样，matlab有三种程序结构：顺序结构、选择结构（或分支结构）、循环结构

1顺序结构，即按照程序语句的顺序逐条运行命令，前面函数作图都是顺序结构，这里不再叙述

2选择结构

例3：if-end语句，例：

cost=10;number=12; % cost和number都是变量名

if number>8 % 如果number>8执行后面的语句，否则跳出

sums=number*0.95*cost; end, sums

例4：if-else-end语句，例：

cost=10;number=5; % 改变number的初值，看结果有何不同 if number>8

sums=number*0.95*cost;

else sums=number*0.5*cost; end, sums

x?2?x?1,?2?x?8?3x,例5：定义函数f(x)??，

4x?5,8?x?20???cosx?sinx,x?20计算f(0.1)，f(1)，f(9)，f(22)，f(2?) 建立函数文件fff.m

function y=fff(x) if x<2

y=x+1;

elseif x>=2 & x<=8

y=3*x;

elseif x>8 & x<=20

y=4*x-5;

else y=cos(x)+sin(x); end

在命令窗口中计算

>> y=[fff(0.1),fff(1),fff(9),fff(22),fff(2*pi)] 结果y =

1.1000 2.0000 31.0000 -1.0088 18.8496

3循环结构

for循环语句（这里的for语句与C语言中的for语句不同，要更简单一些） 例6：一个简单的for循环示例

for i=1:10; % i依次取1,2,?10,. x(i)=2*i; % 对每个i值，重复执行该指令

end; % 表示循环结束，每一个for要对应一个end x % 要求显示运行后数组x的值。 输入后观察结果，体会for语句的作用。

注：在MATLAB里（在C语言中也一样），“?”的作用表示把等号右边的值送给左边的变量，这和数学中相等的含义不同。

while循环语句

例7: Fibonacci 数列：1，1，2，3，5，8，? 即：ak?2?ak?ak?1，（ k?1，2，3?）现要求该数列中第一个大于10000 的元素。

a(1)=1;a(2)=1;i=2;

while a(i)<=10000 % 当a(i)<=10000时执行后面的语句，否则跳

出循环

a(i+1)=a(i-1)+a(i); i=i+1;

end;

i,a(i), %显示i和a(i)

例8：用for 循环语句来寻找Fibonacc 数列中第一个大于10000 的元素。

n=100;a=ones(1,n); % a是一个一行，n列的所有元素为1的矩阵 for i=3:n

a(i)=a(i-1)+a(i-2); if a(i)>=10000

a(i),

break; % 表示跳出循环 end;

end, i

9 练习：1 用循环语句求1到100的整数之和。

2 定义函数3

??(x?1),??2f(x)??1?x,?0,??x??1?1?x?1然后计算f(?2)f(0)f(2)1?x

开始输入a b εf(a)→y0(a+b)/2→x,f(x)→yNYy0*y>0?x→b|b-a|<ε输出x结束x→aN

（提示：建立两个m文件，一个用来定义函数f(x)，一个用来写程序，用while语句，当b?a??时，执行算法，当b?a??时，停止运行，输出结果）

六 插值与拟合

作为Matlab的简单应用，本节我们讲插值与拟合。插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。

1 一维插值

在实际中遇到的函数，有许多我们并不知道其表达式f(x)，而只知道它在某些点上的函数值，如y0?f(x0),y1?f(x1),?,yn?f(xn)，但其它点上的值就不知道了。现在为了获得这些其它点上的函数值，我们设法构造一个比较简单的函数?(x)（例如多项式）来近似代替f(x)，并要求?(x)经过那些已知点(x0,y0),(x1,y1),?,(xn,yn)，即?(x)与f(x)在这些点上的函数值相同。这样我们想要知道f(x)在某点的函数值就可以通过计算?(x)在该点的值来近似。这种方法就叫插值法，最简单的方法就是把各相邻已知点用线段连起来，叫做线性插值，工程上常用的还有二次插值、三次插值、三次样条插值等等，具体方法可参看任何一本数值方法的教材，使用这些插值方法不需编制函数程序，MATLAB自身提供了内部函数interp1供我们直接调用，格式如下：

yi=interp1(x,y,xi,’method’)

对一组点(x,y) 进行插值，计算插值点xi的函数值。x为节点向量值，y为对应的节点函数值。

method用来指定插值的算法。默认为线性算法。其值常用的可以是如下的字符串： nearest 线性最近项插值。 linear 线性插值。

spline 三次样条插值。

cubic 三次插值。

所有的插值方法要求x是单调的（从小到大），但不要求等距。

例1：考虑下列问题，12小时内，一小时测量一次室外温度。数据如下： 时间：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:12; %每隔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 曲线拟合

在科学实验的统计方法研究中,往往要从一组实验数据(xi,yi)中寻找出自变量x 和因变量y之间的函数关系y?f(x)。由于观测数据往往不够准确，因此并不要求y?f(x)经过所有的点 (xi,yi)（前述插值法要求插值函数?(x)必须经过所有已知点）,而只要求在给定点xi上误差?i?f(xi)?yi按照某种标准达到最小，通常采用欧氏范数?2作为误差量

度的标准。这就是所谓的最小二乘法。我们找到的这条曲线f(x)叫做所给数据点(xi,yi)的

拟合曲线，一般先在平面直角坐标系中描出所给数据点(xi,yi)(i?1,?,n)的图形，称为散点图，其次观察散点图，根据经验确定一条曲线，尽可能好的反映x与y的变化关系（通常用多项式），设出其一般形式（如含待定系数的多项式），然后根据已知数据，具体求出这条曲线。

在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行（多项式拟合）。 例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),5) % 上面命令给出的结果是有理数形式，此命令将结果转化为数

值形式，5表示显示5位有效数字

x1=0:0.01:1; % 由此以后三句是画出拟合曲线的图像

y1=polyval(p1,x1); % 此句是在x1这些点处求出多项式的值，送给y1 plot(x,y,'o',x1,y1) 运行结果 p1 =

10.3185 1.4400 ans =

5808773505743561/562949953421312*x+31679/22000 ans =

10.318*x+1.4400

n=2（用二次多项式拟合）时的图形

看来要比用一次函数拟合效果要好一些 n=10时

当n比较大时，如上面n=10，数据点之间出现大的波动（虽然曲线经过所有点）。当企图进行高阶曲线拟合时，这种波动现象经常发生，并不利于我们认识两组数据之间的规律，因此并不是阶数越高越好，实际问题当中，适当选择一个即可。

例3：非线性拟合，这里介绍函数lsqcurvefit来拟合，拟合下面两组数据，t为自变量 t 0.25 0.5 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 1)作散点图

2)观察此图，我们想用指数函数来拟合c?re， 建立函数文件fun01.m用来计算函数值 function c=fun01(x,t) c=x(1)*exp(x(2)*t);

3)建立脚本文件，输入程序 x0=[10, 0.5];

t=[0.25 0.5 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.2420 norm = 1.0659 res =

-0.1568 -0.2152 0.5311 -0.0198 -0.4143 0.4744 0.2394 -0.5006 -0.0889 4)画图： tt=0:0.2:8;

cc=x(1).*exp(x(2).*tt); plot(tt,cc,t,c,'o')

?kt

本文来源：https://www.bwwdw.com/article/4xwo.html