《数学建模讲座一一计算软件的使用.pdf》由会员分享，可在线阅读，更多相关《数学建模讲座一一计算软件的使用.pdf（19页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、数学中国(www.madio.n目专业的数学建模论坛l数学建模讲座一一计算软件的使用司守奎，海军航空工程学院数学教研室Em血1:sishoukui163m，QQ:94339146 一般来说学习数学建模，常用的软件有五种，分别是Matlab、Lingo、Mathematica、SPSS和SAS。实际上使用什么软件都可以，关键是在论文中把数学原理写清楚。如果你搞清了算法，你也可以自己编程从底层一步一步地实现算法。当然如果我们借助于Matl曲的强大工具箱，可以起到事半功倍的作用。数学建模涉及的知识和算法特别多，我认为至少要掌握数学规划，图论，预测方法，评价方法这四大部分知识，其它的数学知识和算法可能
2、就需要数学建模竞赛时现学现实了。这些知识我就不展开讲了，相关的数学建模书上都有这些内容。我新版的数学建模教材数学建模算法及应用c国防工业出版社出版)，第14章把评价方法总结了一下，第15章把预测方法也总结了一下。新版的数学建模教材相对于网上流传的算法大全的电子文档，对一些内容进行了压缩，增加了数字图像处理的知识，并增加了一些抽伽b工具箱的应用。对于数学规划的模型，建议大家使用Lingo软件求解比较方便，对于其它问题，如时间序列模型，你使用什么软件求解都可以，关键看个人的喜好和对某种软件的熟悉程度，例如你可以使用SPSS，SAS或Eviews，R软件等求解时间序列模型。我对Matla如软件比较熟
3、悉，一般问题我都是使用陆tlab软件进行求解。由于时间的关系，我主要讲rMat1ab软件的使用。x=4cos，y=(5+4sin)s z=(5+4sin)sin 其中，0，21。画图的Matlab程序如下。a1pha=0:0.1:2与i;beta=咀:0.1:2、li;x=4*cos(a1pha)*ones(size(b由);y=(5+4、in(a1pha)丰cos(beta);z=(5+4丰sin(a1pha)*sin(be饵);surf1:x品z)画图的脑咀ab程序也可以写成x=(a1pha，b由)4*cos(a1pha);y=(a1pha,beta)(5叫*sin(a1pha)*cos(
4、b眈a);z=(a1pha,beta)(5+4、in(a1pha)*sin(beta);但即你品z)2011年数学中国公益讲座一一数学建模之软件实现(司守奎net 量学中国(WWW.皿)专业的盘学噩摸桔坛l翼叫co.(o:)户.叫E拖到0.(13.).-(5+-4叫哺剧。自】2 ，。圄l例1的图形对于其它的二次曲面，如果可以写成单支的显函数，直接使用命令.，田h或cz，田f画图，杏则必须先化成参数方程.例2绘制工元函数z=SE叫昂，)昂的三维表面图.程序一:r(可，)，叫x*y:+3).1(1+玄A勾;%注意这里函数适用范围要广，对于向量也可逐点计算取值I(uadl仰，1)但U6计算工重积分数
5、学申锢-其中D:x+y耳4.F二(x，y)(4-x.2-y.2).*(x.2+予2斗);%定义被积函数和积分区域合成后的匿名函数I=，Iblq皿d(-2，2，-2，2)%求二重积分的数值解例7计算1=JJJ(x+叫2申忡，其中的ZX2+向x2+y2+Z2日所围。成的区域。l=队y.z)(x.句句:).2.*(咀2+y.2&x.2+y.2+z.2=.刀，I=tripl吨uad(f，-甲啊，町响，-町响1，叭2)，0叭2)3拟合拟舍的准则多数情形下使用的是最小工乘法，即误差平方和最小z一些经济学问题上的拟舍使用最小一乘法，即误差绝对值和最小.3.1 线性最小二乘姑Mat1ab中的线性最小二乘法实际
6、上是解线性方程组，对应的抽血b命令是斜杠飞.对于非齐次线性方程组Ax=b，抽血b的求解命令为x-Ab 无论Ax=b有唯一解，无穷多解还是无解时，陆咀曲总是给出一个解.(J)当Ax=b唯一解时，陆.tl曲给出的就是该唯一解.。当Ax=b有无穷多解时，Mat1ab给出的是最小范数解.当Ax=b无解时，陆.tl曲给出的是最小工乘解.2011年撞学中国告盐讲盛一一数学直模之软件实现司守童数学中国(www.madio.n目专业的数学建模论坛lJ 8 利用表1的数据拟舍函数y=ao+aj码+a2x2中的参数饨，a20 表1.l1 X2 y 5 2 1.5 5 5 2.3 7 10 3.5 8 4 5.2
7、8 6 1.8 3 3 2.9 7 8 3.8 7 3 4.5 2 6 4.7 拟含参数的Matlab程序为数学:中国-m=size(a,l);%求矩阵a的行数目血u=on国(血，1)，a(:，1，2);%构造待定参数作为未知数的线性方程组的系数阵cs=xishua(:,3)%拟合参数求得ao=2.9051,=0.009,a2=0.07630 J 9 利用表l的数据拟合函数中的参数。0叫，饨6电2乌龟x，y=aoe-e 解把函数y=aoe响e电吨两边取对数得lny=lnao+ajxj+a2x2，拟舍的Matlab程序如下clc,cl国ra=5 2 1.5 5 5 2.3 7 10 3.5 8
8、4 5.2 8 6 1.8 3 3 2.9 7 8 3.8 2011年数学中国公益讲座一一数学建模之软件实现(司守奎量学中国(www.皿)专业的盘学噩摸桔坛l7 3 4.5 2 6 4.7;由-size(a，l);回血IF四届(m，I)，叫，1，2);%构造线性方程组的系数阵1rlog(叫，3日;%计算线性方程组的常数项列cs=xishub;%拟舍1og(y)中的未知参数es(IF呻叫1)%求函数y中的参数计算的结果为a，=2.6012，吨=-).0072，叫=0.0415.3.2 用户图形界面解法Matl.b中拟舍一元函数可以使用用户图形界面解法命令cftoo1.拟合二元函数可以使用命奇sf
9、toolo例10(续例引用用户图形界面命令sftool拟合y=aoe响e响中的参数吨，饨，a20操作的步骤如下(J)首先把数据加载到Matlab工作空间中g运行如下程序数掌:;中国-(2)在Matlab命令窗口运行sftool.得到如下界面I.(.皿.t.;./i,.-.1.101 r;t j 1,_,.)二二,.1.-一一飞E侧也1(0.)二二二三呻1.)二二二v回刷.-,伫斗NEY 圄5韧蛐界面曰缸.f.巴1岳二-斗子(3)输入自变量X，Y的数据x1.x2.困变量z的数据Z，并定义要拟舍的函数.0事exp(al*x)*exp(a2句)设置完成以后的界面如下2011年撞学中国告盐讲盛一一数学
10、直模之软件实现司守童量学中国(www.皿)专业的盘学噩摸桔坛l当F唱f1，也时1.,.,fi,I,.G:二二一一如1.，.)二二k勤，由，。由.,.,10.I C.,.!.旦旦曳.l 巾吃些)*exp.(:,1 ,-.6 d I r;，o.，_二.1.矿,固6参戴世置好嗣后曲界面.1.1.巾也町1.咱1;，邸，2d山.。氧M川.tJ.刷.tJ.,J.o!o tl:.,.:.,回%!i阳旦旦兰且坠坦坦旦旦.,t 委|LE吉1It1*exp(tll*叶*exp(:J,0._ 固7拟告结果界面c:;:;:!:I J司.net(5)使用菜单fit中的子菜单Saveto Workspace飞也可以把计算
11、结果输出到Mat1a.b工作空间中.3.3非线性拟合Ma.t1a.b中非线性拟告的命守有lsqcurvefit和nlinfit，用户图形界面命nlintool.下面通过一个例于说明lsqcurvefit的使用.例11利用表2给出的近两个世纪的美国人口统计数据(以百万为单位).建立人口预测模型，最后用它预报200年美国的人口.囊2量国人口统计噩噩1建模与求解2011年撞学中国告盐讲盛一一数学直模之软件实现司守童数学中国(www.madio.n目专业的数学建模论坛l记x(t)为第t年的人口数量，设人口年增长率r(x)为x的线性函数，r(x)=r-sx。自然资源与环境条件所能容纳的最大人口数为xm即
12、当X=Xm时，增长率r(xm)=0，可得r(x)=r(l-王)，建立Logistic人口模型xm 其解为2.参数估计ltt卜巳生纠勾=叮叫=r(1叫呻(1dt x(乌)=xo，x(t)=一一二三旦1+(盐l)e呻句)Xo(1)把表2中的全部数据，包括后面的四个空格，全部保存到纯文本文件也.ta2.txt中。(1)非线性最小二乘估计把表2中的第1个数据作为初始条件，利用余下的数据拟合式(1)中的参数xm和r，F看I.Ml!tla:l担庄姐二f1.1 唱E 钮，cl1-十囚一-www.maO-7aextre町data2.txt);元1原始数据保存在纯文本文乎Fdata2.tit中x=a(2:2:6
13、1,:);%提出人口数据x=nonzeros(x);%去掉后面的零，并变成列向量在=1790:10:20;份=叫1);xl萨=x(I);岛n=(囚，td)四(1).1(1+(臼(1)/xO-l)句xp(-四(2)气时-);%臼(1)=xrn，因(2)=r%注意上面定义的要拟合的匿名函数，自变量要具有广泛使用性，%对向量也可以逐点计算函数值，否则调用Isqcurvefit时出错。=lsqcurvefi胁时m邸，1)，叫2:end)，x(2:end)，zeros(2，1)%非线性问题有时很难求得全局最优解，如果拟合参数不稳定，%耍限制参数向量的下界和上界，计算结果就稳定了。她at=fun(囚，t)
14、%预测己知年代人口x201非击皿(cs，2010)%预测tl2010年的人口求得xm=342.4368,r=0.0274,2010年人口的预测值为282.68百万。(2)线性最小二乘法为了利用简单的线性最小二乘炫估计这个模型的参数r和xm把Logistic方程表示为1 dx r-=r-sx,S=-,(2)x dt xm 2011年数学中国公益讲座一一数学建模之软件实现(司守奎数学中国(www.madio.n目专业的数学建模论坛l利用向后差分，得到差分方程(k)-x(k-1)1&x(k)=r-sx(k),k=2,3,22,(3)其中步长&=10，下面拟合其中的参数r和s0编写如下的Mat1ab程
15、序clc,cl田ra=吃extr四d(也.ta4.txt);%把原始数据保存在纯文本文件也.ta4.txt中x=a(2:2:6矿;x:二nonzeros(x);t=1790:10:201;a=阻四(21，1)，-x(2:end);妒diff(x)./x(2:end)110;cs=a飞b;F臼(1)，XIIFr/cs(2)求得xm=373.5135,r=0.0247 0 也可以利用向前差分，得到差分方程(4)咽在=1790:10:201;a=on四(21，1)，-x(I:end-l);b=diff(x)./x(1:end-l)/10;cs=a飞b;F臼(1)，XllF守ics(2)求得xm=29
16、4.386,r=0.0325 0 从上面的三种拟合方法可以看出，拟合同样的参数，方法不同可能结果相差很大。4 Mat lab图论工具箱4.1 工具箱命令介绍Mat1ab图论工具箱的命令见表30表3Matlab固论工具箱的相关命令命令名功能graphal1sh世恒S陪她s求图中所有顶点对之间的最短距离graphconnc础P找无向图的连通分支，或有向图的强(弱连通分支graphisdag 测试有向图是否含有圃，不含圄返回1，否则返回。graphisomorphism 确定两个图是否同构，同构返回1，否则返回。graphisspsn缸田确定一个图是否是生成树，是返回1，否则返回。graphmaxf
17、low 计算有向图的最太流2011年数学中国公益讲座一一数学建模之软件实现(司守奎量学中国(www.皿)专业的盘学噩摸桔坛lZ叩h函皿组回g甲地m且与，叩.回胆thg叩k句Z甲btravcnc4.2 数据结构在固中找量小生血树把前驱顶点序列变成路径的顶点序列求固中指定的-JIf顶点间的最短距圈和最且路径执行有向无圃固的拓扑排序求从一顶点出盘，所酷擅历固中的JJi点Matlab图论工具箱使用的数据结构是稀疏矩阵.稀疏矩阵在数学上是指矩阵中零元素很多，非零元章很少的矩阵.对于稀疏矩阵，只要存放非零元素的行标、到j标、非零元素的值即可，可以按如下方式存储非零元素的行地址，非零元素的列地址)，非零元素
18、的值.在Matl曲中无向图和有向图邻接矩阵的使用上有很大差异.对于有向圈，只要写出邻接矩阵，直接使用抽血b的命+叩E田命+，就可以把邻接矩阵转化为稀疏矩阵的表示方式.对于无向图，由于邻接矩阵是对称阵，Matlab中只需使用邻接矩阵的下三角元素，即Matlab只存储邻接矩阵下三角元素中的非零元素.稀疏矩阵只是一种存储格式.蛐咀曲中，普通矩阵使用事rso命令变成稀疏矩阵，稀疏矩阵使用full命令变成普通矩阵.4.3应用举倒下面给出图论工具箱在最短路，最小生成树和最大流中的应用例子.例12求圄6中从顶点1到顶点7的最短距离.数学中国-mp旦我 圄8元向圃求解的抽d曲程序为clc;cl国X，a=zcr
19、os(;%邻接矩阵韧始化，下面首先输入邻接矩阵的上三角元素叫1，2)50;a(1，3;叫2，4，5)65，40;%即a(2，4)=司55，a(2，句=40叫3，4，7)52，45;a(4,5:7=50,30,42;叫itoz70，a=o;%上三角矩阵变成下三角矩阵a-甲皿O(a);%构造Matlab使用的稀疏矩阵d，阳也l=graphshortes胆创:a，l，7，Direc时，0)%这里是无向圈，属性D回国的取值为oh-vicwiogra抖卫队，ShowArrows，off，ShowWcigl曲，on)%画出该无向图sot(h.Nod届恫也，Col时，10.40.4)%下面是画出最短路径，参
20、见Ma出b的帮助fowEdgos-g幽静sbynodci他got(h.No由咆础l)，ID);2011年撞学中国告盐讲盛一一数学直模之软件实现司守童数学中国(www.madio.n目专业的数学建模论坛lrevEdges=ge旬dg国:bynodeid(h，get(h.Nod四(tliplr(pa也)，ID);edges=岛wEdges;revEdg四;set(edg田，LineColor，100)set(edg臼，LineWi创1，1.5)save mydata a%把稀疏矩阵保存到町data皿t中，供下面的例子应用例13图8中图的最小生成树。cic;clear;loadmyda阻%把稀疏矩阵
21、a加载到Mat1ab工作壁间T=grapbminspantree(a)%调用求最小生成树的Mat1ab命令例14图9的有向图中弧上的数字表示容量，求从V，到V，的最大流量。v 6 主44 2 2、，V,.4.?v,v3 5 7/亏2 v2 4 V,固9有向固5 随机模拟蒙特卡洛方法也称为计算机随机模拟方法，它源于世界著名的赌城摩纳哥的MonteCarlo(蒙特卡洛)0它是基于对大量事件的统计结果来实现一些确定性问题的计算。使用蒙特卡洛方法必须使用计算机生成相关分布的随机数，Matlab给出了生成各种随机数的命令。5.1 产生随机数的命令(1)随机数矩阵Matlab中可以产生20多种分布的随机数
22、，下面列举部分常用随机数的命令。rand(m，n)产生mXn矩阵，其中的元素是服从0，1上均匀分布的随机数。r扭曲t(m，n，四n，max)产生mXn矩阵，其中的元素是min，max上的随机整数。normrod(mu，si胆战m，n)产生mXn矩阵，其中的元素是服从均值为凹，标准差为S1伊a的正态分布的随机数。exprnd(mu，m，n)产生mXn矩阵，其中的元素是服从均值为皿的指数分布的随机数。2011年数学中国公益讲座一一数学建模之软件实现(司守奎数学中国(www.madio.n目专业的数学建模论坛lpoissmd(mu，m，n)产生mXn矩阵，其中的元素是服从均值为m的泊松(Poisso
23、n)分布的随机数。unifrnd(a丸m，n)产生mXn矩阵，其中的元素是服从区间a，b上均匀分布的随机数。r=mvnrnd(MU，SIGMA，四S田)产生S因对均值向量为MU，协方差阵为SIGMA的多维正态分布的随机数.(2)随机置换r皿d阳m(n)产生1到E的一个随机全排列。阳ms(1叫)产生1到E的所有全排列。5.2 随机事件的模拟(1)贝努里事件的模拟在一次随机试验中，事件A或者发生，或者不发生，事件A发生的概率为0.6。可以用Matlab产生一个服从0，1均匀分布的随机数，若产生的随机数小于等于0.6，我们认为事件A发生，这里实际上是一个儿何概率的问题，它的样本空间的长度为1，区间0
24、,O.6的长度为0.6;当然根据实际需要，我们可以认为产生的随机数大于等于0.4时，事件A发生了。例14事件A发生的概率是0.6，模拟在1000次试验中，事件A发生的次数。Matlab程序如下c1c,c1ear a=rand(1,1000);N=su皿(a(=咀.6)(2)离散随机事件的模拟赌轮法例15事件A，B,C发生的概率见表40主Il.L、-.L版单可随机事件发生的概率唱织子气!4问-tj川，wF 累积概率分布为0.4，0.7,1，如果产生的随机数在0，O.4)上，认为事件A发生，产生的随机数在0.4，O.7)上，认为时间B发生，产生的随机数在0.7,1J上，认为事件C发生。Matlab
25、程序如下c1c,c1ear a=0.4,0.3,0.3J;b=cumsum(a)%求累积概率分布c=rand%产生一个0，1区间上均匀分布的随机数d=(c(b(1)+2*(b(1)(=c&c(b(2)+3*(b(2)(=c&c(=b(3)%模拟结果中l表示A，2表示B，3表示C(3)抓阁的模拟n个人抓阉的模拟使用命令rar巾erm(刑。5.3 单重积分计算方法设区间(a，b)中的随机变量占的概率密度函数为马(吟，g(x)是区间怡，b)上的连续函数，数学期望Eg(q)=g(x)!.2+y.2&x.2+y.2+z.22);%计算落在区域V上的频数Vf1n*4.sqrt(2)%计算体积恤c=h(x品
26、z);%计算被积函数一系列的取值jifenV丰sum(hh.气.2+y.2&x.2+y.2+z.2.2+y.2&x.2+y.2+z.22);%计算落在区域V上的频数Vf1n*16忖伊(2)%计算体积恤c=h(x品z);%计算被积函数一系列的取值jifenV.sur叫hh卢(庐.2+y.2&x.2+y.2+z.2=97&da饱.1 1(作下列假设检验Ho:序列X，平稳，H，:序列X，非平稳(存在上升或下降趋势).Daniel检验方法z对于显著水平，由时间序列X，计算(t，乓)，t=1，2，n的Spearman秩相关系数乱，若ITI乌12(n-2)，则拒绝Ho认为序列非平稳。且当q，0时，认为序列
27、有上升趋势q，X;(L)时拒绝Ho即认为Ii，非白噪声，模型检验未通过:而当矿x;(L-r)时，接受Ho认为Ii，是自噪声，模型通过检验。模型残差通过了LB检验，说明模型是适宜的。计算的Matlab程序如下(12)clc,clear a司目位四d(也饱5.txt);%把表5中的数据保存在纯文本文件也时剧中a=a;a=a(:);a(end-3:end)=;R同e世ank(a)%求原始时间序列的秩n=lengtb(a);仁=l:n;t2=1949:2004;plot(t2,a)。=1-6/(町的-1)*sum(t-Rt).2)%计算Qs的值2011年数学中国公益讲座一一数学建模之软件实现(司守奎数
28、学中国(www.madio.n目专业的数学建模论坛l仁=Qs.sqrt(n-2)/sqrt(1-Qs2)%计算T统计量的值L非也v(O.975，n-2)%计算上alphal2分位数figure,subplot(121),autocorr(a)%画自相关函数图subplot(122),parc田T(a)%画偏相关函数图cs=ar(a,l)%拟合模型res=resid(cs,a);%计算残差向量h=lbqtes叫r四)%进行四(chi2)检验油a仨pre也.ct(阻，a;O)%预测一个值，这里一步预测，最后的一个值是不使用的，a后面可加任意数也atz=xhat:(end)%显示预测的一个值数学中国-2011年数学中国公益讲座一一数学建模之软件实现(司守奎