2011 司守奎（海军航空工程学院数学教研室） 数学中国公益讲座-计算软件的使用

数学中国(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 软件进行求解。 由于时间的关系，我主要讲~r Mat1ab 软件的使用。 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.皿.dio.net) 专业的盘学噩摸桔坛l 翼叫 co.(o:) 户 .~ 叫E拖到 <0. (13.). • - (5 +-4 叫哺剧。自】 2 γ ，。 ‘ 圄 l 例 1 的图形 对于其它的二次曲面，如果可以写成单支的显函数，直接使用命令.."，田h或 cz，田f画 图，杏则必须先化成参数方程. 例 2 绘制工元函数 z =SE叫昂，) 昂' 的三维面图. 程序一 . madio. net :r@(可，) ，叫x*y<呻，)/(x牛呻吟;%定义匿名函数 目四rf(z) %使用匿名函数画图 • a • • 圄 2 倒2 的图形 制 3 绘出椭圆柱面至+zi=l的图形· 9 4 解椭圆柱面的参数万程为 4 2011 年撞学中国告盐讲盛一一数学直模之软件实现〈司守童〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l lx=3cosα y = 2sinα， z=t， t ε R. 编写的抽.tIab 程序如下 x=司副alpha，t) 3*cos(a1pha); %定义参数方程的匿名函数 y=@(alpha,t) 2、in(alpha); z=@(alpha,t) t; ezml四h(x，y，z) 10 s O N -5 plot尔，y2，'*-'); plot(x,y3,'-o'); x=3 。因[叫， y= 28叫时. z= t 。 .1 y .2 -4 x o h=legend('sin($x$)'，'sin($x+\frac{气pi}{3} $)'，'∞s($础)') %1atex格式显示 set(h, 'Interpreter', '1a馆:x') xlabel('$x$'且阳preter','latex') %1au盟格式显示 y1abel('$y$'，'In阳preter','latex') %1au盟格式显示 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 4 量学中国 (www.皿.dio.net) 专业的盘学噩摸桔坛l 3 " 2 " 6到 "' o • 4 2 。 2 , , • 圃 4 制4 国固形 E 积分 Matlab 求符号积分的命令为 int. 求一重数值积分的命令为 quad. quadl，求二重数值 积分的命+为 dblqua.d.求三重数值积分的命奇为 triplequad. .2x+3 例 5 求 I一. . : d旨的数值积分. ,o l+x'" l'=@(X) (2*>:+3).1(1+玄 A勾;%注意这里函数适用范围要广，对于向量也可逐点计算取值 I"'([uadl仰，1) 但U6 计算工重积分 数学申锢-www.madio.net 其中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 & x.^2+y.^2+z.^2<=.刀， I=tripl吨uad(f，-'甲啊，町响，-町响1，'叭2)，0" 叭2)) 3 拟合 拟舍的准则多数情形下使用的是最小工乘法，即误差平方和最小z 一些经济学问题上 的拟舍使用最小一乘法，即误差绝对值和最小. 3.1 线性最小二乘姑 Mat1ab 中的线性最小二乘法实际上是解线性方程组，对应的抽血b命令是斜杠"飞.. 对于非齐次线性方程组Ax =b ， 抽血b 的求解命令为 x-A\b 无论Ax =b有唯一解，无穷多解还是无解时，陆咀曲总是给出一个解. (J)当 Ax =b唯一解时，陆.tl曲给出的就是该唯一解. 。〉当Ax =b有无穷多解时， Mat1ab 给出的是最小范数解. ω 当Ax =b无解时，陆.tl曲给出的是最小工乘解. 2011 年撞学中国告盐讲盛一一数学直模之软件实现〈司守童〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l ~J 8 利用表 1 的数据拟舍函数 y =ao +aj码 +a2x2 中的参数饨，~， a2 0 表 1 .l1 X2 y 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 7 3 4.5 2 6 4.7 拟含参数的Matlab 程序为 数学:中国-www.madio.net m=size(a,l); %求矩阵 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.皿.dio.net) 专业的盘学噩摸桔坛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响 中的参数吨，饨，a2 0 操作的步骤如下 (J)首先把数据加载到 Matlab 工作空间中g 运行如下程序 数掌:;中国-www.madio.net (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.皿.dio.net) 专业的盘学噩摸桔坛l 当F 唱f1 ，也时 1..",.., fi , I , .....' G::::二二一一" · 如 1 <<.，..) 二二k 勤，由， 。 由.. ",,,' "...." ,10< .. I C..,.. !.~旦旦曳 . "l 巾吃些~~) *exp.( :,1 • ,- '" • . 6 d I r;， o.，，_二 .. 1. 矿 ' , , " 固 6 参戴世置好嗣后曲界面 ~ ". .1 <2 (4) 用鼠标点一下界面上 "fit" 按钮，就得到拟合结果. l ih h' 1>.. 1.巾 '也町 1....咱 1; ，，"" ~邸， 2ωd 山..。 氧M川""..tJ..." .. 刷. .tJ ..,,', J .. ... o!o tl: ..,. :., ~回 % !i'阳 旦旦兰且ι」 坠坦坦旦旦 . ,t 委||LE 吉1 I t1Ú*exp (tl l*叶 *exp( :J ', ' 0."_ "' 固 7 拟告结果界面 c::;:;:!:] I J司 . 」 . net (5) 使用菜单 "fit" 中的子菜单 "Save to Workspace'" 飞也可以把计算结果输出到 Mat1a.b 工作空间中. 3.3 非线性拟合 Ma.t1a.b 中非线性拟告的命守有 lsqcurvefit 和 nlinfit，用户图形界面命~ nlintool . 下面通过一个例于说明 lsqcurvefi t 的使用. 例 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叫呻(ο1dt x(乌)=xo， x(t) = 一一二三旦 1+(盐 l)e 呻句) Xo (1) 把表 2 中的全部数据，包括后面的四个空格，全部保存到纯文本文件也.ta2.txt 中。 (1)非线性最小二乘估计 把表 2 中的第 1 个数据作为初始条件，利用余下的数据拟合式(1)中的参数xm和 r ， F看 I.Ml!tla]:l 担庄姐二f 1 .1 唱 E 钮， cl~ 1-十』囚一-www.maOl0.net"'-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 & x(k) = r-sx(k) , k = 2 ,3,…,22 , (3) 其中步长&=10 ，下面拟合其中的参数r和 s 0 编写如下的 Mat1ab 程序 clc, cl田r a=吃extr四d('也.ta4.txt'); %把原始数据保存在纯文本文件也.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) 咽 o.net 在=[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 表 3 Matlab 固论工具箱的相关命令 命令名 功能 graphal1sh世恒S陪她s 求图中所有顶点对之间的最短距离 graphconnc础P 找无向图的连通分支，或有向图的强(弱〉连通分支 graphisdag 测试有向图是否含有圃，不含圄返回1，否则返回。 graphisomorphism 确定两个图是否同构，同构返回 1 ，否则返回。 graphisspsn缸田 确定一个图是否是生成树，是返回 1，否则返回。 graphmaxflow 计算有向图的最太流 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 量学中国 (www.皿.dio.net) 专业的盘学噩摸桔坛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旦我45 adio. net 圄 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'， [1 00]) 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 & 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)上的连续函 数，数学期望 E[g(q)] = [g(x)!<(抽 ω 存在，则积分 (5) 式可用如下方法计算近似值。 设随机变量占的一系列取值为码，码，…，由 y， =g(x，) 形成的随机变量71 =g(幻的取 值数列为Y1'Y2'….显然由此产生的随机抽样y， 是相互独立的，且具有相同分布，若其数 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l 学期望为M ， 则根据太数定理有 叫|骂y;-M 什 =1 可见，只要n充分大，积分 (5) 式有近似值 E=;三g(码) 下面讨论用上述方法计算积分 J=[仰陆 的值。 为此，选择某种密度函数f(x)满足 [f(拙=1. 且能很方便地生成密度函数为f(x) 的随机抽样。同时将积分J写成如下形式 r h(x) r'-.'-L. r = I .:~.~ . f(x沛= I g(x)f(x)dx 品 f(x) 品 在很多情况下，往往取f(x) 为区间 (a， b)上均匀分布的概率密度函数 1-. _1_ , x E (a ,b), f(x)=~b-a 10, 其它. 计算平均值 后=;主h协于是积分 (7) 式的值可近刷刷 J 国 (b-a)最=主二旦艺h;(x) n 丁=1 例 16 计算积分 p. xdx I .JS-4x 解 随机模拟时取区间 (-1，1)上的均匀分布，其概率密度函数 I~ ， XE(-I,I) f(x) = ~ 2 随机模拟的地咀ab 程序如下 c1c, c1国r y=@(x) x.lsqrt(5-俨x);%定义被积函数 n=1000000; %生成随机数的个数 10, 其它 去=unifrnd(-I ， I ， [I ，n]); %生成 n个区间 (-1川上均匀分布的随机数 h=y(x); %计算被积函数的一系列取值 junzhi=sum(h)l同%计算取值的平均值 j江en=2叮unzhi %计算积分的近似值 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 (6) (7) 数学中国(www.madio.n目〉 专业的数学建模论坛 l 4 4 x 苟且fJ PBEJ rrJ 陆 h 算重计多分求积耍重设多假A斗r、" (8) 的值。积分区域Q是有界区域，被积函数f在区域Q 中是有界的。 设g(Xl'… ， x，，) 为区域Q上的概率密度函数，且当 f(x1 ， … ， X.) ;é 0 时 g(Xl'… ， x，，) 亦 不为霉。令 I f(x; … X.) •'""1 ，二一，当g(码，… ， X.) ;é 时h(Xl""'X.)=~ g(码，… ， x，，) lO，当g(叫，…， X.) = 0时. 则积分 (8) 可改写为 J= J- ..Jh帆，'" ,x.)g(x,,"', x.)战 dx. 。 若(XIf-- ， xn)是n维空间区域Q 中的随机变量，概率密度函数为 g(Xl'… ， X，，) ， 则随机变 量h(鸟，… ， X.) 的数学期望为 E[h(元""， X.)] = J… Jh帆，"'， x.)g(x1 ，… ， X，，)C:玩玩 =J. Q 即积分 J 是随机变量 h(X1 ，...，X) 的数学期望。如果选取N 个点 p， (X; ，…，X~)En (i = 1,.. '， N) 服从分布g(Xl' …，x，J ，则根据太数定理，其算术平均值 N J=士Lh(x; ，...，。 嗷攒i学学学毕毕学学毕轩近削哗啤{售但中』且l L一， (矶 ，...， X，，) EQ , g(码，… ， X，.) = 十v ~ '..1 阳， 其它 其中Y表示区域。的体积，则 h(x;,"', x.) = V. f(x1,… ,x ,,) 积分 (8) 的近似值可取为 TT N J=主Lf(儿 ， x!) ω 例 17 (续例7)计算 1= HJ(X+y+Z)2你今宿，其中 Q 为 z;:: x2 + y2 与 。 x2 + y2 +Z2 至 2 所围成的区域。 解 随机模拟时首先要计算Q 的体积，设Q 的体积为V ， 区域。上均匀分布的密度函 数为 |上 ， (x， y， z) ε V， f(x ,y ,z) = ~ V 10, 其它 求。的体积Y 时，在立体区域[-1，1] x [-1,1] x [0，.J2]上产生服从均匀分布的 10. 个随 机点，统计随机点落在Q的频数，则。的体积Y近似为上述立体的体积乘以频率。 仿真一的程序如下 clc, clear h=@(x品z) (x+y+z).^2; %定义被积函数 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l n~ 1 000000; %生成随机数的个数 X9血1量nd(-I ， I ，[1，叶); %生成 n个区间 (-1 ，0 上均匀分布的随机数 y=l皿1企nd(-1 ,1,[1,n]); z=uni量nd(0，sqrt(2)，[1，n]); f=sum(z>~.^2+y.^2 & x.^2+y.^2+z.^2<~2); %计算落在区域V上的频数 V~f1n*4.sqrt(2) %计算体积 恤c=h(x品z);%计算被积函数一系列的取值 jifen~V丰sum(hh.气~.^2+y.^2 & x.^2+y.^2+z.^2<~2))/f 仿真二的程序如下 c1c, c1国r h~@(x品z) (x+y+z).^2; %定义被积函数 n~ 1 000000; %生成随机数的个数 X9血1量nd(-sqrt(2)，sqrt(2)，[I ，n]); %生成 n个区间 (-1 ，0 上均匀分布的随机数 y=unifrnd(叫回(2)，sqrt(2)，[1具]); z=uni量nd(-sqrt(2),sqrt(2),[ 1,n]); f=sum(z>~.^2+y.^2 & x.^2+y.^2+z.^2<~2); %计算落在区域V上的频数 V~f1n*16 忖伊(2)%计算体积 恤c=h(x品z);%计算被积函数一系列的取值 jifen~V.sur叫hh卢(庐~.^2+y.^2 & x.^2+y.^2+z.^2<~))/f 自数据的处理 6.1 Matlab 中的默认数据文件 mat 文件 例 18 把抽tlab 工作空间中的数据矩阵 a， b , c 保存到数据文件 data1.mat 中。 save datal a b c 、 'út; M'l.ltlàb 白跑蝶、 q)(件 mat 文件可以省略后缀名。 唱. 息;;江'马卢王帮宫JJj阳川I~l，俨t~目成成m~t归时们4仨h恤"I咀ma田皿t 6.2 纯文本文件 可以把阳d文档中整行整列的数据粘贴到纯文本文件，然后调入到陆tlab工作空间中。 例 20 把纯文本文件 da匈2.txt 加载到工作空间。 a~load('data2.txt'); 或者是 a~extr叫E也饱2剧')， 例 21 使用 dlm阴阳命令把矩阵b保存到纯文本文件也.ta3.txt。 dlmwrite('da阳3.txt'，b) iJiJ 22 生成服从标准正态分布随机数的100x200 矩阵，然后用伽intf命令保存到纯文 本文件 data4.txt. 解 clc， c1阳 fid~fopen('data4.txt'，'w'); a=fiI町mrnd(0， 1 ， 100，200); 制ntt{fid，'%阳，a'); fclose( fid); 注:对于高维矩阵，用 dlmwrite 构造的纯文本文件， Lingo 软件不识别:为了 Lingo 软 件识别，纯文本文件必须用伽旭tf构造，而且数据之间的分割符为" 6.3 Ex阳cel 文件 iJiJ 23 把一个5xlO矩阵写到U Excel 文件 data5.由表单 Sheet2中 B2 开始的域中. a唱且邸，10); xlswri恒('da旭5.x1s飞a，'Sb田σ，'B2') 例 24 把例 23 生成的 Excel 文件也.ta5.x1s 中表单 Sheet2的域"C3:F6"中的数据赋给 b. b哩lsr四d('data5.x1s'，'Sh田σ，'C3:F6') 6.4 字符串数据 iJiJ 25 统计下列五行字符串中字符 a、 c、 E、 t 出现的频数. 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l l.aggcacg俱aa回cgg回ataacggaggaggacttgg国cggca伽cacggagg 2.cggaggacaaacgggatggcggta忧ggaggtggc部/lctgttcgggga 3.gggacgga饱cgga址ct邵阳acggacgg随aggagga凶cggcggaca:臼ca 4.a:堪回钮acggaaac田ac四.gac扭扭址cgg阳旦回a脑.cagaagc'忧a 5.cggctggcg伊G嗣cggactggcg伊1如回闹剧cggaggaggcggacggaggc 解把上述五行复制到一个纯文本数据文件 shuju.tx.t 中，编写如下程序 c1c fid=fopen('s由ht咱!ÎI且Ltx.t' while (←-f岛eot{f置'id)制〉沙) da饵=fgetl(fid); a=leng也(find(da饱.=97)); b=len脾(find(da饵，=99)); c=leng由(find(由.ta=103)); d=length(find( da也.=116)); e=leng也(find(da幅>=97&da饱.<=122)); t{i,:)=[a b c d e a+b+c+d]; i=1+1; end f, he=sum。 a皿write('pinshu.tx.t'，f); dlmwrite(恤shu.tx.t',he,' -append'); 也:lose叫fid); 其它的一些字符串处理命令有 strc皿Lp， S'园nd 等。 6.5 图像文件 iJiJ 26 把一个比较大的 bmp 图像文件 data6.b呻，转化成比较小的 lPg 文件，命名成 鸣、BPE WU辱示.-h日干r _______ ___ _ _1 者提写回~，蚓6届! 一-www.madio.net r 画ho叫纱 imwrite( a,'data7 .jpg'); figur毡， imshow('data7.jpg') 例 27 生成 10 幅彩色jpg 文件，依次命名成jpq1扣g， .."' jpgl0.jpg. c1c, c1ear f田 i=I:10 end str=[1晖，int2str(i)，' .jpg']; a:(:，:， I)咽硝500); a:(:，:，2)=r阻的则+100; a:(:，:，3)咽且d(500)+200; imwrite(a,str); 6 时间序列 时间序列的建模步骤为 (1)对给定序列做平稳性检验，将非平稳序列转换为平稳序列z (2) 对平稳时间序列X，做自相关函数 (ACF) 和偏自相关函数 (PACF) 图: (3) 根据 ACF 和 PACF 图对X，拟合 ARlMA (p,q ) 模型z (4) 估计模型残差吼: (5) 对残差估计值做检验，判断模型是否是适宜的，若不适宜，要进行修正。 (6) 利用建立的模型进行预测。 考虑芬兰 1949-2004 年的通货膨胀率，数据见表 5. 0.02 0.011 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 0.073 0.092 数学中国(www.madio.n目〉 专业的数学建模论坛 l 0.028 0.075 0.067 0.012 0.078 0.05 0.014 0.2 0.1 0.15 0.05 。 -0 .05 1斟0 1950 回回 四70 回回 回回 2000 2010 圄 10 芬兰 1949-2004 年通货膨胀率 t = 1，2，… ， n 的 Spe町四a相关系数q， ' 有 Vl--l-E(t-乓)2 , n(n" -1) "ti' (10) 构造统计量 T=q，~ - ~l-q; >1 1 ( 作下列假设检验 Ho: 序列X，平稳， H，: 序列X，非平稳(存在上升或下降趋势) . Daniel 检验方法z 对于显著水平α ，由时间序列X，计算 (t，乓)， t = 1，2，… ， n 的 Spearman 秩 相关系数乱，若ITI> 乌12(n-2) ， 则拒绝Ho' 认为序列非平稳。且当 q， >0 时，认为序 列有上升趋势 q， <0时，认为序列有下降趋势。又当 ITI 王 tal2 (n - 2) 时，接受Ho' 可以 认为码是平稳序列， 利用上面的 Daniel 检验方法，在 5%的水平下，认为序列是平稳的。做X， 的自相关函数 和偏自相关函数(见图 11，图 12)，可对X，拟合 ARC l)模型 X, = 0.8422x,_, + a" t = 1950,.. ',2004. 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l Sample Autocorr回国ion Function ' 。.8 6 ••• l~ ↓ 6420 00O ZM 帽甜臣 hnn 蜀，飞@在扣帽回 .{).2 .{).4 0 5 咽 o 15 Lag 圄 11 X， 自相关函数 Sample Partial Autoconelation Function 1 0.8 6420 00o mh 垣斟胆固 H昌，飞回誓，比@在 hh 由 -0.2 20 0.4 0 20 5 唱 o 15 Lag 圄 12 X，偏自相关函数 ARMA模型检验的x2检验 若拟合模型的残差记为乱，它是 Ii， 的估计。例如，对AR(p)序列，设未知参数的估计 是仇，仇，…，仇，则残差 2， =X，一仇X'_l 一-一份'pX，_p' t=1,2"" ,n , (设Xo =X_1 =...=X,_P =0) . 记 其中L为 6，自相关函数的拖尾数， Ljung-Box 的x2检验统计量是 L _2 /=n(n+2)Zj与 检验的假设是 Ho :Pk =0 ， 当 k 王 L; H， : 对某些k~L ， ρ'k "#。 在 Ho成立时，若样本容量n充分大， X2近似于x'{L-r)分布，其中 r是估计的模型参数 个数。 x2检验法z 给定显著性水平α ，查表得上α分位数x;(L-r) ，则当 x2 > 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 年数学中国公益讲座一一数学建模之软件实现(司守奎〉 数学中国(www.madio.n目〉 专业的数学建模论坛 l 仁=Qs.sqrt(n-2)/sqrt(1-Qs^2) %计算 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) %显示预测的一个值 数学中国-www.madio.net 2011 年数学中国公益讲座一一数学建模之软件实现(司守奎〉