古塔變形的數(shù)學模型大學生數(shù)學建模競賽C題全國二等獎

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

2、到嚴肅處理。我們授權(quán)全國大學生數(shù)學建模競賽組委會，可將我們的論文以任何形式進行公開展示（包括進行網(wǎng)上公示，在書籍、期刊和其他媒體進行正式或非正式發(fā)表等）。我們參賽選擇的題號是（從A/B/C/D中選擇一項填寫）： C 我們的參賽報名號為（如果賽區(qū)設(shè)置報名號的話）： 13289002 所屬學校（請?zhí)顚懲暾娜?湖北職業(yè)技術(shù)學院 參賽隊員 (打印并簽名) ：1. 2. 3.指導教師或指導教師組負責人 (打印并簽名)： 數(shù)模組 日期：2013 年 9 月13 日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：高教社杯全國大學生數(shù)學建模競賽編 號 專 用 頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：

3、賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）：古塔變形的數(shù)學模型摘要 古塔被譽為中國古代杰出的高層建筑物，歷史悠久，值得并需要我們的保護。本文研究了關(guān)于古塔變形的問題，古塔的變形與塔身的中心緊密相關(guān)，具體分析了古塔傾斜、彎曲、扭曲的變形情況及趨勢。對于問題1，建立中心位置模型，采用多邊形組合形心的算法，求的結(jié)果是表6，7，8，9中的數(shù)據(jù)；在問題2-1中利用三角函數(shù)相關(guān)知識°°，2009°，2011°對各年份各層的中心點數(shù)據(jù)進行空間直線擬合，采用空間向量法計算兩直線

4、的夾角，結(jié)果是古塔1996年相對于198°，2009年相對于1996年未發(fā)生傾斜，2011年相對于2009°。在問題2-2中，研究塔身的彎曲建立古塔彎曲的數(shù)學模型，采用三次多項式擬合的算法，得到了古塔每年每層的彎曲率。在問題2-3中，研究塔身的扭曲建立古塔扭曲的數(shù)學模型，采用兩點之間的距離公式和扭曲加權(quán)法，得出 ，2011年。在問題3中，采用層次分析法建立古塔變形趨勢的數(shù)學模型，對古塔的變形程度和趨勢做出綜合量化評價，分析得到塔樓變形的趨勢為從1986年到1996年變形程度加深，所以塔樓變形的趨勢為從1996年到2009年變形程度減輕，從2009年到2011年變形程度稍有加

5、深。關(guān)鍵詞： 組合圖形的形心 擬合 中心變形曲線 層次分析法 扭曲加權(quán)一 問題重述古塔由于長時間承受自重、氣溫、風力等各種作用，偶然還要受地震、颶風的影響會產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等。為保護古塔，文物部門需適時對古塔進行觀測，了解各種變形量，以制定必要的保護措施。現(xiàn)有一古塔已上千年歷史，是我國重點保護文物。管理部門委托測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對該塔進行了4次觀測。請根據(jù)題目提供的4次觀測數(shù)據(jù)，討論以下問題：1. 給出確定古塔各層中心位置的通用方法，并列表給出各次測量的古塔各層中心坐標。2. 分析該塔傾斜、彎曲、扭曲等變形情況。3.

6、分析該塔的變形趨勢。二、問題分析塔的變形趨勢與塔身的中心有關(guān)，因此本題的主要目標是計算出塔中心的變化規(guī)律。研究題中數(shù)據(jù)發(fā)現(xiàn)x軸坐標值逐漸偏大，y軸坐標值逐漸偏小，z軸坐標值逐漸偏小，可以由此聯(lián)想到塔正向著x軸偏大的方向、y軸偏小的方向傾斜，且下沉。說明x軸與y軸決定傾斜方向；z軸和y軸數(shù)值變化越大，傾斜程度越大，反之則越小。因此在求解層面中心點的z坐標點取八個點z值的平均。問題1屬于計算多邊形組合形心的數(shù)學問題。因為1986年受自然、社會影響度較小，最接近古塔的原始圖樣，所以圖形較為有參考價值，用MATLAB畫出1986年塔身的大致圖樣為八邊形（圖1-1），圖1-1通過假設(shè)得出正多邊形的形心就

7、是中心，因此可運用組合圖形形心公式計算出古塔各層中心坐標。問題2-0塔的傾斜、彎曲、扭曲都與塔的中心有關(guān)。從數(shù)據(jù)點上看，四年13組中心點的數(shù)據(jù)大致分布在空間的一個平面上（1986年圖2-0-1，1996年圖2-0-2，2009年圖2-0-3，2011年圖2-0-4），所以用空間直線的最小二乘法來擬合中心線。 圖2-0-1 圖2-0-2 圖2-0-3 圖2-0-4問題2-1塔的傾斜程度用傾斜角W表示，塔身擬合后的中心線與地面的夾角為傾斜角（圖2-1-1）。 圖2-1-1問題2-1.1塔自身傾斜情況。觀察圖2-1-1得到RtOPQ，利用三角形正切公式確定自身傾斜角W。1986、1996年題目分別給

8、有4組塔尖坐標，分析數(shù)據(jù)結(jié)合實際生活我們將4組數(shù)據(jù)求平均值，得到1組數(shù)據(jù)用于計算塔的傾斜角W。W可表現(xiàn)塔傾斜幅度。兩直線的方向向量的夾角（小于90O）叫做直線的夾角，現(xiàn)古塔的傾斜度一定小于90O，所以各年份中心線夾角可轉(zhuǎn)化為兩直線的方向向量的夾角，將每年各層中心和塔尖的數(shù)據(jù)用空間直線的最小二乘法進行擬合成中心線。用法向量法得到各中心直線的方向向量，利用兩向量的夾角的余弦公式確定相對傾斜角。問題2-2塔的彎曲與空間直角坐標系的z軸有關(guān)。以塔第一層樓面的中心為原點建立空間直角坐標系，擬合各年各層中心點豎坐標的直線，多次擬合后發(fā)現(xiàn)三次多項式的擬合效果較好（附錄圖2-2-1，圖2-2-2，圖2-2-3

9、），問題2-3塔的扭曲在空間直角坐標系中的x,y軸有關(guān)。扭曲是相鄰年份相對的物理量，以每層中心點偏移的平面距離作為扭曲距離，并用層次分析法對各層的扭曲距離賦予權(quán)重，得到總扭曲距離來衡量不同年份古塔的扭曲程度問題3運用層次分析法對古塔的變形程度做綜合量化評價，并根據(jù)組合權(quán)向量對古塔變形趨勢做出大致分析。3、 模型假設(shè)1，參照中國古塔的塔身樣式假設(shè)塔身水平面是正八角形；2，假設(shè)古塔是均質(zhì)物體，形心就是中心； 3，以地平面為X軸，建立空間直角坐標系X，Y，Z；4，假設(shè)塔尖面積極小，近視為一個點；5，假設(shè)模型4中古塔所在地不存在地震等這種偶然事件；6，假設(shè)古塔所在的地平面是水平面。4 符號說明Ci-第

10、i個三角形的形心（i=1,2,6）；Si-第i個三角形的面積（單位：m2)（i=1,2,6）；Ai-第i個三角形的形心坐標（單位：m）；Wt-第t年塔身的傾斜角(t=1986,1996,2009,2011)；Dt-第t年塔身在地面上的水平投影(t=1986,1996,2009,2011)；Ht-第t年塔上部與下部的相對位移(t=1986,1996,2009,2011)；Lt-第t年塔的擬合直線(t=1986,1996,2009,2011)；R-塔的傾斜率；mt-第t年中心直線x軸的法向量(t=1986,1996,2009,2011)；nt-第t年中心直線y軸的法向量(t=1986,1996,2

11、009,2011)；Jt-第t年的中心直線方向向量(t=1986,1996,2009,2011)；K-塔的曲率；5 模型的建立與求解5.1、模型的準備（1） 給數(shù)據(jù)用圖像形式表現(xiàn)出來題目所給的4次觀測數(shù)據(jù)1986年、1996年塔13層的第五個點數(shù)據(jù)殘缺。參照其10層、11層、12層的數(shù)據(jù)知道第五個點的數(shù)值在第四個點和第六個點的數(shù)值之間，另外分析后發(fā)現(xiàn)這兩組數(shù)據(jù)對題目結(jié)果的影響不大。因此結(jié)合實際和理論，采取取相鄰年份該層的第四個點和第六個點的平均值將兩組數(shù)據(jù)補齊。（2） 用MATLAB軟件將題中所給數(shù)據(jù)及補充的數(shù)據(jù)畫圖表現(xiàn)出來。5.2 問題1 中心位置模型問題一求古塔各層中心坐標，通過假設(shè)知道實

12、質(zhì)就是計算正多邊形的中心，運用組合圖形求形心的方法。將每層正八邊形的塔面都分為六個三角形（圖1-2）S1、S2、S3、S4、S5、S6，圖5-1-2每年各層六個三角形的形心：C1(，)、C2(,)、C3(,)、C4(,)、C5(,)、C6(,)， 代入測量數(shù)據(jù)，運行MATLAB軟件得出：1986年各層分面積的形心塔層坐標(m)A1A2A3A4A5A61565.6707，526.9707 564.3160，524.7823 567.8320，525.6177 234567891011567.9245，519.90551213表11996年各層分面積的形心塔層坐標(m)A1A2A3A4A5A612

13、4345678910111213表22009年各層分面積的形心塔層坐標(m)A1A2A3A4A5A61234567488910111213表32011年各層分面積的形心塔層坐標(m)A1A2A3A4A5A61564.6856，5204421234567891011126513表4每層六個三角形的面積S1、S2、S3、S4、S5、S6，代入測量數(shù)據(jù)，運行MATLAB軟件得出：13層塔中各層劃分后的三角形面積塔層面積(m2)S1S2S3S4S5S612345678910117.4121213表5多邊形的形心（xc,yc,zc)公式：將表1，2，3，4，5中的數(shù)據(jù)代入形心公式中，MATLAB軟件運行

14、的結(jié)果如下： 1986年塔的中心坐標 塔層坐標(m)XcYcZc12345678910111213表61996年塔的中心坐標塔層坐標(m)XcYcZc12345678910111213表72009年塔的中心坐標塔層坐標(m)XcYcZc12345678910111213表82011年塔的中心坐標塔層坐標(m)XcYcZc12345678567.0402910111213表9找出古塔各層中心位置的通用方法：將每層平面多邊形劃分為6個相互獨立、互不相交的三角形，并求出每層的6個三角形的形心和面積，運用組合圖形的形心坐標公式得出古塔各層中心坐標（xc，yc,zc），即中心位置；表6為1986年古塔各

15、層的中心坐標，表7為1996年古塔各層的中心坐標，表8為2009年古塔各層的中心坐標，表9為2011年古塔各層的中心坐標。類別塔身第一層的中心塔尖 坐標(m) 年份xoyoZOXYZ198652199620092011表10各年塔尖平面坐標(x,y）,各年塔身第一層的中心平面坐標(xo,yo)；塔身在地面上的水平投影D:塔上部與下部的相對位移H:在直角三角形OPQ中存在：tanWt即為傾斜度；傾斜角W: 傾斜率： r=arctanWt 代入表10的數(shù)據(jù)運算得出，W1986=1.5308,W1996=1.5558,W2009=1.5564,W2011=1.5339；r1986=2.1334,r1

16、996975,r2009=0.7679,r2011°°，2009°，2011°。5.3-2 問題2-1.2 相對傾斜模型對各年份各層的中心點數(shù)據(jù)（表6、7、8、9）進行空間直線擬合，利用MATLAB軟件擬合各年份塔中心線，結(jié)果如下： L1986 x=0.0028z+567.0229 y L1996 yL2009x=0.0057z+566.9146 yL2011y由于兩平面的交線與這兩平面的法線向量都垂直，所以可取 J1986= = J1996= = J2009= = J2011= 1 = 得出中心直線的方向向量：L1986：（0.0281，-0.0473

17、,1）、L1996：(,-0.0041,1)L2009：(0.0057,-0.0041,1)、L2011：(0.181,-0.0048,1)兩向量的夾角的余弦公式：相對傾斜角：代入1986年和1996年、1996年和2009年、2009年和2011年塔中心線的方向向量的數(shù)據(jù)求得W1=0.049,W2=0,W3=0.0245；即1996年相對于1986年傾斜了°，2009年相對于1996年沒發(fā)生傾斜，2011年相當于2009年傾斜了°。5.3-2 問題2-1 彎曲模型運用表6、7、8中豎坐標的數(shù)據(jù)在MATLAB中進行三次多項式的擬合，得到：y1986x3x2y1996x3x2

18、y2009=x3x2y2011x3x2對擬合后的直線進行一階求導、二階求導，1986x219861996x219962009x220092011x22011曲率計算公式為：MATLAB軟件運行后的結(jié)果：1986年每層曲率：1996年每層曲率：2009年每層曲率：2011年每層曲率：5.3-3 問題2-1 扭曲模型對于古塔的扭曲，我們以古塔每層的中心點作為研究對象，研究其橫縱坐標的變化，通過比較相鄰兩年的中心坐標，算出對應(yīng)的每層中心點偏移的平面距離作為扭曲距離，設(shè)每層扭曲距離的權(quán)重為wi，每層中心點的扭曲距離為mi，則建立古塔總扭曲距離型為：總扭曲距離： 通過matlab軟件算得扭曲距離如下表所

19、示：1986年-1996年1996年-2009年2009年-2011年m1m2m3m4m5m6m7m8m9m10m11m12m13考慮到古塔每層扭曲程度不同，所以對古塔每層的扭曲距離賦予扭曲權(quán)重wi，參考彎曲中每層不同的彎曲情況，用層次分析法對古塔每層的扭曲距離做綜合量化分析，得到相應(yīng)的扭曲權(quán)重，層次結(jié)構(gòu)圖如下表所示：通過層次分析法算得準則層相對于方案層的權(quán)重為w（ 0.6774 0.4552 0.4552 ） 結(jié)合總扭曲距離公式和matlab軟件算得，相應(yīng)年份的總扭曲距離為1996年總扭曲距離11.9399, 2009年總扭曲距離8.4095，2011年總扭曲距離5.4 問題3 變形趨勢模型

20、1.建立遞階層次結(jié)構(gòu) 在古塔變形趨勢問題中，要利用已有年份數(shù)據(jù)描述古塔的變形趨勢，我們采用綜合量化大致分析古塔的變形趨勢，即決策目標是“古塔變形的程度”。 為了實現(xiàn)這一目標，需要考慮的主要準則有三個，即傾斜、彎曲、扭曲。，從相互關(guān)系上分析，這些因素隸屬于主要準則，因此放在下一層次考慮，并且分屬于不同準則。 假設(shè)本問題只考慮這些準則，接下來需要明確為了實現(xiàn)決策目標、在上述準則下可以有哪些方案。根據(jù)題中所述，本問題有四個年份，這四個年份作為措施層元素放在遞階層次結(jié)構(gòu)的最下層。很明顯，這四個年份都與所有準則都相關(guān)。將各個層次的因素按其上下關(guān)系擺放好位置，并將它們之間的關(guān)系用連線連接起來。同時，為了方

21、便后面的定量表示，從上到下用A、B、C、D代表不同層次，同一層次從左到右用1、2、3、4代表不同因素。這樣構(gòu)成的遞階層次結(jié)構(gòu)如下圖。 2. 構(gòu)造判斷矩陣并賦值 根據(jù)遞階層次結(jié)構(gòu)就能很容易地構(gòu)造判斷矩陣。 構(gòu)造判斷矩陣的方法是：每一個具有向下隸屬關(guān)系的元素（被稱作準則）作為判斷矩陣的第一個元素（位于左上角），隸屬于它的各個元素依次排列在其后的第一行和第一列。 重要的是填寫判斷矩陣。填寫判斷矩陣的方法有：大多采取的方法是：向填寫人（專家）反復詢問：針對判斷矩陣的準則，其中兩個元素兩兩比較哪個重要，重要多少，對重要性程度按1-9賦值（重要性標度值見表1） 重要性標度含義表重要性標度含 義1表示兩個元

22、素相比，具有同等重要性3表示兩個元素相比，前者比后者稍重要5表示兩個元素相比，前者比后者明顯重要7表示兩個元素相比，前者比后者強烈重要9表示兩個元素相比，前者比后者極端重要2，4，6，8表示上述判斷的中間值倒數(shù)若元素I與元素j的重要性之比為aij, 則元素j與元素I的重要性之比為aji=1/aij表1 設(shè)填寫后的判斷矩陣為A=(aij)n×n，判斷矩陣具有如下性質(zhì)：(1) aij0(2) aji=1/ aji(3) aii=1 根據(jù)上面性質(zhì)，判斷矩陣具有對稱性，因此在填寫時，通常先填寫aii=1部分，然后再僅需判斷及填寫上三角形或下三角形的n(n-1)/2個元素就可以了。 在特殊情況

23、下，判斷矩陣可以具有傳遞性，即滿足等式：aij*ajk=aikAB1B2B3B111/31/5B2311/2B3521當上式對判斷矩陣所有元素都成立時，則稱該判斷矩陣為一致性矩陣填寫后的判斷矩陣如下：B1C1 C2 C3 C4C111/41/41/2C24114C34114C421/41/41B2C1C2C3C4C111/21/21/5C22111/5C32111/5C45551B3C1C2C3C4C111/81/61/3C28136C361/314C431/61/41CRCIAB1B206B33. 層次單排序（計算權(quán)向量）與檢驗 對于專家填寫后的判斷矩陣，利用一定數(shù)學方法進行層次排序。 層次

24、單排序是指每一個判斷矩陣各因素針對其準則的相對權(quán)重，所以本質(zhì)上是計算權(quán)向量。計算權(quán)向量有特征根法、和法、根法、冪法等，這里簡要介紹和法。 和法的原理是對于一致性判斷矩陣，每一列歸一化后就是相應(yīng)的權(quán)重。對于非一致性判斷矩陣，每一列歸一化后近似其相應(yīng)的權(quán)重，在對這n個列向量求取算術(shù)平均值作為最后的權(quán)重。具體的公式是： 需要注意的是，在層層排序中，要對判斷矩陣進行一致性檢驗。一致性檢驗的步驟如下： 第一步，計算一致性指標C.I.（consistency index）第二步，查表確定相應(yīng)的平均隨機一致性指標R.I.（random index）據(jù)判斷矩陣不同階數(shù)查下表，得到平均隨機一致性指標R.I. 平

25、均隨機一致性指標R.I.表（1000次正互反矩陣計算結(jié)果）矩陣階數(shù)12345678R.I.00矩陣階數(shù)9101112131415R.I.表3 第三步，計算一致性比例C.R.（consistency ratio）并進行判斷 當C.R.<0.1時，認為判斷矩陣的一致性是可以接受的，C.R.>0.1時，認為判斷矩陣不符合一致性要求，需要對該判斷矩陣進行重新修正。 決策：計算權(quán)向量及檢驗 上例計算所得的權(quán)向量及檢驗結(jié)果見下：層次計算權(quán)向量 A權(quán)向量值 B1B2 B3表4 5. 結(jié)果分析從方案層總排序的結(jié)果看，1996年的權(quán)重值大于2011年的權(quán)重值大于2009年的權(quán)重值大于1986年的權(quán)重

26、值，即0.4233>0.2566>0.2497>0.0645,所以塔樓變形的趨勢為從1986年到1996年變形程度加深，從1996年到2009年變形程度減輕，從2009年到2011年變形程度稍有加深。從準則層相對于目標層的權(quán)向量值來看，扭曲對古塔變形程度的影響最大，其次是彎曲，再次是傾斜。 .六、模型的評價與推廣優(yōu)點：1，模型5-3-2采用三次多項式擬合，擬合數(shù)據(jù)精度高；2，文中有大量的圖表，直觀形象；3，模型5-2.2可以簡單快捷的得到結(jié)果；4，對古塔變形的種類進行層次分析法，合理分析出三種變形對塔身整體影響的程度；缺點：1，模型5.2-2.1中運算量大，且聯(lián)系緊密，易發(fā)生

27、計算錯誤； 2，模型5.3-3檢驗誤差； 3，模型5.3-3,5-3-2運用了層次分析法，主觀性強，與實際生活存在誤差；七、參考文獻1 嚴宗元，高等數(shù)學，上海：同濟大學出版社，2009年8月。2 襲揚，空間直線擬合的一種方法，齊齊哈爾大學學報，第25卷第2期：64-68，2009年3月3 同濟大學數(shù)學系，高等數(shù)學第六版下，北京：高等教育出版社，1978年10月 八、附件分割成三角形求解（循環(huán)）A=;B=;C=;AB=B-A;AC=C-A;M=AB;AC;S=1/2*det(M);abs(S) x1=1; x2=2; x3=3; y1=1; y2=2; y3=3; x=(x1+x2+x3)/3

28、y=(y1+y2+y3)/3x1=564.6655;x2=566.0338;x3=567.5067;x4=566.8366;x5=568.3266;x6=569.7207; y1=521.5955;y2=520.8535;y3=521.283;y4=523.4134;y5=523.8619;y6=523.1170; s1=2.2698;s2=5.0618;s3=7.8759;s4=7.7469;s5=5.469;s6=2.2135; Xc=(x1*s1+x2*s2+x3*s3+x4*s4+x5*s5+x6*s6)/(s1+s2+s3+s4+s5+s6) Yc=(y1*s1+y2*s2+y3*s

29、3+y4*s4+y5*s5+y6*s6)/(s1+s2+s3+s4+s5+s6)x=567.1157 567.1355 566.3800 567.6101 567.1858 567.2786 567.0062 566.9845 567.0212 567.0562 567.0363 567.1611 567.3120; y=524.3733 524.2583 524.2672 523.8573 523.9879 523.7846 522.3461 522.5066 522.4868 522.4707 522.3460 522.3893 522.4024; z=1.7874 7.3203 12.7

30、553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8293; F=z;1 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)'86年的古塔圖56567.252;522.244;55.129;xi=560:0.1:580;yi=510:0.1:530;z1i=interp2(x,y,z,xi,yi','nearest');surfc(xi,yi,z1i)x=

31、567.1157 ; y=524.3733 ; z=1.7847 ; for i=2:13 m(i)=x(i)-567.1157 n(i)=y(i)-524.3733 p(i)=z(i)-1.7874 endm = Columns 1 through 12 Column 13 n = Columns 1 through 12 Column 13 p = Columns 1 through 12 Column 13 86年x=1 2 3 4 5 6 7 8 9 10 11 12 13for i=1:13K(i)=abs(0.09*x(i)-0.7706)/(1+(0.045*x(i)2-0.77

32、06*x(i)+6.8809)2)(3/2)end Columns 1 through 12 Column 13 96年x=1 2 3 4 5 6 7 8 9 10 11 12 13for i=1:13K(i)=abs(0.0912*x(i)-0.771)/(1+(0.0456*x(i)2-0.771*x(i)+6.8844)2)(3/2)endK = Columns 1 through 12 Column 13 09年x=1 2 3 4 5 6 7 8 9 10 11 12 13for i=1:13K(i)=abs(0.0912*x(i)-0.77)/(1+(0.0456*x(i)2-0.7

33、7*x(i)+6.8844)2)(3/2)endK = Columns 1 through 12 Column 13 11年x=1 2 3 4 5 6 7 8 9 10 11 12 13for i=1:13K(i)=abs(0.894*x(i)-0.7598)/(1+(0.447*x(i)2-0.7598*x(i)+6.8547)2)(3/2)endK = Columns 1 through 12 Column 13 m=0.1167 0.7556 1.2966 0.4439 0.2235 1.4641 0.1620 0.0417 0.0385 0.1263 0.1321 0.1515;w=0

34、.8845 0.8592 0.6774 0.4522 0.4522 0.6464 1.1508 2.8337 3.0338 1.3334 0.6732 0.4252 0.4252;S=m1*w1+m2*w2+m3*w3+m4*w4+m5*w5+m6*w6+m7*w7+m8*w8+m9*w9+m10*w10+m11*w11+m12*w12x=567.1157 567.1355 566.3800 567.6101 567.1858 567.2786 567.0062 566.9845 567.0212 567.0562 567.0363 567.1611 567.3120; y=524.3733

35、524.2583 524.2672 523.8573 523.9879 523.7846 522.3461 522.5066 522.4868 522.4707 522.3460 522.3893 522.4024; z=1.7874 7.3203 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8293; A=x',y',ones(size(x'); a=Az'x=567.1157 1 2 2 2 5 6 7 7 8 7 4 2; y=524.

36、3733 2 4 6 3 5 1 4 2 7 1 3 1; z=1.7847 1 1 3 5 6 3 5 3 1 6 4 2 1; for i=2:13 endx=1:1:13; y=1.7874 7.3203 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8293; A=polyfit(x,y,1) z=polyval(A,x); plot(x,y,'k+',x,z,'r')z=1.7645 7.309 12.7323 17.0698 21.

37、7094 26.211 29.8246 33.4294 36.8438 40.1368 44.4326 48.6998 52.8184;x=565.8577 566.7785 566.8098 566.8374 566.8662 566.9548 566.9879 567.0402 567.0925 567.1474 567.1901 567.2326 567.241 567.3375;y=520.9166 522.6706 522.6441 522.6197 522.5964 522.5468 522.5246 522.4123 522.2602 522.4057 522.3656 522.

38、3260 522.3260 522.2135;z=1.7633 7.2905 12.7269 17.052 21.7039 26.2045 29.817 33.3366 36.8223 40.2691 44.425 48.6839 52.8131 55.087 ;x=566.7412 567.7195 566.8095 566.8368 566.8656 566.9541 566.9871 567.0394 567.0915 567.1463 567.1897 567.2313 567.2799 567.336;y=522.7005 521.7779 522.6437 522.6234 522

39、.5970 522.5472 522.5305 522.4932 522.4604 522.5895 522.3688 522.0693 522.2822 522.2148;z=1.7645 7.309 12.7323 17.0698 21.7094 26.211 29.8246 33.4294 36.8438 40.1368 44.4326 48.6998 52.8184 55.091;F=z;1 1 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)'x=567.1157 567.1355

40、566.3800 567.6101 567.1858 567.2786 567.0062 566.9845 567.0212 567.0562 567.0363 567.1611 567.3120;y=524.3733 524.2583 524.2672 523.8573 523.9879 523.7846 522.3461 522.5066 522.4868 522.4707 522.3460 522.3893 522.4024;z=1.7874 7.3203 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44

41、.4409 48.7119 52.8293;plot3(x,y,z,'r')x=567.1157 567.1355 566.3800 567.6101 567.1858 567.2786 567.0062 566.9845 567.0212 567.0562 567.0363 567.1611 567.3120 567.2473;y=524.3733 524.2583 524.2672 523.8573 523.9879 523.7846 522.3461 522.5066 522.4868 522.4707 522.3460 522.3893 522.4024 522.243

42、8;z=1.7874 7.3203 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8293 55.1873;F=z;1 1 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)'x1=0.0028*z+567.0229;y1=-0.0473*z+524.5700;z1=z;plot3(x1,y1,z1,'b',x,y,z,'o')x=565.857

43、7 566.7785 566.8098 566.8374 566.8662 566.9548 566.9879 567.0402 567.0925 567.1474 567.1901 567.2326 567.241 567.3375;y=520.9166 522.6706 522.6441 522.6197 522.5964 522.5468 522.5246 522.4123 522.2602 522.4057 522.3656 522.3260 522.3260 522.2135;z=1.7633 7.2905 12.7269 17.052 21.7039 26.2045 29.817

44、33.3366 36.8223 40.2691 44.425 48.6839 52.8131 55.087 ;F=z;1 1 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)'x=565.8577 566.7785 566.8098 566.8374 566.8662 566.9548 566.9879 567.0402 567.0925 567.1474 567.1901 567.2326 567.241 567.3375;y=520.9166 522.6706 522.6441 522.6

45、197 522.5964 522.5468 522.5246 522.4123 522.2602 522.4057 522.3656 522.3260 522.3260 522.2135;z=1.7633 7.2905 12.7269 17.052 21.7039 26.2045 29.817 33.3366 36.8223 40.2691 44.425 48.6839 52.8131 55.087 ;F=z;1 1 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)'x=566.6531 56

46、6.3558 566.7792 566.6888 567.0364 567.0400 566.9986 566.9059 566.8944 567.0710 567.1643 567.2123 567.5044 567.2544;y=524.3733 524.2583 524.2672 523.8573 523.9879 523.7846 522.3461 522.5066 522.4868 522.4707 522.3460 522.3893 522.4024 522.2367;z=1.7645 7.309 12.7323 17.0698 21.7094 26.211 29.8246 33.

47、4294 36.8438 40.1368 44.4326 48.6998 52.8184 55.1198;F=z;1 1 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)'x0=0 3 5 7 9 11 12 13 14 15 ;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ;x=0:0.1:15;y1=lagr1(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,'spline');subplo

48、t(3,1,1)plot(x0,y0,'k+',x,y1,'r')gridtitle('lagrange')subplot(3,1,2)plot(x0,y0,'k+',x,y2,'r')gridtitle('piecewise linear')subplot(3,1,3)plot(x0,y0,'k+',x,y3,'r')gridtitle('spline')hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=inte