数学建模讲座一一计算软件的使用.pdf

数学中国 www madio n目 专业的数学建模论坛l 数学建模讲座一一计算软件的使用 司守奎 海军航空工程学院数学教研室 Em血1 sishoukui 163 m QQ 94339146 一般来说学习数学建模 常用的软件有五种 分别是Matlab Lingo Mathematica SPSS 和SAS 实际上使用什么软件都可以 关键是在论文中把数学原理写清楚 如果你搞清了 算法 你也可以自己编程从底层一步一步地实现算法 当然如果我们借助于Matl曲的强大 工具箱 可以起到事半功倍的作用 数学建模涉及的知识和算法特别多 我认为至少要掌握数学规划 图论 预测方法 评 价方法这四大部分知识 其它的数学知识和算法可能就需要数学建模竞赛时现学现实了 这 些知识我就不展开讲了 相关的数学建模书上都有这些内容 我新版的数学建模教材 数学 建模算法及应用 c国防工业出版社出版 第14章把评价方法总结了一下 第15章把预测 方法也总结了一下 新版的数学建模教材相对于网上流传的算法大全的电子文档 对一些内 容进行了压缩 增加了数字图像处理的知识 并增加了一些抽伽b工具箱的应用 对于数学规划的模型 建议大家使用Lingo软件求解比较方便 对于其它问题 如时间 序列模型 你使用什么软件求解都可以 关键看个人的喜好和对某种软件的熟悉程度 例如 你可以使用SPSS SAS或Eviews R软件等求解时间序列模型 我对Matla如软件比较熟悉 一般问题我都是使用陆tlab软件进行求解 由于时间的关系 我主要讲 rMat1ab软件的使用 x 4cos y 5 4sin s z 5 4sin sin 其中 0 2 1 画图的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 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叫昂 昂 的三维表面图 程序一 madio net r 可 叫x y 3 1 1 玄A勾 注意这里函数适用范围要广 对于向量也可逐点计算取值 I uadl仰 1 但U6计算工重积分 数学申锢 其中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申忡 其中的Z X2 向x2 y2 Z2日所围 成的区域 l 队y z x 句句 2 咀 2 y 2 求矩阵a的行数 目血u on国 血 1 a 1 2 构造待定参数作为未知数的线性方程组的系数阵 cs xishu a 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国r a 5 2 1 5 5 5 2 3 7 10 3 5 8 4 5 2 8 6 1 8 3 3 2 9 7 8 3 8 2011年数学中国公益讲座一一数学建模之软件实现 司守奎 量学中国 www 皿 专业的盘学噩摸桔坛l 7 3 4 5 2 6 4 7 由 size a l 回血IF 四届 m I 叫 1 2 构造线性方程组的系数阵 1rlog 叫 3日 计算线性方程组的常数项列 cs xishu b 拟舍1og y 中的未知参数 es IF呻 叫1 求函数y中的参数 计算的结果为a 2 6012 吨 0072 叫 0 0415 3 2 用户图形界面解法 Matl b中拟舍一元函数可以使用用户图形界面解法命令cftoo1 拟合二元函数可以使 用命奇sftoolo 例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 回 刷 伫斗NE Y 圄5韧蛐界面 曰缸 f 巴 1岳二 斗子 3 输入自变量X Y的数据x1 x2 困变量z的数据Z 并定义要拟舍的函数 0事exp al x exp a2句 设置完成以后的界面如下 2011年撞学中国告盐讲盛一一数学直模之软件实现 司守童 量学中国 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 邸 2 d山 氧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 飞也可以把计算结果输出到 Mat1a b工作空间中 3 3非线性拟合 Ma t1a b中非线性拟告的命守有lsqcurvefit和nlinfit 用户图形界面命 nlintool 下面通过一个例于说明lsqcurvefit的使用 例11利用表2给出的近两个世纪的美国人口统计数据 以百万为单位 建立人口预 测模型 最后用它预报20 0年美国的人口 囊2量国人口统计噩噩 1建模与求解 2011年撞学中国告盐讲盛一一数学直模之软件实现 司守童 数学中国 www madio n目 专业的数学建模论坛l 记x t 为第t年的人口数量 设人口年增长率r x 为x的线性函数 r x r sx 自 然资源与环境条件所能容纳的最大人口数为xm 即当X Xm时 增长率r xm 0 可得 r x r l 王 建立Logistic人口模型 xm 其解为 2 参数估计 ltt 卜巳生 纠勾 叮叫 r 1叫呻 1 dt 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 钮 cl 1 十 囚一 www maO 7a extre町 data2 txt 元1原始数据保存在纯文本文乎Fdata2 tit中 x a 2 2 61 提出人口数据 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 预测己知年代人口 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 把原始数据保存在纯文本文件也 ta4 txt中 x a 2 2 6 矿 x 二nonzeros x t 1790 10 20 1 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 20 1 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 294 386 r 0 0325 0 从上面的三种拟合方法可以看出 拟合同样的参数 方法不同可能结果相差很大 4 Mat lab图论工具箱 4 1 工具箱命令介绍 Mat1ab图论工具箱的命令见表30 表3Matlab固论工具箱的相关命令 命令名功能 graphal1sh世恒S陪她s求图中所有顶点对之间的最短距离 graphconnc础P找无向图的连通分支 或有向图的强 弱 连通分支 graphisdag 测试有向图是否含有圃 不含圄返回1 否则返回 graphisomorphism 确定两个图是否同构 同构返回1 否则返回 graphisspsn缸田确定一个图是否是生成树 是返回1 否则返回 graphmaxflow 计算有向图的最太流 2011年数学中国公益讲座一一数学建模之软件实现 司守奎 量学中国 www 皿 专业的盘学噩摸桔坛l Z叩h函皿 组回 g甲地m且与 叩 回胆th g叩k句 Z甲btravcnc 4 2 数据结构 在固中找量小生血树 把前驱顶点序列变成路径的顶点序列 求固中指定的 JIf顶点间的最短距圈和最且路径 执行有向无圃固的拓扑排序 求从一顶点出盘 所酷擅历固中的JJi点 Matlab图论工具箱使用的数据结构是稀疏矩阵 稀疏矩阵在数学上是指矩阵中零元素很多 非零元章很少的矩阵 对于稀疏矩阵 只 要存放非零元素的行标 到j标 非零元素的值即可 可以按如下方式存储 非零元素的行地址 非零元素的列地址 非零元素的值 在Matl曲中无向图和有向图邻接矩阵的使用上有很大差异 对于有向圈 只要写出邻接矩阵 直接使用抽血b的命 叩E田命 就可以把邻接 矩阵转化为稀疏矩阵的表示方式 对于无向图 由于邻接矩阵是对称阵 Matlab中只需使用邻接矩阵的下三角元素 即 Matlab只存储邻接矩阵下三角元素中的非零元素 稀疏矩阵只是一种存储格式 蛐咀曲中 普通矩阵使用事rso命令变成稀疏矩阵 稀 疏矩阵使用full命令变成普通矩阵 4 3应用举倒 下面给出图论工具箱在最短路 最小生成树和最大流中的应用例子 例12求圄6中从顶点1到顶点7的最短距离 数学中国 mp旦我 圄8元向圃 求解的抽d曲程序为 clc cl国X a zcros 邻接矩阵韧始化 下面首先输入邻接矩阵的上三角元素 叫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回国的取值为o h vicw iogra抖卫队 ShowArrows off ShowWcigl曲 on 画出该无向图 sot h Nod届恫也 Col时 10 40 4 下面是画出最短路径 参见Ma出b的帮助 fowEdgos g幽静sbynodci他got h No由咆础l ID 2011年撞学中国告盐讲盛一一数学直模之软件实现 司守童 数学中国 www madio n目 专业的数学建模论坛l revEdges 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阻 把稀疏矩阵a加载到Mat1ab工作壁间 T grapbminspantree a 调用求最小生成树的Mat1ab命令 例14图9的有向图中弧上的数字表示容量 求从V 到V 的最大流量 v 6 主4 4 2 2 V 4 v v3 5 7 亏 2 v2 4 V 固9有向固 5 随机模拟 蒙特卡洛方法也称为计算机随机模拟方法 它源于世界著名的赌城摩纳哥的Monte Carlo 蒙特卡洛 0它是基于对大量事件的统计结果来实现一些确定性问题的计算 使用蒙 特卡洛方法必须使用计算机生成相关分布的随机数 Matlab给出了生成各种随机数的命令 5 1 产生随机数的命令 1 随机数矩阵 Matlab中可以产生20多种分布的随机数 下面列举部分常用随机数的命令 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目 专业的数学建模论坛l poissmd mu m n 产生mXn矩阵 其中的元素是服从均值为m的泊松 Poisson 分布 的随机数 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 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川 wFa010 net 累积概率分布为0 4 0 7 1 如果产生的随机数在 0 O 4 上 认为事件A发生 产 生的随机数在 0 4 O 7 上 认为时间B发生 产生的随机数在 0 7 1J上 认为事件C发生 Matlab程序如下 c1c c1ear a 0 4 0 3 0 3J b cumsum a 求累积概率分布 c rand 产生一个 0 1 区间上均匀分布的随机数 d c b 1 2 b 1 c 计算落在区域V上的频数 V f1n