函数数值计算(02)：有理函数插值法, 样条函数插值法
有理函数插值法
上一节提到, 当函数在某个区间内的变化行为比较剧烈时, 多项式内插会出现剧烈震荡的行为. 这时候利用有理函数进行内插可能会更加合适一些, 它的构造可以采用一个分式的形式: \[ \Phi_{n,n}(x)=\frac{P_n(x)}{Q_n(x)} \]
其中\(P_n\)和\(Q_n\)都是\(n\)阶多项式, 并且它们互素.
想要完全确定两者, 需要\(2n+2\)个支撑点, 但是在\(\Phi_{n,n}\)中, 两者相除, 实际上造成了一个自由的公倍数, 也就是只需要\(2n+1\)个支撑点, 记为\((x_0,y_0),(x_1,y_1),\cdots,(x_{2n},y_{2n})\). 作为内插函数, \(\Phi_{n,n}\)需要通过所有的支撑点, 可以写成如下形式:
\[ \Phi_{n,n}(x)=\phi_0(x_0)+\frac{x-x_0}{\Phi_{n-1,n-1}(x)} \]
其中, \(\phi_0(x_0)=y_0\)以通过支撑点\(0\). 分母又可以写为:
\[ \Phi_{n-1,n-1}(x)=\phi_1(x_0,x_1)+\frac{x-x_1}{\Phi_{n-2,n-2}(x)} \]
这样, \(\phi_1(x_0,x_1)\)必须满足
\[ \phi_0(x_0)+\frac{x_1-x_0}{\phi_1(x_0,x_1)}=y_1\Rightarrow\phi_1=\frac{x_1-x_0}{y_1-\phi_0(x_0)} \]
为了形式对称性, 把\(y_1\)改写为\(\phi_0(x_1)\), 则
\[ \phi_1(x_0,x_1)=\frac{x_1-x_0}{\phi_0(x_1)-\phi_0(x_0)} \]
对于更一般的情况, 我们不妨倒转递推方向
\[ \Phi_{n-j-1,n-j-1}(x)=\frac{x-x_j}{\Phi_{n-j,n-j}(x)-\phi_j(x_0,x_1,\cdots,x_j)} \]
令\(x=x_{j+1}\), 立即得到
\[ \phi_{j+1}(x_0, x_1, \cdots, x_{j+1})=\frac{x_{j+1}-x_{j}}{\phi_j(x_0,x_1,\cdots,x_{j-1},x_{j+1})-\phi_j(x_0,x_1,\cdots,x_{j-1},x_{j}} \]
它给出了各阶\(\phi\)函数的递推关系. 而\(\Phi_{n,n}\)写成连分数的形式:
\[ \Phi_{n,n}(x)=\phi_0(x_0)+\frac{x-x_0}{\phi_1(x_0,x_1)+\displaystyle\frac{x-x_1}{\phi_2(x_0,x_1,x_2)+\displaystyle\frac{x-x_2}{\ddots+\displaystyle\frac{x-x_{2n-1}}{\phi_{2n}(x_0,\cdots,x_{2n})}}}} \]
如果我们在一个单调的区间插值, 那么上述构造不会发散. 特别是对于 Runge 函数这种有理分式, 把\(x^2\)看作自变量, 只需构造\(\Phi_{1,1}\)形式的有理内插函数, 取三个支撑点就能完成内插, 并且和原函数一模一样, 所谓的 Runge 现象自然也就不再出现.
样条函数插值法
样条函数的基本思想是在每两个支撑点之间用一段多项式函数来表示, 最常用的是三次样条函数, 这里也仅介绍三次样条函数.
基本准备
在区间\([a,b]\)上包括两端一共有\(n+1\)个支撑点: \((x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\). 样条函数记作\(S(x), x\in[a,b]\).
每个子区间都具有一个三次多项式作为样条函数, 也就是\(4n\)个待定参数. 可用的方程又有多少个呢?
每个内部支撑点限制左右两个子区间的函数值, 给出两个方程, 而外部两个支撑点各自只给出一个限制, 一共\(2n\)个方程;
每个内部支撑点限制左右导函数相等, 一共\(n-1\)个方程;
每个内部支撑点限制左右二阶导数相等, 一共\(n-1\)个方程.
我们现在有\(4n-2\)个限制条件可用, 得到的样条函数具有\(2\)个自由度. 自由度的来源无非是端点处边界条件不够完备, 我们只知道边界处的函数值, 不知道边界处的导函数或者是二阶导函数的取值, 因此需要再引入两个独立的约束方程, 具体形式可以研究问题的数学或物理性质得到.
计算参数
定义二阶导数在各个支撑点处的取值为 矩:
\[ M_j=S''(x_j)\ ,\ j =0,1,\cdots,n \]
三次多项式的二阶导数是一次函数, 因而对于每一段\([x_j,x_{j+1}], \ j=0,1,\cdots,n-1\), 不难写出\(S''\)的表达式
\[ S''(x)=\frac{M_j(x_{j+1}-x)+M_{j+1}(x-x_j)}{x_{j+1}-x_j} \]
连续两次积分, 得到
\[ S(x)=\frac{M_j(x_{j+1}-x)^3}{6(x_{j+1}-x_j)}+\frac{M_{j+1}(x-x_j)^3}{6(x_{j+1}-x_j)}+A_j(x-x_j)+B_j\ , x\in[x_j,x_{j+1}] \]
这里\(A, B, M\)一共是\(3n+1\)个参数[1], 考虑\(x_j\)处函数值连续, 有
\[ \frac{M_j(x_{j+1}-x_j)^2}{6}+B_j=\frac{M_j(x_{j}-x_{j-1})^2}{6}+A_{j-1}(x_j-x_{j-1})+B_{j-1}=y_j \]
得到
\[ \boxed{\begin{aligned} B_j&=y_j-\frac{M_j(x_{j+1}-x_j)^2}{6} \newline A_j&=\frac{y_{j+1}-y_j}{x_{j+1}-x_j}-\frac{M_{j+1}-M_j}{6}(x_{j+1}-x_j) \end{aligned}} \]
现在只剩\(n+1\)个待定参数了, 同样地, 内部限制条件只有\(n-1\)个支撑点处的导数连续.
\[ S'(x)=-\frac{M_j(x_{j+1}-x)^2}{2(x_{j+1}-x_j)}+\frac{M_{j+1}(x-x_j)^2}{2(x_{j+1}-x_j)}+A_j\ , x\in[x_j,x_{j+1}] \]
考虑\(x_j\)处的导函数连续, 有
\[ -\frac{M_j(x_{j+1}-x_j)}{2}+A_j=\frac{M_j(x_j-x_{j-1})}{2}+A_{j-1} \]
代入整理, 最终得到
\[ \boxed{ \frac{x_{j+1}-x_j}{6}M_{j+1}+\frac{x_{j+1}-x_j}{3}M_j+\frac{x_j-x_{j-1}}{6}M_{j-1}=\frac{y_{j+1}-y_j}{x_{j+1}-x_j}-\frac{y_j-y_{j-1}}{x_j-x_{j-1}} } \]
这些式子给出全部的支撑点的约束方程, 只要再给边界上[2]引入两个独立约束即可.