本站原创本文是一篇数学建模竞赛论文范文,关于数学建模竞赛相关函授毕业论文,关于确定古塔中心坐标的通用方法MATLAB程序相关毕业论文开题报告范文。适合数学建模竞赛及数学建模及坐标方面的的大学硕士和本科毕业论文以及数学建模竞赛相关开题报告范文和职称论文写作参考文献资料下载。
【摘 要】2013年全国大学生数学建模竞赛C题是借用文物部门的4次观测数据研究某古塔的倾斜、弯曲、扭曲等变形情况.但是要研究其变形情况,必须先确定古塔各层的中心坐标.本文首先对1986年和1996年缺失数据进行补充,其次,利用完整的数据拟合每层各测量点所在平面,最后,将各测量点投影到平面上,得到每层各中心点坐标的通用模型.并且每一步都附有相应的MATLAB计算成程序,这样做极大的减少了计算量,加快了运算速度.
【关 键 词】中心坐标,线性拟合,平面拟合
一、问题提出
2013年全国大学生数学建模竞赛已经结束了,但对竞赛题目的研究还在继续,其中C题是根据附件1提供的4次观测数据研究某古塔的倾斜、弯曲、扭曲等变形情况以及该塔的变形趋势.但是要研究其变形情况,必须先确定古塔各层的中心坐标,那么中心坐标该如何确定,大量的数据计算又该如何处理呢?
二、问题解决
1.缺失数据补充
附件1（参见2013年全国大学生数学建模竞赛C题）给出的数据中,1986年和1996年的第13层第5点的坐标是缺失的,要完整的讨论各层的中心坐标,我们需要补充缺失数据.
首先,提取1986年古塔前12层第5点的坐标,分别作x,y,z与层数t的散点图（如图1所示）.
程序如下：
A等于[567.941517.4071.772,
567.995517.5637.306,
568.048517.71612.741,
568.091517.83817.064,
568.136517.96921.705,
568.18518.09526.189,
568.172518.34629.791,
568.164518.5933.305,
568.156518.83436.809,
568.148519.06840.171],
t等于（1：12）’,x等于A（：,1）,y等于A（：,2）,z等于A（：,3）,
subplot（2,2,1）,
plot（t,x,’*’）,gridon
subplot（2,2,2）,
plot（t,y,’o’）,gridon
subplot（2,2,3）,
plot（t,z,’+’）,gridon
图11986年每层第5点的坐标分量散点图
其次,再观察图1中x坐标发现：x与层数t呈现分段线性关系,我们所求的第13层坐标处于分段函数第三段上,所以我们只需对第三段用MATLAB进行数据拟合[1],程序如下：
t1等于[ones（3,1）,（10：12）'],x等于A（10：12,1）,
[b1,bint1,r1,rint1,stats1]等于regress（y,t1）
由该程序可得x坐标与层数t的第三段线性关系为：
x等于568.6932-0.0545t（t≥10）（1）
第三,观察图1中y坐标与z坐标发现：y,z与层数t是线性关系,所以我们用MATLAB对坐标数据进行线性拟合,程序如下：
t2等于[ones（12,1）,（1：12）'],y等于A（：,2）,
[b2,bint2,r2,rint2,stats2]等于regress（y,t2）
t3等于[ones（11,1）,（2：12）’],z等于A（2：12,3）,
[b3,bint3,r3,rint3,stats3]等于regress（z,t3）
由该程序可得y和z坐标与层数t的线性关系为：
y等于517.1185+0.1880t（2）
z等于0.8997+4.0044t（由于第1点是异常点,所以剔除掉）（3）
最后,我们在（1）-（3）式中分别令t等于13可得1986年第13层第5点的坐标为：（567.9847,519.5625,52.9569）.同理得1996年第13层第5点的坐标为：（567.9903,519.5547,52.9540）.
2.平面拟合
数据补充完整后,因各层基本处于同一平面内,可先拟合各层所在平面.利用MATLAB统计工具箱中的regress命令拟合各层所在平面[2]：
x等于S（：,1）,y等于S（：,2）,z等于S（：,3）,%S为1986年的观测数据构成的矩阵
X等于[ones（8,1）,x,y],
b,bint,r,rint,stats]等于regress（z,X）
由以上程序可得拟合平面方程为：
z等于b（1）+b（2）x+b（3）y（4）
3.求中心坐标
得到各层拟合平面方程后,将各测量点投影到拟合平面内,然后再用均匀物体的重心公式计算中心坐标[3].
为了计算方便,将（4）式转化为Ax+BY+Cz+D等于0的形式,即有：
b（2）x+b（3）y-z+b（1）等于0（5）过点（xi,yi,zi）作平面的垂线,垂足即为点（xi,yi,zi）在拟合平面上的投影,垂线的直线参数方程为：
将其代入（5）式可得：
于是点（xi,yi,zi）的投影点（xi,yi,zi）为：
（）
投影点的平均值即为古塔各层的中心坐标：
其中：
4.MATLAB程序
将以上过程归纳总结,我们得到一个MATLAB程序,应用该程序可一次性得到同一年份各层的中心坐标.
我们以1986年的观测数据为例,介绍方法如下：
首先,在MATLAB中建立一个M函数文件并保存,程序为：
function[b,X,Y,Z]等于zxzb（S）%S为1986年的观测数据构成的矩阵
x等于S（：,1）,y等于S（：,2）,z等于S（：,3）,
PX等于[ones（8,1）,x,y],
[b,bint,r,rint,stats]等于regress（z,PX）,
SS等于mean（S）,
px等于SS（1）,
py等于SS（2）,
pz等于SS（3）,
u等于b（2）^2+b（3）^2+（-1）^2,
v等于b（1）+b（2）*px+b（3）*py-pz,
X等于px-b（2）/u*v,
Y等于py-b（3）/u*v,
Z等于pz+1/u*v,
其次,将1986年的数据构成的矩阵存入MATLAB中的workspace空间,并重命名为S.
最后,在MATLAB的编辑窗口中输入以下程序：
[b1,X1,Y1,Z1]等于b_zx（A（1：8,：））,
[b2,X2,Y2,Z2]等于b_zx（A（9：16,：））,
[b3,X3,Y3,Z3]等于b_zx（A（17：24,：））,
[b4,X4,Y4,Z4]等于b_zx（A（25：32,：））,
[b5,X5,Y5,Z5]等于b_zx（A（33：40,：））,
[b6,X6,Y6,Z6]等于b_zx（A（41：48,：））,
[b7,X7,Y7,Z7]等于b_zx（A（49：56,：））,
[b8,X8,Y8,Z8]等于b_zx（A（57：64,：））,
[b9,X9,Y9,Z9]等于b_zx（A（65：72,：））,
[b10,X10,Y10,Z10]等于b_zx（A（73：80,：））,
[b11,X11,Y11,Z11]等于b_zx（A（81：88,：））,
[b12,X12,Y12,Z12]等于b_zx（A（89：96,：））,
[b13,X13,Y13,Z13]等于b_zx（A（97：104,：））,
zxzb等于[X1,Y1,Z1,
X2,Y2,Z2,
X3,Y3,Z3,
X4,Y4,Z4,
X5,Y5,Z5,
X6,Y6,Z6,
X7,Y7,Z7,
X8,Y8,Z8,
X9,Y9,Z9,
X10,Y10,Z10,
X11,Y11,Z11,
X12,Y12,Z12,
X13,Y13,Z13]
运行以上程序即可得1986年各层的中心坐标,而要得到1996年、2009年和2011年各层的中心坐标只需改变存入workspace空间的数据即可.
三、结束语
使用以上方法我们得到了古塔各层的中心坐标,为检测古塔倾斜、弯曲、扭曲等各种变形做好了准备,而且MATLAB程序的使用,极大的减少了计算量,加快了运算速度.该方法也可以推广到其他建筑物中心坐标的计算中.