paper.m

这是用于工程测量进行大坝变形分析的软件

clear all
clc
[filename1,pathname1]=uigetfile('*.txt','txt文本');
fit1=fopen(strcat(pathname1,filename1),'rt');
if(fit1==-1)
    msgbox('Input File or Path is not correct','Warning','warn');
    break;
end

[filename2,pathname2]=uiputfile('*.txt','txt文本');
fit2=fopen(strcat(pathname2,filename2),'wt');
if(fit2==-1)
    msgbox('Error by Opening Output File','Warning','warn');
    break;
end
No=fscanf(fit1,'%d',1);  %读取点的个数
j=0;  %是日期的编号
while (1) 
    temp=fscanf(fit1,'%d',1);
    if temp==''
        break;
    end
    j=j+1;
    year(j)=temp;
    month(j)=fscanf(fit1,'%d',1);
    day(j)=fscanf(fit1,'%d',1);
    for i=1:No                                     
        x(j,i)=fscanf(fit1,'%f',1);

        y(j,i)=fscanf(fit1,'%f',1);
        if j~=1 
            u(j,i)=x(j,i)-x(j-1,i);
            v(j,i)=y(j,i)-y(j-1,i);
        end
    end  
end                      %数据读取完成
time=j;                   %把观测的次数送给time


for j=2:time

	for i=2:No
        A((i-1)*2-1,1)=x(j,i)-x(j,1);
        A((i-1)*2-1,2)=0;
        A((i-1)*2-1,3)=y(j,i)-y(j,1);
        %A((i-1)*2-1,4)=y(j,i)-y(j,1);
        
        A((i-1)*2,1)=0;
        A((i-1)*2,2)=y(j,i)-y(j,1);
        A((i-1)*2,3)=x(j,i)-x(j,1);
        %A((i-1)*2,4)=x(j,i)-x(j,1); 
        
        L((i-1)*2-1)=u(j,i)-u(j,1);
        L((i-1)*2)=v(j,i)-v(j,1);
        
	end
	temp=inv(A'*A)*A'*L'
	E11(j)=temp(1);
	E22(j)=temp(2);
	%E12(j)=temp(3);
	w(j)=temp(3);
end



%画出E11各点
for k=2:time
    hold on
    plot(k,E11(k))
end



%对E11进行直线拟合
%设直线为y=ax+b
A=[];
L=[];
for k=2:time
    A(k-1,1)=k;
    A(k-1,2)=1;
    L(k-1)=E11(k);
end
    temp=inv(A'*A)*A'*L';
    a=temp(1);
    b=temp(2);
%算出E11与直线拟合值的差
for k=2:time
    dE11(k)=a*k+b-E11(k);
end
%画出拟合好的直线
plot([2:time],a.*[2:time]+b) 



%对E11进行三次样条拟合 
%设样条函数为a1*x^3+a2*x^2+a3*x+a4=y
A=[];
L=[];
for k=2:time
    A(k-1,1)=k^3;
    A(k-1,2)=k^2;
    A(k-1,3)=k;
    A(k-1,4)=1;
    
    L(k-1)=E11(k);
end
    temp=inv(A'*A)*A'*L';
    a1=temp(1);
    b1=temp(2);
    c1=temp(3);
    d1=temp(4);
    
%算出E11与直线拟合值的差
for k=2:time
    dE11(k)=a1*k^3+b1*k^2+c1*k+d1-E11(k);
end
%画出拟合好的直线
plot([2:time],a1.*[2:time].^3+b1.*[2:time].^2+c1.*[2:time]+d1) 



%对E11进行任意段三次样条拟合
n=input('请输入要拟合的样条函数分的段数值，（不得大于观次数－2）')

num=fix((time-1)/n);   %平均分的话，每段所占有的点数
A=[];
L=[];
linepub=1;
for i=2:time
    %for j=1:n*4
        A(i-1,linepub*4-3)=i^3;
        A(i-1,linepub*4-2)=i^2;
        A(i-1,linepub*4-1)=i;
        A(i-1,linepub*4)=1;
        
        
    
        L(i-1)=E11(i);
        
        if ((i-1)/num)==fix((i-1)/num) & (i-1)<n*num %判断如果点号可以整除每段共有的点数，而且不是最后一段
            linepub=linepub+1;
        end
end

for k=1:(n-1)
    C(k,k*4-3)=(k*num+1+1)^3;  %连续条件   (注：因为第一个点是0所以每3点实际上是第4点
    C(k,k*4-2)=(k*+num+1+1)^2;
    C(k,k*4-1)=(k*+num+1+1);
    C(k,k*4)=1;    
    C(k,k*4+1)=-(k*num+1+1)^3;
    C(k,k*4+2)=-(k*+num+1+1)^2;
    C(k,k*4+3)=-(k*+num+1+1);
    C(k,k*4+4)=-1;     
    
    C(k+(n-1),k*4-3)=3*(k*num+1+1)^2;            %一阶导数连续
    C(k+(n-1),k*4-2)=2*(k*num+1+1);
    C(k+(n-1),k*4-1)=1;
    C(k+(n-1),k*4)=0;
    C(k+(n-1),k*4+1)=3*(k*num+1+1)^2;
    C(k+(n-1),k*4+2)=2*(k*num+1+1);
    C(k+(n-1),k*4+3)=1;
    C(k+(n-1),k*4+4)=0;
    
    C(k+2*(n-1),k*4-3)=6*(k*num+1+1);            %二阶导数连续
    C(k+2*(n-1),k*4-2)=2;
    C(k+2*(n-1),k*4-1)=0;
    C(k+2*(n-1),k*4)=0;
    C(k+2*(n-1),k*4+1)=6*(k*num+1+1);
    C(k+2*(n-1),k*4+2)=2;
    C(k+2*(n-1),k*4+3)=0;
    C(k+2*(n-1),k*4+4)=0;  
       
end

C(1+3*(n-1),1)=6*2;   %端点条件　从第2点开始
C(1+3*(n-1),2)=2;
C(1+3*(n-1),3)=0;
C(1+3*(n-1),4)=0;

C(2+3*(n-1),(n-1)*4+1)=6*time;   %端点条件　最后点
C(2+3*(n-1),(n-1)*4+2)=2;
C(2+3*(n-1),(n-1)*4+3)=0;
C(2+3*(n-1),(n-1)*4+4)=0;

W(3*(n-1)+2)=0;

Na=A'*A;
U=A'*L';

temp=inv([Na C';C zeros(3*(n-1)+2,3*(n-1)+2)])*[U;W']

for i=1:n
    a2(i)=temp(i*4-3);
    b2(i)=temp(i*4-2);
    c2(i)=temp(i*4-1);
    d2(i)=temp(i*4);
end
%画出拟合好的直线
for i=1:n
    plot([2:i*num+1+1],a2(i).*[2:i*num+1+1].^3+b2(i).*[2:i*num+1+1].^2+c2(i).*[2:i*num+1+1]+d2(i)) ;

end



% for k=
%     lamda(k)=h(k+1)/(h(k)+h(k+1));
%     mu(k)=1-lamda(k);
%     d(k)=6/(h(k)+h(k+1))*((E11(k+1)-e11(k))/h(k+1)-(E11(k)-E(k-1))/h(k))
% end
% A=[];
% L=[];
% 
% A(1,1)=2;
% A(1,2)=1;
% L(1)=d(1);
% for i=2:time-1
%     A(i,i-1)=mu(i-1);
%     A(i,i)=2;
%     A(i,i+1)=lamda(i-1);
%     L(i)=d(i);
%     
% end
% 
% A(time,time-1)=1;
% A(time,time)=2;
% L(time)=d(time);
% 
% temp1=inv(A)*L;