splint1.m

来自「用matlab写的一些数值算法」· M 代码 · 共 39 行

function  p = splint1(x,f,df)
% Compute coefficients for interpolating cubic spline
% Call    p = splint1(x,f)     % Natural spline
% or      p = splint1(x,f,df)  % Correct boundary conds.
% Input
% x:  Vector with knots.  Length n+1.
% f:  Vector with nodal function values.  Length n+1.
% df: If present, then a vector with endpoint derivatives,
%     df(1) = f'(x(1)),  df(2) = f'(x(n+1))
% Output
% p : n*4 array with p(i,:) = [ai bi ci di].

% Version 3.06.2004.  INCBOX 

n = length(x) - 1;     % number of knot intervals
h = diff(x);           % n-vector with knot spacings
v = diff(f)./h;        % n-vector with divided differences
% Set up matrix A and right hand side r
A = zeros(n+1,n+1);   r = zeros(n+1,1);
for  i = 1 : n-1
  A(i+1,i:i+2) = [h(i+1)  2*(h(i)+h(i+1))  h(i)];
  r(i+1) = 3*(h(i+1)*v(i) + h(i)*v(i+1));
end
if  nargin == 3                  % Correct boundary conds.
  A(1,1) = 1;  r(1) = df(1);
  A(n+1,n+1) = 1;  r(n+1) = df(2);
else                             % Natural spline.  
  A(1,1:2) = [2 1];  r(1) = 3*v(1);
  A(n+1,n:n+1) = [1 2];  r(n+1) = 3*v(n);
end
ds = A\r;                        % Solve  A*ds = r
% Compute coefficients
p = zeros(n,4);
for  i = 1:n
  p(i,1) = f(i);
  p(i,2) = h(i)*ds(i);
  p(i,3) = 3*(f(i+1) - f(i)) - h(i)*(2*ds(i) + ds(i+1));
  p(i,4) = 2*(f(i) - f(i+1)) + h(i)*(ds(i) + ds(i+1));
end