13. SciPy埃尔米特插值

埃尔米特插值是另一类插值问题，这类插值在给定的节点处，不但要求插值多项式的函数值与原函数值相同。同时还要求在节点处，插值多项式的一阶直至指定阶的导数值，也与被插函数的相应阶导数值相等，这样的插值称为埃尔米特($Hermite$)插值。

13.1 Hermite插值基本原理

$Hermite$插值在不同的节点，提出的差值条件个数可以不同，若在某节点$x_i$，要求插值函数多项式的函数值，一阶导数值，直至$\;m_i - 1\;$阶导数值均与被插函数的函数值相同及相应的导数值相等。我们称为重插值点节，因此，$Hermite$插值应给出两组数，一组为插值点{$\;x_i, \; i \in [0, n]\;$}节点，另一组为相应的重数标号{$m_i, \; i \in [0, n]$}。

对于$f(x)$在$[a, b]$上存在$n + 1$个互异点即$x_i$有$a \leq x_0 < x_1 < \cdots < x_n \leq b$,有$y_j = f(x_j)$, $m_j = f'(x_j)$对于$j$取值于$0,1,2,\cdots,n$，求插值多项式$H(x)$满足： $$H(x_j) = y_j,\; H'(x_j) = m_j, (j=0,1,2,\cdots,n)$$ 这里共需要$2n+2$个插值条件可以唯一确定一个$2n + 1$的多项式$H_{2n+1} = H(x)$,其形式为： $$ H_{2n + 1} = a_0 + a_1x+\cdots+a_{2n + 1}x^{2n + 1} $$ 先求出$\;2n+2\;$个插值基函数$\;\alpha_j(x)\;$和$\;\beta_j(x)\;$此处$\;j\;$取值$\;0,1,\cdots,n\;$，每个基函数都是$\;2n+1\;$次多项式，且满足： $$ \alpha_j(x_k)=\delta_{jk} = 0, \quad j\neq k\quad \alpha'_j(x_k) = 0 $$

$$ \alpha_j(x_k)=\delta_{jk} =1, \quad j = k\quad \alpha'_j(x_k) = 0 $$

$$ \beta_j(x_k) = 0, \quad \beta_j{'}(x_k) = \delta_{jk} $$

$$ j,k = 0,1,2,\cdots,n $$

此时，$H(x)$可以改写为： $$ H_{2n + 1}(x) = \sum_{j = 0}^{n}(y_j\alpha_j(x) + m_j\beta_j(x)) $$ 接下来求解这些插值基函数$\;\alpha_j(x)\;$和$\;\beta_j(x)\;$，利用朗日插值基函数$\;P_j(x)\;$，令： $$ \alpha_j(x) = (ax+b)P_j^2(x) $$ 有： $$ \alpha_j(x_j) = (ax_j+b)P_j^2(x_j) = 1 $$

$$ P_j(x_j) = 1 $$

对$\;\alpha_j(x)\;$求导得到$\;\alpha'_j(x)\;$即： $$ \alpha'_j(x) = aP_j^2(x_j) + 2(ax_j + b)P_j(x_j)P'_j(x_j) = 0 $$ 即： $$ a + 2(ax_j + b)P'_j(x_j) = 0 $$ 整理以上两式得到： $$ ax_j + b = 1 $$

$$ a + 2P'_j(x_j) = 0 $$

解之得： $$ a = -2P'_j(x_j) $$

$$ b = 1 + 2x_jP'_j(x_j) $$

已知 $$ P_j(x) = \prod_{i =0,i \neq j}^{n}\frac{x - x_i}{x_j - x_i} $$ 求导可得： $$ P_j^{'}(x_j) = \frac{1}{x_j - x_0}\times \frac{x_j - x_1}{x_j - x_1}\times \cdots\times \frac{x_j - x_n}{x_j - x_n} + \cdots + \frac{x_j - x_0}{x_j - x_0}\times \cdots \times \frac{x_j - x_{n-1}}{x_j - x_{n-1}}\times \frac{1}{x_j - x_n} =\sum_{i =0,i \neq j}^{n}\frac{1}{x_j - x_i} $$ 将$a = -2P_j^{'}(x_j)$和$b = 1 + 2x_jP_j^{'}(x_j)$代入公式$\;\alpha_j(x) = (ax+b)P_j^2(x)\;$得： $$ a_j(x) = (1-2(x - x_j)\sum_{i =0,i \neq j}^{n}\frac{1}{x_j - x_i})P_j^2(x) $$ 同理，$\beta_j(x)$为： $$ \beta_j(x) = (x - x_j)P_j^2(x) $$

13.2 Hermite插值计算

下面以经过$f(x) = x^3 + 1$的点(0,1)、(1,2)和导数$f'(0)=0$、$f'(1)=3$来求$Hermite$插值多项式。

这里有$x_0 =0, x_1 = 1, y_0 = 1, y_1 = 2$便可求得$\alpha_j(x)$和$\beta_j(x)$了。

首先计算$P_j(x)$: $$ P_0(x) = \frac{x - x_1}{x_0 - x_1} = \frac{x-1}{0-1} = 1-x $$

$$ P_1(x) = \frac{x - x_0}{x_1 - x_0} = \frac{x-0}{1-0} = x $$

接着计算$\alpha_j(x)$： $$ \alpha_0 = (1 - 2(x -0)\frac{1}{0-1})(1-x)^2 = (1+2x)(1-x)^2= 1-3x^2 + 2x^3 $$

$$ \alpha_1 = (1-2(x - 1)\frac{1}{1-0})x^2 = 3x^2-2x^3 $$

计算$\beta_j(x)$: $$ \beta_0(x) = (x - 0)(1-x)^2= x - 2x^2 + x^3 $$

$$ \beta_1 = (x - 1)x^2 = x^3 - x^2 $$

计算$H_{2n + 1}(x)$: $$ H_3(x) = \sum_{j = 0}^{1}(y_j\alpha_j(x) + m_j\beta_j(x)) = y_0\alpha_0(x) + m_0\beta_0(x)+y_1\alpha_1(x) + m_1\beta_1(x) $$ 代入上边求得的各值： $$ H_3(x) = 1\times(1-3x^2 + 2x^3) +0\times(x-2x^2+x^3) + 2\times(3x^2 -2x^3) + 3\times(x^3 - x^2) = x^3 + 1 $$ 即最后求得： $$ H_3(x) = x^3 + 1 $$ 与之前的$f(x)$表达式一模一样！

13.3 Hermite插值编程

SciPy里的Hermite插值实现函数是scipy.interpolate模块下的KroghInterpolator函数，其形参$x$和$y$必须同等长度，可以给出若干个坐标点作为形参，也可以给出每个$x_i$对应的$y_i$和导数。

13.3.1 坐标点作为形参

下面的程序代码是用四个坐标点作为函数的参数求得插值。KroghInterpolator函数要求$x_i$是递增序列。

from scipy.interpolate import KroghInterpolator
import numpy as np
import numpy as np, matplotlib.pyplot as plt

x = np.linspace(-5, 5, 100)
nodes = np.array([0,1,2,3])
def f(x):
    return x**3 + 1
xi = np.array([0,1,2,3])
yi = np.array([1,2,9,28])
interpolant = KroghInterpolator(xi, yi)
plt.figure()

plt.subplot(121)
plt.plot(x, interpolant(x), 'b--',label='Hermite Interpolation')
plt.plot(nodes, f(nodes), 'ro', label='nodes')
plt.legend(loc=9)
plt.xlim(-10.5, 10.5)
plt.title('$f(x) = x^3 + 1$')

plt.subplot(122)
plt.plot(x, f(x), 'g--', label='orginal')
plt.plot(xi, f(xi), 'ro', label='nodes')
plt.legend(loc=9)
plt.title('$f(x) = x^3 + 1$')
plt.xlim(-10.5, 10.5)
plt.show()

程序执行结果如下：

13.3.2 坐标和导数作为形参

由于KroghInterpolator要求$x$是递增序列，所以形参$x$里会多次出项某一个x轴坐标$x_i$，如果$x$里有多个同样的$x_i$,那么第二个参数$y$对应值分别为$x_i$的$f(x_i)$值、$f'(x_i)$一阶导数值、$f''(x_i)$二阶导数值$\cdots$$f^{(n)}(x_i)$n阶导数值。

from scipy.interpolate import KroghInterpolator
import numpy as np
import numpy as np, matplotlib.pyplot as plt

x = np.linspace(-5, 5, 100)
nodes = np.array([0,1])
def f(x):
    return x**3 + 1
xi = np.array([0,0,1,1])
yi = np.array([1,0,2,3])
interpolant = KroghInterpolator(xi, yi)
plt.figure()
plt.subplot(121)
plt.plot(x, interpolant(x), 'r--',label='Hermite Interpolation')
plt.plot(xi, yi, 'go', label='nodes', markersize=8)
plt.legend(loc=9)
plt.xlim(-10.5, 10.5)
plt.title('$f(x) = x^3 + 1$')

plt.subplot(122)
plt.plot(x, f(x), 'b--', label='orginal')
plt.plot(xi, f(xi), 'go', label='nodes', markersize=8)
plt.legend(loc=9)
plt.xlim(-10.5, 10.5)
plt.title('$f(x) = x^3 + 1$')
plt.show()

程序执行结果：

需要说明的函数的形参$xi$和$yi$:

xi = np.array([0,0,1,1])
yi = np.array([1,0,2,3])

在$xi$里有两个坐标$x_0= 0$和$x_1 = 1$,在$yi$里1是$f(0)$、0是$f'(0)$导数值、2是$f(1)$的计算值、3则是$f'(1)$的导数值。

感谢Klang(金浪)智能数据看板klang.org.cn鼎力支持！