数学建模数据之简单处理技巧(Mathematica)

1、数学建模数据之简单处理技巧人们在生产实践与科学研究中经常会得到一系列的数据，然后通过这些数据得到某种内在规律，这就叫数据处理（Adjustment of Data）。科学家开发了许多方法来处理这个问题，最初由Gauss发展起来，用于彗星轨道（Orbits of Comets）的计算以及三角测量术中。主要方法有：最小二乘平方法、平均误差及误差延伸法则、直接测量的处理、以及一个函数用较简单函数表示的问题。数据拟合(Fit)就是其中的一种。假设已经得到数据列data1 = y1, y2, y3，yn, 现在需要寻找此数据列所满足的规律。Mathematica系统提供了拟合命令Fit，使用的格式如下，

2、例如：fx = Fit data1, 1, x, x2, x3 , x 表示用最小误差平方法去拟合数据data1，而且指明用构成的函数基，线性表出拟合函数fx。此处，得到的拟合函数fx 按x = j, f j = yj (data1中第j个数据)处理数据；一般地，假设有2维数据 data2 = x1, y1 , x2, y2 , , 则命令Fit data2, 1, f1x, f2x, , x 表示用最小误差平方法去拟合数据data2，而且指明用一元函数列 1, f1x, f2x, 去线性表出拟合函数Fx。假设有3维数据 data3 = x1, y1, z1 , x2, y2, z2 , ,

3、则命令fx, y = Fit data3, 1,f1x,y,f2x,y,x,y 表示用最小误差平方法去拟合数据data3，而且指明用2元函数列 1, f1x, y, f2x, y, 去线性表出拟合函数fx, y。数据拟合典型例子 d = 1, 1, 2, -2 , 3, 3 , 4, -4 , 5, 5 , 6, 6 ;g1 = ListPlot d, PlotStyle -> Hue 0 , PointSize .03 f1 = Fit d, 1, x, x2, x3, x4 , x ;Print“f1 = ”, f1g2 = Plot f1, x, 1, 10 , PlotStyle

4、 -> Hue .6 f2 = Fit d, 1, x, x2, x3, x4, x5, x ;Print“f2 = ”, f2g3 = Plot f2, x, 1, 10 , PlotStyle -> GrayLevel 0 , Dashing .03 Show g1, g2, g3 得到结果：图1-1-52图1-1-53图1-1-54图1-1-55例1 数据拟合示例。 要在山区修一条公路，测量员测得地形高程数据如下：平面区域的海拔（单位：百米）高程表X=0400800120016002000Y=037.04.705.506.006.706.904005.106.207.308.

5、008.508.708006.507.608.809.7010.2010.5012007.408.8010.8011.3012.5012.8016008.309.8011.8013.2014.5014.2020008.8010.6012.3013.9015.0015.00要从(0,0)到(2000,2000)修一条公路，问：如何选址？当必须经过某一指定点时，又如何选址？解：假设测点选择合理，即山坡表面是充分光滑的曲面，测点之间无山沟与山谷。第一步，将数据送到data，画图模拟：Mathematica命令如下：data = 省略 ;ListPlot3Ddata,ViewPoint->*,*

6、,*，AxesLabel->x,y,z。例如，可以选视点ViewPoint->-1,-1,2 或 -1,1,2画图做参考。第二步，画三角剖分面构成的拟合曲面。打开子程序包<<DiscreteMathComputationalGeometryTriangularSurfacePlotdata 表示为data画一个由三角剖分面构成的拟合曲面（图像略）。第三步，拟合。ff = Fit data, 1, x, x2, x3, x4, y,y2, y3, y4, x*y, x*y2, x*y3, x*y4,x2*y, x2*y2, x2*y3, x2*y4, x3*y, x3*y

7、2, x3*y3, x3*y4,x4*y, x4*y2, x4*y3, x4*y4, x, y 程序执行后，得到一个拟合函数ff。再画 ff 的图，效果更好。Plot3Dff,x,0,2000,y,0,2000（图像略）。例2 估计水箱流量（美国大学生MCM1991-A题）某些州的用水管理机构需估计公众用水速度（单位：加仑/小时）和每日总用水量的数据。许多地方没有测量流入或流出市政水箱流量的设备，而只能测量水箱中的水位（误差不超过0.5%）。当水箱水位某最低水位L时，水泵抽水灌入水箱，直到水位达到最高水位H为止。但是，也无法测量水泵的流量。因此，在水泵开动时，无法立即将水箱中的水位和用水量联系

8、起来。这种情形一天发生一次或两次，每次约为2小时。估计所有时刻，包括水泵抽水期间流出水箱的流量f(t)，并估计一天总用水量。下表给出某天某小镇的真实数据：时间（秒）033166635106191393717921212402522328543水 位（英尺）31.7531.1030.5429.9429.4728.9228.5027.9527.523228435923393323943543118466364995353936577546057426.97水泵开动水泵开动35.5034.4533.5032.6031.6730.8730.1264554685357185475020792548264

9、985948899539327029.2728.4227.6726.97水泵开动水泵开动34.7533.9733.40表中给出距开始测量的时间及即时水位。水箱是高40英尺、直径57英尺的正圆柱。通常水位低于L=27英尺时，水泵开始抽水，高于H=35.5英尺时，水泵停止工作。参考解答：（1） 由于水泵开动时，没有水位数据，所以，分段画图观察：Clear data data = 0,31.35,3316,31.10,6635,30.54,10619,29.94,13937,29.47,17921,28.92,21240,28.50,25223,27.95,28543,27.52,32284,26.

10、97,39435,35.50,43118,34.45,46636,33.50,49953,32.60,53936,31.67,57754,30.87,60574,30.12,64554,29.27,68535,28.42,71854,27.67,75020,26.97,85948,34.75,89953,33.97,93270,33.40 ;ListPlot data 得到结果：图1-1-56从图上看出，水位函数h(t)是分段线性函数。分段拟合如下：拟合后，得到结果：图1-1-57图1-1-58图1-1-59再画h(t) 的图与数据点图比较，将所画4张图放在一起，可见，拟合的非常好。图1-1-

11、60（2） 水箱中的水量 英尺3，流量 英尺3/秒于是：得到：f1 = 36.7628 英尺3/秒；f2 = 60.2972 英尺3/秒；f3 = 47.1388 英尺3/秒；一天的总用水量大约是：V = 36.7658*32284 + 60.2972*( 85948 39435 ) + 47.1388*( 93270 82649 )= 英尺3。（3） 再进一步就f1、f2、f3的数值，讨论管理机构应该如何操作（略）。数据处理建模实验例题演示例1台湾灾难性地震的数学模型及预测 台湾处于地震多发地带，其大大小小的地震每年都要发生十几次，但由于大多数地震发生在台湾的外海地区，给台湾造成的损失不大，

12、媒体报道，台湾地震专家认为“台湾中部、南部及北部在未来10年或数十年可能会相对平静，但花莲、台东地区发生大地震的可能性越来越高，已经进入地震警戒期”。资料显示，1900年以来给台湾带来灾难性损失的地震主要有以下几次：日期震级1904,11,66.51906,3,177.11935,4,117.11941,12,177.11964,1,186.51986,11,156.81999,9,217.62001,12,186.72002,3,307.5我们很想知道下一次台湾发生灾难性地震的日期与震级。特做以下数学模型给出预测，所有计算都用Mathematica软件包计算。1. 震级的预测 将震级看作是一

13、个时间序列：z = 6.5, 7.1, 7.1, 7.1, 6.5, 6.8, 7.6, 6.7, 7.5 。用Mathematica软件包画图观察，程序如下： z = 6.5, 7.1, 7.1, 7.1, 6.5, 6.8, 7.6, 6.7, 7.5 ; g1=ListPlotz,PlotJoined->True得到图形：可以看出图形有震荡，规律性不强。对数据做累加处理再画图： z = 6.5, 7.1, 7.1, 7.1, 6.5, 6.8, 7.6, 6.7, 7.5 ; zz = 6.5, 6.5+7.1, 6.5+7.1+7.1, 6.5+7.1+7.1+7.1, 6.5+

14、7.1+7.1+7.1+6.5, 6.5+7.1+7.1+7.1+6.5+6.8,6.5+7.1+7.1+7.1+6.5+6.8+7.6, 6.5+7.1+7.1+7.1+6.5+6.8+7.6+6.7, 6.5+7.1+7.1+7.1+6.5+6.8+7.6+6.7+7.5 ; g2=ListPlotzz,PlotJoined->True图形具有很强的直线性。于是，用线性函数模拟并画图如下： z = 6.5, 7.1, 7.1, 7.1, 6.5, 6.8, 7.6, 6.7, 7.5 ; zz = 6.5, 6.5+7.1, 6.5+7.1+7.1, 6.5+7.1+7.1+7.1,

15、 6.5+7.1+7.1+7.1+6.5, 6.5+7.1+7.1+7.1+6.5+6.8, 6.5+7.1+7.1+7.1+6.5+6.8+7.6, 6.5+7.1+7.1+7.1+6.5+6.8+7.6+6.7,6.5+7.1+7.1+7.1+6.5+6.8+7.6+6.7+7.5 ; g2=ListPlotzz,PlotJoined->True fz = Fitzz, 1, x, x; Print“fz = ”,fz g3=Plotfz,x,0,9，PlotStyle->RGBColor1,0,0 Showg2,g3得到zz的图像：得到拟合函数fz：fz = - 0.4694

16、44 + 7.005x 。拟合函数fz的图像：把这两个图放在一起：可见拟合得非常好。求出预测震级如下： u = fz/.x->10; u ( 6.5+7.1+7.1+7.1+6.5+6.8+7.6+6.7+7.5 )得到结果：6.68056即，台湾下一次灾难性地震的震级大约为 6.7 级。2. 地震日期的预测 首先，求出地震日期之间的间隔天数作为时间序列。打开Mathematica的日历子程序软件包： <<MiscellaneousCalendar然后计算： a = DaysBetween1904,11,6,1906,3,17, DaysBetween1906,3,17,19

17、35,4,11,DaysBetween1935,4,11,1941,12,17, DaysBetween1941,12,17,1964,1,18,DaysBetween1964,1,18,1986,11,15,DaysBetween1986,11,15,1999,9,21,DaysBetween1999,9,21,2001,12,18,DaysBetween2001,12,18,2002,3,30; Print“a = ”,a ListPlota,PlotJoined->True得到结果： a = 496,10617,2442,8067,8337,4693,819,102图形显示数据的规

18、律性也不强。因为数据都是正的，对数据做累加处理： a = 496,10617,2442,8067,8337,4693,819,102； b = 496,496+10617, 496+10617+2442, 496+10617+2442+8067, 496+10617+2442+8067+8337, 496+10617+2442+8067+8337+4693,496+10617+2442+8067+8337+4693+819, 496+10617+2442+8067+8337+4693+819+102; g1=ListPlotb,PlotJoined->True得到结果：图形较平缓的地方说

19、明地震间隔较短，对应灾害地震频发期，图形较陡峭的地方说明地震间隔较长，对应灾害地震休眠期。模型1 数据一次累加模型假设：台湾已经进入灾害地震休眠期在这种假设下，图形应该向陡峭发展。于是做以下拟合： 得到如下结果：58779.7这说明拟合函数ffx在x=9的值等于58779.7。再做计算 58779.7 (496+10617+2442+8067+8337+4693+819+102 ) = 23207 （天）即：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后23207天，台湾有可能发生6.7级地震。打开Mathematica软件包的日历功能： <<Miscellane

20、ousCalendar DaysPlus2002,3,30,23207得到结果： 2065,10,12 结论1：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后23207天，也就是2065年10月12日左右，台湾有可能发生6.7级地震。模型2 数据3次累加模型假设：台湾已经进入灾害地震休眠期在这种假设下，图形应该向陡峭发展。将数据做3次累加，然后用3次多项式做以下拟合。程序如下：得到结果如下：（1）3次累加数据b3的图像：（2）通过拟合得到的3次多项式及其图像：（3）把两个图像合并起来：可见，拟合的比较好。（4）数据还原计算：得到aa如下：于是，打开Mathematica软件包

21、的日历功能： <<MiscellaneousCalendar DaysPlus2002,3,30,4508得到结果： 2014,8,2 结论2：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后4508天，也就是2014年8月2日左右，台湾有可能发生6.7级地震。模型3 假设台湾任然是灾害地震频发期 在这种假设下，图形应该走向平缓。于是做以下拟合： <<StatisticsNonlinearFit 得到如下结果：35692.5119.532可以看出图形拟合得较好。另外，求出拟合函数yyx在x = 9的值等于35692.5。而且，在台湾任然是灾害地震频发期的

22、假设下， 2002年3月30日之后大约120( 35692.5 35573 = 119.5)天，台湾有可能发生6.7级地震。打开Mathematica软件包的日历功能： <<MiscellaneousCalendar DaysPlus2002,3,30,120得到结果： 2002,7,28 结论3：在台湾任然是灾害地震频发期的假设下， 2002年3月30日之后120天，也就是2002年7月28日左右，台湾有可能发生6.7级地震。总结如下：结论1：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后23207天，也就是2065年10月12日左右，台湾有可能发生6.7级地震

23、。结论2：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后4508天，也就是2014年8月2日左右，台湾有可能发生6.7级地震。结论3：在台湾任然是灾害地震频发期的假设下， 2002年3月30日之后120天，也就是2002年7月28日左右，台湾有可能发生6.7级地震。例2 艾滋病患者人数及时间预测。问题：某城市新闻媒体报道：医疗机构发现本城市有9人发作艾滋病，过了3个月，又报道发现11人发作艾滋病，再过3个月又报道有7人发作艾滋病。请预测一下此城市一共有多少人将发作艾滋病，大约需要多长时间？我们将尝试用不同的方法为这个问题建立模型。模型1 用微积分的方法建立模型。假设（1）以上

24、每批病人都在同一天发病，且每月按30天计算；（2）第t天发病人数为y(t)，且y(t)是连续函数；由上面假设，y(0) = 9, y(90) = 11, y(180) = 7。在Mathematica软件包中画图观察：data = 0, 9 , 90, 11 , 180, 7 ;ListPlot data, PlotJoined -> True 执行后得到图像如下：通过观察图像，我们进一步假设（3）函数y(t)是一个2次多项式函数。于是，我们通过数据拟合求函数y(t)。data = 0, 9 , 90, 11 , 180, 7 ;y = Fit data, 1, t, t2 , t ;P

25、rint “y = ”, y 执行后得到：y = 9 + 0.0555556t 0.00037037t2画一下y的图像： Plot y, t, -200, 300 得到图像如下：把这个函数图像与数据图像放在一起观察：求函数y(t)与坐标轴的两个交点，这两个交点分别代表有病人发作的最初时间和最后结束的时间。Solve y = 0, t 执行后得到： t -> - 97.9884 , t -> 247.988 。负的时间说明最初报道有9人发病之前就有人发病。这两个时间相减就得到所有病人发病总的持续时间。 247.988 ( - 97.9884 ) = 345.976 （天），即，大约3

26、46天。 总的发病病人数 = = 2556.38 （人）结论1：此城市大约有2557人将发病，发病持续时间大约346天，换句话说，最后一个病人发病距离最初报道9人发病大约248天。评注：这个模型简单易懂，但假设不一定合理，需要改进。改进后的模型如下：模型2 仍然用微积分的方法建立模型。（1）假设第t天发病人数为y(t)，且y(t)是连续函数。因为3个数据只能确定3个未知数，所以令函数y(t)是一个2次多项式函数：y(t) = a x2 +b x + c ；（2）假设报道的病人数都是在一段时间内发现的。即，从第一个病人发作到第一次报道时间3个月，共有9人发作，后3个月共有11人发作，再3个月共有7人发作。即于是在Mathematica软件包中编程如下：运行后得到结果：。求y与x轴的交点，再编程序如下：a = -