数值方法笔记4：插值、近似和拟合

多项式方程如$f(x)=1+x;f(x)=x^2+x^4$
非多项式方程如$f(x)=e^{x}f(x)=\sin (x)\ln (x)$

多项式方程对于计算机好计算，非多项式方程则不然，可以用多项式方程来计算非多项式方程吗？怎么用一个多项式方程来表示一个非多项式方程，是我们要解决的问题。

1. 插值

1.1 插值的一些概念

1.1.1 插值的定义

插值：给定一组点$(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$,构造一个函数$y=f(x)$,满足$y_i=f(x_i),i=1,2,...,n$

例子：

1.1.2 插值的存在性

设多项式方程为：
$$\begin{gathered} P_n(x)=a_0+a_{1}x+a_2{x}^2+\dots \dots +a_n{x}^{\mathrm{n}} \\ {P}_n\left(x_i\right)=y_i,\mathrm{i}=0,1,2,\dots ,\mathrm{n},x_i\ne x_{\mathrm{j}} \end{gathered}$$

有：

$$\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^n\\ 1& x_1& x_1^2& \cdots & x_1^n\\ 1& x_2& x_2^2& \cdots & x_2^n\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& x_n& x_n^2& \cdots & x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_n\end{array}\right]$$
也即$\mathbf{Aa}=\mathbf{y}$的形式。

要求$\mathbf{A}$行列式不等于0，根据vandermonde行列式，我们知道行列式：

$$\left|\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^n\\ 1& x_1& x_1^2& \cdots & x_1^n\\ 1& x_2& x_2^2& \cdots & x_2^n\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& x_n& x_n^2& \cdots & x_n^n\end{array}\right|=\prod _{n\ge i>j\ge 0}\left(x_i-x_j\right)$$

1.1.3 插值的误差分析

一般来说，如果点$(x_i,y_i)$是任意选择的，我们无法进行误差分析。

但如果这些点是按顺序从一个已知的函数$f(x)$中抽取的，这意味着$(x_i,y)=(x_i,f(x_i))$，误差$R(x)=|f(x)-P_n(x)|$可以被预测到。

1.2 拉格朗日插值(Lagrange Interpolation)

给定一组点$(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$

$$\begin{array}{l}y=P_N(x)=y_0{L}_{N,0}(x)+y_1{L}_{N,1}(x)+\dots +y_N{L}_{N,N}(x)\\ ={\sum }_{k=0}^N{y}_k{L}_{N,k}(x)\end{array}$$

$$\begin{array}{rl}& L_{N,k}(x)=\frac{\left(x-x_0\right)\dots \left(x-x_{k-1}\right)\left(x-x_{k+1}\right)\dots \left(x-x_N\right)}{\left(x_k-x_0\right)\dots \left(x_k-x_{k-1}\right)\left(x_k-x_{k+1}\right)\dots \left(x_k-x_N\right)}=\frac{{\prod }_{\begin{array}{c}j=0j\ne k\end{array}}^N\left(x-x_j\right)}{{\prod }_{\begin{array}{c}j=0j\ne k\end{array}}^N\left(x_k-x_j\right)}\\ & L_{N,k}\left(x_0\right)=0,\dots ,L_{N,k}\left(x_1\right)=0,L_{N,k}\left(x_k\right)=1,L_{N,k}\left(x_{k+1}\right)=0,\dots ,L_{N,k}\left(x_N\right)=0\end{array}$$

1.2.1 拉格朗日插值误差分析

定理：假设$f(x)$在$[a,b]$上是$(n+1)$阶可微，$L_n(x)$是拉格朗日函数，那么误差估计如下

$$R_n(x)=f(x)-P_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}\prod _{k=0}^n\left(x-x_k\right),\xi \in \left[x_0,x_n\right]$$

证明：
$$\begin{array}{rl}设& w(x)=\prod _{k=0}^n\left(x-x_k\right)\\ 令& R(t)=f(t)-\frac{f(x)-P_n(x)}{w(x)}w(t)-P_n(t)\\ 有& R\left(x_i\right)=f\left(x_i\right)-\frac{f(x)-P_n(x)}{w(x)}w\left(x_i\right)-P_n\left(x_i\right)=0\\ & R(x)=f(x)-\frac{f(x)-P_n(x)}{w(x)}w(x)-P_n(x)=0\end{array}$$

于是$R(t)$有$n+2$个根：$x_0,x_1,\cdots ,x_n,x$，使用罗尔定理，$R$的$n+1$阶导数$R^{(n+1)}(x)$至少有一个根，使得$R^{(n+1)}(\xi )=0\xi \in [x_0,x_n]$,于是：

$$\begin{array}{rl}& R^{(n+1)}(t)=f^{(n+1)}(t)-\frac{f(x)-P_n(x)}{w(x)}{w}^{(n+1)}(t)-P_n^{(n+1)}(t)\\ & =f^{(n+1)}(t)-\frac{f(x)-P_n(x)}{w(x)}(n+1)!\\ & R^{(n+1)}(\xi )=f^{(n+1)}(\xi )-\frac{f(x)-P_n(x)}{w(x)}(n+1)!=0\\ & f(x)-P_n(x)=f^{(n+1)}(\xi )w(x)/(n+1)!\end{array}$$

可以导出的其他性质：

$$\begin{array}{rl}& w(x)=\prod _{k=0}^n\left(x-x_k\right)\\ & w^{\prime }(x)=\left(x-x_1\right)\left(x-x_2\right)\dots \left(x-x_n\right)+\left(x-x_0\right)\left(x-x_2\right)\dots +\dots \\ & L_{n,k}(x)=\left(x-x_0\right)\dots \left(x-x_{k-1}\right)\left(x-x_{k+1}\right)\dots \left(x-x_n\right)/w^{\prime }\left(x_k\right)\end{array}$$

于是有：

$$\sum _{k=0}^n{L}_{n,k}(x)=\sum _{k=0}^n\left(x-x_0\right)\dots \left(x-x_{k-1}\right)\left(x-x_{k+1}\right)\dots \left(x-x_n\right)/w^{\prime }\left(x_k\right)\equiv 1$$

证明： $$\text{ 令}f(x)=1,给定一系列点\left(x_{\mathrm{i}},1\right)\text{, i.e. }{y}_{\mathrm{i}}=1,P_n(x)=\sum _{i=0}^n{L}_{n,i}(x),R(x)\equiv 0$$

如何预测误差的边界

分别找到$f^{n+1}$和$w(x)$的上界有：

$$\begin{array}{rl}& \left|f^{(n+1)}(x)\right|<M_n\\ & |w(x)|=\left|\prod _{k=0}^n\left(x-x_k\right)\right|<W_{\mathrm{n}}\end{array}$$

于是有：

$$\left|R_n(x)\right|=\left|\frac{f^{(n+1)}(\xi )}{(n+1)!}\prod _{k=0}^n\left(x-x_k\right)\right|<\frac{M_n{W}_n}{(n+1)!}$$

如果$x_k=x_0+kh,$ 则$$W_1=\frac{h^2}{4},W_2=\frac{2h^3}{3\sqrt{3}},W_3=h^4$$

1.3 Newton多项式插值

给定一组点$(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$

我们观察这样一个递推序列：

$$P_N(x)=P_{N-1}(x)+a_N\left(x-x_0\right)\left(x-x_1\right)\dots \left(x-x_{N-1}\right)$$

怎么求其中的系数呢？例子：给两个点和三个点的插值如下：

我们把三点的插值的$a_2$计算一下，可以发现Newton的系数和微分很像！

于是我们可以给出多个点的递推公式：
$$\begin{array}{l}{P}_N(x)=a_0+a_1\left(x-x_0\right)+a_1\left(x-x_0\right)\left(x-x_1\right)+a_1\left(x-x_0\right)\left(x-x_1\right)\left(x-x_2\right)+\\ \dots +a_N\left(x-x_0\right)\left(x-x_1\right)\dots \left(x-x_{N-1}\right)\\ a_k=f\left[x_0,x_1,\dots ,x_k\right]\end{array}$$

其中

可以验算一下$P_2$

差分计算和微分也很相似，二次经过差分变一次。$P(x)$是$n$次多项式，$P^{\prime }(x)$是$n-1$次多项式。

性质：如果$f(x)$是$n$次多项式，那么$f[x,x_0,x_1,\cdots ,x_n]=0$

1.3.1 Newton多项式插值误差分析

如果$f(x)$是$n$次多项式，误差为：

$$R_n(x)=f(x)-P_n(x)=f\left[x,x_0,\dots ,x_n\right]\left(x-x_0\right)\left(x-x_1\right)\dots \left(x-x_n\right)$$

对于拉格朗日插值我们又知道：

$$R_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}\left(x-x_0\right)\left(x-x_1\right)\dots \left(x-x_n\right)$$

于是有：

$$f\left[x,x_0,\dots ,x_n\right]=\frac{f^{(n+1)}(\xi )}{(n+1)!}$$

当$n$的0，就变成中值定理。
于是我们知道拉格朗日插值和Newton插值可以互换。只是形式不同。

1.4 Chebyshev多项式确定插值点

根据前面的拉格朗日插值我们知道

我们怎么选择最佳的$(x_0,x_1,\dots ,x_N)$以最小化$|R(x)|$呢？

答案是，使用Chebyshev多项式的根。

Chebyshev多项式一般表达式：


注意$T_n(x)$最高次项的系数是$2^{n-1}$, 对于$T_n$有$n$个零点。

1.4.1 Chebyshev多项式性质

定理：
$$\begin{array}{rl}& \underset{-1\le x\le 1}{\max }\left|\frac{T_n(x)}{2^{n-1}}\right|\le \underset{-1\le x\le 1}{\max }|w(x)|\\ & \underset{-1\le x\le 1}{\max }\left|\frac{T_n(x)}{2^{n-1}}\right|=\frac{1}{2^{n-1}}\end{array}$$
使用反证法证明：
假设$$\underset{-1\le x\le 1}{\max }\left|\frac{T_n(x)}{2^{n-1}}\right|>\underset{-1\le x\le 1}{\max }|w(x)|$$
具体证明如下：
对于至少$n-1$次多项式(给定了$n+1$个点)，$w(x)-\frac{T_n(x)}{2^{n-1}}$有$n$个根，矛盾。
这是书里摘出的证明

$T_n(x)$最高次项的系数是$2^{n-1}$的证明。

正交性质

注意$x_k$ 是$T_{n+1}(x)=0$, 不是$T_n(x)=0$的根.

证明：

以下几个图的$f(x_k)$和$P_n(x_k)$是一个意思

1.5 有理插值

其中


例子：

1.5.1 有理插值误差分析

$$E(x)=f(x)-\frac{P(x)}{Q(x)}=\frac{\left(x-x_0\right)\dots \left(x-x_{m+n}\right)}{(m+n+1)!Q^2(x)}{\left[f(x)Q^2(x)\right]}_{x=\xi }^{(m+n+1)}$$

例子：

1.6 分段插值

单段拟合会出现Runge现象

1.6.1 分段线性插值

上图的$h=x_1-x_0$

1.6.2 分段Hermite插值

保证拟合曲线导数的光滑性。

归一化 Hermite Problem 形式

1.6.2.1 Hermite插值误差

对于$H_3$

Hermite赫米特多项式误差介于泰勒级数和拉格朗日之间

总结

拉格朗日插值和Newton插值形式不同，但是结果是一样的。可以互相转换。把拉格朗日插值和Newton插值和泰勒展式比较一下：

2. 近似

2.1 近似的一些概念

给出一个一般函数$y=f(x)$,求一个多项式函数$y=P(x)$,满足$|f(x)-P(x)|<\epsilon (x)$

例子：

2.2 帕德近似 (有理近似)

$$\begin{array}{rl}& R_{n,m}(x)=\frac{P_n(x)}{Q_m(x)}\\ & P_n(x)=p_0+p_{1}x+p_2{x}^2+\dots +p_n{x}^n\\ & Q_m(x)=1+q_{1}x+q_2{x}^2+\dots +q_m{x}^m\end{array}$$

共有$n+m$个未知数。

给定$n+m$的大小，当$n=m$或$n=m+1$近似误差最小。

例子：

连续形式：

3. 拟合

拟合：问题：给定一个曲线类型和一组点，在该类型的曲线中找出一条与这些点距离最小的曲线。

问题：如何评估误差？
有以下几种形式

3.1 最小二乘拟合

3.1.1 线性拟合

拟合类型：
$$y=f(x)=Ax+B$$

给定点：
$$\left(x_{\mathrm{i}},y_{\mathrm{i}}\right),i=1,2,\dots ,n$$

使用最小二乘法

$$\begin{array}{rl}d(A,B)& =\sum _{i=1}^n{\left[f\left(x_i\right)-y_i\right]}^2=\sum _{i=1}^n{\left[A{x}_i+B-y_i\right]}^2\\ & =\sum _{i=1}^n\left[A^2{x}_i^2+B^2+y_i^2+2AB{x}_i-2A{x}_i{y}_i-2B{y}_i\right]\\ & =\left(\sum _{i=1}^n{x}_i^2\right)A^2+n{B}^2+\left(2\sum _{i=1}^n{x}_i\right)AB-2\left(\sum _{i=1}^n{x}_i{y}_i\right)A-2\left(\sum _{i=1}^n{y}_i\right)B+\left(\sum _{i=1}^n{y}_i^2\right)\end{array}$$

$$\begin{array}{rl}\frac{\partial d(A,B)}{\partial A}& =2\left(\sum _{i=1}^n{x}_i^2\right)A+2\left(\sum _{i=1}^n{x}_i\right)B-2\sum _{i=1}^n{x}_i{y}_i=0\\ \frac{\partial d(A,B)}{\partial B}& =2nB+2\left(\sum _{i=1}^n{x}_i\right)A-2\sum _{i=1}^n{y}_i=0\end{array}$$

求解上面的线性方程组可以得到$A$和$B$。

3.1.2 指数拟合

拟合类型
$$y=f(x)=A{x}^M$$其中$M$是已知的。

给定点：
$$\left(x_{\mathrm{i}},y_{\mathrm{i}}\right),i=1,2,\dots ,n$$

应用最小二乘法：
$$\begin{array}{rl}& E(A)=\sum _{i=1}^n{\left(A{x}_i^M-y_i\right)}^2\to \min \\ & E^{\prime }(A)=2\sum _{i=1}^n\left[\left(A{x}_i^M-y_i\right)x_i^M\right]=2\sum _{i=1}^n\left(A{x}_i^{2M}-x_i^M{y}_i\right)\\ & =2\left(\sum _{i=1}^n{x}_i^{2M}\right)A-2\sum _{i=1}^n{x}_i^M{y}_i=0\\ & A=\frac{{\sum }_{i=1}^n{x}_i^M{y}_i}{{\sum }_{i=1}^n{x}_i^{2M}}\end{array}$$

另一种形式的指数拟合：
拟合形式：
$$y=C{e}^{Ax}$$

给定点：
$$\left(x_{\mathrm{i}},y_{\mathrm{i}}\right),i=1,2,\dots ,n$$

应用最小二乘法：
$$\begin{array}{rl}& E(A,C)=\sum _{i=1}^n{\left(C{e}^{A{x}_i}-y_i\right)}^2\to \min \\ & \frac{\partial E}{\partial A}=2\sum _{i=1}^n\left[\left(C{e}^{A{x}_i}-y_i\right)C{e}^{A{x}_i}{x}_i\right]=2\sum _{i=1}^n\left(C^2{e}^{2A{x}_i}{x}_i-C{e}^{A{x}_i}{x}_i{y}_i\right)=0\\ & \frac{\partial E}{\partial C}=2\sum _{i=1}^n\left[\left(C{e}^{A{x}_i}-y_i\right)e^{A{x}_i}\right]=2\sum _{i=1}^n\left(C{e}^{2A{x}_i}-e^{A{x}_i}{y}_i\right)=0\end{array}$$

这是非线性方程不好解得$A$和$C$。求解有两个思路：一是使用优化的方法比如单纯形法，二是使用Newton法迭代，另外可以考虑线性化：
$$\begin{array}{rl}& y=C{e}^{Ax}\\ & \ln (y)=Ax+\ln (C)\\ & Y=AX+B\end{array}$$

其中
$$Y=\ln (y),X=x,B=\ln (C)$$

3.1.3 其他形式拟合（线性化思想）

3.1.4 一般线性最小二乘法形式

3.1.5 非线性最小二乘法

在一般最小二乘的基础上添加权重。

由此导出：Levenberg-Marquardt方法

3.1.5.1 Levenberg-Marquardt方法

回顾我们的残差函数

下面推导Gradient descent Method和Gauss-Newton Method下得到的优化方法需要的步长$h$。

Gauss-Newton方法推导步长。

各种方法下的步长对比：

（$\rho (h)的表达式$的推导的第二行是把上面的$h_{lm}$带进去消掉$r$）

先看下面红色框框的了解迭代的思路：

巧妙得到$J$的递推式。

(这一页没有讲清楚)

3.1.6 多项式拟合

按照3.1.4. 多项式拟合取

$$y=f(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\dots +a_{1}x+a_0$$

转化为解线性方程组求解$a_j$

3.1.7 拟合的问题

欠拟合和过拟合

多项式摆动

3.2 基于优化的拟合

逻辑回归基于优化的方法（查看本文的3.4.1逻辑回归）：

3.3 基于概率分布的模型

3.4 非线性曲线拟合

3.4.1 逻辑回归

优化方法逻辑回归参考3.2 基于优化的方法

3.5 三角拟合

有周期性的数据，可以用三角拟合。

当数据的周期不等于$2\pi$，进行延拓。

$$\{\begin{array}{c}{\bar{x}}_i=2\pi x_i/P\\ {\bar{y}}_i=y_i\end{array},i=0,1,2,\dots ,n$$

对于拟合确保$n>2m$

三角拟合其实是使用了傅里叶级数，下面展示了连续傅里叶级数和离散的傅里叶级数。

3.5.1 三角拟合的正交性质

我们规定$x_k$是$[-\pi ,\pi ]$之中均匀取点，第一个点是$-\pi$，最后一个点是$\pi$,$x_k=-\pi +\frac{2k\pi }{n},k=1,\cdots ,n$

最小二乘法求系数。

偶函数和基函数下可以简化

3.6 分段三次样条曲线拟合

过点、0阶、1阶、2阶连续，四个条件约束的三次曲线。
和Hermite插值比较：





例子：

3.7 Bezier曲线拟合

3.8 Matlab函数