数模中常用的matlab语言和方法.docx
《数模中常用的matlab语言和方法.docx》由会员分享，可在线阅读，更多相关《数模中常用的matlab语言和方法.docx（23页珍藏版）》请在冰豆网上搜索。
数模中常用的matlab语言和方法
数模培训
一、矩阵的输入
二、循环语句
三、M文件及M函数
四、优化问题
五、回归分析
一矩阵的输入
常用的矩阵输入方法有三种:
1.直接输入法2下载法3.编程输入法
下面以例解释下载法
例1．现要根据瑞士地图计算其国土面积。
于是对地图作如下的测量：
以西东方向为横轴，以南北方向为纵轴。
（选适当的点为原点）将国土最西到最东边界在x轴上的区间划取足够多的分点xi，在每个分点处可测出南北边界点的对应坐标y1，y2。
用这样的方法得到下表
x
7.0
10.5
13.0
17.5
34.0
40.5
44.5
48.0
56.0
y1
44
45
47
50
50
38
30
30
34
y2
44
59
70
72
93
100
110
110
110
x
61.0
68.5
76.5
80.5
91.0
96.0
101.0
104.0
106.5
y1
36
34
41
45
46
43
37
33
28
y2
117
118
116
118
118
121
124
121
121
x
111.5
118.0
123.5
136.5
142.0
146.0
150.0
157.0
158.0
y1
32
65
55
54
52
50
66
66
68
y2
121
122
116
83
81
82
86
85
68
根据地图比例知18mm相当于40km，试由上表计算瑞士国土的近似面积。
（精确值为41288km2）。
解题思路：
数据实际上表示了两条曲线，实际上我们要求由两曲线所围成的图形的面积。
解此问题的方法是数值积分的方法。
具体解时我们遇到两个问题：
1。
数据如何输入；2。
没有现成的命令可用。
解：
对于第一个问题，我们可把数据考备成M文件（或纯文本文件）。
然后，利用数据绘制平面图形。
键入
loadmianji.txt
a=mianji;
A=[a（1:
3,1:
9）,a（4:
6,1:
9）,a（7:
9,1:
9）]';
plot（A（:
1）,A（:
2）,'r',A（:
1）,A（:
3）,'g'）
接下来可以计算面积。
键入：
a1=trapz（A（:
1）*40/18,A（:
2）*40/18）;
a2=trapz（A（:
1）*40/18,A（:
3）*40/18）;
d=a2-a1
d=4.2414e+004
至此，问题可以说得到了解决。
之所以说还有问题，是我们觉得误差较大。
但计算方法的理论给了我们更精确计算方法。
只是MATLAB没有相应的命令。
想得到更理想的结果，我们可以自己设计解决问题的方法。
（可以编写辛普森数值计算公式的程序，或用拟合的方法求出被积函数，再利用MATLAB的命令quad,quad8）
2.下面以例解释编程输入法
若矩阵的各元素之间存有某种关系，那么可利用存有的关系，通过编写程序来建立矩阵；还可以利用MATLAB的矩阵命令来建立矩阵。
例：
输入矩阵
fori=1:
5%变量i的步长为1，从1变到5
forj=1:
5%变量j的步长为1,从1变到5
H（i,j）=1/（i+j-1）;%将值1/（i+j-1）赋予矩阵H的第i行第j列的元素
end
end
循环语句
MATLAB也有控制流语句,用于控制程序的流程.主要有for循环、while循环、if和break三种控制语句.虽然语句很少,但功能很强.
ⅰfor循环语句
for循环语句的一般表达形式为：
fori=表达式
可执行语句1
……
可执行语句n
end
例：
求S=1+2+3+…+50,可编程如下：
s=0;
fork=1:
50
s=s+k;
end
ⅱwhile循环
while循环语句用来控制一个或一组语句在某逻辑条件下重复预先确定或不确定的次数.
while循环语句的一般表达形式为：
while表达式
循环体语句
end
例：
对于上面同样的问题,可编程如下：
S=0;k=0;
whilek<51
S=S+k;k=k+1;%当条件k〈51时,反复执行语句S=S+k,k=k+1
end
以上这个例子事先已知循环次数是51,下面再看一个预先不能确定循环次数的例子.
例：
用迭代法
xk+1=e
k=0,1,2,…
求解方程x-e-x=0的根.初始值x0=0.5,相对误差限ε=10-8,编程如下：
ep=10^（-8）;dx=1;x0=0.5;k=0;
whiledx>ep
k=k+1;
x=exp（-x0）;
dx=abs（x-x0）/（1+abs（x））;
x0=x;
end
ⅲif和break语句
MATLAB中if和break语句的作用与使用方式同其它编程语言一样,用来将控制流程
进行分流与中断退出.
例：
可以把上面的解方程的例子中的循环语句改写成:
x0=0;
while
（1）
k=k+1;
x=exp（-x0）;
dx=abs（x-x0）/（1+abs（x））;
ifdx<=epsbreak;end
x0=x;
end
程序中while
（1）说明循环条件总是真,直到满足dx<=eps条件跳出循环体.
ⅳif–else–end分支结构
分支结构有三种形式：
if表达式
执行语句
end
如果表达式的值非0,则执行下面的语句.否则执行end后面的语句.
if表达式
执行语句1
else
执行语句2
end
if表达式1
执行语句1
elseif表达式2
执行语句2
elseif表达式3
执行语句3
……
else（此句可以省略）
执行语句n
end
例：
对函数
可编如下的程序：
ifx<1
y=x;
elseifx>=1&x<10
y=2*x-1;
else
y=3*x-11;
end
练习:
求H矩阵的每一行元素的和,最大值,最小值.
求H矩阵的每一行元素的和,最大值,最小值.
求H矩阵的主对角线元素的和,最大值,最小值.
创建M文件
创建M文件是MATLAB中的非常重要的内容.事实上,正是由于在MATLAB工具箱中存放着大量的M文件,使得MATLAB在应用起来显得简单、方便,且功能强大.M文件有两种形式：
命令文件和函数文件当用户要运行的命令较多时,如果直接在命令窗口中逐条输入和运行,有诸多不便.此时可通过编写命令文件来解决这个问题.另外,从前面的许多例子可以看到：
MATLAB的许多命令,需要用户通过编写函数文件来执行.
ⅰ命令文件的创立
进入MATLAB命令窗口后,选择“file”下拉式菜单中的“new”进入编辑/调试器（Editer/Debugger）,在编辑/调试器中,编写符合语法规则的命令.编写完命令文件后,选择“file”下拉式菜单中的“save”项,然后依提示输入一个文件名.至此,完成了命令文件的创建.
ⅱ函数文件的创立
函数文件的创立方法与命令文件的创立方法完全一样,只是函数文件的第一句可执行语句是以function引导的定义语句,并且输入文件名时要与定义语句中的函数名相同.
建立了函数文件或命令文件后,只要在命令窗口键入命令文件名或函数名,就可执行M文件中所包含的所有命令.
下面分别创建并运行一个命令文件和一个函数文件,以了解M文件的创建和运行的全过程.
计算所有小于1000的Fibonnaci数.
命令文件的创建和运行：
⑴在MATLAB的命令窗口点击“新建”工具栏或在“file”下拉菜单中选“New”中的
“M-file”项,进入编辑/调试器.
⑵在编辑/调试器中,输入以下命令：
%计算小于1000的Fibonnaci数
f=[1,1];
i=1;
whilef（i）+f（i+1）<1000
f（i+2）=f（i）+f（i+1）;
i=i+1;
end
f,i
⑶在“file”下拉菜单中选“Save”项,依提示输入文件名“fibno”至此,完成了命令文件fibno.m的创建.
⑷执行fibno.
在MATLAB窗口中输入fibno并敲回车键,计算机依次执行fibno中的各条命令后显示如下的结果：
ans=
Columns1through12
1123581321345589144
Columns13through16
233377610987
函数文件的创建和运行：
⑴在MATLAB的命令窗口点击“新建”工具栏或在“file”下拉菜单中选“New”中的“M-file”项,进入编辑/调试器.
⑵在编辑/调试器中,输入以下命令：
functionf=ffibno（n）
f=[1,1];
i=1;
whilef（i）+f（i+1）f（i+2）=f（i）+f（i+1）;
i=i+1;
end
f
⑶在“ffile”下拉菜单中选“Save”项,依提示输入文件名“ffibno”至此,完成了函数文件ffibno.m的创建.
⑷执行ffibno.
在MATLAB窗口中输入ffibno（1000）并敲回车键即可.
将上面的关于H矩阵的运算编成M文件和M函数.
6：
微分方程数值解（单摆问题）
单摆问题的数学模型是
在初始角度不大时,问题可以得到很好地解决,但如果初始角较大,此方程无法求出解析解.现问题是当初始角为100和300时,求出其解,画出解的图形进行比较。
解:
若θ0较小,则原方程可用
来近似.其解析解为θ（t）=θ0cosωt,
.
若不用线性方程来近似,那么有两个模型:
取g=9.8,l=25,100=0.1745,300=0.5236.用MATLAB求这两个模型的数值解,先要作如下的处理:
令x1=θ,x2=θ’,则模型变为
再编函数文件
functionxdot=danbai（t,x）
xdot=zeros（2,1）;
xdot
（1）=x
（2）;xdot
（2）=-9.8/25*sin（x
（1））;
在命令窗口键入
[t,x]=ode45（‘danbai’,[0:
0.1:
20],[0.1745,0]）;
[t,y]=ode45（‘danbai’,[0:
0.1:
20],[0.5236,0]）;
plot（t,x（:
1）,’r’,t,y（:
1）,’k’）;
参考书:
数学实验,高等教育出版社
7．解优化问题
线性规划有约束极小问题
模型
用命令
[x,fval]=linprog（f,A,b,A1,b1,lb,ub）
例1：
Findxthatminimizes
f（x）=-5x1-4x2-6x3
subjectto
x1-x2+x3≦20
3x1+2x2+4x3≦42
3x1+2x2≦30
0≦x1,0≦x2,0≦x3
First,enterthecoefficients:
f=[-5;-4;-6]
A=[1-11
324
320];
b=[20;42;30];
lb=zeros（3,1）;
Next,callalinearprogrammingroutine:
[x,fval,exitflag,output,lambda]=linprog（f,A,b,[],[],lb）；
Enteringx,fval,lambda.ineqlin,andlambda.lowergets
x=
0.0000
15.0000
3.0000
fval=
-78.0000
和其它信息。
例2：
解问题
把问题极小化并将约束标准化
键入c=[-2,-3,5];a=[-2,5,-1];b=-10;a1=[1,1,1];b1=7;LB=[0,0,0];
[x，y]=linprog（c,a,b,a1,b1,LB）得当X=（6.4286,0.5714,0.0000）时,z=-14.5714最大.
例3:
解问题
解:
键入
c=[-2,-1,1];a=[1,4,-1;2,-2,1];b=[4;12];a1=[1,1,2];b1=6;lb=[0;0;-inf];ub=[inf;inf;5];
[x,z]=linprog（c,a,b,a1,b1,lb,ub）得当X=（4.6667,0.0000,0.6667）时,z=-8.6667最小.
非线性规划有约束极小问题
模型1
模型2：
MATLAB求解此问题的命令是：
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon（‘fun’,x0,A,b,A1,b1,LB,UB,’nonlcon’,
options,p1,p2,…）
fun是目标函数的m_文件名.nonlcon是约束函数C（x）和C1（x）的m_文件名.文件输出为[C,C1].
例:
求解最优化问题
第1步:
建立目标函数和非线性约束的m_文件.
functiony=e1511（x）%目标函数的m_文件
y=exp（x
（1））*（4*x
（1）^2+2*x
（2）^2+4*x
（1）*x
（2）+2*x
（2）+1）;
function[c1,c2]=e1511b（x）%非线性约束的m_文件
c1=[1.5+x
（1）*x
（2）-x
（1）-x
（2）;-x
（1）*x
（2）-10];
c2=[];
第2步:
运行程序.键入
x0=[-1,1];a1=[1,1];b1=0;
[x,f,exitflab,output]=fmincon（‘e1511’,x0,[],[],a1,b1,[],[],’e1511b’）
得结果.
输出结果的意义:
经过4次迭代（iterations:
4）收敛到了（exitfag=1）最优解x
（1）=-1.2247,x
（2）=1.2247,目标函数最优值为1.8951.
8．回归分析
前面我们曾学过拟合。
但从统计的观点看，对拟合问题还需作回归分析。
例如：
有描述问题甲和问题乙的两组数据（x,y）和（x,z）。
设
x=[01234];y=[1.01.31.5,2.02.3];z=[0.6,1.95,0.9,2.85,1.8];
如果在平面上画出散点图，
plot（x,y,'r.',x,z,'b*'）
那么问题甲的5个点基本在一条直线上而问题乙的5个点却很散乱。
如果都用命令
c1=polyfit（x,y,1）,
c2=polyfit（x,z,1）
来拟合，将得到同一条直线。
自然对问题甲的信任程度会高于对问题乙的信任程度。
所以有必要对所得结果作科学的评价分析。
回归分析就是解决这种问题的科学方法。
x=[01234]';y=[1.01.31.5,2.02.3]';z=[0.6,1.95,0.9,2.85,1.8]';
x=[ones（5,1）,x];
[b,bint,r,rint,stats]=regress（y,x,0.05）;
b=
0.9600
0.3300
bint=
0.76381.1562
0.24990.4101
输出参数的意义：
第一个参数的置信区间：
[0.76381.1562]
第二个参数的置信区间：
[0.24990.4101]
stats=
0.9829171.94740.00100.0063
R2=0.9829F=171.9474,p=0.0010.
R是衡量y与x的相关程度的指标，称为相关系数。
R越大，x与y关系越密切。
通常R大于0.9才认为相关关系成立。
F是一统计指标
p是与F对应的概率，当p<0.05时，回归模型成立。
此例中p=0<0.0010<0.05，所以，所得回归模型成立。
再输入
rcoplot（r,rint）
得图形,说明第3个数据应剔除,将会得到更准确的模型.
对于第2条曲线做同样的工作.
[b,bint,r,rint,stats]=regress（z,x,0.05）;
b
b=
0.9600
0.3300
bint=
-1.11423.0342
-0.51681.1768
置信区间大很多.
stats=
0.33891.53810.30310.7080
说明模型不可靠.
下面结合三个具体的例子介绍MATLAB实现回归分析的命令。
1．合金强度y与其中含碳量x有密切关系，如下表
x
0.10
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.20
0.21
0.23
y
42.0
41.5
45.0
45.5
45.0
47.5
49.0
55.0
50.0
55.0
55.5
60.5
根据此表建立y（x）。
并对结果作可信度进行检验、判断x对y影响是否显著、检查数据中有无异常点、由x的取值对y作出预测。
解1：
在x-y平面上画散点图，直观地知道y与x大致为线性关系。
用命令polyfit（x,y,1）可得y=140.6194x+27.0269。
作回归分析用命令
[b,bint,r,rint,ststs]=regress（y,x,alpha）可用help查阅此命令的具体用法
残差及置信区间可以用rcoplot（r,rint）画图
设回归模型为
y=β0+β1x,
在MATLAB命令窗口中键入下列命令进行回归分析
x=0.1:
0.01:
0.18;x=[x,0.2,0.21,0.23]';
y=[42,41.5,45,45.5,45,47.5,49,55,50,55,55.5,60.5]';
X=[ones（12,1）,x];
[b,bint,r,rint,stats]=regress（y,X,0.05）;
b,bint,stats,rcoplot（r,rint）
得结果和图
b=
27.0269
140.6194
bint=
22.322631.7313
111.7842169.4546
stats=
0.9219118.06700.0000
结果含义为
β0=27.0269β1=140.6194
β0的置信区间是[22.3226,31.7313]
β1的置信区间是[111.7842,169.4546]
R2=0.9219F=118.0670,p<10-4.
R是衡量y与x的相关程度的指标，称为相关系数。
R越大，x与y关系越密切。
通常R大于0.9才认为相关关系成立。
F是一统计指标
p是与F对应的概率，当p<0.05时，回归模型成立。
此例中p=0<10-4<0.05，所以，所得回归模型成立。
观察所得残差分布图，看到第8个数据的残差置信区间不含零点，此点视为异常点，剔除后重新计算。
此时键入：
X（8,:
）=[];
y（8）=[];
[b,bint,r,rint,stats]=regress（y,X）;
b,bint,stats,rcoplot（r,rint）
得：
b=
27.0992
137.8085
bint=
23.856330.3421
117.8534157.7636
stats=
0.9644244.05710.0000
可以看到：
置信区间缩小；R2、F变大，所以应采用修改后的结果。
2．将17至29岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄（x）对这种运动能力（y）的影响。
现得到一组数据如下表
年龄
17
19
21
23
25
27
29
第一人
20．48
25．13
26．15
30．0
26．1
20．3
19．35
第二人
24．35
28．11
26．3
31．4
26．92
25．7
21．3
试建立关系y（x），并作必要的统计分析。
解2：
在x-y平面上画散点图，直观地知道y与x大致为二次函数关系。
x=17:
2:
29;x=[x,x];
y=[20.48,25.13,26.1530,26.1,20.3,19.35,24.35,28.11,26.3,31.4,26.92,25.7,21.3];
plot（x,y,'*'）
设模型为y=a1x2+a2x+a3
此问题可以利用命令polyfit（x,y,2）来解，也可以像上题一样求解。
polyfit（x,y,2）
ans=
-0.20038.9782-72.2150
即
t=17:
0.2:
29;
d=-0.2003*t.^2+8.9782*t-72.215;
plot（t,d）
下面介绍用命令polytool来解。
首先在命令窗口键入
polytool（x,y,2）
得到一个交互式窗口
窗口中绿线为拟合曲线、红线为y的置信区间、可通过移动鼠标的十字线或通过在窗口下方输入来设定x值，窗口左边则输出与x对应的y值及y的置信区间。
通过左下方的Export下拉菜单可输出回归系数等。
更详细的解释可通过help查阅。
3．某厂生产的某产品的销售量与竞争对手的价格x1和本厂的价格x2有关。
下表是该产品在10个城市的销售记录。
x1
120
140
190
130
155
175
125
145
180
150
x2
100
110
90
150
210
150
250
270
300
250
y（个）
102
100
120
77
46
93
26
69
65
85
试建立关系y（x1，x2），对结果进行检验。
若某城市本厂产品售价160（元），对手售价170（元），预测此产品在该城市的销售量。
解3.这是一个多元回归问题。
若设回归模型是线性的，即设y=β0+β1x1+β2x2
那么依然用regress（y,x,alpha）求回归系数。
键入
x1=[120,140,190,130,155,175,125,145,180,150];
x2=[100,110,90,150,210,150,250,270,300,250];
y=[102,100,120,77,46,93,26,69,65,85]';
x=[ones（10,1）,x1',x2'];
[b,bint,r,rint,stats]=regress（y,x）;b,bint,stats,
b=
66.5176
0.4139
-0.2698
bint=
-32.5060165.5411
-0.20181.0296
-0.4611-0.0785
stats=
0.65276.57860.0247351.0445
得
b=
66.5176
0.4139
-0.2698
bint=
-32.5060165.5411
-0.20
升级会员 