第一种三次样条函数的Matlab求解

来源：李庆杨等第五版《数值分析》P44页例题。目的：回顾一下matlab的基本用法。

首先根据书中给出的数学表达式编写求解弯矩（书上这么叫的）函数：

function [M,A,d,H] = SolveSpline_1(x,y,diff_y0,diff_yn)
    %Function：求解第一种三次样条函数插值
    %变量说明：
	%Input varibles:
		% x：自变量数据（以列的形式输入）
		% y：因变量数据（以列的形式输入）
		% diff_y0：左端点导数值
		% diff_yn：右端点到数值
	%Output varibles:
		%M：三次样条函数中的弯矩
		%A：三对角矩阵
		%d：右端向量
		%H：节点步长

    xx = x(2:end);
    a = x(1);%作为x0点
    [n,m] = size(xx);

    A = 2*eye(n+1,n+1);%初始化三对角矩阵

    H = diff(x);
    %对角参数设置
    mu = zeros(n,1);
    mu(n) = 1;
    for k = [1:n-1]
        mu(k) = H(k)/(H(k) + H(k+1));
    end

    lambda = zeros(n,1);
    lambda(1) = 1;
    for k = [2:n]
        lambda(k) = H(k)/(H(k-1) + H(k)); 
    end

    d = zeros(n+1,1);
    %temp = (y(2) - y(1))/(x(2) - x(1))
    d(1) = (6/H(1))*((y(2) - y(1))/(x(2) - x(1)) - diff_y0);
    %temp = diff_yn - ((y(end) - y(end-1))/(x(end) - x(end-1)))
    d(n+1) = (6/H(n))*(diff_yn - ((y(end) - y(end-1))/(x(end) - x(end-1))));
    for k = [2:n]
        %temp = (y(k+1) - y(k))/(x(k+1) - x(k))
        %temp = (y(k) - y(k-1))/(x(k) - x(k-1))
        d(k) = 6*((y(k+1) - y(k))/(x(k+1) - x(k)) - (y(k) - y(k-1))/(x(k) - x(k-1)))/(x(k+1) - x(k-1));
    end
    
    %给出三对角矩阵
    A = A + diag(mu,-1) + diag(lambda,1);
    %求解M
    M = A^(-1)*d;
end

%这样求解完 M 后，便可将 M 带入到三次样条函数表达式中去，从而便得出了小区间上的三次样条插值函数

再编写主函数并绘制三次样条插值函数的图像，如下：

clc,clear all
x = [27.7,28,29,30]';
y = [4.1,4.3,4.1,3.0]';
diff_y0  = 3.0;
diff_yn = -4.0;
[M,A,d,H] = SolveSpline_1(x,y,diff_y0,diff_yn);

%利用分段函数的写法给出三次样条插值函数
[n,m] = size(x);
n = n-1;
t = [27.7:0.001:30];
S = 0;
for k = [1:n]

%{
    temp = (x(k+1) - t).^3./(6*H(k))
    temp = (t - x(k)).^3./(6*H(k))
    temp = y(k) - M(k)*H(k)^2/6
    temp = (x(k+1) - t)/H(k)
    temp = y(k+1) - M(k+1)*H(k)^2/6
    temp = x(k+1) - t
    temp = (t - x(k+1))./H(k) 
%}
    
    S = S + (M(k)*((x(k+1) - t).^3./(6*H(k))) + M(k+1)*((t - x(k)).^3./(6*H(k))) ... 
        + (y(k) - M(k)*H(k)^2/6).*((x(k+1) - t)/H(k)) + (y(k+1) - M(k+1)*H(k)^2/6).*((t - x(k))./H(k))).*(t > x(k) & t < x(k+1));
end

figure
title('Demo about spline');
hold on
plot(t,S);
scatter(x,y)
hold off
xlabel('X');
ylabel('Y');
legend('spline function','previous points')
grid on;

最后，简单的绘图结果如下：