Matlab上机学习指导.docx
《Matlab上机学习指导.docx》由会员分享,可在线阅读,更多相关《Matlab上机学习指导.docx(13页珍藏版)》请在冰豆网上搜索。
![Matlab上机学习指导.docx](https://file1.bdocx.com/fileroot1/2023-1/22/87594627-f95e-454f-bb19-fe2999d53c83/87594627-f95e-454f-bb19-fe2999d53c831.gif)
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')