自行编程第一种边界条件三次样条插值matlab

下面是一个使用Matlab编写的三次样条插值程序，它使用第一种边界条件（即在第一个和最后一个点上的一阶导数分别为给定值）：

function [s] = cubic_spline(x,y,dy1,dyn)
% x,y: 数据点的横纵坐标
% dy1,dyn: 第一种边界条件下，第一个和最后一个点上的一阶导数值

n = length(x);
h = diff(x);
delta = diff(y) ./ h;
d = [3*(delta(1)-dy1); 3*(delta(2:n-1)-delta(1:n-2)); 3*(dyn-delta(n-1))];
A = spdiags([h(1:n-2) 2*(h(1:n-2)+h(2:n-1)) h(2:n-1)], [-1 0 1], n-2, n-2);
m = A \ d;

% 构造插值函数
s = cell(1,n-1);
syms t;
for i = 1:n-1
    if i == 1
        s{i} = simplify(y(1) + delta(1)*(t-x(1)) + (3*m(1)-2*delta(1)-m(2))*(t-x(1))^2/h(1) + (m(1)-delta(1)-2*m(2))/h(1)*(t-x(1))^3);
    elseif i == n-1
        s{i} = simplify(y(n) + delta(n-1)*(t-x(n)) + (3*m(n-2)-2*delta(n-1)-m(n-3))*(t-x(n))^2/h(n-1) + (m(n-2)-delta(n-1)-2*m(n-3))/h(n-1)*(t-x(n))^3);
    else
        s{i} = simplify(y(i) + delta(i-1)*(t-x(i)) + (3*m(i-1)-2*delta(i-1)-m(i))*(t-x(i))^2/h(i-1) + (m(i-1)-delta(i-1)-2*m(i))/h(i-1)*(t-x(i))^3);
    end
end

% 画出插值函数
fplot(s{1}, [x(1) x(n)]);
hold on;
for i = 2:n-1
    fplot(s{i}, [x(i) x(i+1)]);
end
scatter(x, y);
end

该程序使用了Matlab内置的spdiags函数来构造系数矩阵，并使用了syms函数来构造插值函数。在程序中，输入数据点的横纵坐标存储在向量x和y中，第一个和最后一个点上的一阶导数分别为dy1和dyn。程序返回一个cell数组s，其中每个元素都是一个符号表达式，表示在对应区间上的插值函数。