数学建模数据之简单处理技巧(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 - H

4、ue .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.008.508.7080

5、06.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-*,*,*，AxesLabel-x,y

6、,z。例如，可以选视点ViewPoint-1,-1,2 或 -1,1,2画图做参考。第二步，画三角剖分面构成的拟合曲面。打开子程序包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+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、+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, 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、.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.469444 + 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

9、 )得到结果：6.68056即，台湾下一次灾难性地震的震级大约为 6.7 级。2. 地震日期的预测 首先，求出地震日期之间的间隔天数作为时间序列。打开Mathematica的日历子程序软件包： True得到结果： a = 496,10617,2442,8067,8337,4693,819,102图形显示数据的规律性也不强。因为数据都是正的，对数据做累加处理： a = 496,10617,2442,8067,8337,4693,819,102； b = 496,496+10617, 496+10617+2442, 496+10617+2442+8067, 496+10617+2442+8067+

10、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得到结果：图形较平缓的地方说明地震间隔较短，对应灾害地震频发期，图形较陡峭的地方说明地震间隔较长，对应灾害地震休眠期。模型1 数据一次累加模型假设：台湾已经进入灾害地震休眠期在这种假设下，图形应该向陡峭发展。于是做以下拟合： 得到如下结果：58779.7这说明拟合函数ffx在x=9的值等于58779.7。再做计算 5877

11、9.7 (496+10617+2442+8067+8337+4693+819+102 ) = 23207 （天）即：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后23207天，台湾有可能发生6.7级地震。打开Mathematica软件包的日历功能： MiscellaneousCalendar DaysPlus2002,3,30,23207得到结果： 2065,10,12 结论1：在台湾已经进入灾害地震休眠期的假设下， 2002年3月30日之后23207天，也就是2065年10月12日左右，台湾有可能发生6.7级地震。模型2 数据3次累加模型假设：台湾已经进入灾害地震休眠期在这

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

13、在这种假设下，图形应该走向平缓。于是做以下拟合： StatisticsNonlinearFit 得到如下结果：35692.5119.532可以看出图形拟合得较好。另外，求出拟合函数yyx在x = 9的值等于35692.5。而且，在台湾任然是灾害地震频发期的假设下， 2002年3月30日之后大约120( 35692.5 35573 = 119.5)天，台湾有可能发生6.7级地震。打开Mathematica软件包的日历功能： True 执行后得到图像如下：通过观察图像，我们进一步假设（3）函数y(t)是一个2次多项式函数。于是，我们通过数据拟合求函数y(t)。data = 0, 9 , 90, 1

14、1 , 180, 7 ;y = Fit data, 1, t, t2 , t ;Print “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

15、( - 97.9884 ) = 345.976 （天），即，大约346天。 总的发病病人数 = = 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）假设报道的病人数都是在一段时间内发现的。即，从第一个病人发作到第一次

16、报道时间3个月，共有9人发作，后3个月共有11人发作，再3个月共有7人发作。即于是在Mathematica软件包中编程如下：运行后得到结果：。求y与x轴的交点，再编程序如下：a = - 1/243000; b = 2/2025; c = 1/15;y = a t2 + b t + c ;Solve y = 0, t 运行后得到结果：这两个根代表最初病人的发作时刻和最后病人的发作时刻，二者相减就是发病总时间：（天）。求发病的病人总数，编程如下：a = - 1/243000; b = 2/2025; c = 1/15;y = a t2 + b t + c ;运行后得到结果：29.3707（人）。结

17、论2：此城市大约有30人将发病，发病持续时间大约350天，换句话说，最后一个病人发病距离最初报道9人发病大约天。评注：这个模型简单易懂，但假设与结果都不一定合理，需要改进。进一步的模型如下：模型3 应用微分方程方法建立模型。（1）假设t时刻发病人数为y(t)，且y(t)是连续函数；（2）假设此城市共有K个人感染艾滋病毒，则在t时刻已经感染病毒但没有发作的人数是K y(t)。（3）发病速度与y(t)和K y(t)成正比，假设比例系数为a；（4）以上每批病人都在同一天发病，且每月按30天计算；由上面假设，y(0) = 9, y(90) = 11, y(180) = 7。依据假设建立模型如下：应用分离变量法求解此微分方程。令，得到方程组如下：从后两个方程解出x：于是，我们得到：解得：c = -1, c = 。将c = -1代入得到x = 1，推出a K = 0，不合题意，舍去。所