数值分析：插值与拟合

1 插值

定义设函数 $y=f(x)$ 在区间 $[a,b]$ 上的 $n+1$ 个点， $a\leq x_0,x_1,...,x_n\leq b$ 上的函数值为 $y_0,y_1,...,y_n$ ，若粗在函数 $P(x)$ ，使 $P(x_i)=y_i,i=0,1,...n$ 成立，则称函数 $P(x)$ 为 $f(x)$ 的 插值函数， $f(x)$ 称为被插值函数，点 $x_0,x_1,...,x_n$ 称为 插值节点，包含插值节点的区间 $[a,b]$ 称为 插值区间

在这里我们进讨论 多项式插值 与 分段插值

多项式插值

插值多项式是否存在，若存在是否唯一
怎么推导插值多项式
如何估计其逼近程度
如何应用

定理1 在 $n+1$ 个相异插值节点 $x_0,x_1,...,x_n$ 处取给定值 $y_0,y_1,...,y_n$ 的次数不高于 $n$ 的插值多项式存在且唯一（可以通过非齐次线性方程组的范德蒙德行列式证明）

2 拉格朗日插值

2.1 拉格朗日插值多项式

线性插值
过曲线 $y=f(x)$ 上两点 $M_0(x_0,y_0),M_1(x_1,y_1)$ 做直线 $L_1(x)$ ，由两点式： $L_1(x)=y_0\frac{x-x_1}{x_0-x_1}+y_1\frac{x-x_0}{x_1-x_0}$ 显然 $L_1(x)$ 为插值多项式
抛物线插值
过曲线 $y=f(x)$ 上三点 $M_0(x_0,y_0),M_1(x_1,y_1),M_2(x_2,y_2)$ 做抛物线 $L_2(x)$ ，同样不难得到： $L_2(x)=y_0\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+y_1\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}+y_2\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$ 显然 $L_2(x)$ 也为插值多项式

用类似的推导方式可以证明， $n+1$ 个节点的拉格朗日插值多项式应定义为如下形式。
定义 $y=f(x)$ 在插值节点 $x_0,x_1,...,x_n$ 上的 拉格朗日插值多项式 $L_n(x)$ 为 $L_n(x)=y_0l_0(x)+y_1l_1(x)+...+y_nl_n(x)=\sum\limits_{i=0}^ny_il_i(x)$ 其中 $l_i(x)=\frac{x-x_0}{x_i-x_0}\frac{x-x_1}{x_i-x_1}...\frac{x-x_n}{x_i-x_n},i=0,1,...n$ $l_i(x)$ 称为 插值基函数

2.2 插值余项

若在 $[a,b]$ 上使用 $L_n(x)$ 近似 $f(x)$ ，则其截断误差为 $R_n(x)=f(x)-L_n(x)$ ，也称为 插值多项式的余项。关于插值余项估计有如下定理

定理2 设 $f^{(n)}(x)$ 在区间 $[a,b]$ 上连续， $f^{(n+1)}(x)$ 在区间 $(a,b)$ 内存在， $L_n(x)$ 是以 $a\le x_0<x_1<...<x_n\le b$ 为节点的拉格朗日插值多项式，则对任何 $x\in [a,b]$ $R_n(x)=f(x)-L_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)$ $\omega_{n+1}(x)=(x-x_0)(x-x_1)...(x-x_n)$
这里 $\xi\in(a,b)$ 且依赖于 $x$

2.3 算法步骤

(1) 输入 $x,x_i,y_i,(i=0,1,...,n)$
(2) 对 $i=0,1,...,n$ 计算 $l_i=\prod\limits_{j=o,j\neq i}^n\frac{x-x_j}{x_i-x_j}$
(3) 计算 $L=\sum\limits_{i=0}^ny_il_i$
(4) 输出 $L \approx f(x)$ ，结束

2.4 注意事项

$L_n(x)$ 为次数不高于 $n$ 次的多项式，其次数可能小于 $n$
若 $f(x)$ 是次数不超过 $n$ 次的多项式，那么以 $n+1$ 个点为节点的插值多项式就一定是其本身。特别是当取 $f(x)\equiv 1$ 时，得恒等式 $\sum\limits_{i=0}^nl_i(x)\equiv 1$ 这一等式为我们提供了验证插值基函数的一种方法
$R_n(x)$ 用于误差估计。若 $|f^{(n+1)}(x)|\le M$ ，则 $|R_n(x)|\le \frac{M}{(n+1)!}|\omega_{n+1}(x)|$
$L_n(x)$ 只与节点及函数 $f(x)$ 在节点处的值有关，与 $f(x)$ 无关；而 $R_n(x)$ 与 $f(x)$ 关系最为密切。

3 差商与牛顿插值

3.1 差商

定义设已给插值节点 $x_0,x_1,...,x_n$ 以及相应的函数值 $f(x_0),f(x_1),...,f(x_n)$ 称 $f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}$ 为函数 $y=f(x)$ 关于点 $x_0,x_1$ 的一阶差商，称 $f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}$ 为函数 $y=f(x)$ 关于点 $x_0,x_1,x_2$ 的二阶差商，称 $f[x_0,x_1,...,x_k]=\frac{f[x_1,x_2,...,x_k]-f[x_0,x_1,...,x_{k-1}]}{x_k-x_0}$ 为函数 $y=f(x)$ 关于点 $x_0,x_1,x_2$ 的 $k$ 阶差商

3.2 差商的性质

线性性
$k$ 阶差商 $f[x_0,x_1,...,x_k]$ 是函数值 $f(x_0,f(x_1),...,f(x_k)$ 的线性组合，即 $f[x_0,x_1,...,x_k]=\sum\limits_{i=0}^k{\frac{f(x_i)}{\omega'_{k+1}(x_i)}}$ 其中 $\omega_{k+1}(x)=\prod\limits_{i=0}^k(x-x_i)$
对称性
差商 $f[x_0,x_1,...,x_k]$ 是 $x_0,x_1,...,x_k$ 的对称函数。
差商与导数的关系
设函数 $f(x)$ 在区间 $[a,b]$ 上存在 $n$ 阶导数，且 $x_0,x_1,...,x_k\in [a,b]$ ，则存在 $\xi \in [a,b]$ 使得 $f[x_0,x_1,...,x_k]=\frac{f^{(n)}(\xi)}{n!}$

3.3 牛顿插值多项式

由差商的定义可导出 $f(x)=N_n(x)+R_n(x)$ 其中 $\begin{align} N_n(x)=&f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]+...\\ &+(x-x_0)(x-x_1)...(x-x_{n-1})f[x_0,x_1,...,x_n]\\ \end{align}$ 称为牛顿插值多项式 $R_n(x)=f[x,x_0,x_1,...,x_n]\omega_{n+1}(x)$ 称为牛顿插值多项式的余项

3.4 算法步骤

(1) 输入 $x,x_i,y_i(i=0,1,...,n)$
(2) 对 $i=0,1,...,n$ 计算 $f_i=y_i$
(3)对 $i=1,2,...,n$ 计算 $f_k=\frac{f_{k-1}-f_k}{x_{k-1}-x_k}(k=n,n-1,...,i)$
(3) 计算 $p=f_0+\sum\limits_{k=0}^n{f_k(\prod\limits_{j=0}^{k-1}(x-x_j))}$
(4) 输出 $p \approx f(x)$ ，结束

4 Hermite插值

4.1 Hermite插值多项式及其余项

Hermite插值多项式: $H_{2n+1}(x)=\sum\limits_{j=0}^n[y_j\alpha_j(x)+y_j'\beta_j(x)]$ 其中 $\alpha_j(x)=\left[1-2(x-x_j)\sum\limits_{k=0,k\neq j}^n\frac{1}{x_j-x_k}\right]l^2_j(x)$ $\beta_j(x)=(x-x_j)l_j^2(x)$ 余项: $R_{2n+1}(x)=\frac{f^{(2n+2)}}{2n+2}\omega^2_{n+1}(x)$ 其中 $\omega_{n+1}(x)=\prod\limits_{i=0}^n(x-x_i),\xi\in [a,b]$

4.2 两点三次Hermite插值多项式

已知 $y=f(x)$ 在 $[a,b]$ 上的节点 $x_0,x_1$ 上的函数值 $y_0,y_1$ 及一阶导数值 $y_0',y_1'$ ，则三次Hermite插值多项式 $\begin{align} H_3(x)=&y_0\alpha_0(x)+y_1\alpha_1(x)+y_0'\beta_0(x)+y_1'\beta_1(x)\\ =&y_0\left(1-2\frac{x-x_0}{x_0-x_1}\right)\left(\frac{x-x_1}{x_0-x_1}\right)^2+y_1\left(1-2\frac{x-x_1}{x_1-x_0}\right)\left(\frac{x-x_1}{x_0-x_1}\right)^2\\ &+y_0'(x-x_0)\left(\frac{x-x_1}{x_0-x_1}\right)^2+y_1'(x-x_1)\left(\frac{x-x_0}{x_1-x_0}\right)^2\\ \end{align}$ 余项 $R_3(x)=\frac{f^{(4)}(\xi)}{4!}(x-x_0)^2(x-x_i)^2,\xi\in[a,b]$

5 分段低次插值

5.1 龙格现象

随着插值多项式次数的增大，插值函数 $L_n(x)$ 在两端会发生激烈的震荡，这就是龙格现象

5.2 分段线性插值

对给定区间 $[a,b]$ 做分割： $a=x_0<x_1<...<x_n=b$ ，在每个小区间 $[x_i,x_{i+1}]$ 上以 $x_i,x_{i+1}$ 为节点作 $f(x)$ 的线性插值： $S_i(x)=\frac{x-x_{i+1}}{x_i-x_{i+1}}f(x_i)+\frac{x-x_i}{x_{i+1}-x_i}f(x_{i+1}),x\in[x_i,x_{i+1}]$ 把每个小区间上的线性插值函数连起来，就得到了分段线性插值函数 $S(x)$
误差公式为 $|f(x)-S(x)|\leq\frac{M_2}{8}(x_i-x_{i+1})^2,M_2=\max\limits_{a\leq x\leq b}|f''(x)|$

5.3 分段三次Hermite插值

对给定区间 $[a,b]$ 做分割： $a=x_0<x_1<...<x_n=b$ ，在每个小区间 $[x_i,x_{i+1}]$ 上以 $x_i,x_{i+1}$ 为节点作分段三次Hermite插值： $\begin{align} S_i(x)=&f(x_i)\left(1-2\frac{x-x_i}{x_i-x_{i+1}}\right)\left(\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)^2\\ &+f(x_{i+1})\left(1-2\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)\left(\frac{x-x_i}{x_{i+1}-x_{i}}\right)^2\\ &+f'(x_i) (x-x_i)\left(\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)^2\\ &+f'(x_{i+1}) (x-x_{i+1})\left(\frac{x-x_i}{x_{i+1}-x_{i}}\right)^2\\ \end{align}$ 误差公式为 $|f(x)-S(x)|\leq\frac{h^4}{384}\max\limits_{x_i\leq x \leq x_{i-1}}|f^{(4)}(x)|,h=\max\limits_{0\leq i \leq n-1}|x_{i+1}-x_i|$

5 最小二乘拟合

定义确定 $a_k(k=0,1,...,m$ 使得 $R(a_0,a_1,...,a_m)=\sum\limits_{k=1}^n(\varphi(x_k)-y_k)^2$ 最小问题称为观测数据的最小二乘拟合问题。若函数系 $\{b_k(x)\}$ 为幂函数系 $\{x^k\}$ ，此时的到的 $\varphi(x)$ 称为回归曲线

5.1 最小二乘法

由于 $R({{a}_{0}},{{a}_{1}},...,{{a}_{m}})\ge 0$ ，且为连续函数，故一定存在一组 ${{a}_{k}}(k=0,1,...,m)$ 使其达到极小值，这时我们只需要求出驻点即可。
因此我们假设 $y=f(x)$ 的观测数据为 $({{x}_{k}},{{y}_{k}})(k=1,2,...,n)$ ，对于 $\varphi(x)$ 为 $m$ 次回归曲线，通过求解正规方程组可以求出 ${{a}_{k}}(k=0,1,...,m)$ 的值，可求得近似函数 $f(x)$ 。 $\left( \begin{matrix} n & \sum\limits_{k=1}^n x_k & \sum\limits_{k=1}^n x^2_k & \cdots & \sum\limits_{k=1}^n x^m_k\\ \sum\limits_{k=1}^n x_k & \sum\limits_{k=1}^n x^2_k & \sum\limits_{k=1}^n x^3_k & \cdots & \sum\limits_{k=1}^n x^{m+1}_k\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \sum\limits_{k=1}^n x^m_k & \sum\limits_{k=1}^n x^{m+1}_k & \sum\limits_{k=1}^n x^{m+2}_k& \cdots & \sum\limits_{k=1}^n x^{2m}_k\\ \end{matrix} \right) \left( \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_m \\ \end{matrix} \right)= \left( \begin{matrix} \sum\limits_{k=1}^n y_k \\ \sum\limits_{k=1}^n x_ky_k \\ \vdots \\ \sum\limits_{k=1}^n x_k^my_k \\ \end{matrix} \right)$

解题时，先写出正规方程组所需数据的数据表，再求解即可得到参数

5.2 最小二乘法的病态现象

衡量一个线性方程组是否“病态”及其病态程度，可通过矩阵的条件数理论来完成。若线性方程组系数矩阵 $A$ 的条件数 $Cond(A)\gg 1$ ，则该方程组“病态”，并且系数矩阵的条件数越大，方程组的“病态”程度就越严重，其解随系数矩阵或自由项的变化（灵敏度）就越大。