Matlab实现——Clamped Cubic Spline

欢迎前往个人博客 驽马点滴 和视频空间 哔哩哔哩-《挨踢日志》

%csfit.m

%Program 5.3 (Clamped Cubic Spline ).

function S=csfit(X,Y,dx0,dxn)

%Input - X is the 1xn abscissa vector

%           - Y is the 1xn ordinate vector

%       - dxo = S'(x0) first derivative boundary condition

%       - dxn = S'(xn) first derivative boundary condition

%Output     - S: rows of S are the coefficients for the  

%                     cubic interpolants

N=length(X)-1;

H=diff(X);

D=diff(Y)./H;

A=H(2:N-1);

B=2*(H(1:N-1)+H(2:N));

C=H(2:N);

U=6*diff(D);

%Clamped spline endpoint constraints

B(1)=B(1)-H(1)/2;

U(1)=U(1)-3*(D(1)-dx0);

B(N-1)=B(N-1)-H(N)/2;

U(N-1)=U(N-1)-3*(dxn-D(N));

for k=2:N-1

   temp=A(k-1)/B(k-1);   % multiple

   B(k)=B(k)-temp*C(k-1); % elimination

   U(k)=U(k)-temp*U(k-1); % elimination 

end

M(N)=U(N-1)/B(N-1);   % back substitution

for k=N-2:-1:1

   M(k+1)=(U(k)-C(k)*M(k+2))/B(k);  

end

M(1)=3*(D(1)-dx0)/H(1)-M(2)/2;

M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;

for k=0:N-1

   S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));

   S(k+1,2)=M(k+1)/2;

   S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;

   S(k+1,4)=Y(k+1);

end

欢迎前往个人博客 驽马点滴 和视频空间 哔哩哔哩-《挨踢日志》