MATLA在拟合与插值中的应用在大量的应用领域中，人们经常面临用一个解析函数描述数据 （通常是测量值）的任务比如在土木工程中对实验梁的应力应变 （6 -- £）曲线的数据进行拟合，从而得出钢筋混凝土的弹性模量的计算式在这里讨论的方法是曲线拟合与插值 其中包括曲线拟合， 一维插值，二维插值以及如何解决插值中求值时的单调性问题曲线拟合曲线拟合涉及回答两个基本问题： 最佳拟合 意味着什么？应该用什么样的曲线？可用许 多不同的方法定义 最佳拟合，并存在无穷数目的曲线 我们将最佳拟合解释为在数据点的最 小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的先看看图 1在MATLAB中，函数polyfit求解最小二乘曲线拟合问题简单阐述这个函数的用法, 让我们以上面图11.1中的数据开始 x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 .11];? y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];为了用polyfit，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶 次或度如果我们选择 n=1作为阶次，得到最简单的线性近似。
通常称为线性回归如果我们选择n=2作为阶次，得到一个 2阶多项式 n=2; % polyno mial order? p=polyfit（x, y, n）p =-9.8108 20.1293 -0.0317polyfit 的输出是一个多项式系数的行向量其解是 y = － 9.8108x 2 ＋ 20.1293x － 0.0317为了将曲线拟合解与数据点比较，把二者都绘成图 xi=linspace(0, 1, 100); % x-axis data for plotting? z=polyval(p, xi);为了计算在 xi 数据点的多项式值，调用 MATLAB 的函数 polyval plot(x, y, ' o ' , x, y, xi, z, ' : ' )画出了原始数据x和y，用'o'标出该数据点，在数据点之间，再用直线重画原始数据， 并用点 ' : ' 线，画出多项式数据 xi 和 zxlabel(' x '), ylabel(' y=f(x) '), title(' Second Order Curve Fitting ')将图作标志这些步骤的结果表示于前面的图 1 中 多项式阶次的选择是有点任意的。
两点决定一直线或一阶多项式 三点决定一个平方或2 阶多项式按此进行， n+1 数据点唯一地确定 n 阶多项式于是，在上面的情况下，有 11 个数据点，我们可选一个高达 10 阶的多项式然而，高阶多项式给出很差的数值特性，我 们不应选择比所需的阶次高的多项式此外，随着多项式阶次的提高，近似变得不够光滑， 因为较高阶次多项式在变零前，可多次求导不妨选一个 10 阶多项式? pp=polyfit(x, y, 10) ;? format short e? pp.' 则 ans =-4.6436e+0052.2965e+006-4.8773e+0065.8233e+006-4.2948e+0062.0211e+006-6.0322e+0051.0896e+005-1.0626e+0044.3599e+002-4.4700e-001要注意在现在情况下，多项式系数的规模与前面的 2 阶拟合的比较还要注意在最小(-4.4700e-001)和最大(5.8233e+006)系数之间有7个数量级的幅度差将这个解作图，并把此 图与原始数据及 2 阶曲线拟合相比较 zz=polyval(pp, xi); plot(x, y, ' o ' , xi, z, ' : ' , xi, zz)xlabel(' x '), ylabel(' y=f(x) '), title(' 2nd and 10th Order curve Fitti ng ')在下面的图11.2中，原始数据标以'o' , 2阶曲线拟合是虚线，10阶拟合是实线。
注意, 在10阶拟合中，在左边和右边的极值处，数据点之间出现大的纹波当企图进行高阶曲线 拟合时，这种纹波现象经常发生根据图 2，显然，‘越多就越好’在这里不适用图2 2阶和10阶曲线拟合维插值正如曲线拟合所描述的那样，插值定义为对数据点之间函数的估值方法， 这些数据点是由某些集合给定当我们不能很快地求出所需中间点的函数值时， 插值是一个有价值的工具例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况举例一维插值，考虑下列问题， （由于手头没有我们专业相关的实验数据，故采用老师经常提到的测温的例子，只是数据名称不一样） 12小时内，一小时测量一次室外温度数据存储在两个 MATLAB 变量中hours=1:12; % in dex for hour data was recordedtemps=[5 8 9 15 25 29 31 30 22 25 27 24]; % recorded temperaturesplot(hours, temps, hours, temps,' + ') % view temperaturestitle(' Temperature ')xlabel(' Hour '), ylabel(' Degrees Celsius ')Temperature35 1 1 1 sppec seproeD图3性插值下室外温度曲线% estimate temperature at hour=9.3% estimate temperature at hour=4.7t=in terp1(hours, temps, [3.26.5 7.1 11.7])% find temp at many points!正如图3看到的，MATLAB画出了数据点线性插值的直线。
为了计算在任意给定时间 的温度，人们可试着对可视的图作解释另外一种方法，可用函数 in terp1 t=interp1(hours, temps, 9.3) t =22.9000? t=in terp1(hours, temps, 4.7) t =2210.200030.000030.900024.9000若不采用直线连接数据点， 我们可采用某些更光滑的曲线来拟合数据点 最常用的方法 是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模， 每个3次多项式的 头两个导数与该数据点相一致 这种类型的插值被称为 3次样条或简称为样条函数interp1 也能执行3次样条插值这在工程中经常用到）? t=interp1(hours, temps, 9.3, ' spline ') t =21.8577% estimate temperature at hour=9.3? t=interp1(hours, temps, 4.7, ' spline ') t =22.3143% estimate temperature at hour=4.7? t=interp1(hours, temps, [3.2 6.5 7.1 11.7], ' spline ')t =9.673430.042731.175525.3820样条插值得到的结果， 与上面所示的线性插值的结果不同。
因为插值是一个估计或猜测 的过程，其意义在于，应用不同的估计规则导致不同的结果一个最常用的样条插值是对数据平滑 也就是， 给定一组数据， 使用样条插值在更细的 间隔求值例如，? h=1:0.1:12; % estimate temperature every 1/10 hour? t=interp1(hours, temps, h, ' spline ') ;? plot(hours, temps, ' - ' , hours, temps, ' + ' , h, t) % plot comparative results? title(' Springfield Temperature ')? xlabel(' Hour '), ylabel(' Degrees Celsius ')在图 4中，虚线是线性插值，实线是平滑的样条插值，标有 ' + ' 的是原始数据如要求在时间轴上有更细的分辨率， 并使用样条插值， 我们有一个更平滑、 但不一定更精确地对温 度的估计尤其应注意，在数据点，样条解的斜率不突然改变作为这个平滑插值的回报， 3 次样条插值要求更大量的计算，因为必须找到 3 次多项式以描述给定数据之间的特征。
S^Rrec See—QBSpringfield Temperature105 0 52 2122 4 6 8 10Hour图4在不同插值下室外温度曲线二维插值二维插值是基于与一维插值同样的基本思想 然而，正如名字所隐含的，二维插值是对两变量的函数z=f（x, y）进行插值（比如钢筋混凝土实验中的正应力和剪应力都对挠度产生 影响）这里依然考虑温度问题数据由课件中改动而成） 设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅采集了下列的数据:width=1:5;depth=1:3;% in dex for width of plate (i.e.,the x-dime nsi on)% in dex for depth of plate (i,e,,the y-dime nsion)? temps=[82 81808284; 79 63temperature datatemps =82 8180828479 6361658184 8482858661 65 81; 84 84 82 85 86] %如同在标引点上测量一样，矩阵 temps表示整个平板的温度分布temps的列与下标depth或y-维相联系，行与下标 width或x-维相联系（见图5）。
为了估计在中间点的温度，我们必须对它们进行辨识wi=1:0.2:5;% estimate across width of plate? d=2; % at a depth of 2? zli near= in terp2(width, depth, temps, wi, d) ; % lin ear in terpolati on? zcubic= in terp2(width, depth, temps, wi,d, ' cubic ') ; % cubic in terpolati on? plot(wi, zli near, ' - ' , wi, zcubic) % plot results? xlabel(' Width of Plate '), ylabel(' Degrees Celsius ')? title( [' Temperature at Depth = ' nu m2str(d)])另一种方法，我们可以在两个方向插值 先在三维坐标画出原始数据， 看一下该数据的粗糙程度(见图6) mesh(width, depth, temps) % use mesh plot? xlabel(' Width of Plate '), ylabel(' Depth of Plate ')? zlabel(' Degrees Celsius '), axis(' ij '), gridsnrRrec sere6Temperature at Depth = 2o 58 7562 3 4 5Width of Plate图5在深度d=2处的平板温度90807042331Depth of PlateWidth of Plate平板温度601susDC seekyB然后在两个方向上插值，以平滑数据。
di=1:0.2:3; % choose higher resoluti on for depth? wi=1:0.2:5; % choose higher resoluti on for width? zcubic= interp2(width, depth, temps, wi, di, ' cubic ') ; % cubic? mesh(wi, di, zcubic)? xlabel(' Width of Plate '), ylabel(' Depth of Plate ')? zlabel(' Degrees Celsius '), axis(' ij '), grid该例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持跟踪 in terp2的基本形式是interp2(x, y,乙xi, yi, method)这里x和y是两个独立变量，z是一个应变量 矩阵x和y对z的关系是z(i, :) = f(x, y(i)) 和 z(:, j) = f(x(j), y).也就是，当x变化时，z的第i行与y的第i个元素y(i)相关，当y变化时，z的第j列 与x的第j个元素x(j)相关，。
xi是沿x-轴插值的一个数值数组； yi是沿y-轴插值的一个数值数组90 f80 r70、423Depth of PlateWidth of Plate6013 1图7二维插值后的平板温度虽然对于许多应用，函数interp1和interp2是很有用的，但它们限制为对 单调向量进 行插值在某些情况，这个限制太严格例如，考虑下面的插值:x=li nspace(0, 5);y=1-exp(-x).*si n(2*pi*x); plot(x, y)if all(above = = 0) | all(below = = 0), % handle simplest case1.61.4 - "1.2 - -1 -0.8 - -0.6 -0.4 - -0.2 1 ■ J 1 '0 1 2 3 4 5图 8 函数 1-exp(-x).*sin(2*pi*x)的曲线函数interp1可用来在任何值或 x的值上估计y值 yi=interp1(x, y, 1.8) yi =1.1556然而，interp1不能找出对应于某些 y值的x值例如，如在图8上所示，考虑寻找y=1.1 处的x值：1.61.41.210.80.60.40.22 3 4 5图8 给y值在函数曲线上求 x的值plot(x, y, [0, 5], [1.11.1])从图 8 上，我们看到有四个交点。
使用 interp1 ，我们得到：? xi=interp1(y, x, 1.1)??? Error using ==> table1 First column of the table must be monotonic.这个函数 interp1 失败 ，由于 y 不是单调的如何消除了单调性的要求 (我尝试搜索了一些资料，对此问题可由如下解答)? table=[x; y].' ; % create column oriented table from data? xi=mminterp(table, 2, 1.1)xi =0.52811.10000.95801.10001.58251.10001.88471.1000这里使用了线性插值，函数 mminterp 估计了 y=1.1 处的四个点由于函数 mminterp 的一般性质，要插值的数据是由面向列矩阵给出，在上面的例子中称作为表 (table )第二个输入参量是被搜索矩阵 table 的列，第三个参量是要找的值函数的主体由下面给出：function y=mminterp(tab, col, val)% MMINTERP 1-D Table Search by Linear Interpolation.% Y=MMINTERP(TAB,COL,V AL) linearly interpolates the table% TAB searching for the scalar value V AL in the column COL.% All crossings are found and TAB(:,COL) need not be monotonic. % Each crossing is returned as a separate row in Y and Y has as % many columns as TAB.Naturally,the column COL of Y contains % the value VAL. If V AL is not found in the table,Y=[].[rt, ct]=size(tab);if length(val) > 1, error(' VAL must be a scalar. '), endif col>ct|col < 1, error(' Chosen column outside table width. '), end if rt < 2, error(' Table too small or not oriented in columns. '), endabove=tab(: , col) > val; %below=tab(: , col) < val; equal=tab(: , col) = = val;True where > VAL % True where < VAL % True where = VALendy=tab(find(equal), : ); returnpslope=find(below(1:rt-1)&above(2:rt));% indices where slope is +nslope=find(below(2:rt)&above(1:rt-1));% indices where slope is -ib=sort([pslope; nslope+1]); ia=sort([nslope; pslope+1]); ie=find(equal);% put indices below in order % put indices above in order % indices where equal to val[tmp,ix]=sort( [ib, ie] ); ieq=ix > length(ib); ry=length(tmp);y=zeros(ry, ct);% find where equals fit in result % True where equals values fit % # of rows in result y% poke data into a zero matrixalpha=(val-tab(ib,col))./(tab(ia,col)-tab(ib,col));alpha=alpha(: , ones(1, ct)); % duplicate for all columnsy(~ieq, : )=alpha.*tab(ia, : )+(1-alpha).*tab(ib, : ); % interpolated valuesy(ieq, : )=tab(ie, : ); % equal valuesy( : , col)=val*ones(ry, 1); % remove roundoff error正如所见的， mminterp 利用了 find 和 sort 函数、逻辑数组和数组操作技术。
没有 For 循环和 While 循环 不论用其中哪一种技术来实现将使运行变慢， 尤其对大的表 mminterp 与含有大于或等于 2 的任意数列的表一起工作，如同函数 interp1 一样而且，在这种情况 下，插值变量可以是 任意的列 例如，? z=sin(pi*x);add more data to table% same interpolation as earlier? table=[x; y; z].' ;? t=mminterp(table, 2, 1.1) t =0.52811.10000.99300.95801.10000.13141.58251.1000-0.96391.88471.1000-0.3533minterp(table, 3, -.5)%1.16690.7316-0.50001.83291.1377-0.50003.16710.9639-0.5000t =second third column now3.83311.0187 -0.5000这些最后的结果估计了 x和y在z= -0.5处的值小结曲线的插值和拟合是一个很复杂的工作，但在 MATLAB中能由几句轻松的命令来实现，为工程技术人员和科研工作者带来极大的方便， 让人不禁感叹它的强大，实为工科学生必备之工具！！F面的表总结了在 MATLAB中所具有的曲线拟合和插值函数。
可供同学们参考曲线拟合和插值函数polyfit(x, y, n)对描述n阶多项式y=f(x)的数据 进行最小二乘曲线拟合in terp1(x, y, xo)1维线性插值in terp1(x, y, xo, ' spli ne ')1维3次样条插值in terp1(x, y, xo, ' cubic ')1维3次插值in terp2(x, y, Z, xi, yi)2维线性插值in terp2(x, y, Z, xi, yi, ' cubic ')2维3次插值in terp2(x, y, Z, xi, yi, ' n earest ')2维最近邻插值参考文献《精通Matlab综合与指南》。