數(shù)學(xué)建模數(shù)據(jù)之簡單處理技巧(Mathematica)

1、數(shù)學(xué)建模數(shù)據(jù)之簡單處理技巧人們在生產(chǎn)實(shí)踐與科學(xué)研究中經(jīng)常會得到一系列的數(shù)據(jù)，然后通過這些數(shù)據(jù)得到某種內(nèi)在規(guī)律，這就叫數(shù)據(jù)處理（Adjustment of Data）。科學(xué)家開發(fā)了許多方法來處理這個(gè)問題，最初由Gauss發(fā)展起來，用于彗星軌道（Orbits of Comets）的計(jì)算以及三角測量術(shù)中。主要方法有：最小二乘平方法、平均誤差及誤差延伸法則、直接測量的處理、以及一個(gè)函數(shù)用較簡單函數(shù)表示的問題。數(shù)據(jù)擬合(Fit)就是其中的一種。假設(shè)已經(jīng)得到數(shù)據(jù)列data1 = y1, y2, y3，yn, 現(xiàn)在需要尋找此數(shù)據(jù)列所滿足的規(guī)律。Mathematica系統(tǒng)提供了擬合命令Fit，使用的格式如下，

2、例如：fx = Fit data1, 1, x, x2, x3 , x 表示用最小誤差平方法去擬合數(shù)據(jù)data1，而且指明用構(gòu)成的函數(shù)基，線性表出擬合函數(shù)fx。此處，得到的擬合函數(shù)fx 按x = j, f j = yj (data1中第j個(gè)數(shù)據(jù))處理數(shù)據(jù)；一般地，假設(shè)有2維數(shù)據(jù) data2 = x1, y1 , x2, y2 , , 則命令Fit data2, 1, f1x, f2x, , x 表示用最小誤差平方法去擬合數(shù)據(jù)data2，而且指明用一元函數(shù)列 1, f1x, f2x, 去線性表出擬合函數(shù)Fx。假設(shè)有3維數(shù)據(jù) data3 = x1, y1, z1 , x2, y2, z2 , ,

3、則命令fx, y = Fit data3, 1,f1x,y,f2x,y,x,y 表示用最小誤差平方法去擬合數(shù)據(jù)data3，而且指明用2元函數(shù)列 1, f1x, y, f2x, y, 去線性表出擬合函數(shù)fx, y。數(shù)據(jù)擬合典型例子 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 得到結(jié)果：圖1-1-52圖1-1-53圖1-1-54圖1-1-55例1 數(shù)據(jù)擬合示例。 要在山區(qū)修一條公路，測量員測得地形高程數(shù)據(jù)如下：平面區(qū)域的海拔（單位：百米）高程表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)修一條公路，問：如何選址？當(dāng)必須經(jīng)過某一指定點(diǎn)時(shí)，又如何選址？解：假設(shè)測點(diǎn)選擇合理，即山坡表面是充分光滑的曲面，測點(diǎn)之間無山溝與山谷。第一步，將數(shù)據(jù)送到data，畫圖模擬：Mathematica命令如下：data = 省略 ;ListPlot3Ddata,ViewPoint->*,*

6、,*，AxesLabel->x,y,z。例如，可以選視點(diǎn)ViewPoint->-1,-1,2 或 -1,1,2畫圖做參考。第二步，畫三角剖分面構(gòu)成的擬合曲面。打開子程序包<<DiscreteMathComputationalGeometryTriangularSurfacePlotdata 表示為data畫一個(gè)由三角剖分面構(gòu)成的擬合曲面（圖像略）。第三步，擬合。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 程序執(zhí)行后，得到一個(gè)擬合函數(shù)ff。再畫 ff 的圖，效果更好。Plot3Dff,x,0,2000,y,0,2000（圖像略）。例2 估計(jì)水箱流量（美國大學(xué)生MCM1991-A題）某些州的用水管理機(jī)構(gòu)需估計(jì)公眾用水速度（單位：加侖/小時(shí)）和每日總用水量的數(shù)據(jù)。許多地方?jīng)]有測量流入或流出市政水箱流量的設(shè)備，而只能測量水箱中的水位（誤差不超過0.5%）。當(dāng)水箱水位某最低水位L時(shí)，水泵抽水灌入水箱，直到水位達(dá)到最高水位H為止。但是，也無法測量水泵的流量。因此，在水泵開動(dòng)時(shí)，無法立即將水箱中的水位和用水量聯(lián)系

8、起來。這種情形一天發(fā)生一次或兩次，每次約為2小時(shí)。估計(jì)所有時(shí)刻，包括水泵抽水期間流出水箱的流量f(t)，并估計(jì)一天總用水量。下表給出某天某小鎮(zhèn)的真實(shí)數(shù)據(jù)：時(shí)間（秒）033166635106191393717921212402522328543水 位（英尺）31.7531.1030.5429.9429.4728.9228.5027.9527.523228435923393323943543118466364995353936577546057426.97水泵開動(dòng)水泵開動(dòng)35.5034.4533.5032.6031.6730.8730.1264554685357185475020792548264

9、985948899539327029.2728.4227.6726.97水泵開動(dòng)水泵開動(dòng)34.7533.9733.40表中給出距開始測量的時(shí)間及即時(shí)水位。水箱是高40英尺、直徑57英尺的正圓柱。通常水位低于L=27英尺時(shí)，水泵開始抽水，高于H=35.5英尺時(shí)，水泵停止工作。參考解答：（1） 由于水泵開動(dòng)時(shí)，沒有水位數(shù)據(jù)，所以，分段畫圖觀察：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 得到結(jié)果：圖1-1-56從圖上看出，水位函數(shù)h(t)是分段線性函數(shù)。分段擬合如下：擬合后，得到結(jié)果：圖1-1-57圖1-1-58圖1-1-59再畫h(t) 的圖與數(shù)據(jù)點(diǎn)圖比較，將所畫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） 再進(jìn)一步就f1、f2、f3的數(shù)值，討論管理機(jī)構(gòu)應(yīng)該如何操作（略）。數(shù)據(jù)處理建模實(shí)驗(yàn)例題演示例1臺灣災(zāi)難性地震的數(shù)學(xué)模型及預(yù)測 臺灣處于地震多發(fā)地帶，其大大小小的地震每年都要發(fā)生十幾次，但由于大多數(shù)地震發(fā)生在臺灣的外海地區(qū)，給臺灣造成的損失不大，

12、媒體報(bào)道，臺灣地震專家認(rèn)為“臺灣中部、南部及北部在未來10年或數(shù)十年可能會相對平靜，但花蓮、臺東地區(qū)發(fā)生大地震的可能性越來越高，已經(jīng)進(jìn)入地震警戒期”。資料顯示，1900年以來給臺灣帶來災(zāi)難性損失的地震主要有以下幾次：日期震級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我們很想知道下一次臺灣發(fā)生災(zāi)難性地震的日期與震級。特做以下數(shù)學(xué)模型給出預(yù)測，所有計(jì)算都用Mathematica軟件包計(jì)算。1. 震級的預(yù)測 將震級看作是一

13、個(gè)時(shí)間序列：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得到圖形：可以看出圖形有震蕩，規(guī)律性不強(qiáng)。對數(shù)據(jù)做累加處理再畫圖： 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圖形具有很強(qiáng)的直線性。于是，用線性函數(shù)模擬并畫圖如下： 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的圖像：得到擬合函數(shù)fz：fz = - 0.4694

16、44 + 7.005x 。擬合函數(shù)fz的圖像：把這兩個(gè)圖放在一起：可見擬合得非常好。求出預(yù)測震級如下： u = fz/.x->10; u ( 6.5+7.1+7.1+7.1+6.5+6.8+7.6+6.7+7.5 )得到結(jié)果：6.68056即，臺灣下一次災(zāi)難性地震的震級大約為 6.7 級。2. 地震日期的預(yù)測 首先，求出地震日期之間的間隔天數(shù)作為時(shí)間序列。打開Mathematica的日歷子程序軟件包： <<MiscellaneousCalendar然后計(jì)算： 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得到結(jié)果： a = 496,10617,2442,8067,8337,4693,819,102圖形顯示數(shù)據(jù)的規(guī)

18、律性也不強(qiáng)。因?yàn)閿?shù)據(jù)都是正的，對數(shù)據(jù)做累加處理： 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得到結(jié)果：圖形較平緩的地方說

19、明地震間隔較短，對應(yīng)災(zāi)害地震頻發(fā)期，圖形較陡峭的地方說明地震間隔較長，對應(yīng)災(zāi)害地震休眠期。模型1 數(shù)據(jù)一次累加模型假設(shè)：臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期在這種假設(shè)下，圖形應(yīng)該向陡峭發(fā)展。于是做以下擬合： 得到如下結(jié)果：58779.7這說明擬合函數(shù)ffx在x=9的值等于58779.7。再做計(jì)算 58779.7 (496+10617+2442+8067+8337+4693+819+102 ) = 23207 （天）即：在臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期的假設(shè)下， 2002年3月30日之后23207天，臺灣有可能發(fā)生6.7級地震。打開Mathematica軟件包的日歷功能： <<Miscellane

20、ousCalendar DaysPlus2002,3,30,23207得到結(jié)果： 2065,10,12 結(jié)論1：在臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期的假設(shè)下， 2002年3月30日之后23207天，也就是2065年10月12日左右，臺灣有可能發(fā)生6.7級地震。模型2 數(shù)據(jù)3次累加模型假設(shè)：臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期在這種假設(shè)下，圖形應(yīng)該向陡峭發(fā)展。將數(shù)據(jù)做3次累加，然后用3次多項(xiàng)式做以下擬合。程序如下：得到結(jié)果如下：（1）3次累加數(shù)據(jù)b3的圖像：（2）通過擬合得到的3次多項(xiàng)式及其圖像：（3）把兩個(gè)圖像合并起來：可見，擬合的比較好。（4）數(shù)據(jù)還原計(jì)算：得到aa如下：于是，打開Mathematica軟件包

21、的日歷功能： <<MiscellaneousCalendar DaysPlus2002,3,30,4508得到結(jié)果： 2014,8,2 結(jié)論2：在臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期的假設(shè)下， 2002年3月30日之后4508天，也就是2014年8月2日左右，臺灣有可能發(fā)生6.7級地震。模型3 假設(shè)臺灣任然是災(zāi)害地震頻發(fā)期 在這種假設(shè)下，圖形應(yīng)該走向平緩。于是做以下擬合： <<StatisticsNonlinearFit 得到如下結(jié)果：35692.5119.532可以看出圖形擬合得較好。另外，求出擬合函數(shù)yyx在x = 9的值等于35692.5。而且，在臺灣任然是災(zāi)害地震頻發(fā)期的

22、假設(shè)下， 2002年3月30日之后大約120( 35692.5 35573 = 119.5)天，臺灣有可能發(fā)生6.7級地震。打開Mathematica軟件包的日歷功能： <<MiscellaneousCalendar DaysPlus2002,3,30,120得到結(jié)果： 2002,7,28 結(jié)論3：在臺灣任然是災(zāi)害地震頻發(fā)期的假設(shè)下， 2002年3月30日之后120天，也就是2002年7月28日左右，臺灣有可能發(fā)生6.7級地震?？偨Y(jié)如下：結(jié)論1：在臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期的假設(shè)下， 2002年3月30日之后23207天，也就是2065年10月12日左右，臺灣有可能發(fā)生6.7級地震

23、。結(jié)論2：在臺灣已經(jīng)進(jìn)入災(zāi)害地震休眠期的假設(shè)下， 2002年3月30日之后4508天，也就是2014年8月2日左右，臺灣有可能發(fā)生6.7級地震。結(jié)論3：在臺灣任然是災(zāi)害地震頻發(fā)期的假設(shè)下， 2002年3月30日之后120天，也就是2002年7月28日左右，臺灣有可能發(fā)生6.7級地震。例2 艾滋病患者人數(shù)及時(shí)間預(yù)測。問題：某城市新聞媒體報(bào)道：醫(yī)療機(jī)構(gòu)發(fā)現(xiàn)本城市有9人發(fā)作艾滋病，過了3個(gè)月，又報(bào)道發(fā)現(xiàn)11人發(fā)作艾滋病，再過3個(gè)月又報(bào)道有7人發(fā)作艾滋病。請預(yù)測一下此城市一共有多少人將發(fā)作艾滋病，大約需要多長時(shí)間？我們將嘗試用不同的方法為這個(gè)問題建立模型。模型1 用微積分的方法建立模型。假設(shè)（1）以上

24、每批病人都在同一天發(fā)病，且每月按30天計(jì)算；（2）第t天發(fā)病人數(shù)為y(t)，且y(t)是連續(xù)函數(shù)；由上面假設(shè)，y(0) = 9, y(90) = 11, y(180) = 7。在Mathematica軟件包中畫圖觀察：data = 0, 9 , 90, 11 , 180, 7 ;ListPlot data, PlotJoined -> True 執(zhí)行后得到圖像如下：通過觀察圖像，我們進(jìn)一步假設(shè)（3）函數(shù)y(t)是一個(gè)2次多項(xiàng)式函數(shù)。于是，我們通過數(shù)據(jù)擬合求函數(shù)y(t)。data = 0, 9 , 90, 11 , 180, 7 ;y = Fit data, 1, t, t2 , t ;P

25、rint “y = ”, y 執(zhí)行后得到：y = 9 + 0.0555556t 0.00037037t2畫一下y的圖像： Plot y, t, -200, 300 得到圖像如下：把這個(gè)函數(shù)圖像與數(shù)據(jù)圖像放在一起觀察：求函數(shù)y(t)與坐標(biāo)軸的兩個(gè)交點(diǎn)，這兩個(gè)交點(diǎn)分別代表有病人發(fā)作的最初時(shí)間和最后結(jié)束的時(shí)間。Solve y = 0, t 執(zhí)行后得到： t -> - 97.9884 , t -> 247.988 。負(fù)的時(shí)間說明最初報(bào)道有9人發(fā)病之前就有人發(fā)病。這兩個(gè)時(shí)間相減就得到所有病人發(fā)病總的持續(xù)時(shí)間。 247.988 ( - 97.9884 ) = 345.976 （天），即，大約3

26、46天。 總的發(fā)病病人數(shù) = = 2556.38 （人）結(jié)論1：此城市大約有2557人將發(fā)病，發(fā)病持續(xù)時(shí)間大約346天，換句話說，最后一個(gè)病人發(fā)病距離最初報(bào)道9人發(fā)病大約248天。評注：這個(gè)模型簡單易懂，但假設(shè)不一定合理，需要改進(jìn)。改進(jìn)后的模型如下：模型2 仍然用微積分的方法建立模型。（1）假設(shè)第t天發(fā)病人數(shù)為y(t)，且y(t)是連續(xù)函數(shù)。因?yàn)?個(gè)數(shù)據(jù)只能確定3個(gè)未知數(shù)，所以令函數(shù)y(t)是一個(gè)2次多項(xiàng)式函數(shù)：y(t) = a x2 +b x + c ；（2）假設(shè)報(bào)道的病人數(shù)都是在一段時(shí)間內(nèi)發(fā)現(xiàn)的。即，從第一個(gè)病人發(fā)作到第一次報(bào)道時(shí)間3個(gè)月，共有9人發(fā)作，后3個(gè)月共有11人發(fā)作，再3個(gè)月共有7人發(fā)作。即于是在Mathematica軟件包中編程如下：運(yùn)行后得到結(jié)果：。求y與x軸的交點(diǎn)，再編程序如下：a = -