五、插值与拟合

  如果只知道函数$f(x)$在区间[a,b]上的一系列点上的值：
$$y_i=f(x_i),i=0,1,...,n$$
当需要求其他点上的值时，用一个函数$\varphi (x_i)$代替$f(x)$。
插值： 插值函数$\varphi (x_i)$要满足：$$\varphi (x_i)=y_i,i=0,1,...,n$$
拟合： 求近似函数而不要求过已知数据点，只要求在某种意义下在这些点上的总偏差最小。

5.1 插值

5.1.1 分段线性插值

将相邻点用直线连接即是分段线性插值，插值函数$I_n(x_i)=y_i$，且在每个区间$[x_i,x_{i+1}](i=0,1,...,n-1)$为线性函数。

5.1.2 拉格朗日插值多项式

5.1.3 样条插值

许多问题对插值函数的光滑性又较高要求，要求曲线具有较高光滑程度，要连续且有连续曲率。

1. 样条函数

2. 三次样条插值

利用样条函数进行插值，即为样条插值。
已知函数$f(x)$在区间$[a,b]$上的$n+1$个节点：
$$a=x_0<x_1<...<x_{n-1}<x_n=b$$
上的值$y_i=f(x_i)$，求插值函数$S(x)$，使得：
（1）$S(x_i)=y_i$
（2）在给个区间上为三次多项式，记为$S_i(x)$
（3）在区间$[a,b]$上二阶连续可微。
由条件（2），
$$S(x)=S_i(x),x\in [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$$
$a_i,b_i,c_i,d_i$为待定参数，共4n个。
由条件（3），
$$\begin{gathered} S_i(x_{i+1})=S_{i+1}(x_{i+1}), \\ {S}_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1}), \\ {S}_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1}), \\ i=0,1,...,n-2 \end{gathered}$$
为确定4n个参数，还需两个边界条件。

常用边界条件：
（1）$S^{\prime }(a)=y_0^{\prime },S^{\prime }(b)=y_n^{\prime }$这种边界条件建立的样条插值函数称为完备三次样条插值函数。
（2）$S^{\prime \prime }(a)=y_0^{\prime \prime },S^{\prime \prime }(b)=y_n^{\prime \prime }$。$y_0^{\prime \prime }=y_n^{\prime \prime }=0$时称为自然边界条件。
（3）$S^{\prime }(a+0)=S^{\prime }(b-0)，S^{\prime \prime }(a+0)=S^{\prime \prime }(b-0{)}^{\prime \prime }$ 此条件称为周期条件。

5.1.4 Matlab插值工具

1.一维插值函数

y=interp1(x0,y0,x,'method')

其中：method指定插值方法，默认线性插值。

‘nearest’	最近项插值
‘linear’	线性插值
‘spline’	立方样条插值
‘cubic’	立方插值

要求x0单调。当x0等距时可以使用快速插值法：‘*nearest’…
2.三次样条插值
如无边界条件，常用非扭结(not-a-knot)条件，该条件强迫第1个和第2个三次多项式的三阶导数相等，最后一个和倒数第2个三次多项式也同样。

y=interp1(x0,y0,x,'spline')
y=spline(x0,y0,x);
pp=csape(x0,y0,conds)
pp=csape(x0,y0,conds,valconds);y=fnval(pp,x);

x0,y0是已知点，x是插值点，y是插值点函数值。
conds是插值边界条件：

‘complete’	边界为一阶导数，导数值在valconds中
‘not-a-knot’	非扭结条件
‘periodic’	周期条件
‘second’	边界为二阶导数，导数值在valconds中，默认[0,0]
‘variational’	设置边界二阶导数值为[0,0]

对于一些特殊边界条件，将conds设为1*2的矩阵，元素值为0，1，2，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件。

clc,clear
x0=[0 3 5 7 9 11 12 13 14 15];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x=0:0.1:15;
y1=interp1(x0,y0,x);
y2=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0);
y3=fnval(pp1,x);
pp2=csape(x0,y0,'second');
y4=fnval(pp2,x);
[x',y1',y2',y3',y4']
subplot(1,3,1)
plot(x0,y0,'+',x,y1)
title('Pircewise linear')
subplot(1,3,2)
plot(x0,y0,'+',x,y2)
title('Splinel')
subplot(1,3,3)
plot(x0,y0,'+',x,y3)
title('Spline2')
dx=diff(x);
dy=diff(y3);
dy_dx=dy./dx;
dy_dx0=dy_dx(1)
ytemp=y3(131:151);
ymin=min(ytemp);
index=find(y3==ymin);
xmin=x(index);
[xmin,ymin]

3.二维插值
1）插值点为网格节点
已知$m\times n$个节点$(x_i,y_j,z_{ij})(i=1,2,...,m;j=1,2,...,n),且x_1<...<x_m;y_1<...<y_n$
求点$(x,y)$处的插值。

z=interp2(x0,y0,z0,x,y,'method')

其中，x0,y0为m维和n维向量，z0维n*m矩阵。x,y为一维数组，表示插值点，且是方向不同的向量，即一个行向量，一个列向量；z为矩阵，行数为y的维数，列数为x的维数；'method‘的用法同一维插值。

%三次样条插值：
pp=csape({x0,y0},z0,conds,valconds)
z=fnval(pp,{x,y})

其中，x0,y0为m维和n维向量，z0维m*n矩阵。z为矩阵，行数为x的维数，列数为y的维数；

2）插值点为散乱节点
已知n个节点$(x_i,y_i,z_i)(i=1,2,...,n),$
求点$(x,y)$处的插值$z$。

ZI=griddata(x,y,z,XI,YI)

其中：x,y,z均为n维向量；XI和YI是方向不同的向量，即一个行向量，一个列向量；

5.2 曲线拟合

5.2.1 线性最小二乘法

5.2.2 Matlab实现最小二乘法

1.解方程组法

clc,clear
x=[19;25;31;38;44];
y=[19.0;32.3;49.0;73.3;97.8];
r=[ones(5,1),x.^2];
a=r\y;
xi=19:0.1:44;
yi=a(1)+a(2)*xi.^2;
plot(x,y,'o',xi,yi,'r')

2.多项式拟合法
取{r1(x),...,rm+1(x)}={1,x,...,xm}\{r_1(x),...,r_{m+1}(x)\}=\{1,x,...,x^m\}{r1​(x),...,rm+1​(x)}={1,x,...,xm}

a=polyfit(x0,y0,m)

a为拟合多项式$y=a(1)x^m+...+a(m+1)$的系数。
多项式在x处的值用函数计算：

y=polyval(a,x)

5.3 最小二乘优化

5.3.1 lsqlin函数

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

5.3.2 lsqcurvefit函数

x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)

5.3.3 lsqnonlin函数

x=lsqnonlin(fun,x0,lb,ub,options)

5.3.4 lsqnoneg函数

x=lsqnonneg(C,d,options)

5.3.5 Matlab曲线拟合图形界面

cftool

1. 数据导入工作空间
2. 运行cftool，打开用户图形界面窗口
3. 选择适当的模型进行拟合
4. 生成相关统计量，进行预测

5.4 曲线拟合与函数逼近

syms x
base=[1,x^2,x^4];
y1=base.'*base
y2=cos(x)*base.'
r1=int(y1,-pi/2,pi/2)
r2=int(y2,-pi/2,pi/2)
a=r1\r2
xinshu1=double(a) %符号数据转为数值数据
xinshu2=vpa(a,6) %符号数据转为小数型的符号数据