校内培训插值拟合与

插值与拟合§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米…）处的水温。解决这个问题，可以通过构造一个与给定数据相适应的函数来解决，这是一个被称为插值的问题。深度x(m)46671495014221634水温y(C)7.044.283.402.542.13解决这个问题，可以通过构造一个与给定数据相适应的函数来解决，这是一个被称为插值的问题。

插值问题的基本提法：对于给定的函数表其中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）。xx0x1…xny=f(x)y0y1…yn为便于叙述，通常称区间[a,b]为插值区间，称点x0，x1，…，xn为插值节点，称函数类{P(x)}为插值函数类，称式(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)的图象的一条拟合曲线。时间t（分）12345678浓度y1034.006.408.008.809.229.509.709.86时间t（分）910111213141516浓度y10310.0010.2010.3210.3210.5010.5510.5810.60

数据拟合问题的基本提法：对于给定的函数表其中f(x)在区间[a,b]上连续，x0，x1，…，xn为[a,b]上n1个互不相同的点，要求找一个简单合理的函数近似表达式

(x)，使

(x)与f(x)在某种准则下最为接近，这就是最基本的数据拟合问题（见图7.1.2）。通常，我们称

(x)为给定数据点的拟合函数。xx0x1…xny=f(x)y0y1…yn图7.1.2数据拟合问题示意图

§7.1.3插值与数据拟合的基本理论依据插值方法与数据拟合的基本理论依据，就是数学分析中的Weierstrass定理：设函数f(x)在区间[a,b]上连续，则对

>0，存在多项式P(x)，使得即：有界区间上的连续函数被多项式一致逼近。§7.1.4实际应用中两种方法的选择在实际应用中，究竟选择哪种方法比较恰当？总的原则是根据实际问题的特点来决定采用哪一种方法。具体说来，可从以下两方面来考虑：1．如果给定的数据是少量的且被认为是严格精确的，那么宜选择插值方法。采用插值方法可以保证插值函数与被插函数在插值节点处完全相等。2．如果给定的数据是大量的测试或统计的结果，并不是必须严格遵守的，而是起定性地控制作用的，那么宜选用数据拟合的方法。这是因为，一方面测试或统计数据本身往往带有测量误差，如果要求所得的函数与所给数据完全吻合，就会使所求函数保留着原有的测量误差；另一方面，测试或统计数据通常很多，如果采用插值方法，不仅计算麻烦，而且逼近效果往往较差。§7.2一维数据的基本插值方法简介

插值函数类的取法很多，可以是代数多项式，也可以是三角多项式或有理函数；可以是[a,b]上任意光滑函数，也可以是分段光滑函数。在此介绍最基本、最常用的两种插值方法：分段多项式插值与三次样条插值，及其Matlab实现。§7.2.1一维数据的分段多项式插值对于给定的一维数据分段多项式插值就是求一个分段（共n段）多项式P(x)，使其满足P(xi)=yi（i=0,1,…,n）或更高的要求。一般地，分段多项式插值中的多项式都是低次多项式（不超过三次）。xx0x1…xny=f(x)y0y1…yn1．分段线性插值

分段线性插值函数P1(x)是一个分段一次多项式（分段线性函数）。在几何上就是用折线代替曲线，如图7.2.1，故分段线性插值亦称为折线插值。其插值公式为

其中x[xi,xi+1]

图7.2.1分段线性插值示意图

2．分段二次插值分段二次插值函数P2(x)是一个分段二次多项式。在几何上就是分段抛物线代替曲线y=f(x)，故分段二次插值又称为分段抛物插值。其插值公式其中x[xi-1

,xi+1]

3．三次Hermite插值三次Hermite插值问题的基本提法一：已知一维数据求一个三次多项式P3(x)，使之满足P3(xi)=yi，P3(xi)=mi，i=0,1(7.2.3)xx0x1y=f(x)y0y1y

=f(x)m0m1下面的(7.2.5)、(7.2.6)两式构成里三次Hermite插值基本提法一的插值公式P3(x)=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)xx0x1x2y=f(x)y0y1y2y

=f(x)m1下面的(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一维数据的三次样条插值上述介绍的分段多项式插值，其优点为计算简单、稳定性好、收敛性有保证，且易于在计算机上实现。但它也明显存在着缺陷。它只能保证在每个小区间段[xi,xi+1]内光滑，在各小区间连接点xi

处连续，却不能保证整条曲线的光滑、光顺性，难以满足某些工程的要求。对于象高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即有二阶连续导数。而由60年代开始，首先起源与航空、造船业等工程设计的实际需要而发展起来的样条插值，既保留了分段多项式插值的各种优点，又提高了插值函数的光滑度。在此，仅介绍应用最广且具有二阶连续导数的三次样条插值方法。

1．三次样条插值问题的基本提法对于给定的一维数据求一个三次多项式S(x)满足条件（1）S(xi)=yi，i=0,1,…,n；（2）S(x)具有二阶连续导数，特别在节点xi

上应满足连续性要求，即对i=0,1,…,n有xx0x1…xny=f(x)y0y1…yn

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)为三次样条插值函数。对应于划分的三次样条插值函数的表达式为其中

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)=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——为被插值点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）三次样条插值函数的调用格式有两种：yi=interp1(x,y,xi,‘spline’)yi=spline(x,y,xi)它们是等价的。§7.3二维数据的基本插值方法简介

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

二维数据插值的方法也有很多。在此，针对给定数据有规律分布和散乱分布两种情形，简单介绍双三次样条插值方法和改进的Shepard方法（反距离平方法）的基本概念和基本思想，及其Matlab实现。§7.3.1双三次样条插值双三次样条插值方法，是用来解决规则区域上给定数据有规律分布的插值问题的常用方法。实际上，双三次样条函数是由两个一维三次样条函数作直积产生的。对任意固定的y0[c,d]，S(x,y0)是关于x的三次样条函数，同理，对任意固定的x0[a,b]，S(x0,y)是关于y的三次样条函数。从而，根据一维三次样条函数的算法可以设计出S(x,y)的具体算法。§7.3.2改进的Shepard方法改进的Shepard方法，也称反距离加权平均法，这是解决规则区域上给定数据散乱、随机分布的插值问题的一个常用的方法。问题：设T=[a,b][c,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。一个典型的容易想到的方法是“反距离加权平均”方法，又称之为Shepard方法。这方法的基本思想是，在非给定数据的点处，定义其函数值由已知数据点与该点距离的近或远作加权平均决定。按照上述的思想，可从给定的数据恢复出规则分布点上的数据，接下来就可应用双三次样条插值或其它的二维数据插值方法来处理。§7.3.3二维数据插值的Matlab实现1．规则区域上给定数据有规律分布的二维插值数据形式为：y1y2…ynx1x11z12…z1nx2z21z22…z2n……………xmzm1zm2…zmn

插值函数为：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)z1z2…zn

插值函数为：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．实例例7.3.1气旋变化情况的可视化表7.3.1是气象学家测量得到的气象资料，它们分别表示在南半球地区按不同纬度、不同月份的平均气旋数字。根据这些数据，绘制出气旋分布曲面图形。0—1010—2020—3030—4040—5050—6060—7070—8080—901月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.4一维数据拟合的最小二乘法简介

数据拟合问题的形式多种多样，解决的方法也有许多。在此，我们只简单介绍以最小二乘法为准则的一维数据拟合方法。§7.4.1数据拟合的最小二乘法简介对于给定的一组测量数据设y=

(x)为一拟合函数，记i

=(xi)

yi(i=1,2,…,n)，(7.4.1)则称i

为拟合函数(x)在xi

点处的偏差或残量。xx1x2…xny=f(x)y1y2…yn为使(x)在整体上尽可能与给定数据的函数f(x)近似，我们常采用偏差的平方和达到最小，即来保证每个偏差的绝对值|i|都很小。这一原则称为最小二乘原则，根据最小二乘原则确定拟合函数(x)的方法称为最小二乘法。

1.线性最小二乘拟合我们知道，函数系{xk|k=0,1,…,m}的线性组合(x)=a0

a1xa2x2

…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)式给出的广义多项式，使得(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)=(x,a0,a1,a2,…,am)的待定参数a0,a1,a2,…,am不能全部以线性形式出现，如指数拟合函数等，这便是非线性最小二乘拟合问题。一般地，非线性最小二乘拟合问题是一个非线性函数的极小化问题，可用非线性优化方法求解。

3．最小二乘拟合函数的选择最小二乘法中，拟合函数的选择是很重要的。可以通过对给定数据的分析来选择，也可以直接由实际问题给定。最常用的是多项式和样条函数，尤其是当不知道该选择什么样的拟合函数时，通常可以考虑选择样条函数。另外，对同一问题，也可选择不同的函数进行最小二乘拟合，比较各自误差的大小，从中选出误差较小的作为拟合函数。§7.4.2一维数据拟合的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,xdata,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)拟合运算首先，分别用二、三、六次多项式拟合，计算得输出参数分别为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.04451.0711x

4.3252x2(2)(x)=0.00600.1963x

2.1346x2

2.5952x3(3)(x)=0.0004x

0.0103x2

0.1449x3

1.1395x4

4.9304