五 程序设计
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)