第九讲 数据插值与拟合

1、第九讲第九讲 曲线拟合与插值曲线拟合与插值 在工程实践和科学实验中，常常需要从一组实验观测数据 niyx ii , 1 , 0),(揭示自变量x与因变量y之间的关系， 一般可以用一个近似的函数关系式yf（x）来表示 通常可以采用两种方法：曲线拟合和插值 拟合主要是考虑到观测数据受随机误差的影响，寻求整体 误差最小、较好反映观测数据的近似函数，并不保证所得 到的函数一定满足 )( ii xfy 曲线拟合的目的是根据实验获得的数据去建立因变量 与自变量之间有效的经验函数关系，为进一步的深入 研究提供线索 插值函数一般是已知函数的线性组合或者称为加权平 均插值在工程实践和科学实验中有着非常广泛而又十

2、 分重要的应用，例如，信息技术中的图像重建、图像放 大中为避免图像的扭曲失真的插值补点、建筑工程的外 观设计。化学工程实验数据与模型的分析、天文观测数 据、地理信息数据的处理如(天气预报）以及社会经济 现象的统计分析等等 插值则要求函数在每个观测点处一定要满足 )( ii xfy 1、船在该海域会搁浅吗船在该海域会搁浅吗 在某海域测得一些点(x，y)处的水深z（单位：英尺）由 下表给出，水深数据是在低潮时测得的船的吃水深度 为5英尺，问在矩形区域(75,200)*(-50，150）里的哪些 地方船要避免进人 一、实例及其模型 分析分析 由于测量点是散乱分布的，先在平面上作出测量点的分 布图，再

3、利用二维插值方法补充一些点的水深，然后作 出海底曲面图和等高线图，并求出水深小于5的海域范 围 在化学反应中，为研究某化合物的浓度随时间的变化规律， 测得一组数据如表 2、浓度的变化规律、浓度的变化规律 表中的数据反映了浓度随时间变化的函数关系，它是一 种离散关系若需要推断20，40分钟时的浓度值，能否用 一个显函数y=f(t)来拟合表中的离散数据，然后再计算浓 度值f(20), f（40）? 问题分析问题分析 （1）首先将这些离散数据分布在直角坐标系下，由此可 发现浓度与时间之间呈现什么规律这种数据分布在 直角坐标系下的图形被称为散点图； （2）根据散点图，判段它接近于哪类函数曲线， 即确定

4、函数形式 （3）函数形式确定以后，关键是要确定函数中含有的 待定参数。 最常用的确定待定系数的方法是，曲线拟合的最小二乘法 二、二、 插值与拟合插值与拟合 1、插值方法、插值方法 （1）分段线性插值）分段线性插值 分段线性插值的提法如下： （2）分段三次埃尔米特插值分段三次埃尔米特插值 在插值问题中，如果除了插值节点的函数值给定外，还 要求在节点的导数值为给定值，即插值问题变为 相当于在每一小段上应满足四个条件（方程），可以确 定四个待定参数三次多项式正好有四个系数，所以可 以考虑用三次多项式函数作为插值函数，这就是分段三分段三 次埃尔米特插值次埃尔米特插值，它与分段线性插值一起都称为分段多

5、项式插值 （3）三次样条插值）三次样条插值 2、曲线拟合的最小二乘法、曲线拟合的最小二乘法 给定平面上的点, 2 , 1),(niyxi 进行曲线拟合有多种方法，最小二乘法是解决曲线拟 合最常用的一种方法 最小二乘法的原理是求f(x),使 n i ii n i i yxf 1 2 1 2 )( 达到最小 简单地说，最小二乘法准则就是使所有散点到曲线的距 离平方和最小 线性最小二乘法线性最小二乘法 拟合函数可由一些简单的“基函数”（例如幂函数，三 角 函数等等） )(,),(),( 10 xxx n 来线性表示 )()()()( 1100 xcxcxcxf mm 现在要确定系数 , 10m cc

6、c 使达到极小为此 三、插值的matlab实现 MATLAB中的插值函数为interp1，其调用格式为 ) ,( 1intmethodxiyxerpyi 其中x，y为插值点，yi为在被插值点xi处的插值结果， x, y为向量。 注意：所有的插值方法都要求x是单调的，并且xi不 能够超过x的范围。 linear spline cubic nearest method ) ,( 1intmethodxiyxerpyi MATLAB提供的插值方法有几种 表示采用的插值方法 ：分段线性插值 pchip：三次Hermite插值（立方插值） ：三次分段样条插值 ：最近点等值方式 缺省时表示线性插值 例例1

7、 在一 天24小时内，从零点开始每间隔2小时测得的环 境温度数据分别为 12，9，9，1，0，18 ，24，28，27，25，20，18，15，13， 推测中午（即13点）时的温度 x=0：2：24； y=12 9 9 10 18 24 28 27 25 20 18 15 13； x113 ； y1interp1（x，y，x1，spline） 若要得到一天24小时的温度曲线 x=0：2：24； y=12 9 9 10 18 24 28 27 25 20 18 15 13 xi0：13600：24； yi=interp1(x,y,xi,spline ); plot(x, y, o, xi, yi

8、) 2、高维插值、高维插值 N维插值函数interpN（） 其中N可以为2，3，如N2为二维插值，调用格式为 ) ,(2intmethodyixizyxerpzi 其中 x，y，z为插值节点，zi为被插值点(xi，yi)处的插值结果 且， xi， yi为被插值节点构成的新的网格数据 methods代表的意思和可选择的插值方法和前面一样 注意：注意：所有的插值方法都要求x和y是单调的网格，x和 y可以 是等距的也可以是不等距的 (1) 网格数据插值问题 例例2 气旋变化情况的可视化 下表是气象学家测量得到的气象资料，它们分别表示在南 半球地时按不同纬度。不同月份的平均气旋数字根据这 些数据，绘制

9、出气旋分布曲面图形 y=5:10:85;x=1:12; x,y=meshgrid(x,y); plot(x,y,*); pause z=2.4,1.6,2.4,3.2,1.0,0.5,0.4,0.2,0.5,0.8,2.4,3.6; 18.7 21.4 16.2 9.2 2.8 1.7 1.4 2.4 5.8 9.2 10.3 16; 20.8 18.5 18.2 16.6 12.9 10.1 8.3 11.2 12.5 21.1 23.9 25.5; 22.1 20.1 20.5 25.1 29.2 32.6 33.0 31.0 28.6 32.0 28.1 25.6; 37.3 28.8

10、27.8 37.2 40.3 41.7 46.2 39.9 35.9 40.3 38.2 43.4; 48.2 36.6 35.5 40 37.6 35.4 35 34.7 35.7 39.5 40 41.9; 25.6 24.2 25.5 24.6 21.1 22.2 20.2 21.2 22.6 28.5 25.3 24.3; 5.3 5.3 5.4 4.9 4.9 7.1 5.3 7.3 7 8.6 6.3 6.6; 0.3,0,0,0.3,0,0,0.1,0.2,0.3,0,0.1,0.3; figure surf(x,y,z) pause xi,yi=meshgrid(1:12,5:

11、1:85); zi=interp2(x,y,z,xi,yi,spline ); figure mesh(xi,yi,zi) xlabel(月份), ylabel(纬度), zlabel(气旋), axis(0 12 0 90 0 50) title(南半球气旋可视化图形) （2）、一般二维分布的数据插值）、一般二维分布的数据插值 在实际应用问题中，大部分的数据以实测的多组 （xi,yi,zi）给出，所以不能直接使用interp2（）函数。 Matlab中提供了另一个函数griddata( )，用来专 门解决这类问题。其调用格式如下 Z=griddata(x,y,z,x0,y0,method)

12、x,y,z是已知样本点的坐标，可以是任意分布的。 X0,y0是期望的插值位置，即被插值节点, 可以是单点， 向量或者网格型矩阵 插值方法，除了上面的 方法外，还有一个是4.0版本提供 的一个插值方法，选项为v4 四、曲线拟合的四、曲线拟合的matlab实现实现 1、已知函数原型的、已知函数原型的 （1）多项式拟合）多项式拟合 假设已知函数原型为 11 nn n axaxay Matlab提供的拟合函数为 a=polyfit(xdata,ydata,n) 其中n表示多项式的最高阶数，xdata，ydata为将要拟合 的数据，它是用数组的方式输入 输出参数a为拟合多项式 的系数, 11 nn aa

13、aa 注：注：多项式在x处的值y可用下面程序计算 y=polyval(a,x) T=19.1 25.0 30.1 36.0 40.0 45.1 50.0; R=76.30 77.80 79.25 80.80 82.35 83.90 85.10; PR=polyfit(T,R,1); t=10:60; r=polyval(PR,t); plot(T,R,*,t,r) 解：Matlab程序 （2）一般函数线性组合的曲线拟合）一般函数线性组合的曲线拟合 假设已知函数原型为 )()()()( 1100 xcxcxcxf mm 通过求解线性方程可得待定系数，一般方法： X= %已知数据x的列向量 Y=

14、%已知数据y的列向量 A=f1(X),f2(X),fm(X) %系数矩阵，fm（）为基函数 c=Ay 解：matlab程序 X=0 0.1 0.2 0.3 0.4 0.5 0.6 Y=2 2.20254 2.40715 2.61592 2.83096 3.05448 3.28876 A=ones(size(X),exp(X),exp(-X); c1=AY; C=c1 x=0:0.05:1; y=C(1)+C(2)*exp(x)+C(3)*exp(-x); plot(X,Y,*,x,y) （3）一般的曲线拟合）一般的曲线拟合 假设已知函数原型是一般的函数，可以是多项式，可以 是线性，也可以是非线

15、性的，一般情况下用这个来求解 非线性情况 Matlab在优化工具箱中提供的求解一般的曲线拟合函 数lsqcurvefit(),其调用格式如下 其中Fun表示函数Fun(p,data)的M函数文件，p0表示 函数的初值.。 p=lsqcurvefit(Fun,p0,xdata,ydata) 注：若要求解点x处的函数值可用程序f=Fun(p,x)计算 (1)函数原型m文件 function y=fname(a,t) y=a(1)*exp(-a(2)*t); 注：注：因为后面可能要用到计算函数在一些点上的值， 因此写函数原型表达式时，记得用点运算点运算 tk=0.2 0.3 0.4 0.5 0.6

16、0.7 0.8; Ik=3.16 2.38 1.75 1.34 1.00 0.74 0.56; a=lsqcurvefit(fname,1,1,tk,Ik) x=0:0.05:1; y=fname(a,x); plot(tk,Ik,*,x,y) 运行结果：a = 5.6361 2.8906 (2)拟合程序 2、函数原型未知、函数原型未知 已知一组数据，用什么样的曲线拟合最好呢？可以根据 散点图进行直观判断，在此基础上，选择几种曲线分别 拟合，然后观察哪条曲线的最小二乘指标最小。 图(a)，数据接近于直线，故宜采用线性函数y=a+bx拟合； 解解：先将表中的数据用曲线表示 X=1.1052 1.2214 1.3499 1.4918 1.6487 1.8221 2.0138 2.2255 2.4596 2.7183 3.6693; Y=0.6795 0.6006 0.5309 0.4693 0.4148 0.3666 0.3241 0.2865 0.2532 0.2238