matlab三次样条插值 代码

好的，以下是 MATLAB 的三次样条插值代码：

function [spline] = cubic_spline(x, y)
% x: 自变量向量
% y: 因变量向量
% spline: 三次样条插值函数

n = length(x);
h = diff(x);
alpha = zeros(n-2, 1);
for i = 2:n-1
    alpha(i-1) = 3/h(i-1)*(y(i+1)-y(i)) - 3/h(i)*(y(i)-y(i-1));
end

l = zeros(n, 1);
mu = zeros(n, 1);
z = zeros(n, 1);
l(1) = 1;
mu(1) = 0;
z(1) = 0;

for i = 2:n-1
    l(i) = 2*(x(i+1)-x(i-1))-h(i-1)*mu(i-1);
    mu(i) = h(i)/l(i);
    z(i) = (alpha(i-1)-h(i-1)*z(i-1))/l(i);
end

l(n) = 1;
z(n) = 0;
c = zeros(n, 1);
b = zeros(n-1, 1);
d = zeros(n-1, 1);

for j = n-1:-1:1
    c(j) = z(j)-mu(j)*c(j+1);
    b(j) = (y(j+1)-y(j))/h(j)-h(j)*(c(j+1)+2*c(j))/3;
    d(j) = (c(j+1)-c(j))/(3*h(j));
end

spline = @(xq) interp_spline(x, y, b, c, d, xq);

end

function [yq] = interp_spline(x, y, b, c, d, xq)
% x: 自变量向量
% y: 因变量向量
% b, c, d: 样条插值系数
% xq: 插值点向量
% yq: 插值点处的函数值

n = length(x);
nq = length(xq);
yq = zeros(nq,1);

for i = 1:nq
    j = find(x <= xq(i), 1, 'last');
    if j == n
        j = n-1;
    end
    h = xq(i)-x(j);
    yq(i) = y(j)+b(j)*h+c(j)*h^2+d(j)*h^3;
end

end