www.pudn.com > naknot.rar > naknot.m, change:2009-10-25,size:1446b


%%%%%% div为等分数目 
function naknot(X,Y,div) 
%%%%%% X,Y中的下标从1开始 
n=length(X); 
a=(X(2)-X(1))/(X(3)-X(2)); 
b=(X(n)-X(n-1))/(X(n-1)-X(n-2)); 
%%%%%% 计算u,v,d 
for k=2:(n-1) 
    u(k)=(X(k)-X(k-1))/(X(k+1)-X(k-1)); 
    v(k)=1-u(k); 
    d(k)=((Y(k+1)-Y(k-1))/(X(k+1)-X(k-1))-(Y(k)-Y(k-1))/(X(k)-X(k-1)))/(X(k+1)-X(k)); 
end 
d=d'; 
%%%%%% 计算p,p为三对角阵 
p(1,1)=2+a; 
p(1,2)=1-a; 
p(n-2,n-3)=1-b; 
p(n-2,n-2)=2+b; 
for k=2:(n-3) 
    p(k,k-1)=u(k+1); 
    p(k,k)=2; 
    p(k,k+1)=v(k+1); 
end 
%%%%%% 计算M 
M(2:n-1)=inv(p)*(d(2:n-1).*6); 
M(1)=(M(2)-u(2)*M(3))/v(2); 
M(n)=(M(n-1)-v(n-1)*M(n-2))/u(n-1); 
%%%%%% 计算每段样条系数s 
for k=1:(n-1) 
    s(k,1)=(M(k+1)-M(k))/(6*(X(k+1)-X(k))); 
    s(k,2)=M(k)/2; 
    s(k,3)=(Y(k+1)-Y(k))/(X(k+1)-X(k))-(2*M(k)+M(k+1))*(X(k+1)-X(k))/6; 
    s(k,4)=Y(k); 
end 
%%%%%% 绘图 
x=linspace(X(1),X(n),(div+1)); 
%%%%%% deltax为间隔,beginn数组的起始,endn数组的结束,deltan数组每次增加的数量 
deltax=(X(n)-X(1))/(div); 
beginn=1; 
deltan=fix((X(2)-X(1))/deltax)+1; 
endn=deltan; 
y(beginn:endn)=s(1,1).*((x(beginn:endn)-X(1)).^3)+s(1,2).*((x(beginn:endn)-X(1)).^2)+s(1,3).*((x(beginn:endn)-X(1)))+s(1,4); 
for k=2:(n-1) 
    beginn=endn+1; 
    deltan=fix((X(k+1)-X(1))/deltax)-fix((X(k)-X(1))/deltax); 
    endn=beginn+deltan-1; 
    y(beginn:endn)=s(k,1).*((x(beginn:endn)-X(k)).^3)+s(k,2).*((x(beginn:endn)-X(k)).^2)+s(k,3).*((x(beginn:endn)-X(k)))+s(k,4); 
end 
figure 
plot(X,Y,'o',x,y)