高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽C題論文冷磊

1、古塔變形、問(wèn)題分析:通過(guò)對(duì)數(shù)據(jù)的初步觀(guān)察，知道個(gè)問(wèn)題其實(shí)就是一個(gè)通過(guò)古塔的觀(guān)測(cè)數(shù)據(jù)找到古塔變形的規(guī)律進(jìn)而推測(cè)其未來(lái)發(fā)展的問(wèn)題。我認(rèn)為題目的關(guān)鍵是對(duì)這些數(shù)據(jù)進(jìn)行擬合，然后根據(jù)函數(shù)的變化規(guī)律來(lái)進(jìn)行預(yù)測(cè)。 先對(duì)古塔大致形態(tài)進(jìn)行分析，我先用 Matlab將原始數(shù)據(jù)用三維圖的 形式描繪出來(lái)，這樣有助于我對(duì)問(wèn)題進(jìn)一步分析與求解：由于前兩次數(shù)據(jù)有缺失，故選擇2009年的數(shù)據(jù)：古塔側(cè)視圖:515 560古塔俯視圖:通過(guò)對(duì)這個(gè)大致圖形的觀(guān)察，我掌握了古塔的一些基本信息，這 是一個(gè)十三層的八角古塔，每一層都是由正八邊形（由于測(cè)量誤差或 變形可能不那么正，但古人肯定是這么設(shè)計(jì)的），于是可以大膽假設(shè) 古塔的各層都是正

2、八邊形。二、模型假設(shè)：1）假設(shè)古塔的各層質(zhì)量分布均勻，即重心即為古塔幾何中心2）假設(shè)古塔的各層都是正八邊形，即便有些數(shù)據(jù)誤差，不過(guò)可以忽略二、符號(hào)申明：三、模型求解:由于1986年和1996年的數(shù)據(jù)有缺失（如下圖13層第5點(diǎn)）:這對(duì)于利用幾何法求中心是有阻礙的,于是先將它們補(bǔ)齊。根據(jù)假設(shè)100101J02l.l IBJ525. 05232.雉1w uni u agae. 3M2525. 335752.S&4271fi523. 61632. ETfi2364. T222523. 609852. S773”也冏521. 552.瞅3el.4.872521.514852. 889+505. 31b.

3、9. H&3884b. J1G151.L 39&B豆.104%lu55C9l 701521.0552. 7naf.56 g. 7072521.14/852. 696V5白九S9?523L LBS52. ?沁7569, 9032523. 131732, 789cufc rmuc _i n mU n CfTWhrr n iS ITTfcOrrr fT U FTm tiTi+ TT2:古塔各層都是正八邊形，故這里用幾何法進(jìn)行推算由于是正八邊形，如圖，E為所求的點(diǎn):由幾何關(guān)系應(yīng)該有：k ahk efk abk ed用坐標(biāo)表示出來(lái)：即為:yaybyeyf TOC o 1-5 h z HYPERLINK

4、l bookmark14 o Current Document %XbXeXf HYPERLINK l bookmark16 o Current Document Vayhyey HYPERLINK l bookmark18 o Current Document XaXhXeX縱坐標(biāo)可直接利用其余各點(diǎn)平均值求得。利用Matlab求解這個(gè)方程，分別將1986和1996的數(shù)據(jù)帶入方程,解得缺失的中心點(diǎn)的坐標(biāo)分別為：年份XYZ19861996于是現(xiàn)在每一次數(shù)據(jù)都是齊全的了，下面來(lái)求各層的中心坐標(biāo)點(diǎn)：首先給出各層中心點(diǎn)求法的通用方法：首先，先將各層都擬合到一個(gè)平面，由于豎直面上各層差異不大，于 是采用

5、對(duì)各層各點(diǎn)z坐標(biāo)取均值的方法將古塔的每一層擬合到同一 個(gè)平面上，然后再求中心坐標(biāo)，根據(jù)假設(shè)1:塔的各層分布均勻，故幾 何中心即為重心，在網(wǎng)上找到一篇論文就是研究關(guān)于求多邊形幾何中 心的方法，這里不妨直接拿來(lái)用一下：定理2任邊形一出一也的頂點(diǎn)4國(guó)*工fi=12 -按逆時(shí)的方向排列一則n邊形的重心坐標(biāo)為在這里，只需將n設(shè)置為8,利用Matlab求解公式即可得到各層的中心坐標(biāo)。得到的結(jié)果，如下:1986 年:566.6648522.71051.787375566.7196522.66847.32025566.7735522.627312.75525566.8161522.594417.0782556

6、6.8621522.559121.7205566.9084522.524426.23513566.9468522.508129.83688566.9843522.492433.35088567.0218522.476436.85488567.0569522.462440.172131996 年:566.665522.71021.783566.7205522.66747.314625566.7751522.625612.75075566.8183522.592217.07513566.8649522.556321.716566.9118522.52126.2295566.9506522.5042

7、29.83225566.9884522.488133.34538567.0265522.471436.84825567.062522.457240.167632009 年:566.7268522.70151.7645566.764522.66937.309566.8001522.638412.73225566.8293522.613217.06975566.8604522.586621.70938566.9471522.534226.211566.9792522.512329.82463567.0305522.479733.33988567.0816522.446636.84375567.13

8、7522.393740.161132011 年:566.727522.70141.76325566.7642522.6697.2905566.8004522.638712.72688566.8297522.612717.052566.861522.58621.70388566.9478522.533526.2045566.98522.511529.817567.0313522.478833.33663567.0825522.445736.82225567.1381522.392640.14413散點(diǎn)連線(xiàn)圖:-I將四年的放在一起看:現(xiàn)在大致可以知道1996到2009年這個(gè)時(shí)間段古塔變形較大，下面

9、進(jìn)行進(jìn)一步量化分析：（一）、傾斜程度：對(duì)中心點(diǎn)進(jìn)行線(xiàn)性擬合，求出擬合直線(xiàn)與豎坐標(biāo)軸的夾角， 夾角的余弦值可作為傾斜程度的標(biāo)識(shí)：Ay即:cos aXi|n |l |vxi22Yi2Zi對(duì)各個(gè)年份的中心點(diǎn)進(jìn)行擬合，采用線(xiàn)性插值法得:放大后: ; / J;:丁產(chǎn) 果+:i=i=的奘j:申+ 194119961jJ!a41R41B41E4LB4J i J!ii：ii：Zill11killI:;:J2QC 2011:i;:1 n ： I ：J_ .1 _一二 LL._ . 1IH1IBQ解得的數(shù)據(jù)為:年份余弦值角度值正切值余弦值1986199620092011有表格可以看出，塔身與地面的夾角從變化到，逐

10、漸靠近地面，下面利用Matlab對(duì)其擬合分析：原始散點(diǎn)圖：年份-角度用三次樣條插值后:1.5604-IIIIN,：: I M k,；9iii1.5602155581,G5S51.55S41 5992W951990199520002005201020151,559由模擬的函數(shù)圖可知，古塔的傾斜度每年都有所增加，但是在2000年到2010年這個(gè)時(shí)間段內(nèi)變化特別明顯，這說(shuō)明在這個(gè)時(shí)間段內(nèi)自然環(huán)境可能有比較劇烈的變化，聯(lián)想到 2008年我國(guó)的汶川地震，可 以知道那段時(shí)間我國(guó)的地殼運(yùn)動(dòng)是比較活躍的，由此猜測(cè)也許與地殼活動(dòng)劇烈有關(guān)，到了 2010年以后，變化速度急劇變慢，回到了 90年 代的速度，由此我大

11、膽的做出一個(gè)與本題無(wú)關(guān)的推測(cè)， 趨勢(shì)變緩意味 著地殼運(yùn)動(dòng)變緩，也就是說(shuō)在未來(lái)的一段日子里我國(guó)應(yīng)該不會(huì)發(fā)生強(qiáng) 烈的地震災(zāi)害，所以大家可以安穩(wěn)的生活。同時(shí)，我對(duì)古塔未來(lái)一段 時(shí)間的傾斜度變化做出一個(gè)預(yù)測(cè)：變化存在，但很緩慢，與地面的夾 角逐漸在變小，由圖形走向來(lái)看，可能趨向停止（如果沒(méi)有不可抗力 因素，如地震狂風(fēng)的影響）。（二）、彎曲程度: 由于彎曲程度是指中心軸線(xiàn)的彎曲度， 將各個(gè)中心點(diǎn)連接起來(lái)，依次求出相鄰點(diǎn)的之間的夾角，然后取平均值，作為彎曲度的度量：用MatLab解得為：年份余弦角度（弧度）1986199620092011用三次樣條擬合后得到的圖像:彎曲程度在1986到2007年之前都基本

12、保持線(xiàn)性，但是 2008到2010 突然變?yōu)橹笖?shù)級(jí)增長(zhǎng)，估計(jì)為自然災(zāi)害影響。所以可以預(yù)測(cè)，2011年 以后，若無(wú)自然災(zāi)害或其他不可抗力的因素的影響，彎曲程度將與1986年左右持平，每年大約彎曲6.1 10 7rad ,其實(shí)這個(gè)只是很小的，至少要等上千年才能看到較明顯的變化。分析這個(gè)問(wèn)題我還采用了另外一種做法，就是將各層坐標(biāo)投影到xoz平面上，然后通過(guò)觀(guān)測(cè)曲線(xiàn)曲率的變化來(lái)分析彎曲程度的變化：具體做法是，分別對(duì)每年都做一次三次樣條插值， 然后在曲線(xiàn)上抽樣一些 點(diǎn)，取這些點(diǎn)的曲率的平均值作為該年度塔的彎曲度的特征值，然后通過(guò)對(duì)改特征值的分析來(lái)找到彎曲度變化的規(guī)律。對(duì)曲率擬合后的到的結(jié)果得到的結(jié)果與前

13、面基本相同（三）、扭曲程度：平面間的旋轉(zhuǎn)角度可作為扭曲的度量標(biāo)識(shí)， 先將各層中心與塔角的向 量與底部的同位置向量求余弦，然后將所有的點(diǎn)都求一次，取平均值 作為該年度塔的彎曲程度。原理示意圖:cosXi X2YiY2:22222VXiYi X2Y2將這些數(shù)據(jù)利用MatLab求解得到四年的數(shù)據(jù)如下:年份余弦值1986199620092011現(xiàn)在對(duì)余弦值進(jìn)行分析，進(jìn)行三次樣插值后得到的結(jié)果:由圖中可以看出，其實(shí)在1996年以前幾乎是沒(méi)有多少扭曲的，但是 在1996到2009年這段時(shí)間里變化十分巨大，由此估計(jì)在這段時(shí)間發(fā) 生過(guò)一些比較特別的事，比如地震，這和之前的論斷不謀而合，但到 了 2011年以后

14、，角度變化又趨于平緩，可以預(yù)測(cè)，在未來(lái)的幾年內(nèi) 如果沒(méi)有不可抗力因素的影響，古塔的是不會(huì)發(fā)生太大的扭曲的?，F(xiàn)在來(lái)觀(guān)察一下古塔扭曲的方向：我的做法是先將各個(gè)平面都投影到地面， 然后平移他們，使他們幾何 中心重合到一點(diǎn)，然后將他們邊角上的 按順序相連和擬合得到古塔 的扭曲形狀。坐標(biāo)變換公式為：(x, y) (x,y) (xo,y。)(x,y)為原坐標(biāo)(。為中心點(diǎn)坐標(biāo)(x, y)為平移后的坐標(biāo)對(duì)邊沿點(diǎn)的擬合曲線(xiàn)為:6.4: i I ) I ：； 之一g塔邊沿?cái)M臺(tái)線(xiàn)n_二:一；J- : : : :： ：n _ _；?L!jS4kfr；j二 ,；*；：/：t；N：i金：0i：斗；:JT !：I ：;G

15、111111O-6-4-20246理論上如果古塔沒(méi)有扭曲，那么邊沿?cái)M合線(xiàn)應(yīng)該是直線(xiàn)，但現(xiàn)在邊沿 線(xiàn)呈現(xiàn)出了螺旋的形狀，并且還呈現(xiàn)出了逆時(shí)針扭曲的趨勢(shì)， 至于為 什么會(huì)出現(xiàn)逆時(shí)針扭曲趨勢(shì)，其實(shí)可以用地轉(zhuǎn)偏向力來(lái)解釋?zhuān)剞D(zhuǎn)偏 向力理論根據(jù)大地的磁場(chǎng)提出，北半球的物體在運(yùn)動(dòng)時(shí)會(huì)產(chǎn)生地轉(zhuǎn)偏 向力，這個(gè)里的方向可以有左手理論提出：地轉(zhuǎn)偏向力示意:即古塔（肯定在北半球）在偏移過(guò)程收到地轉(zhuǎn)偏向力的作用，他的偏 轉(zhuǎn)方向會(huì)偏向運(yùn)動(dòng)方向的左邊，在整體上表現(xiàn)就是扭曲方向（俯視） 呈現(xiàn)逆時(shí)針。至此，三種情況個(gè)數(shù)據(jù)都已分析完全，對(duì)古塔1986到2011這段時(shí)間 內(nèi)的變化分析可以得到如下結(jié)論，1986到1996年間無(wú)論是

16、扭曲，傾 斜還是彎曲表現(xiàn)得都很微弱，但這以后，三者的變動(dòng)極大，估計(jì)這段 時(shí)間自然環(huán)境發(fā)生過(guò)一些大的變化或者是遭遇過(guò)自然的災(zāi)害（初步估 計(jì)為地震）。到了 2009年以后，變化速度急劇減小，特別是扭曲速度,基本減到了 0,傾斜速率也減小到即為微弱，但是古塔的趨勢(shì)是與地面夾角在逐漸減小。五、附件：源代碼：限于篇幅，值附上部分代碼%對(duì)古塔形狀進(jìn)行大致描繪point=xlsread( C:Usersdsd-dellDesktopC);grid on;hold on ; for i=0 : 1 : 13 if i 13 x=point(i*8+1 : i*8+8 ,1); y=point(i*8+1 :

17、i*8+8,2); z=point(i*8+1 : i*8+8,3); plot3(x,y,z,red-*);%meshz(x,y,z);plot3(x(1);x(8),y(1);y(8),z(1);z(8),red-*);else if i = 13x=point(i*8+1 : i*8+4 , 1);y=point(i*8+1 : i*8+4,2);z=point(i*8+1 : i*8+4,3);plot3(x,y,z,red-*);plot3(x(1);x(4),y(1);y(4),z(1);z(4),red-*);endendend%?u ?t?D a ?o ?DD? TOC o 1

18、-5 h z cen=xlsread( C:Usersdsd-dellDesktopC);plot3(cen(1:13,1),cen(1:13,2),cen(1:13,3),red-*);%幾何法補(bǔ)充缺失坐標(biāo)point=xlsread( C:Usersdsd-dellDesktopC);syms cen xg yg cen=zeros(13,3); for i=0:1:12x=point(i*8+1 : i*8+8 , 1);y=point(i*8+1 : i*8+8,2);z=point(i*8+1 : i*8+8,3);xg=0；yg=o；for j=1:1:8xg=xg+x(j);yg=

19、yg+y(j);endcen(i+1,1)=xg/8;cen(i+1,2)=yg/8;cen(i+1,3)=z(1);endxlswrite( C:Usersdsd-dellDesktopC,cen);%對(duì)傾斜程度進(jìn)行分析%?ao?3?D ? e TOC o 1-5 h z point=xlsread( C:Usersdsd-dellDesktopC);p1=xlsread(C:Usersdsd-dellDesktopC);p2=xlsread(C:Usersdsd-dellDesktopC);p3=xlsread(C:Usersdsd-dellDesktopC);p4=xlsread(C:U

20、sersdsd-dellDesktopC);%xi=:;%yi=interp1(C,xi,linear);%mesh(xi,yi,zi);%plot3(point(1:13,1),point(1:13,2),point(1:13,3),r-*);gridon ;holdon ;%plot(xi,yi,black-*);%xi=522:523;%yi=interp1(point(1:13,2),point(1:13,3),xi,linear);%plot(xi,yi,red-*);x=zeros(4,4);p=polyfit(p1(1:13,1),p1(1:13,3),13);f=polyval

21、(p,p1(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(1,1)=1986;x(1,2)=(f(13)-f(1)/(p1(13,1)-p1(1,1);x(1,3)=atan(x(1,2);x(1,4)=cos(x(1,3);p=polyfit(p2(1:13,1),p2(1:13,3),13); f=polyval(p,p2(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(2,1)=1996;x(2,2)=(f(13)-f(1)

22、/(p2(13,1)-p2(1,1);x(2,3)=atan(x(2,2);x(2,4)=cos(x(2,3);p=polyfit(p3(1:13,1),p3(1:13,3),13);f=polyval(p,p3(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(3,1)=2009;x(3,2)=(f(13)-f(1)/(p3(13,1)-p3(1,1);x(3,3)=atan(x(3,2);x(3,4)=cos(x(3,3);p=polyfit(p4(1:13,1),p4(1:13,3),13);f=polyva

23、l(p,p4(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(4,1)=2011;x(4,2)=(f(13)-f(1)/(p4(13,1)-p4(1,1);x(4,3)=atan(x(4,2);x(4,4)=cos(x(4,3); TOC o 1-5 h z xlswrite( C:Usersdsd-dellDesktopC,x);%對(duì)扭曲程度進(jìn)行分析cen1=xlsread(C:Usersdsd-dellDesktopC)cen2=xlsread(C:Usersdsd-dellDesktopC)cen3=xls

24、read(C:Usersdsd-dellDesktopC)cen4=xlsread(C:Usersdsd-dellDesktopC)p1=xlsread(C:Usersdsd-dellDesktopC);p2=xlsread(C:Usersdsd-dellDesktopC);p3=xlsread(C:Usersdsd-dellDesktopC);p4=xlsread(C:Usersdsd-dellDesktopC);syms c t%?DD? ?6?a?6u?土？&6?for i=1:1:13a=cen1(i,1),cen1(i,2);for j=1:1:8p1(j+(i-1)*8,1)=p1

25、(j+(i-1)*8,1)-a；p1(j+(i-1)*8,2)=p1(j+(i-1)*8,2)-a(2);endendfor i=1:1:13a=cen2(i,1),cen2(i,2);for j=1:1:8p2(j+(i-1)*8,1)=p2(j+(i-1)*8,1)-a(1);p2(j+(i-1)*8,2)=p2(j+(i-1)*8,2)-a(2);endendfor i=1:1:13a=cen3(i,1),cen3(i,2);for j=1:1:8p3(j+(i-1)*8,1)=p3(j+(i-1)*8,1)-a(1);p3(j+(i-1)*8,2)=p3(j+(i-1)*8,2)-a(

26、2);endendfor i=1:1:13a=cen4(i,1),cen4(i,2);for j=1:1:8p4(j+(i-1)*8,1)=p4(j+(i-1)*8,1)-a(1);p4(j+(i-1)*8,2)=p4(j+(i-1)*8,2)-a(2);endend% theta=zeros(4,2);theta(1:1:4,1)=1986,1996,2009,2011;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p1(j,1)*p1(i*8+j,1)+(p1(j,2)*p1(i*8+j,2)/(p1(j,1)A2+p1(j,2)A2)A*(p1(i*8+j,1)

27、A2+p1(i*8+j,2)A2)A)/8;endc=c+t/12;endtheta(1,2)=c;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p2(j,1)*p2(i*8+j,1)+(p2(j,2)*p2(i*8+j,2)/(p2(j,1)A2+p2(j,2)A2)A*(p2(i*8+j,1)A2+p2(i*8+j,22)A)/8;endc=c+t/12;endtheta(2,2)=c;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p3(j,1)*p3(i*8+j,1)+(p3(j,2)*p3(i*8+j,2)/(p3(j,1)A2+p3

28、(j,2)A2)A*(p3(i*8+j,1)A2+p3(i*8+j,2)A2)A)/8;endc=c+t/12;endtheta(3,2)=c;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p4(j,1)*p4(i*8+j,1)+(p4(j,2)*p4(i*8+j,2)/(p4(j,1)A2+p4(j,2)A2)A*(p4(i*8+j,1)A2+p4(i*8+j,2)A2)A)/8;endc=c+t/12;endtheta(4,2)=c;%， %xlswrite( C:Usersdsd-dellDesktopC,theta);hold on;grid on;p=zer

29、os(13,2);for i=1:1:8for j=0:1:12p(j+1,1:2)=p1(8*j+i,1),p1(8*j+i,2); endplot(p(1:13,1),p(1:13,2),black-* );end plot(0,0, black-* );%對(duì)彎曲程度進(jìn)行擬合 TOC o 1-5 h z p1=xlsread(C:Usersdsd-dellDesktopC);p2=xlsread(C:Usersdsd-dellDesktopC);p3=xlsread(C:Usersdsd-dellDesktopC);p4=xlsread(C:Usersdsd-dellDesktopC);p=xlsread( C:Usersdsd-dellDesktopC);xi=522:523; yi=interp1(p1(1:13,2),p1