数值分析——埃尔米特（Hermit）插值

创作时间:

作者:

@小白创作中心

数值分析——埃尔米特（Hermit）插值

引用

CSDN

https://blog.csdn.net/qq_40619699/article/details/140842711

埃尔米特（Hermit）插值是一种在插值点处不仅要求函数值相等，还要求导数值相等的插值方法。本文将详细介绍两种常见的埃尔米特插值形式：三点三次Hermit插值和两点三次Hermit插值，并提供MATLAB实现代码和测试用例。

理论推导

1.三点三次Hermit插值

考虑满足条件
$$
P(x_i) = H_3(x_i) = f(x_i) \quad (i=0,1,2)
$$
以及
$$
H_3'(x_k) = f'(x_k) \quad (k=0或1或2)
$$
因此可得三次Hermit插值函数有以下形式：
$$
H_3(x) = f(x_0) + fx_0, x_1 + fx_0, x_1, x_2(x-x_1) + A(x-x_0)(x-x_1)(x-x_2) \tag{1}
$$
其中A为待定系数，可通过条件
$$
H_3'(x_k) = f'(x_k) \quad (k=0或1或2)
$$
确定：
$$
A = \frac{f'(x_k) - f[x_0, x_1] - (x_1-x_0)f[x_0, x_1, x_2]}{(x_1-x_0)(x_1-x_2)} \tag{2}
$$

2.两点三次Hermit插值

考虑满足条件
$$
H_3(x_k) = y_k, \quad H_3(x_{k+1}) = y_{k+1}, \quad H_3'(x_k) = m_k, \quad H_3'(x_{k+1}) = m_{k+1}
$$
因此可得三次Hermit插值函数有以下形式：
$$
H_3(x) = \alpha_k(x)y_k + \alpha_{k+1}(x)y_{k+1} + \beta_k(x)m_k + \beta_{k+1}(x)m_{k+1} \tag{3}
$$
其中$\alpha_k(x), \alpha_{k+1}(x), \beta_k(x), \beta_{k+1}(x)$是插值点$x_k, x_{k+1}$的三次Hermit插值基函数：
$$
\begin{matrix}
\alpha_{k}(x_{k})=1 & \alpha_{k+1}(x_{k})=0 & \beta_{k}(x_{k})=0 & \beta_{k+1}(x_{k})=0 \
\alpha_{k}(x_{k+1})=0 & \alpha_{k+1}(x_{k+1})=1 & \beta_{k}(x_{k+1})=0 & \beta_{k+1}(x_{k+1})=0 \
\alpha'{k}(x{k})=0 & \alpha'{k+1}(x{k})=0 & \beta'{k}(x{k})=1 & \beta'{k+1}(x{k})=0 \
\alpha'{k}(x{k+1})=0 & \alpha'{k+1}(x{k+1})=0 & \beta'{k}(x{k+1})=0 & \beta'{k+1}(x{k+1})=1
\end{matrix} \tag{4}
$$
由$\alpha_k(x_{k+1}) = 0, \alpha'k(x{k+1}) = 0$可设：
$$
\alpha_k(x) = \left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 (a+bx) \tag{5}
$$
其中$a, b$是待定系数，可通过$\alpha_k(x_k) = 1, \alpha'k(x_k) = 0$求得：
$$
a = \frac{-2}{x_k-x{k+1}}, \quad b = 1+\frac{2x_k}{x_k-x_{k+1}} \tag{6}
$$
于是
$$
\alpha_k(x) = \left( 1+2\frac{x-x_k}{x_{k+1}-x_{k}} \right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 \tag{7}
$$
同理可得
$$
\begin{align}
\alpha_{k+1}(x) &= \left( 1+2\frac{x-x_{k+1}}{x_{k}-x_{k+1}} \right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^2 \
\beta_{k}(x) &= \left( x-x_k \right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 \
\beta_{k+1}(x) &= \left( x-x_{k+1} \right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^2
\end{align} \tag{8}
$$
因此两点三次Hermit插值函数：
$$
\begin{align}
H_3(x) = &\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 \left( 1+2\frac{x-x_k}{x_{k+1}-x_{k}} \right)f_k+ \left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^2\left( 1+2\frac{x-x_{k+1}}{x_{k}-x_{k+1}} \right)f_{k+1} \
&+\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2\left( x-x_k \right)f'k+\left(\frac{x-x{k}}{x_{k+1}-x_{k}}\right)^2\left( x-x_{k+1} \right)f'_{k+1}
\end{align} \tag{9}
$$

MATLAB实现

1.三点三次Hermit插值多项式

实现代码如下：

% 三点三次Hermit插值
function [aArr]=func_3dotsCubicHermitInterpolation(xArr,yArr,dy)
    % ************************************************************
    % 识别误差
    % ************************************************************
    xLength=length(xArr);
    yLength=length(yArr);
    % 插值数组不匹配
    if xLength~=yLength
        error('Error:The length of xArr=%d and yArr=%d are not equal in length.',xLength,yLength);
    end
    % ************************************************************
    % 计算插值多项式
    % ************************************************************
    syms x;
    A=(dy-func_DifferenceQuotient2(xArr(1:1:2),yArr(1:1:2),2)-...
        (xArr(2)-xArr(1))*func_DifferenceQuotient2(xArr(1:1:3),yArr(1:1:3),3))/...
        ((xArr(2)-xArr(1))*(xArr(2)-xArr(3)));
    aArr=sym2poly(poly2sym(func_NewtonInterpolation(xArr,yArr,xLength),x)+A*...
        (x-xArr(1))*(x-xArr(2))*(x-xArr(3)));
end

2.两点三次Hermit插值多项式

实现代码如下：

% 两点三次Hermit插值
function [aArr]=func_2dotsCubicHermitInterpolation(xArr,yArr,dy1,dy2)
    % ************************************************************
    % 识别误差
    % ************************************************************
    xLength=length(xArr);
    yLength=length(yArr);
    % 插值数组不匹配
    if xLength~=yLength
        error('Error:The length of xArr=%d and yArr=%d are not equal in length.',xLength,yLength);
    end
    % ************************************************************
    % 计算插值多项式
    % ************************************************************
    syms x;
    poly1=(1+2*(x-xArr(1))/(xArr(2)-xArr(1)))*((x-xArr(2))/(xArr(1)-xArr(2)))^2*yArr(1);
    poly2=(1+2*(x-xArr(2))/(xArr(1)-xArr(2)))*((x-xArr(1))/(xArr(2)-xArr(1)))^2*yArr(2);
    poly3=(x-xArr(1))*((x-xArr(2))/(xArr(1)-xArr(2)))^2*dy1;
    poly4=(x-xArr(2))*((x-xArr(1))/(xArr(2)-xArr(1)))^2*dy2;
    aArr=sym2poly(poly1+poly2+poly3+poly4);
end

3.三次Hermit插值多项式测试

针对函数$f(x) = x^{\frac{3}{2}}$的均差表如下：

基于上述均差表编写测试程序：

% 测试三次Hermit差分函数——func_3dotsCubicHermitInterpolation、func_2dotsCubicHermitInterpolation
clc;
clear;
close all;
x=[1/4,1,9/4];
y=[1/8,1,27/8];
A1=func_3dotsCubicHermitInterpolation(x,y,3/2)
A2=func_2dotsCubicHermitInterpolation(x(1:1:2),y(1:1:2),3/4,3/2)
polyval(A1,0.36)
polyval(A2,0.36)
0.36^(3/2)
polyval(A1,0.49)
polyval(A2,0.49)
0.49^(3/2)
polyval(A1,0.64)
polyval(A2,0.64)
0.64^(3/2)

测试结果如下：