插值与拟合.docx
《插值与拟合.docx》由会员分享，可在线阅读，更多相关《插值与拟合.docx（16页珍藏版）》请在冰豆网上搜索。
插值与拟合
实验十一：
曲线的插值与拟合（设计性实验）
一、实验目的
1、掌握用MATLAB计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目，对三种插值结果进行初步分析；
2、掌握matlab多项式拟合及最小二乘拟合命令的用法；
3、通过实例学习用曲线插值和拟合解决实际问题。
二、实验原理
在大量应用中，人们常面临用一个解析函数描述数据间的关系问题。
解决这个问题有两种方法。
函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。
而面对一个实际问题，究竟用插值还是拟合，有时容易确定，有时则并不明显。
在插值法里，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。
在曲线拟合或回归方法里，是设法找出某条光滑曲线，使它最佳地拟合数据，而不必经过任何数据点。
1、一维插值
插值问题的提法是，已知
个节点
，其中
互不相同，不妨设
，求任一插值点
处的插值
。
可以看成是由某个函数
产生的，
的解析表达式可能十分复杂，或不存在封闭形式。
也可以未知。
求解的基本思路是，构造一个相对简单的函数
，使
通过全部节点，即
，再由
计算插值，即
。
MATLAB中用interp1（）函数来进行一维插值计算，其调用格式为：
y=interp1（x0,y0,x,’method’）
其中x0,y0是同维数据向量，分别表示插值节点的横、纵坐标，x是待求函数值的插值节点向量。
’method’为可选项，说明插值使用的方法。
对于一维插值interp1来说，matlab提供的可选的方法有：
nearest,linear,spline,cubic，它们分别表示最近插值，线性插值，三次样条插值和三次插值。
命令返回值y是插值曲线在节点向量x（横坐标）处的纵坐标向量。
三次样条插值命令：
y=interp1（x0,y0,x,'spline'）
或
y=spline（x0,y0,x）
其中输入x0,y0,x和输出y的意义同上。
为说明一维插值，考虑下列问题，12小时内每一小时测一次室外温度。
数据存储在两个变量中。
hours=1:
12;%indexforhourdatawasrecorded
temps=[5,8,9,15,25,29,31,30,22,25,7,24];%recordedtemperatures
plot（hours,temps,hours,temps,'+'）%viewtemperatures
title（'Temperature'）
xlabel（'Hour'）,ylabel（'DegreesCelsius'）
为计算任给时间的温度，可用函数interp1作插值运算。
t=interp1（hours，temps，9.3）%estimatetemperatureathour=9.3
t=22.9000
t=interp1（hours,temps,4.7）
t=22
t=interp1（hours,temps,[3.26.57.111.7]）
t=10.2000
30.0000
30.9000
24.9000
若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。
最常用的方法是用3次样条插值。
t=interp1（hours,temps,9.3,'spline'）%estimatetemperatureathour=9.3
t=21.8577
t=interp1（hours,temps,4.7,'spline'）%estimatetemperatureathour=4.7
t=22.3143
t=interp1（hours,temps,[3.26.57.111.7],'spline'）
t=9.6734
30.0427
31.1755
25.3820
注意，样条插值得到的结果，与上面所示的线性插值的结果不同。
因为插值是一个估计或猜测的过程，其意义在于应用不同的估计规则导致不同的结果。
.2、二维插值
二维插值是基于与一维插值同样的基本思想，是对两个变量的函数z=f（x，y）进行插值。
MATLAB中用函数interp2来拟合二维网格（X，Y）上的数据Z，语法是：
z1=interp2（x,y,z,x1,y1,’method’）
其中（x,y,z）是已给的数据点的横、纵、竖坐标，（x1，y1）是插值点的横、纵坐标，z1为插值点的竖坐标。
‘method’为插值方法，主要有
'linear'线性插值，默认
'pchip'逐段三次Hermite插值
'spline'逐段三次样条函数插值
其中最后一种插值的曲面比较平滑
例：
设有一平板，在均匀分布的格栅上采集温度值。
数据如下：
width=1:
5;%indexforwidthofplate（i.e.thex-dimension）
depth=1:
3;%indexfordepthofplate（i.e.they-dimension）
temps=[8281808284;7963616581;8484828586]%temperaturedata
temps=
8281808284
7963616581
8484828586
如同在标引点上测量一样，矩阵temps表示整个平板的温度分布。
temps的列与下标depth或y-维相联系，行与下标width或x-维相联系。
为了估计在中间点的温度，我们必须对它们进行辨识。
wi=1:
0.2:
5;%estimateacrosswidthofplate
d=2;%atadepthof2
zlinear=interp2（width,depth,temps,wi,d）;%linearinterpolation
zcubic=interp2（width,depth,temps,wi,d,'cubic'）;%cubicinterpolation
plot（wi,zlinear,'-',wi,zcubic）%plotresults
xlabel（'WidthofPlate'），ylabel（'DegreesCelsius'）
title（['TemperatureatDepth='num2str（d）]）
3、曲线拟合
设
是直角平面坐标系下给出的一组数据，设
，且已知它们满足某一函数y=f（a,x）,其中a是待定的参数，曲线拟合就是要确定这些参数，使得拟合值与真实值之间的误差达到最小。
最常用的判断曲线拟合效果好坏的方法就是最小二乘法，即使得目标函数sum（yi-f（a,xi））^2）达到最小。
当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。
（1）多项式拟合
MATLAB中进行多项式拟合的函数为polyfit（）。
调用格式：
p=polyfit（x,y,n）
[p,s]=polyfit（x,y,n）
说明：
x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。
矩阵s用于生成预测值的误差估计。
如
x=[0.1.2.3.4.5.6.7.8.91];
y=[-.4471.9783.286.167.087.347.669.569.489.3011.2];
p=polyfit（x,y,2）
p=-9.810820.1293-0.0317
可得其解为y=-9.8108x2＋20.1293x-0.0317。
为了将曲线拟合解与数据点比较，可将二者都绘成图。
xi=linspace（0,1,100）;%x-axisdataforplotting
z=polyval（p,xi）;%为计算在xi数据点的多项式值，调用MATLAB的函数polyval
plot（x,y,'o',x,y,xi,z,':
'）
多项式阶次的选择是任意的。
两点决定一直线或一阶多项式。
三点决定一曲线或2阶多项式。
按此进行，n+1数据点唯一地确定n阶多项式。
于是对于有11个数据点的上例情况，可选一个高达10阶的多项式。
然而，高阶多项式给出很差的数值特性，故不应选择比所需阶次高的多项式。
拟合曲线阶次"越多就越好"的观念是不适用的。
（2）最小二乘拟合
在Matlab的最优化工具箱中提供了lsqcurvefit（）函数，可以解决最小二乘曲线拟合的问题。
该函数的调用格式为：
[a,fval]=lsqcurvefit（fun,a0,x,y）
其中：
fun为原型函数的Matlab表示，可以是M文件定义的函数或inline（）函数，a0为最优化的初值，x,y为原始输入输出数据向量，a为返回的待定系数向量，fval为在此待定系数下的目标函数的值。
例：
由函数
生成一组数据x和y
x=0:
.1:
10;
y=0.5*exp（-0.1*x）+0.8*exp（-0.3*x）.*sin（2.5*x）;
假设以上数据满足原型为
，其中，
为待定系数。
采用最小二乘曲线拟合的目的就是获得这些待定系数的值，使得目标函数的值最小。
根据已知的函数原型，可以编写出如下函数：
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）;
xx'
res
绘制出拟合曲线与样本点的图形：
x1=0:
0.01:
10;y1=f（xx,x1）;plot（x1,y1,x,y,'o'）
三、实验内容
1、一维插值方法的实现。
2、二维插值方法的实现。
3、多项式拟合命令的使用方法。
4、最小二乘拟合命令的使用方法。
四、实验报告
实验十一：
曲线的插值与拟合
实验名称：
实验日期：
年月日
姓名：
班级学号：
成绩：
一、实验目的
1、掌握用MATLAB计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目，对三种插值结果进行初步分析；
2、掌握matlab多项式拟合及最小二乘拟合命令的用法；
3、通过实例学习用曲线插值和拟合解决实际问题。
二、实验内容及步骤
1、已知数据：
x
0
.1
.2
.3
.4
.5
.6
.7
.8
.9
1
y
.3
.5
1
1.4
1.6
1.9
.6
.4
.8
1.5
2
画出用线性、三次样条和三次多项式插值所得[0,1]区间内的曲线图，并求当xi=0.25、0.35、0.45时的yi的值。
（1）画曲线图：
程序：
x=0:
0.1:
1;
y=[0.30.511.41.61.90.60.40.81.52];
subplot（2,2,1）
plot（x,y,'b+'）
title（'散点图'）
subplot（2,2,2）
x2=0:
0.05:
1;
y2=interp1（x,y,x2）;
plot（x2,y2）
title（'线性插值'）
subplot（2,2,3）
x3=0:
0.05:
1;
y3=interp1（x,y,x3,'spline'）;
plot（x3,y3）
title（'三次样条插值'）
subplot（2,2,4）
x4=0:
0.05:
1;
p=polyfit（x,y,3）;
y4=polyval（p,x4）;
plot（x4,y4）
title（'三次多项式插值'）
运行结果：
（2）求xi=0.25、0.35、0.45时的yi的值
程序：
x=0:
0.1:
1;
y=[0.30.511.41.61.90.60.40.81.52];
z=interp1（x,y,[0.250.350.45]）
运行结果：
z=
1.20001.50001.7500
2、已知某处山区地形选点测量坐标数据为：
x=00.511.522.533.544.55
y=00.511.522.533.544.555.56
海拔高度数据为：
z=8990878592919693908782
9296989995918986848284
9698959290888584838185
8081828995969392898686
8285879899969788858283
8285899495939291868488
8892939495898786838192
9296979896939584828184
8585818280808185909395
8486819899989796958487
8081858283848790958688
8082818485868382818082
8788899899979698949287
（1）画出其地貌图：
程序：
x=[00.511.522.533.544.55];
y=[00.511.522.533.544.555.56];
z=[8990878592919693908782;
9296989995918986848284;
9698959290888584838185;
8081828995969392898686;
8285879899969788858283;
8285899495939291868488;
8892939495898786838192;
9296979896939584828184;
8585818280808185909395;
8486819899989796958487;
8081858283848790958688;
8082818485868382818082;
8788899899979698949287];
[xi,yi]=meshgrid（0:
0.5:
5,0:
0.5:
6）;
mesh（xi,yi,z）
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）;
title（'地貌图'）
运行结果：
（2）对数据插值加密形成地貌图，原始数据用小圆圈标出。
（将程序补充完整）
程序：
x=0:
0.5:
5;
y=0:
0.5:
6;
z=[8990878592919693908782
9296989995918986848284
9698959290888584838185
8081828995969392898686
8285879899969788858283
8285899495939291868488
8892939495898786838192
9296979896939584828184
8585818280808185909395
8486819899989796958487
8081858283848790958688
8082818485868382818082
8788899899979698949287];
[x,y]=meshgrid（x,y）;
plot3（x,y,z,'bo'）
holdon
xi=linspace（0,5,50）;
yi=linspace（0,6,80）;
[x1,y1]=meshgrid（xi,yi）;
z1=interp2（x,y,z,x1,y1,'cubic'）;
mesh（x1,y1,z1）
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）;
title（'地貌图'）
运行结果：
3、由离散数据
x
0
.1
.2
.3
.4
.5
.6
.7
.8
.9
1
y
.3
.5
1
1.4
1.6
1.9
.6
.4
.8
1.5
2
用3阶多项式拟合数据，并将原始曲线与拟合曲线进行比较。
（1）程序：
x=0:
0.1:
1;
y=[0.30.511.41.61.90.60.40.81.52];
p=polyfit（x,y,3）;
x1=0:
0.01:
1;
y1=polyval（p,x1）;
plot（x,y,'b',x1,y1,'r'）
legend（'原始数据','拟合数据'）
运行结果：
（2）从图像上观察拟合的效果。
从上面两幅图来看，用三次多项式拟合的效果并不是很好。
4、已知数据可能满足
，求满足数据的最小二乘解a,b,c,d的值，并将原始曲线与拟合曲线进行比较。
0.1
0.2
0.3
0.4
0.5
2.3201
2.6470
2.9707
3.2885
3.6008
0.6
0.7
0.8
0.9
1.0
3.9090
4.2147
4.5191
4.8232
5.1275
（1）建立函数文件e11f1.m，用来存储函数
。
程序：
ellfl=inline（'a
（1）*x+a
（2）*x.^2.*exp（-a（3）*x）+a（4）','a','x'）;
（2）用最小二乘拟合函数拟合数据。
程序：
xi=0.1:
0.1:
1;
yi=[2.32012.64702.97073.28853.60083.90904.21474.51914.82325.1275];
[xx,res]=lsqcurvefit（ellfl,[1,1,1,1],xi,yi）;
xx'
res
运行结果：
ans=
2.9957
1.0151
2.2540
2.0197
res=
1.9194e-004
所以用最小二乘拟合的函数为
。
在此拟合函数下的目标函数的值为
。
5、在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表：
氮施肥量（公斤/公顷）
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
根据上表数据分别给出土豆产量与氮、磷肥的关系式。
（1）画出土豆产量与氮施肥量的散点图，观察它们的大致图形确定多项式的阶数。
程序：
x1=[03467101135202259336404471];
y1=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];
x2=[024497398147196245294342];
y2=[33.4632.4736.0637.9641.0440.0941.2642.1740.3642.73];
subplot（1,2,1）;plot（x1,y1,'o'）;
xlabel（'氮施肥量'）;
ylabel（'土豆产量'）;holdon;
subplot（1,2,2）;plot（x2,y2,'o'）;
xlabel（'磷施肥量'）;
ylabel（'土豆产量'）;
运行结果：
结论：
从散点图判断土豆产量与氮、磷肥的关系式应该采用的模型。
从上面的散点图来看，土豆产量与氮、磷肥的关系式都应该采用二次多项式模型。
（2）利用MATLAB对数据进行拟合。
程序：
x1=[03467101135202259336404471];
y1=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];
x2=[024497398147196245294342];
y2=[33.4632.4736.0637.9641.0440.0941.2642.1740.3642.73];
p=polyfit（x1,y1,2）
p2=polyfit（x2,y2,2）
运行结果：
p=
-0.00030.197114.7416
p2=
-0.00010.071932.9161
结论：
用二次多项式去拟合土豆产量与氮肥的关系式为
。
用二次多项式去拟合土豆产量与磷肥的关系式为
。
升级会员 