插值、拟合与MATLAB编程

1、插值、拟合与MATLAB编程相关知识在生产和科学实验中，自变量x与因变量y间的函数关系y f（x）有时不能写 出解析表达式，而只能得到函数在若干点的函数值或导数值，或者表达式过于 复杂而需要较大的计算量。当要求知道其它点的函数值时，需要估计函数值在 该点的值。为了完成这样的任务，需要构造一个比较简单的函数y （X），使函数在观测 点的值等于已知的值，或使函数在该点的导数值等于已知的值，寻找这样的函 数y （X）有很多方法。根据测量数据的类型有以下两类处理观测数据的方法。（1）测量值是准确的，没有误差，一般用插值。（2）测量值与真实值有误差，一般用曲线拟合。在MATLAB中，无论是插值还是拟合，

2、都有相应的函数来处理。一、插值1、一维插值：已知离散点上的数据集（Xi,yi）,（X2,y2）,（Xn,yn），即已知在点集 X二Xl，X2，,Xn上的函数值 Y二yi，y2,，yn，构造一个解析函数（其图形为一 曲线）通过这些点，并能够求出这些点之间的值，这一过程称为一维插值。MATLAB 命令：yi=interp1（X, Y, Xi, method）该命令用指定的算法找出一个一元函数 y f（x），然后以f（x）给出x处的值。 xi可以是一个标量，也可以是一个向量，是向量时，必须单调，method可以下列方法之一： nearest最近邻点插值，直接完成计算； spline三次样条函数插值；

3、linear线性插值(缺省方式)，直接完成计算； cubic '次函数插值；对于minxi,maxxi外的值，MATLAB使用外推的方法计算数值。例1:已知某产品从1900年到2010年每隔10年的产量为：75.995, 91.972, 105.711, 123.203, 131.699, 150.697, 179.323, 203.212, 226.505, 249.633,256.344, 267.893,计算出1995年的产量，用三次样条插值的方法，画出每隔一年的插值 曲线图形，同时将原始的数据画在同一图上。解：程序如下year=1900:10:2010;product二75.9

4、95, 91.972, 105.711, 123.203, 131.699, 150.697, 179.323, 203.212, 226.505, 249.633, 256.344, 267.893p1995=i nterp1(year,product,1995) x=1900:2010;y=in terp1(year,product,x,'cubic');plot(year,product,'o',x,y);计算结果为：p1995=252.98852、二维插值已知离散点上的数据集(冷,弓),&22忆2),(召,召)，即已知在点集(x1,y1),(x2

5、,y2),(Xn,%)上的函数值乙厶，Zn，构造一个解析函数(其图形为一 曲面）通过这些点，并能够求出这些已知点以外的点的函数值，这一过程称为 二维插值。MATLAB 函数：Zi=interp2（X,Y ,Z,Xi,Yi,method）该命令用指定的算法找出一个二元函数 z f（x，y），然后以f（x，y）给出（x,y）处 的值。返回数据矩阵乙，Xi，Yi是向量，且必须单调，乙和meshgrid（Xi,Yi）是 同类型的。method可以下列方法之一： nearest最近邻点插值，直接完成计算； spline三次样条函数插值；linear线性插值（缺省方式），直接完成计算； cubic 

6、9;次函数插值；例2:已知1950年到1990年间每隔10年，服务年限从10年到30年每隔 10年的劳动报酬表如下：表：某企业工作人员的月平均工资（元）服务年限年份1020301950150.697169.592187.6521960179.323195.072250.2871970203.212239.092322.7671980226.505273.706426.7301990249.633370.281598.243试计算1975年时，15年工龄的工作人员平均工资解：程序如下：years=1950:10:1990;service=10:10:30;wage二150.697 169.592

7、 187.652179.323 195.072 250.287203.212 239.092 322.767226.505 273.706 426.730249.633 370.281 598.243;mesh(service,years,wage) % 绘原始数据图w=interp2(service,years,wage,15,1975); %求点(15,1975)处的值计算结果为：235.6288例3:设有数据x=1,2,3,4,5,6, y=1,2,3,4,在由x,y构成的网格上，数据为:12,10,11,11,13,1518,21,26,32,28,25求通过这些点的插值曲面解：程序为

8、：x=1:6;y=1:4;t=12,10,11,11,13,1516,22,28,35,27,2018,21,26,32,28,25;20,25,30,33,32,20subplot(1,2,1)mesh(x,y,t)x1=1:0.1:6;y1=1:0.1:4;x2,y2=meshgrid(x1,y1);t1=interp2(x,y,t,x2,y2,'cubic');subplot(1,2,2)mesh(x1,y1,t1);结果如右图。作业： 已知某处山区地形选点测量坐标数据为：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

9、.5 3 3.5 4 4.5 5 5.5 6海拔高度数据为：z=89 90 87 85 92 91 96 93 90 87 8292 96 98 99 95 91 89 86 84 82 8496 98 95 92 90 88 85 84 83 81 8580 81 82 89 95 96 93 92 89 86 8682 85 87 98 99 96 97 88 85 82 8382 85 89 94 95 93 92 91 86 84 8888 92 93 94 95 89 87 86 83 81 9292 96 97 98 96 93 95 84 82 81 84 85 85 81 82

10、 80 80 81 85 90 93 9584 86 81 98 99 98 97 96 95 84 8780 81 85 82 83 84 87 90 95 86 8880 82 81 84 85 86 83 82 81 80 8287 88 89 98 99 97 96 98 94 92 871、画出原始数据图；2、画出加密后的地貌图，并在图中标出原始数据1、拟合曲线拟合已知离散点上的数据集（冷）"22），（Xn，yn），即已知在点集石必，,召上的函数值%,齐，构造一个解析函数（其图形为一曲线）使f（x）在原离散点x上尽可能接近给定的yi值，这一过程称为曲线拟合。最常用的曲线拟合

11、方法是n2Mf（xj y,最小二乘法，该方法是寻找函数f（x）使得 i 1最小。MATLA函数：p=polyfit（x,y,n）p,s= polyfit(x,y, n)说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低 的多项式系数向量p。x必须是单调的。矩阵s用于生成预测 值的误差估计。（见下一函数polyval）多项式曲线求值函数：polyval（）调用格式：y=polyval(p,x) y,DELTA二polyval(p,x,s)说明：y=polyval（p,x）为返回对应自变量x在给定系数P的多项式的值 y,DELTA=polyval(p,x,s)使用polyfit函数的选项输

12、出s得出误差估计Y DELTA。 它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则Y DELTA 将至少包含50%的预测值。例5:求如下给定数据的拟合曲线，x=0.5,1.0,1.5,2.0,2.5,3.0, y二1.75,2.45,3.81,4.80,7.00,8.60。解：MATLAB程序如下：x=0.5,1.0,1.5,2.0,2.5,3.0;y二1.75,2.45,3.81,4.80,7.00,8.60;p=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'

13、;-b')计算结果为：p =0.5614 0.8287 1.1560此结果表示拟合函数为：f(x) 0.5614x2+0.8287x+1.1560,用此函数拟合数据的效果如图所示例2:由离散数据x _0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52拟合出多项式。程序：x=0:.1:1;y=.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2n=3;p=polyfit(x,y, n)xi=li nspace(0,1,100);z=polyval(p,xi); %多项式求值plot(x,y, ' o' ,xi,z,

14、' k: ' ,x,y,' b')legend(原始数据',阶曲线')结果：p =16.7832 -25.7459 10.9802 -0.0035多项式为：16.7832x3-25.7459x2+10.9802x-0.0035曲线拟合图形：也可由函数给出数据。例 3： x=1:20,y=x+3*sin(x)程序：x=1:20;y=x+3*s in (x);p=polyfit(x,y,6) xi=linspace(1,20,100);z=polyval(p,xi);%?d ?e ?o ?卩 o色 yplot(x,y, 'o' ,xi

15、,z, 'k:',x,y,'b' )结果：p =0.0000 -0.0021 0.0505 -0.5971 3.6472 -9.7295 11.3304再用 10 阶多项式拟合程序： x=1:20;y=x+3*sin(x); p=polyfit(x,y,10) xi=linspace(1,20,100);z=polyval(p,xi); plot(x,y,'o',xi,z,'k:',x,y,'b') 结果： p =Columns 1 through 70.0000 -0.0000 0.0004 -0.0114 0.1814 -1.8065 11.2360Columns 8 through 11-42.0861 88.5907 -92.8155 40.267可用不同阶的多项式来拟合数据，但也不是阶数越高拟合的越好作业：1.已知 x二0.1,0.8,1.3,1.925,3.1, y二1.2,1.627,2.0,1.3,0.5，利用其中的部 分数据，分别用线性函数插值， 3 次函数插值，求 x=2.0 处的值。2.已知二元函数 zf (x, y) 在点集 D( x,y) |x 0,1,2,3,4; y 0,1,2,3,4 上