古塔的變形模型(共46頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上2013高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽承 諾 書我們仔細(xì)閱讀了中國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽的競(jìng)賽規(guī)則.我們完全明白，在競(jìng)賽開始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、我網(wǎng)上咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競(jìng)賽規(guī)則的, 如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競(jìng)賽規(guī)則，以保證競(jìng)賽的公正、公平性。如有違反競(jìng)賽規(guī)則的行為，我們將受到嚴(yán)肅處理。我們參賽選擇的題號(hào)是（從A/B/C/D中選擇一項(xiàng)填寫）： C 我們

2、的參賽報(bào)名號(hào)為（如果賽區(qū)設(shè)置報(bào)名號(hào)的話）： sm063 所屬學(xué)校（請(qǐng)?zhí)顚懲暾娜?山東現(xiàn)代職業(yè)學(xué)院 參賽隊(duì)員 (打印并簽名) ：1. 馬昱軒 2. 李倩 3. 梁帥健 指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)： 宋祖芳 日期：2013年 9月 16 日賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：sm06303專心-專注-專業(yè)2013高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽編 號(hào) 專 用 頁(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)）：古塔的變形模型摘要某古塔是我

3、國(guó)重點(diǎn)保護(hù)文物,已有上千年歷史。由于長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生諸如傾斜、彎曲、扭曲等各種變形。為保護(hù)古塔，文物部門需適時(shí)對(duì)古塔進(jìn)行觀測(cè)，了解各種變形量，以制定必要的保護(hù)措施。管理部門委托測(cè)繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對(duì)該塔進(jìn)行了4次觀測(cè)?；诟郊?提供的4次觀測(cè)數(shù)據(jù)：對(duì)于問題1，1986年和1996年的觀測(cè)數(shù)據(jù)中都缺少13層一個(gè)點(diǎn)的數(shù)據(jù)（因MATLAB程序中用的是循環(huán)語(yǔ)句，所以計(jì)算時(shí)賦予0值），其它各層均給出8個(gè)點(diǎn)的觀測(cè)值，為使所得數(shù)據(jù)更具真實(shí)性，確定古塔各層中心位置的方法更具適用性，本文將每層所給

4、8點(diǎn)構(gòu)成的圖形看做不規(guī)則八邊形，用中垂線求交點(diǎn)法求得古塔各層中心坐標(biāo)。對(duì)于問題2，結(jié)合問題1的分析，采用垂直投影法1求古塔的傾斜度，根據(jù)所得數(shù)據(jù)，分析古塔的傾斜程度（因1986年和1996年13層賦予值后所得數(shù)據(jù)偏差較大，為使所得數(shù)據(jù)更具真實(shí)性，所以本問題起1986年和1996年13層數(shù)據(jù)予以舍棄）；彎曲是建立在二維平面上的一條曲線，通過截取古塔過軸、軸的界面，求出每層古塔的傾斜度，從而分析得出古塔塔身在軸、軸的界面的彎曲程度。同理，也可分析軸、軸的界面的彎曲程度；扭曲同樣是采用垂直投影法1求古塔每層的傾斜度，建立三維立體空間，根據(jù)所得數(shù)據(jù)，分析古塔塔身的扭曲程度。對(duì)于問題3，利用問題1、2所

5、得數(shù)據(jù)，進(jìn)行合理的分析與猜想，進(jìn)而分析出古塔塔身的傾斜、彎曲、扭曲等變化趨勢(shì)。本文的模型解決了題目給出的問題，計(jì)算過程中充分尊重觀測(cè)數(shù)據(jù)，給出更符合實(shí)際的結(jié)果。本文所得結(jié)果大部分由圖表給出，結(jié)合圖像，較為直觀地表現(xiàn)出古塔變形情況。結(jié)果表明，采用數(shù)學(xué)軟件可得出可靠結(jié)論。關(guān)鍵詞：古塔變形監(jiān)測(cè) 垂直投影法 傾斜度一、問題的重述由于長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等。為保護(hù)古塔，文物部門需適時(shí)對(duì)古塔進(jìn)行觀測(cè)，了解各種變形量，以制定必要的保護(hù)措施。某古塔已有上千年歷史，是我國(guó)重點(diǎn)保護(hù)文物。管理部門委托測(cè)繪公司先后于1986年7月、

6、1996年8月、2009年3月和2011年3月對(duì)該塔進(jìn)行了4次觀測(cè)。請(qǐng)你們根據(jù)附件1提供的4次觀測(cè)數(shù)據(jù)，討論以下問題：1. 給出確定古塔各層中心位置的通用方法，并列表給出各次測(cè)量的古塔各層中心坐標(biāo)。2. 分析該塔傾斜、彎曲、扭曲等變形情況。3. 分析該塔的變形趨勢(shì)。二、問題的分析2.1問題的背景某古塔是我國(guó)重點(diǎn)保護(hù)文物。已有上千年歷史，由于長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生諸如傾斜、彎曲、扭曲等各種變形。對(duì)現(xiàn)存歷史文物古塔的保護(hù)，掌握古塔的形變，顯得尤為重要，文物部門需適時(shí)對(duì)古塔進(jìn)行觀測(cè)，了解各種變形量，以制定必要的保護(hù)措施。管理部門委托測(cè)繪公司先后于

7、1986年7月、1996年8月、2009年3月和2011年3月對(duì)該塔進(jìn)行了4次觀測(cè)。2.2問題的分析 基于附件1提供的4次觀測(cè)數(shù)據(jù)： 2.2.1問題一的分析 本題建立模型是利用中垂線（即中垂線相交于一點(diǎn)為中點(diǎn))法，求解古塔各層中心的坐標(biāo)，以確保古塔中心坐標(biāo)的準(zhǔn)確性。 2.2.2問題二的分析 根據(jù)層與層之間的中心點(diǎn)，塔體傾斜變化采用垂直投影法1(即垂線法，利用相對(duì)位移量進(jìn)行監(jiān)測(cè))；彎曲是建立在二維平面上的一條曲線，通過截取古塔過軸、軸的界面，求出每層古塔的傾斜度，從而分析得出古塔塔身在軸、軸的界面的彎曲程度。同理，也可分析軸、軸的界面的彎曲程度；扭曲同樣是采用垂直投影法求古塔每層的傾斜度，建立三

8、維立體空間，根據(jù)所得數(shù)據(jù)，分析古塔塔身的扭曲程度。 2.2.3問題三的分析 我們?cè)趯?duì)古塔進(jìn)行安全性監(jiān)測(cè)方面要滿足兩個(gè)方面：塔體的傾斜變化；塔身各段的變化1。針對(duì)上述幾個(gè)方面就采取對(duì)應(yīng)的監(jiān)測(cè)方法。3、 問題的假設(shè)1. 假設(shè)求解古塔各層中心坐標(biāo)時(shí)縱坐標(biāo)先不予以考慮；2. 假設(shè)古塔的地基土抗壓能力強(qiáng)、韌性大、不易收縮，對(duì)古塔的變形不產(chǎn)生任何影響；3. 假設(shè)對(duì)古塔觀測(cè)數(shù)據(jù)的缺失可忽略不計(jì)；4、 符號(hào)及文字說(shuō)明表示建立古塔坐標(biāo)軸的橫坐標(biāo)表示古塔各層的第(=1,2,.,8）個(gè)點(diǎn)的橫坐標(biāo)表示建立古塔坐標(biāo)軸的縱坐標(biāo)表示古塔各層的第（=1,2，.,8)個(gè)點(diǎn)的縱坐標(biāo)表示建立古塔坐標(biāo)軸豎坐標(biāo)表示古塔各層的第(=1,

9、2,.,8)個(gè)點(diǎn)的豎坐標(biāo)表示古塔的第條直線的斜率表示兩個(gè)點(diǎn)之間的中點(diǎn)表示兩個(gè)點(diǎn)之間的中點(diǎn)表示古塔中心點(diǎn)到水平面的距離表示古塔中心點(diǎn)的位移量表示古塔的傾斜度表示古塔的第(=1,2,.,13)層數(shù)（MATLAB程序中）表示古塔各層的第(=1,2,.,8）個(gè)點(diǎn)（MATLAB程序中）表示關(guān)于橫坐標(biāo)的矩陣（MATLAB程序中）表示關(guān)于縱坐標(biāo)的矩陣（MATLAB程序中）表示各層每?jī)牲c(diǎn)之間中垂線的斜率（MATLAB程序中）表示每層8條中垂線兩兩相交所得交點(diǎn)的關(guān)于橫坐標(biāo)的矩陣（MATLAB程序中）表示每層8條中垂線兩兩相交所得交點(diǎn)的關(guān)于縱坐標(biāo)的矩陣（MATLAB程序中）五、模型的建立與求解（1） 問題一的模型

10、建立與求解1. 通過中垂線的方法對(duì)古塔各層中心坐標(biāo)進(jìn)行求解，根據(jù)附件1表格數(shù)據(jù)，1986年和1996年的觀測(cè)數(shù)據(jù)中都缺少13層一個(gè)點(diǎn)的數(shù)據(jù)計(jì)算時(shí)賦予0值。然后觀察數(shù)據(jù)，古塔的每層八個(gè)點(diǎn)的高度變化不大，假設(shè)把古塔的每一層的點(diǎn)看作是在一個(gè)平面上，求解古塔各層中心坐標(biāo)時(shí)縱坐標(biāo)先不予以考慮，之后用平均值的方法計(jì)算出。 分析過程如下圖1-1：（用MATLAB軟件作圖（附件1，程序1） 1-1 根據(jù)上圖，做出任意兩邊的中垂線，交與一點(diǎn)，如下圖所示1-2（在上圖上用畫圖工具作圖）： 1-2由圖可知，點(diǎn)坐標(biāo)為，點(diǎn)坐標(biāo)。根據(jù)公式，得: ，；， 根據(jù)直線的斜率公式，求出直線的斜率。同理，求出直線的斜率。 再根據(jù)中

11、垂線定理，求出中垂線，的斜率，； 設(shè)中垂線的直線方程為： 設(shè)中垂線的直線方程為： 通過方程、建立方程組: 解得： ； 再把、式中所求的數(shù)據(jù)代入中，即可解出兩中垂線交點(diǎn)坐標(biāo)。同理，通過數(shù)學(xué)軟件MATLAB編程求解（附錄1，程序2）可得出每層八條中垂線兩兩相交的八個(gè)交點(diǎn)坐標(biāo)，得出數(shù)據(jù)并對(duì)這八個(gè)中垂線的交點(diǎn)求平均值，因1986年和1996年的觀測(cè)數(shù)據(jù)中都缺少13層一個(gè)點(diǎn)的數(shù)據(jù)計(jì)算時(shí)賦予0值得出的數(shù)據(jù)偏差太大，此時(shí)予以舍棄。塔頂數(shù)據(jù)是塔頂?shù)乃膫€(gè)點(diǎn)求平均值而得。各層中心坐標(biāo)如下表：表一： 古塔層次各層中心坐標(biāo)1566.522.2566.522.3566.522.4566.522.5566.522.656

12、6.522.7566.522.8566.522.9567.522.10567.522.11567.522.12567.522.13-塔頂567.24725522.24375 而每層古塔都有一定的高度差，因此，假設(shè)每層古塔的中心高度為每層八點(diǎn)的平均高度，即：，因此，古塔每層的中心坐標(biāo)是：表二：1986年古塔層次古塔各層中心坐標(biāo)1566.522.1.2566.522.7.320253566.522.12.755254566.522.17.078255566.522.21.72056566.522.26.7566.522.29.8566.522.33.9567.522.36.10567.522.40

13、.11567.522.44.12567.522.48.13-塔頂567.24725522.2437555.表三：1996年古塔層次古塔各層中心坐標(biāo)1566.522.1.7832566.522.7.3566.522.12.750754566.522.17.5566.522.21.766566.522.26.22957566.522.29.8566.522.33.9567.522.36.10567.522.4760940.11567.522.44.12567.522.48.13-塔頂567.25435522.2366555.11975表四：2009年古塔層次古塔各層中心坐標(biāo)1566.522.1.7

14、6452566.522.7.3093566.522.12.732254566.522.17.069755566.522.21.6566.522.26.2117566.991522.29.8567.522.33.9567.522.36.8437510567.522.40.11567.522.44.12567.522.48.6997513567.522.52.塔頂567.336522.214855.091表五：2011年古塔層次古塔各層中心坐標(biāo)1566.522.1.763252566.522.7.29053566.81276522.12.4566.522.17.0255566.522.21.656

15、6.522.26.20457566.9918522.29.8178567.522.33.9567.522.36.8222510567.522.40.11567.522.44.12567.522.48.13567.522.52.塔頂567.3375522.213555.087（2） 問題二的模型建立與求解1，根據(jù)問題一求出的每層的中心坐標(biāo)數(shù)據(jù)，通過分析古塔的傾斜、彎曲、扭曲及層與層之間的中心點(diǎn)，確定了塔體傾斜變化采用垂直投影法1(即垂線法，利用相對(duì)位移量進(jìn)行監(jiān)測(cè));彎曲、扭曲即為塔身各段的變化采用求傾斜度的方法。用MATLAB數(shù)學(xué)軟件編程（附件1，程序3）和畫圖工具畫圖，如下圖2-1： 2-1

16、通過對(duì)上圖的觀察與分析，每層的中心大致在一條傾斜的直線上，采用垂直投影的方法，對(duì)古塔的傾斜角進(jìn)行求解，如下圖通過MATLAB數(shù)學(xué)軟件編程和畫圖工具并用作圖（附件1，程序4），建立的數(shù)學(xué)模型，如下圖2-2； 2-2 通過對(duì)圖形的觀察，點(diǎn)為古塔塔頂?shù)闹行淖鴺?biāo)，點(diǎn)為古塔第一層的中心坐標(biāo)，點(diǎn)是點(diǎn)在點(diǎn)所在水平面上的投影， 根據(jù)傾斜度公式,3計(jì)算出的值： 根據(jù)以上公式，得出以下數(shù)據(jù)，如表六：傾斜度198619962009201171.70.69.69.比較數(shù)據(jù)： 分析數(shù)據(jù)，的值越小，傾斜程度越嚴(yán)重，古塔隨著年代的增加而越來(lái)越傾斜。2. 對(duì)于彎曲的情況，彎曲是建立在二維平面上的一條曲線，在解決問題的過程中，

17、為了對(duì)彎曲便以分析，用軸、軸所在的平面建立二維平面的數(shù)學(xué)模型，（由于古塔的彎曲方向具有不確定性，此方法具有一定的局限性）通過畫圖工具畫出簡(jiǎn)易二維圖形2-3； 2-3 根據(jù)傾斜度公式：3，求出每層的傾斜度，每層塔的的高度： 中心點(diǎn)的位移距離： 根據(jù)上述公式計(jì)算出1986、1996、2009、2011年的每層古塔的傾斜度，由于1986、1996年的第十三層丟失數(shù)據(jù)，因此舍去第十三層的斜率，用MATLAB數(shù)學(xué)軟件編程（附件1，程序5），計(jì)算結(jié)果，如下表：層次年代1986199620092011傾斜度197.95.159.159.297.96.164.164.398.97.161.159.498.98

18、.165.166.592.90.49.49.6108.107.107.107.7106.105.66.66.8101.99.66.65.9101.99.59.59.1075.75.97.97.1176.76.97.97.1276.75.82.82.13-49.49. 根據(jù)表格中所得的數(shù)據(jù)，可以分析出，每層古塔的傾斜度不同，塔身的彎曲程度也不同，傾斜度越小，塔身的彎曲程度越嚴(yán)重，通過這四個(gè)年代每層對(duì)應(yīng)的傾斜度的比較，得出古塔的塔身的彎曲程度每年都在加深。3. 扭曲等情況，需要按著上述方法計(jì)算出每層古塔的傾斜角，通過比較每層的傾 斜角，可得出古塔的扭曲情況。下圖2-4為古塔每年的中心坐標(biāo)圖（通過M

19、ATLAB數(shù)學(xué)軟件作圖（附件1，程序6），旋轉(zhuǎn)選取適當(dāng)?shù)囊暯牵?2-4 通過對(duì)上圖的觀察，每層都有一定的扭曲程度，采用第1小題中的方法求出每層古塔的傾斜度，如圖（用畫圖工具作簡(jiǎn)易圖）： 根據(jù)傾斜度公式：，求出每層的傾斜度，每層塔的的高度： 中心點(diǎn)的位移距離： 根據(jù)上述公式計(jì)算出1986、1996、2009、2011年的每層古塔的傾斜度，由于1986、1996年的第十三層丟失數(shù)據(jù)，因此舍去第十三層的斜率，用MATLAB數(shù)學(xué)軟件編程（附件1，程序7），計(jì)算結(jié)果，如下表：層次年代1986199620092011傾斜度181.80.123.122.281.80.127.128.382.81.123.

20、120.482.81.127.127.578.76.44.44.696.94.90.90.794.93.57.57.891.89.57.56.993.92.42.42.1061.61.71.71.1162.61.71.71.1235.35.60.60.13-26.26.觀察圖表中數(shù)據(jù)，（3） 問題三的模型建立與求解 古塔存在時(shí)間久遠(yuǎn)，經(jīng)過長(zhǎng)時(shí)間大自然的侵襲和破壞，并且還受到人們生產(chǎn)生活的影響。塔的內(nèi)部和外部都受到一定的影響。根據(jù)建筑力學(xué)得知，當(dāng)發(fā)生地震時(shí)，地基結(jié)構(gòu)發(fā)生了改變使得古塔發(fā)生不均勻沉降，導(dǎo)致古塔發(fā)生了變形。經(jīng)過四次觀測(cè)數(shù)據(jù)得出古塔每次都發(fā)生了很微小的變形，如建筑物超過自重，會(huì)使得抗外

21、界干擾能力及承載力降低從而造成破壞。就風(fēng)壓而言，相當(dāng)于是對(duì)古塔產(chǎn)生額外的與重力相同的外力作用，會(huì)增加內(nèi)部古塔柱梁板等構(gòu)件的內(nèi)力，可能造成構(gòu)件內(nèi)力超過其承載能力，從而造成古塔變形。就風(fēng)吸而言，建筑上一般不考慮其對(duì)古塔減輕重力的有利一面，而且需要古塔在風(fēng)的吸力作用下是否會(huì)造成風(fēng)吸力大于建筑物屋面重力的情況，而造成屋面被掀起。這樣的例子，主要表現(xiàn)在建筑物外圍雨棚的設(shè)計(jì)中，對(duì)于一般屋面，其外形一般產(chǎn)生的只是風(fēng)壓。因此會(huì)出現(xiàn)古塔的傾、彎曲、扭曲等各種變形。根據(jù)以上問題得出的數(shù)據(jù)，分析得出的古塔的傾斜、彎曲、扭曲等變化情況，如果不受外力影響（除重力外），塔身會(huì)隨著時(shí)間的變化會(huì)繼續(xù)傾斜、彎曲、扭曲等變化，如

22、圖2-4，根據(jù)四次測(cè)量的數(shù)據(jù)及分析，猜想古塔會(huì)繼續(xù)向某些方向傾斜、彎曲、扭曲等變化；如果考慮到外力對(duì)塔身的影響，塔身會(huì)發(fā)生不規(guī)則傾斜、彎曲、扭曲等變化。如果古塔到了一定的年代，塔身的變形程度過大，古塔重心不足以維持塔身的平衡，古塔將會(huì)出現(xiàn)倒塌一部分，或者完全倒塌。 6、 模型的評(píng)價(jià)與改進(jìn)方向6.1模型的優(yōu)缺點(diǎn)（1） 本模型綜合利用了、中垂線法、垂直投影法等數(shù)學(xué)方法，數(shù)學(xué)推導(dǎo)嚴(yán)謹(jǐn)，理論性強(qiáng)，所得數(shù)據(jù)真實(shí)可靠。（2）模型具有一般性，在一定條件下適用于其他古建筑物的變形情況監(jiān)測(cè)。（3） 由于數(shù)據(jù)較多，計(jì)算量過大，對(duì)于問題一，算得八條垂直平分線的八個(gè)交點(diǎn)后，采用極限思想的計(jì)算會(huì)使數(shù)據(jù)更準(zhǔn)確一些。（4）

23、 問題二中的彎曲，扭曲在已得的數(shù)據(jù)上用擬合方法會(huì)更優(yōu)化一些。6.2模型的改進(jìn)放向我們的模型通過一些合理的假設(shè)，使得問題的數(shù)學(xué)描述較為簡(jiǎn)單直觀，但實(shí)際上還有很多因素與要討論的問題密切相關(guān)，而且是應(yīng)該加以考慮的。對(duì)于問題一的求解方法，如果多次計(jì)算每次求出的數(shù)據(jù)，根據(jù)極限思想，也就是利用求極限方法求各層中心坐標(biāo)，最終的結(jié)果更無(wú)限接近于真實(shí)值。如問題二，在解決彎曲變化時(shí)，應(yīng)多方位考慮彎曲方向等問題，這些都是在研究最優(yōu)策略時(shí)應(yīng)該認(rèn)真考慮的問題，因?yàn)闀r(shí)間關(guān)系，我們只討論了比基礎(chǔ)、較簡(jiǎn)單的建模方法，還有待于完善和改進(jìn)。7、 參考文獻(xiàn)1梁?？?古塔變形測(cè)量方法探討.城市勘探2011年第03期2軟件編程；清華大

24、學(xué)出版社3韓煊,李寧.地鐵施工引起的建筑物扭曲變形分析.土木工程學(xué)報(bào)2010年1月8、 附錄附件1程序1：x=565.454562.058561.39563.782567.941571.255571.938569.5565.454;y=528.012525.544521.447518.108517.407519.857523.953527.356528.012;plot(x,y)程序2：（1986年）xa=load(xa.txt);ya=load(ya.txt);X=xa;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y

25、,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1

26、(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data1.xls,xx,1)xlswrite(data2.xls,yy,1)（1996年）xa=load(xb.txt);ya=load(yb.txt);X=xa;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)

27、/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(da

28、ta3.xls,xx,1)xlswrite(data4.xls,yy,1)（2009年）xa=load(xc.txt);ya=load(yc.txt);X=xa;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=z

29、eros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data5.xls,xx,1)xlswrite(data6.xls,yy,1)（2011年）xa=load(xd.txt);ya=load(yd.txt);X=x

30、a;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1

31、(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data7.xls,xx,1)xlswrite(data8.xls,yy,1)程序3：x=566.566.566.566.566.566.566.566.567.567.567.567.567.24725;y=522.522.522.522.522.522.522.522.522.522.522.522.522.243

32、75;z=1.7.3202512.7552517.0782521.720526.29.33.36.40.44.48.55.;plot3(x,y,z,r.)程序4：（1）x=565.454562.058561.39563.782567.941571.255571.938569.5565.48562.238561.663564.001567.995571.165571.801569.414565.506562.415561.931564.216568.048571.076571.666569.33565.526562.555562.144564.387568.091571.005571.55856

33、9.263565.548562.706562.373564.571568.136570.929571.443569.191565.57562.854562.6564.752568.18570.857571.333569.121565.671563.132562.883564.949568.172570.679571.094568.994565.77563.403563.158565.141568.164570.506570.862568.87565.868563.674563.433565.333568.156570.333570.63568.747565.961563.927563.6935

34、65.516568.148570.171570.408568.631566.078564.193563.958565.649568.094570.013570.236568.615566.195564.459564.224565.782568.039569.854570.063568.598566.308564.716564.481565.91569.701569.897568.582567.255567.235567.247567.252;y=528.012525.544521.447518.108517.407519.857523.953527.356527.764525.364521.4

35、2518.226517.563519.961523.908527.141527.52525.188521.394518.343517.716520.063523.864526.93527.327525.047521.373518.435517.838520.144523.829526.762527.119524.896521.351518.534517.969520.232523.791526.581526.915524.748521.329518.632518.095520.315523.755526.406526.652524.585521.356518.846518.346520.441

36、523.672526.167526.397524.427521.382519.055518.59520.564523.591525.933526.141524.268521.408519.263518.834520.686523.51525.701525.9524.12521.433519.462519.068520.801523.433525.482525.628523.95521.463519.607519.242520.885523.35525.259525.355523.78521.492519.753519.415520.969523.268525.037525.092523.616

37、521.521519.893521.05523.188524.822522.238522.242522.251522.244;z=1.7921.8181.7831.7691.7721.771.7941.8017.3267.3517.3147.3017.3067.3047.3247.33612.76112.78612.74912.73612.74112.7412.75812.77117.08417.10917.07217.05917.06417.06317.08117.09421.72621.75121.71421.70121.70521.70821.72321.73626.26726.3092

38、6.30826.26426.18926.13626.16426.24429.86929.91129.9129.86629.79129.73729.76529.84633.38333.42533.42433.3833.30533.25133.27933.3636.88736.92936.92836.88436.80936.75536.78336.86440.20140.21440.24440.22340.17140.03840.12940.15744.47244.48544.50544.48644.44244.30944.444.42848.74348.75648.77648.75748.713

39、48.5848.67148.69952.86652.87852.89752.8852.70352.79452.82255.12855.10855.12855.129;plot3(x,y,z,b.)（2）z=-1:1;x=565.454562.058561.39563.782567.941571.255571.938569.5565.454;y=528.012525.544521.447518.108517.407519.857523.953527.356528.012;z=000000000; plot3(x,y,z,r)程序5：（1986年）i=1:13;xa=566.566.566.566

40、.566.566.566.566.567.567.567.567.567.24725;za=1.7.3202512.7552517.0782521.720526.29.33.36.40.44.48.55.;x=xa;z=za;Q=zeros(1,13);for i=1:12 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,1)（1996年）i=1:13;xa=566.566.566.566.566.566.566.566.567.567.567.567.567.25435;za=1.7837.12.7507517.21.7626.2

41、29529.33.36.40.44.48.55.11975;x=xa;z=za;Q=zeros(1,13);for i=1:12 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,2)（2009年）i=1:13;xa=566.566.566.566.566.566.566.991567.567.567.567.567.567.567.336;za=1.76457.30912.7322517.0697521.26.21129.33.36.8437540.44.48.6997552.55.091;x=xa;z=za;Q=zeros(1,1

42、3);for i=1:13 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,3)（2011年）i=1:13;xa=566.566.566.81276566.566.566.566.9918567.567.567.567.567.567.567.3375;za=1.763257.290512.17.02521.26.204529.81733.36.8222540.44.48.52.55.087;x=xa;z=za;Q=zeros(1,13);for i=1:13 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=

43、Qxlswrite(bb.xls,Q,4)程序6：（1986年）x=566.566.566.566.566.566.566.566.567.567.567.567.567.24725;y=522.522.522.522.522.522.522.522.522.522.522.522.522.24375;z=1.7.3202512.7552517.0782521.720526.29.33.36.40.44.48.55.;plot3(x,y,z)（1996年）x=566.566.566.566.566.566.566.566.567.567.567.567.567.25435;y=522.522.

44、522.522.522.522.522.522.522.522.47609522.522.522.23665;z=1.7837.12.7507517.21.7626.229529.33.36.40.44.48.55.11975;plot3(x,y,z)（2009年）x=566.566.566.566.566.566.566.991567.567.567.567.567.567.567.336;y=522.522.522.522.522.522.522.522.522.522.522.522.522.522.2148;z=1.76457.30912.7322517.0697521.26.21129.33.36.8437540.44.48.6997552.55.091;plot3(x,y,z)（2011年）x=566.566.566.81276566.566.566.566.9918567.567.567.567.567.567.567.3375;y=522.522.522.522.522.