插值与拟合

实际问题中，仅能观测到函数f(x)在某区间[a,b]上一系列点的值$y_i=f(x_i),i=0,1,\cdots ,n$。想得知观测点之间的函数值时，可用插值或拟合进行近似。

• 插值：用较简单的、满足一定条件的函数$\phi (x)$代替$f(x)$，且插值函数满足$\phi (x_i)=y_i,i=0,1,\cdots ,n$
• 拟合：拟合的近似函数不要求过已知数据点，只要求在某种意义下在这些点上的总偏差最小

1. 插值方法

1.1 分段线性插值
将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记为$I_n(x)$，满足$I_n(x_i)=y_i$，且$I_n(x)$在每个小区间$[x_i,x_{i+1}]$上是线性函数。$I_n(x)$可以表示为$I_n(x)={\sum }_{i=0}^n{y}_i{l}_i(x)$，其中$$l_i(x)=\{\begin{array}{l}\frac{x-x_{i-1}}{x_i-x_{i-1}},x\in [x_{i-1},x_i],i\ne 0\\ \frac{x-x_{i+1}}{x_i-x_{i+1}},x\in [x_i,x_{i+1}],i\ne n\\ 0,其他\end{array}$$用$I_n(x)$计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关。但n越大，分段越多，插值误差越小。

1.2 拉格朗日插值多项式
拉格朗日插值的基函数为$$l_i(x)=\prod _{\genfrac{}{}{0ex}{}{j=0}{j\ne i}}^n\frac{x-x_j}{x_i-x_j},i=0,1,\cdots ,n$$其中$l_i(x)$是n次多项式，满足$$l_i(x_j)=\{\begin{array}{l}0,j\ne i\\ 1,j=i\end{array}$$拉格朗日插值函数为$$L_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)=\sum _{i=0}^n{y}_i(\prod _{\genfrac{}{}{0ex}{}{j=0}{j\ne i}}^n\frac{x-x_j}{x_i-x_j})$$

1.3 样条插值
许多工程问题对插值函数的光滑性有较高的要求，要求曲线不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。

1.3.1 样条函数的概念
数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间[a,b]的一个划分$$\Delta :a=x_0<x_1<\cdots <x_{n-1}<x_n=b$$如果函数S(x)满足：

• 在每个小区间$[x_i,x_{i+1}]$上S(x)是m次多项式
• S(x)在[a,b]上具有m-1阶连续导数

则称S(x)为关于划分$\Delta$的m次样条函数，其图形为m次样条曲线。折线是一次样条曲线。

1.3.2 三次样条插值
取插值函数为样条函数，即为样条插值。记函数S(x)为f(x)的三次样条插值函数，且有 S ( x ) = { S i ( x )   , x ∈ [ x i , x i + 1 ]   , i = 0 , 1 , ⋯   , n − 1 } S i ( x ) = a i x 3 + b i x 2 + c i x + d i S(x)=\left\{S_i(x)\ ,x\in [x_i,x_{i+1}]\ ,i=0,1,\cdots,n-1\right\}\\ S_i(x)=a_ix^3+b_ix^2+c_ix+d_i S(x)={Si​(x) ,x∈[xi​,xi+1​] ,i=0,1,⋯,n−1}Si​(x)=ai​x3+bi​x2+ci​x+di​式中$a_i,b_i,c_i,d_i$为待定系数，共4n个。由二阶连续可微的性质有$$\{\begin{array}{l}{S}_i(x_{i+1})=S_{i+1}(x_{i+1})\\ S{’}_i(x_{i+1})=S{’}_{i+1}(x_{i+1})\\ S’{’}_i(x_{i+1})=S’{’}_{i+1}(x_{i+1})\end{array},i=0,1,\cdots ,n-2$$联立$S(x_i)=y_i(i=0,1,\cdots ,n)$共有4n-2个方程，为确定S(x)的4n个待定参数，需再给出两个边界条件。
常用的三次样条函数的边界条件由3种类型：

1. $S’(a)=y{’}_0,S’(b)=y{’}_n$。由这种边界条件建立的样条插值函数称为f(x)的完备三次样条插值函数。
   如果f’(x)不知道，可要求S’(x)与f’(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’(a)=N_a^{\prime }(a),S’(b)=N{’}_b(b)$。由这种边界条件建立的三次样条称为f(x)的Lagrange三次样条插值函数。
2. $S’’(a)=y’{’}_0,S’’(b)=y’{’}_n$。特别地，$y’{’}_0=y’{’}_n=0$时，称为自然边界条件。
3. $S’(a+0)=S’(b-0),S’’(a+0)=S’’(b-0)$，此条件称为周期条件。

1.4 Matlab插值工具箱
1.4.1 一维插值函数

y=interp1(x0,y0,x,’method’)		%其中x0,y0为已知数据点,x为插值点，y为插值点的函数值，method指定插值的方法，默认为线性插值，其值可以为：
	‘nearest’最近项插值	‘linear’线性插值	‘spline’立方样条插值	‘cubic’立方插值	所有的插值方法要求x0是单调的

1.4.2 三次样条插值
如果三次样条插值没有边界条件，最常用的方法就是采用非扭结（not-a-knot）条件。该条件强迫第一个和第二个三次多项式的三阶导数相等，最后一个和倒数第二个三次多项式的三阶导数相等。
Matlab中三次样条插值由如下函数：

y=interp1(x0,y0,x,’spline’)
y=spline(x0,y0,x)
pp=csape(x0,y0)	%提倡使用该函数，csape的返回值是pp形式，要得到插值点的函数值，需调用函数fnval
	y=fnval(pp,x)，默认边界条件为Lagrange边界条件
pp=csape(x0,y0,conds,valconds)	%conds指定插值的边界条件，当边界设为二阶导数时，valconds给定二阶导数的值

1.4.3 二维插值
若节点是二维的，插值函数就是二元函数，即曲面。

• 插值节点为网格节点

z=interp2(x0,y0,z0,x,y,’method’)	%x0,y0分别为m维和n维向量，表示节点；z0为n*m矩阵，表示节点值
%x,y为一维数组，表示插值点；z为矩阵，其行数为y的维数，列数为x的维数，表示得到的插值
pp=csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})	%三次样条插值

• 插值节点为散乱节点

ZI=griddata(x,y,z,XI,YI)	%x,y,z均为n维向量，指明所给数据点的横、纵、竖坐标
%向量XI,YI为给定网格点的横坐标和纵坐标，ZI为网格(XI,YI)处的函数值

2. 曲线拟合的线性最小二乘法

2.1 线性最小二乘法
已知一组数据，即平面上n个点$(x_i,y_i)$，要寻找一个函数f(x)使曲线拟合最好。线性最小二乘法是解决曲线拟合最常用的方法，基本思路是，令$$f(x)=a_1{r}_1(x)+a_2{r}_2(x)+\cdots +a_m{r}_m(x)$$其中$r_k(x)$为事先选定的一组线性无关的函数，$a_k$为待定系数，m<n。拟合准则是使$y_i$与$f(x_i)$的距离${\delta }_i$的平方和最小，称为最小二乘准则。

2.1.1 系数$a_k$的确定
记$$J(a_1,\cdots ,a_m)=\sum _{i=1}^n{\delta }_i^2=\sum _{i=1}^n[f(x_i)-y_i{]}^2$$要使J达到最小，即要$\frac{\partial J}{\partial a_j}=0$，可得到关于$a_1,\cdots .a_m$的线性方程组，矩阵化后为$${\mathbf{R}}^{\mathbf{T}}\mathbf{R}\mathbf{A}={\mathbf{R}}^{\mathbf{T}}\mathbf{Y}$$式中$$\mathbf{R}={\left[\begin{array}{ccc}{r}_1(x_1)& \cdots & r_m(x_1)\\ ⋮& ⋮& ⋮\\ r_1(x_n)& \cdots & r_m(x_n)\end{array}\right]}_{n\times m},\mathbf{A}=[a_1,\cdots ,a_m{]}^T,\mathbf{Y}=[y_1,\cdots ,y_n{]}^T$$当 { r 1 ( x ) , ⋯   , r m ( x ) } \left\{r_1(x),\cdots,r_m(x)\right\} {r1​(x),⋯,rm​(x)}线性无关时，$\mathbf{R}$列满秩，${\mathbf{R}}^{\mathbf{T}}\mathbf{R}$可逆，于是方程组式有唯一解$$\mathbf{A}=({\mathbf{R}}^{\mathbf{T}}\mathbf{R}{)}^{-1}{\mathbf{R}}^{\mathbf{T}}\mathbf{Y}$$
2.1.2 函数$r_k(x)$的选取
若通过机理分析，能够知道y和x之间的函数关系，则$r_1(x),\cdots ,r_m(x)$容易确定。若无法知道y和x之间的关系，通常可以将数据$(x_i,y_i),i=1,2,\cdots ,n$作图，直观地判断应该用什么样的曲线取作拟合。常用的曲线有：

• 直线：$y=a_{1}x+a_2$
• 多项式：$y=a_1{x}^m+\cdots +a_{m}x+a_{m+1}$
• 双曲线：$y=\frac{a_1}{x}+a_2$
• 指数曲线：$y=a_1{e}^{a_{2}x}$，用指数曲线拟合前需作变量代换，化为对$a_1,a_2$的线性函数

2.2 最小二乘法的Matlab实现
2.2.1 解方程组方法
Matlab中的线性最小二乘标准型为$\underset{A}{\min }||\mathbf{R}\mathbf{A}-\mathbf{Y}|{|}_2^2$，求解命令为$\mathbf{A}=\mathbf{R}\setminus \mathbf{Y}$。

2.2.2 多项式拟合方法
若取 { r 1 ( x ) , ⋯   , r m + 1 ( x ) } = { 1 , x , ⋯   , x m } \left\{r_1(x),\cdots,r_{m+1}(x)\right\}=\left\{1,x,\cdots,x^m\right\} {r1​(x),⋯,rm+1​(x)}={1,x,⋯,xm}，即用m次多项式拟合给定数据，Matlab中有现成的函数

a=polyfit(x0,y0,m)		%x0,y0为要拟合的数据，m为拟合多项式的次数，a为拟合多项式的系数向量
y=polyval(a,x)		%计算多项式在x出的值

3. 最小二乘优化

在非线性规划问题中，若目标函数由若干个函数的平方和构成，把极小化这类函数的问题称为最小二乘优化问题$$\min F(x)=\sum _{i=1}^m{f}_i^2(x)$$在Matlab优化工具箱中，用于求解最小二乘优化问题的函数由lsqlin、lsqcurvefit、lsqnonlin、lsqnonneg。

1. lsqlin函数：求解$$\begin{gathered} \underset{x}{\min }\frac{1}{2}||\mathbf{C}\cdot \mathbf{x}-\mathbf{d}|{|}_2^2 \\ \text{s.t.}\{\begin{array}{l}\mathbf{A}\cdot \mathbf{x}⩽\mathbf{b}\\ Aeq\cdot \mathbf{x}=beq\\ lb⩽\mathbf{x}⩽ub\end{array} \end{gathered}$$式中$\mathbf{C},\mathbf{A},Aeq$为矩阵，$\mathbf{d},\mathbf{x},\mathbf{b},beq,lb,ub$为向量。
   MATLAB中的函数为

x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)

2. lsqcurvefit函数：给定输入输出数列xdata，ydata，求参量x，使得$$\underset{x}{\min }||F(x,xdata)-ydata|{|}_2^2=\sum _i(F(x,xdat{a}_i)-ydat{a}_i{)}^2$$MATLAB中的函数为

x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)	%fun为定义函数F(x,xdata)的M文件

3. lsqnonlin函数：已知函数向量$F(\mathbf{x})=[f_1(\mathbf{x}),\cdots ,f_k(\mathbf{x}){]}^T$，求x使得$$\underset{x}{\min }||F(\mathbf{x})|{|}_2^2$$MATLAB中的函数为

x = lsqnonlin(fun,x0,lb,ub,options)		%fun为定义向量函数F(x)的M文件

4. lsqnonneg函数：求解非负的x，使得$$\underset{x}{\min }||\mathbf{C}\cdot \mathbf{x}-\mathbf{d}|{|}_2^2$$MATLAB中的函数为

x = lsqnonneg(C,d,options)

4. 曲线拟合与函数逼近

曲线拟合是选择一个较简单函数f(x)在一定准则下最接近一组离散数据。函数逼近是选择一个较简单函数f(x)在一定准则下最接近一个较复杂的连续函数y(x),[a,b]。与曲线拟合的最小二乘准则对应，函数逼近常用最小平方逼近，即$$J={\int }_a^b[f(x)-y(x){]}^{2}dx$$达到最小。与曲线拟合一样，选一组函数 { r k ( x ) , k = 1 , ⋯   , m } \left\{r_k(x),k=1,\cdots,m\right\} {rk​(x),k=1,⋯,m}构造f(x)，即令$$f(x)=a_1{r}_1(x)+\cdots +a_m{r}_m(x)$$利用极值必要条件可得$$\left[\begin{array}{ccc}\left(r_1,r_1\right)& \cdots & \left(r_1,r_m\right)\\ ⋮& ⋮& ⋮\\ \left(r_m,r_1\right)& \cdots & \left(r_m,r_m\right)\end{array}\right]\left[\begin{array}{c}{a}_1\\ ⋮\\ a_m\end{array}\right]=\left[\begin{array}{c}\left(y,r_1\right)\\ ⋮\\ \left(y,r_m\right)\end{array}\right]$$其中$(g,h)={\int }_a^{b}g(x)h(x)dx$，当方程组的系数矩阵非奇异时，有唯一解。