插值法与数据拟合法

《插值法与数据拟合法》由会员分享，可在线阅读，更多相关《插值法与数据拟合法（56页珍藏版）》请在装配图网上搜索。

1、第七讲插值方法与数据拟合 7.1引言在工程和科学实验中，常常需要从一组实验观测数据(Xi , yi ) (i = 1,2,n)揭示自变量x与因变量y之间的关系，一般可以用一个近似的函数关系式y = f (x)来表示。函数f (x)的产生办法因观测数据与要求的不同而异，通常可采用两种方法：插值与数据拟合。 7.1.1 插值方法1 .引例1已经测得在北纬32.3 海洋不同深度处的温度如下表:表 7.1.1深度x (m)46671495014221634水温y (C )7.044.283.402.542.13根据这些数据，我们希望能合理地估计出其它深度(如500 米、600 米、1000米)处的水温

2、。解决这个问题，可以通过构造一个与给定数据相适应的函数来解决，这是一个被称为插值的问题。2 插值问题的基本提法对于给定的函数表xX0X1Xny = f (x)yoy1-yn其中f(x)在区间a,b上连续，xo,xi,，xn为a, b上n 1个互不相同的点，要求在一个性质优良、便于计算的函数类P(x)中，选出一个使P(xi ) = yi, i = 0, 1,，n(7.1.1)成立的函数P(x)作为f (x)的近似，这就是最基本的插值问题(见图 7.1.1 )。为便于叙述，通常称区间a,b为插值区间，称点xo,X1,，xn为插值节点，称函数类P(x)为插值函数类，称式(7.1.1)为插值条件，称函

3、数P(x)为插值函数，称f (x)为被插函数。求插值函数 P(x) 的方法称为插值法。 7.1.2 数据拟合1 .引例2在某化学反应中，已知生成物的浓度与时间有关。今测得一组数据如下：表 7.1.2时间t (分)12345678浓度y 10 34.006.408.008.809.229.509.709.86时间t (分)910111213141516浓度y 10 310.0010.2010.3210.3210.5010.5510.5810.60根据这些数据，我们希望寻找一个 y = f (t)的近似表达式(如建立浓度 y与时间t之间的经验公式等)。 从几何上看，就是希望根据给定的一组点(1,4

4、.00 ),(16, 10.60 ),求函数y = f (t)的图象的一条拟合曲线。2 .数据拟合问题的基本提法xX0X1Xny = f (x)yoy1-yn对于给定的函数表其中f （x）在区间a, b上连续，xo, xi,，xn为a, b上n 1个互不相同的点，要求找一个简单合理的函数近似表达式（x）,使 （x）与f （x）在某种准则下最为接近，这就是最基本的数据拟合问题（见图 7.1.2 ）。通常，我们称（x）为给定数据点的拟合函数。图7.1.1插值问题示意图图7.1.2数据拟合问题示意图 7.1.3插值方法与数据拟合的基本理论依据插值方法与数据拟合的基本理论依据，就是数学分析中的 Wei

5、erstrass定理：设函数f （x）在区间a,b上连续，则对 0，存在多项式P（x），使得maxjx） p（x）即：有界区间上的连续函数被多项式一致逼近。 7.1.4实际应用中两种方法的选择在实际应用中，究竟选择哪种方法比较恰当？总的原则是根据实际问题的特点来决定采用哪一种方 法。具体说来，可从以下两方面来考虑：1 如果给定的数据是少量的且被认为是严格精确的，那么宜选择插值方法。采用插值方法可以保证 插值函数与被插函数在插值节点处完全相等。2 如果给定的数据是大量的测试或统计的结果，并不是必须严格遵守的，而是起定性地控制作用的，那么宜选用数据拟合的方法。这是因为，一方面测试或统计数据本身往往

6、带有测量误差，如果要求所得的 函数与所给数据完全吻合，就会使所求函数保留着原有的测量误差；另一方面，测试或统计数据通常很多，如果采用插值方法，不仅计算麻烦，而且逼近效果往往较差。 7.2一维数据的基本插值方法简介插值函数类的取法很多，可以是代数多项式，也可以是三角多项式或有理函数；可以是a, b上任意光滑函数，也可以是分段光滑函数。在此介绍最基本、最常用的两种插值方法：分段多项式插值与三次 样条插值，及其Matlab实现。 7.2.1 一维数据的分段多项式插值对于给定的一维数据xxoX1Xny = f(x)yoyiyn分段多项式插值就是求一个分段(共n段)多项式P(x),使其满足P(xi )

7、= yi (i = 0, 1,，n)或更高的要求。一般地，分段多项式插值中的多项式都是低次多项式(不超过三次)1 .分段线性插值分段线性插值函数 Pi (x)是一个分段一次多项式(分段线性函数)。在几何上就是用折线代替曲线，如图7.2.1 ,故分段y线性插值亦称为折线插值。其插值公式为X XiXXi 1f/j , XXi , Xi +1 P(X)i 1yXi 1XiXiXi 1(7.2.1)0 xo X1Xn-1Xn2 .分段二次插值图7.2.1分段线性插值示意图分段二次插值函数 P2 (x)是一个分段二次多项式。在几何上就是分段抛物线代替曲线y = f (x)，故分段二次插值又称为分段抛物插

8、值。其插值公式为P2(X)(X Xi )(xXi 1)yi 1(Xi 1 Xi )(Xi 1Xi 1)(XXi 1)(x Xi1)(XiXi 1)(XiXi 1) yi(X Xi1)(x Xi)yi 1(Xi 1Xi1)(Xi1 Xi)Xi -1 , Xi +1 (7.2.2)i 1 i 1 X xjyk, xk i 1 j i 1 XkXjj k3 .三次Hermite 插值三次Hermite插值问题的基本提法一：已知一维数据XxoX1y = f (x)yoy1(723)mom i求一个三次多项式 P3(X)，使之满足P3 (xi ) = yi, P3 (xi ) = mi, i = 0,

9、1i(x)，使之满足构造三次插值基函数o(x)，i(x)，o(x)，i(Xj)o，iji, iji (Xj)0,i,j o,ii(Xj)0,o, i (Xj)i,(7.2.4)利用这四个插值基函数，取三次多项式P3 (X)P3 (x) =o(x) yoi(x) yio(x) m oi (x) m i(7.2.5)将插值条件(7.2.3)式代入，可推得:o(x)X XoXiXoi(x)2 XXiXoXi2Xi2XXoXi Xo2xo(7.2.6)o(x)(xX。)xiXoXii(x)(x2、XXoXi)XiXo(725)、(7.2.6)两式构成了三次 Hermite插值基本提法一的插值公式。三次

10、Hermite插值问题的基本提法二：已知一维数据xoxiX2y = f (x)yoyiy2求一个三次多项式 P3(X)，使之满足(727)P3 (Xi ) = yi, i = 0, 1,2 , P3 (xi ) = mi构造三次插值基函数0(X), 1(X),2 (X), 1(X),使之满足0, i j(7.2.8)J, i (Xi)0, i,j 0,1,21, i j1(xJ 0,1 (X1)1利用这四个插值基函数，取三次多项式P3 (X)为P3 (x) =0(x) y01 (x) y12(x) y21(x) m1(7.2.9)将插值条件(7.2.7)式代入，可推得:0(X)(XX1)2(x

11、X2)(X。X1)2(x。X2)1(X)(XX0)(XX2) 1 (x x) 1 1(X1X0)2(X11 (xx1)X2)X1X0X1X2(7.2.10)3(X)(XX)(xX1)2(X2X)(X2X1)21(X)(Xx)(x xj(x X2)(X1 X)(X1X2)(729)、(7.2.10)两式构成了三次 Hermite插值基本提法二的插值公式。 7.2.2一维数据的三次样条插值上述介绍的分段多项式插值，其优点为计算简单、稳定性好、收敛性有保证，且易于在计算机上实现。 但它也明显存在着缺陷。它只能保证在每个小区间段Xi , Xi +1内光滑，在各小区间连接点 Xi处连续，却60年代开始，

12、首先起源与航空、造船业等工程不能保证整条曲线的光滑、光顺性，难以满足某些工程的要求。对于象高速飞机的机翼形线，船体放样等 型值线往往要求有二阶光滑度，即有二阶连续导数。而由设计的实际需要而发展起来的样条插值，既保留了分段多项式插值的各种优点，又提高了插值函数的光滑 度。在此，仅介绍应用最广且具有二阶连续导数的三次样条插值方法。1.三次样条插值问题的基本提法对于给定的一维数据XxoX1Xny = f (x)yoy1-Yn求一个三次多项式S(x)满足条件(1) S(xi) = yi, i = 0, 1,，n ；(2) S(x)具有二阶连续导数，特别在节点Xi上应满足连续性要求，即对i = 0, 1

13、,n有S(Xi 0)S(xi0)S0)S(Xi0)S。0)S(xi 0)2 .三次样条插值函数给定区间a, b的一个划分 ：a = xo X1 =466&x1(i)714 &x1(i)950&x1(i)=1422 a4=1,x1(i),(x1(i)A2)/2,(x1(iF3)/6,(x1(i)-x(2)A3)/6,(x1(i)-x(3)A3)/6,0;y1(i)=a4*alpha;elsea5=1,x1(i),(x1(i)A2)/2,(x1(i)A3)/6,(x1(i)-x(2)A3)/6,(x1(i)-x(3)A3)/6,(x1(i)-x (4) )A3)/6;y1(i)=a5*alpha;

14、endendy1 7.3 二维数据的基本插值方法简介对于二维数据的插值，首先要考虑两个问题：一是二维区域是任意区域还是规则区域，二是给定的数 据是有规律分布的还是散乱的、随机分布的。第一个问题比较容易处理。目前的插值方法基本上是基于规则区域的，对于不规则区域，只需将其， 划分为规则区域或扩充为规则区域来讨论即可。对于第二个问题，当给定的数据是有规律分布时，方法较 多也较成熟；而给定的数据是散乱的、随机分布时，没有固定的方法，但一般的处理思想是：从给定的数 据出发，依据一定的规律恢复出规则分布点上的数据，转化为数据分布有规律的情形来处理。二维数据插值的方法也有很多。在此，针对给定数据有规律分布和

15、散乱分布两种情形，简单介绍双三 次样条插值方法和改进的 Shepard 方法(反距离平方法)的基本概念和基本思想，及其 Matlab 实现。 7.3.1 双三次样条插值双三次样条插值方法，是用来解决规则区域上给定数据有规律分布的插值问题的常用方法。设R： a, b c, d是xy平面上的一个矩形区域。在x轴和y轴上分别取定分割x： a = xo xi Xnb , y : c = yo yi ym 0，令由于（）是可微函数，使得如下定义的F(x,y)274R0,F(x，y)在性能上有所改善rk 0NWk(x,y)fk, rk 0，k 1其中 Wk(x,y)(厲)N(rk)k 1(7.3.3)按照

16、上述的思想，可从给定的数据恢复出规则分布点上的数据，接下来就可应用双三次样条插值或其 它的二维数据插值方法来处理。 7.3.3 二维数据插值的 Matlab实现1 规则区域上给定数据有规律分布的二维插值数据形式为：y1y 2-y nX1X11Z12Z1nx 2Z21Z22-Z2nx mZm1Zm2Zmn插值函数为：interp2()。其调用格式为zi = in terp2(x, y, z, xi, yi, methos),其中x , y , z 为插值节点，均为向量;zi为被插值点(xi, yi)处的插值结果；methos为米用的插值方法：nearest:表示最临近插值，linear：表示双线

17、性插值，cubic:表示双二次插值，spline：表示双三次样条插值。注意：上述methos中所有的插值方法都要求 x和y是单调的网格，x和y可以是等距的也可以是不等距的。2 规则区域上给定数据散乱或随机分布的二维插值数据形式为：g y1)(X2, y2)(Xn, yn)Z1Z2Zn插值函数为：eOlsef和e01sff ,。通常两者配合使用，其调用格式为fnodes, a, rnw, b, c = e01sef(x, y, z)pf(l, j), ifail = e01sff(x, y, z, rnw, fnodes, px(j), py(i)其中x , y , z为插值节点，均为向量；px

18、(j) , py(i)为被插值节点;Pf(i, j)为被插值点(px(j), py(i)处的插值结果；其它输出参数涉及插值算法，可以不用了解。eOlsef的输出fnodes和rnw为确定插值的参数，它们是eOlsff需要的输入参数，因此两函数需配合使用。3 例例1气旋变化情况的可视化表7.3.1是气象学家测量得到的气象资料，它们分别表示在南半球地区按不同纬度、不同月份的平均气旋数字。根据这些数据，绘制出气旋分布曲面图形。表 7.3.10 1010 2020 3030 4040 5050 6060 7070 8080 901月2.418.720.822.137.348.225.65.30.32月

19、1.621.418.520.128.836.624.25.303月2.416.218.220.527.835.525.55.404月3.29.216.625.137.24024.64.90.35月1.02.812.929.240.337.621.14.906月0.51.710.132.641.735.422.27.107月0.41.48.333.046.23520.25.30.18月0.22.411.231.039.934.721.27.30.28月0.55.812.528.625.935.722.670.310月0.89.221.132.040.339.528.58.6011月2.410.3

20、23.928.138.24025.36.30.112月3.61625.525.643.441.924.36.60.3解 下面分别用最邻近插值、双线性插值、双三次插值和双三次样条插值，给出不同月份按纬度变化 的气旋值（插值结果），并作出可视化图形如下。最邻近插值的南半就气旋可视化图飛双线性插值的南半球W族可视化图略双三次插值的南半球岂能可视化图形双三次禅条掘11的南半球宅能可视化S惑纬S0 0 月份图732 四种插值方法的可视化图形最邻近插值的Matlab程序为：x=1:12;y=5:10:85;z=2.4, 18.7, 20.8, 22.1,37.3, 48.2, 25.6, 5.3, 0.3

21、1.6, 21.4, 18.5, 20.1,28.8, 36.6, 24.2, 5.3, 0 2.4, 16.2, 18.2, 20.5, 27.8, 35.5, 25.5, 5.4, 03.2, 9.2, 16.6, 25.1, 37.2, 40, 24.6, 4.9, 0.3 1.0, 2.8, 12.9, 29.2, 40.3, 37.6, 21.1, 4.9, 00.5, 1.7, 10.1, 32.6, 41.7, 35.4, 22.2, 7.1, 00.4, 1.4, 8.3, 33.0, 46.2, 35, 20.2, 5.3, 0.10.2, 2.4, 11.2, 31.0,

22、39.9, 34.7, 21.2, 7.3, 0.20.5, 5.8, 12.5, 28.6, 25.9, 35.7, 22.6, 7, 0.30.8, 9.2, 21.1, 32.0, 40.3, 39.5, 28.5, 8.6, 02.4, 10.3, 23.9, 28.1, 38.2, 40, 25.3, 6.3, 0.13.6, 16, 25.5, 25.6, 43.4, 41.9, 24.3, 6.6, 0.3;xi,yi=meshgrid(1:12,5:1:85);zi=interp2(x,y,z,xi,yi,nearest);figuresurf(xi,yi,zi)xlabel(

23、月份),ylabel(纬度),zlabel(气旋),axis(0 12 0 90 0 50)title( 南半球气旋可视化图形 )双线性插值、双三次插值、双三次样条插值的 Matlab 程序为：分别将最邻近线性插值程序中的zi=interp2(x,y,z,xi,yi,nearest)改写为zi=interp2(x,y,z,xi,yi,linear)zi=interp2(x,y,z,xi,yi,cubic)zi=interp2(x,y,z,xi,yi,spline)例 2 水道测量数据（ AMCM86A 题）在某海域测得一些点（X, y）处的水深z （单位：英尺）由表 732给出，水深数据是在低

24、潮时测得的。船的吃水深度为 5英尺，问在矩形区域（75, 200） （ 50, 150）里的哪些地方船要避免进入。表7.3.2 水道水深测量数据（单位：英尺）X129.0140.0108.588.0185.5195.0105.5y7.5141.528.0147.022.5137.585.5z4868688X157.5107.577.081.0162.0162.0117.5y6.581.03.056.566.584.038.5z9988949解 （1 ）假设由题目给出的信息是很少的，除了14个位置的水深之外一无所知。显然，题目要求我们找出水深不到5英尺的区域。为了讨论方便，下面三个假设是合理的：

25、 所给数据是精确的； 讨论区域的海底曲面是光滑的，更确切地说，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底,不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形。 水深是一个按区域来划分的变量，在某个位置的水深与其周围区域的水深是相互依赖的，但这种 依赖作用随距离的增大而减小。就我们讨论的问题来说，每一个给定数据点影响周围的每一个未知点，一 个给定数据点离未知点越近，作用就越大。（2 ）问题分析根据假设，海底曲面是连续光滑的，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形，因而很自然的想法就是用某种光滑的

26、拟合曲面去逼近已知的14 个数据点或以 14 个已知的数据点为基础， 利用二维插值补充一些点的水深，以求得水深不超过 5 米的区域。在此，我们采用二维插值方法，应用 Matlab 程序，作出矩形区域 （75 ， 200） （ 50 ，150） 范围内 的海底地形图、 水深不超过 5 米的危险区域的平面图以及水深不超过 5 米的危险区域的海底地貌图， 并求 出水深不超过 5 米的危险海域范围。（ 3）问题求解采用改进的 Shepard 方法，利用 Matlab 软件作出作出矩形区域 （75 ，200） （ 50 ，150） 范围内的 海底地形图、水深不超过 5 米的危险区域的平面图以及水深不超过

27、 5 米的危险区域的海底地貌图（见图 7.3.3 ），并求出水深不超过 5 米的危险海域范围为： 115 ，200 3 ，119 。（ 4）求解的 Matlat 程序如下：clear;x=129,140,108.5,88,185.5,195,105.5,157.5,107.5,77,81,162,162,117.5; y=7.5,141.5,28,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-38.5;subplot（2,2,1）,plot（x,y,+）,title（ 测量点分布图 ）; z=-4,-8,-6,-8,-6,-8,-8,-9,-9,-

28、8,-8,-9,-4,-9;fnodes,minnq,rnw,rnq,ifail=e01sef（x,y,z）;nx=100;px=linspace(75,200,nx);(75i 200(-60 150)范围内的海底地形图150；-, +100+十50. +*0+ +50100. +1十2DO50100150 2W水深不超辻占米的危脸反域的平面图水深不超讨吕米的危脸区域的海底地粽图150501001502D0图图7.3.3ny=200;py=li nspace(-50,150 ,n y);for i=1: nyfor j=1: nxpf(i,j),ifail=e01sff(x,y,z,r nw

29、,fnodes,px(j),py(i);endendsubplot(2,2,2),meshz(px,py,pf+5),title(75, 200)x(-50 , 150)范围内的海底地形图);a,b=fi nd(pf=-5);amin=min(a);amax=max(a);bmin=min(b);bmax=max(b);xmin=75+(200-75)/100)*bminxmax=75+(200-75)/100)*bmaxymin=-50+(150+50)/200)*aminymax=-50+(150+50)/200)*amaxfor k=1:length(b)i0(k)=75+(200-75

30、)/100)*b(k);endfor k=1:length(a)j0(k)=-50+(150+50)/200)*a(k);endsubplot(2,2,3),plot(i0,j0,+),title(水深不超过 5 米的危险区域的平面图 );i1,j1=find(pfo, b o ）。（2 ）拟合运算首先，分别用二、三、六次多项式拟合，计算得输出参数分别为p1 = 0.0445 , 1.0711 , 4.3252p2 = 0.00600.1963 , 2.1346 , 2.5952p3 = 0.0000 , 0.0004 ,0.0103 , 0.1449 ,1.1395 , 4.9604 , 0

31、.0498拟合函数分别为(X)=0.04451.0711 x 4.3252 x2(x) = 0.0060(x) = 0.0004 x 0.0103 x20.1963 x 2.1346 x22.5952 x30.1449 x3 1.1395 x4 4.9304 x5 0.0498 x6;x其次，再用有理分式 （x）拟合，计算得输出参数分别为ax bp = 0.0841，0.1392拟合函数为最后，用指数函数(x)(x)t0.0841t0.1392aex拟合，计算得输出参数分别为p = 11.3578，1.0873拟合函数为1.0873(x)11.3578e三种方式五个种函数的拟合曲线见图7.4.2 7.4.4 。图7.4.2多项式函数拟合曲线图图743 有理分式函数拟合曲线图图744 指数函数拟合曲线图(3 )误差分析和给定的16组数据比较，三种方式五个函数拟合的误差见下表：表7.4.1五个函数拟合的误差表偏差平方