Spline3_inter.m

% 函数功能：三次样条插值函数，第三种边界条件：已知函数是以T=b-a为周期的周期函数
% 输入：n+1个插值点列向量x（从0到n）和y （这里首尾都输入了，实际上重复了）
% 输出：分段插值函数S
function S = Spline3_inter(x, y)
%% 三次样条插值的M值的方程组的导出
n = length(x);
h = diff(x); % 这里定义的h1 = x2 - x1
h = [h; h(1)];
y = [y; y(2)]; % 这两行是第三种边界条件特地加的
u = h(1: end - 1) ./ (h(1: end - 1) + h(2: end)); % n-1维列向量
la = 1 - u; % n-1维列向量
f2 = (diff(y(2: end)) ./ h(2: end) - diff(y(1: end - 1)) ./ h(1: end - 1)) ./ (h(1: end - 1) + h(2: end)); % 二阶差商
d = 6 * f2;

%% 求关于M的方程组的系数矩阵
Mu = eye(n - 1); % 存放u元素的矩阵
Mu(Mu == 1) = u;
Mu = Mu(:, [2: end, 1]);
L = eye(n - 1); % 存放la元素的矩阵
L(L == 1) = la;
L = L(:, [end, 1: end - 1]);
A = 2 * eye(n - 1) + Mu + L;

%% 不是三对角方程组标准形式，用LU分解求解（因为不想为一个特殊矩阵编写特殊算法，牺牲了速度）
M = LU_equ(A, d); % LU分解法
M = [M(end); M]; % 周期函数的首尾元素相等，补上第一个元素值

%% 代入S的方程组，S(x)为分段函数，段数不确定，用循环表示
S = cell(n - 1, 2); % 第一列存放插值表达式，第二列存放区间端点
for i = 1: n - 1
    %% 代入S的表达式
    S{i, 1} = @(xx)( ...
        (x(i + 1) - xx) .^3 ./ (6 .* (x(i + 1) - x(i))) .* M(i) ...
        + (xx - x(i)) .^3 ./ (6 .* (x(i + 1) - x(i))) .* M(i + 1) ...
        + (y(i) - (x(i + 1) - x(i)) .^2 / 6 .* M(i)) .* (x(i + 1) - xx) ./ (x(i + 1) - x(i)) ...
        + (y(i + 1) - (x(i + 1) - x(i)) .^2 ./ 6 .* M(i + 1)) .* (xx - x(i)) ./ (x(i + 1) - x(i)) ...
        ) .* (xx >= x(i) & xx <= x(i + 1)); % xx的定义域
    S{i, 2} = [x(i), x(i + 1)];
    
    %% 绘图
    fplot(S{i, 1}, S{i, 2}); % 函数句柄中的运算符、逻辑符需要转换成向量适用的(./ .* &)
    hold on
end

end