2013年数学建模大赛C题.doc

.\ 2013高教社杯全国大学生数学建模竞赛 承 诺 书 我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。 我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。 我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。 我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。 我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）：石家庄职业技术学院 参赛队员 (打印并签名) ：1.魏鹏飞 2.邢磊 3.刘力恒 指导教师或指导教师组负责人 (打印并签名)：陈佩宁 （论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。） 日期：2013年9月16日 赛区评阅编号（由赛区组委会评阅前进行编号）： 2013高教社杯全国大学生数学建模竞赛 编 号 专 用 页 赛区评阅编号（由赛区组委会评阅前进行编号）： 赛区评阅记录（可供赛区评阅时使用）： 评 阅 人 评 分 备 注 全国统一编号（由赛区组委会送交全国前编号）： 全国评阅编号（由全国组委会评阅前进行编号）： C题：古塔的变形 摘要 古塔由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。 对于第一个问题，求中心点坐标，采用的是均值法，由于前两次测量中第13层第5个点没有数据，要是采用均值法求中心坐标，会产生较大的误差，所以在求第13层中心坐标，采用的是拟合法。 对于第二个问题，分析古塔倾斜、弯曲、扭曲等变形情况。这个问题可以分三个小问题考虑。 1、 分析古塔的倾斜情况，先用Matlab软件绘制出，古塔的俯视图，观察古塔的倾斜情况，大致的倾斜方向，再用三角函数求出古塔的倾斜角度，再把四次算的倾斜角，做一下比较，观察古塔的倾斜状况。 2、 分析古塔的弯曲情况，首先观察X-Z坐标系中心点坐标，用Matlab软件把X-Z坐标系中的中心坐标拟合成一条曲线，求出这条曲线的曲率，然后按照上述方法求出Y-Z坐标系中心点坐标的曲线方程，求出这条线的曲率，分别观察古塔在X轴方向的弯曲情况，和Y轴方向的弯曲情况。 3、 分析古塔的扭曲情况，由于时间关系，没有分析古塔的扭曲。 对于第三个问题，分析古塔的变形趋势，可以根据第二问中的倾斜角，弯曲情况，进行简单的分析。 关键词：Matlab拟合，Matlab绘图，均值法，Matlab curve fitting软件，Matlab编程 一、 问题重述 古塔由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。 某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。 请你们根据附件1提供的4次观测数据，讨论以下问题： 1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。 2. 分析该塔倾斜、弯曲、扭曲等变形情况。 3. 分析该塔的变形趋势。 二、 问题分析 本文我们是根据相关人员对古塔的观测的数据来分析该古塔的倾斜与弯曲的程度，并且分析出未来古塔的趋势走向。首先，我们利用均值法求出各层的中心点并且拟合出图形，然后我们对该塔构建了三角形，利用三角函数求出该塔的倾斜角度，并且利用曲率算出弯曲的程度。最后，我们利用所求出的数据以及图表进行分析得到该塔未来的发展趋势。 三、 模型假设 1、 假设古塔每层都是正八边形。 2、 假设题目中提供的数据真实可靠。 3、 假设地面平整。 4、 假设每层的测量点在一个平面内。 四、 符号说明 S 塔身长度 古塔测量的数据x坐标 中心点的x坐标 古塔测量的数据y坐标 中心点的y坐标 古塔测量的数据z坐标 中心点的z坐标 N 古塔每层的测量点的个数 塔的倾斜角 五、 模型的建立与求解 5.1 问题1模型的建立与求解 正八边形的重心等于中心，所以可以用均值法求每个面的中心点，公式如下： 根据每个面内点的坐标(Xi , Yi ,Zi)，可求得平面的中心坐标： 由于每个面都有八个测量点，所以在这里N=8。 在求第13层中心点时，由于缺失数据，用均值法得出的中心坐标有很大的偏差，所以在求13层中心点改用拟合法。下面以求1996年古塔的中心点为例。 因为古塔的每层测量点都在一个平面内，所以13层的Z轴坐标为7个测量点Z轴坐标的平均值。 Matlab curvefitting软件对古塔X-Z坐标系12层中心点和一个塔顶坐标进行拟合，结果如下图： Figure 1 Linear model Poly1: Coefficients (with 95% confidence bounds): p1 = 93.06 (90.11, 96.01) p2 = -5.273e+04 (-5.44e+04, -5.106e+04) 拟合的公式为：= 93.06* -52730 将代入拟合公式求得 同上用Matlab curvefitting软件对古塔Y-Z坐标系12层中心点和一个塔顶坐标进行拟合，结果如下图： Figure 2 Linear model Poly1: f(x) = p1*x + p2 Coefficients (with 95% confidence bounds): p1 = -114.5 (-133, -96.11) p2 = 5.987e+04 (5.025e+04, 6.95e+04) 拟合公式：= -114.5* +59870 将代入拟合公式求得 依据此方法，求出1986年的13层中心坐标。 古塔各层的中心坐标，如下表： 塔的各层中心点坐标 年份 层次 X坐标 Y坐标 Z坐标 年份 层次 X坐标 Y坐标 Z坐标 1986 1 566.6648 522.7105 1.787375 1996 1 566.665 522.7102 1.783 2 566.7196 522.6684 7.32025 2 566.7205 522.6674 7.314625 3 566.7735 522.6273 12.75525 3 566.7751 522.6256 12.75075 4 566.8161 522.5944 17.07825 4 566.8183 522.5922 17.07513 5 566.8621 522.5591 21.7205 5 566.8649 522.5563 21.716 6 566.9084 522.5244 26.23513 6 566.9118 522.521 26.2295 7 566.9468 522.5081 29.83688 7 566.9506 522.5042 29.83225 8 566.9843 522.4924 33.35088 8 566.9884 522.4881 33.34538 9 567.0218 522.4764 36.85488 9 567.0265 522.4714 36.84825 10 567.0569 522.4624 40.17213 10 567.062 522.4572 40.16763 11 567.1045 522.423 44.44088 11 567.1102 522.4173 44.43538 12 567.1518 522.3836 48.71188 12 567.1578 522.3775 48.70738 13 567.204 522.286 52.83429 13 567.1951 522.228 52.83 塔尖 567.2473 522.2438 55.12325 塔尖 567.2544 522.2367 55.11975 表格 1 年份 层次 X坐标 Y坐标 Z坐标 年份 层次 X坐标 Y坐标 Z坐标 2009 1 566.7268 522.7015 1.7645 2011 1 566.727 522.7014 1.76325 2 566.764 522.6693 7.309 2 566.7642 522.669 7.2905 3 566.8001 522.6384 12.73225 3 566.8004 522.6387 12.72688 4 566.8293 522.6132 17.06975 4 566.8297 522.6127 17.052 5 566.8604 522.5866 21.70938 5 566.861 522.586 21.70388 6 566.9471 522.5342 26.211 6 566.9478 522.5335 26.2045 7 566.9792 522.5123 29.82463 7 566.98 522.5115 29.817 8 567.0305 522.4797 33.33988 8 567.0313 522.4788 33.33663 9 567.0816 522.4466 36.84375 9 567.0825 522.4457 36.82225 10 567.137 522.3937 40.16113 10 567.1381 522.3926 40.14413 11 567.1799 522.3547 44.43263 11 567.181 522.3535 44.42488 12 567.2225 522.316 48.69975 12 567.2238 522.3147 48.68388 13 567.2712 522.2715 52.81838 13 567.2725 522.2701 52.81313 塔尖 567.336 522.2148 55.091 塔尖 567.3375 522.2135 55.087 续表1 5.2问题2模型的建立 对于第二个问题，可以分成三个问题，倾角问题，弯曲问题，和扭曲问题。 5.2.1古塔的倾斜角 图 1 图 2 图 3 图 4 图 5 图 6 图 7 图 8 从上述八幅图中，不难看出，古塔已经向着X轴的正方向，Y轴的负方向发生了倾斜， 从专业角度来讲倾斜是指基础两端点倾斜方向的沉降差与其距离的比值，所以我们先取第一层的中心和塔尖进行倾斜角的测量。 首先测出两点间的空间距离，两点间的空间距离就是古塔的塔身长度S，计算公式： 然后，求塔身与垂直方向的夹角的余弦值，公式为： 进而求得倾斜角的值： 把四年的测量数据分别带入公式可以得到古塔四次测量时的倾斜角，见下表 年份 1986 1996 2009 2011 倾斜角 0.8018 0.8121 0.8377 0.8397 5.2.2古塔的弯曲 Matlab curvefitting软件对古塔进行X-Z坐标系中心坐标曲线拟合，Y-Z坐标系中心坐标曲线拟合，结果如下图。 Figure 3 General model Gauss1: f(x) = a1*exp(-((x-b1)/c1)^2) Coefficients (with 95% confidence bounds): a1 = 54.24 (50.25, 58.22) b1 = 567.3 (567.2, 567.3) c1 = 0.4135 (0.3546, 0.4723) 拟合公式：然后代入曲率公式 Figure 4 General model Gauss1: f(x) = a1*exp(-((x-b1)/c1)^2) Coefficients (with 95% confidence bounds): a1 = 55.34 (53.18, 57.49) b1 = 522.3 (522.3, 522.3) c1 = 0.2665 (0.2421, 0.291) 拟合公式：然后代入曲率公式 5.2.3古塔的扭曲 5.3分析古塔的变形趋势 图 9 如图9所示，古塔的右下方的等高线比较密集，古塔左上方的等高线比较稀疏，说明古塔在X轴的正方向，Y轴的负方向已经有倾斜，在未来忽略不可抗力，古塔会一直沿着这个方向倾斜。 六、 模型的推广与改进 本模型简单易懂。 本模型解决了，古塔的各层中心点的确定，古塔倾斜角的求解。 改进建议，通过对倾斜角的取值，应该可以预测古塔的倒塌时间。 七、 参考文献 1、 石宁 刘竞 刘青桂 高等数学 中国水利水电出版社 2010年7月 2、 梁国业 廖建平 数学建模 冶金工业出版社 2004年9月 八、 附录 1、求古塔倾斜角matlab程序 function angle = angleacos(x,y,z) %UNTITLED5 Summary of this function goes here % Detailed explanation goes here s=sqrt((x(14)-x(1))^2+(y(14)-y(1))^2+(z(14)-z(1))^2); angle=acos((z(14)-z(1))/s); angle=180/pi*angle; end 2、matlab画出古塔的3维图程序 hold on; plot3(tower1986x, tower1986y, tower1986z); plot3(tower1986x, tower1986y, tower1986z,g.); plot3(towercentre1986x, towercentre1986y, towercentre1986z,r*); xlabel(X); ylabel(Y); zlabel(Z); grid on; hold off