全國(guó)大學(xué)生高教杯數(shù)學(xué)建模C題古塔變形論文

1、2013高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽承 諾 書(shū)我們仔細(xì)閱讀了全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽章程和全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽參賽規(guī)則（以下簡(jiǎn)稱(chēng)為“競(jìng)賽章程和參賽規(guī)則”，可從全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽網(wǎng)站下載）。我們完全明白，在競(jìng)賽開(kāi)始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢(xún)等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問(wèn)題。我們知道，抄襲別人的成果是違反競(jìng)賽章程和參賽規(guī)則的，如果引用別人的成果或其他公開(kāi)的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競(jìng)賽章程和參賽規(guī)則，以保證競(jìng)賽的公正、公平性。如有違反競(jìng)賽章程和參賽

2、規(guī)則的行為，我們將受到嚴(yán)肅處理。我們授權(quán)全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽組委會(huì)，可將我們的論文以任何形式進(jìn)行公開(kāi)展示（包括進(jìn)行網(wǎng)上公示，在書(shū)籍、期刊和其他媒體進(jìn)行正式或非正式發(fā)表等）。我們參賽選擇的題號(hào)是（從A/B/C/D中選擇一項(xiàng)填寫(xiě)）： C 我們的參賽報(bào)名號(hào)為（如果賽區(qū)設(shè)置報(bào)名號(hào)的話）： 所屬學(xué)校（請(qǐng)?zhí)顚?xiě)完整的全名）： 長(zhǎng)春工業(yè)大學(xué) 參賽隊(duì)員 (打印并簽名) ：1. 武太彬 2. 賈光輝 3. 牛文正 指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)： 李純凈 （論文紙質(zhì)版與電子版中的以上信息必須一致，只是電子版中無(wú)需簽名。以上內(nèi)容請(qǐng)仔細(xì)核對(duì)，提交后將不再允許做任何修改。如填寫(xiě)錯(cuò)誤，論文可能被取消評(píng)獎(jiǎng)資格

3、。） 日期： 2013 年 9 月_16_日賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：2013高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽編 號(hào) 專(zhuān) 用 頁(yè)賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：賽區(qū)評(píng)閱記錄（可供賽區(qū)評(píng)閱時(shí)使用）：評(píng)閱人評(píng)分備注全國(guó)統(tǒng)一編號(hào)（由賽區(qū)組委會(huì)送交全國(guó)前編號(hào)）：全國(guó)評(píng)閱編號(hào)（由全國(guó)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：古塔的變形摘 要本文對(duì)古塔的變形問(wèn)題建立數(shù)學(xué)模型，它實(shí)質(zhì)上是一個(gè)空間解析幾何問(wèn)題。首先建立空間解析幾何模型，并利用這個(gè)模型對(duì)問(wèn)題1進(jìn)行求解，然后對(duì)模型進(jìn)行數(shù)據(jù)處理和圖形分析，精確地給出了中心位置坐標(biāo)；之后對(duì)于問(wèn)題2，在問(wèn)題1的基礎(chǔ)上我們計(jì)算出中心坐標(biāo)的擬合直線，計(jì)算出曲率的

4、值；最后對(duì)于問(wèn)題3，運(yùn)用AR自回歸模型給出古塔的變形趨勢(shì)。在問(wèn)題1中，通過(guò)分析這四年古塔的每一層中的8 個(gè)離散點(diǎn)，運(yùn)用MATLAB數(shù)學(xué)軟件建立觀測(cè)數(shù)據(jù)模擬圖，依據(jù)此圖作出每一層相應(yīng)散點(diǎn)的投影平面，經(jīng)過(guò)計(jì)算基本近似正八邊形，得出正八邊形的中心，并運(yùn)用空間直線擬合模型，求出古塔的一條軸心擬合直線；并且運(yùn)用空間平面擬合模型，求出每一層的8個(gè)散點(diǎn)擬合的平面，那么直線與每一平面相交的點(diǎn)，即為每一層的中心坐標(biāo)，以所求第一層中心坐標(biāo)為例，1986年第一層中心坐標(biāo)為，1996年第一層中心坐標(biāo)為，2009年第一層中心坐標(biāo)為，2011年第一層中心坐標(biāo)，其它（算上塔尖）14層見(jiàn)正文表5。在問(wèn)題2中，首先運(yùn)用MATL

5、AB軟件畫(huà)出每一年的俯視圖，即各年的平面圖，可以看出這四年的古塔是逐漸傾斜和彎曲的，且在第五層開(kāi)始發(fā)生了扭曲。建立曲率數(shù)學(xué)模型，運(yùn)用軟件求解出相應(yīng)的曲率值。在問(wèn)題3中，首先建立帶季節(jié)項(xiàng)的AR自回歸模型，運(yùn)用問(wèn)題2中所求的曲率作為自變量，代入到模型中，運(yùn)用時(shí)間序列分析的統(tǒng)計(jì)軟件SPSS，得出古塔變形的趨勢(shì)，即隨年代的增長(zhǎng)，古塔的傾斜和彎曲將更加嚴(yán)重，且在第五層由于空間中心的改變，將更加扭曲。在本文最后，對(duì)模型的優(yōu)缺點(diǎn)及改進(jìn)之處進(jìn)行分析。關(guān)鍵詞：空間解析幾何 最小二乘法擬合 曲率 AR預(yù)測(cè)模型 MATLAB7.0軟件1問(wèn)題的重述與分析由于許多外在原因，古塔會(huì)產(chǎn)生各種變形，傾斜，彎曲，扭曲等。為了了

6、解各種變形量，測(cè)繪公司先后于1986年7月，1996年8月，2009年3月和2011年3月對(duì)該塔進(jìn)行了四次觀測(cè)。討論3個(gè)問(wèn)題，問(wèn)題1，給出確定古塔各層中心位置的通用方法，列表給出各次測(cè)量的古塔各層中心坐標(biāo)。本文首先對(duì)古塔的變形情況進(jìn)行分析，可以獲取位置信息，且只有四次的觀測(cè)數(shù)據(jù)信息。建立適當(dāng)?shù)淖鴺?biāo)系，先研究一層塔的八個(gè)點(diǎn)的投影所得的八個(gè)點(diǎn)的坐標(biāo)點(diǎn)，然后再確定各中心的坐標(biāo)，從而找到各層的中心點(diǎn)，可以擬合成一條直線。并且運(yùn)用空間擬合平面模型，求出每一層的8個(gè)散點(diǎn)擬合的平面，那么直線與每一平面相交的點(diǎn)，即為每一層的中心坐標(biāo)。 問(wèn)題2，分析該塔傾斜，彎曲，扭曲等變形情況。首先運(yùn)用MATLAB軟件畫(huà)出每

7、一年的俯視圖，即各年的平面圖，可以看出四年的古塔是傾斜還是彎曲或是扭曲，建立相應(yīng)的數(shù)學(xué)模型。問(wèn)題3，分析該塔的變形趨勢(shì)。首先建立帶季節(jié)項(xiàng)的AR自回歸模型，運(yùn)用問(wèn)題2中所求的結(jié)果作為自變量，代入到模型中，得出古塔變形的趨勢(shì)。2 基本假設(shè)1. 觀測(cè)的數(shù)據(jù)準(zhǔn)確；2. 以古塔的地面為水平面；3. 假設(shè)人在觀測(cè)時(shí)是圍繞塔的中心進(jìn)行八個(gè)不同的方面進(jìn)行大量的測(cè)量；4. 兩點(diǎn)的距離作到小數(shù)點(diǎn)后百分位即為等長(zhǎng)，后面忽略；5. 不考慮異常點(diǎn)；3 符號(hào)說(shuō)明已知每個(gè)點(diǎn)的橫坐標(biāo)；已知每個(gè)點(diǎn)的縱坐標(biāo)；已知每個(gè)點(diǎn)的豎坐標(biāo)；直線方程的待估系數(shù)；空間直線的對(duì)應(yīng)的方向向量；：分別是邊和邊由其原始位置旋轉(zhuǎn)的角度；: 變化率；4 模

8、型的建立與求解4.1 問(wèn)題1的求解 首先根據(jù)已知每一年數(shù)據(jù)中的每一層的8個(gè)散點(diǎn)坐標(biāo)運(yùn)用MATLAB軟件，得到如下每年觀測(cè)數(shù)據(jù)的模擬圖。圖1 1986年觀測(cè)數(shù)據(jù)模擬圖圖2 1996年觀測(cè)數(shù)據(jù)模擬圖圖3 2009年觀測(cè)數(shù)據(jù)模擬圖圖4 2011年觀測(cè)數(shù)據(jù)模擬圖由圖1到圖4的模擬圖可以直觀地看到，四次測(cè)量古塔的顯著變化，相同坐標(biāo)上塔的上部越來(lái)越小，把圖2和圖3疊加到一個(gè)圖。圖5 古塔扭曲變化對(duì)比圖圖5中，紅線代表1996年的塔圖，黃線代表2009年的塔圖。兩個(gè)塔圖疊加在一起對(duì)比更加顯明，而在第五層出現(xiàn)扭曲，從第五層往上越來(lái)越小，到2011年成了一個(gè)點(diǎn)的塔尖。 經(jīng)過(guò)以上的圖形分析，本文以空間解析幾何為背

9、景，運(yùn)用最小二乘擬合理論給出空間直線與平面的擬合。下面建立空間內(nèi)曲線的最小二乘問(wèn)題，以空間直線為例研究曲線擬合的方法。 已知空間直線的標(biāo)準(zhǔn)方程整理得直線射影式方程 ，其中;這樣直線可以看作是用這 2 個(gè)方程表示的平面相交的直線，所以可以分別對(duì) 2 個(gè)方程進(jìn)行數(shù)據(jù)擬合。 設(shè)表示按擬合方程求得的近似值。一般地，它不同于實(shí)測(cè)值兩者之差，同理可得 當(dāng)Q取最小值時(shí)的值即為方程的系數(shù)，即滿(mǎn)足下列方程時(shí)Q值最小有令 (1)方程組(1)可寫(xiě)成 其中根據(jù)m組數(shù)據(jù)點(diǎn)解方程組就可以求得的值。圖6 1996年古塔各層中心點(diǎn)擬合直線分析圖圖7 2009年古塔各層中心點(diǎn)擬合函數(shù)分析圖從圖6，圖7中可以得出：由1996年古

10、塔各層中心點(diǎn)擬合直線與2009年古塔各層中心點(diǎn)擬合直線，方程為，由一次方程的系數(shù)k看出1996年為0.0105,2009年為0.0112，可以得出古塔的中心發(fā)生了一次傾斜，且2009年比1986年傾斜角度變大。運(yùn)用附件中MATLAB的程序經(jīng)過(guò)調(diào)試，得到如下結(jié)果： 1986年 1996年 2009年 2011年接下來(lái)運(yùn)用最小二乘法擬合求空間平面方程：設(shè)空間平面方程為，其中為待估參數(shù)。設(shè)古塔的每層8個(gè)點(diǎn)，設(shè)點(diǎn)坐標(biāo)，考慮到數(shù)據(jù)在三個(gè)方向均存在誤差，得矩陣形式：即。通常采用矩陣奇異值分解解算待定參數(shù)的整體最小二乘解。其中則參數(shù)的整體最小二乘估計(jì)為：殘差矩陣為：。運(yùn)用MATLAB軟件得到古塔各層平面方程

11、，列表為：表1 1986年的擬合平面方程z=(0.003417)*x+(-0.000831)*y+(0.471956) z=(0.003629)*x+(-0.000818)*y+(5.887189) z=(0.003726)*x+(-0.000883)*y+(11.308475) z=(0.003838)*x+(-0.000935)*y+(15.602646) z=(0.003941)*x+(-0.000874)*y+(20.156847) z=(0.005443)*x+(-0.017498)*y+(33.310426) z=(0.005668)*x+(-0.018744)*y+(37.502

12、066) z=(0.005887)*x+(-0.020009)*y+(41.620122) z=(0.006115)*x+(-0.021455)*y+(45.825845) z=(0.000887)*x+(-0.021682)*y+(52.003938) z=(0.001868)*x+(-0.022189)*y+(56.048295) z=(0.002175)*x+(-0.023884)*y+(61.121714) z=(0.004664)*x+(-0.027610)*y+(66.050359)表2 1996年的擬合平面方程z=(0.003715)*x+(-0.001488)*y+(0.684

13、400)z=(0.003782)*x+(-0.000493)*y+(5.617203)z=(0.003735)*x+(-0.001266)*y+(11.515906)z=(0.004045)*x+(0.000715)*y+(14.555943)z=(0.003941)*x+(-0.001286)*y+(20.385763)z=(0.005639)*x+(-0.017135)*y+(32.997144)z=(0.005880)*x+(-0.018358)*y+(37.168214)z=(0.005877)*x+(-0.020489)*y+(41.891314)z=(0.006366)*x+(-0

14、.021014)*y+(45.437916)z=(0.000873)*x+(-0.022225)*y+(52.314084)z=(0.001853)*x+(-0.022770)*y+(56.380505)z=(0.002159)*x+(-0.024511)*y+(61.481252)z=(0.005234)*x+(-0.027989)*y+(65.963401)表3 2009年的擬合平面方程z=(-0.003770)*x+(-0.002373)*y+(5.079985)z=(-0.000194)*x+(-0.003972)*y+(9.661222)z=(-0.003697)*x+(-0.002

15、317)*y+(15.977842)z=(-0.000065)*x+(-0.004433)*y+(19.616420)z=(-0.000946)*x+(-0.004248)*y+(24.611974)z=(-0.020936)*x+(-0.004708)*y+(39.819633)z=(-0.018595)*x+(-0.006811)*y+(43.402884)z=(-0.021604)*x+(-0.006530)*y+(48.330004)z=(-0.022140)*x+(-0.006879)*y+(52.311658)z=(-0.022761)*x+(-0.001523)*y+(52.91

16、5196)z=(-0.022023)*x+(-0.002515)*y+(57.363039)z=(-0.023948)*x+(-0.002804)*y+(62.798782)z=(-0.027663)*x+(-0.005448)*y+(70.356472)表4 2011年的擬合平面方程z=(-0.003740)*x+(-0.002333)*y+(5.040576)z=(-0.002530)*x+(-0.002555)*y+(10.061014)z=(-0.004408)*x+(-0.002633)*y+(16.522721)z=(-0.002736)*x+(-0.003593)*y+(20.5

17、18390)z=(-0.001351)*x+(-0.004259)*y+(24.823896)z=(-0.021377)*x+(-0.004712)*y+(40.045928)z=(-0.018205)*x+(-0.007014)*y+(43.305695)z=(-0.021560)*x+(-0.006476)*y+(48.273096)z=(-0.025300)*x+(-0.004743)*y+(52.729697)z=(-0.026175)*x+(0.000449)*y+(53.562983)z=(-0.024393)*x+(-0.002167)*y+(58.395736)z=(-0.01

18、8681)*x+(-0.004368)*y+(60.918705)z=(-0.028368)*x+(-0.005161)*y+(70.556460)聯(lián)立直線與平面方程，即得到表5的各年各層的中心坐標(biāo)。表5 各年各層的中心坐標(biāo)1986年1996年2009年2011年1層2層3層4層5層6層7層8層9層10層11層12層13層4.2 問(wèn)題2的求解 問(wèn)題2是為了分析古塔的變形情況，根據(jù)問(wèn)題1的圖1-圖5，運(yùn)用MATLAB軟件，得到如下關(guān)于軸的俯視圖：圖8 古塔模擬俯視圖圖9 1986年觀測(cè)數(shù)據(jù)俯視模擬圖圖10 1996年觀測(cè)數(shù)據(jù)模擬俯視圖圖11 2009年觀測(cè)數(shù)據(jù)模擬俯視圖 圖8至圖11分別給出了整

19、體和各個(gè)年份的俯視圖，從四個(gè)圖中可以看到并不是均勻的八邊形環(huán)，有些邊環(huán)密，有些邊環(huán)稀，這就直觀地說(shuō)明了古塔發(fā)生了傾斜和彎曲，加上問(wèn)題1的第五層的中心坐標(biāo)的改變，說(shuō)明古塔在第五層上也發(fā)生了扭曲現(xiàn)象。本文在第一問(wèn)最后得出每層各中心坐標(biāo)的擬合直線，運(yùn)用到第二問(wèn)中古塔的傾斜問(wèn)題上，傾斜是指基礎(chǔ)兩端點(diǎn)傾斜方向的沉降差與其距離的比值建筑中心線或其墻、柱等，在不同高度的點(diǎn)對(duì)其相應(yīng)底部點(diǎn)的偏移現(xiàn)象。運(yùn)用公式兩條直線夾角公式：彎曲，即不直?？梢苑譃樾巫儚澢翱臻g彎曲。當(dāng)桿件受到與桿軸線垂直的外力或在軸線平面內(nèi)的力偶作用時(shí)，桿的軸線由原來(lái)的直線變成曲線，這種變形叫彎曲變形。曲率處處不為零的空間稱(chēng)為彎曲空間。物體發(fā)

20、生彎曲時(shí)產(chǎn)生的形變叫做“彎曲形變”。物體彎曲得越厲害，產(chǎn)生的彈力就越大。例如，將弓拉得越滿(mǎn)，箭就射得越遠(yuǎn)。把一個(gè)物體放在支持物上，物體越重，支持物被壓彎曲得越厲害，支持力就越大。運(yùn)用曲率公式：扭曲物體因外力作用而扭轉(zhuǎn)變形,也用于比喻數(shù)學(xué)的一種改變物體形狀的方法，扭曲變形是在彎曲變形的基礎(chǔ)上，旋轉(zhuǎn)，扭成螺旋狀態(tài)，叫扭曲變形。問(wèn)題1中得出扭曲的發(fā)生主要在第五層。首先,計(jì)算方向的扭曲變形。由于邊和的旋轉(zhuǎn)角度一般都很小(對(duì)于一個(gè)實(shí)際問(wèn)題,一般因此可以認(rèn)為約等于,因此:式中:分別是邊和邊由其原始位置旋轉(zhuǎn)的角度。如果，則說(shuō)明微元的一邊產(chǎn)生均勻豎向位移,則變形后各點(diǎn)仍保持在一個(gè)平面內(nèi),意味著這個(gè)微元發(fā)生的是

21、剛性的旋轉(zhuǎn)。否則,當(dāng),即在兩條邊之間存在一個(gè)旋轉(zhuǎn)角的差，就說(shuō)明產(chǎn)生了扭曲變形。則在變形后的面上,扭曲變形的大小可以采用這兩個(gè)角度之差隨著兩個(gè)對(duì)邊之間距離的變化率來(lái)表示:或者:如果應(yīng)用以上分析過(guò)程考慮y方向的扭曲變形,即首先計(jì)算微元另外兩條對(duì)邊和的旋轉(zhuǎn)角度的差,將會(huì)得到相同的結(jié)果。最終運(yùn)用MATLAB軟件得出4.3 問(wèn)題3的求解在問(wèn)題1的基礎(chǔ)上，運(yùn)用MATLAB軟件，給出1996年與2009年的中心擬合直線圖來(lái)。時(shí)間序列分析的目標(biāo)就是通過(guò)分析要素（變量）隨時(shí)間變化的歷史過(guò)程，揭示其變化發(fā)展規(guī)律，并對(duì)未來(lái)狀態(tài)進(jìn)行分析預(yù)測(cè)。如在變形測(cè)量中，可以采用時(shí)間序列分析方法對(duì)觀測(cè)數(shù)據(jù)進(jìn)行分析，以便建立變形體的

22、動(dòng)態(tài)變形預(yù)測(cè)模型，并對(duì)其變形趨勢(shì)進(jìn)行預(yù)測(cè)。而自回歸( AR) 模型的參數(shù)估計(jì)是時(shí)間序列分析的基本問(wèn)題，是在模型結(jié)構(gòu)及階次已確定的條件下，對(duì)模型參數(shù)進(jìn)行估計(jì)，使所建立的模型是實(shí)際時(shí)間序列的“最佳”擬合模型。但在實(shí)際的觀測(cè)中，觀測(cè)值是由一定觀測(cè)手段得到的，不可避免地含有隨機(jī)誤差。理論上講AR( p) 參數(shù)的最小二乘估計(jì)也就是線性最小二乘估計(jì)，所以可以用現(xiàn)在廣泛應(yīng)用的整體最小二乘估計(jì)方法進(jìn)行估計(jì)，即不僅考慮自身觀測(cè)值的誤差，同時(shí)考慮與其有關(guān)的自身前一個(gè)或前幾個(gè)時(shí)刻的觀測(cè)值的誤差，從而進(jìn)行參數(shù)估計(jì)，得到的預(yù)測(cè)觀測(cè)數(shù)據(jù)更準(zhǔn)確。但AR 模型是典型的預(yù)測(cè)模型，而最小二乘估計(jì)是常用的數(shù)據(jù)分析工具，在預(yù)測(cè)數(shù)據(jù)上

23、還是常應(yīng)用最小二乘估計(jì)比較準(zhǔn)確，反而最小二乘估計(jì)的數(shù)據(jù)不夠準(zhǔn)確和穩(wěn)定。自回歸模型，子樣觀察值，白噪聲序列表示為，回歸系數(shù)用表示，則可以得到模型： ，模型參數(shù)的最小二乘估計(jì)。設(shè)樣本觀測(cè)值，記,則模型可表示為有最小二乘原理可得到模型參數(shù)的估計(jì)為。那么根據(jù)最小二乘估計(jì)值可以得到噪聲的估計(jì)值為噪聲方差為，由此模型可以求出未來(lái)的變形趨勢(shì)。5模型的評(píng)價(jià)及推廣5.1模型的優(yōu)點(diǎn)(1). 在問(wèn)題的求解中，充分運(yùn)用了數(shù)據(jù)和圖形，即表格，使結(jié)果明了清晰。(2). 本文采用了多種專(zhuān)業(yè)軟件對(duì)模型進(jìn)行求解，如Matlab，SPSS等提高了模型 的精確度。(3). 本文所有模型建立均完全基于實(shí)際的統(tǒng)計(jì)數(shù)據(jù)，擬合成近似的線進(jìn)

24、行處理， 科學(xué)合理。(4). 本文建立多種模型，多種模型的對(duì)比更能體現(xiàn)不同模型的特點(diǎn)及優(yōu)勢(shì)。(5). 應(yīng)用相關(guān)幾何知識(shí)，將空間復(fù)雜的點(diǎn)分布簡(jiǎn)化成直線和面，應(yīng)用幾何問(wèn)題求 空間點(diǎn)。5.2模型的缺點(diǎn)(1). 背景資料的篩選方法有待進(jìn)一步優(yōu)化和改進(jìn)；(2). 模型建立中采用的數(shù)據(jù)不太準(zhǔn)確，利用擬合估計(jì)值，對(duì)模型產(chǎn)生了相對(duì)較大 的誤差。(3). 建立的模型不夠完善，忽略了許多的因素，過(guò)于簡(jiǎn)單。(4). 由于古塔在現(xiàn)實(shí)生活中受到各種影響，理論值難免與實(shí)際情況有所偏差。5.3模型的改進(jìn) 由于都是求出每個(gè)點(diǎn)都是一個(gè)估計(jì)值，從而有較大的誤差，需要在模型的建立中進(jìn)行改進(jìn)，使估計(jì)的所有的數(shù)據(jù)更加準(zhǔn)確。 5.4模型

25、的推廣用中心位置的模型方法，不難預(yù)測(cè)全國(guó)所有古塔或高樓由于外力和內(nèi)力的作用在以后的變化趨勢(shì)。6.參考文獻(xiàn)1 姚澤清,鄭旭東,全國(guó)大學(xué)生數(shù)學(xué)建模賽題優(yōu)秀論文評(píng)析,北京:國(guó)防工業(yè)出版社，2012: 198-218.2 周培德，計(jì)算幾何M. 北京：清華大學(xué)出版社，2011：100-105.3 襲楊，空間直線擬合的一種方法J.齊齊哈爾大學(xué)學(xué)報(bào)，2009,25(2):64-68.4 王能超. 數(shù)值分析簡(jiǎn)明教程M. 北京：高等教育出版社，2003：36-37. 5 呂林根，許子道. 解析幾何M. 北京：高等教育出版社，2004：156-157.6 于曉秋，李欣，張宏禮. MATLAB 數(shù)值實(shí)驗(yàn)M. 沈陽(yáng)：

26、遼寧科學(xué)技術(shù)出版社，2004：34-40. 7 王正東. 數(shù)學(xué)軟件與數(shù)學(xué)實(shí)驗(yàn)M. 北京：科學(xué)出版社，2004：88-89. 8 吳懷宇 時(shí)間序列分析與綜合M 武漢: 武漢大學(xué)出版社， 20049 俞錦成 關(guān)于整體最小二乘的可解性J 南京師范大學(xué)學(xué)報(bào)( 自然科學(xué)版) ， 1996， 19( 1) : 13-167 附錄直線擬合： x=566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9843 567.0218 567.0569 567.1045 567.1518; y=522.7105 522.6684 522.

27、6273 522.5944 522.5591 522.5244 522.5081 522.4924 522.4764 522.4624 522.423 522.3836; z=1.7874 7.3202 12.7552 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A = 0.0104 566.6411B = -0.0066 522.7125 /1986年 x=566.

28、665 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.062 567.1102 567.1578; y=522.7102 522.6674 522.6256 522.5922 522.5563 522.521 522.5042 522.4881 522.4714 522.4572 522.4173 522.3775; z=1.783 7.3146 12.7507 17.0751 21.716 26.2295 29.8323 33.3454 36.8483 40.1676 44.4354 4

29、8.7074; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A = 0.0105 566.6411B = -0.0067 522.7124 /1996年 x=566.727 566.7642 566.8004 566.8297 566.861 566.9478 566.98 567.0313 567.0825 567.1381 567.181 567.2238; y=522.7014 522.669 522.6387 522.6127 522.586 522.5335 522.5115 522.4788 522

30、.4457 522.3926 522.3535 522.3147; z=1.7632 7.2905 12.7269 17.052 21.7039 26.2045 29.817 33.3366 36.8222 40.1441 44.4249 48.6839; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A =0.0112 566.6643B =-0.0084 522.7436 /2011年 x=566.7268 566.764 566.8001 566.8293 566.8603 566.9471 566.979

31、2 567.0305 567.0816 567.137 567.1799 567.2225; y=522.7015 522.6693 522.6384 522.6132 522.5866 522.5342 522.5123 522.4797 522.4466 522.3937 522.3547 522.316; z=1.7645 7.309 12.7323 17.0697 21.7094 26.211 29.8246 33.3399 36.8438 40.1611 44.4326 48.6998; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y

32、; A=(MN) B=(MO) A =0.0112 566.6642B = -0.0084 522.7436 /2009年平面擬合公式：86.1 z=(0.003417)*x+(-0.000831)*y+(0.471956) z=(0.003629)*x+(-0.000818)*y+(5.887189) z=(0.003726)*x+(-0.000883)*y+(11.308475) z=(0.003838)*x+(-0.000935)*y+(15.602646) z=(0.003941)*x+(-0.000874)*y+(20.156847) z=(0.005443)*x+(-0.01749

33、8)*y+(33.310426) z=(0.005668)*x+(-0.018744)*y+(37.502066) z=(0.005887)*x+(-0.020009)*y+(41.620122) z=(0.006115)*x+(-0.021455)*y+(45.825845) z=(0.000887)*x+(-0.021682)*y+(52.003938) z=(0.001868)*x+(-0.022189)*y+(56.048295) z=(0.002175)*x+(-0.023884)*y+(61.121714) z=(0.004664)*x+(-0.027610)*y+(66.0503

34、59)96.1 z=(0.003715)*x+(-0.001488)*y+(0.684400) z=(0.003782)*x+(-0.000493)*y+(5.617203)z=(0.003735)*x+(-0.001266)*y+(11.515906) z=(0.004045)*x+(0.000715)*y+(14.555943) z=(0.003941)*x+(-0.001286)*y+(20.385763) z=(0.005639)*x+(-0.017135)*y+(32.997144) z=(0.005880)*x+(-0.018358)*y+(37.168214) z=(0.0058

35、77)*x+(-0.020489)*y+(41.891314) z=(0.006366)*x+(-0.021014)*y+(45.437916)z=(0.000873)*x+(-0.022225)*y+(52.314084) z=(0.001853)*x+(-0.022770)*y+(56.380505) z=(0.002159)*x+(-0.024511)*y+(61.481252) z=(0.005234)*x+(-0.027989)*y+(65.963401)99.1 z=(-0.003770)*x+(-0.002373)*y+(5.079985) z=(-0.000194)*x+(-0

36、.003972)*y+(9.661222) z=(-0.003697)*x+(-0.002317)*y+(15.977842) z=(-0.000065)*x+(-0.004433)*y+(19.616420) z=(-0.000946)*x+(-0.004248)*y+(24.611974) z=(-0.020936)*x+(-0.004708)*y+(39.819633) z=(-0.018595)*x+(-0.006811)*y+(43.402884) z=(-0.021604)*x+(-0.006530)*y+(48.330004) z=(-0.022140)*x+(-0.006879

37、)*y+(52.311658) z=(-0.022761)*x+(-0.001523)*y+(52.915196) z=(-0.022023)*x+(-0.002515)*y+(57.363039) z=(-0.023948)*x+(-0.002804)*y+(62.798782) z=(-0.027663)*x+(-0.005448)*y+(70.356472)11.1 z=(-0.003740)*x+(-0.002333)*y+(5.040576) z=(-0.002530)*x+(-0.002555)*y+(10.061014) z=(-0.004408)*x+(-0.002633)*y

38、+(16.522721)z=(-0.002736)*x+(-0.003593)*y+(20.518390) z=(-0.001351)*x+(-0.004259)*y+(24.823896) z=(-0.021377)*x+(-0.004712)*y+(40.045928) z=(-0.018205)*x+(-0.007014)*y+(43.305695) z=(-0.021560)*x+(-0.006476)*y+(48.273096)z=(-0.025300)*x+(-0.004743)*y+(52.729697) z=(-0.026175)*x+(0.000449)*y+(53.5629

39、83) z=(-0.024393)*x+(-0.002167)*y+(58.395736) z=(-0.018681)*x+(-0.004368)*y+(60.918705) z=(-0.028368)*x+(-0.005161)*y+(70.556460)平面擬合： x=xlsread(D:data.xls,C4:C11);y=xlsread(D:data.xls,D4:D11);z=xlsread(D:data.xls,E4:E11);plot3(x,y,z);scatter3(x,y,z,filled);hold on;X = ones(8,1) x y; b = regress(z,X

40、);hold on;xfit = min(x):0.1:max(x); yfit = min(y):0.1:max(y);XFIT,YFIT= meshgrid (xfit,yfit); ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT;mesh (XFIT,YFIT,ZFIT);title(sprintf( z=(%f)*x+(%f)*y+(%f),b(3), b(2),b(1);兩直線夾角：function a=JiaJiao(x,y)% 求兩條直線夾角% x,y 是已知三點(diǎn)的橫坐標(biāo)和縱坐標(biāo)% eg: x=1 2 3;y=4 1 5;if x(2)=x(1)

41、k1=(y(2)-y(1)/(x(2)-x(1);endif x(3)=x(2)k2=(y(3)-y(2)/(x(3)-x(2);endif x(2)=x(1) & x(3)=x(2)a=0;elseif x(3)=x(2) a=pi/2-atan(abs(k1);elseif x(1)=x(2)a=pi/2-atan(abs(k2);elseif 1+k1*k2=0a=pi/2;elsea=atan(abs(k2-k1)/(1+k2*k1); % 夾角enda=a*360/(2*pi); % 轉(zhuǎn)化為角度制向量的夾角：A=1 2 9;B=0 1 0;acos(dot(A,B)/(norm(A)

42、*norm(B)*180/pi曲率:clc; clear all; close all;x0 = linspace(0, 1);y0 = sin(x0).*cos(x0);h = abs(diff(x0(2), x0(1);% 模擬一階導(dǎo)figure; box on; hold on;ythe1 = cos(x0).2 - sin(x0).2; %理論一階導(dǎo)yapp1 = gradient(y0, h); %matlab數(shù)值近似plot(x0, ythe1, .);plot(x0, yapp1, r);legend(理論值, 模擬值);title(模擬一階導(dǎo));% 模擬二階導(dǎo)figure; box on; hold on;ythe2 = (-4)*cos(x0).*sin(x0); %理論二階導(dǎo)yapp2 = 2*2*del2(y0, h); %matlab數(shù)值近似plot(x0, ythe2,.);plot(x0, yapp2,r);legend(理論值, 模擬值);title(模擬二階導(dǎo));% 模擬曲率syms x yy = sin(x)*cos(x);yd2 = diff(y, 2);yd1 = diff(y, 1);k = abs(yd2)/(1+yd12)(3/2);k1 = subs(k, x, x0);k2 = abs(yapp2)./(