数值分析——埃尔米特（Hermit）插值
关键字：Matalb；曲线拟合；导数条件
系列文章目录
数值分析——拉格朗日插值
数值分析——牛顿插值多项式
数值分析——分段低次插值
文章目录
- 系列文章目录
- 前言
- 一、理论推导
- 1.三点三次Hermit插值
- 2.两点三次Hermit插值
- 二、MATLAB实现
- 1.三点三次Hermit插值多项式
- 2.两点三次Hermit插值多项式
- 3.三次Hermit插值多项式测试
- 总结
- 参考文献
前言
在之前提到的插值方法——线性插值、抛物线插值、拉格朗日插值、牛顿插值中，插值函数需要满足的条件都是在插值点与原函数数值相同。在部分实际问题中要求在插值点的一阶导数甚至高阶导数也相等，满足这种条件的插值函数就是今天的主角——埃尔米特（Hermit）插值
本文对一些复杂的流程与概念进行了简化与省略，如果想详细了解的小伙伴可以参考李庆扬老师的《数值分析 -第5版》的第2章（需要电子版请私戳）
一、理论推导
1.三点三次Hermit插值
考虑满足条件 P ( x i ) = H 3 ( x i ) = f ( x i ) ( i = 0 , 1 , 2 ) P\left( x_i \right)=H_3\left( x_i \right)=f\left( x_i \right)\left( i=0,1,2 \right) P(xi)=H3(xi)=f(xi)(i=0,1,2)以及 H 3 ′ ( x k ) = f ′ ( x k ) ( k = 0 或 1 或 2 ) H^{'}_3\left( x_k \right)=f^{'}\left( x_k \right)\left( k=0或1或2 \right) H3′(xk)=f′(xk)(k=0或1或2)，因此可得三次Hermit插值函数有以下形式：
H 3 ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) + A ( x − x 0 ) ( x − x 1 ) ( x − x 2 ) (1) H_3\left( x \right)=f\left( x_0 \right)+f\left[ x_0,x_1 \right]\left(x-x_0\right)+f\left[ x_0,x_1,x_2 \right]\left(x-x_0\right)\left(x-x_1\right)+A\left(x-x_0\right)\left(x-x_1\right)\left(x-x_2\right) \tag{1} H3(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+A(x−x0)(x−x1)(x−x2)(1)其中 A A A为待定系数，可通过条件 H 3 ′ ( x k ) = f ′ ( x k ) ( k = 0 或 1 或 2 ) H^{'}_3\left( x_k \right)=f^{'}\left( x_k \right)\left( k=0或1或2 \right) H3′(xk)=f′(xk)(k=0或1或2)确定：
A = 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 ) (2) A=\frac{f^{'}\left( x_k \right)-f\left[ x_0,x_1 \right]-\left(x_1-x_0\right)f\left[ x_0,x_1,x_2 \right]}{\left(x_1-x_0\right)\left(x_1-x_2\right)} \tag{2} A=(x1−x0)(x1−x2)f′(xk)−f[x0,x1]−(x1−x0)f[x0,x1,x2](2)
2.两点三次Hermit插值
考虑满足条件 H 3 ( x k ) = y k 、 H 3 ( x k + 1 ) = y k + 1 、 H 3 ′ ( x k ) = m k 、 H 3 ′ ( x k + 1 ) = m k + 1 H_3\left( x_k \right)=y_k、H_3\left( x_{k+1} \right)=y_{k+1}、H^{'}_3\left( x_k \right)=m_k、H^{'}_3\left( x_{k+1} \right)=m_{k+1} H3(xk)=yk、H3(xk+1)=yk+1、H3′(xk)=mk、H3′(xk+1)=mk+1，因此可得三次Hermit插值函数有以下形式：
H 3 ( x ) = α k ( x ) y k + α k + 1 ( x ) y k + 1 + β k ( x ) m k + β k + 1 ( x ) m k + 1 (3) H_3\left( x \right)=\alpha_k\left( x \right)y_k+\alpha_{k+1}\left( x \right)y_{k+1}+\beta_k\left( x \right)m_k+\beta_{k+1}\left( x \right)m_{k+1} \tag{3} H3(x)=αk(x)yk+αk+1(x)yk+1+βk(x)mk+βk+1(x)mk+1(3)其中 α k ( x ) 、 α k + 1 ( x ) 、 β k ( x ) 、 β k + 1 ( x ) \alpha_k\left( x \right)、\alpha_{k+1}\left( x \right)、\beta_k\left( x \right)、\beta_{k+1}\left( x \right) αk(x)、αk+1(x)、βk(x)、βk+1(x)是插值点 x k 、 x k + 1 x_{k}、x_{k+1} xk、xk+1的三次Hermit插值基函数：
α k ( x k ) = 1 α k + 1 ( x k ) = 0 β k ( x k ) = 0 β k + 1 ( x k ) = 0 α k ( x k + 1 ) = 0 α k + 1 ( x k + 1 ) = 1 β k ( x k + 1 ) = 0 β k + 1 ( x k + 1 ) = 0 α k ′ ( x k ) = 0 α k + 1 ′ ( x k ) = 0 β k ′ ( x k ) = 1 β k + 1 ′ ( x k ) = 0 α k ′ ( x k + 1 ) = 0 α k + 1 ′ ( x k + 1 ) = 0 β k ′ ( x k + 1 ) = 0 β k + 1 ′ ( x k + 1 ) = 1 (4) \begin{matrix} \alpha_{k}\left( x_{k} \right)=1 & \alpha_{k+1}\left( x_{k} \right)=0 & \beta_{k}\left( x_{k} \right)=0 & \beta_{k+1}\left( x_{k} \right)=0 \\ \alpha_{k}\left( x_{k+1} \right)=0 & \alpha_{k+1}\left( x_{k+1} \right)=1 & \beta_{k}\left( x_{k+1} \right)=0 & \beta_{k+1}\left( x_{k+1} \right)=0 \\ \alpha^{'}_{k}\left( x_{k} \right)=0 & \alpha^{'}_{k+1}\left( x_{k} \right)=0 & \beta^{'}_{k}\left( x_{k} \right)=1 & \beta^{'}_{k+1}\left( x_{k} \right)=0 \\ \alpha^{'}_{k}\left( x_{k+1} \right)=0 & \alpha^{'}_{k+1}\left( x_{k+1} \right)=0 & \beta^{'}_{k}\left( x_{k+1} \right)=0 & \beta^{'}_{k+1}\left( x_{k+1} \right)=1 \end{matrix} \tag{4} αk(xk)=1αk(xk+1)=0αk′(xk)=0αk′(xk+1)=0αk+1(xk)=0αk+1(xk+1)=1αk+1′(xk)=0αk+1′(xk+1)=0βk(xk)=0βk(xk+1)=0βk′(xk)=1βk′(xk+1)=0βk+1(xk)=0βk+1(xk+1)=0βk+1′(xk)=0βk+1′(xk+1)=1(4)由 α k ( x k + 1 ) = 0 、 α k ′ ( x k + 1 ) = 0 \alpha_{k}\left( x_{k+1} \right)=0、\alpha^{'}_{k}\left( x_{k+1} \right)=0 αk(xk+1)=0、αk′(xk+1)=0可设：
α k ( x ) = ( x − x k + 1 x k − x k + 1 ) 2 ( a + b x ) (5) \alpha_k\left( x \right)=\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 \left( a+bx \right) \tag{5} αk(x)=(xk−xk+1x−xk+1)2(a+bx)(5)其中 a 、 b a、b a、b是待定系数，可通过 α k ( x k ) = 1 、 α k ′ ( x k ) = 0 \alpha_{k}\left( x_{k} \right)=1、\alpha^{'}_{k}\left( x_{k} \right)=0 αk(xk)=1、αk′(xk)=0求得：
a = − 2 x k − x k + 1 , b = 1 + 2 x k x k − x k + 1 (6) a=\frac{-2}{x_k-x_{k+1}},b=1+\frac{2x_k}{x_k-x_{k+1}} \tag{6} a=xk−xk+1−2,b=1+xk−xk+12xk(6)于是
α k ( x ) = ( 1 + 2 x − x k x k + 1 − x k ) ( x − x k + 1 x k − x k + 1 ) 2 (7) \alpha_k\left( x \right)=\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} αk(x)=(1+2xk+1−xkx−xk)(xk−xk+1x−xk+1)2(7)同理可得
α k + 1 ( x ) = ( 1 + 2 x − x k + 1 x k − x k + 1 ) ( x − x k x k + 1 − x k ) 2 β k ( x ) = ( x − x k ) ( x − x k + 1 x k − x k + 1 ) 2 β k + 1 ( x ) = ( x − x k + 1 ) ( x − x k x k + 1 − x k ) 2 (8) \begin{align} \alpha_{k+1}\left( x \right)&=\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}\left( x \right)&=\left( x-x_k \right)\left(\frac{x-x_{k+1}}{x_{k}-x_{k+1}}\right)^2 \\ \beta_{k+1}\left( x \right)&=\left( x-x_{k+1} \right)\left(\frac{x-x_{k}}{x_{k+1}-x_{k}}\right)^2 \end{align} \tag{8} αk+1(x)βk(x)βk+1(x)=(1+2xk−xk+1x−xk+1)(xk+1−xkx−xk)2=(x−xk)(xk−xk+1x−xk+1)2=(x−xk+1)(xk+1−xkx−xk)2(8)因此两点三次Hermit插值函数：
H 3 ( x ) = ( x − x k + 1 x k − x k + 1 ) 2 ( 1 + 2 x − x k x k + 1 − x k ) f k + ( x − x k x k + 1 − x k ) 2 ( 1 + 2 x − x k + 1 x k − x k + 1 ) f k + 1 + ( x − x k + 1 x k − x k + 1 ) 2 ( x − x k ) f k ′ + ( x − x k x k + 1 − x k ) 2 ( x − x k + 1 ) f k + 1 ′ (9) \begin{align} H_3\left( x \right)= &\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} H3(x)=(xk−xk+1x−xk+1)2(1+2xk+1−xkx−xk)fk+(xk+1−xkx−xk)2(1+2xk−xk+1x−xk+1)fk+1+(xk−xk+1x−xk+1)2(x−xk)fk′+(xk+1−xkx−xk)2(x−xk+1)fk+1′(9)
二、MATLAB实现
1.三点三次Hermit插值多项式
实现代码如下：
% 三点三次Hermit插值
function [aArr]=func_3dotsCubicHermitInterpolation(xArr,yArr,dy)% ************************************************************% 识别误差% ************************************************************xLength=length(xArr);yLength=length(yArr);% 插值数组不匹配if xLength~=yLengtherror('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~=yLengtherror('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 3 2 f\left( x \right)=x^{\frac{3}{2}} f(x)=x23的均差表如下：
基于上述均差表编写测试程序：
% 测试三次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)
测试结果如下：
总结
以上在本文中对Hermit插值的两种常见形式的推导进行了简单的介绍，并提供了MATLAB的实现代码与测试用例。
参考文献
[1]李庆扬,王能超,易大义. 数值分析.第5版[M]. 清华大学出版社, 2008.