MTALAB实现三次样条插值及埃尔米特（Hermite）插值

问题描述：

本文将采用三次样条插值及赫米特（Hermite）插值解决这一问题。

Sec 1

三次样条插值的数学原理可自行百度

function Sanci(N)
%程序使用三次样条插值
%N代表节点个数，数组形式

n=length(N);
for i=1:n
    jiedian=linspace(-2*pi,2*pi,N(i));%插值节点的确认
    %三次样条插值
    %计算函数值
    for j=1:N(i)
        F(j)=0.5*jiedian(j)-cos(jiedian(j));
    end
    
    %计算差商
    deltx=diff(jiedian);
    delty=diff(F);
    firstorder=delty./deltx ;%一阶差商
    for k=1:length(jiedian)-2
        delt2x(k)=jiedian(k+2)-jiedian(k);
    end
    delt2y=diff(firstorder);
    D=delt2y./delt2x; %二阶差商
    
    %构造系数矩阵B
    CN=length(jiedian);
    for k=1:CN
        for j=1:CN
            if j==k
                B(k,j)=2;
                if k==1
                    B(k,k+1)=1;
                end
                if k==CN
                    B(k,k-1)=1;
                end
                if k>=2
                    B(k,j-1)=0.5;
                    B(k,j+1)=0.5;
                end
                if j~=k&&j~=k-1&&j~=k+1
                    B(k,j)=0;
                end
            end
        end
    end
    %求解待定参数M
    %根据手算得出边界条件对应D值得表达式D=(6*(f_i-f_i-1)-3*h_i)/h_i^2
    D_1=(6*firstorder(1)-12*pi/N(i))/(4*pi/N(i))^2;
    D_Last=(6*firstorder(CN-1)-12*pi/N(i))/(4*pi/N(i))^2;
    C=[D_1,D,D_Last];%附加边界条件的新系数矩阵
    M=(B'*C')';          
    for j=1:N(i)-1
        x1=jiedian(j);x2=jiedian(j+1);
        x_mid(j)=(x1+x2)/2;
        %构造插值函数
        S1(j)=(((x2-x_mid(j))^3)*M(j)+((x_mid(j)-x1)^3)*M(j+1))/(24*pi/N(i))+((x2-x_mid(j))*F(j)+(x_mid(j)-x1)*F(j+1))/(4*pi/N(i))-(4*pi*((x2-x_mid(j))*M(j)+(x_mid(j)-x1)*M(j+1)))/(6*N(i));
    end
    save('三次插值中点値.mat','S1');
    save('z','x_mid');
    subplot(3,2,i);
    plot(x_mid,S1,'k-')%绘图
    hold on;
    plot(x_mid,S1,'ko');%三次样条插值点
    hold on;
    %绘制f=0.5*x-cosx图形
    sym x;
    f='0.5*x-cos(x)';
    h=ezplot(f,[-2*pi,2*pi]);
    set(h,'color','r');
    bt=strcat('节点=',num2str(N(i)));% 字符串连接
    title(bt);
end

在命令行窗口输入：

Sanci([7,14,21,56,112])

可得到下面的图片，可以看到，插值节点越多，插值曲线与原曲线拟合程度越高

Sec 2

赫米特插值

function Hermite(N)
%程序使用赫米特插值
%N代表节点个数，数组形式

n=length(N);
for i=1:n
    jiedian=linspace(-2*pi,2*pi,N(i));%插值节点的确认
    %Hermite 插值绘图
    for j=1:N(i)-1
        x1=jiedian(j);x2=jiedian(j+1);
        x_mid(j)=(x1+x2)/2;
        h1(j)=0.5*x1-cos(x1); %插值节点x1的函数值
        h2(j)=0.5*x2-cos(x2); %插值节点x2的函数值
        hx1(j)=0.5+sin(x1); %插值节点x1的导数值
        hx2(j)=0.5+sin(x2); %插值节点x2的导数值
        %插值函数值的计算
        H(j)=(1+2*(x_mid(j)-x1)/(x2-x1))*(((x2-x_mid(j))/(x2-x1))^2)*h1(j)+(1+2*(x2-x_mid(j))/(x2-x1))*(((x_mid(j)-x1)/(x2-x1))^2)*h2(j)+((x_mid(j)-x1)*((x2-x_mid(j))^2)/((x2-x1)^2))*hx1(j)-((x2-x_mid(j))*((x_mid(j)-x1)^2)/((x2-x1)^2))*hx2(j);
    end
    save('Hermite插值中点値.mat','H');
    subplot(3,2,i);
plot(x_mid,H,'b-');%绘图
hold on;
    plot(x_mid,H,'b*');%赫米特插值点
    %绘制f=0.5*x-cosx图形
    sym x;
    f='0.5*x-cos(x)';
    h=ezplot(f,[-2*pi,2*pi]);
    set(h,'color','r');
    grid on;
    bt=strcat('节点=',num2str(N(i)));% 字符串连接
    title(bt);
end

在命令行窗口输入：

>>Hermite([7,14,28,56,112])