分段三次Hermite插值的原理、matlab代码和Python代码网站首页 技术杂谈

分段三次Hermite插值的原理、matlab代码和Python代码

天下弈星~ 2023-06-08 12:00:02

简介python埃尔米特插值

一、分段三次Hermite插值

已知 $y_{i}=f(x_{i}),m_{i}=f^{'}(x_{i})(i=0,1,...,n) ,a=x_{0}<x_{1}<...<x_{n}=b$ ，求一个分段函数H(x)，使其满足：
（1） $H(x_{i})=y_{i},H^{'}(x_{i})=m_{i}(i=0,1,...,n)$

（2）在每个子区间 $[x_{i},x_{i+1}]$ 上，H(x)是次数不超过3的多项式。

称满足上述条件的函数H(x)为分段三次Hermite插值多项式。

二、分段三次Hermite插值多项式的推导

采用基函数的方法来构造 $H_{3}(x)$ 。

将 $H_{3}(x)$ 表示为：

$H_{3}(x)=y_{0}alpha _{0}(x)+y_{1}alpha _{1}(x)+m_{0}eta _{0}(x)+m_{1}eta _{1}(x)$

其中 $alpha _{0}(x),alpha _{1}(x),eta _{0}(x),eta _{1}(x)$ 为插值函数，且均为次数不超过3的多项式。为满足插值条件，它们应满足：

$alpha _{i}(x_{j})=delta _{ij}=left{egin{matrix} 0,i eq j & \ 1,i=j & end{matrix} ight.$ , $alpha_{i}^{'}(x_{j})=0$

$eta _{i}(x_{j})=0,eta _{i}^{'}=delta _{ij}$

$(i,j=0,1)$

$H_{3}(x_{i})=y_{i},H_{3}^{'}(x_{i})=m_{i},(i=0,1)$

$H_{3}(x)=y_{0}alpha _{0}(x)+y_{1}alpha _{1}(x)+m_{0}eta _{0}(x)+m_{1}eta _{1}(x)$

由于 $alpha _{0}(x_{1})=alpha _{0}^{'}(x_{1})=0$ ,故 $alpha _{0}(x)$ 含有 $(x-x_{1})^{2}$ 因子。可设

$alpha _{0}(x)=(a+b(x-x_{0}))(x-x_{1})^{2}$

其中a,b为待定系数。

由 $alpha _{0}(x_{0})=1$ ，可得 $a=frac{1}{(x_{0}-x_{1})^{2}}$ 。

由 $alpha _{0}^{'}(x_{0})=0$ ，可得 $b=frac{-2}{(x_{0}-x_{1})^{3}}$ 。

将a,b代入得

$a_{0}(x)=(1+2frac{x-x_{0}}{x_{1}-x_{0}})(frac{x-x_{1}}{x_{0}-x_{1}})^{2}$

类似地，将 $x_{0},x_{1}$ 互换，可得：

$a_{1}(x)=(1+2frac{x-x_{1}}{x_{0}-x_{1}})(frac{x-x_{0}}{x_{1}-x_{0}})^{2}$

由于 $eta _{0}(x_{0})=eta _{0}(x_{1})=eta _{0}^{'}(x_{1})=0$ ，故 $eta _{0}(x)$ 含有 $(x-x_{0})(x-x_{1})^{2}$ 因子。可设

$eta _{0}=c(x-x_{0})(x-x_{1})^{2}$

其中c为待定系数。

由 $eta _{0}^{'}(x_{0})=1$ ，可得 $c=frac{1}{(x_{0}-x_{1})^{2}}$ 。

$eta _{0}(x)=(x-x_{0})(frac{x-x_{1}}{x_{0}-x_{1}})^{2}$

类似地，将 $x_{0},x_{1}$ 互换，可得：

$eta _{1}(x)=(x-x_{1})(frac{x-x_{0}}{x_{1}-x_{0}})^{2}$

综上所述，三次Hermite插值多项式 $H_{3}(x)$ 的表达式为：

$H_{3}(x)=y_{0}alpha _{0}(x)+y_{1}alpha _{1}(x)+m_{0}eta _{0}(x)+m_{1}eta _{1}(x)$

$a_{0}(x)=(1+2frac{x-x_{0}}{x_{1}-x_{0}})(frac{x-x_{1}}{x_{0}-x_{1}})^{2}$

$a_{1}(x)=(1+2frac{x-x_{1}}{x_{0}-x_{1}})(frac{x-x_{0}}{x_{1}-x_{0}})^{2}$

$eta _{0}(x)=(x-x_{0})(frac{x-x_{1}}{x_{0}-x_{1}})^{2}$

$eta _{1}(x)=(x-x_{1})(frac{x-x_{0}}{x_{1}-x_{0}})^{2}$

可以证明，其余项为：

$R_{3}(x)=f(x)-H_{3}(x)=frac{f^{(4)}(xi )}{4!}(x-x_{0})^{2}(x-x_{1})^{2}$

其中， $xi$ 介于 $x_{0},x_{1}$ 之间。

注：

作为多项式插值，三次已经是较高的次数，次数再高就有可能发生Runge现象，因此，对有n+1节点的插值问题，我们可以使用分段两点三次Hermite插值。

三、分段三次Hermite插值的matlab实现

例：已知f(x)在节点1，2处的函数值为f(1)=2,f(2)=3,f(x)在节点1，2处的导数值为 $f^{'}(1)=0,f^{'}(2)=-1$ 。求f(x)的两点三次插值多项式及f(x)在x=1.5,1.7处的函数值。

function Hermite(A,B,C)
%输入  A代表分段三次Hermite插值所对应的节点
%      B代表节点对应的函数值
%      C代表节点对应的一阶导数值
%输出  f代表分段三次Hermite插值多项式
syms x a1 a2 b1 b2 df f
a1=(1+2*(x-A(1))/(A(2)-A(1)))*((x-A(2))/(A(1)-A(2)))^2
a2=(1+2*(x-A(2))/(A(1)-A(2)))*((x-A(1))/(A(2)-A(1)))^2
b1=(x-A(1))*((x-A(2))/(A(1)-A(2)))^2
b2=(x-A(2))*((x-A(1))/(A(2)-A(1)))^2
df=B(1)*a1+B(2)*a2+C(1)*b1+C(2)*b2
f=collect(df,x)

在命令行窗口中输入：
>> A=[1,2];
>> B=[2,3];
>> C=[0,-1];
>> Hermite(A,B,C)

最后得到的结果如下：

a1 =
 
(2*x - 1)*(x - 2)^2
 
 
a2 =
 
-(2*x - 5)*(x - 1)^2
 
 
b1 =
 
(x - 1)*(x - 2)^2
 
 
b2 =
 
(x - 1)^2*(x - 2)
 
 
df =
 
2*(2*x - 1)*(x - 2)^2 - (x - 1)^2*(x - 2) - 3*(2*x - 5)*(x - 1)^2
 
 
f =
 
- 3*x^3 + 13*x^2 - 17*x + 9

因此，f(x)的两点三次插值多项式是： $-3x^{3}+13x^{2}-17x+9$

四、分段三次Hermite插值的Python实现

同样采用上面的例子。

from sympy import *
def Hermite(A,B,C):
    x=Symbol('x')
    a1=Symbol('a1')
    a2=Symbol('a2')
    b1=Symbol('b1')
    b2=Symbol('b2')
    df=Symbol('df')
    f=Symbol('f')
    a1=(1+2*(x-A[0])/(A[1]-A[0]))*((x-A[1])*(A[0]-A[1]))**2
    a2=(1+2*(x-A[1])/(A[0]-A[1]))*((x-A[0])/(A[1]-A[0]))**2
    b1=(x-A[0])*((x-A[0])/(A[1]-A[0]))**2
    b2=(x-A[1])*((x-A[0])*(A[1]-A[0]))**2
    df=B[0]*a1+B[1]*a2+C[0]*b1+C[1]*b2
    f=factor(df)
    return(f)

A=[1,2]
B=[2,3]
C=[0,-1]
f=Hermite(A,B,C)
print(f)

结果如下：