插值与数据拟合.ppt

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

2、来解决，这是一个被称为插值的问题。,解决这个问题，可以通过构造一个与给定数据相适应的函数来解决，这是一个被称为插值的问题。 插值问题的基本提法：对于给定的函数表 其中 f(x) 在区间 a, b 上连续，x0，x1，xn为 a, b 上 n1个互不相同的点，要求在一个性质优良、便于计算的函数类 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) 为插值函数类，称式

3、(7.1.1) 为插值条件，称函数 P(x) 为插值函数，称 f(x) 为被插函数。求插值函数 P(x) 的方法称为插值法。,图 7.1.1 插值问题示意图,引例 7.1.2 在某化学反应中，已知生成物的浓度与时间有关。今测得一组数据如下： 表7.1.2 根据这些数据，我们希望寻找一个 y = f(t) 的近似表达式（如建立浓度y与时间 t 之间的经验公式等）。从几何上看，就是希望根据给定的一组点（1, 4.00），,（16, 10.60），求函数 y = f(t) 的图象的一条拟合曲线。,数据拟合问题的基本提法：对于给定的函数表 其中 f(x) 在区间 a, b 上连续，x0，x1，xn为

4、a, b 上 n1个互不相同的点，要求找一个简单合理的函数近似表达式 (x)，使 (x) 与 f(x) 在某种准则下最为接近，这就是最基本的数据拟合问题（见图7.1.2）。 通常，我们称 (x) 为给定数据点的拟合函数。,图7.1.2 数据拟合问题示意图, 7.1.3 插值与数据拟合的基本理论依据 插值方法与数据拟合的基本理论依据，就是数学分析中的 Weierstrass 定理：设函数 f(x) 在区间 a, b 上连续，则对 0，存在多项式P(x)，使得 即：有界区间上的连续函数被多项式一致逼近。, 7.1.4 实际应用中两种方法的选择 在实际应用中，究竟选择哪种方法比较恰当？总的原则是根据

5、实际问题的特点来决定采用哪一种方法。具体说来，可从以下两方面来考虑： 1如果给定的数据是少量的且被认为是严格精确的，那么宜选择插值方法。采用插值方法可以保证插值函数与被插函数在插值节点处完全相等。,2如果给定的数据是大量的测试或统计的结果，并不是必须严格遵守的，而是起定性地控制作用的，那么宜选用数据拟合的方法。这是因为，一方面测试或统计数据本身往往带有测量误差，如果要求所得的函数与所给数据完全吻合，就会使所求函数保留着原有的测量误差；另一方面，测试或统计数据通常很多，如果采用插值方法，不仅计算麻烦，而且逼近效果往往较差。,7.2 一维数据的基本插值方法简介,插值函数类的取法很多，可以是代数多项

6、式，也可以是三角多项式或有理函数；可以是 a, b 上任意光滑函数，也可以是分段光滑函数。在此介绍最基本、最常用的两种插值方法：分段多项式插值与三次样条插值，及其 Matlab 实现。, 7.2.1 一维数据的分段多项式插值 对于给定的一维数据 分段多项式插值就是求一个分段（共 n 段）多项式 P(x)，使其满足 P(xi) = yi（i = 0, 1, , n）或更高的要求。一般地，分段多项式插值中的多项式都是低次多项式（不超过三次）。,1分段线性插值 分段线性插值函数 P1(x) 是一个分段一次多项式（分段线性函数）。在几何上就是用折线代替曲线，如图 7.2.1，故分段线性插值亦称为折线插

7、值。其插值公式为 其中 xxi, xi +1,图 7.2.1 分段线性插值示意图,2分段二次插值 分段二次插值函数 P2(x) 是一个分段二次多项式。在几何上就是分段抛物线代替曲线 y = f(x)，故分段二次插值又称为分段抛物插值。其插值公式 其中 xxi -1 , xi +1,3三次 Hermite 插值 三次 Hermite 插值问题的基本提法一：已知一维数据 求一个三次多项式 P3(x)，使之满足 P3(xi) = yi，P3(xi) = mi，i = 0, 1 (7.2.3),下面的 (7.2.5)、 (7.2.6) 两式构成里三次 Hermite插值基本提法一的插值公式 P3(x)

8、 = 0(x)y0 1(x)y1 0(x)m0 1(x)m1 (7.2.5),三次 Hermite 插值问题的基本提法二：已知一维数据 求一个三次多项式 P3(x)，使之满足 P3(xi) = yi，i = 0, 1, 2，P3(x1) = mi (7.2.3),下面的 (7.2.9)、 (7.2.10) 两式构成里三次Hermite 插值基本提法二的插值公式 P3(x) = 0(x)y0 1(x)y1 0(x)m0 1(x)m1 (7.2.9), 7.2.2 一维数据的三次样条插值 上述介绍的分段多项式插值，其优点为计算简单、稳定性好、收敛性有保证，且易于在计算机上实现。但它也明显存在着缺陷

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

10、二阶连续导数，特别在节点 xi 上应满足连续性要求，即对 i = 0, 1, , n 有,2三次样条插值函数 给定区间 a, b 的一个划分：a = x0 x1 xn = b，设函数 y = f(x) 在节点 xi 上的值为 yi = f(xi)，i = 0, 1, , n。如果 S(x) 于 a, b 有二阶连续导数，且在每个小区间 xi , xi+1 上是三次多项式，则称 S(x) 是节点 x0，x1，xn上的三次样条函数。如果 S(x) 在节点 xi上还满足插值条件 S(xi) = yi，i = 0, 1, , n， (7.2.11) 则称 S(x) 为三次样条插值函数。,对应于划分 的

11、三次样条插值函数的表达式为 其中,3边界条件 在式 (7.2.12) 给出的三次多项式中，共含有 n3 个待定系数。而由插值条件 (7.2.11) 式，可列出 n1个方程，方程组中未知数的个数比方程个数多 2，还需附加 2 个条件才能进行求解。通常可在区间端点 x0 = a 和 xn = b 处各附加一个条件（称为边界条件或边值条件）去确定 S(x)。,边界条件类型很多，较基本而又常见的有三类： (1) 第一边值条件，即给出边界点的一阶导数值 S (x0) = y0， S(xn) = yn (7.2.13) (2) 第二边值条件，即给出边界点的二阶导数值 S (x0) = y0， S(xn)

12、= yn (7.2.14) 特别地，当 S(x0) = S(xn) = 0 时，称为自然边界条件。满足自然边界条件的三次样条插值函数称为自然样条插值函数。,(3) 第三边值条件（混合边值条件） 其中 1、2、1、2、1、2 为定数。当 1、2 为零时，则为第一边值条件，当 1、2 为零时，则为第二边值条件。, 7.2.3 一维数据插值的 Matlab 实现 Matlab 中一维数据的插值函数为：interp1( )。其调用格式为 yi = interp1(x, y, xi, methos), 其中，x，y 为插值节点，均为向量； xi 任取的被插值点可以是一个数值，也可以是 一个向量； yi

13、为被插值点xi处的插值结果； methos 为采用的插值方法： nearest：表示最临近插值， linear：表示线性插值， cubic：表示三次插值， spline：表示三次样条插值。,注意： （1）上述 methos 中所有的插值方法都要求 x 是单调的，并且 xi 不能超过 x 的取值范围； （2）三次样条插值函数的调用格式有两种： yi = interp1(x, y, xi, spline) yi = spline(x, y, xi) 它们是等价的。,注意： （1）上述 methos 中所有的插值方法都要求 x 是单调的，并且 xi 不能超过 x 的取值范围； （2）三次样条插值函数

14、的调用格式有两种： yi = interp1(x, y, xi, spline) yi = spline(x, y, xi) 它们是等价的。 例 7.2.1 引例7.1.1的求解（见教材）,7.3 二维数据的基本插值 方法简介,对于二维数据的插值，首先要考虑两个问题：一是二维区域是任意区域还是规则区域，二是给定的数据是有规律分布的还是散乱的、随机分布的。,第一个问题比较容易处理。目前的插值方法基本上是基于规则区域的，对于不规则区域，只需将其，划分为规则区域或扩充为规则区域来讨论即可。对于第二个问题，当给定的数据是有规律分布时，方法较多也较成熟；而给定的数据是散乱的、随机分布时，没有固定的方法，

15、但一般的处理思想是：从给定的数据出发，依据一定的规律恢复出规则分布点上的数据，转化为数据分布有规律的情形来处理。,二维数据插值的方法也有很多。在此，针对给定数据有规律分布和散乱分布两种情形，简单介绍双三次样条插值方法和改进的 Shepard方法（反距离平方法）的基本概念和基本思想，及其 Matlab 实现。, 7.3.1 双三次样条插值 双三次样条插值方法，是用来解决规则区域上给定数据有规律分布的插值问题的常用方法。 实际上，双三次样条函数是由两个一维三次样条函数作直积产生的。对任意固定的 y0c, d，S(x, y0) 是关于 x 的三次样条函数，同理，对任意固定的 x0a, b，S(x0,

16、 y) 是关于 y 的三次样条函数。从而，根据一维三次样条函数的算法可以设计出 S(x, y) 的具体算法。, 7.3.2 改进的 Shepard 方法 改进的 Shepard 方法，也称反距离加权平均法，这是解决规则区域上给定数据散乱、随机分布的插值问题的一个常用的方法。 问题：设 T = a, bc, d 上散乱分布 N 个点 V1，V2，VN ，其中 Vk = (xk，yk) 处给出数据 fk ，k = 1, 2, , N。要寻求 T 上的二元函数 F(x，y)，使 F(xk，yk) = fk ，k = 1, 2, , N。,一个典型的容易想到的方法是“反距离加权平均”方法，又称之为 S

17、hepard 方法。这方法的基本思想是，在非给定数据的点处，定义其函数值由已知数据点与该点距离的近或远作加权平均决定。 按照上述的思想，可从给定的数据恢复出规则分布点上的数据，接下来就可应用双三次样条插值或其它的二维数据插值方法来处理。, 7.3.3 二维数据插值的 Matlab 实现 1规则区域上给定数据有规律分布的二维插值 数据形式为：,插值函数为：interp2( )。其调用格式为 zi = interp2(x, y, z, xi, yi, methos), 其中 x，y，z 为插值节点，均为向量； zi 为被插值点 (xi, yi) 处的插值结果； methos 为采用的插值方法： n

18、earest：表示最临近插值， linear：表示双线性插值， cubic：表示双三次插值， spline：表示双三次样条插值。 注意：上述 methos 中所有的插值方法都要求 x 和 y是单调的网格，x 和 y 可以是等距的也可以是不等距的。,2规则区域上给定数据散乱或随机分布的二维插值 数据形式为：,插值函数为：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

19、 为插值节点，均为向量； px(j)，py(i) 为被插值节点； pf(i, j) 为被插值点 (px(j), py(i) 处的插值结果； 它输出参数涉及插值算法，可以不用了解。e01sef 的输出 fnodes 和 rnw 为确定插值的参数，它们是 e01sff 需要的输入参数，因此两函数需配合使用。,3实例 例 7.3.1 气旋变化情况的可视化 表 7.3.1 是气象学家测量得到的气象资料，它们分别表示在南半球地区按不同纬度、不同月份的平均气旋数字。根据这些数据，绘制出气旋分布曲面图形。,解：下面分别用最邻近插值、双线性插值、双三次插值和双三次样条插值，给出不同月份按纬度变化的气旋值（插值

20、结果），并作出可视化图形如下。,7.4 一维数据拟合的最小二乘法简介,数据拟合问题的形式多种多样，解决的方法也有许多。在此，我们只简单介绍以最小二乘法为准则的一维数据拟合方法。, 7.4.1 数据拟合的最小二乘法简介 对于给定的一组测量数据 设y = (x) 为一拟合函数，记 i = (xi) yi (i = 1, 2, , n)， (7.4.1) 则称 i 为拟合函数 (x) 在 xi 点处的偏差或残量。,为使 (x) 在整体上尽可能与给定数据的函数 f(x) 近似，我们常采用偏差的平方和达到最小，即 来保证每个偏差的绝对值 |i| 都很小。这一原则称为最小二乘原则，根据最小二乘原则确定拟合

21、函数 (x) 的方法称为最小二乘法。,1. 线性最小二乘拟合 我们知道，函数系 xk | k = 0, 1, , m 的线性组合 (x) = a0 a1x a2x2 amxm 为 m 次多项式。一般地，若函数系 k(x) | k = 0, 1, , m 是线性无关的，则其线性组合 称为函数系 k(x) | k = 0, 1, , m 的广义多项式。如三角多项式 就是函数系 1, cosx, sinx, cos2x, sin2x, , cosmx, sinmx 的广义多项式。,设 k(x) | k = 0, 1, , m 为一线性无关的函数系，取拟合函数为 (7.4.3) 式给出的广义多项式，使

22、得 (7.4.2) 成立。由于 (x) 的待定系数a0, a1, a2, , am 全部以线性形式出现，故我们称之为线性最小二乘拟合。 在式 (7.4.2) 中，目标函数 S 是关于参数 a0, a1, a2 , , am 的多元函数，由多元函数取得最小值的必要条件知，欲使 S 达到极小，须满足,即 亦即 其中 k = 0, 1, , m，式 (7.4.4) 是关于 a0, a1, a2 , , am 的线性方程组，称为正规方程组。 从正规方程组 (7.4.4) 中解出 a0, a1, a2 , , am，于是就得到了最小二乘拟合函数 (x)。,2非线性最小二乘拟合 如果拟合函数 (x) =

23、(x, a0, a1, a2 , , am) 的待定参数 a0, a1, a2 , , am 不能全部以线性形式出现，如指数拟合函数 等，这便是非线性最小二乘拟合问题。 一般地，非线性最小二乘拟合问题是一个非线性函数的极小化问题，可用非线性优化方法求解。,3最小二乘拟合函数的选择 最小二乘法中，拟合函数的选择是很重要的。可以通过对给定数据的分析来选择，也可以直接由实际问题给定。最常用的是多项式和样条函数，尤其是当不知道该选择什么样的拟合函数时，通常可以考虑选择样条函数。 另外，对同一问题，也可选择不同的函数进行最小二乘拟合，比较各自误差的大小，从中选出误差较小的作为拟合函数。, 7.4.2 一

24、维数据拟合的 Matlab 实现 1. 多项式函数拟合 拟合函数的命令为：polyfit( )，其调用格式为：a = polyfit(xdata, ydata, m)，其中 m 表示多项式的最高阶数； xdata，ydata 为将要拟合的数据，它们都 是以数组方式输入； a 输出参数，为拟合多项式的系数 a = a0, a1, a2 , , am。 多项式在 x 处的值 y 可用如下命令格式计算： y = polyval(a, x)。,2. 一般的曲线拟合 拟合函数的命令为：curvefit( )，或lsqcurvefit( )，其调用格式为 p = curvefit(Fun, p0, xda

25、ta, ydata)，或 p = lsqcurvefit(Fun, p0, xdata, ydata)， 其中 Fun 表示函数Fun (p, xdta) 的M文件； P0 为函数的初值。 若要求点 x 处的函数值 y，可用程序 f = Fun(p,x) 计算。,3实例 例 7.4.1 求解7.1.2 引例7.1.2。 解：(1) 选择拟合曲线 作出给定数据的散点图如下：,通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。具有这种特性的曲线很多，我们选取如下三种函数来拟合： 多项式 (x) = a0 a1x amxm， m 为适当选取的正整数； ； （a 0, b 0）。,(2)

26、拟合运算 首先，分别用二、三、六次多项式拟合，计算得输出参数分别为 p1 = 0.0445，1.0711，4.3252 p2 = 0.0060，0.1963，2.1346，2.5952 p3 = 0.0000，0.0004，0.0103，0.1449， 1.1395，4.9604，0.0498 拟合函数分别为 (1)(x) = 0.0445 1.0711x 4.3252x2 (2)(x) = 0.0060 0.1963x 2.1346x2 2.5952x3 (3)(x) = 0.0004x 0.0103x2 0.1449x3 1.1395x4 4.9304x5 0.0498x6；,其次，再用有

27、理分式 拟合，计算得输出参数分别为 p = 0.0841，0.1392 拟合函数为,最后，用指数函数 拟合，计算得输出参数分别为 p = 11.3578，1.0873 拟合函数为 三种方式五个种函数的拟合曲线见图 7.4.27.4.4。,图 7.4.2 多项式函数拟合曲线图,图7.4.3 有理分式函数拟合曲线图,图7.4.4 指数函数拟合曲线图,(3) 误差分析 和给定的 16 组数据比较，三种方式五个函数拟合的误差见下表： 表7.4.1 五个函数拟合的误差表 其中：偏差平方和以及平均偏差平方和为：,从上表中求得的误差情况来看，好似六次多项式函数拟合的最好，指数函数拟合次之，然后分别是有理分式

28、函数拟合、三次多项式函数拟合和二次多项式函数拟合。但是，就这个实际问题的本质来说，化学反应中生成物的浓度到一定时间后应基本稳定，即当 t 时， f(t) 常数。 而我们有,三个多项式函数拟合曲线的趋势也可从图 7.4.5 看出。 图 7.4.5,(4) 结论 通过以上的计算和分析，我们得出如下结论：本问题可用指数函数或有理分式函数来拟合，其拟合函数分别为,拟合误差为 表7.4.2 (5) Matlab 程序 略。,7.5 范例,水道测量数据（AMCM 86A题）。 本问题由加州海军研究生院数学系的Richard Franke 提供，问题如下： 在某海域测得一些点（x, y）处的水深 z（单位：英尺）由表 7.5.1 给出，水深数据是在低潮时测得的。船的吃水深度为 5 英尺，问在矩形区域 (75，200)(50，150) 里的哪些地方船要避免进入。,表7.5.1 水道水深测量数据（单位：英尺）,解：