matlab在数学建模中的应用.docx
《matlab在数学建模中的应用.docx》由会员分享，可在线阅读，更多相关《matlab在数学建模中的应用.docx（14页珍藏版）》请在冰豆网上搜索。
matlab在数学建模中的应用
Matlab在数学建模中的应用
它大致分为模型准备、模型假设、模型构成、模型求解、模型分析、模型检验及应用等步骤。
这一过程往往需要对大量的数据进行分析、处理、加工，建立和求解复杂的数学模型，这些都是手工计算难以完成的，往往在计算机上实现。
在目前用于数学建模的软件中，matlab强大的数值计算、绘图以及多样化的工具箱功能，能够快捷、高效地解决数学建模所涉及的众多领域的问题，倍受数学建模者的青睐。
1Matlab在数学建模中的应用
下面将联系数学建模的几个环节，结合局部实例，介绍matlab在数学建模中的应用。
1.1模型准备阶段
模型准备阶段往往需要对问题中的给出的大量数据或图表等进行分析，此时matlab的数据处理功能以及绘图功能都能得到很好的应用。
确定变量间关系
例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
-859.579151516612
即
由
stats=
0.990850141482672920.4761130081070
知z的99.085%可由模型确定，F远超过F检验的临界值，p远小于α=0.05.
bint=
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
>>
求数字特征
例250个数据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模型的求解分析与检验
拟合数据做预测
〔单位：
百万〕
根据分析，第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百万
绘制误差条图
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））
对模型进行模拟
对于一些没有给出数据的实际问题，建立模型后往往需要找一组随机数据进行模拟，从而检验模型的优劣。
例4一栋大厦有9部电梯，上下班顶峰期和非顶峰期上下电梯的人数有显著的差异，为节约用电，试建立数学模型进行电梯的调试。
题中没有给出等电梯的人数，在建立完数学模型后，就可以利用matlab模拟一组各时间段等电梯的人数带入模型求解和检验。
由概率知识知道，到达电梯的人数呈正态分布且在上班之前的某一刻和下班之后的某一刻到达峰值，可以使用
X=normrnd（mu,sigma,1,n）
来生成均值为mu,方差为sigma的一组（n个）随机数来模拟。
2实例分析
实例1〔身高问题〕
学校随机抽取100名学生，测量他们的身高，得一组数据。
1〕根据这些数据对全校学生的平均身高作出估计，并给出估计的误差范围；2〕学校10年前作过普查，学生的平均身高为167.8cm，试根据这次抽查的数据，对学生的平均身高有无显著提高作出结论。
身高为h=[161175172172175175180179172174164170158176178178178171168169179163182174160163170160168176163170178178174172170170172180169171170168171179156158171171162175170170154175170168166164170168173162163160160172170172174172175160168170170158169173167164168170171176173169164167167168172163172164172168165160]
〔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这一数学原理相接近。
而峰度
④平均身高的估计及误差范围：
此即需由样本去推断总体，由数理统计知识，需对总体均值
和标准差
进行点估计和区间估计。
[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]
语言需几十行的程序才能解决的问题，防止了繁杂的数值计算和复杂的程序设计，能使数学建模者将主要的精力放在问题的分析、模型的建立、算法研究等方面，既节约了时间，大大提高了数学建模的效率，又有利于提高数学建模的质量和人们解决实际问题的能力。
另外，如本例，其先进的数据可视化功能，能将一组大规模的杂乱无章的数据通过图形的方式表现出来，根据几何直观，数学建模者能快速
实例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）里的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）
>>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+2500*y）
x1=
17/10+1/2*y
x=
3.65003.60003.55003.50003.45003.4000
f3=
1.0e+003*
2.11251.80001.51251.25001.01250.8000
>>plot（y,f3）
>>
hights'
ans=
11301250128012301040900500700
13201450142014001300700900850
139015001500140090011001060950
15001200110013501450120011501010
15001200110015501600155013801070
15001550160015501600160016001550
1480150015501510143013001200980
>>x=1200:
400:
4000;
y=1200:
400:
3600;
>>mesh（x,y,hights'）
>>xi=1200:
100:
4000;
>>yi=1200:
100:
3600;
>>zi=interp2（x,y,hights',xi',yi,'cubic'）;
>>mesh（xi,yi,zi）
>>
升级会员 