【笔记】大一下数值分析碎碎念——插值
\(\newcommand\op[1]{\operatorname{#1}}\)
插值
给定数据点 \((x_i,y_i)\),要求找到函数满足 \(f(x_i)=y_i\)。
线性插值:全局信息维护,光滑性(求导),积分都不太好搞。但是原理简单。
多项式?
指数?变化快。
三角函数?周期性。
多项式插值
Weierstrass 逼近定理:设 \(f\in \op C[a,b]\),则 \(\forall \epsilon > 0\),存在多项式 \(p(x)\) 满足 \(|f(x)-p(x)| < \epsilon\) over \([a,b]\)
过 \(n\) 个点 \((x_i,y_i)\),假设多项式是 \(n-1\) 次的,则正好可以唯一解出所有系数。注意到只要有 \(x_i\) 两两不等,则必然有解,因为系数矩阵是范德蒙矩阵。
Vandermonde 矩阵的行列式是 \(\prod_{j < i}(x_i-x_j)\)
证明:
数学归纳法,\(n=1,n=2\) 显然。
\(n>2\) 时考虑每一行减去上一行的 \(x_1\) 倍,再按第一行展开,提取公因式,由此可知。
Lagrange 插值
对于每一个点 \((x_i,y_i)\),构造出多项式函数 \(l_i(x)\) 满足在 \(x_j\)(\(j \ne i\))时取 \(0\),\(x_i\) 处取 \(1\)。然后我们把所有函数各自乘上 \(y_i\) 加起来 \(f(x) = \sum l_i(x)y_i\) 即得到插值函数。
显然存在构造:
逐次线性插值(Neville 算法)
记 \(I_{i_1,\cdots,i_n}(x)\) 为满足 \((x_{i_j},y_{i_j})\) 的多项式,那么我们可以递归计算:
验证一下:显然在 \(x=x_0,\cdots,x_{k}\) 正确。在 \(x=x_l\) 处显然也正确!每步次数最多 \(+1\)。
但是直接实现复杂度爆炸,把指标换一下改成:
这下指标一定是连续子段,所以计算过程中被用到的 \(I\) 不会超过 \(O(k^2)\) 级别。
单点求值时,如需增加一个点,复杂度为 \(O(k)\),空间复杂度维持 \(O(k)\)。
优势:可以逐步比较精度。
Newton 插值
求出一组系数 \(c_n\) 使得:
写一些值看看:
\(P_n(x_0)= y_0 = c_0\)
\(P_n(x_1)=y_1 = c_0 + c_1(x_1-x_0)\) 得到 \(c_1 = \frac{y_1 - y_0}{x_1-x_0}\)
\(P_n(x_2) = y_2 = c_0 + c_1(x_2-x_0)+c_2(x_2-x_0)(x_2-x_1)\) 得到 \(c_2 = \frac{\frac{y_2-y_0}{x_2-x_0}- \frac{y_1-y_0}{x_1-x_0}}{x_2-x_1}\)
我们先记
于是我们就可以递归计算:
发现 \(c_k = f[x_0,\cdots,x_k]\)。
所以我们需要想办法计算 \(f[x_0,\cdots,x_k]\)。
可以数学归纳法证明:
于是可知顺序不会影响差商,故我们可以控制参数一定是一个区间,那么又可以计算了。
另一种证明 \(c_k = f[x_0,\cdots,x_k]\) 的方法:
首先我们定义 \(f[x_0,\cdots,x_n]\) 过这些点的多项式最高次项的系数,最后我们证明这个东西就是我们上面定义的差商。接下来我们证明:
差分
\(\Delta f_k = f_{k+1} - f_k\)
\(\nabla f_k = f_k - f_{k-1}\)
\(\delta f_k = f_{k + \frac{1}{2}} - f_{k-\frac{1}{2}}\)
误差
内差。
龙格现象
Hermite 插值
过 \((x_1,y_1),\cdots,(x_n,y_n)\) 且告诉导数 \(y_1',\cdots,y_n'\),则得到 \(2n-1\) 阶插值多项式为 Hermite 插值多项式。构造:
检查一下满足:
- \(H_j(x_i) = \delta_{i,j}\)
- \(H_j'(x_i)=0\)
- \(\hat H_j(x_i) = 0\)
- \(\hat H _j(x_i) = \delta_{i,j}\)
分段插值
三次样条插值
分段,每一段是三次函数,且满足最后二阶连续可导。
分段处有 \(3n-3\) 个条件,过各点有 \(n+1\) 个条件,还有两个条件需要自行补充。
- 自然样条:补充两个端点处二阶导为 \(0\)。
- 曲率调整样条:直接指定两个端点处的二阶导。
- 钳制样条:指定两个断点处的一阶导。
一定有解:考虑写成矩阵,发现对角占优,所以可以求逆!