分段三次多项式插值

分段三次多项式插值是一种应用最广泛的多项式插值方法，背景是：给定$n+1$个节点及函数值$(x_i, f_i), i = 0,1,2,\cdots,n$，要求过这些点的一个分段三次多项式$p(x)$，它满足在这些点上具有二阶可微的性质。

分析

首先我们来分析这样的多项式$p(x)$的存在唯一性，分析若不加任何限制，函数$p(x)$被分为了$n$个三次多项式，因此它的自由度为$4n$，在$n$个节点处，需要满足二阶可微的条件，因此：节点$(x_i, f_i), i = 0,1,2,\cdots,n$函数值相等、内部节点$(x_i, f_i), i = 1,2,\cdots,n-1$一二阶左右导数值相等，得到$3n-2$个线性无关的限制条件，还需要另外两个条件就可以唯一确定多项式$p(x).$

为了给出这两个条件，我们可以追加边界条件：

1. （第一类）$p'(x_0) = f'_0, p'(x_n) = f'_n;$
2. （第二类）$p''(x_0) = f''_0, p''(x_n) = f''_n.$

以上等式右端均为常数。

三弯矩方法

以第一类边界条件为例，在证明$p(x)$存在唯一性时也将解构造了出来。引入记号$I_i = [x_{i-1}, x_i], h_i = x_i - x_{i-1}, i = 1,2,\cdots,n.$

首先，函数$p''(x)$在每个区间$I_i$上是一次函数，那么有

$$p_i''(x) = \dfrac{M_{i-1}}{h_i}(x_i - x) + \dfrac{M_i}{h_i} (x - x_{i-1}) , \quad x \in I_i.$$

其中$M_i = p''(x_i), i = 0,1,2,\cdots,n.$是待定常数。

积分两次，注意到$p(x_i) = f_i, i = 0,1,2,\cdots,n$，因此

$${\begin{aligned}p_{i}(x)&={\dfrac {M_{i-1}}{6h_{i}}}(x_{i}-x)^{3}+{\dfrac {M_{i}}{6h_{i}}}(x-x_{i-1})^{3}\\&\quad +{\dfrac {f_{i-1}-M_{i-1}}{6}}(x_{i}-x)+{\dfrac {f_{i}-M_{i}}{6}}(x-x_{i-1}),\quad \forall x\in I_{i}.\end{aligned}}$$

注意到$p'(x)$在$\bigcup_{i=1}^n I_i$上连续，那么

$$p'_i(x_i^-) = p'_{i+1} (x_i^+), i = 1,2,\cdots,n-1.$$

对式#A1求一次导数，代入对应的$x_i$得到

$$\dfrac{h_i}{6}M_{i-1} + \dfrac{h_i + h_{i+1}}{3} M_i + \dfrac{h_{i+1}}{6}M_{i+1} = \dfrac{f_{i+1} - f_i}{h_{i+1}} - \dfrac{f_i - f_{i-1}}{h_i}, i = 1,2,\cdots,n-1.$$

引入记号

$$\lambda_i = \dfrac{h_{i+1}}{h_i + h_{i+1}}, \quad \mu_i = \dfrac{h_i}{h_i + h_{i+1}}.$$

即有递推公式

$$\mu _{i}M_{i-1}+2M_{i}+\lambda _{i}M_{i+1}=:d_{i}.$$

给出的边界条件是第一类的，带入边界条件可获得初始迭代值以及终止迭代值满足

$$\begin{cases} -\dfrac{h_1}{3} M_0 - \dfrac{h_1}{6} M_1 + \dfrac{f_1 - f_0}{h_1} &= f'_0, \\ \dfrac{h_n}{3} M_{n-1} + \dfrac{h_n}{6} M_n + \dfrac{f_n - f_{n-1}}{h_n} &= f'_n. \end{cases}$$

化简得到

$$\begin{cases} 2M_0 + M_1 &= \dfrac{6}{h_1} \left( \dfrac{f_1 - f_0}{h_1} - f'_0 \right), \\ M_{n-1} + 2M_0 &= \dfrac{6}{h_n} \left( f'_n \dfrac{f_n - f_{n-1}}{h_n} \right). \end{cases}$$

将上述条件以及递推公式写成矩阵方程

$$\begin{pmatrix} 2 & 1 \\ \mu_1 & 2 & \lambda_1 \\ &\mu_2 & 2 & \lambda_2 \\ && \ddots & \ddots & \ddots \\ &&&\mu_{n-1} & 2 & \lambda_{n-1} \\ &&&& 1 & 2 \end{pmatrix} \begin{pmatrix} M_0 \\ M_1 \\ M_2 \\ \vdots \\ M_{n-1} \\ M_n \end{pmatrix} = \begin{pmatrix} d_0 \\ d_1 \\ d_2 \\ \vdots \\ d_{n-1} \\ d_n \end{pmatrix}$$

左端的行列式不为零，方程有唯一解，解出$M_i$代入$p(x)$后就得到所求多项式。

三转角方法

记$m_{i} = p^{\prime}\left(x_{i}\right), i=0,1,2,\cdots,n$，将其设为待定参数。利用两点 Hermite 插值公式可得$p(x)$的分段表达式

$$\begin{align} p(x) =& m_{i-1} \frac{\left(x_{i}-x\right)^{2}\left(x-x_{i-1}\right)}{h_{i}^{2}}-m_{i} \frac{\left(x-x_{i-1}\right)^{2}\left(x_{i}-x\right)}{h_{i}^{2}} \\ & \quad + f_{i-1} \frac{\left(x_{i}-x\right)^{2}\left[2\left(x-x_{i-1}\right)+h_{i}\right]}{h_{i}^{3}} + f_{i} \frac{\left(x-x_{i-1}\right)^{2}\left[2\left(x_{i}-x\right)+h_{i}\right]}{h_{i}^{3}}, \quad x \in I_i. \end{align}$$

对$p(x)$求二阶导数, 有

$$\begin{align} p''(x) &= -2 m_{i-1} \frac{2 x_{i}+x_{i-1}-3 x}{h_{i}^{2}}-2 m_{i} \frac{2 x_{i-1}+x_{i}-3 x}{h_{i}^{2}} \\ &\quad + 6 \frac{f_{i}-f_{i-1}}{h_{i}^{3}} (x_{i}+x_{i-1}-2 x). \end{align}$$

显然，$p(x)$和$p'(x)$在区间$\bigcup_{i=1}^n I_i$上连续，得

$$p''(x_i^+) = p''(x_i^-), \quad i = 1,2,\cdots,n-1.$$

因此有

$$\begin{cases} p''(x_i^-) = \dfrac{2 m_{i-1}}{h_{i}}+\dfrac{4 m_{i}}{h_{i}}-6 \dfrac{f_{i}-f_{i-1}}{h_{i}^{2}}, \\ p''(x_i^+) = -\dfrac{4 m_{i-1}}{h_{i}}-\dfrac{2 m_{i}}{h_{i}}+6 \dfrac{f_{i}-f_{i-1}}{h_{i}^{2}}. \end{cases}$$

利用上式得到

$$\lambda_i m_{i-1} + 2 m_i + \mu_i m_{i+1} =: c_i, \quad i = 1,2,\cdots,n-1.$$

其中$\lambda_i, \mu_i$假设同上一小节。利用附加边界条件有

$$m_0 = f_0', \quad m_n = f_n'.$$

将上述条件以及递推公式写成矩阵方程

$$\begin{pmatrix} 1 \\ \mu_1 & 2 & \lambda_1 \\ &\mu_2 & 2 & \lambda_2 \\ && \ddots & \ddots & \ddots \\ &&&\mu_{n-1} & 2 & \lambda_{n-1} \\ &&&&& 1 \end{pmatrix} \begin{pmatrix} m_0 \\ m_1 \\ m_2 \\ \vdots \\ m_{n-1} \\ m_n \end{pmatrix} = \begin{pmatrix} f_0' \\ c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ f_n' \end{pmatrix}$$

左端的行列式不为零，方程有唯一解，解出$m_i$代入$p(x)$后就得到所求多项式。

参考资料

1. 黄云清, 《数值计算方法》, 科学出版社, 北京, 2012-06, ISBN 978-7-0302-3428-5.