2013年国赛高教杯数学建模
C题 古塔的变形
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。
某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。
请你们根据附件1提供的4次观测数据，讨论以下问题：
1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。
2. 分析该塔倾斜、弯曲、扭曲等变形情况。
3. 分析该塔的变形趋势。
整体求解过程概述(摘要)
本文要求根据测绘公司对古塔的4次测量数据，给出确定古塔各层中心位置的通用方法，并分析古塔的变形情况及其变形趋势。为了计算的精度，我们首先对各变形量进行了合理的数学定义，并对附录的缺失数据进行合理的赋值。
对于问题一，我们通过最小二乘法拟合出观测点所在平面，再建立优化模型，在拟合平面上寻找到各观测点距离的平方和最小的点作为古塔该层的中心点。利用MATLAB编程求解，得到了每次观测古塔各层中心坐标的通用方法及各层的中心点坐标。
对于问题二，我们将古塔的倾斜、弯曲和扭曲等变形情况，分别给予合理的数学描述。对于倾斜变形，我们定义了倾斜角α，即塔尖与底层中心的水平距离与塔高的比值；对于弯曲变形，我们定义了弯曲率K，即用中心点所拟合出的空间曲线的曲率来描述古塔各处弯曲率；对于扭曲变形，我们定义了相对扭曲度θ，利用坐标的旋转变换角度描述古塔的扭曲变形情况。利用空间曲线拟合、坐标变换等方法以及MATLAB程序，分别求出了三个变形刻画量的量化指标。
对于问题三，我们考虑通过古塔的倾斜、弯曲及扭曲程度来分析古塔的变形趋势。由于数据量较少，我们建立灰色预测模型分析这三种变形因素的变化趋势，利用相应的MATLAB程序，得到了倾斜角、弯曲率以及相对扭曲度的预测函数和误差检验，验证了模型的可靠性，并继而分析古塔的变形趋势。
本文巧妙地将各种变形量给予了合理的数学描述及模型，并运用最小二乘法、曲线投影拟合、坐标变换等数学方法实现了求解，并利用灰色预测对未来变形趋势进行了预测，具有较好的实用性和可推广性。
模型假设：
1.由于中国古代建筑物多为对称图形，假设古塔是对称的。
2.假设每次古塔的测量点选取是固定的。
3.假设测量数据都是准确可靠的。
4.假设古塔的变形只由倾斜、弯曲和扭曲变形造成，不考虑其他因素。
问题分析：
问题一要求确定古塔各层中心位置的通用方法。根据建筑变形测量规范，在建筑物变形测量中，为更好地测量出建筑物变形程度的各个指标，我们假设每次测量应选取固定的测量点，且在同一层所选取的测量点在未变形前处于同一个水平面上。而经过对各层观测点三维散点图的绘制发现，各层的八个点近似对称地分布在一个平面上，只是因为年代久远发生变形导致了些许偏差。因此为了更准确地找出各层中心点，我们考虑先利用最小二乘法拟合出各层观测点所在的平面方程，再建立优化模型在该平面上寻找一点使其到各观测点距离的平方和最小，以此确立古塔各层中心坐标。
问题二要求分析古塔的各种变形情况。根据《中华人民共和国行业标准建筑变形测量规范(JGJ8—2007)》知，变形是建筑的地基、基础、上部结构及其场地受各种作用力而产生的形状或位置变化现象。在本问中，我们主要分析古塔三种主要的变形情况：倾斜、弯曲、扭曲。对于倾斜变形，我们定义倾斜角α进行描述 ，其正切值等于塔尖与底层中心的水平距离与塔高的比值，即tanα=d/h ；对于弯曲变形，我们首先通过投影法拟合出古塔各层中心点所在空间曲线的参数方程，再利用空间曲线的曲率来刻画古塔的弯曲度K；对于扭曲变形，考虑到扭曲变形实际为古塔水平面的旋转产生，因此我们采用二维坐标( x, y)旋转的矩阵变换，通过各观测量点前后的坐标确定古塔的旋转角度θ，以此刻画古塔的扭曲度。但是，实际中水平面坐标( x, y)不仅发生了旋转变换，还受到倾斜弯曲变形等所引起的平移变化的影响，因此我们在考虑坐标变换的时候加入了平移量( p,q )使其更加准确合理。
问题三为分析古塔的变形情况。本文中，我们认为建筑物变形由建筑物的倾斜、弯曲、扭曲等因素共同造成。由于附录只给出了四次统计的数据，而我们的目标是分析古塔未来多年的变化趋势，因此我们采用信息不完全、不充分的预测系统——灰色预测对古塔未来的变形趋势进行预测。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
clc,clear
x0=[0.0141,0.0142,0.0146,0.0147];%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
clc
clear
x1=[0.000141404 0.000121639 0.000089860 0.000056555];
x2=[0.000141405 0.000121641 0.000089920 0.000056556 ];
x3=[0.000141406 0.000121642 0.000089977 0.000056556 ];
x4=[0.000141407 0.000121643 0.000090030 0.000056557 ];
x5=[0.000141408 0.000121644 0.000090089 0.000056557 ];
x6=[0.000141408 0.000121645 0.000090149 0.000056558 ];
x7=[0.000141409 0.000121646 0.000090200 0.000056558 ];
x8=[0.000141409 0.000121646 0.000090250 0.000056558 ];
x9=[0.000141409 0.000121647 0.000090301 0.000056559 ];
x10=[0.000141409 0.000121647 0.000090352 0.000056559 ];
x11=[0.000141408 0.000121648 0.000090417 0.000056559 ];
x12=[0.000141408 0.000121648 0.000090486 0.000056560 ];
x13=[0.000141408 0.000121648 0.000090555 0.000056560 ];
x14=[0.000141407 0.000121648 0.000090599 0.000056560 ];x0=x1(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a1(1,:)=x0;
a1(2,:)=x0_hat;
a1(3,:)=epsilon;
a1(4,:)=delta;x0=x2(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a2(1,:)=x0;
a2(2,:)=x0_hat;
a2(3,:)=epsilon;
a2(4,:)=delta;x0=x3(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a3(1,:)=x0;
a3(2,:)=x0_hat;
a3(3,:)=epsilon;
a3(4,:)=delta;x0=x4(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a4(1,:)=x0;
a4(2,:)=x0_hat;
a4(3,:)=epsilon;
a4(4,:)=delta;x0=x5(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a5(1,:)=x0;
a5(2,:)=x0_hat;
a5(3,:)=epsilon;
a5(4,:)=delta;x0=x6(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a6(1,:)=x0;
a6(2,:)=x0_hat;
a6(3,:)=epsilon;
a6(4,:)=delta;x0=x7(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a7(1,:)=x0;
a7(2,:)=x0_hat;
a7(3,:)=epsilon;
a7(4,:)=delta;x0=x8(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a8(1,:)=x0;
a8(2,:)=x0_hat;
a8(3,:)=epsilon;
a8(4,:)=delta;x0=x9(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a9(1,:)=x0;
a9(2,:)=x0_hat;
a9(3,:)=epsilon;
a9(4,:)=delta;x0=x10(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a10(1,:)=x0;
a10(2,:)=x0_hat;
a10(3,:)=epsilon;
a10(4,:)=delta;x0=x11(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a11(1,:)=x0;
a11(2,:)=x0_hat;
a11(3,:)=epsilon;
a11(4,:)=delta;x0=x12(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a12(1,:)=x0;
a12(2,:)=x0_hat;
a12(3,:)=epsilon;
a12(4,:)=delta;x0=x13(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a13(1,:)=x0;
a13(2,:)=x0_hat;
a13(3,:)=epsilon;
a13(4,:)=delta;x0=x14(1,:);%原始数据序列
n=length(x0);
a_x0=diff(x0)';%求1次累减序列，即1阶向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘拟合参数
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二阶微分方程的符号解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知数据点1次累加序列的预测值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知数据点的预测值
epsilon=x0-x0_hat%求残差
delta=abs(epsilon./x0)%求相对误差
a14(1,:)=x0;
a14(2,:)=x0_hat;
a14(3,:)=epsilon;
a14(4,:)=delta;
