2013年國賽高教杯數(shù)學(xué)建模

C題 古塔的變形

??由于長時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等。為保護(hù)古塔，文物部門需適時(shí)對(duì)古塔進(jìn)行觀測(cè)，了解各種變形量，以制定必要的保護(hù)措施。
??某古塔已有上千年歷史，是我國重點(diǎn)保護(hù)文物。管理部門委托測(cè)繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對(duì)該塔進(jìn)行了4次觀測(cè)。
??請(qǐng)你們根據(jù)附件1提供的4次觀測(cè)數(shù)據(jù)，討論以下問題：
??1. 給出確定古塔各層中心位置的通用方法，并列表給出各次測(cè)量的古塔各層中心坐標(biāo)。
??2. 分析該塔傾斜、彎曲、扭曲等變形情況。
??3. 分析該塔的變形趨勢(shì)。

整體求解過程概述(摘要)

??本文要求根據(jù)測(cè)繪公司對(duì)古塔的4次測(cè)量數(shù)據(jù)，給出確定古塔各層中心位置的通用方法，并分析古塔的變形情況及其變形趨勢(shì)。為了計(jì)算的精度，我們首先對(duì)各變形量進(jìn)行了合理的數(shù)學(xué)定義，并對(duì)附錄的缺失數(shù)據(jù)進(jìn)行合理的賦值。
??對(duì)于問題一，我們通過最小二乘法擬合出觀測(cè)點(diǎn)所在平面，再建立優(yōu)化模型，在擬合平面上尋找到各觀測(cè)點(diǎn)距離的平方和最小的點(diǎn)作為古塔該層的中心點(diǎn)。利用MATLAB編程求解，得到了每次觀測(cè)古塔各層中心坐標(biāo)的通用方法及各層的中心點(diǎn)坐標(biāo)。
??對(duì)于問題二，我們將古塔的傾斜、彎曲和扭曲等變形情況，分別給予合理的數(shù)學(xué)描述。對(duì)于傾斜變形，我們定義了傾斜角α，即塔尖與底層中心的水平距離與塔高的比值；對(duì)于彎曲變形，我們定義了彎曲率K，即用中心點(diǎn)所擬合出的空間曲線的曲率來描述古塔各處彎曲率；對(duì)于扭曲變形，我們定義了相對(duì)扭曲度θ，利用坐標(biāo)的旋轉(zhuǎn)變換角度描述古塔的扭曲變形情況。利用空間曲線擬合、坐標(biāo)變換等方法以及MATLAB程序，分別求出了三個(gè)變形刻畫量的量化指標(biāo)。
??對(duì)于問題三，我們考慮通過古塔的傾斜、彎曲及扭曲程度來分析古塔的變形趨勢(shì)。由于數(shù)據(jù)量較少，我們建立灰色預(yù)測(cè)模型分析這三種變形因素的變化趨勢(shì)，利用相應(yīng)的MATLAB程序，得到了傾斜角、彎曲率以及相對(duì)扭曲度的預(yù)測(cè)函數(shù)和誤差檢驗(yàn)，驗(yàn)證了模型的可靠性，并繼而分析古塔的變形趨勢(shì)。
??本文巧妙地將各種變形量給予了合理的數(shù)學(xué)描述及模型，并運(yùn)用最小二乘法、曲線投影擬合、坐標(biāo)變換等數(shù)學(xué)方法實(shí)現(xiàn)了求解，并利用灰色預(yù)測(cè)對(duì)未來變形趨勢(shì)進(jìn)行了預(yù)測(cè)，具有較好的實(shí)用性和可推廣性。

模型假設(shè)：

??1.由于中國古代建筑物多為對(duì)稱圖形，假設(shè)古塔是對(duì)稱的。
??2.假設(shè)每次古塔的測(cè)量點(diǎn)選取是固定的。
??3.假設(shè)測(cè)量數(shù)據(jù)都是準(zhǔn)確可靠的。
??4.假設(shè)古塔的變形只由傾斜、彎曲和扭曲變形造成，不考慮其他因素。

問題分析：

??問題一要求確定古塔各層中心位置的通用方法。根據(jù)建筑變形測(cè)量規(guī)范，在建筑物變形測(cè)量中，為更好地測(cè)量出建筑物變形程度的各個(gè)指標(biāo)，我們假設(shè)每次測(cè)量應(yīng)選取固定的測(cè)量點(diǎn)，且在同一層所選取的測(cè)量點(diǎn)在未變形前處于同一個(gè)水平面上。而經(jīng)過對(duì)各層觀測(cè)點(diǎn)三維散點(diǎn)圖的繪制發(fā)現(xiàn)，各層的八個(gè)點(diǎn)近似對(duì)稱地分布在一個(gè)平面上，只是因?yàn)槟甏眠h(yuǎn)發(fā)生變形導(dǎo)致了些許偏差。因此為了更準(zhǔn)確地找出各層中心點(diǎn)，我們考慮先利用最小二乘法擬合出各層觀測(cè)點(diǎn)所在的平面方程，再建立優(yōu)化模型在該平面上尋找一點(diǎn)使其到各觀測(cè)點(diǎn)距離的平方和最小，以此確立古塔各層中心坐標(biāo)。

??問題二要求分析古塔的各種變形情況。根據(jù)《中華人民共和國行業(yè)標(biāo)準(zhǔn)建筑變形測(cè)量規(guī)范(JGJ8—2007)》知，變形是建筑的地基、基礎(chǔ)、上部結(jié)構(gòu)及其場地受各種作用力而產(chǎn)生的形狀或位置變化現(xiàn)象。在本問中，我們主要分析古塔三種主要的變形情況：傾斜、彎曲、扭曲。對(duì)于傾斜變形，我們定義傾斜角α進(jìn)行描述 ，其正切值等于塔尖與底層中心的水平距離與塔高的比值，即tanα=d/h ；對(duì)于彎曲變形，我們首先通過投影法擬合出古塔各層中心點(diǎn)所在空間曲線的參數(shù)方程，再利用空間曲線的曲率來刻畫古塔的彎曲度K；對(duì)于扭曲變形，考慮到扭曲變形實(shí)際為古塔水平面的旋轉(zhuǎn)產(chǎn)生，因此我們采用二維坐標(biāo)( x, y)旋轉(zhuǎn)的矩陣變換，通過各觀測(cè)量點(diǎn)前后的坐標(biāo)確定古塔的旋轉(zhuǎn)角度θ，以此刻畫古塔的扭曲度。但是，實(shí)際中水平面坐標(biāo)( x, y)不僅發(fā)生了旋轉(zhuǎn)變換，還受到傾斜彎曲變形等所引起的平移變化的影響，因此我們?cè)诳紤]坐標(biāo)變換的時(shí)候加入了平移量( p,q )使其更加準(zhǔn)確合理。

??問題三為分析古塔的變形情況。本文中，我們認(rèn)為建筑物變形由建筑物的傾斜、彎曲、扭曲等因素共同造成。由于附錄只給出了四次統(tǒng)計(jì)的數(shù)據(jù)，而我們的目標(biāo)是分析古塔未來多年的變化趨勢(shì)，因此我們采用信息不完全、不充分的預(yù)測(cè)系統(tǒng)——灰色預(yù)測(cè)對(duì)古塔未來的變形趨勢(shì)進(jìn)行預(yù)測(cè)。

模型的建立與求解整體論文縮略圖

全部論文請(qǐng)見下方“ 只會(huì)建模 QQ名片” 點(diǎn)擊QQ名片即可

clc,clear
x0=[0.0141,0.0142,0.0146,0.0147];%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差

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,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a1(1,:)=x0;
a1(2,:)=x0_hat;
a1(3,:)=epsilon;
a1(4,:)=delta;x0=x2(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a2(1,:)=x0;
a2(2,:)=x0_hat;
a2(3,:)=epsilon;
a2(4,:)=delta;x0=x3(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a3(1,:)=x0;
a3(2,:)=x0_hat;
a3(3,:)=epsilon;
a3(4,:)=delta;x0=x4(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a4(1,:)=x0;
a4(2,:)=x0_hat;
a4(3,:)=epsilon;
a4(4,:)=delta;x0=x5(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a5(1,:)=x0;
a5(2,:)=x0_hat;
a5(3,:)=epsilon;
a5(4,:)=delta;x0=x6(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a6(1,:)=x0;
a6(2,:)=x0_hat;
a6(3,:)=epsilon;
a6(4,:)=delta;x0=x7(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a7(1,:)=x0;
a7(2,:)=x0_hat;
a7(3,:)=epsilon;
a7(4,:)=delta;x0=x8(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a8(1,:)=x0;
a8(2,:)=x0_hat;
a8(3,:)=epsilon;
a8(4,:)=delta;x0=x9(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a9(1,:)=x0;
a9(2,:)=x0_hat;
a9(3,:)=epsilon;
a9(4,:)=delta;x0=x10(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a10(1,:)=x0;
a10(2,:)=x0_hat;
a10(3,:)=epsilon;
a10(4,:)=delta;x0=x11(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a11(1,:)=x0;
a11(2,:)=x0_hat;
a11(3,:)=epsilon;
a11(4,:)=delta;x0=x12(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a12(1,:)=x0;
a12(2,:)=x0_hat;
a12(3,:)=epsilon;
a12(4,:)=delta;x0=x13(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a13(1,:)=x0;
a13(2,:)=x0_hat;
a13(3,:)=epsilon;
a13(4,:)=delta;x0=x14(1,:);%原始數(shù)據(jù)序列
n=length(x0);
a_x0=diff(x0)';%求1次累減序列，即1階向前差分
B=[-x0(2:end)',ones(n-1,1)];
u=B\a_x0%最小二乘擬合參數(shù)
x=dsolve('D2x+a*Dx=b','x(0)=c1,Dx(0)=c2');%求二階微分方程的符號(hào)解
x=subs(x,{'a','b','c1','c2'},{u(1),u(2),x0(1),x0(1)});
yuce=subs(x,'t',0:n-1)%求已知數(shù)據(jù)點(diǎn)1次累加序列的預(yù)測(cè)值
x=vpa(x,6)
x0_hat=[yuce(1),diff(yuce)]%求已知數(shù)據(jù)點(diǎn)的預(yù)測(cè)值
epsilon=x0-x0_hat%求殘差
delta=abs(epsilon./x0)%求相對(duì)誤差
a14(1,:)=x0;
a14(2,:)=x0_hat;
a14(3,:)=epsilon;
a14(4,:)=delta;

全部論文請(qǐng)見下方“ 只會(huì)建模 QQ名片” 點(diǎn)擊QQ名片即可