数学建模入门笔记(3) 插值与拟合

插值与拟合

插值和拟合的区别：拟合不要求过每一个已知点，而插值要求过每一个已知点，因而插值可以看作过每一个点的拟合。

插值适用于补全缺失值，因为使用一般拟合就有可能使已知值偏移，不符合需求。据说PS用某种样条插值，放大的时候最大程度的保留连续性，因此显得不是那么模糊.

数学建模笔记1.3——插值 - Infty的文章 - 知乎
https://zhuanlan.zhihu.com/p/390028714

1 插值

1.1 分段线性插值

从几何上解释就是用直线把已知点相连形成折线

数学表示

分段线性插值函数记为 \(I_n(x)\)

\[I_n(x)=\sum_{i=0}^n y_i l_i(x) \]

其中基函数记为 \(l_i(x)\)

\[l_i(x)=\left\{\begin{array}{l} \frac{x-x_{i-1}}{x_i-x_{i-1}}, x \in\left[x_{i-1}, x_i\right], i \neq 0, \\ \frac{x-x_{i+1}}{x_i-x_{i+1}}, x \in\left[x_i, x_{i+1}\right], i \neq n, \\ 0, \text { 其他 } \end{array}\right. \]

\(I_n(x)\) 有良好的收敛性, 即对于 \(x \in[a, b]\), 有

\[\lim _{n \rightarrow \infty} I_n(x)=f(x) \]

特点

计算量小，用 \(I_n(x)\) 计算 \(x\) 点的插值时，只用到 \(x\) 左右的两个节点， 计算量与节点个数 \(n\) 无关。但 \(n\) 越大插值误差越小。

不够光滑。

1.2 拉格朗日插值

数学表示

拉格朗日 (Lagrange) 插值的基函数为

\[\begin{aligned} l_i(x) & =\frac{\left(x-x_0\right) \cdots\left(x-x_{i-1}\right)\left(x-x_{i+1}\right) \cdots\left(x-x_n\right)}{\left(x_i-x_0\right) \cdots\left(x_i-x_{i-1}\right)\left(x_i-x_{i+1}\right) \cdots\left(x_i-x_n\right)} \\ & =\prod_{\substack{j=0 \\ j \neq i}}^n \frac{x-x_j}{x_i-x_j}, i=0,1, \cdots, n \end{aligned} \]

\(l_i(x)\) 是 \(n\) 次多项式, 满足

\[l_i\left(x_j\right)=\left\{\begin{array}{l} 0, j \neq i, \\ 1, j=i 。 \end{array}\right. \]

拉格朗日插值函数

\[L_n(x)=\sum_{i=0}^n y_i l_i(x)=\sum_{i=0}^n y_i\left(\prod_{\substack{j=0 \\ j \neq i}}^n \frac{x-x_j}{x_i-x_j}\right) \]

特点

\((1)\) 增加或改变已知节点则需要重新计算生成拉格朗日插值函数，计算量大。

\((2)\) 会出现龙格现象（Runge）：
指的是对于某些函数，使用均匀节点构造高次多项式差值时，在插值区间的边缘的误差可能很大的现象。它是由Runge在研究多项式差值的误差时发现的，这一发现很重要，因为它表明，并不是插值多项式的阶数越高，效果就会越好。

龙格现象(Runge Phenomenon) - sslchi的文章 - 知乎
https://zhuanlan.zhihu.com/p/138747068

1.3 牛顿插值

数学表示

差商

自变量之差与因变量之差之比叫差商
定义: 函数 \(y=f(x)\) 在区间 \(\left[x_i, x_{i+1}\right]\) 上的平均变化率

\[f\left[x_i, x_{i+1}\right]=\frac{f\left(x_{i+1}\right)-f\left(x_i\right)}{x_{i+1}-x_i} \]

称为 \(f(x)\) 关于 \(x_i, x_{i+1}\) 的一阶差商，并记为 \(f\left[x_i, x_{i+1}\right]\)
二阶差商:

\[f\left[x_i, x_{i+1}, x_{i+2}\right]=\frac{f\left[x_{i+1}, x_{i+2}\right]-f\left[x_i, x_{i+1}\right]}{x_{i+2}-x_i} \]

\(\mathrm{m}\) 阶差商:

\[f\left[x_0, x_1, \cdots, x_m\right]=\frac{f\left[x_1, x_2, \cdots, x_m\right]-f\left[x_0, x_1, \cdots, x_{m-1}\right]}{x_m-x_0} \]

牛顿插值公式

\[\begin{gathered} a_0=f\left(x_0\right) \\ a_1=f\left[x_0, x_1\right] \\ a_2=f\left[x_0, x_1, x_2\right] \end{gathered} \]

其中一般式:

\[a_k=f\left[x_0, x_1, \cdots, x_k\right] \quad(k=0,1, \cdots, n) \]

将求得系数代入多项式中即可得到n次牛顿插值公式

\[N_n(x)=f\left(x_0\right)+f\left[x_0, x_1\right]\left(x-x_0\right)+\cdots+f\left[x_0, x_1, \cdots, x_n\right]\left(x-x_0\right)\left(x-x_1\right) \cdots\left(x-x_n\right) \]

其余项为

\[\begin{gathered} R_n(x)=f\left[x_0, x_1, \cdots, x_n, x\right]\left(x-x_0\right)\left(x-x_1\right) \cdots\left(x-x_n\right) \\ =f\left[x_0, x_1, \cdots, x_n, x\right] \prod_{i=0}^n\left(x-x_i\right)=\frac{f^{(n+1)}(\xi)}{(n+1) !} \prod_{i=0}^n\left(x-x_i\right) \\ f\left[x_0, x_1, \cdots, x_n\right]=\frac{f^{(n+1)}(\xi)}{(n+1) !} \end{gathered} \]

特点

\((1)\) 添加新节点计算量小

\((2)\) 仍会出现龙格现象

http://t.csdnimg.cn/ksi89

牛顿插值的几何解释是怎么样的？ - 马同学的回答 - 知乎
https://www.zhihu.com/question/22320408/answer/141973314

1.4 三次样条插值

样条插值

即取插值函数为样条函数, 称为样条插值。例如分段线性插值是一次样条插值。

三次样条插值的数学表示

定义

即已知函数 \(y=f(x)\) 在区间 \([a, b]\) 上的 \(n+1\) 个节点

\[a=x_0<x_1<\cdots<x_{n-1}<x_n=b \]

上的值 \(y_i=f\left(x_i\right)(i=0,1, \cdots, n)\), 求插值函数 \(S(x)\), 使得
\((1)\) \(S\left(x_i\right)=y_i(i=0,1, \cdots, n)\) 。
\((2)\) 在每个小区间 \(\left[x_i, x_{i+1}\right](i=0,1, \cdots, n-1)\) 上 \(S(x)\) 是三次多项式, 记为 \(S_i(x)\) 。
\((3)\) \(S(x)\) 在 \([a, b]\) 上二阶连续可微。

函数 \(S(x)\) 称为 \(f(x)\) 的三次样条插值函数。

固定条件

由条件(2), 不妨记

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

式中: \(a_i, b_i, c_i, d_i\) 为待定系数, 共 \(4 n\) 个。
由条件 (3), 有

\[\left\{\begin{array}{l} S_i\left(x_{i+1}\right)=S_{i+1}\left(x_{i+1}\right), \\ S_i^{\prime}\left(x_{i+1}\right)=S_{i+1}^{\prime}\left(x_{i+1}\right), i=0,1, \cdots, n-2 \\ S^{\prime \prime}{ }_i\left(x_{i+1}\right)=S^{\prime \prime}{ }_{i+1}\left(x_{i+1}\right), \end{array}\right. \]

容易看出, 式 (5.1) 和式 (5.2) 共含有 \(4 n-2\) 个方程, 为确定 \(S(x)\) 的 \(4 n\) 个待定参数, 尚需再给出两个边界条件。

边界条件

常用的三次样条函数的边界条件有 3 种类型:
\((1)\) \(S^{\prime}(a)=y_0^{\prime}, S^{\prime}(b)=y_n^{\prime}\) 。 由这种边界条件建立的样条插值函数称为 \(f(x)\) 的完备三次样条插值函数。
特别地, \(y^{\prime}{ }_0=y_n^{\prime}=0\) 时, 样条曲线在端点处呈水平状态。
如果 \(f^{\prime}(x)\) 不知道, 我们可以要求 \(S^{\prime}(x)\) 与 \(f^{\prime}(x)\) 在端点处近似相等。这时以 \(x_0, x_1\), \(x_2, x_3\) 为节点作一个三次 Newton 插值多项式 \(N_a(x)\), 以 \(x_n, x_{n-1}, x_{n-2}, x_{n-3}\) 作一个三次 Newton 插值多项式 \(N_b(x)\),要求

\[S^{\prime}(a)=N_a^{\prime}(a), S^{\prime}(b)=N_b^{\prime}(b) 。 \]

由这种边界条件建立的三次样条称为 \(f(x)\) 的 Lagrange 三次样条插值函数。
\((2)\) \(S^{\prime \prime}(a)=y_0^{\prime \prime}, S^{\prime \prime}(b)=y_n^{\prime \prime}\) 。特别地, \(y^{\prime \prime}{ }_0=y_n^{\prime \prime}=0\) 时, 称为自然边界条件。
\((3)\) \(S^{\prime}(a+0)=S^{\prime}(b-0), S^{\prime \prime}(a+0)=S^{\prime \prime}(b-0)\), 此条件称为周期条件。

2 拟合

曲线拟合问题的提法是, 已知一组 (二维) 数据, 即平面上的 \(n\) 个点 \(\left(x_i, y_i\right), i=1\), \(2, \cdots, n, x_i\) 互不相同, 寻求一个函数 (曲线) \(y=f(x)\), 使 \(f(x)\) 在某种准则下与所有数据点最为接近,即曲线拟合得最好。

2.1 线性最小二乘法

线性最小二乘法是解决曲线拟合最常用的方法。

基本思路

令

\[f(x)=a_1 r_1(x)+a_2 r_2(x)+\cdots+a_m r_m(x) \]

其中:

\(r_k(x)\) 为事先选定的一组线性无关的函数; \(a_k\) 为待定系数 \((k=1,2, \cdots, m ; m<n)\) 。
拟合准则是使 \(y_i(i=1,2, \cdots, n)\) 与 \(f\left(x_i\right)\) 的距离 \(\delta_i\) 的平方和最小, 称为最小二乘准则。

步骤

1. 函数 \(r_k(x)\) 的选取

面对一组数据 \(\left(x_i, y_i\right), i=1,2, \cdots, n\), 用线性最小二乘法作曲线拟合时, 首要的也是关键的一步是恰当地选取 \(r_1(x), \cdots, r_m(x)\) 。

如果通过机理分析, 能够知道 \(y\) 与 \(x\) 之间的函数关系, 则 \(r_1(x), \cdots, r_m(x)\) 容易确定。

若无法知道 \(y\) 与 \(x\) 之间的关系, 通常可以将数据 \(\left(x_i, y_i\right), i=1,2, \cdots, n\) 作图, 直观地判断应该用什么样的曲线去作拟合。

常用的曲线有：
(1) 直线 \(y=a_1 x+a_2\) 。
(2) 多项式 \(y=a_1 x^m+\cdots+a_m x+a_{m+1}\) (一般 \(m=2, 3\), 不宜太高 \() 。\)
(3) 双曲线 (一支) \(y=\frac{a_1}{x}+a_2\) 。
(4) 指数曲线 \(y=a_1 \mathrm{e}^{a_2 x}\) 。
对于指数曲线, 拟合前需作变量代换, 化为对 \(a_1, a_2\) 的线性函数。

已知一组数据, 用什么样的曲线拟合最好, 可以在直观判断的基础上, 选几种曲线分别拟合,然后比较,看哪条曲线的最小二乘指标 \(J\) 最小。

2. 系数 \(a_k\) 的确定

记

\[J\left(a_1, \cdots, a_m\right)=\|\boldsymbol{R A}-\boldsymbol{Y}\|_{2}^2=\sum_{i=1}^n \delta_i^2=\sum_{i=1}^n\left[f\left(x_i\right)-y_i\right]^2, \]

为求 \(a_1, \cdots, a_m\) 使 \(J\) 达到最小, 只需利用极值的必要条件 \(\frac{\partial J}{\partial a_j}=0(j=1, \cdots, m)\), 得到关于 \(a_1, \cdots, a_m\) 的线性方程组

\[\sum_{i=1}^n r_j\left(x_i\right)\left[\sum_{k=1}^m a_k r_k\left(x_i\right)-y_i\right]=0, j=1, \cdots, m, \]

即

\[\sum_{k=1}^m a_k\left[\sum_{i=1}^n r_j\left(x_i\right) r_k\left(x_i\right)\right]=\sum_{i=1}^n r_j\left(x_i\right) y_i, j=1, \cdots, m, \]

记

\[\begin{aligned} \boldsymbol{R} & =\left[\begin{array}{ccc} r_1\left(x_1\right) & \cdots & r_m\left(x_1\right) \\ \vdots & \vdots & \vdots \\ r_1\left(x_n\right) & \cdots & r_m\left(x_n\right) \end{array}\right]_{n \times m}, \\ \boldsymbol{A} & =\left[a_1, \cdots, a_m\right]^{\mathrm{T}}, \boldsymbol{Y}=\left[y_1, \cdots, y_n\right]^{\mathrm{T}}, \end{aligned} \]

方程组式 (5.5) 可表为

\[\boldsymbol{R}^{\mathrm{T}} \boldsymbol{R} \boldsymbol{A}=\boldsymbol{R}^{\mathrm{T}} \boldsymbol{Y} \]

当 \(\left\{r_1(x), \cdots, r_m(x)\right\}\) 线性无关时, \(\boldsymbol{R}\) 列满秩, \(\boldsymbol{R}^{\mathrm{T}} \boldsymbol{R}\) 可逆, 于是方程组式 (5.6) 有唯一解

\[\boldsymbol{A}=\left(\boldsymbol{R}^{\mathrm{T}} \boldsymbol{R}\right)^{-1} \boldsymbol{R}^{\mathrm{T}} \boldsymbol{Y} \]

3. Matlab解法

Matlab 中的线性最小二乘的标准型为

\[\min _A\|\boldsymbol{R A}-\boldsymbol{Y}\|_2^2 \]

命令为 \(\boldsymbol{A}=\boldsymbol{R} \backslash \boldsymbol{Y}\)

2.2 最小二乘优化

优化问题中目标函数由若干函数的平方和组成且求其最小，则属于最小二乘优化问题

可以用Matlab优化工具箱解决

2.3 曲线拟合与函数逼近

拟合： 离散\(\Rightarrow\) 连续函数，符合最小二乘准则

逼近： 复杂连续函数\(\Rightarrow\) 简单连续函数，符合最小平方逼近准则

最小平方逼近准则原理类似最小二乘准则，可看作连续的版本