分段插值

插值

张瑞
中国科学技术大学数学科学学院

分段插值

Runge现象

插值的次数是越高越好吗?

例 1. 函数$f(x)=\frac1{1+x^2}, x\in[-5,5]$做等距节点插值,给出最大误差

. 在节点$x_i=-5+\frac{i}{10}$, $i=0,1,\cdots,n$上插值,编程计算

=============== n= 5
Max error at point  0.0  with error  0.4326923076923076
=============== n= 10
Max error at point  -4.7  with error  1.9156430502192496
=============== n= 20
Max error at point  4.870  with error  59.768327839903954
=============== n= 40
Max error at point  4.950  with error  104371.93698015573

可以看到,次数越高,误差越大

$n=10$时,等距节点的插值作图

intp-runge

可以看到,插值多项式在端点附近抖动最大。

定义 1.
等距节点下,高次插值多项式在端点附近的抖动现象,称为Runge现象

. 是否高次多项式一定会有Runge现象呢?

用Chebyshev点,可以得到

=============== n= 5
Max error at cheby point  0.0  with error  0.5559113388123956
=============== n= 10
Max error at cheby point  -0.7800000000000002  with error  0.10914672464976671
=============== n= 20
Max error at cheby point  -1.1099999999999999  with error  0.01533291731815517
=============== n= 40
Max error at cheby point  0.9500000000000002  with error  0.00028938783946064195

可以看到,插值多项式越高,误差越小。

. 插值的效果与插值节点的选取有关。

  • 要得到插值效果好的高次插值多项式,需要插值节点取的很特殊。
  • 但是这样的节点在实际中应用不太方便,操作复杂。
  • 实际应用中,等距节点最普遍,操作简单,但效果不好。

分段插值

在节点组$\{x_i\}_{i=0}^n$上构造它的分段线性插值多项式。

在区间$[x_{i},x_{i+1}]$上,构造线性插值多项式,得到

\[p_i(x)=f(x_i)\frac{x-x_{i+1}}{x_i-x_{i+1}}+f(x_{i+1})\frac{x-x_i}{x_{i+1}-x_i} \]

这样,整个插值曲线是个分段函数

\[g(x)=\begin{cases} p_i(x), \quad x\in[x_i,x_{i+1}] \\ \end{cases} \]

intp-piecewise-linear

  • 整个函数是分段函数。
  • 每一个小区间内,是一个次数不超过1次的多项式。
  • 在节点上,只能保证连续性,不能保证可导性。

误差

在每个小区间$[x_i, x_{i+1}]$内,有误差表达式

\[f(x)-p_i(x)=\frac{f''(\xi_i)}{2!}(x-x_i)(x-x_{i+1}) , \xi_i\in(x_i,x_{i+1}) \]

因此

\[\begin{aligned} |f(x)-g(x)|\leq &\max_{i=0,\cdots,n-1} |\frac{f''(\xi_i)}{2!}(x-x_i)(x-x_{i+1})| \\ \leq &\max_{i=0,\cdots,n-1} |\frac{f''(\xi_i)}{2!}|\left(\frac{x_{i+1}-x_i}2\right)^2 \end{aligned} \]

$M_2=\max(|f''(x)|)$, $h=\max(x_{i+1}-x_i)$,则有

\[|f(x)-p(x)|\leq \frac{M_2}8 h^2 \]

分段线性插值的误差满足

\[R=|f(x)-p(x)|\leq \frac{M_2}8 h^2 \]

可以看到$\displaystyle\lim_{h\to0}R=0$。也就是说,节点越密,误差越小

问题. 如何以最小的代价,提高整条插值曲线的光滑性?

样条插值

定义 2.
若分段函数在每个小区间上是次数不超过3次的多项式,且具有2阶连续导数,则称曲线为3次样条函数

满足插值条件的3次样条函数则称为3次样条插值函数

$S(x)$在区间$[x_i, x_{i+1}]$上的表达式$s_i(x)$

\[\begin{aligned} s_i(x) = a_i + b_i x+c_i x^2+d_i x^3,& x\in[x_i, x_{i+1}] , \\ & i=0,1,\cdots,n-1 \end{aligned} \]

需要求出$\{a_i, b_i, c_i, d_i\}$$4n$个系数。

由函数特性和插值条件,可得

\[\begin{aligned} s_i(x_i)=f(x_i), &i=0,1,\cdots, n-1 \\ s_i(x_{i+1})=f(x_{i+1}),& i=0,1,\cdots,n-1 \\ s'_i(x_i)=s'_{i-1}(x_i),& i=1,2,\cdots,n-1 \\ s''_i(x_i)=s''_{i-1}(x_i),& i=1,2,\cdots,n-1 \\ \end{aligned} \]
  • 可以看到,共有$4n$个系数待求,但只有$4n-2$个方程。
  • 还需要增加2个条件,才可能唯一确定这个插值函数。增加的条件通常会在边界处给出,也称为边界条件
  • 这样的做法需要解含$4n$个未知量的方程,比较复杂

m关系式

记每个节点处的1阶导数为$m_i, i=0,1,\cdots,n$, 则区间$[x_i, x_{i+1}]$上的表达式$s_i(x)$满足

\[\begin{cases} s_i(x_i)=f(x_i),& s_i(x_{i+1})=f(x_{i+1}) \\ s'_i(x_i)=m_i,& s'_i(x_{i+1})=m_i \end{cases} \]

这是一个4个条件的Hermite插值,正好可以确定一个3次多项式。

可以得到

\[\begin{aligned} s_i(x)=&\frac1{h_i^3}(h_i+2(x-x_i))(x-x_{i+1})^2f(x_i) \\ &+\frac1{h_i^3}(h_i-2(x-x_{i+1}))(x-x_i)^2f(x_{i+1}) \\ &+\frac1{h_i^2}(x-x_i)(x-x_{i+1})^2m_i \\ &+\frac1{h_i^2}(x-x_{i+1})(x-x_i)^2m_{i+1} \\ \end{aligned} \]

由3次样条函数具有2阶导数,因此

\[s''_i(x_i)=s''_{i-1}(x_i), i=1,2,\cdots,n \]

\[\begin{aligned} s''_i(x)=&\frac6{h_i^3}(h_i+2(x-x_{i+1}))f(x_i) \\ &+\frac6{h_i^3}(h_i-2(x-x_i))f(x_{i+1}) \\ &+\frac2{h_i^2}(3(x-x_{i+1})+h_i)m_i \\ &+\frac2{h_i^2}(3(x-x_i)-h_i)m_{i+1} \\ \end{aligned} \]
\[\begin{aligned} s''_i(x_i)=&\frac6{h_i^3}(h_i+2(x_i-x_{i+1}))f(x_i) \\ &+\frac6{h_i^3}(h_i-2(x_i-x_i))f(x_{i+1}) \\ &+\frac2{h_i^2}(3(x_i-x_{i+1})+h_i)m_i \\ &+\frac2{h_i^2}(3(x_i-x_i)-h_i)m_{i+1} \\ =&\frac6{h_i}\frac{f(x_{i+1})-f(x_i)}{h_i} \\ &+\frac2{h_i}(-2m_i-m_{i+1}) \\ \end{aligned} \]
\[\begin{aligned} s''_i(x_{i+1})=&\frac6{h_i^3}(h_i+2(x_{i+1}-x_{i+1}))f(x_i) \\ &+\frac6{h_i^3}(h_i-2(x_{i+1}-x_i))f(x_{i+1}) \\ &+\frac2{h_i^2}(3(x_{i+1}-x_{i+1})+h_i)m_i \\ &+\frac2{h_i^2}(3(x_{i+1}-x_i)-h_i)m_{i+1} \\ =&-\frac6{h_i}\frac{f(x_{i+1})-f(x_i)}{h_i} \\ &+\frac2{h_i}(m_i+2m_{i+1}) \\ \end{aligned} \]

整理可得,对$i=1,2,\cdots,n-1$,有

\[\begin{aligned} \frac{h_{i}}{h_i+h_{i-1}}m_{i-1}+&2m_i+\frac{h_{i-1}}{h_i+h_{i-1}}m_{i+1} \\ =3\bigg(&\frac{h_{i}}{h_i+h_{i-1}}\frac{f(x_{i})-f(x_{i-1})}{h_{i-1}} \\ &+\frac{h_{i-1}}{h_i+h_{i-1}}\frac{f(x_{i+1})-f(x_{i})}{h_{i}} \bigg) \end{aligned} \]

\[\begin{aligned} \displaystyle\lambda_i=\frac{h_{i}}{h_i+h_{i-1}}>0, \quad \mu_i=1-\lambda_i \\ c_i=3(\lambda_if[x_{i-1},x_i]+\mu_if[x_i, x_{i+1}]) \end{aligned} \]

可得线性方程组

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

这是一个$n+1$个未知量,$n-1$个方程的线性方程组。还需要2个条件就有可能得到唯一的解。

不防设$x_0<x_1<\cdots<x_n$,常见的三种边界条件

  1. $x_0$$x_n$处的1阶导数$m_0$, $m_n$已知,则有
    \[\begin{pmatrix} 2 & \mu_1 \\ \lambda_2 & 2 & \mu_2 \\ & \ddots & \ddots & \ddots \\ & & \lambda_{n-1} & 2 \end{pmatrix} \begin{pmatrix} m_1 \\ m_2 \\ \vdots \\ m_{n-1} \end{pmatrix} =\begin{pmatrix} c_1-\lambda_1m_0 \\ c_2 \\ \vdots \\ c_{n-1}-\mu_nm_n \end{pmatrix} \]
    很明显,$\lambda_i, \mu_i\in(0,1)$,因此矩阵是一个$(n-1)\times(n-1)$的对角占优的三对角阵, 因此,解存在唯一。
  1. $x_0$$x_n$处的2阶导数给定,分别是$M_0$, $M_n$,则有
    \[\begin{aligned} M_n=&s''_{n-1}(x_n)=-\frac{6}{h_{n-1}}f[x_{n-1},x_n]+\frac2{h_{n-1}}(m_{n-1}+2m_n) \\ M_0=&s''_{0}(x_0)=\frac{6}{h_{0}}f[x_{0},x_1]-\frac2{h_{0}}(m_{1}+2m_0) \end{aligned} \]

整理后,得到条件

\[\begin{aligned} m_{n-1}+2m_n =& 3f[x_{n-1},x_n]+\frac{h_{n-1}}2M_n \\ 2m_0+m_{1} = &3f[x_0,x_1]-\frac{h_0}2M_0 \end{aligned} \]
\[\begin{pmatrix} 2 & 1 & \\ \lambda_1 & 2 & \mu_1 \\ & \lambda_2 & 2 & \mu_2 \\ & & \ddots & \ddots & \ddots & \\ & & & \lambda_{n-1} & 2 & \mu_{n-1} \\ & & & & 1 & 2 \end{pmatrix} \begin{pmatrix} m_0 \\ m_1 \\ m_2 \\ \vdots \\ m_{n-1} \\ m_n \end{pmatrix} =\begin{pmatrix} 3f[x_0,x_1]-\frac{h_0}2M_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{n-1} \\ 3f[x_{n-1},x_n]+\frac{h_{n-1}}2M_n \end{pmatrix} \]

得到的仍然是一个对角占优的三对角的线性方程组。

  1. 周期边界。 原函数$f(x)$是一个周期函数。 若函数的周期正好是$x_n-x_0$,则此时应该有$f(x_0)=f(x_n)$, $f'(x_0)=f'(x_n)$, $f''(x_0)=f''(x_n)$。 然后,由
    \[\begin{aligned} m_0=&m_n \\ s''_0(x_0)=&s''_{n-1}(x_n) \end{aligned} \]
    可以得到周期边界的表达式。

具体的作为作业

M关系式

类似m关系式,假定每个节点上的2阶导数是$M_i$, $i=0,1,\cdots,n$。 则在每个区间$[x_{i}, x_{i+1}]$上的函数$s_i(x)$满足

\[\begin{aligned} s_i(x_i)=f(x_i), & s_i(x_{i+1})=f(x_{i+1}) \\ s''_i(x_i)=M_i, & s''_i(x_{i+1})=M_{i+1} \end{aligned} \]

4个条件,仍然有可能确定一个3次多项式。

  • 可以假定
    \[s_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3 \]
    代入条件,求出系数$a_i, b_i, c_i, d_i$

还有更直观的方法可以使用

注意到$s_i(x)$是3次多项式,求二阶导后,得到1次多项式。 由$s''_i(x)$满足的插值条件,可得

\[s''_i(x)=M_i\frac{x-x_{i+1}}{x_i-x_{i+1}}+M_{i+1}\frac{x-x_i}{x_{i+1}-x_i} \]

积分2次后,有

\[s_i(x)=\frac{M_{i+1}}{6h_i}(x-x_i)^3+\frac{M_i}{6h_i}(x_{i+1}-x)^3+ax+b \]

$a$, $b$为积分常数。利用插值条件来确定$a$$b$

为计算方便,可以写成

\[\begin{aligned} s_i(x)=&\frac{M_{i+1}}{6h_i}(x-x_i)^3+\frac{M_i}{6h_i}(x_{i+1}-x)^3 \\ &+c(x_{i+1}-x)+d(x-x_i) \end{aligned} \]

由3次样条函数1阶导数的存在性,有

\[s'_i(x_i)=s'_{i-1}(x_i) , i=1,2,\cdots,n-1 \]

整理后,得到

\[\mu_iM_{i-1}+2M_i+\lambda_iM_{i+1}=d_i, i=1,2,\cdots,n-1 \]

其中$d_i=6f[x_{i-1},x_i,x_{i+1}]$

边界条件的表达式请参考教材,这里不再细述。

作业: 推导m关系式和M关系式在周期边界条件($f(x)$的周期是$x_n-x_0$)下的表达式。

3次样条函数插值的特点:

  1. 每个小区间的多项式次数不高,因此不存在Runge现象。
  2. 可以得到随着点越来越密,误差也会越来越小。
  3. 3次样条函数插值具有“整体性”。如果一个节点处的函数值发生了变化,所有区间上的系数都会产生变化。

intp-spline

. 程序使用了 python 的 scipy 库来计算样条。

谢谢

vertical slide 2

目录

本节读完

例 2.

[#ex9-1-0].