matlab在数学建模中的应用.docx

1、matlab在数学建模中的应用Matlab在数学建模中的应用 它大致分为模型准备、模型假设、模型构成、模型求解、模型分析、模型检验及应用等步骤。这一过程往往需要对大量的数据进行分析、处理、加工，建立和求解复杂的数学模型，这些都是手工计算难以完成的，往往在计算机上实现。在目前用于数学建模的软件中，matlab 强大的数值计算、绘图以及多样化的工具箱功能，能够快捷、高效地解决数学建模所涉及的众多领域的问题，倍受数学建模者的青睐。1 Matlab在数学建模中的应用下面将联系数学建模的几个环节，结合局部实例，介绍matlab在数学建模中的应用。1.1 模型准备阶段模型准备阶段往往需要对问题中的给出的大

2、量数据或图表等进行分析，此时matlab的数据处理功能以及绘图功能都能得到很好的应用。 确定变量间关系例1 某地连续20年的实际投资额、国民生产总值、物价指数的统计数据见表，由这些数据建立一个投资额模型，根据对未来国民生产总值及物价指数的估计，预测未来的投资额。表1 实际投资额、国民生产总值、物价指数的统计表记该地区第t年的投资为z(t)，国民生产总值为x(t)，物价指数为y(t)。赋值：z=90.9 97.4 113.5 125.7 122.8 133.3 149.3 144.2 166.4 195 229.8 228.7 206.1 257.9 324.1 386.6 423 401.9

3、474.9 424.5x=596.7 637.7 691.1 756 799 873.4 944 992.7 1077.6 1185.9 1326.4 1434.2 1549.2 1718 1918.3 2163.9 2417.8 2631.6 2954.7 3073y=0.7167 0.7277 0.7436 0.7676 0.7906 0.8254 0.8679 0.9145 0.9601 1 1.0575 1.1508 1.2579 1.3234 1.4005 1.5042 1.6342 1.7842 1.9514 2.0688先观察x与z之间，y与z之间的散点图plot(x,z,*)p

4、lot(y,z,*)由散点图可以看出，投资额和国民生产总值与物价指数都近似呈线性关系，因此可以建立多元线性回归模型直接利用统计工具箱直接计算b,bint,r,rint,stats=regress(z,X,alpha) 输入z：n维数据向量X:ones(20,1) x y，这里的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.7

5、56305635088 -859.579151516612即由stats = 0.990850141482672 920.476113008107 0知z的99.085%可由模型确定，F远超过F检验的临界值 ，p远小于=0.05 .bint = 0.477375412990184 0.759657810478151 -1121.49331646023 -597.664986572995b的置信区间不包含零点，x,y对z影响都是显著的。 z=90.9 97.4 113.5 125.7 122.8 133.3 149.3 144.2 166.4 195 229.8 228.7 206.1 257.

6、9 324.1 386.6 423 401.9 474.9 424.5;x=596.7 637.7 691.1 756 799 873.4 944 992.7 1077.6 1185.9 1326.4 1434.2 1549.2 1718 1918.3 2163.9 2417.8 2631.6 2954.7 3073;y=0.7167 0.7277 0.7436 0.7676 0.7906 0.8254 0.8679 0.9145 0.9601 1 1.0575 1.1508 1.2579 1.3234 1.4005 1.5042 1.6342 1.7842 1.9514 2.0688; X=

7、ones(20,1) x y; b,bint,r,rint,stats=regress(z,X,0.05) b = 322.7563 0.6185 -859.5792bint = 1.0e+003 * 0.2244 0.4211 0.0005 0.0008 -1.1215 -0.5977r = 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

8、.6605rint = -8.7701 39.0405 -19.9490 31.4118 -23.6775 28.6173 -30.8377 21.1539 -39.6068 10.4712 -44.0093 3.6652 -37.0101 14.3956 -32.8144 19.8691 -24.2139 29.0382 -28.3542 25.0022 -30.0489 21.3453 -18.4680 34.6097 -16.3235 29.1283 -15.2378 35.4362 -6.1337 43.5015 -4.5227 41.3519 -13.6047 32.6417 -38

9、.9498 9.1828 -22.0553 26.0461 -38.2783 -3.0427stats = 0.9909 920.4761 0 161.5988 求数字特征例2 50个数据x=451.42 43.895 27.185 312.69 12.863 383.97 683.1 292.842 35.338 612.4 608.54 15.76 16.355 190.07 586.92 57.581 367.57 631.45 717.63 692.67 84.079 454.36 441.83 353.25 153.61 675.64 699.21 727.51 478.38 554

10、.84 121.05 450.75 715.88 892.84 273.1 254.77 865.6 232.35 804.87 908.4 231.89 239.31 49.754 78.384 640.82 190.89 843.87 173.9 170.79 994.3，计算其数字特征。输入数据，利用以下min(x): 向量x的元素的最小值 max(x): 向量x的元素的最大值 mean(x): 向量x的元素的算术平均值 geomean(x):向量x的元素的几何平均值n个正数的连乘积的n次算术根叫做这n个数的几何平均数median(x): 向量x的元素的中位数 var(x):向量x的元素

11、的方差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.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4 28

12、1.4y=log(x);f=polyfit(t,y,1),得到r=0.2022，=6.045x(22)=516.770百万 绘制误差条图t=0:21,x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4 281.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

13、) 对模型进行模拟对于一些没有给出数据的实际问题，建立模型后往往需要找一组随机数据进行模拟，从而检验模型的优劣。例4 一栋大厦有9部电梯，上下班顶峰期和非顶峰期上下电梯的人数有显著的差异，为节约用电，试建立数学模型进行电梯的调试。题中没有给出等电梯的人数，在建立完数学模型后，就可以利用matlab模拟一组各时间段等电梯的人数带入模型求解和检验。由概率知识知道，到达电梯的人数呈正态分布且在上班之前的某一刻和下班之后的某一刻到达峰值，可以使用X=normrnd(mu,sigma ,1 ,n)来生成均值为mu,方差为sigma的一组(n个)随机数来模拟。2 实例分析实例1 身高问题学校随机抽取100

14、名学生，测量他们的身高，得一组数据。1根据这些数据对全校学生的平均身高作出估计，并给出估计的误差范围；2学校10年前作过普查，学生的平均身高为167.8cm，试根据这次抽查的数据，对学生的平均身高有无显著提高作出结论。身高为h=161 175 172 172 175 175 180 179 172 174 164 170 158 176 178 178 178 171 168 169 179 163 182 174 160 163 170 160 168 176 163 170 178 178 174 172 170 170 172 180 169 171 170 168 171 179 15

15、6 158 171 171 162 175 170 170 154 175 170 168 166 164 170 168 173 162 163 160 160 172 170 172 174 172 175 160 168 170 170 158 169 173 167 164 168 170 171 176 173 169 164 167 167 168 172 163 172 164 172 168 165 1602模型的建立与求解：作学生身高的直方图和频数表，对学生身高作直观考察hist(h) 作出身高直方图N,X=hist(h) 作学生的身高频数表由结果可以使我们对这所学校学生的身

16、高有这样的一些粗略认识：近70%学生身高在165至175之间，平均约为169，身高的分布大致呈中间高、两端低的钟形，故可以假设为正态分布N().对分布作假设检验：采用正态概率图纸法检验，matlab统计工具箱中提供的是Q-Q图检验：normplot(h)由图可知，样本点在一条直线附近，故可得学生身高服从正态分布这一结论。考察样本统计量所反映的数据特征：mean(h) 计算样本均值median(h) 计算中位数std(h) 计算标准差range(h) 计算极差skewness(h) 计算偏度kurtosis(h) 计算峰度均值中位数标准差极差偏度峰度169.441705.946428-0.324

17、22.6849标准差为s=5.9464，说明数据与均值偏离程度不算太大，偏度,这与正态分布是对称的，偏度接近于0这一数学原理相接近。而峰度平均身高的估计及误差范围：此即需由样本去推断总体，由数理统计知识，需对总体均值和标准差进行点估计和区间估计。mu sigma muci sigmaci=normfit(h,0.01)可得到全校学生平均身高，标准差的点估计和区间估计显著性水平为0.01)解决平均身高是否有显著提高的问题：由数理统计知识知，此即需要对总体均值进行假设检验：。由于总体标准差未知，故用t检验，取显著性水平。H,p,ci=ttest(h,167.8,0.01,1)得h=1表示拒绝,p=

18、0.0035,ci=168.0339,inf语言需几十行的程序才能解决的问题，防止了繁杂的数值计算和复杂的程序设计，能使数学建模者将主要的精力放在问题的分析、模型的建立、算法研究等方面，既节约了时间，大大提高了数学建模的效率，又有利于提高数学建模的质量和人们解决实际问题的能力。另外，如本例，其先进的数据可视化功能，能将一组大规模的杂乱无章的数据通过图形的方式表现出来，根据几何直观，数学建模者能快速实例2 价格竞争问题位于同一条公路旁的甲、乙两个加油站彼此竞争剧烈。当甲站突然宣布降价后，乙站根据甲站的售价应如何调整自己的售价，使得既能和甲站竞争，又可以获得尽可能高的利润？解 1问题分析：加油站的

19、利润主要来自汽油的销售价和销售量。这场价格战中，乙加油站汽油降价销售主要受以下3个因素影响：甲加油站汽油降价的幅度；乙加油站汽油降价的幅度；两站之间汽油销售价之差。2模型假设：汽油的正常销售价格保持常数不变；1中的3个因素对乙加油站销售量的影响是线性的。3模型的建立：P：汽油的正常销售价格元/升L：降价前乙加油站的销售量升/日W：汽油的本钱价格元/升a：因素对乙加油站汽油销售量影响的比例常数b: 因素对乙加油站汽油销售量影响的比例常数c: 因素对乙加油站汽油销售量影响的比例常数x:乙加油站的销售价格元/升y:甲加油站的销售价格元/升这里的a,b,c0.(4)模型的求解：以上是建立的数学模型，下

20、面用matlab求解：syms L P W a b c x yf=(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)里的xf2=subs(f1,L,P,W,a,b,c,2000,4,3,1000,1000,4000

21、) 取具体的数据代入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*yf3=(-13/10+1/2*y).*(-6500+2500*y) 求乙站相应的利润plot(y,f3) t=0:21; x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4 281.4; y=log(x); f=polyfit(t,y,1); f=polyfit(t,y

22、,1)f = 0.2022 1.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=161 175 172 172 175 175 180 179 172 174 164 170 158 176 178 178 178 171 168 169 179 163 182 174 160 163 170 160 168 176 163 170 178 178 174 172 170 170 172 180 16

23、9 171 170 168 171 179 156 158 171 171 162 175 170 170 154 175 170 168 166 164 170 168 173 162 163 160 160 172 170 172 174 172 175 160 168 170 170 158 169 173 167 164 168 170 171 176 173 169 164 167 167 168 172 163 172 164 172 168 165 160; hist(h) N,X=hist(h) N = 2 3 9 11 13 18 20 13 8 3X = Columns 1

24、 through 7 155.4000 158.2000 161.0000 163.8000 166.6000 169.4000 172.2000 Columns 8 through 10 175.0000 177.8000 180.6000 normplot(h) mean(h)ans = 169.4400 std(h)ans = 5.9464 mu sigma muci sigmaci=normfit(h,0.01)mu = 169.4400sigma = 5.9464muci = 167.8782 171.0018sigmaci = 5.0187 7.2549 H,p,ci=ttest(

25、h,167.8,0.01,1)H = 1p = 0.0035ci = 168.0339 Inf syms L P W a b c x yf=(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

26、 = 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

27、*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*yf3=(-13/10+1/2*y).*(-6500+2500*y) x1 = 17/10+1/2*y x = 3.6500 3.6000 3.5500 3.5000 3.4500 3.4000f3 = 1

28、.0e+003 * 2.1125 1.8000 1.5125 1.2500 1.0125 0.8000 plot(y,f3) hightsans = 1130 1250 1280 1230 1040 900 500 700 1320 1450 1420 1400 1300 700 900 850 1390 1500 1500 1400 900 1100 1060 950 1500 1200 1100 1350 1450 1200 1150 1010 1500 1200 1100 1550 1600 1550 1380 1070 1500 1550 1600 1550 1600 1600 1600 1550 1480 1500 1550 1510 1430 1300 1200 980 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)