关键字：Matalb；曲线拟合；导数条件
系列文章目录
数值分析——拉格朗日插值
数值分析——牛顿插值多项式
数值分析——分段低次插值
数值分析——三次样条插值
文章目录
前言
在之前提到的插值方法——线性插值、抛物线插值、拉格朗日插值、牛顿插值中，插值函数需要满足的条件都是在插值点与原函数数值相同。在部分实际问题中要求在插值点的一阶导数甚至高阶导数也相等，满足这种条件的插值函数就是今天的主角——埃尔米特（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~=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
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.
6767
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
