插值与拟合.docx

1、插值与拟合实验十一：曲线的插值与拟合（设计性实验）一、实验目的1、掌握用MATLAB计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目，对三种插值结果进行初步分析；2、掌握matlab多项式拟合及最小二乘拟合命令的用法；3、通过实例学习用曲线插值和拟合解决实际问题。二、实验原理在大量应用中，人们常面临用一个解析函数描述数据间的关系问题。解决这个问题有两种方法。函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟用插值还是拟合，有时容易确定，有时则并不明显。在插值法里，数据假定是正确的，要求以某种方法描述

2、数据点之间所发生的情况。在曲线拟合或回归方法里，是设法找出某条光滑曲线，使它最佳地拟合数据，而不必经过任何数据点。1、一维插值 插值问题的提法是，已知个节点，其中互不相同，不妨设，求任一插值点处的插值。可以看成是由某个函数产生的，的解析表达式可能十分复杂，或不存在封闭形式。也可以未知。求解的基本思路是，构造一个相对简单的函数，使通过全部节点，即，再由计算插值，即。MATLAB中用interp1()函数来进行一维插值计算，其调用格式为：y=interp1(x0,y0,x,method)其中x0,y0是同维数据向量，分别表示插值节点的横、纵坐标，x是待求函数值的插值节点向量。method为可选项，

3、说明插值使用的方法。对于一维插值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; % index for hour data was recordedtemps=5,8,

4、9,15,25,29, 31,30,22,25,7,24; % recorded temperaturesplot(hours,temps,hours,temps, + ) %view temperaturestitle( Temperature )xlabel( Hour ),ylabel( Degrees Celsius )为计算任给时间的温度，可用函数interp1作插值运算。t=interp1(hours，temps，9.3) %estimate temperature at hour=9.3t =22.9000t=interp1(hours, temps, 4.7)t =22t=in

5、terp1(hours, temps,3.2 6.5 7.1 11.7) t =10.200030.000030.900024.9000若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法是用3次样条插值。t=interp1(hours, temps,9.3, spline ) % estimate temperature at hour=9.3t=21.8577t=interp1(hours, temps, 4.7, spline ) % estimate temperature at hour=4.7t=22.3143 t=interp1(hours, temps,

6、 3.2 6.5 7.1 11.7, spline ) t= 9.673430.042731.175525.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为插值点的竖坐标

7、。 method为插值方法，主要有 linear 线性插值，默认 pchip 逐段三次Hermite插值 spline 逐段三次样条函数插值其中最后一种插值的曲面比较平滑例：设有一平板，在均匀分布的格栅上采集温度值。数据如下：width=1:5; %index for width of plate (i.e. the x-dimension)depth=1:3; %index for depth of plate (i.e. the y-dimension)temps=82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86 %temperature dat

8、atemps =82 81 80 82 8479 63 61 65 8184 84 82 85 86如同在标引点上测量一样，矩阵temps表示整个平板的温度分布。temps的列与下标depth或y-维相联系，行与下标width或x-维相联系。为了估计在中间点的温度，我们必须对它们进行辨识。wi=1:0.2:5; %estimate across width of plated=2; % at a depth of 2zlinear=interp2(width,depth,temps,wi,d);%linear interpolationzcubic=interp2(width, depth,

9、temps, wi,d, cubic ) ; % cubic interpolationplot(wi, zlinear, - , wi, zcubic) % plot resultsxlabel( Width of Plate )， ylabel( Degrees Celsius )title( Temperature at Depth = num2str(d) )3、 曲线拟合设是直角平面坐标系下给出的一组数据，设，且已知它们满足某一函数y=f(a,x),其中a 是待定的参数，曲线拟合就是要确定这些参数，使得拟合值与真实值之间的误差达到最小。最常用的判断曲线拟合效果好坏的方法就是最小二乘法

10、，即使得目标函数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 .9 1;y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30

11、11.2;p=polyfit(x, y, 2)p=-9.8108 20.1293 -0.0317可得其解为y=-9.8108x220.1293x-0.0317。为了将曲线拟合解与数据点比较，可将二者都绘成图。xi=linspace(0, 1, 100); % x-axis data for plottingz=polyval(p, xi);%为计算在xi数据点的多项式值，调用MATLAB的函数polyvalplot(x, y, o , x, y, xi, z, : )多项式阶次的选择是任意的。两点决定一直线或一阶多项式。三点决定一曲线或2阶多项式。按此进行，n+1数据点唯一地确定n阶多项式。于

12、是对于有11个数据点的上例情况，可选一个高达10阶的多项式。然而，高阶多项式给出很差的数值特性，故不应选择比所需阶次高的多项式。拟合曲线阶次越多就越好的观念是不适用的。（2）最小二乘拟合在Matlab的最优化工具箱中提供了lsqcurvefit()函数，可以解决最小二乘曲线拟合的问题。该函数的调用格式为：a,fval=lsqcurvefit(fun,a0,x,y)其中：fun为原型函数的Matlab表示，可以是M文件定义的函数或inline()函数，a0为最优化的初值，x,y为原始输入输出数据向量，a为返回的待定系数向量，fval为在此待定系数下的目标函数的值。例：由函数生成一组数据x和yx=

13、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); xxres绘制出拟合曲线与样本点的图形： x1=0:0.01:10; y1=

14、f(xx,x1); plot(x1,y1,x,y,o)三、实验内容1、一维插值方法的实现。2、二维插值方法的实现。3、多项式拟合命令的使用方法。4、最小二乘拟合命令的使用方法。四、实验报告实验十一：曲线的插值与拟合实验名称： 实验日期： 年 月 日姓名： 班级学号： 成绩： 一、实验目的1、掌握用MATLAB计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目，对三种插值结果进行初步分析；2、掌握matlab多项式拟合及最小二乘拟合命令的用法；3、通过实例学习用曲线插值和拟合解决实际问题。二、实验内容及步骤1、已知数据：x0.1.2.3.4.5.6.7.8.91y.3.511.41.

15、61.9.6.4.81.52画出用线性、三次样条和三次多项式插值所得0,1区间内的曲线图，并求当xi=0.25、0.35、0.45时的yi的值。（1） 画曲线图：程序：x=0:0.1:1;y=0.3 0.5 1 1.4 1.6 1.9 0.6 0.4 0.8 1.5 2;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,y

16、3)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.3 0.5 1 1.4 1.6 1.9 0.6 0.4 0.8 1.5 2;z=interp1(x,y,0.25 0.35 0.45)运行结果：z = 1.2000 1.5000 1.75002、已知某处山区地形选点测量坐标数据为：x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5y

17、=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6海拔高度数据为：z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95

18、84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87（1） 画出其地貌图：程序：x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5; y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6;z=89 90 87 85 92 91 96 93 90 87 82;

19、 92 96 98 99 95 91 89 86 84 82 84; 96 98 95 92 90 88 85 84 83 81 85; 80 81 82 89 95 96 93 92 89 86 86; 82 85 87 98 99 96 97 88 85 82 83; 82 85 89 94 95 93 92 91 86 84 88; 88 92 93 94 95 89 87 86 83 81 92; 92 96 97 98 96 93 95 84 82 81 84; 85 85 81 82 80 80 81 85 90 93 95; 84 86 81 98 99 98 97 96 95

20、84 87; 80 81 85 82 83 84 87 90 95 86 88; 80 82 81 84 85 86 83 82 81 80 82; 87 88 89 98 99 97 96 98 94 92 87;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=89 90 87 85 92 91 96 93 90 87 82 92 96 9

21、8 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 8

22、2 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87;x,y=meshgrid(x,y);plot3(x,y,z,bo)hold onxi=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、由离散数据x0.1.2.3.4.5.

23、6.7.8.91y.3.511.41.61.9.6.4.81.52用3阶多项式拟合数据，并将原始曲线与拟合曲线进行比较。（1）程序： x=0:0.1:1;y=0.3 0.5 1 1.4 1.6 1.9 0.6 0.4 0.8 1.5 2;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.10.

24、20.30.40.52.32012.64702.97073.28853.60080.60.70.80.91.03.90904.21474.51914.82325.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.3201 2.6470 2.9707 3.2885 3.6008 3.9090 4.2147 4.5191 4.8232 5.1275;xx,res=lsqcurvefit(ellfl,1,1,

25、1,1,xi,yi); xxres运行结果：ans = 2.9957 1.0151 2.2540 2.0197res = 1.9194e-004所以用最小二乘拟合的函数为。在此拟合函数下的目标函数的值为。5、在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表：氮施肥量(公斤/公顷)03467101135202259336404471土豆产量(公斤)15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75磷施肥量(公斤/公顷)024497398147196245294342土豆产量(公斤)

26、33.4632.4736.0637.9641.0440.0941.2642.1740.3642.73根据上表数据分别给出土豆产量与氮、磷肥的关系式。(1)画出土豆产量与氮施肥量的散点图，观察它们的大致图形确定多项式的阶数。程序：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;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

27、.17 40.36 42.73;subplot(1,2,1);plot(x1,y1,o); xlabel(氮施肥量); ylabel(土豆产量); hold on;subplot(1,2,2);plot(x2,y2,o);xlabel(磷施肥量); ylabel(土豆产量);运行结果：结论：从散点图判断土豆产量与氮、磷肥的关系式应该采用的模型。从上面的散点图来看，土豆产量与氮、磷肥的关系式都应该采用二次多项式模型。（2）利用MATLAB对数据进行拟合。程序：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;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;p=polyfit(x1,y1,2)p2=polyfit(x2,y2,2)运行结果：p = -0.0003 0.1971 14.7416p2 = -0.0001 0.0719 32.9161结论：用二次多项式去拟合土豆产量与氮肥的关系式为。用二次多项式去拟合土豆产量与磷肥的关系式为。