sf_spline.m

来自「基于matlab的三次样条插值算法程序」· M 代码 · 共 51 行

function [S] = SF_Spline(x,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 三次样条 插值多项式的程序
%
% x    插值节点 向量
% y    插值点上的函数值 向量
%
% 返回 分段 3 次 Hermite 插值多项式 Sn, 系数存储在矩阵 S 的每一行中
%
% Oct 2, 2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           
n = size(x,2);              % 插值点个数, 注意: 不是书上的 n
for j = 1 : n-1
    h(j) = x(j+1) - x(j);             % 步长
    f(j) = (y(j+1) - y(j)) / h(j);    % 一阶差商
end    
for j = 2 : n-1
    lambda(j) = h(j)/(h(j-1) + h(j));
    mu(j) = h(j-1)/(h(j-1) + h(j));
    g(j) = 3 * (lambda(j) * f(j-1) + mu(j) * f(j));
end

A = 2 * eye(n,n);
for i = 2 : n-1
    A(i,i+1) = mu(i);
    A(i,i-1) = lambda(i);
end

b = zeros(n,1);
for i = 2 : n-1
    b(i) = g(i);
end

% 2 阶导数边界条件
b(1) = 3*f(1);        
A(1,2) = 1;
b(n) = 3*f(n-1);       
A(n,n-1) = 1;

m = A \ b; % 求解线性方程组, 得到各节点的一阶导数值
S = [];
for j = 1 : n-1
    s = (y(j)/(h(j)^3))*conv([1 -2*x(j+1) x(j+1)^2],[2 h(j)-2*x(j)]) + ...
    (y(j+1)/(h(j)^3))*conv([1 -2*x(j) x(j)^2],[-2 h(j)+2*x(j+1)]) + ...
    (m(j)/(h(j)^2))*conv([1 -2*x(j+1) x(j+1)^2],[1 -x(j)]) + ...
    (m(j+1)/(h(j)^2))*conv([1 -2*x(j) x(j)^2],[1 -x(j+1)]);
    S = [S;s]; 
end