插值法与数据拟合法.doc

第七讲 插值方法与数据拟合 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插值问题的基本提法对于给定的函数表xx0x1xny = f (x)y0y1yn其中f (x) 在区间 a, b 上连续，x0，x1，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 为插值区间，称点x0，x1，xn为插值节点，称函数类 P(x) 为插值函数类，称式 (7.1.1) 为插值条件，称函数 P(x) 为插值函数，称f (x) 为被插函数。求插值函数 P(x) 的方法称为插值法。 7.1.2 数据拟合1 引例2 在某化学反应中，已知生成物的浓度与时间有关。今测得一组数据如下： 表7.1.2时间t（分）12345678浓度y10-34.006.408.008.809.229.509.709.86时间t（分）910111213141516浓度y10-310.0010.2010.3210.3210.5010.5510.5810.60根据这些数据，我们希望寻找一个y = f (t) 的近似表达式（如建立浓度y与时间t之间的经验公式等）。从几何上看，就是希望根据给定的一组点（1, 4.00），,（16, 10.60），求函数y = f (t) 的图象的一条拟合曲线。2 数据拟合问题的基本提法对于给定的函数表xx0x1xny = f (x)y0y1yn其中f (x) 在区间 a, b 上连续，x0，x1，xn为 a, b 上n + 1个互不相同的点，要求找一个简单合理的函数近似表达式 j (x)，使 j (x) 与f (x) 在某种准则下最为接近，这就是最基本的数据拟合问题（见图7.1.2）。 通常，我们称 j (x) 为给定数据点的拟合函数。 图7.1.1 插值问题示意图图 7.1.2 数据拟合问题示意图 7.1.3 插值方法与数据拟合的基本理论依据插值方法与数据拟合的基本理论依据，就是数学分析中的Weierstrass定理：设函数f (x) 在区间 a, b 上连续，则对 e 0，存在多项式P(x)，使得。即：有界区间上的连续函数被多项式一致逼近。 7.1.4 实际应用中两种方法的选择在实际应用中，究竟选择哪种方法比较恰当？总的原则是根据实际问题的特点来决定采用哪一种方法。具体说来，可从以下两方面来考虑：1如果给定的数据是少量的且被认为是严格精确的，那么宜选择插值方法。采用插值方法可以保证插值函数与被插函数在插值节点处完全相等。2如果给定的数据是大量的测试或统计的结果，并不是必须严格遵守的，而是起定性地控制作用的，那么宜选用数据拟合的方法。这是因为，一方面测试或统计数据本身往往带有测量误差，如果要求所得的函数与所给数据完全吻合，就会使所求函数保留着原有的测量误差；另一方面，测试或统计数据通常很多，如果采用插值方法，不仅计算麻烦，而且逼近效果往往较差。 7.2 一维数据的基本插值方法简介插值函数类的取法很多，可以是代数多项式，也可以是三角多项式或有理函数；可以是 a, b 上任意光滑函数，也可以是分段光滑函数。在此介绍最基本、最常用的两种插值方法：分段多项式插值与三次样条插值，及其Matlab实现。 7.2.1 一维数据的分段多项式插值对于给定的一维数据xx0x1xny = f (x)y0y1yn分段多项式插值就是求一个分段（共n段）多项式P(x)，使其满足P(xi ) = yi（i = 0, 1, , n）或更高的要求。一般地，分段多项式插值中的多项式都是低次多项式（不超过三次）。1 分段线性插值 y 分段线性插值函数P1 (x) 是一个分段一次多项式（分段线 f (x)性函数）。在几何上就是用折线代替曲线，如图7.2.1，故分段线性插值亦称为折线插值。其插值公式为 P(x) ，xxi , xi +1 (7.2.1) 0 x0 x1 xn-1 xn x 2分段二次插值 图7.2.1 分段线性插值示意图分段二次插值函数P2 (x) 是一个分段二次多项式。在几何上就是分段抛物线代替曲线y = f (x)，故分段二次插值又称为分段抛物插值。其插值公式为， xxi -1 , xi +1 (7.2.2) 3三次Hermite插值三次Hermite插值问题的基本提法一：已知一维数据xx0x1y = f (x)y0y1y = f (x)m0m1求一个三次多项式P3 (x)，使之满足P3 (xi ) = yi，P3 (xi ) = mi，i = 0, 1 (7.2.3)构造三次插值基函数a0(x)，a1(x)，b0(x)，b1(x)，使之满足 (7.2.4)利用这四个插值基函数，取三次多项式P3 (x) 为P3 (x) = a0(x) y0 + a1(x) y1 + b0(x) m0 + b1(x) m1 (7.2.5)将插值条件 (7.2.3) 式代入，可推得： (7.2.6)(7.2.5)、 (7.2.6) 两式构成了三次Hermite插值基本提法一的插值公式。三次Hermite插值问题的基本提法二：已知一维数据xx0x1x2y = f (x)y0y1y2y = f (x)m1求一个三次多项式P3 (x)，使之满足P3 (xi ) = yi，i = 0, 1, 2，P3 (x1 ) = mi (7.2.7)构造三次插值基函数a0(x)，a1(x)，a2 (x)，b1(x)，使之满足 (7.2.8)利用这四个插值基函数，取三次多项式P3 (x) 为P3 (x) = a0(x) y0 + a1(x) y1 + a2(x) y2 + b1(x) m1 (7.2.9)将插值条件 (7.2.7) 式代入，可推得： (7.2.10)(7.2.9)、 (7.2.10) 两式构成了三次Hermite插值基本提法二的插值公式。 7.2.2 一维数据的三次样条插值上述介绍的分段多项式插值，其优点为计算简单、稳定性好、收敛性有保证，且易于在计算机上实现。但它也明显存在着缺陷。它只能保证在每个小区间段 xi , xi +1 内光滑，在各小区间连接点xi 处连续，却不能保证整条曲线的光滑、光顺性，难以满足某些工程的要求。对于象高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即有二阶连续导数。而由60年代开始，首先起源与航空、造船业等工程设计的实际需要而发展起来的样条插值，既保留了分段多项式插值的各种优点，又提高了插值函数的光滑度。在此，仅介绍应用最广且具有二阶连续导数的三次样条插值方法。1 三次样条插值问题的基本提法对于给定的一维数据xx0x1Xny = f (x)y0y1Yn求一个三次多项式S(x) 满足条件 （1）S(xi) = yi，i = 0, 1, , n； （2）S(x) 具有二阶连续导数，特别在节点xi上应满足连续性要求，即对i = 0, 1, , n有2 三次样条插值函数给定区间 a, b 的一个划分D：a = x0 x1 =466&x1(i)714&x1(i)950&x1(i)=1422 a4=1,x1(i),(x1(i)2)/2,(x1(i)3)/6,(x1(i)-x(2)3)/6,(x1(i)-x(3)3)/6,0; y1(i)=a4*alpha; else a5=1,x1(i),(x1(i)2)/2,(x1(i)3)/6,(x1(i)-x(2)3)/6,(x1(i)-x(3)3)/6,(x1(i)-x(4)3)/6; y1(i)=a5*alpha; endendy1 7.3 二维数据的基本插值方法简介对于二维数据的插值，首先要考虑两个问题：一是二维区域是任意区域还是规则区域，二是给定的数据是有规律分布的还是散乱的、随机分布的。第一个问题比较容易处理。目前的插值方法基本上是基于规则区域的，对于不规则区域，只需将其，划分为规则区域或扩充为规则区域来讨论即可。对于第二个问题，当给定的数据是有规律分布时，方法较多也较成熟；而给定的数据是散乱的、随机分布时，没有固定的方法，但一般的处理思想是：从给定的数据出发，依据一定的规律恢复出规则分布点上的数据，转化为数据分布有规律的情形来处理。二维数据插值的方法也有很多。在此，针对给定数据有规律分布和散乱分布两种情形，简单介绍双三次样条插值方法和改进的Shepard方法（反距离平方法）的基本概念和基本思想，及其Matlab实现。 7.3.1 双三次样条插值 双三次样条插值方法，是用来解决规则区域上给定数据有规律分布的插值问题的常用方法。设R：a, bc, d 是xy平面上的一个矩形区域。在x轴和y轴上分别取定分割Dx：a = x0 x1 xnb，Dy：c = y0 y1 ym 0，令由于w(g) 是可微函数，使得如下定义的F(x，y) 在性能上有所改善， 其中。 (7.3.3) 按照上述的思想，可从给定的数据恢复出规则分布点上的数据，接下来就可应用双三次样条插值或其它的二维数据插值方法来处理。 7.3.3 二维数据插值的Matlab实现1规则区域上给定数据有规律分布的二维插值数据形式为：y1y 2y nx1x11z12z1nx 2z21z22z2nx mzm1zm2zmn 插值函数为：interp2( )。其调用格式为zi = interp2(x, y, z, xi, yi, methos),其中 x，y，z 为插值节点，均为向量； zi 为被插值点 (xi, yi) 处的插值结果； methos 为采用的插值方法：nearest：表示最临近插值， linear：表示双线性插值， cubic：表示双三次插值， spline：表示双三次样条插值。注意：上述 methos 中所有的插值方法都要求x和y是单调的网格，x和y可以是等距的也可以是不等距的。2规则区域上给定数据散乱或随机分布的二维插值数据形式为：(x1, y1)(x2, y2)(xn, yn)z1z2zn插值函数为：e01sef和e01sff，。通常两者配合使用，其调用格式为 fnodes, a, rnw, b, c = e01sef(x, y, z)；pf(I, j), ifail = e01sff(x, y, z, rnw, fnodes, px(j), py(i)；其中 x，y，z 为插值节点，均为向量； px(j)，py(i) 为被插值节点； pf(i, j) 为被插值点 (px(j), py(i) 处的插值结果；其它输出参数涉及插值算法，可以不用了解。e01sef的输出fnodes和rnw为确定插值的参数，它们是e01sff需要的输入参数，因此两函数需配合使用。 3例例1 气旋变化情况的可视化表7.3.1是气象学家测量得到的气象资料，它们分别表示在南半球地区按不同纬度、不同月份的平均气旋数字。根据这些数据，绘制出气旋分布曲面图形。 表7.3.1010102020303040405050606070708080901月2.418.720.822.137.348.225.65.30.32月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.323.928.138.24025.36.30.112月3.61625.525.643.441.924.36.60.3解 下面分别用最邻近插值、双线性插值、双三次插值和双三次样条插值，给出不同月份按纬度变化的气旋值（插值结果），并作出可视化图形如下。图7.3.2 四种插值方法的可视化图形最邻近插值的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 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, 0 3.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, 0 0.5, 1.7, 10.1, 32.6, 41.7, 35.4, 22.2, 7.1, 0 0.4, 1.4, 8.3, 33.0, 46.2, 35, 20.2, 5.3, 0.1 0.2, 2.4, 11.2, 31.0, 39.9, 34.7, 21.2, 7.3, 0.2 0.5, 5.8, 12.5, 28.6, 25.9, 35.7, 22.6, 7, 0.3 0.8, 9.2, 21.1, 32.0, 40.3, 39.5, 28.5, 8.6, 0 2.4, 10.3, 23.9, 28.1, 38.2, 40, 25.3, 6.3, 0.1 3.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); figure surf(xi,yi,zi) xlabel(月份),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（单位：英尺）由表7.3.2给出，水深数据是在低潮时测得的。船的吃水深度为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.5y-6.5-81.03.056.5-66.584.0-38.5z9988949解 （1）假设 由题目给出的信息是很少的，除了14个位置的水深之外一无所知。显然，题目要求我们找出水深不到5英尺的区域。为了讨论方便，下面三个假设是合理的： 所给数据是精确的； 讨论区域的海底曲面是光滑的，更确切地说，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形。 水深是一个按区域来划分的变量，在某个位置的水深与其周围区域的水深是相互依赖的，但这种依赖作用随距离的增大而减小。就我们讨论的问题来说，每一个给定数据点影响周围的每一个未知点，一个给定数据点离未知点越近，作用就越大。 （2）问题分析 根据假设，海底曲面是连续光滑的，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形，因而很自然的想法就是用某种光滑的拟合曲面去逼近已知的14个数据点或以14个已知的数据点为基础，利用二维插值补充一些点的水深，以求得水深不超过5米的区域。在此，我们采用二维插值方法，应用Matlab程序，作出矩形区域 (75，200) (-50，150) 范围内的海底地形图、水深不超过5米的危险区域的平面图以及水深不超过5米的危险区域的海底地貌图，并求出水深不超过5米的危险海域范围。 （3）问题求解采用改进的Shepard方法，利用Matlab软件作出作出矩形区域 (75，200) (-50，150) 范围内的海底地形图、水深不超过5米的危险区域的平面图以及水深不超过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,-8,-8,-9,-4,-9;fnodes,minnq,rnw,rnq,ifail=e01sef(x,y,z);nx=100;px=linspace(75,200,nx);图图7.3.3ny=200;py=linspace(-50,150,ny);for i=1:ny for j=1:nx pf(i,j),ifail=e01sff(x,y,z,rnw,fnodes,px(j),py(i); endendsubplot(2,2,2),meshz(px,py,pf+5),title(75，200)x(-50，150) 范