张瑞
中国科学技术大学数学科学学院
rui@ustc.edu.cn
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$ 时,等距节点的插值作图
可以看到,插值多项式在端点附近抖动最大。
定义 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
可以看到,插值多项式越高,误差越小。
要得到插值效果好的高次插值多项式,需要插值节点取的很特殊。
但是这样的节点在实际中应用不太方便,操作复杂。
实际应用中,等距节点最普遍,操作简单,但效果不好。
如何克服Runge现象 ,在节点越多的情形下,得到插值误差越小的插值函数呢?
使用低阶的多项式做插值,以避免高次多项式的Runge现象
使用分段处理的方法,每次处理少数的点,最终得到分段的插值函数
分段插值
在节点组$\{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}
\]
整个函数是分段函数。
每一个小区间内,是一个次数不超过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$ 。也就是说,节点越密,误差越小
可以通过增加每次插值使用的区间数(即节点数)来提高每个区间上的多项式的次数,
但不能改变插值函数整体只有连续性,没有可导性 。
如果想保证整体具有1阶导数,可选的办法是定义每个节点处的1阶导数值,
这样在每个区间上可以构造3次的Hermite插值。
问题 . 如何以最小的代价,提高整条插值曲线的光滑性?
样条插值
定义 2 .
若分段函数在每个小区间上是次数不超过3次的多项式,且具有2阶连续导数,则称曲线为3次样条函数 。
满足插值条件的3次样条函数则称为3次样条插值函数 。
有节点$(x_i, f(x_i))$ , $i=0,1,\cdots,n$ 。如何构造3次样条函数$S(x)$ 满足
\[S(x_i)=f(x_i), \quad i=0,1,\cdots,n
\]
设$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$ ,常见的三种边界条件 是
若$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)$ 的对角占优的三对角阵,
因此,解存在唯一。
若$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}
\]
得到的仍然是一个对角占优的三对角的线性方程组。
周期边界。 原函数$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次样条函数插值的特点:
每个小区间的多项式次数不高,因此不存在Runge现象。
可以得到随着点越来越密,误差也会越来越小。
3次样条函数插值具有“整体性”。如果一个节点处的函数值发生了变化,所有区间上的系数都会产生变化。
注 .
程序使用了 python 的 scipy
库来计算样条。