Matlab上机学习指导.docx
《Matlab上机学习指导.docx》由会员分享，可在线阅读，更多相关《Matlab上机学习指导.docx（13页珍藏版）》请在冰豆网上搜索。
Matlab上机学习指导
五程序设计
MATLAB作为一种高级语言，它不仅可以如前几节所介绍的那样，以一种人机交互式的命令行的方式工作，还可以像BASIC、FORTRAN、C等其他高级计算机语言一样进行控制流的程序设计，即编制一种以.m为扩展名的MATLAB程序,简称M文件,所谓M文件就是由MATLAB语言编写的可在MATLAB语言环境下运行程序源代码文件。
由于商用的MATLAB软件是用C语言编写而成。
因此，M文件的语法与C语言十分相似。
对广大参加建模竞赛且学过C语言的同学来说，M文件的编写是相当容易的。
8.4.1M文件可以分为脚本文件（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）函数名和文件名必须相同。
例1:
编写函数
计算
（1）打开MATLAB文件编辑器，即点图标
，输入
functiony=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=[431];
y=[-216];
fun2（x,y）
运行结果为：
L1=
5.0990
L2=
6.4031
nei=
1
wai=
17-2610
8.4.2程序结构
和任何计算机语言一样，matlab有三种程序结构：
顺序结构、选择结构（或分支结构）、循环结构
1顺序结构，即按照程序语句的顺序逐条运行命令，前面函数作图都是顺序结构，这里不再叙述
2选择结构
例3：
if-end语句，例：
cost=10;number=12;%cost和number都是变量名
ifnumber>8%如果number>8执行后面的语句，否则跳出
sums=number*0.95*cost;
end,
sums
例4：
if-else-end语句，例：
cost=10;number=5;%改变number的初值，看结果有何不同
ifnumber>8
sums=number*0.95*cost;
elsesums=number*0.5*cost;
end,
sums
例5：
定义函数
，
计算
，
，
，
，
建立函数文件fff.m
functiony=fff（x）
ifx<2
y=x+1;
elseifx>=2&x<=8
y=3*x;
elseifx>8&x<=20
y=4*x-5;
elsey=cos（x）+sin（x）;
end
在命令窗口中计算
>>y=[fff（0.1）,fff
（1）,fff（9）,fff（22）,fff（2*pi）]
结果y=
1.10002.000031.0000-1.008818.8496
3循环结构
for循环语句（这里的for语句与C语言中的for语句不同，要更简单一些）
例6：
一个简单的for循环示例
fori=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，…即：
，（
1，2，3…）现要求该数列中第一个大于10000的元素。
a
（1）=1;a
（2）=1;i=2;
whilea（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的矩阵
fori=3:
n
a（i）=a（i-1）+a（i-2）;
ifa（i）>=10000
a（i）,
break;%表示跳出循环
end;
end,i
9练习：
1用循环语句求1到100的整数之和。
2定义函数
然后计算
3
（提示：
建立两个m文件，一个用来定义函数f（x），一个用来写程序，用while语句，当
时，执行算法，当
时，停止运行，输出结果）
六插值与拟合
作为Matlab的简单应用，本节我们讲插值与拟合。
插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。
随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。
1一维插值
在实际中遇到的函数，有许多我们并不知道其表达式
，而只知道它在某些点上的函数值，如
，但其它点上的值就不知道了。
现在为了获得这些其它点上的函数值，我们设法构造一个比较简单的函数
（例如多项式）来近似代替
，并要求
经过那些已知点
，即
与
在这些点上的函数值相同。
这样我们想要知道
在某点的函数值就可以通过计算
在该点的值来近似。
这种方法就叫插值法，最简单的方法就是把各相邻已知点用线段连起来，叫做线性插值，工程上常用的还有二次插值、三次插值、三次样条插值等等，具体方法可参看任何一本数值方法的教材，使用这些插值方法不需编制函数程序，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=[589152529313022252724];
t=interp1（hours,temps,[3.2,4.7]）%一阶线性插值，如果只估计一个点的值，则无须加方括号
运行结果
t=
10.200022.0000
输入如下程序，画出插值曲线
hours=1:
12;
temps=[589152529313022252724];
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曲线拟合
在科学实验的统计方法研究中,往往要从一组实验数据
中寻找出自变量x和因变量y之间的函数关系
。
由于观测数据往往不够准确，因此并不要求
经过所有的点
（前述插值法要求插值函数
必须经过所有已知点）,而只要求在给定点
上误差
按照某种标准达到最小，通常采用欧氏范数
作为误差量度的标准。
这就是所谓的最小二乘法。
我们找到的这条曲线
叫做所给数据点
的拟合曲线，一般先在平面直角坐标系中描出所给数据点
的图形，称为散点图，其次观察散点图，根据经验确定一条曲线，尽可能好的反映x与y的变化关系（通常用多项式），设出其一般形式（如含待定系数的多项式），然后根据已知数据，具体求出这条曲线。
在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行（多项式拟合）。
例2：
两组数据如下：
x=[0.1.2.3.4.5.6.7.8.91];
y=[-.4471.9783.286.167.087.347.669.569.489.3011.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.31851.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）观察此图，我们想用指数函数来拟合
，
建立函数文件fun01.m用来计算函数值
functionc=fun01（x,t）
c=x
（1）*exp（x
（2）*t）;
3）建立脚本文件，输入程序
x0=[10,0.5];
t=[0.250.511.523468];
c=[19.2118.1515.3614.1012.899.327.455.243.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.21520.5311-0.0198-0.41430.47440.2394-0.5006-0.0889
4）画图：
tt=0:
0.2:
8;
cc=x
（1）.*exp（x
（2）.*tt）;
plot（tt,cc,t,c,'o'）
升级会员 