数据插值与拟合
插值与插值函数：已知由(可能未知或非常复杂)产生的一批离散数据，且个互异插值节点，在插值区间内寻找一个相对简单的函数，使其满足下列插值条件：
再利用已求得的计算任一非插值节点的近似值，这就是插值。其中称为插值函数，称为被插函数。
?最小二乘拟合：已知一批离散的数据，互不相同，寻求一个拟合函数，使与的误差平方和在最小二乘意义下最小。在最小二乘意义下确定的称为最小二乘拟合函数。
1）Lagrange插值法
????a．待定系数法：假设插值多项式，利用待定系数法即可求得满足插值条件的插值函数。关键在于确定待定系数。
????b．利用基函数的构造方法首先构造个满足条件：的次插值基函数，再将其线性组合即可得如下的Lagrange插值多项式：
其中
????c．Lagrange插值余项
注：上述两种构造方法所得的Lagrange插值多项式是一样的，即满足插值条件的Lagrange插值多项式是唯一的。
2）分段线性插值
????作分段线性插值的目的在于克服Lagrange插值方法可能发生的不收敛性缺点。所谓分段线性插值就是利用每两个相邻插值节点作线性插值，即可得如下分段线性插值函数：
其中
特点：插值函数序列具有一致收敛性，克服了高次Lagrange插值方法的缺点，故可通过增加插值节点的方法提高其插值精度。但存在于节点处不光滑、插值精度低的缺点。
3）三次样条插值
????三次样条插值的目的在于克服Lagrange插值的不收敛性和提高分段线性插值函数在节点处的光滑性。所谓三次样条插值方法就是在满足下列条件：
????a．
????b．在每个子区间上是三次多项式的三次样条函数中寻找满足如下插值条件：
以及形如等边界条件的插值函数的方法。
????特点：三次样条插值函数序列一致收敛于被插函数，因此可通过增加节点的方法提高插值的精度。
4）插值方法的Matlab实现
一维数据插值
MATLAB中用函数interp1来拟合一维数据，语法是
YI=INTERP1(X,Y,XI,方法)
其中（X，Y）是已给的数据点，XI是插值点，
其中方法主要有
''linear''-线性插值，默认
''pchip''-逐段三次Hermite插值
''spline''-逐段三次样条函数插值
其中最后一种插值的曲线比较平滑
例：
x=0:.12:1;x1=0:.02:1;
y=(x.^2-3x+5).exp(-5x).sin(x);
plot(x,y,''o'');holdon;
y1=interp1(x,y,x1,''spline'');
plot(x1,y1,'':'')
如果要根据样本点求函数的定积分，而函数又是比较光滑的，则可以用样条函数进行插值后再积分，在MATLAB中可以编写如下程序：
functiony=quadspln(x0,y0,a,b)
f=inline(‘interp1(x0,y0,x,’’spline’’)’,’x’,’x0’,’y0’);
y=quadl(f,a,b,1e-8,[],x0,y0);
现求six(x)在区间[0，pi]上的定积分，只取5点
x0=[0,0.4,1,2,pi];
y0=sin(x0);
I=quadspln(x0,y0,0,pi)
结果得到的值为2.01905，精确值为2
二元函数插值：
MATLAB中用函数interp2来拟合二维网格（X，Y）上的数据Z，语法是
YI=INTERP2(X,Y,Z，XI,YI，方法)
其中（X，Y，Z）是已给的数据点，（XI，YI）是插值点坐标，
其中方法主要有
''linear''-线性插值，默认
''pchip''-逐段三次Hermite插值
''spline''-逐段三次样条函数插值
其中最后一种插值的曲面比较平滑
例：[x,y]=meshgrid(-3:.6:3,-2:.4:2);
z=(x.^2-2x).exp(-x.^2-y.^2-x..y);
[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);%生成网格，x1和y1均为同样size的矩阵
z1=interp2(x,y,z,x1,y1,’spline’);%z1是矩阵，size和x1，y1相同
surf(x1,y1,z1);
axis([-3,3,-2,2,-0.7,1.5]);
如果数据不是在网格上取的，则可用函数griddata来解决
语法是
YI=griddata(X,Y,Z，XI,YI，‘v4’)
其中（X，Y，Z）是已给的数据点，（XI，YI）是插值点坐标，
其中除了方法‘v4’外还有
''linear''-线性插值，默认
''cublc''-逐段三次Hermite插值
''nearest''
其中‘v4’方法比较好
例x=-3+6rand(200,1);%生成随机点的x坐标向量x
y=-2+4rand(200,1);%生成随机点的y坐标向量y
z=(x.^2-2x).exp(-x.^2-y.^2-x.y);%上述点的样本值向量z
[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);%生成网格，x1和y1均为同样size的矩阵
z1=griddata(x,y,z,x1,y1,’v4’);
surf(x1,y1,z1);
axis([-3,3,-2,2,-0.7,1.5]);
生成的图类似上图。
最小二乘拟合
假设有一组数据，xi,yi,且已知它们满足某一函数y=f(a,x),其中a是待定的参数，最小二乘拟合就是要确定这些参数，使得目标函数sum(yi-f(a,xi))^2)达到最小。在MATLAB的最优化工具箱中有函数lsqcurvefit,语法如下
[a，J]=lsqcurvefit(fun,a0,x,y,Lb,Ub,options),其中fun是内联函数，也可以是m函数，但要写为@fun的形式，a0是参数的初值，x,y为数据向量。J=min(sum(yi-f(a,xi))^2))
Lb和Ub是a的下上限
例
x=0:.1:10;
y=0.12exp(-0.213x)+0.54exp(-0.17x).sin(1.23x);
f=inline(''a(1)exp(-a(2)x)+a(3)exp(-a(4)x).sin(a(5)x)'',''a'',''x'');
[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y)，
运行得a向量为0.11970.21250.54040.17021.2300
最小二乘为7.1637e-007
如果要提高精度，可把上面程序的最后一句改为以下语句。首先自定义选项
Ff=optimset;
ff.TolFun=1e-20;ff.TolX=1e-15,
[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff);
题目：拟合问题
????在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表：
氮施肥量(公斤/公顷) 0 34 67 101 135 202 259 336 404 471 土豆产量(公斤) 15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75 磷施肥量(公斤/公顷) 0 24 49 73 98 147 196 245 294 342 土豆产量(公斤) 33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73 根据上表数据分别给出土豆产量与氮、磷肥的关系式。
建模分析：
????使用Matlab语言首先画出土豆产量与氮施肥量的散点图，从图可看出土豆产量与氮肥量的关系是二次函数关系，因此可选取拟合函数为：
，
其中x和y分别为氮肥量和土豆产量，a、b和c为待定系数。再画出磷肥量与土豆产量的散点图，从图可看出从0到98、从98到342之间分别呈明显的线性关系。由此可选取所求拟合函数为一分段的线性函数，换言之，用前5点作一线性拟合函数，再用后6个点也作一线性拟合函数，最后用两个线性函数求出其分界点即可得分段线性函数。
数学模型：
????对氮肥的拟合函数为：
????对磷肥的拟合函数为：
x1=[0,34,67,101,135,202,259,336,404,471];
y1=[15.18,21.36,25.72,32.29,34.03,39.45,43.15,43.46,40.83,30.75];
plot(x1,y1,''r+'')
aa=polyfit(x1,y1,2)
xx=0:471;
yy=aa(1)xx.xx+aa(2)xx+aa(3);
plot(xx,yy,x1,y1,''r+'')
%Nextproblem
x2=[0,24,49,73,98,147,196,245,294,342];
y2=[33.46,32.47,36.06,37.96,41.04,40.09,41.26,42.17,40.36,42.73];
plot(x2,y2,''r+'')
a1=polyfit(x2(1:5),y2(1:5),1)
a2=polyfit(x2(5:10),y2(5:10),1)
x0=(a2(2)-a1(2))/(a1(1)-a2(1))
xx1=0:x0;yy1=a1(1)xx1+a1(2);
xx2=x0:342;yy2=a2(1)xx2+a2(2);
plot(x2,y2,''r+'',xx1,yy1,xx2,yy2)
程序运行结果
8
|
|
