poly_nihe.m

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

function [a,b,c,d,cigema_squre]=poly_nihe(x,y,phasenum)
%设样条函数为a1*x^3+a2*x^2+a3*x+a4=y
num=fix((length(x)-1)/phasenum);   %平均分的话，每段所占有的点数
rest_num=(length(x)-1)-phasenum*num;%平均分剩下的点;
for i=1:phasenum
    if rest_num~=0
        phase_num(i)=num+1;%每段享有的点数
        rest_num=rest_num-1;
    else 
        phase_num(i)=num;
    end
end
if phasenum==1%只有一段样条＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
		A=[];
		L=[];
		for k=1:length(x)
            A(k,1)=x(k)^3;
            A(k,2)=x(k)^2;
            A(k,3)=x(k);
            A(k,4)=1;
            L(k)=y(k);
		end
            temp=inv(A'*A)*A'*L';
            a=temp(1);
            b=temp(2);
            c=temp(3);
            d=temp(4);
elseif phasenum>1%n次样条＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
	A=[];
	L=[];
	linepub=1;
    sum_num=0;
	for i=1:length(x)      
        if i==sum_num+phase_num(linepub)+1 & i~=length(x) %判断该点在哪一段内
            sum_num=sum_num+phase_num(linepub);
            linepub=linepub+1;
        end
        A(i,linepub*4-3)=x(i)^3;
        A(i,linepub*4-2)=x(i)^2;
        A(i,linepub*4-1)=x(i);
        A(i,linepub*4)=1;
        L(i)=y(i);
	end
%==============================写条件式=================================
	for k=1:(phasenum-1)
        total=0;
        for t=1:k
            total=total+phase_num(t);
        end
        continue_point_num=total+1;
        C(k,k*4-3)=x(continue_point_num)^3;  %连续条件  
        C(k,k*4-2)=x(continue_point_num)^2;
        C(k,k*4-1)=x(continue_point_num);
        C(k,k*4)=1;    
        C(k,k*4+1)=-x(continue_point_num)^3;
        C(k,k*4+2)=-x(continue_point_num)^2;
        C(k,k*4+3)=-x(continue_point_num);
        C(k,k*4+4)=-1;     
        
        C(k+(phasenum-1),k*4-3)=3*x(continue_point_num)^2;  %一阶导数连续
        C(k+(phasenum-1),k*4-2)=2*x(continue_point_num);
        C(k+(phasenum-1),k*4-1)=1;
        C(k+(phasenum-1),k*4)=0;
        C(k+(phasenum-1),k*4+1)=-3*x(continue_point_num)^2;
        C(k+(phasenum-1),k*4+2)=-2*x(continue_point_num);
        C(k+(phasenum-1),k*4+3)=-1;
        C(k+(phasenum-1),k*4+4)=0;
        
        C(k+2*(phasenum-1),k*4-3)=6*x(continue_point_num);            %二阶导数连续
        C(k+2*(phasenum-1),k*4-2)=2;
        C(k+2*(phasenum-1),k*4-1)=0;
        C(k+2*(phasenum-1),k*4)=0;
        C(k+2*(phasenum-1),k*4+1)=-6*(continue_point_num);
        C(k+2*(phasenum-1),k*4+2)=-2;
        C(k+2*(phasenum-1),k*4+3)=0;
        C(k+2*(phasenum-1),k*4+4)=0;  
	end
	C(1+3*(phasenum-1),1)=6*x(1);   %端点条件　从第1点开始
	C(1+3*(phasenum-1),2)=2;
	C(1+3*(phasenum-1),3)=0;
	C(1+3*(phasenum-1),4)=0;
	
 	C(2+3*(phasenum-1),(phasenum-1)*4+1)=6*x(length(x));   %端点条件　最后点
 	C(2+3*(phasenum-1),(phasenum-1)*4+2)=2;
 	C(2+3*(phasenum-1),(phasenum-1)*4+3)=0;
 	C(2+3*(phasenum-1),(phasenum-1)*4+4)=0;
	
	W(3*(phasenum-1)+2)=0;
	Na=A'*A;
	U=A'*L';
	temp=inv([Na C';C zeros(3*(phasenum-1)+2,3*(phasenum-1)+2)])*[U;W'];
	for i=1:phasenum
        a(i)=temp(i*4-3);
        b(i)=temp(i*4-2);
        c(i)=temp(i*4-1);
        d(i)=temp(i*4);
	end
end
if phasenum>1
    V=A*temp(1:(length(temp)-3*(phasenum-1)-2))-L';
else
    V=A*temp(1:(length(temp)-3*(phasenum-1)))-L';
end
    cigema_squre=sqrt(V'*V/(length(x)))*1000