前言

  本篇为插值法专栏第四篇内容讲述，此章主要讲述 Hermite（埃尔米特）插值法及matlab代码，其中也给出详细的例子让大家更好的理解Hermite插值法
提示 之前已经介绍牛顿插值法和三次样条插值，如果没看过前两篇的可以点击以下链接阅读

1. 数值分析（一）牛顿插值法
2. 数值分析（二）三次样条插值法
3. 数值分析（二续） 三次样条插值二类边界完整matlab代码

一、Hermite插值

  在许多实际应用中，不仅要求函数值相等，而且要求若干阶导数也相等。

1. Hermite定理

  定义：满足函数值相等且导数也相等的插值方法 $f(x)\approx \phi (x)$，φ(xi)=f(xi)(i=0,1,…,n)\\varphi ({x_i}) = f({x_i}) (i=0,1,\\ldots ,n)φ(xi​)=f(xi​)(i=0,1,…,n)φ′(xi)=f′(xi)\\varphi '({x_i}) = f'({x_i})φ′(xi​)=f′(xi​)φ(2)(xi)=f(2)(xi)\\varphi ^{(2)} ({x_i}) = f ^{(2)} ({x_i})φ(2)(xi​)=f(2)(xi​)$$⋮$$φ(m)(xi)=f(m)(xi)\\varphi ^{(m)} ({x_i}) = f ^{(m)} ({x_i})φ(m)(xi​)=f(m)(xi​)

定理：设 f(x)∈Cn[a,b],x0,...,xnf(x) \\in {C^n}[a,b],{x_0},...,{x_n}f(x)∈Cn[a,b],x0​,...,xn​ 为 $[a,b]$ 上的互异节点，则 f[x0,...,xn]f[{x_0},... ,{x_n}]f[x0​,...,xn​] 是其变量的连续函数

2. 重节点差商

  注意： 差商知识不清楚的话可以看之前这篇 数值分析（一）牛顿插值法 中有对差商详细的讲解
f[x0,x0]=lim⁡x1→x0f[x0,x1]=f(x1)−f(x0)x1−x0=f′(x0)f[{x_0},{x_0}] = \\mathop {\\lim }\\limits_{{x_1} \\to {x_0}} f[{x_0},{x_1}] = \\frac{{f({x_1}) - f({x_0})}}{{{x_1} - {x_0}}} = f'({x_0})f[x0​,x0​]=x1​→x0​lim​f[x0​,x1​]=x1​−x0​f(x1​)−f(x0​)​=f′(x0​)
f[x0,x0,x0]=lim⁡x1→x0x2→x0f[x0,x1,x2]=12!f′′(x0)f[{x_0},{x_0},{x_0}] = \\mathop{\\lim }\\limits_{{x_1} \\to {x_0}\\\\{{x_2} \\to {x_0}}}f[{x_0},{x_1},{x_2}] = \\frac{1}{{2!}}f''({x_0})f[x0​,x0​,x0​]=x1​→x0​x2​→x0​lim​f[x0​,x1​,x2​]=2!1​f′′(x0​)
一般 $n$ 阶重节点差商定义为：
f[x0,x0,⋯,x0]=lim⁡xi→x0f[x0,x1,⋯,xn]=1n!f(n)(x0)f[{x_0},{x_0},\\cdots,{x_0}] = \\mathop{\\lim }\\limits_{{x_i} \\to {x_0}}f[{x_0},{x_1},\\cdots,{x_n}] = \\frac{1}{{n!}}f^{(n)}({x_0})f[x0​,x0​,⋯,x0​]=xi​→x0​lim​f[x0​,x1​,⋯,xn​]=n!1​f(n)(x0​)

3. 重节点Newton插值

  在Newton插值公式中，令xi→x0,i=1,⋯,nx_i \\to x_0, i = 1, \\cdots, nxi​→x0​,i=1,⋯,n，则Nn(x)=f(x0)+f[x0,x1](x−x0)+⋯+f[x0,x1,⋯,xn]∏i=1n−1(x−xi)=f(x0)+f′(x0)(x−x0)+⋯+f(n)(x0)n!(x−x0)nN_n(x) = f(x_0)+f[x_0, x_1](x - x_0)+ \\cdots + f[x_0, x_1, \\cdots, x_n]\\prod_{i=1 }^{n-1} (x-x_i) \\\\ = f(x_0) + f'(x_0)(x-x_0)+ \\cdots + \\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n Nn​(x)=f(x0​)+f[x0​,x1​](x−x0​)+⋯+f[x0​,x1​,⋯,xn​]i=1∏n−1​(x−xi​)=f(x0​)+f′(x0​)(x−x0​)+⋯+n!f(n)(x0​)​(x−x0​)n

担心有些同学看不懂推导，解释一下如何转化为第二个等式（即，Taylor 插值多项式）：

将1.1中讲述的重节点差商 带入Newton插值公式中：lim⁡xi→x0f[x0,x1,⋯,xn]=1n!f(n)(x0)\\mathop{\\lim }\\limits_{{x_i} \\to {x_0}}f[{x_0},{x_1},\\cdots,{x_n}]=\\frac{1}{{n!}}f^{(n)}({x_0})xi​→x0​lim​f[x0​,x1​,⋯,xn​]=n!1​f(n)(x0​)所有的xi→x0x_i \\to x_0xi​→x0​，所以∏i=1n−1(x−xi)\\prod_{i=1 }^{n-1} (x-x_i)∏i=1n−1​(x−xi​) 就变成了(x−x0)n(x-x_0)^n(x−x0​)n

插值余项Rn(x)=f(n+1)(ξ)(n+1)!(x−x0)n+1R_n(x) = \\frac{f^{(n+1)}(\\xi)}{(n+1)!}(x-x_0)^{n+1}Rn​(x)=(n+1)!f(n+1)(ξ)​(x−x0​)n+1

4. Hermite 插值公式

  一般来说，给定 $m+1$ 个插值条件，就可以构造出一个 $m$ 次Hermite 插值多项式，接下来介绍两个典型的Hermite插值：三点三次 Hermite插值 和 两点三次 Hermite 插值

4.1 三点三次 Hermite插值

  插值节点：x0，x1，x2x_0，x_1， x_2x0​，x1​，x2​
  插值条件：P(xi)=f(xi)，i=0,1,2，P′(x1)=f′(x1)P(x_i) = f(x_i)，i=0,1,2，P'(x_1) = f'(x_1)P(xi​)=f(xi​)，i=0,1,2，P′(x1​)=f′(x1​)

  设：P(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+A(x−x0)(x−x1)(x−x2)P(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)P(x)=f(x0​)+f[x0​,x1​](x−x0​)+f[x0​,x1​,x2​](x−x0​)(x−x1​)+A(x−x0​)(x−x1​)(x−x2​)
将P′(x1)=f′(x1)P'(x_1) = f'(x_1)P′(x1​)=f′(x1​)代入可得A=f′(x1)−f[x0,x1]−f[x0,x1,x2](x1−x0)(x1−x0)(x1−x2)A = \\frac{f'(x_1)-f[x_0,x_1] - f[x_0,x_1,x_2](x_1-x_0)}{(x_1-x_0)(x_1-x_2)}A=(x1​−x0​)(x1​−x2​)f′(x1​)−f[x0​,x1​]−f[x0​,x1​,x2​](x1​−x0​)​
由于x0,x1,x2x_0,x_1,x_2x0​,x1​,x2​是$R(x)$的零点，且x1x_1x1​是二重零点，故可设，余项公式：R(x)=f(x)−P(x)=k(x)(x−x0)(x−x1)2(x−x2)R(x) = f(x) - P(x) = k(x)(x-x_0)(x-x_1)^2(x-x_2)R(x)=f(x)−P(x)=k(x)(x−x0​)(x−x1​)2(x−x2​)与 Lagrange 插值余项公式的推导过程类似，可得
R(x)=f(4)(ξx)4!(x−x0)(x−x1)2(x−x2)R(x) = \\frac{f^{(4)}(\\xi_x)}{4!}(x-x_0)(x-x_1)^2(x-x_2)R(x)=4!f(4)(ξx​)​(x−x0​)(x−x1​)2(x−x2​)
其中ξx\\xi_xξx​ 位于由 x0,x1,x2x_0, x_1, x_2x0​,x1​,x2​和$x$所界定的区间内

4.2 两点三次 Hermite插值

  插值节点：x0，x1x_0，x_1x0​，x1​
  插值条件：P(xi)=f(xi)=yi，P′(xi)=f′(xi)=mi，i=0,1P(x_i) = f(x_i)=y_i，P'(x_i) = f'(x_i)=m_i，i=0,1P(xi​)=f(xi​)=yi​，P′(xi​)=f′(xi​)=mi​，i=0,1

  模仿Lagrange 多项式的思想，设
H3(x)=y0α0(x)+y1α1(x)+m0β0(x)+m1β1(x)H_3(x) = y_0\\alpha_0(x)+y_1\\alpha_1(x)+m_0\\beta_0(x)+m_1\\beta_1(x)H3​(x)=y0​α0​(x)+y1​α1​(x)+m0​β0​(x)+m1​β1​(x)其中α0(x),α1(x),β0(x),β1(x)\\alpha_0(x), \\alpha_1(x), \\beta_0(x), \\beta_1(x)α0​(x),α1​(x),β0​(x),β1​(x)均为3次多项式，且满足αj(xi)=δij,αj′(xi)=0,βj(xi)=0,βj′(xi)=δji(i,j=0,1)\\alpha_j(x_i) = \\delta_{ij}, \\alpha'_j(x_i)=0,\\\\ \\beta_j(x_i)=0, \\beta'_j(x_i) = \\delta_{ji} \\\\ (i, j =0, 1)αj​(xi​)=δij​,αj′​(xi​)=0,βj​(xi​)=0,βj′​(xi​)=δji​(i,j=0,1)

上述中的δij\\delta_{ij}δij​表达式在Lagrange多项式中提到过，怕有些同学忘了这里帮忙回顾一下
lk(xi)=δki={1(i=k)0(i≠k)l_k(x_i) = \\delta_{ki} = \\left \\{\\begin{matrix}1&(i=k) \\\\ 0&(i \\ne k)\\end{matrix}\\right.lk​(xi​)=δki​={10​(i=k)(i=k)​

将插值条件代入立即可得 α0(x)\\alpha_0(x)α0​(x)求解过程

  Step 1：根据前面可以得到： α0(x0)=1,α0′(x0)=0,α0(x1)=0,α0′(x1)=0\\alpha_0(x_0) = 1, \\alpha'_0(x_0)=0, \\alpha_0(x_1) = 0, \\alpha'_0(x_1)=0α0​(x0​)=1,α0′​(x0​)=0,α0​(x1​)=0,α0′​(x1​)=0

  Step 2：由上面提到过α0(x),α1(x),β0(x),β1(x)\\alpha_0(x), \\alpha_1(x), \\beta_0(x), \\beta_1(x)α0​(x),α1​(x),β0​(x),β1​(x)都是3次多项式，α0(x1)=0,α0′(x1)=0\\alpha_0(x_1) = 0, \\alpha'_0(x_1)=0α0​(x1​)=0,α0′​(x1​)=0，根据以上条件可以构造函数为： α0(x)=(ax+b)(x−x1x0−x1)2\\alpha_0(x) = (ax + b)(\\frac{x-x_1}{x_0-x_1})^2α0​(x)=(ax+b)(x0​−x1​x−x1​​)2
  Step 3：将 α0(x0)=1,α0′(x0)=0\\alpha_0(x_0) = 1, \\alpha'_0(x_0) = 0α0​(x0​)=1,α0′​(x0​)=0 代入上式，可以求解得到 a=−2x0−x1,b=3x0−x1x0−x1=1+2x0x0−x1a = -\\frac{2}{x_0-x_1}, b = \\frac{3x_0-x_1}{x_0-x_1} = 1+ \\frac{2x_0}{x_0-x_1}a=−x0​−x1​2​,b=x0​−x1​3x0​−x1​​=1+x0​−x1​2x0​​
  Step 4：将求得出来的 $a,b$ 代入到构造的函数α0(x)\\alpha_0(x)α0​(x) 中，即得：α0(x)=(1+2x−x0x1−x0)(x−x1x0−x1)2\\alpha_0(x) = (1+2\\frac{x-x_0}{x_1-x_0})(\\frac{x-x_1}{x_0-x_1})^2α0​(x)=(1+2x1​−x0​x−x0​​)(x0​−x1​x−x1​​)2

注意：我在推公式的过程感觉是2前面是个减号(很容易出现自己推导和给出的不一致)，我仔细一看发现是这么一回事，那么我来把省略的一些步骤给他详细写一下
ax+b=−2xx0−x1+1+2x0x0−x1=1−2x−x0x0−x1=1+2x−x0x1−x0ax+b = -\\frac{2x}{x_0-x_1}+1+ \\frac{2x_0}{x_0-x_1} = 1 - 2\\frac{x-x_0}{x_0-x_1} = 1 + 2\\frac{x-x_0}{x_1-x_0}ax+b=−x0​−x1​2x​+1+x0​−x1​2x0​​=1−2x0​−x1​x−x0​​=1+2x1​−x0​x−x0​​才发现将-放到了分母

同理： α1(x)=(1+2x−x1x0−x1)(x−x0x1−x0)2\\alpha_1(x) = (1+2\\frac{x-x_1}{x_0-x_1})(\\frac{x-x_0}{x_1-x_0})^2α1​(x)=(1+2x0​−x1​x−x1​​)(x1​−x0​x−x0​​)2

类似得可以得到： β0(x)=(x−x0)(x−x1x0−x1)2\\beta_0(x) = (x-x_0)(\\frac{x-x_1}{x_0-x_1})^2β0​(x)=(x−x0​)(x0​−x1​x−x1​​)2 β1(x)=(x−x1)(x−x0x1−x0)2\\beta_1(x) = (x-x_1)(\\frac{x-x_0}{x_1-x_0})^2β1​(x)=(x−x1​)(x1​−x0​x−x0​​)2

满足插值条件 P(x0)=f(x0)=y0，P′(x0)=f′(x0)=m0P(x1)=f(x1)=y1，P′(x1)=f′(x1)=m1P(x_0) = f(x_0)=y_0，P'(x_0) = f'(x_0)=m_0 \\\\ P(x_1) = f(x_1)=y_1，P'(x_1) = f'(x_1)=m_1P(x0​)=f(x0​)=y0​，P′(x0​)=f′(x0​)=m0​P(x1​)=f(x1​)=y1​，P′(x1​)=f′(x1​)=m1​ 的三次Hermite插值多项式为
H3(x)=y0(1+2x−x0x1−x0)(x−x1x0−x1)2+y1(1+2x−x1x0−x1)(x−x0x1−x0)2+m0(x−x0)(x−x1x0−x1)2+m1(x−x1)(x−x0x1−x0)2H_3(x) = y_0(1+2\\frac{x-x_0}{x_1-x_0})(\\frac{x-x_1}{x_0-x_1})^2+y_1 (1+2\\frac{x-x_1}{x_0-x_1})(\\frac{x-x_0}{x_1-x_0})^2 \\\\ +m_0(x-x_0)(\\frac{x-x_1}{x_0-x_1})^2+m_1(x-x_1)(\\frac{x-x_0}{x_1-x_0})^2H3​(x)=y0​(1+2x1​−x0​x−x0​​)(x0​−x1​x−x1​​)2+y1​(1+2x0​−x1​x−x1​​)(x1​−x0​x−x0​​)2+m0​(x−x0​)(x0​−x1​x−x1​​)2+m1​(x−x1​)(x1​−x0​x−x0​​)2

余项R3(x)=f(4)(ξx)4!(x−x0)2(x−x1)2R_3(x) = \\frac{f^{(4)}(\\xi_x)}{4!}(x-x_0)^2(x-x_1)^2 R3​(x)=4!f(4)(ξx​)​(x−x0​)2(x−x1​)2其中ξx∈(x0,x1)\\xi_x \\in (x_0, x_1)ξx​∈(x0​,x1​)

4.3 $2n+1$次Hermite插值多项式

  由上面 4.2 中给出了 α0(x),α1(x),β0(x),β1(x)\\alpha_0(x), \\alpha_1(x), \\beta_0(x), \\beta_1(x)α0​(x),α1​(x),β0​(x),β1​(x) 的表达式，那么依次类推推导出 $2n+1$ 次Hermite插值多项式，且条件为H(xi)=f(xi)，H′(xi)=f′(xi)，(i=0,1,⋯,n)H(x_i) = f(x_i)，H'(x_i) = f'(x_i)，(i=0,1,\\cdots,n)H(xi​)=f(xi​)，H′(xi​)=f′(xi​)，(i=0,1,⋯,n)
  求解：αj(x)\\alpha_j(x)αj​(x) 和 βj(x)\\beta_j(x)βj​(x) $(j=0,1,\cdots ,n)$ αj(x)=[1−2(x−xj)∑k=0k≠jn1xj−xk]lj2(x)\\alpha_j(x) = [1-2(x-x_j)\\sum_{\\begin{matrix} k = 0\\\\ k\\ne j \\end{matrix}}^{n}\\frac{1}{x_j-x_k}]l_j^2(x)αj​(x)=[1−2(x−xj​)k=0k=j​∑n​xj​−xk​1​]lj2​(x)
βj(x)=(x−xj)lj2(x)\\beta_j(x) = (x-x_j)l_j^2(x)βj​(x)=(x−xj​)lj2​(x)
H2n+1(x)=∑j=0n[1−2(x−xj)∑k=0k≠jn1xj−xk]lj2(x)f(xj)+∑j=0n(x−xj)lj2(x)f′(xj)H_{2n+1}(x) = \\sum_{j=0}^{n}[1-2(x-x_j)\\sum_{\\begin{matrix} k = 0\\\\ k\\ne j \\end{matrix}}^{n}\\frac{1}{x_j-x_k}]l_j^2(x)f(x_j)+\\sum_{j=0}^{n}(x-x_j)l_j^2(x)f'(x_j)H2n+1​(x)=j=0∑n​[1−2(x−xj​)k=0k=j​∑n​xj​−xk​1​]lj2​(x)f(xj​)+j=0∑n​(x−xj​)lj2​(x)f′(xj​)

上述中的 lj(x)l_j(x)lj​(x) 在 Lagrange 多项式中提到过不太清楚得同学可以去回顾一下

二、Hermite插值算法及matlab代码

Matlab 版本号 2022b

1. $2n+1$次Hermite插值matlab代码实现

function y=Hermitezi(X,Y,Y1,x)
% 2n+1次Hermite插值函数
%   X为已知数据点的x坐标
%   Y为已知数据点的y坐标
%   Y1为数据点y值导数
%   x0为插值点的x坐标
if(length(X) == length(Y))if(length(X) == length(Y1))n=length(X);elsedisp('y和y的导数维数不相等');renturn;endelsedisp('x和y的维数不相等');return;
end
%以上为输入判断和确定“n”的值
m=length(x);
for t=1:mz=x(t);s=0.0;for i=1:nh=1.0;a=0.0;for j=1:nif(j~= i)h=h*(z-X(j))^2/((X(i)-X(j))^2);%求得值为(li(x))^2a=a+1/(X(i)-X(j));   %求得ai（x）表达式之中的累加部分endends=s+h*((X(i)-z)*(2*a*Y(i)-Y1(i))+Y(i));end
y(t)=s;
end
end