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

古塔变形古塔变形一、问题分析：通过对数据的初步观察，知道个问题其实就是一个通过古塔的观测数据找到古塔变形的规律进而推测其未来发展的问题。我认为题目的关键是对这些数据进行拟合，然后根据函数的变化规律来进行预测。先对古塔大致形态进行分析，我先用Matlab将原始数据用三维图的形式描绘出来，这样有助于我对问题进一步分析与求解：由于前两次数据有缺失，故选择2009年的数据：古塔侧视图：古塔俯视图：通过对这个大致图形的观察，我掌握了古塔的一些基本信息，这是一个十三层的八角古塔，每一层都是由正八边形（由于测量误差或变形可能不那么正，但古人肯定是这么设计的），于是可以大胆假设古塔的各层都是正八边形。二、模型假设：1）假设古塔的各层质量分布均匀，即重心即为古塔几何中心2）假设古塔的各层都是正八边形，即便有些数据误差，不过可以忽略二、符号申明：三、模型求解：由于1986年和1996年的数据有缺失（如下图13层第5点）：这对于利用几何法求中心是有阻碍的，于是先将它们补齐。根据假设2：古塔各层都是正八边形，故这里用几何法进行推算:由于是正八边形，如图，E为所求的点：由几何关系应该有：用坐标表示出来：即为：纵坐标可直接利用其余各点平均值求得。利用Matlab求解这个方程，分别将1986和1996的数据带入方程，解得缺失的中心点的坐标分别为：年份XYZ1986568.1644519.625352.834291996568.1704519.619152.83于是现在每一次数据都是齐全的了，下面来求各层的中心坐标点：首先给出各层中心点求法的通用方法：首先，先将各层都拟合到一个平面，由于竖直面上各层差异不大，于是采用对各层各点z坐标取均值的方法将古塔的每一层拟合到同一个平面上,然后再求中心坐标，根据假设1：塔的各层分布均匀，故几何中心即为重心，在网上找到一篇论文就是研究关于求多边形几何中心的方法，这里不妨直接拿来用一下：在这里，只需将n设置为8，利用Matlab求解公式即可得到各层的中心坐标。得到的结果，如下：1986年：1996年：2009年：2011年：散点连线图：将四年的放在一起看：现在大致可以知道1996到2009年这个时间段古塔变形较大，下面进行进一步量化分析：（一）、倾斜程度：对中心点进行线性拟合，求出拟合直线与竖坐标轴的夹角，夹角的余弦值可作为倾斜程度的标识：即：对各个年份的中心点进行拟合，采用线性插值法得：放大后：解得的数据为：年份余弦值角度值正切值余弦值19860.0106541.56014293.852230.01065419960.0107751.56002192.804840.01077520090.0115791.55921786.353880.01157920110.0116041.55919286.170690.011604有表格可以看出，塔身与地面的夹角从1.560142变化到1.559192，逐渐靠近地面，下面利用Matlab对其拟合分析：原始散点图：年份-角度用三次样条插值后：由模拟的函数图可知，古塔的倾斜度每年都有所增加，但是在2000年到2010年这个时间段内变化特别明显，这说明在这个时间段内自然环境可能有比较剧烈的变化，联想到2008年我国的汶川地震，可以知道那段时间我国的地壳运动是比较活跃的，由此猜测也许与地壳活动剧烈有关，到了2010年以后，变化速度急剧变慢，回到了90年代的速度，由此我大胆的做出一个与本题无关的推测，趋势变缓意味着地壳运动变缓，也就是说在未来的一段日子里我国应该不会发生强烈的地震灾害，所以大家可以安稳的生活。同时，我对古塔未来一段时间的倾斜度变化做出一个预测：变化存在，但很缓慢，与地面的夹角逐渐在变小，由图形走向来看，可能趋向停止（如果没有不可抗力因素，如地震狂风的影响）。（二）、弯曲程度：由于弯曲程度是指中心轴线的弯曲度，将各个中心点连接起来，依次求出相邻点的之间的夹角，然后取平均值，作为弯曲度的度量：用MatLab解得为：年份余弦角度（弧度）1986-9.2E-051.5708891996-0.000191.5709832009-0.000291.5710872011-0.00041.571193用三次样条拟合后得到的图像：弯曲程度在1986到2007年之前都基本保持线性，但是2008到2010突然变为指数级增长，估计为自然灾害影响。所以可以预测，2011年以后，若无自然灾害或其他不可抗力的因素的影响，弯曲程度将与1986年左右持平，每年大约弯曲,其实这个只是很小的，至少要等上千年才能看到较明显的变化。分析这个问题我还采用了另外一种做法，就是将各层坐标投影到xoz平面上，然后通过观测曲线曲率的变化来分析弯曲程度的变化：具体做法是，分别对每年都做一次三次样条插值，然后在曲线上抽样一些点，取这些点的曲率的平均值作为该年度塔的弯曲度的特征值，然后通过对改特征值的分析来找到弯曲度变化的规律。对曲率拟合后的到的结果得到的结果与前面基本相同。（三）、扭曲程度：平面间的旋转角度可作为扭曲的度量标识，先将各层中心与塔角的向量与底部的同位置向量求余弦，然后将所有的点都求一次，取平均值作为该年度塔的弯曲程度。原理示意图:余弦公式：将这些数据利用MatLab求解得到四年的数据如下：年份余弦值19860.99853319960.99853320090.99891620110.998916现在对余弦值进行分析,进行三次样插值后得到的结果：由图中可以看出，其实在1996年以前几乎是没有多少扭曲的，但是在1996到2009年这段时间里变化十分巨大，由此估计在这段时间发生过一些比较特别的事，比如地震，这和之前的论断不谋而合，但到了2011年以后，角度变化又趋于平缓，可以预测，在未来的几年内如果没有不可抗力因素的影响，古塔的是不会发生太大的扭曲的。现在来观察一下古塔扭曲的方向：我的做法是先将各个平面都投影到地面，然后平移他们，使他们几何中心重合到一点，然后将他们边角上的 按顺序相连和拟合得到古塔的扭曲形状。坐标变换公式为：为原坐标为中心点坐标为平移后的坐标对边沿点的拟合曲线为：理论上如果古塔没有扭曲，那么边沿拟合线应该是直线，但现在边沿线呈现出了螺旋的形状，并且还呈现出了逆时针扭曲的趋势，至于为什么会出现逆时针扭曲趋势，其实可以用地转偏向力来解释，地转偏向力理论根据大地的磁场提出，北半球的物体在运动时会产生地转偏向力，这个里的方向可以有左手理论提出：地转偏向力示意：即古塔（肯定在北半球）在偏移过程收到地转偏向力的作用，他的偏转方向会偏向运动方向的左边，在整体上表现就是扭曲方向（俯视）呈现逆时针。至此，三种情况个数据都已分析完全，对古塔1986到2011这段时间内的变化分析可以得到如下结论，1986到1996年间无论是扭曲，倾斜还是弯曲表现得都很微弱，但这以后，三者的变动极大，估计这段时间自然环境发生过一些大的变化或者是遭遇过自然的灾害（初步估计为地震）。到了2009年以后，变化速度急剧减小，特别是扭曲速度，基本减到了0，倾斜速率也减小到即为微弱，但是古塔的趋势是与地面夹角在逐渐减小。五、附件：源代码：限于篇幅，值附上部分代码1.Tempshape.m %对古塔形状进行大致描绘point=xlsread(C:Usersdsd-dellDesktopC2009.xlsx);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 : 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 = 13 x=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-*); end endend%cen=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);plot3(cen(1:13,1),cen(1:13,2),cen(1:13,3),red-*);2.Center.m %几何法补充缺失坐标point=xlsread(C:Usersdsd-dellDesktopC2011nihe.xlsx);syms cen xg ygcen=zeros(13,3);for i=0:1:12 x=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=0; for j=1:1:8 xg=xg+x(j); yg=yg+y(j); end cen(i+1,1)=xg/8; cen(i+1,2)=yg/8; cen(i+1,3)=z(1);endxlswrite(C:Usersdsd-dellDesktopC2011center.xlsx,cen); 3.QinxieNihe.m %对倾斜程度进行分析%point=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);p1=xlsread(C:Usersdsd-dellDesktopC1986center.xlsx);p2=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);p3=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);p4=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx);%xi=522.4:0.001:522.6; %yi=interp1(C,xi,linear);%mesh(xi,yi,zi);%plot3(point(1:13,1),point(1:13,2),point(1:13,3),r-*);grid on;hold on;%plot(xi,yi,black-*);%xi=522:0.01: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(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)/(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=polyval(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); xlswrite(C:Usersdsd-dellDesktopCjiaoduqinxie.xlsx,x); 4.NiuquNihe.m %对扭曲程度进行分析cen1=xlsread(C:Usersdsd-dellDesktopC1986center.xlsx);cen2=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);cen3=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);cen4=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx); p1=xlsread(C:Usersdsd-dellDesktopC1986nihe.xlsx);p2=xlsread(C:Usersdsd-dellDesktopC1996nihe.xlsx);p3=xlsread(C:Usersdsd-dellDesktopC2009nihe.xlsx);p4=xlsread(C:Usersdsd-dellDesktopC2011nihe.xlsx); syms c t %for i=1:1:13 a=cen1(i,1),cen1(i,2); for j=1:1:8 p1(j+(i-1)*8,1)=p1(j+(i-1)*8,1)-a(1); p1(j+(i-1)*8,2)=p1(j+(i-1)*8,2)-a(2); endend for i=1:1:13 a=cen2(i,1),cen2(i,2); for j=1:1:8 p2(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); endend for i=1:1:13 a=cen3(i,1),cen3(i,2); for j=1:1:8 p3(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(2); endend for i=1:1:13 a=cen4(i,1),cen4(i,2); for j=1:1:8 p4(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:12 t=0; for j=1:1:8 t=t+(p1(j,1)*p1(i*8+j,1)+(p1(j,2)*p1(i*8+j,2)/(p1(j,1)2+p1(j,2)2)0.5)*(p1(i*8+j,1)2+p1(i*8+j,2)2)0.5)/8; end c=c+t/12;endtheta(1,2)=c; c=0;for i=1:1:12 t=0; for j=1:1:8 t=t+(p2(j,1)*p2(i*8+j,1)+(p2(j,2)*p2(i*8+j,2)/(p2(j,1)2+p2(j,2)2)0.5)*(p2(i*8+j,1)2+p2(i*8+j,2)2)0.5)/8; end c=c+t/12;endtheta(2,2)=c; c=0;for i=1:1:12 t=0; for j=1:1:8 t=t+(p3(j,1)*p3(i*8+j,1)+(p3(j,2)*p3(i*8+j,2)/(p3(j,1)2+p3(j,2)2)0.5)*(p3(i*8+j,1)2+p3(i*8+j,2)2)0.5)/8; end c=c+t/12;endtheta(3,2)=c; c=0;for i=1:1:12 t=0; for j=1:1:8 t=t+(p4(j,1)*p4(i*8+j,1)+(p4(j,2)*p4(i*8+j,2)/(p4(j,1)2+p4(j,2)2)0.5)*(p4(i*8+j,1)2+p4(i*8+j,2)2)0.5)/8; end c=c+t/12;endtheta(4,2)=c;%xlswrite(C:Usersdsd-dellDesktopCniuqu1.xlsx,theta); hold on;grid on;p=zeros(13,2);for i=1:1:8 for j=0:1:12 p(j+1,1:2)=p1(8*j+i,1),p1(8*j+i,2); end plot(p(1:13,1),p(1:13,2),black-*);end plot(0,0,black-*); 5.WanquNihe.m %对弯曲程度进行拟合p1=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);p2=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);p3=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);p4=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx);p=xlsread(C:Usersdsd-dellDesktopCniuqu1.xlsx);xi=522:0.02:523;yi=interp1(p1(1:13,2),p1(1:13,3),xi,cubic);plot(xi,yi,black-*);%plot(point(1:4,1),point(1:4,3),-*);xi=1986:1:2011;yi=interp1(p(1:4,1),p(1:4,2),xi,cubic);plot(xi,yi,black-*);grid on;hold on;plot(p(1:4,1),p(1:4,2),r-*);6.QingxieDu.m 对倾斜度的计算