matlab三次样条插值代码

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

function [spline] = cubicSpline(x, y, xq)
% x: 源数据的x坐标
% y: 源数据的y坐标
% xq: 查询点的x坐标
% spline: 查询点的y坐标

% 计算差分
n = length(x);
h = diff(x);
delta = diff(y)./h;

% 计算三对角矩阵的上下对角线和右侧的向量
u = [0; h(1:n-2)];
v = [h(2:n-1); 0];
d = [1; 2*(h(1:n-2)+h(2:n-1)); 1];
b = [0; 3*(delta(2:n-1).*h(1:n-2)+delta(1:n-2).*h(2:n-1)); 0];

% 解三对角矩阵
c = ThomasAlgorithm(u,d,v,b);

% 计算系数
a = y(1:n-1);
b = delta;
d = (c(2:n)-c(1:n-1))./(3*h);

% 计算样条插值
spline = zeros(size(xq));
for i = 1:length(xq)
    if xq(i)<x(1)
        spline(i) = y(1);
    elseif xq(i)>x(n)
        spline(i) = y(n);
    else
        j = find(xq(i)>=x,1,'last');
        spline(i) = a(j)+b(j)*(xq(i)-x(j))+c(j)*(xq(i)-x(j))^2+d(j)*(xq(i)-x(j))^3;
    end
end
end

% 解三对角矩阵的函数
function [x] = ThomasAlgorithm(a,b,c,d)
% a: 下对角线
% b: 主对角线
% c: 上对角线
% d: 右侧向量
% x: 解

n = length(d);
c(1) = c(1)/b(1);
d(1) = d(1)/b(1);

for i = 2:n-1
    c(i) = c(i)/(b(i)-a(i)*c(i-1));
    d(i) = (d(i)-a(i)*d(i-1))/(b(i)-a(i)*c(i-1));
end

d(n) = (d(n)-a(n)*d(n-1))/(b(n)-a(n)*c(n-1));

x(n) = d(n);
for i = n-1:-1:1
    x(i) = d(i)-c(i)*x(i+1);
end
end

使用方法：

假设有以下源数据：

x = [1 2 3 4 5];
y = [3 4 2 1 2];

要在查询点`xq`处进行插值，可以使用以下代码：

xq = linspace(1,5,100);
spline = cubicSpline(x,y,xq);
plot(x,y,'o',xq,spline);

这将绘制出原始数据点和使用三次样条插值算法得到的插值曲线。