2021年数学建模大赛C题

承诺书
我们仔细阅读了?全国大学生数学建模竞赛章程?和?全国大学生数学建模竞赛参赛规那么?〔以下简称为“竞赛章程和参赛规那么〞，可从全国大学生数学建模竞赛网站下载〕。

我们完全明白，在竞赛开场后参赛队员不能以任何方式〔包括、电子邮件、网上咨询等〕与队外的任何人〔包括指导老师〕研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违背竞赛章程和参赛规那么的，假如引用别人的成果或其他公开的资料〔包括网上查到的资料〕，必须按照规定的参考文献的表述方式在正文引用途和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规那么，以保证竞赛的公正、公平性。

如有违背竞赛章程和参赛规那么的行为，我们将受到严肃处理。

我们受权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进展公开展示〔包括进展网上公示，在书籍、期刊和其他媒体进展正式或非正式发表等〕。

我们参赛选择的题号是〔从A/B/C/D中选择一项填写〕： C
我们的参赛报名号为〔假如赛区设置报名号的话〕：
所属学校〔请填写完好的全名〕：石家庄职业技术学院
参赛队员(打印并签名) ：1.魏鹏飞
2.邢磊
3.刘力恒
指导老师或指导老师组负责人(打印并签名)：陈佩宁
〔论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。

以上内容请仔细核对，提交后将不再允许做任何修改。

如填写错误，论文可能被取消评奖资格。

〕
日期：2013年9月16日赛区评阅编号〔由赛区组委会评阅前进展编号〕：
编号专用页
赛区评阅编号〔由赛区组委会评阅前进展编号〕：
全国统一编号〔由赛区组委会送交全国前编号〕：全国评阅编号〔由全国组委会评阅前进展编号〕：
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月、2021年3月和2021年3月对该塔进展了4次观测。

请你们根据附件1提供的4次观测数据，讨论以下问题：
1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。

2. 分析该塔倾斜、弯曲、扭曲等变形情况。

3. 分析该塔的变形趋势。

二、问题分析
本文我们是根据相关人员对古塔的观测的数据来分析该古塔的倾斜与弯曲的程度，并且分析出将来古塔的趋势走向。

首先，我们利用均值法求出各层的中心点并且拟合出图形，然后我们对该塔构建了三角形，利用三角函数求出该塔的倾斜角度，并且利用曲率算出弯曲的程度。

最后，我们利用所求出的数据以及图表进展分析得到该塔将来的开展趋势。

三、模型假设
1、假设古塔每层都是正八边形。

2、假设题目中提供的数据真实可靠。

3、假设地面平整。

4、假设每层的测量点在一个平面内。

四、符号说明
X<<古塔测量的数据x坐标S 塔身长度(1i108)
i
Y<<古塔测量的数据y坐标X中心点的x坐标(1i108)
i
Y中心点的y坐标(1i108)
X<<古塔测量的数据z坐标
i
Z中心点的z坐标N 古塔每层的测量点的个数
α塔的倾斜角
五、 模型的建立与求解
5.1 问题1模型的建立与求解
正八边形的重心等于中心，所以可以用均值法求每个面的中心点，公式如下： 根据每个面内点的坐标(
)，可求得平面的中心坐标：
1
N
i
i X
X N ==
∑ 1
N
i
i Y
Y N ==
∑
1
N
i
i Z
Z N ==
∑ 由于每个面都有八个测量点，所以在这里N=8。

在求第13层中心点时，由于缺失数据，用均值法得出的中心坐标有很大的偏向，所以在求13层中心点改用拟合法。

下面以求1996年古塔的中心点为例。

因为古塔的每层测量点都在一个平面内，所以13层的Z 轴坐标为7个测量点Z 轴坐标的平均值。

1352.83
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) 拟合的公式为：Z = 93.06*X -52730
将1352.83Z =代入拟合公式求得13567.1951X =
同上用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) 拟合公式：Z= *Y+59870
将
1352.83
Z=代入拟合公式求得
13522.2280
Y=
根据此方法，求出1986年的13层中心坐标。

塔的各层中心点坐标
年份层
次
X坐标Y坐标Z坐标年份层次X坐标Y坐标Z坐标
1986 1
1996
1
2 2
3 3 522.6256
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
塔
尖
塔尖
年份层
次
X坐标Y坐标Z坐标年份层次X坐标Y坐标Z坐标
2021 1
2021
1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13 塔
尖
塔尖
续表1
对于第二个问题，可以分成三个问题，倾角问题，弯曲问题，和扭曲问题。

古塔的倾斜角
图 1
图 2 图 3
图 4 图 5
图 6 图7
图 8
从上述八幅图中，不难看出，古塔已经向着X 轴的正方向，Y 轴的负方向发生了倾从专业角度来讲倾斜是指根底两端点倾斜方向的沉降差与其间隔 的比值，所以我们先取第一层的中心和塔尖进展倾斜角的测量。

首先测出两点间的空间间隔 ，两点间的空间间隔 就是古塔的塔身长度S ，计算公式：
222141141141()(y y )(z z )S x x =-+-+-
然后，求塔身与垂直方向的夹角的余弦值cos α，公式为：
141cos S
α=
进而求得倾斜角α的值：
141arccos z z
S
α-=
年份 1986 1996 2021 2021 倾斜角α
古塔的弯曲
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)
拟合公式：
2
567.3
()
0.415
54.24*
X
Z e
-
-
=然后代入曲率公式
''
3
'2
2
(1)
Z
K
Z
=
+
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)
拟合公式：
2
522.3
()
0.2665
55.34*
Y
Z e
-
-
=然后代入曲率公式
''
3
'2
2
(1)
Z
K
Z
=
+
古塔的扭曲
图9
如图9所示，古塔的右下方的等高线比拟密集，古塔左上方的等高线比拟稀疏，说明古塔在X轴的正方向，Y轴的负方向已经有倾斜，在将来忽略不可抗力，古塔会一直沿着这个方向倾斜。

六、模型的推广与改良
本模型简单易懂。

本模型解决了，古塔的各层中心点确实定，古塔倾斜角的求解。

改良建议，通过对倾斜角的取值，应该可以预测古塔的倒塌时间。

七、参考文献
1、石宁刘竞刘青桂高等数学中国水利水电出版社2021年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。