数值计算04-插值与拟合教材

1数值计算第四章插值与拟合

我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算法。此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。

在实际中，常常要处理由实验或测量所得到的一些离散数据。例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值．%interp_1hours=1:12;temps=[589152529313022252724];plot(hours,temps,'rs‘)

在实际中，常常要处理由实验或测量所得到的一些离散数据。例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值．

在实际中，常常要处理由实验或测量所得到的一些离散数据。例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值．

在实际中，常常要处理由实验或测量所得到的一些离散数据。例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值．插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为插值问题。（不需要函数表达式）

如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近似函数的方法称为数据拟合。（必须有函数表达式）近似函数不一定（曲线或曲面）通过所有的数据点。9第四章插值与拟合4.1插值的含义4.2插值的Matlab实现4.3拟合问题104.1插值的含义一维插值:已知函数在n+1个节点的值,求任一点的函数值114.1插值的含义一维插值:已知函数在n+1个节点的值,求任一点的函数值二维插值:已知函数在个节点的值,求任一点的函数值一维插值的定义已知n+1个节点其中互不相同，不妨设求任一插值点处的插值

节点可视为由产生，表达式复杂，或无封闭形式，或未知。

构造一个(相对简单的)函数通过全部节点,即再用计算插值，即

2.2插值的Matlab实现基本格式：yc=interp1(x,y,cx,’method’)%x,y分别表示已知数据点的横、纵坐标向量，x必须单调；%cx为需要插值的横坐标数据%method为插值方法，有‘nearest’最近邻点插值‘linear’线性插值（默认）‘spline’三次样条插值‘cubic’三次插值注：Lagrange,Newton插值法需自编程序。线性插值主要参数

15参数名称说明特点nearest邻近点插值法。根据已知两点间的插值点与这两点之间的位置远近插值。当插值点距离前点近时，取前点的值，否则取后点的值速度最快，但平滑性差linear线性插值。把相邻的数据点用直线连接，按所生成的曲线进行插值，是默认的插值方法占有的内存较邻近点插值方法多，运算时间也稍长，与邻近点插值不同，其结果是连续的，但在顶点处的斜率会改变spline三次样条插值。用已知数据求出样条函数后，按照样条函数插值运算时间长，但内存的占有较立方插值方法要少，三次样条插值的平滑性很好，但如果输入的数据不一致或数据点过近，可能出现很差的插值结果cubic立方插值法，也称三次多项式插值。用已知数据构造出三次多项式进行插值需要较多的内存和运算时间，平滑性很好bicubic双立方插值法。利用已知的数据点拟合一个双立方曲面，然后根据插值点的坐标插值，每个插值点的值由该点附近的六个点的坐标确定二维插值函数独有。插值点处的值和该点值的导数都连续例2.9：对函数

进行11个点的三次样条插值：x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];y=[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26];xi=-5:0.01:5;yi=interp1(x,y,xi);plot(xi,yi,x,y,'o')holdonf=inline('1/(x^2+1)');fplot(f,[-5,5],'r')

例2.9：对函数

进行11个点的三次样条插值：x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];y=[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26];xi=-5:0.01:5;yi=interp1(x,y,xi,'spline');plot(xi,yi,x,y,'o')holdonf=inline('1/(x^2+1)');fplot(f,[-5,5],'r')

例2.9：对函数

进行11个点的三次样条插值：x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];y=[1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26];xi=-5:0.01:5;yi=interp1(x,y,xi,‘cubic');plot(xi,yi,x,y,'o')holdonf=inline('1/(x^2+1)');fplot(f,[-5,5],'r')

二维插值的定义

xyO第一种（网格节点）：

已知m

n个节点其中互不相同，不妨设

构造一个二元函数通过全部已知节点,即再用计算插值，即第二种（散乱节点）：

yx0已知n个节点其中互不相同，

构造一个二元函数通过全部已知节点,即再用计算插值，即

注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

将四个插值点（矩形的四个顶点）处的函数值依次简记为：

分片线性插值xy

(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足

双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O要求x0,y0单调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。z=interp2(x0,y0,z0,x,y,’method’)被插值点插值方法用MATLAB作网格节点数据的插值插值节点被插值点的函数值‘nearest’

最邻近插值‘linear’

双线性插值‘cubic’

双三次插值缺省时,双线性插值例2.10：测得平板表面3*5网格点处的温度分别为：828180828479636165818484828586试作出平板表面的温度分布曲面z=f(x,y)的图形。输入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.%输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)2．以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值,画出插值后的温度分布曲面图.32已知某处山区地形选点测量坐标数据为：x=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5y=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6海拔高度数据为：z=8990878592919693908782

9296989995918986848284

9698959290888584838185

8081828995969392898686

8285879899969788858283

8285899495939291868488

8892939495898786838192

9296979896939584828184

8585818280808185909395

8486819899989796958487

8081858283848790958688

8082818485868382818082

8788899899979698949287例2.11画出该山区地形地貌图.33插值前

插值后Tomountain.m最近邻点插值双线性插值双三次插值

插值函数griddata格式为:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散点数据的插值计算

要求cx取行向量，cy取为列向量．被插值点插值方法插值节点被插值点的函数值‘nearest’最邻近插值‘linear’

双线性插值‘cubic’

双三次插值'v4'-MATLAB提供的插值方法缺省时,双线性插值例2.13:在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）×（-50，150）里的哪些地方船要避免进入．4.作出水深小于5的海域范围,即z=5的等高线.2.在矩形区域(75,200)×(-50,150)进行插值。1.输入插值基点数据3.作海底曲面图步骤：%程序一：插值并作海底曲面图x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5];z=[48686889988949];x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,'v4');meshc(x1,y1,z1)海底曲面图%程序二：插值并作出水深小于5的海域范围。x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,'v4');%插值z1(z1>=5)=nan;%将水深大于5的置为nan，这样绘图就不会显示出来meshc(x1,y1,z1)水深小于5的海域范围

在实际应用中，经常遇到下列数据处理问题：已知函数在m个点上的数据表，寻求其近似函数。

2.3曲线拟合问题所谓”曲线拟合”,是指根据给定的数据表,寻找一个简单的表达式来”拟合”该组数据,此处的”拟合”的含义为:不要求该表达式对应的近似曲线完全通过所有的数据点,只要求该近似曲线能够反映数据的基本变化趋势.

已知函数在m个点上的数据表，寻求其近似函数。设的近似函数为其中是某函数族中的已知线性无关函数。(1)最小二乘多项式拟合称为残差向量寻求一组常数，要求的2-范数达到最小。如果m=n，且以及即多项式插值记则得到最小二乘问题：上述问题的解也称为方程组的最小二乘解。当时称之为超定（或矛盾）方程组。例2.14：考察某种纤维的强度与其拉伸倍数的关系.下表是实际测定的24个纤维样品的强度与相应的拉伸倍数的数据记录:编号拉伸倍数强度编号拉伸倍数强度11.91.41355.5221.3145.2532.11.81565.542.52.5166.36.452.72.8176.5662.72.5187.15.373.531986.583.52.72087944218.98.51043.52298114.54.2239.58.1124.63.524108.1可以看出,纤维强度随拉伸倍数增加而增加并且24个点大致分布在一条直线附近该直线称为这一问题的数学模型。因此可认为强度与拉伸倍数之间的主要关系是线性关系怎样确定a,b,使得直线能较好地反映所给数据的基本“变化趋势”?采用最小二乘的思想令问题转化为求参数使

达到最小值。这种求线性函数y=a+bx的过程称为线性拟合。(2)非线性拟合已知函数在若干个点上的数据表，确定参数和利用经验函数拟合某组数据：某些非线性拟合问题可转化为线性拟合问题线性化处理：令则

由线性拟合方法可得到和，从而得到和。如何选取拟合函数？1.通过机理分析建立数学模型来确定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.将数据(xi,yi)i=1,…n作图，通过直观判断确定f(x)：用MATLAB解拟合问题1、线性最小二乘拟合2、非线性最小二乘拟合用MATLAB作线性最小二乘拟合1.作多项式f(x)=a1xm+…+amx+am+1拟合,可利用已有程序:a=polyfit(x,y,m)2.对超定方程组可得最小二乘意义下的解。，用3.多项式在x处的值y可用以下命令计算：

y=polyval（a，x）输出拟合多项式系数a=[a1,…am,

am+1](数组)）输入同长度的数组X，Y拟合多项式次数即要求出二次多项式:中的使得:例2.15对下面一组数据作二次多项式拟合1）输入以下命令：x=0:0.1:1;y=[-04471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(10,1)]；

A=R\y'解法1．用解超定方程的方法2）计算结果：Ａ=-9.810820.1293-0.03171）输入以下命令：

x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出数据点和拟合曲线的图形2）计算结果：Ａ=-9.810820.1293-0.0317解法2．用多项式拟合的命令1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非线性最小二乘拟合Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的.

lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的参变量x(向量),使得

输入格式为:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一个事先建立的定义函数F(x,xdata)

的M-文件,自变量为x和xdata说明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知数据点选项见无约束优化

lsqnonlin用以求含参量x（向量）的向量值函数

f(x)=(f1(x),f2(x),…,fn(x))T

中的参量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知数据点：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）输入格式为：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin

（‘fun’，x0，options）；

3）x=lsqnonlin

（‘fun’，x0，options，‘grad’）；

4）[x，options]=lsqnonlin

（‘fun’，x0，…）；

5）[x，options，funval]=lsqnonlin

（‘fun’，x0，…）；说明：x=lsqnonlin

（‘fun’，x0，options）；fun是一个事先建立的定义函数f(x)的M-文件，自变量为x迭代初值选项见无约束优化

例2.16用下面一组数据拟合中的参数a，b，k该问题即解最优化问题：

1）编写M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）输入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];

x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a，b，k)解法1.用命令lsqcurvefit3）运算结果为：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）结论：a=0.0063,b=-0.0034,k=0.2542

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）编写M-文件curvefun2.m

functionf=curvefun2(x)tdata=100:100:1000;cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata2）输入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdatatdata的值写在curvefun2.m中3）运算结果为

f=1.0e-003*(0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792x=0.0063-0.00340.2542可以看出,两个命令的计算结果是相同的.4）结论：即拟合得a=0.0063b=-0.0034k=0.2542如何预报人口的增长人口的增长是当前世界上引起普遍关注的问题，并且我们会发现在不同的刊物预报同一时间的人口数字不相同，这显然是由于用了不同的人口模型计算的结果。我国是世界第一人口大国，基本上地球每九个人中就有一个中国人。有效地控制我国人口的增长是使我过全面进入小康社会、到21世纪中叶建成富强民主文明的社会主义国家的需要。而有效控制人口增长的前提是要认识人口数量的变化规律，建立人口模型，作出较准确的预报。

例:如何预报人口的增长例2.17：1949年—1994年我国人口数据资料如下：

年份xi1949195419591964196919741979198419891994人口数yi5.46.06.77.08.19.19.810.311.311.8建模分析我国人口增长的规律,预报1999年我国人口数。模型一：假设人口随时间线性地增加模型：参数估计观测值的模型：

拟合的精度:误差平方和。可以算出：a=-283.2320b=0.1480模型：y=–1.93+0.146x

则可看成是线性方程,用polyfit命令计算得：模型二：指数增长模型可变为YA=+BXa=2.33,b=0.0179则所求模型为:%程序nihe_2.m如下：x=[1949195419591964196919741979198419891994];y=[5.46.06.77.08.19.19.810.311.311.8];a=polyfit(x,y,1);x1=[1949:10:1994];y1=a(2)+a(1)*x1;b=polyfit(x,log(y),1);y2=exp(b(2))*exp(b(1)*x1);plot(x,y,'*')holdonplot(x1,y1,'--r')holdonplot(x1,y2,'-k')legend('原曲线','模型一曲线','模型二曲线')77例3：体重约70kg的某人在短时间内喝下2瓶啤酒后，隔一定时间测量他的血液中酒精含量（mg/100ml），得到数据如表所示。试用所给数据用函数

进行拟合，并求出未知常数。时间t/h0.250.50.7511.522.533.544.55酒精含量hmg/100ml306875828277686858515041时间t/h678910111213141516酒精含量hmg/100ml3835282518151210774

78分析：由于拟合函数形式已经确定，且不是多项式函数，但若取对数可得这样对参数

是线性的。因此可考虑先对数据h进行对数变换，再调用lsqcurvefit()函数进行拟合。79clearall;closeall;t=[0.250.50.7511.522.533.544.55678910111213141516];h=[3068758282776868585150413835282518151210774];h1=log(h);%对数变换f=inline(‘a(1)+a(2).*log(t)+a(3).*t’,‘a’,‘t’);[x,r]=lsqcurvefit(f,[1,0.5,-0.5],t,h1)%求参数lna,b,c的拟合值x=4.48340.4709-0.2663r=0.4097输出结果：80

近年来我国的电信事业发展迅速，现已成世界第一电信大国。据统计某市的在过去近９年中通信工具的拥有量如下表（单位为万台）：

应用案例1：通信工具的发展趋势

问题的提出：使用matlab语句plot(x0,y0,’r.’,’Markersize’,10)81Oxy..........246810108642.........

应用案例