【数模学习笔记】插值算法和拟合算法

声明：以下笔记中的图片以及内容
均整理自“数学建模学习交流”清风老师的课程资料，仅用作学习交流使用

插值算法

实际上本栏重点只有三次Hermite插值和三次样条插值的两行调用代码，其他的全是废（原）话（理）

定义

设函数$y=f(x)$在区间$[a,b]$上有定义，且已知在点
$a\le x_0<x_1<\cdots <x_n\le b$
上的值分别为：$y_0,y_1,\cdots ,y_n$，
若存在一简单函数$P（x）$，使
$$P(x_i)=y_i(i=0,1,2\cdots ,n)$$
则称$P(x)$为$f(x)$的插值函数，点$x_0,x_1,\cdots ,x_n$称为插值节点，包含插值节点的区间$[a,b]$称为插值区间，求插值函数$P(x)$的方法称为插值法。

三个类型插值举例

插值多项式

• 示例：已知函数$y=f(x)$在点$x_0=0$，$x_1=1$，$x_2=2$上的值分别为$f(0)=1$，$f(1)=2$，$f(2)=5$。我们可以构造一个二次插值多项式$P(x)$来逼近$f(x)$。
  • 设$P(x)=a_0+a_{1}x+a_2{x}^2$，将已知点代入可得方程组：
    • $\{\begin{array}{l}{a}_0+a_1\times 0+a_2\times 0^2=1\\ a_0+a_1\times 1+a_2\times 1^2=2\\ a_0+a_1\times 2+a_2\times 2^2=5\end{array}$
  • 解这个方程组可得$a_0=1$，$a_1=0$，$a_2=1$，所以插值多项式$P(x)=1+x^2$。在区间$[0,2]$上，这个多项式可以用来近似原函数$f(x)$。

分段插值

• 示例：假设我们要对函数$f(x)=\frac{1}{1+25x^2}$在区间$[-1,1]$上进行插值。如果我们只用一个高次多项式进行插值，可能会出现龙格现象（即在区间端点附近出现剧烈振荡）。这时可以采用分段插值，比如将区间$[-1,1]$等分成若干个子区间，如$[-1,-0.5],[-0.5,0],[0,0.5],[0.5,1]$。在每个子区间上分别进行低次多项式插值（如线性插值或二次插值）。
  • 在区间$[-1,-0.5]$上，已知端点值$f(-1)$和$f(-0.5)$，通过线性插值得到该区间上的插值函数$P_1(x)$，使得$P_1(-1)=f(-1)$，$P_1(-0.5)=f(-0.5)$，并且在该区间内$P_1(x)$近似$f(x)$。
  • 同理，在其他子区间上也进行类似的操作，得到相应的插值函数$P_2(x)$、$P_3(x)$、$P_4(x)$。这样就构成了整个区间$[-1,1]$上的分段插值函数。

三角插值

• 示例：对于一个周期为$2\pi$的函数$f(x)$，我们在区间$[0,2\pi ]$上取$n$个等距节点$x_k=\frac{2k\pi }{n}$，$k=0,1,\cdots ,n-1$，已知这些节点上的函数值$f(x_k)$。
  • 三角插值多项式可以表示为$S_n(x)=\frac{a_0}{2}+{\sum }_{k=1}^n(a_k\cos kx+b_k\sin kx)$，其中系数$a_k$和$b_k$通过以下公式计算：
    • $a_k=\frac{1}{\pi }{\int }_0^{2\pi }f(x)\cos kxdx$（$k=0,1,\cdots ,n$）
    • $b_k=\frac{1}{\pi }{\int }_0^{2\pi }f(x)\sin kxdx$（$k=1,2,\cdots ,n$）
  • 例如，假设$f(x)$是一个简单的周期函数，在$[0,2\pi ]$上$f(x)=\{\begin{array}{ll}x,& 0\le x<\pi \\ 2\pi -x,& \pi \le x<2\pi \end{array}$，通过计算上述积分可以得到相应的三角插值多项式$S_n(x)$，用它来近似原函数$f(x)$在整个周期上的行为。

这些插值方法在不同的场景下各有优势，分段插值可以避免高次多项式的振荡问题，插值多项式计算相对简单，三角插值适用于周期函数等情况，它们都为函数的近似和数据的处理提供了有效的手段。

一般插值多项式

原理

• 定理：
  • 设有$n+1$个互不相同的节点$(x_i,y_i)$ $(i=0,1,2,\cdots ,n)$
  • 则存在唯一的多项式：
    $L_n(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$
  • 使得 $L_n(x_j)=y_j,(j=0,1,2,\cdots ,n)$ 成立
    （即就是仅存在唯一多项式使其过给定的$n+1$个点）
• 证明：（利用范德蒙行列式）
  • 构造方程组
    $\{\begin{array}{l}{a}_0+a_1{x}_0+a_2{x}_0^2+\cdots +a_n{x}_0^n=y_0\\ a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n=y_1\\ \cdots \cdots \cdots \\ a_0+a_1{x}_n+a_2{x}_n^2+\cdots +a_n{x}_n^n=y_n\end{array}$
  • 令：
    $A=\left[\begin{array}{cccc}1& x_0& \cdots & x_0^n\\ 1& x_1& \cdots & x_1^n\\ ⋮& ⋮& \cdots & ⋮\\ 1& x_n& \cdots & x_n^n\end{array}\right]$
    $X={\left[\begin{array}{c}{a}_0{a}_1\cdots a_n\end{array}\right]}^T$
    $Y={\left[\begin{array}{c}{y}_0{y}_1\cdots y_n\end{array}\right]}^T$
  • 方程组的矩阵形式：$AX=Y$
  • 由于$|A|={\prod }_{i=1}^n{\prod }_{j=0}^{i-1}(x_i-x_j)\ne 0$，所以方程组有唯一解，从而$L_n(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$唯一存在.
• 注：
  • 只要$n+1$个节点互异，满足上述插值条件的多项式是唯一存在的。
  • 如果不限制多项式的次数，插值多项式并不唯一。

拉格朗日插值法

在数值分析中，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。

• 两个点:$(x_0,y_0),(x_1,y_1)$
  $f(x)=\frac{x-x_1}{x_0-x_1}{y}_0+\frac{x-x_0}{x_1-x_0}{y}_1$

• 三个点:$(x_0,y_0),(x_1,y_1),(x_2,y_2)$
  $f(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}{y}_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}{y}_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}{y}_2$

• 拉格朗日插值多项式一般形式
  记 ${\omega }_{n+1}(x)=(x-x_0)(x-x_1)\cdots (x-x_n)$
  则 ${\omega }_{n+1}^{\prime }(x_k)=(x_k-x_0)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n).$
  记 $l_i(x)=\frac{(x-x_0)\cdots (x-x_{i-1})(x-x_{i+1})\cdots (x-x_n)}{(x_i-x_0)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}$
  $=\frac{(x-x_0)\cdots (x-x_{i-1})(x-x_i)(x-x_{i+1})\cdots (x-x_n)}{(x-x_i)(x_i-x_0)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}$
  $=\frac{{\omega }_{n+1}(x)}{(x-x_i){\omega }_{n+1}^{\prime }(x_i)}$
  则拉格朗日插值多项式为 $L_n(x)={\sum }_{i=0}^n{y}_i{l}_i(x)$
  即 $L_n(x)={\sum }_{k=0}^n{y}_k\frac{{\omega }_{n+1}(x)}{(x-x_k){\omega }_{n+1}^{\prime }(x_k)}.$

但实际上，拉格朗日插值法实际情况下并不常用，因为有一个很大的问题——龙格现象。

龙格现象

龙格现象是指在使用多项式插值逼近函数时，当插值节点等距分布且插值多项式的次数较高时，在插值区间的两端会出现剧烈振荡的现象，导致插值结果在区间端点附近与原函数偏差较大，不能很好地逼近原函数。
所以在不清楚曲线具体运动趋势的情况下，不要轻易使用高次插值。

分段线性插值

为了避免龙格现象带来的影响，故我们使用分段低次插值。
这里以分段二次插值为例：
选取跟节点$x$最近的三个节点$x_{i-1},x_i,x_{i+1}$进行二次插值。
即在每一个区间$[x_{i-1},x_{i+1}]$上，取：
$$f(x)\approx L_2(x)=\sum _{k=i-1}^{i+1}[y_k\prod _{\begin{array}{c}j=i-1\\ j\ne k\end{array}}^{i+1}\frac{(x-x_j)}{(x_k-x_j)}]$$

这种分段的低次插值称为分段二次插值，在几何上就是用分段抛物线代替$y=f(x)$，故分段二次插值又称为分段抛物插值。

牛顿插值法

牛顿插值法每次插值只和前n项的值有关，这样每次只要在原来的函数上添加新的项，就能够产生新的函数。

$f(x)=f(x_0)+f[x_0,x_1](x-x_0)$
$+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdots$
$+f[x_0,x_1,\cdots ,x_{n-2},x_{n-1}](x-x_0)(x-x_1)\cdots (x-x_{n-3})(x-x_{n-2})$
$+f[x_0,x_1,\cdots ,x_{n-1},x_n](x-x_0)(x-x_1)\cdots (x-x_{n-2})(x-x_{n-1})$

• 差商定义:
  称$f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}$为函数$f(x)$关于点$x_0,x_k$的一阶差商（亦称均差）.

• 二阶差商:
  $f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}$

• K阶差商:（$x_k,x_{k-1}$可以不相邻）
  $f[x_0,x_1,\cdots ,x_k]=\frac{f[x_1,\cdots ,x_{k-1},x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_0}$

与拉格朗日插值法相比，牛顿插值法的计算过程具有继承性。但是牛顿插值也存在龙格现象的问题。同时，两种插值仅仅要求插值多项式在插值节点处与被插函数有相等的函数值，而这种插值多项式却不能全面反映被插值函数的性态。然而在许多实际问题中，不仅要求插值函数与被插值函数在所有节点处有相同的函数值，它也需要在一个或全部节点上插值多项式与被插函数有相同的低阶甚至高阶的导数值。

Hermite埃尔米特插值

原理

• 设函数$f(x)$在区间$[a,b]$上有$n+1$个互异节点 $a=x_0<x_1<x_2<\cdots <x_n=b$，
• 定义在$[a,b]$上函数$f(x)$在节点上满足：
  $f(x_i)=y_i,f^{\prime }(x_i)=y_i^{\prime }(i=0,1,2,\cdots ,n)$ （$2n+2$个条件）
• 可唯一确定一个次数不超过$2n+1$的多项式 $H_{2n+1}(x)=H(x)$满足：
  $H(x_j)=y_j,H^{\prime }(x_j)=m_j(j=0,1,\cdots ,n)$.
• 其余项为：
  $R(x)=f(x)-H(x)=\frac{f^{(2n+2)}(\xi )}{(2n+2)!}{\omega }_{2n+2}(x)$

可以看到，直接使用Hermite插值得到的次数较高，也存在龙格现象，因此在实际应用中，往往使用分段三次Hermite插值多项式(pchip).

分段三次埃尔米特插值

构造

• 已知给定$n+1$个节点$x_0,x_1,\cdots ,x_n$，以及对应的函数值$y_0,y_1,\cdots ,y_n$和导数值$m_0,m_1,\cdots ,m_n$。

• 函数值基函数$h_i(x)$：

  • 对于节点$x_i$，构造函数值基函数$h_i(x)$，使其满足$h_i(x_j)={\delta }_{ij}$（其中${\delta }_{ij}$是克罗内克符号，当$i=j$时，${\delta }_{ij}=1$；当$i\ne j$时，${\delta }_{ij}=0$）。
  • 一种常见的构造方式是：
    $h_i(x)=\left(1+2\frac{x-x_i}{x_{i+1}-x_i}\right){\left(\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)}^2$，当$i=0,1,\cdots ,n-1$；
    $h_n(x)=\left(1+2\frac{x-x_n}{x_{n-1}-x_n}\right){\left(\frac{x-x_{n-1}}{x_n-x_{n-1}}\right)}^2$。

• 导数值基函数${\hat{h}}_i(x)$：

  • 构造导数值基函数${\hat{h}}_i(x)$，使其满足${\hat{h}}_i^{\prime }(x_j)={\delta }_{ij}$。
  • 例如：
    ${\hat{h}}_i(x)=(x-x_i){\left(\frac{x-x_{i+1}}{x_i-x_{i+1}}\right)}^2$，当$i=0,1,\cdots ,n-1$；
    ${\hat{h}}_n(x)=(x-x_n){\left(\frac{x-x_{n-1}}{x_n-x_{n-1}}\right)}^2$。

• 构造三次埃尔米特插值多项式$H(x)$
  $H(x)={\sum }_{i=0}^n[y_i{h}_i(x)+m_i{\hat{h}}_i(x)]$

应用

Matlab有内置的函数（直接调用就行了）：
p = pchip(x,y,new_x)
x,y是已知的样本点的坐标，new_x是插值横坐标

x = -pi:pi; y = sin(x);
new_x = -pi:0.1:pi;
p = pchip(x,y,new_x);
plot(x, y, 'o', new_x, p, 'r-')

三次样条插值

定义

设$y=f(x)$在点$x_0,x_1,x_2,\cdots x_n$的值为$y_0,y_1,y_2,\cdots y_n$，若函数$S(x)$满足下列条件

• $S(x_i)=f(x_i)=y_i$，$i=0,1,2,\cdots ,n$
• 在每个子区间$[x_i,x_{i+1}](i=0,1,2,\cdots ,n-1)$上$S(x)$是三次多项式
• $S(x)$在$[a,b]$上二阶连续可微。

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

• $S(x)$除了满足基本插值条件$s(x_i)=f_i$外还应具有如下形式：
  $S(x)=\{\begin{array}{ll}{S}_0(x),& x\in [x_0,x_1],\\ S_1(x),& x\in [x_1,x_2],\\ ⋮& \\ S_{n-1}(x),& x\in [x_{n-1},x_n];\end{array}$
  $S_i(x)\in C^3([x_i,x_{i+1}])$.
  并且满足条件：
  $\{\begin{array}{l}{S}_{i-1}(x_i)=S_i(x_i),\\ S_{i-1}^{\prime }(x_i)=S_i^{\prime }(x_i),\\ S_{i-1}^{\prime \prime }(x_i)=S_i^{\prime \prime }(x_i),\end{array}$

应用

Matlab有内置的函数：
p = spline (x,y, new_x)
x,y是已知的样本点的坐标，new_x是要插入处对应的横坐标

x = -pi:pi;
y = sin(x);
new_x = -pi:0.1:pi;
p1 = pchip(x,y,new_x);    %分段三次埃尔米特插值
p2 = spline(x,y,new_x);   %三次样条插值
plot(x,y,'o',new_x,p1,'r-',new_x,p2,'b-')
Legend('样本点','三次埃尔米特插值','三次样条插值','Location','SouthEast')    %标注显示在东南方向

三次Hermite插值和三次样条插值的对比

三次样条比Hermite更光滑，实际应用中两种都可以使用。

n维数据的插值

这个了解即可。
p = interpn(x1,x2,...,xn, y, new_x1,new_x2,...,new_xn, method)
x1,x2,...,xn是已知的样本点的横坐标
y是已知的样本点的纵坐标
new_x1,new_x2,...,new_xn是要插入点的横坐标
method是要插值的方法
‘linear’: 线性插值 (默认算法);
‘cubic’: 三次插值;
‘spline’: 三次样条插值法; (最为精准)
‘nearest’: 最邻近插值算法。

x = -pi:pi; y = sin(x);
new_x = -pi:0.1:pi;
p = interpn (x, y, new_x,'spline');
%等价于p = spline(x, y, new_x);
plot(x, y, 'o', new_x, p, 'r-')

---

省流：本栏最有用的两行
p = pchip(x,y,new_x)
p = spline (x,y, new_x)

---

拟合算法

插值算法中，得到的多项式$f(x)$要经过所有样本点。但是如果样本点太多，那么这个多项式次数过高，会造成龙格现象。
比如说这个例子

尽管我们可以选择分段的方法避免龙格现象，但是更多时候我们更倾向于得到一个确定的曲线，尽管这条曲线不能经过每一个样本点，但只要保证误差足够小即可，这就是拟合的思想。（拟合的结果是得到一个确定的较为简单的曲线）

最小二乘法

高中学的 懒得解释。


解释一下$$k,b=\underset{\begin{array}{c}k,b\end{array}}{\mathrm{argmin}}(f(x))$$指的是求使得$f(x)$取最小值时的参数$k,b$的值

• 推导过程：
  
• 结论：

$$\hat{k}=\frac{n{\sum }_{i=1}^n{x}_i{y}_i-{\sum }_{i=1}^n{y}_i{\sum }_{i=1}^n{x}_i}{n{\sum }_{i=1}^n{x}_i^2-{\sum }_{i=1}^n{x}_i{\sum }_{i=1}^n{x}_i}$$
$$\hat{b}=\frac{{\sum }_{i=1}^n{x}_i^2{\sum }_{i=1}^n{y}_i-{\sum }_{i=1}^n{x}_i{\sum }_{i=1}^n{x}_i{y}_i}{n{\sum }_{i=1}^n{x}_i^2-{\sum }_{i=1}^n{x}_i{\sum }_{i=1}^n{x}_i}$$

拟合评价

• 三个基本概念

  • 总体平方和 $SST$ (Total sum of squares)：$SST={\sum }_{i=1}^n(y_i-\bar{y}{)}^2$

  • 误差平方和 $SSE$ (The sum of squares due to error)：$SSE={\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2$

  • 回归平方和 $SSR$ (Sum of squares of the regression)：$SSR={\sum }_{i=1}^n({\hat{y}}_i-\bar{y}{)}^2$

• 可以证明：$SST=SSE+SSR$

  • 证明：

• 拟合优度：$0\le R^2=\frac{SSR}{SST}=\frac{SST-SSE}{SST}=1-\frac{SSE}{SST}\le 1$
  $R^2$越接近$1$，说明误差平方和 $SSE$ 越接近$0$，说明拟合的越好。

• 注意：$R^2$只能用于拟合函数是线性函数时，拟合结果的评价.
  线性函数和其他函数（如复杂的指数函数）比较拟合的好坏，直接看$SSE$即可

  • 线性函数的界定：这里我们说的线性函数指的是对参数为线性，即在函数中，参数仅以一次方出现，并且不能乘以或除以其他任何参数，也不能出现参数的复合函数形式。
    形如 $y=a+b{x}^2$ 是线性函数，这里的参数就是$a,b$.
    形如 $y=a+b^{3}x$、$y=a+bx+bc{x}^2$、$y=a(x-b{)}^2$、$y=a\sin (b+cx)$都不是线性函数，不能用$R^2$.

实际应用中，常使用多个不同的拟合函数去试，然后比较这几个拟合函数的拟合优度，取最优的一个拟合函数。
曲线复杂度与拟合误差之间的取舍：次数越高（或者指数函数越复杂），误差越小，但是产生龙格现象的概率越大，且函数曲线越复杂。（一个理解：在插值算法中，曲线经过每一个点，那么曲线的$SSE$即为0，如果放在拟合算法中来讲，这种曲线的拟合优度极高，但是曲线本身十分复杂，这就违背了拟合的初衷）
因此实际上，如果两个拟合函数的$SSE$相差不大，但是曲线复杂度相差较大，那么就取明显更为简单的那个函数，哪怕它的$SSE$相较另一个更大一点。

y\_hat = k*x+b;  % y的拟合值
SSR = sum((y\_hat-mean(y)).^2)  % 回归平方和 注:mean()是求均值的函数。
SSE = sum((y\_hat-y).^2)  % 误差平方和
SST = sum((y-mean(y)).^2)  % 总体平方和
SST-SSE-SSR
R\_2 = SSR / SST

Matlab自带拟合工具箱cftool

• 拟合方法

  • Custom Equations：用户自定义的函数类型 注意修改范围Fit options
  • Exponential：指数逼近，有2种类型， $a\ast exp(b\ast x)$ 、 $a\ast exp(b\ast x)+c\ast exp(d\ast x)$
  • Fourier：傅立叶逼近，有7种类型，基础型是 $a0+a1\ast cos(x\ast w)+b1\ast sin(x\ast w)$
  • Gaussian：高斯逼近，有8种类型，基础型是 $a1\ast exp(-((x-b1)/c1{)}^2)$
  • Interpolant：插值逼近，有4种类型，linear、nearest neighbor、cubic spline、shape-preserving
  • Polynomial：多形式逼近，有9种类型，linear ~、quadratic ~、cubic ~、4-9th degree ~
  • Power：幂逼近，有2种类型，$a\ast x^b、a\ast x^b+c$
  • Rational：有理数逼近，分子、分母共有的类型是linear ~、quadratic ~、cubic ~、4-5th degree ~；此外，分子还包括constant型
  • Smoothing Spline：平滑逼近
  • Sum of Sin Functions：正弦曲线逼近，有8种类型，基础型是 $a1\ast sin(b1\ast x+c1)$
  • Weibull：只有一种，$a\ast b\ast x^{(b-1)}\ast exp(-a\ast x^b)$

• 拟合工具还可以直接生成代码Generate code，可放在附录里。

• 小tips：导出图像可以在文件-导出设置-渲染-修改分辨率，使图更清晰