高教社杯全国大学生数学建模竞赛C题论文冷磊

1、古塔变形、问题分析:通过对数据的初步观察，知道个问题其实就是一个通过古塔的观测数据找到古塔变形的规律进而推测其未来发展的问题。我认为题目的关键是对这些数据进行拟合，然后根据函数的变化规律来进行预测。 先对古塔大致形态进行分析，我先用 Matlab将原始数据用三维图的 形式描绘出来，这样有助于我对问题进一步分析与求解：由于前两次数据有缺失，故选择2009年的数据：古塔侧视图:515 560古塔俯视图:通过对这个大致图形的观察，我掌握了古塔的一些基本信息，这 是一个十三层的八角古塔，每一层都是由正八边形（由于测量误差或 变形可能不那么正，但古人肯定是这么设计的），于是可以大胆假设 古塔的各层都是正

2、八边形。二、模型假设：1）假设古塔的各层质量分布均匀，即重心即为古塔几何中心2）假设古塔的各层都是正八边形，即便有些数据误差，不过可以忽略二、符号申明：三、模型求解:由于1986年和1996年的数据有缺失（如下图13层第5点）:这对于利用几何法求中心是有阻碍的,于是先将它们补齐。根据假设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:古塔各层都是正八边形，故这里用几何法进行推算由于是正八边形，如图，E为所求的点:由几何关系应该有：k ahk efk abk ed用坐标表示出来：即为: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纵坐标可直接利用其余各点平均值求得。利用Matlab求解这个方程，分别将1986和1996的数据带入方程,解得缺失的中心点的坐标分别为：年份XYZ19861996于是现在每一次数据都是齐全的了，下面来求各层的中心坐标点：首先给出各层中心点求法的通用方法：首先，先将各层都拟合到一个平面，由于竖直面上各层差异不大，于 是采用

5、对各层各点z坐标取均值的方法将古塔的每一层拟合到同一 个平面上，然后再求中心坐标，根据假设1:塔的各层分布均匀，故几 何中心即为重心，在网上找到一篇论文就是研究关于求多边形几何中 心的方法，这里不妨直接拿来用一下：定理2任边形一出一也的顶点4国*工fi=12 -按逆时的方向排列一则n边形的重心坐标为在这里，只需将n设置为8,利用Matlab求解公式即可得到各层的中心坐标。得到的结果，如下: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散点连线图:-I将四年的放在一起看:现在大致可以知道1996到2009年这个时间段古塔变形较大，下面

9、进行进一步量化分析：（一）、倾斜程度：对中心点进行线性拟合，求出拟合直线与竖坐标轴的夹角， 夹角的余弦值可作为倾斜程度的标识：Ay即:cos aXi|n |l |vxi22Yi2Zi对各个年份的中心点进行拟合，采用线性插值法得:放大后: ; / J;:丁产 果+:i=i=的奘j:申+ 194119961jJ!a41R41B41E4LB4J i J!ii：ii：Zill11killI:;:J2QC 2011:i;:1 n ： I ：J_ .1 _一二 LL._ . 1IH1IBQ解得的数据为:年份余弦值角度值正切值余弦值1986199620092011有表格可以看出，塔身与地面的夹角从变化到，逐

10、渐靠近地面，下面利用Matlab对其拟合分析：原始散点图：年份-角度用三次样条插值后:1.5604-IIIIN,：: I M k,；9iii1.5602155581,G5S51.55S41 5992W951990199520002005201020151,559由模拟的函数图可知，古塔的倾斜度每年都有所增加，但是在2000年到2010年这个时间段内变化特别明显，这说明在这个时间段内自然环境可能有比较剧烈的变化，联想到 2008年我国的汶川地震，可 以知道那段时间我国的地壳运动是比较活跃的，由此猜测也许与地壳活动剧烈有关，到了 2010年以后，变化速度急剧变慢，回到了 90年 代的速度，由此我大

11、胆的做出一个与本题无关的推测， 趋势变缓意味 着地壳运动变缓，也就是说在未来的一段日子里我国应该不会发生强 烈的地震灾害，所以大家可以安稳的生活。同时，我对古塔未来一段 时间的倾斜度变化做出一个预测：变化存在，但很缓慢，与地面的夹 角逐渐在变小，由图形走向来看，可能趋向停止（如果没有不可抗力 因素，如地震狂风的影响）。（二）、弯曲程度: 由于弯曲程度是指中心轴线的弯曲度， 将各个中心点连接起来，依次求出相邻点的之间的夹角，然后取平均值，作为弯曲度的度量：用MatLab解得为：年份余弦角度（弧度）1986199620092011用三次样条拟合后得到的图像:弯曲程度在1986到2007年之前都基本

12、保持线性，但是 2008到2010 突然变为指数级增长，估计为自然灾害影响。所以可以预测，2011年 以后，若无自然灾害或其他不可抗力的因素的影响，弯曲程度将与1986年左右持平，每年大约弯曲6.1 10 7rad ,其实这个只是很小的，至少要等上千年才能看到较明显的变化。分析这个问题我还采用了另外一种做法，就是将各层坐标投影到xoz平面上，然后通过观测曲线曲率的变化来分析弯曲程度的变化：具体做法是，分别对每年都做一次三次样条插值， 然后在曲线上抽样一些 点，取这些点的曲率的平均值作为该年度塔的弯曲度的特征值，然后通过对改特征值的分析来找到弯曲度变化的规律。对曲率拟合后的到的结果得到的结果与前

13、面基本相同（三）、扭曲程度：平面间的旋转角度可作为扭曲的度量标识， 先将各层中心与塔角的向 量与底部的同位置向量求余弦，然后将所有的点都求一次，取平均值 作为该年度塔的弯曲程度。原理示意图:cosXi X2YiY2:22222VXiYi X2Y2将这些数据利用MatLab求解得到四年的数据如下:年份余弦值1986199620092011现在对余弦值进行分析，进行三次样插值后得到的结果:由图中可以看出，其实在1996年以前几乎是没有多少扭曲的，但是 在1996到2009年这段时间里变化十分巨大，由此估计在这段时间发 生过一些比较特别的事，比如地震，这和之前的论断不谋而合，但到 了 2011年以后

14、，角度变化又趋于平缓，可以预测，在未来的几年内 如果没有不可抗力因素的影响，古塔的是不会发生太大的扭曲的。现在来观察一下古塔扭曲的方向：我的做法是先将各个平面都投影到地面， 然后平移他们，使他们几何 中心重合到一点，然后将他们边角上的 按顺序相连和拟合得到古塔 的扭曲形状。坐标变换公式为：(x, y) (x,y) (xo,y。)(x,y)为原坐标(。为中心点坐标(x, y)为平移后的坐标对边沿点的拟合曲线为:6.4: i I ) I ：； 之一g塔边沿拟台线n_二:一；J- : : : :： ：n _ _；?L!jS4kfr；j二 ,；*；：/：t；N：i金：0i：斗；:JT !：I ：;G

15、111111O-6-4-20246理论上如果古塔没有扭曲，那么边沿拟合线应该是直线，但现在边沿 线呈现出了螺旋的形状，并且还呈现出了逆时针扭曲的趋势， 至于为 什么会出现逆时针扭曲趋势，其实可以用地转偏向力来解释，地转偏 向力理论根据大地的磁场提出，北半球的物体在运动时会产生地转偏 向力，这个里的方向可以有左手理论提出：地转偏向力示意:即古塔（肯定在北半球）在偏移过程收到地转偏向力的作用，他的偏 转方向会偏向运动方向的左边，在整体上表现就是扭曲方向（俯视） 呈现逆时针。至此，三种情况个数据都已分析完全，对古塔1986到2011这段时间 内的变化分析可以得到如下结论，1986到1996年间无论是

16、扭曲，倾 斜还是弯曲表现得都很微弱，但这以后，三者的变动极大，估计这段 时间自然环境发生过一些大的变化或者是遭遇过自然的灾害（初步估 计为地震）。到了 2009年以后，变化速度急剧减小，特别是扭曲速度,基本减到了 0,倾斜速率也减小到即为微弱，但是古塔的趋势是与地面夹角在逐渐减小。五、附件：源代码：限于篇幅，值附上部分代码%对古塔形状进行大致描绘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-*);%几何法补充缺失坐标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);%对倾斜程度进行分析%?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);%对扭曲程度进行分析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-* );%对弯曲程度进行拟合 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