Matlab三次样条插值程序

以下是Matlab的三次样条插值程序示例：

function [S] = cubic_spline(x, y, xx)
% x: vector of x coordinates of data points
% y: vector of y coordinates of data points
% xx: vector of x coordinates of points to be interpolated
% S: vector of interpolated y values at points xx

n = length(x);
h = zeros(n-1,1);
alpha = zeros(n-1,1);
l = zeros(n,1);
mu = zeros(n-1,1);
z = zeros(n,1);
c = zeros(n,1);
b = zeros(n-1,1);
d = zeros(n-1,1);

% calculate h
for i = 1:n-1
    h(i) = x(i+1) - x(i);
end

% calculate alpha
for i = 2:n-1
    alpha(i) = 3/h(i)*(y(i+1)-y(i)) - 3/h(i-1)*(y(i)-y(i-1));
end

% calculate l, mu, z
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) - h(i-1)*z(i-1))/l(i);
end

l(n) = 1;
z(n) = 0;
c(n) = 0;

% calculate c
for i = n-1:-1:1
    c(i) = z(i) - mu(i)*c(i+1);
end

% calculate b and d
for i = 1:n-1
    b(i) = (y(i+1)-y(i))/h(i) - h(i)*(c(i+1)+2*c(i))/3;
    d(i) = (c(i+1)-c(i))/(3*h(i));
end

% interpolate at xx
S = zeros(size(xx));
for j = 1:length(xx)
    if xx(j) < x(1) || xx(j) > x(n)
        error('xx(j) is out of range')
    end
    i = 1;
    while xx(j) > x(i+1)
        i = i+1;
    end
    S(j) = y(i) + b(i)*(xx(j)-x(i)) + c(i)*(xx(j)-x(i))^2 + d(i)*(xx(j)-x(i))^3;
end
end

在上述代码中，我们首先计算了样条插值所需的中间变量（如$l, \mu, z$）和系数（如$b, c, d$），然后使用这些变量和系数进行插值。要使用此函数进行插值，请将数据点的x坐标存储在一个向量x中，y坐标存储在向量y中，要插值的点的x坐标存储在向量xx中，然后调用函数cubic_spline(x, y, xx)即可。