三次样条插值二类边界完整matlab代码

以下是三次样条插值二类边界的完整 MATLAB 代码示例：

% 三次样条插值二类边界完整matlab代码
% 输入数据：x, y，其中 x 和 y 分别为插值节点的横纵坐标
% 输出数据：pp，为一个结构体，包含三次样条插值函数的所有信息

% 数据准备
x = [0, 1, 2, 3, 4];
y = [1, 2.5, 3, 2.5, 1];

% 计算差分商
n = length(x);
h = diff(x);
delta = diff(y)./h;
l = [1, 2*h(1:end-1), 1];
mu = [0, h(2:end-1), 0];
z = [0, tridiag(l, mu, l, delta), 0];

% 计算三次样条插值函数的系数矩阵
M = zeros(n);
M(1, 1:2) = [1, 0];
M(n, n-1:n) = [0, 1];
for i = 2:n-1
    M(i, i-1:i+1) = [h(i-1), 2*(h(i-1)+h(i)), h(i)];
end

% 求解线性方程组，得到三次样条插值函数的系数向量
a = M\z';

% 构造三次样条插值函数的结构体
pp.form = 'pp';
pp.breaks = x;
pp.coefs = [a(1:n-1)', diff(a)./h(1:end-1)', 3*delta - 2*a(1:end-1)' - a(2:end)', (a(1:end-1)' + a(2:end)' - 2*delta)./h(1:end-1)'.^2];
pp.pieces = n-1;
pp.order = 4;
pp.dim = 1;
end

% 三对角线性方程组求解函数
function x = tridiag(a,b,c,r)
n = length(b);
x = zeros(n,1);
c(1) = c(1)/b(1);
r(1) = r(1)/b(1);
for i = 2:n-1
    temp = b(i) - a(i)*c(i-1);
    c(i) = c(i)/temp;
    r(i) = (r(i)-a(i)*r(i-1))/temp;
end
r(n) = (r(n)-a(n)*r(n-1))/(b(n)-a(n)*c(n-1));
x(n) = r(n);
for i = n-1:-1:1
    x(i) = r(i) - c(i)*x(i+1);
end
end

使用示例：

pp = spline2classbound(x, y);
xx = linspace(0, 4, 100);
yy = ppval(pp, xx);
plot(x, y, 'o', xx, yy);