matlab在数学建模中的应用.docx
《matlab在数学建模中的应用.docx》由会员分享，可在线阅读，更多相关《matlab在数学建模中的应用.docx（18页珍藏版）》请在冰点文库上搜索。
matlab在数学建模中的应用
Matlab在数学建模中的应用
数学建模是通过对实际问题的抽象和简化，引入一些数学符号、变量和参数，用数学语言和方法建立变量参数间的内在关系，得出一个可以近似刻画实际问题的数学模型，进而对其进行求解、模拟、分析检验的过程。
它大致分为模型准备、模型假设、模型构成、模型求解、模型分析、模型检验及应用等步骤。
这一过程往往需要对大量的数据进行分析、处理、加工，建立和求解复杂的数学模型，这些都是手工计算难以完成的，往往在计算机上实现。
在目前用于数学建模的软件中，matlab强大的数值计算、绘图以及多样化的工具箱功能，能够快捷、高效地解决数学建模所涉及的众多领域的问题，倍受数学建模者的青睐。
1Matlab在数学建模中的应用
下面将联系数学建模的几个环节，结合部分实例，介绍matlab在数学建模中的应用。
1.1模型准备阶段
模型准备阶段往往需要对问题中的给出的大量数据或图表等进行分析，此时matlab的数据处理功能以及绘图功能都能得到很好的应用。
1.1.1确定变量间关系
例1已知某地连续20年的实际投资额、国民生产总值、物价指数的统计数据（见表），由这些数据建立一个投资额模型，根据对未来国民生产总值及物价指数的估计，预测未来的投资额。
表1实际投资额、国民生产总值、物价指数的统计表
记该地区第t年的投资为z（t），国民生产总值为x（t），物价指数为y（t）。
赋值：
z=[90.997.4113.5125.7122.8133.3149.3144.2166.4195229.8228.7206.1257.9324.1386.6423401.9474.9424.5]'
x=[596.7637.7691.1756799873.4944992.71077.61185.91326.41434.21549.217181918.32163.92417.82631.62954.73073]'
y=[0.71670.72770.74360.76760.79060.82540.86790.91450.960111.05751.15081.25791.32341.40051.50421.63421.78421.95142.0688]'
先观察x与z之间，y与z之间的散点图
plot（x,z,'*'）
plot（y,z,'*'）
由散点图可以看出，投资额和国民生产总值与物价指数都近似呈线性关系，因此可以建立多元线性回归模型
直接利用统计工具箱直接计算
[b,bint,r,rint,stats]=regress（z,X,alpha）
输入
z：
n维数据向量
X:
[ones（20,1）xy]，这里的1是个向量，元素全为常数1，即为ones（n,1）
Alpha:
置信水平,一般为0.05
输出
b：
β的估计值
bint:
b的置信区间
r:
残差向量z-Xb
rint:
r的置信区间
Stats:
检验统计量
F,p
代入上述公式
[b,bint,r,rint,stats]=regress（z,X,0.05）
有b=
322.756305635088
0.618516611734168
-859.579151516612
即
由
stats=
0.9972920.4761130081070
知z的99.085%可由模型确定，F远超过F检验的临界值，p远小于α=0.05.
bint=
224.4022221134421.110389156777
0.4773754129901840.759657810478151
-1121.49331646023-597.664986572995
b的置信区间不包含零点，x,y对z影响都是显著的。
z=[90.997.4113.5125.7122.8133.3149.3144.2166.4195229.8228.7206.1257.9324.1386.6423401.9474.9424.5]';
x=[596.7637.7691.1756799873.4944992.71077.61185.91326.41434.21549.217181918.32163.92417.82631.62954.73073]';
y=[0.71670.72770.74360.76760.79060.82540.86790.91450.960111.05751.15081.25791.32341.40051.50421.63421.78421.95142.0688]';
>>X=[ones（20,1）xy];
>>[b,bint,r,rint,stats]=regress（z,X,0.05）
b=
322.7563
0.6185
-859.5792
bint=
1.0e+003*
0.22440.4211
0.00050.0008
-1.1215-0.5977
r=
15.1352
5.7314
2.4699
-4.8419
-14.5678
-20.1721
-11.3072
-6.4726
2.4121
-1.6760
-4.3518
8.0709
6.4024
10.0992
18.6839
18.4146
9.5185
-14.8835
1.9954
-20.6605
rint=
-8.770139.0405
-19.949031.4118
-23.677528.6173
-30.837721.1539
-39.606810.4712
-44.00933.6652
-37.010114.3956
-32.814419.8691
-24.213929.0382
-28.354225.0022
-30.048921.3453
-18.468034.6097
-16.323529.1283
-15.237835.4362
-6.133743.5015
-4.522741.3519
-13.604732.6417
-38.94989.1828
-22.055326.0461
-38.2783-3.0427
stats=
0.9909920.47610161.5988
>>
1.1.2求数字特征
例2已知50个数据x=[451.4243.89527.185312.6912.863383.97683.1292.84235.338612.4608.5415.7616.355190.07586.9257.581367.57631.45717.63692.6784.079454.36441.83353.25153.61675.64699.21727.51478.38554.84121.05450.75715.88892.84273.1254.77865.6232.35804.87908.4231.89239.3149.75478.384640.82190.89843.87173.9170.79994.3]，计算其数字特征。
输入数据，利用下列提供的函数可以求得各数字特征。
min（x）:
向量x的元素的最小值
max（x）:
向量x的元素的最大值
mean（x）:
向量x的元素的算术平均值
geomean（x）:
向量x的元素的几何平均值
（n个正数的连乘积的n次算术根叫做这n个数的几何平均数）
median（x）:
向量x的元素的中位数
var（x）:
向量x的元素的方差
std（x）:
向量x的元素的标准差
diff（x）:
向量x的相邻元素的差
sort（x）:
对向量x的元素进行排序（Sorting）
length（x）:
向量x的元素个数
sum（x）:
向量x的元素总和
prod（x）:
向量x的元素总乘积
1.2模型的求解分析与检验
1.2.1拟合数据做预测
例3以下是美国1790年至2000年的人口统计数据（单位：
百万），建立人口发展模型并预测2010年美国的人口数目。
根据分析，第t年的人口x满足
（指数增长模型）
将上式两边取对数，得
，
，
由t=0:
21,x=[3.95.37.29.612.917.123.231.438.650.262.97692106.5123.2131.7150.7179.3204226.5251.4281.4]
y=log（x）;f=polyfit（t,y,1）,得到
r=0.2022，
=
=6.045
x（22）=516.770百万
1.2.2绘制误差条图
将模型得出的结果与真实结果作比较，绘制出对比图和误差条图，反应模型与实际的吻合程度。
如上例，模型结果与实际人口数的对比图以及误差条图可由命令
t=0:
21,x=[3.95.37.29.612.917.123.231.438.650.262.97692106.5123.2131.7150.7179.3204226.5251.4281.4];
plot（t,x,'*',t,6.0448*exp（0.2022*t）,'o',t,6.0448*exp（0.2022*t））;errorbar（[1790:
10:
2000],ones（1,22）,x-6.0448*exp（0.2022*t））
1.2.3对模型进行模拟
对于一些没有给出数据的实际问题，建立模型后往往需要找一组随机数据进行模拟，从而检验模型的优劣。
例4已知一栋大厦有9部电梯，上下班高峰期和非高峰期上下电梯的人数有显著的差别，为节约用电，试建立数学模型进行电梯的调试。
题中没有给出等电梯的人数，在建立完数学模型后，就可以利用matlab模拟一组各时间段等电梯的人数带入模型求解和检验。
由概率知识知道，到达电梯的人数呈正态分布且在上班之前的某一刻和下班之后的某一刻达到峰值，可以使用
X=normrnd（mu,sigma,1,n）
来生成均值为mu,方差为sigma的一组（n个）随机数来模拟。
2实例分析
实例1（身高问题）
学校随机抽取100名学生，测量他们的身高，得一组数据。
1）根据这些数据对全校学生的平均身高作出估计，并给出估计的误差范围；2）学校10年前作过普查，学生的平均身高为167.8cm，试根据这次抽查的数据，对学生的平均身高有无显著提高作出结论。
身高为h=[161175172172175175180179172174164170158176178178178171168169179163182174160163170160168176163170178178174172170170172180169171170168171179156158171171162175170170154175170168166164170168173162163160160172170172174172175160168170170158169173167164168170171176173169164167167168172163172164172168165160]
解：
（1）分析与假设：
这是一个直接根据观测数据组建模型的问题，需用部分去推断整体，属于随机性数学模型，运用数理统计的方法可解决这一问题。
受测量工具、观测方法等因素的影响，上述测量数据可能会有不同程度的误差，不尽准确，但作为我们分析的基点，可假设这些测量数据是准确的。
（2）模型的建立与求解：
①作学生身高的直方图和频数表，对学生身高作直观考察
hist（h）作出身高直方图
[N,X]=hist（h）作学生的身高频数表
由结果可以使我们对这所学校学生的身高有这样的一些粗略认识：
近70%学生身高在165至175之间，平均约为169，身高的分布大致呈中间高、两端低的钟形，故可以假设为正态分布N（
）.
②对分布作假设检验：
采用正态概率图纸法检验，matlab统计工具箱中提供的是Q-Q图检验：
normplot（h）
由图可知，样本点在一条直线附近，故可得学生身高服从正态分布这一结论。
③考察样本统计量所反映的数据特征：
mean（h）计算样本均值
median（h）计算中位数
std（h）计算标准差
range（h）计算极差
skewness（h）计算偏度
kurtosis（h）计算峰度
均值
中位数
标准差
极差
偏度
峰度
169.44
170
5.9464
28
-0.3242
2.6849
标准差为s=5.9464，说明数据与均值偏离程度不算太大，偏度
这与正态分布是对称的，偏度接近于0这一数学原理相接近。
而峰度
比正态分布的峰度3稍小一些，考虑到样本抽取方法和测量误差，可以认为这一模型是比较合理的，与实际情况比较相符。
④平均身高的估计及误差范围：
此即需由样本去推断总体，由数理统计知识，需对总体均值
和标准差
进行点估计和区间估计。
[musigmamucisigmaci]=normfit（h,0.01）
可得到全校学生平均身高
，标准差
的点估计和区间估计（显著性水平为0.01）
⑤解决平均身高是否有显著提高的问题：
由数理统计知识知，此即需要对总体均值进行假设检验：
。
由于总体标准差
未知，故用t检验，取显著性水平
。
[H,p,ci]=ttest（h,167.8,0.01,1）
得h=1表示拒绝
p=0.0035,ci=[168.0339,inf]
根据这一命令结果可知全校学生平均身高有显著提高。
由该例题可以看出matlab在数学建模中的巨大优势，充分显示了它超强的数值计算、数据处理和图形处理功能，无论是在实际问题的分析阶段，数学模型的建立阶段，还是模型求解、分析阶段，matlab都有其他语言无法比拟的方便、快捷、高效的运用，不论是数值计算，还是图形的描绘，matlab只需要一两个命令就能解决诸如C，C++语言需几十行的程序才能解决的问题，避免了繁杂的数值计算和复杂的程序设计，能使数学建模者将主要的精力放在问题的分析、模型的建立、算法研究等方面，既节约了时间，大大提高了数学建模的效率，又有利于提高数学建模的质量和人们解决实际问题的能力。
另外，如本例，其先进的数据可视化功能，能将一组大规模的杂乱无章的数据通过图形的方式表现出来，根据几何直观，数学建模者能快速而轻易地提到有意义的特征和结果，探索、发现规律，进而较快地找到数学建模的方法，丰富了数学建模的方法和手段，有力地促进了问题的解决。
实例2价格竞争问题
位于同一条公路旁的甲、乙两个加油站彼此竞争激烈。
当甲站突然宣布降价后，乙站根据甲站的售价应如何调整自己的售价，使得既能和甲站竞争，又可以获得尽可能高的利润？
解
（1）问题分析：
加油站的利润主要来自汽油的销售价和销售量。
这场价格战中，乙加油站汽油降价销售主要受以下3个因素影响：
①甲加油站汽油降价的幅度；②乙加油站汽油降价的幅度；③两站之间汽油销售价之差。
（2）模型假设：
①汽油的正常销售价格保持常数不变；②
（1）中的3个因素对乙加油站销售量的影响是线性的。
（3）模型的建立：
引入符号：
P：
汽油的正常销售价格（元/升）
L：
降价前乙加油站的销售量（升/日）
W：
汽油的成本价格（元/升）
a：
因素①对乙加油站汽油销售量影响的比例常数
b:
因素②对乙加油站汽油销售量影响的比例常数
c:
因素③对乙加油站汽油销售量影响的比例常数
x:
乙加油站的销售价格（元/升）
y:
甲加油站的销售价格（元/升）
根据问题的分析和模型的假设，可得乙加油站的利润函数为：
这里的a,b,c>0.
（4）模型的求解：
以上是建立的数学模型，下面用matlab求解：
symsLPWabcxy
f=（x-W）*（L-a*（P-y）+b*（P-x）-c*（x-y））
df=diff（f,x）%求导
x0=solve（'L-a*（P-y）+b*（P-x）-c*（x-y）+（x-W）*（-b-c）=0','x'）%求驻点
x0=1/2*（L-a*P+a*y+b*P+c*y+W*b+W*c）/（b+c）
即当甲加油站把汽油的销售价格降到y元时，乙加油站把汽油的销售价格定为x0时可以使乙加油站获得最高的利润。
（5）模型检验
f1=subs（f,x,x0）将x0代入求乙站的利润函数里的x
f2=subs（f1,{L,P,W,a,b,c},{2000,4,3,1000,1000,4000}）取具体的数据代入
x1=subs（x0,{L,P,W,a,b,c},{2000,4,3,1000,1000,4000}）
y=3.9:
-0.1:
3.4;
x=17/10+1/2*y
f3=（-13/10+1/2*y）.*（-6500+2500*y）求乙站相应的利润
plot（y,f3）绘制乙站的利润函数
由该例题可以看出，matlab具有强大的符号计算功能，这也是其他语言所不具有的。
在一般的数学模型中，往往有很多未知参数，此时要求解该数学模型不得不进行符号运算，这对其他语言来说，编程是非常困难和繁琐的，而matlab则只需几个简单的命令就可以解决问题。
另外，matlab先进的数据可视化功能和方便的绘图功能，可以将数学模型的求解结果用可视化、动态化的形式表现出来，使数学建模者能以视觉图像方式对模型的结果进行观察，作出解释和评价，有利于加深数学建模都对问题本质的进一步认识，进而进一步修改和完善数学模型，使之更加符合实际。
>>t=0:
21;
>>x=[3.95.37.29.612.917.123.231.438.650.262.97692106.5123.2131.7150.7179.3204226.5251.4281.4];
>>y=log（x）;
>>f=polyfit（t,y,1）;
>>f=polyfit（t,y,1）
f=
0.20221.7992
>>plot（t,x,'*',t,6.0448*exp（0.2022*t）,'o',t,6.0448*exp（0.2022*t））
>>errorbar（[1790:
10:
2000],ones（1,22）,x-6.0448*exp（0.2022*t））
>>h=[161175172172175175180179172174164170158176178178178171168169179163182174160163170160168176163170178178174172170170172180169171170168171179156158171171162175170170154175170168166164170168173162163160160172170172174172175160168170170158169173167164168170171176173169164167167168172163172164172168165160];
>>hist（h）
>>[N,X]=hist（h）
N=
239111318201383
X=
Columns1through7
155.4000158.2000161.0000163.8000166.6000169.4000172.2000
Columns8through10
175.0000177.8000180.6000
>>normplot（h）
>>mean（h）
ans=
169.4400
>>std（h）
ans=
5.9464
>>[musigmamucisigmaci]=normfit（h,0.01）
mu=
169.4400
sigma=
5.9464
muci=
167.8782
171.0018
sigmaci=
5.0187
7.2549
>>[H,p,ci]=ttest（h,167.8,0.01,1）
H=
1
p=
0.0035
ci=
168.0339Inf
>>symsLPWabcxy
f=（x-W）*（L-a*（P-y）+b*（P-x）-c*（x-y））
df=diff（f,x）%求导
x0=solve（'L-a*（P-y）+b*（P-x）-c*（x-y）+（x-W）*（-b-c）=0','x'）%求驻点
x0=1/2*（L-a*P+a*y+b*P+c*y+W*b+W*c）/（b+c）
f=
（x-W）*（L-a*（P-y）+b*（P-x）-c*（x-y））
df=
L-a*（P-y）+b*（P-x）-c*（x-y）+（x-W）*（-b-c）
x0=
1/2*（L-a*P+a*y+b*P+c*y+W*b+W*c）/（b+c）
x0=
（1/2*L-1/2*a*P+1/2*a*y+1/2*b*P+1/2*c*y+1/2*W*b+1/2*W*c）/（b+c）
>>f1=subs（f,x,x0）
f1=
（（1/2*L-1/2*a*P+1/2*a*y+1/2*b*P+1/2*c*y+1/2*W*b+1/2*W*c）/（b+c）-W）*（L-a*（P-y）+b*（P-（1/2*L-1/2*a*P+1/2*a*y+1/2*b*P+1/2*c*y+1/2*W*b+1/2*W*c）/（b+c））-c*（（1/2*L-1/2*a*P+1/2*a*y+1/2*b*P+1/2*c*y+1/2*W*b+1/2*W*c）/（b+c）-y））
>>f2=subs（f1,{L,P,W,a,b,c},{2000,4,3,1000,1000,4000}）
f2=
（-13/10+1/2*y）*（-6500+2500*y）
>>x1=subs（x0,{L,P,W,a,b,c},{2000,4,3,1000,1000,4000}）
y=3.9:
-0.1:
3.4;
x=17/10+1/2*y
f3=（-13/10+1/2*y）.*（-6500
升级会员 