分段三次Hermite插值Matlab实现

function  [m_matrix]=hermite(x,y,y0,yn,x_value)
%% 输入值分配，x_input,y_input均为数组,y_0,y_n为x_0,x_n分别对应的一阶导数值
x_input = x;
y_input = y;
y_0 = y0;
y_n = yn;
%% 
[~,number] = size(x_input); %获取输入x_input的大小1×number
delta_h = zeros(1,number-1); %给delta_h分配数组大小1×(number-1)，并全部初始化为0
delta_f = zeros(1,number-1); %给delta_f分配数组大小1×(number-1)，并全部初始化为0
lambda_ = zeros(1,number-2); %给lambda_分配数组大小1×(number-2)，并全部初始化为0
miu = zeros(1,number-2); %给miu_分配数组大小1×(number-2)，并全部初始化为0
e = zeros(1,number-2); %给e分配数组大小1×(number-2)，并全部初始化为0
% 计算delta_h、delta_f的值
for i = 1:(number-1)
    delta_h(i) = x_input(i+1) - x_input(i);
    delta_f(i) = (y_input(i+1) - y_input(i))/ delta_h(i);
end

%%计算lambda,miu,e
for i=1:number-2
    lambda_(1,i) = delta_h(1,i+1) / (delta_h(1,i+1) + delta_h(1,i));
     miu(1,i) = 1 - lambda_(1,i);
     e(1,i) = 3*(lambda_(1,i)*delta_f(1,i) + miu(1,i)*delta_f(1,i+1));
end

A = zeros(number-2,number-2); %初始化系数矩阵A，(n-2)*(n-2)
B = zeros(number-2,1); %初始化系数矩阵B，(n-2)*1
%当i=1时
A(1,1) = 2;
A(1,2) = miu(1,1);
B(1,1) = e(1,1) - lambda_(1,1) * y_0;

%当i=2:n-2时
for i = 2:number-3
    B(i,1) = e(1,i);
    A(i,i-1) = lambda_(1,i);
    A(i,i) = 2;
    A(i,i+1) = miu(1,i);
end
%当i=n-1时
A(number-2,number-3) = lambda_(1,number-2);
A(number-2,number-2) = 2;
B(number-2,1) = e(1,number-2) - miu(1,number-2)*y_n;

%% 计算A的逆，A*B
m_matrix = A\B;
m = zeros(1,number);
m(1,1) = y_0;
m(1,number) = y_n;
for i = 2:number-1
    m(1,i) = m_matrix(i-1,1);
end

for i =1:number-1
    % 获取相邻两点间的插值函数
    x_ = linspace(x_input(1,i),x_input(1,i+1));
    s1 = y_input(1,i).*((x_-x_input(1,i+1)).^2).*(delta_h(1,i) + 2.*(x_ - x_input(1,i)))./(delta_h(1,i).^3);
    s2 = y_input(1,i+1).*((x_-x_input(1,i)).^2).*(delta_h(1,i) + 2.*(x_input(1,i+1) - x_))./(delta_h(1,i).^3);
    s3 = m(1,i).*((x_ - x_input(1,i+1)).^2).*(x_ - x_input(1,i))./(delta_h(1,i).^2);
    s4 = m(1,i+1).*((x_ - x_input(1,i)).^2).*(x_ - x_input(1,i+1))./(delta_h(1,i).^2);
    s = s1 + s2 + s3 +s4;
    % 判断输入的x属于哪个插值区间，满足则计算对应的f(x)的值
    if x_value>x_input(1,i)&&x_value<x_input(1,i+1)
        s1_ = y_input(1,i)*((x_value-x_input(1,i+1)).^2)*(delta_h(1,i) + 2.*(x_value - x_input(1,i)))./(delta_h(1,i).^3);
        s2_ = y_input(1,i+1)*((x_value-x_input(1,i)).^2)*(delta_h(1,i) + 2.*(x_input(1,i+1) - x_value))./(delta_h(1,i).^3);
        s3_ = m(1,i)*((x_value - x_input(1,i+1)).^2)*(x_value - x_input(1,i))./(delta_h(1,i).^2);
        s4_ = m(1,i+1)*((x_value - x_input(1,i)).^2)*(x_value - x_input(1,i+1))./(delta_h(1,i).^2);
        s_ = s1_ + s2_ + s3_ +s4_;
        % 绘制计算点的结果
        plot(x_value,s_,'ro');
    end
    % 绘制每一段插值函数图像
    plot(x_,s,'b');
    hold on;
end
% 显示网格
plot(x_input,y_input,'*');
axis([0,10,0,6]);
grid on

验证数据：

已知数据

x	0	1	2	3	4	5	6	7	8	9	10
y	2.51	3.30	4.04	4.70	5.22	5.54	5.78	5.40	5.57	5.70	5.80

边界条件为y'0 = 0.8,y'n = 0.2.

测试结果：