第X章 MATLAB在拟合与插值中的应用Word格式文档下载.docx

1、 plot（x, y, o , x, y, xi, z, : ） 画出了原始数据x和y，用o标出该数据点，在数据点之间，再用直线重画原始数据，并用点线，画出多项式数据xi和z。 xlabel（ x ）, ylabel（ y=f（x） ）, title（ Second Order Curve Fitting ） 将图作标志。这些步骤的结果表示于前面的图1中。 多项式阶次的选择是有点任意的。两点决定一直线或一阶多项式。三点决定一个平方或2阶多项式。按此进行，n+1数据点唯一地确定n阶多项式。于是，在上面的情况下，有11个数据点，我们可选一个高达10阶的多项式。然而，高阶多项式给出很差的数值特性，我

2、们不应选择比所需的阶次高的多项式。此外，随着多项式阶次的提高，近似变得不够光滑，因为较高阶次多项式在变零前，可多次求导。不妨选一个10阶多项式 pp=polyfit（x, y, 10） ; format short e pp.则 ans = -4.6436e+005 2.2965e+006 -4.8773e+006 5.8233e+006 -4.2948e+006 2.0211e+006 -6.0322e+005 1.0896e+005 -1.0626e+004 4.3599e+002 -4.4700e-001要注意在现在情况下，多项式系数的规模与前面的2阶拟合的比较。还要注意在最小（-4.4

3、700e-001）和最大（5.8233e+006）系数之间有7个数量级的幅度差。将这个解作图，并把此图与原始数据及2阶曲线拟合相比较。 zz=polyval（pp, xi）; , xi, z, , xi, zz） 2nd and 10th Order curve Fitting 在下面的图11.2中，原始数据标以，2阶曲线拟合是虚线，10阶拟合是实线。注意，在10阶拟合中，在左边和右边的极值处，数据点之间出现大的纹波。当企图进行高阶曲线拟合时，这种纹波现象经常发生。根据图2，显然， 越多就越好 在这里不适用。图2 2阶和10阶曲线拟合一维插值 正如曲线拟合所描述的那样，插值定义为对数据点之间函

4、数的估值方法，这些数据点是由某些集合给定。当我们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况。 举例一维插值，考虑下列问题，（由于手头没有我们专业相关的实验数据，故采用老师经常提到的测温的例子，只是数据名称不一样）12小时内，一小时测量一次室外温度。数据存储在两个MATLAB变量中。 hours=1:12; % index for hour data was recorded temps=5 8 9 15 25 29 31 30 22 25 27 24; % recorded temperatures plot

5、（hours, temps, hours, temps, + ） % view temperatures title（ Temperature Hour Degrees Celsius 图3 在线性插值下室外温度曲线 正如图3看到的，MATLAB画出了数据点线性插值的直线。为了计算在任意给定时间的温度，人们可试着对可视的图作解释。另外一种方法，可用函数interp1。 t=interp1（hours, temps, 9.3） % estimate temperature at hour=9.3 t = 22.9000 t=interp1（hours, temps, 4.7） % estimat

6、e temperature at hour=4.7 22 t=interp1（hours, temps, 3.2 6.5 7.1 11.7） % find temp at many points! 10.2000 30.0000 30.9000 24.9000 若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的头两个导数与该数据点相一致。这种类型的插值被称为3次样条或简称为样条。函数interp1也能执行3次样条插值。（这在工程中经常用到） t=interp1（hours, temps,

7、9.3, spline ） % estimate temperature at hour=9.3 21.8577 t=interp1（hours, temps, 4.7, ） % estimate temperature at hour=4.7 22.3143 t=interp1（hours, temps, 3.2 6.5 7.1 11.7, ） 9.6734 30.0427 31.1755 25.3820样条插值得到的结果，与上面所示的线性插值的结果不同。因为插值是一个估计或猜测的过程，其意义在于，应用不同的估计规则导致不同的结果。 一个最常用的样条插值是对数据平滑。也就是，给定一组数据，使

8、用样条插值在更细的间隔求值。例如， h=1:0.1: % estimate temperature every 1/10 hour t=interp1（hours, temps, h, ） ; plot（hours, temps, - , hours, temps, , h, t） % plot comparative results Springfield Temperature 在图4中，虚线是线性插值，实线是平滑的样条插值，标有的是原始数据。如要求在时间轴上有更细的分辨率，并使用样条插值，我们有一个更平滑、但不一定更精确地对温度的估计。尤其应注意，在数据点，样条解的斜率不突然改变。作为这

9、个平滑插值的回报，3次样条插值要求更大量的计算，因为必须找到3次多项式以描述给定数据之间的特征。图4 在不同插值下室外温度曲线二维插值 二维插值是基于与一维插值同样的基本思想。然而，正如名字所隐含的，二维插值是对两变量的函数z=f（x, y） 进行插值（比如钢筋混凝土实验中的正应力和剪应力都对挠度产生影响）。这里依然考虑温度问题。（数据由课件中改动而成）设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅。 采集了下列的数据： width=1:5; % index for width of plate （i.e.,the x-dimension） depth=1:3; %

10、 index for depth of plate （i,e,the y-dimension） temps=82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86 % temperature data temps = 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 如同在标引点上测量一样，矩阵temps表示整个平板的温度分布。temps的列与下标depth或y-维相联系，行与下标width或x-维相联系（见图5）。为了估计在中间点的温度，我们必须对它们进行辨识。 wi=1:0.2: % estimate acros

11、s width of plate d=2; % at a depth of 2 zlinear=interp2（width, depth, temps, wi, d） ; % linear interpolation zcubic=interp2（width, depth, temps, wi,d, cubic % cubic interpolation plot（wi, zlinear, , wi, zcubic） % plot results Width of Plate title（ Temperature at Depth = num2str（d） ） 另一种方法，我们可以在两个方向插

12、值。先在三维坐标画出原始数据，看一下该数据的粗糙程度（见图6）。 mesh（width, depth, temps） % use mesh plot Depth of Plate zlabel（）, axis（ ij ）, grid图5 在深度d=2处的平板温度图6 平板温度 然后在两个方向上插值，以平滑数据。 di=1: % choose higher resolution for depth % choose higher resolution for width zcubic=interp2（width, depth, temps, wi, di, % cubic mesh（wi, di

13、, zcubic） 该例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持跟踪。interp2的基本形式是interp2（x, y, z, 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-轴插值的一个数值数组。图7 二维插值后的平板温度 虽然对于许多应用，函数interp1和interp2是很

14、有用的，但它们限制为对单调向量进行插值。在某些情况，这个限制太严格。例如，考虑下面的插值： x=linspace（0, 5）; y=1-exp（-x）.*sin（2*pi*x）; plot（x, y）图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值：图8 给y值在函数曲线上求x的值 plot（x, y, 0, 5, 1.1 1.1 ） 从图8上，我们看到有四个交点

15、。使用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.5281 1.1000 0.9580 1.1000 1.5825 1.1000 1.8847

16、1.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,VAL） linearly interpolates the table % TAB searching

17、 for the scalar value VAL 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 VAL is not found in the table,Y=. rt, ct=size（

18、tab）; if length（val） 1, error（ VAL must be a scalar. ）, end if colct|col Chosen column outside table width. if rt val; % True where VAL below=tab（: , col） % True where length（ib）; % True where equals values fit ry=length（tmp）; % # of rows in result y y=zeros（ry, ct）; % poke data into a zero matrix a

19、lpha=（val-tab（ib,col）./（tab（ia,col）-tab（ib,col）; alpha=alpha（: , ones（1, ct）; % duplicate for all columns y（ieq, : ）=alpha.*tab（ia, : ）+（1-alpha）.*tab（ib, : % interpolated values y（ieq, : ）=tab（ie, : % equal values y（ : , col）=val*ones（ry, 1）; % remove roundoff error 正如所见的，mminterp利用了find和sort函数、逻辑数

20、组和数组操作技术。没有For循环和While循环。不论用其中哪一种技术来实现将使运行变慢，尤其对大的表。mminterp与含有大于或等于2的任意数列的表一起工作，如同函数interp1一样。而且，在这种情况下，插值变量可以是任意的列。 z=sin（pi*x）; % add more data to table y; z. t=mminterp（table, 2, 1.1） % same interpolation as earlier t = 0.5281 1.1000 0.9930 0.9580 1.1000 0.1314 1.5825 1.1000 -0.9639 1.8847 1.100

21、0 -0.3533 t=mminterp（table, 3, -.5） % second third column now 1.1669 0.7316 -0.5000 1.8329 1.1377 -0.5000 3.1671 0.9639 -0.5000 3.8331 1.0187 -0.5000这些最后的结果估计了x和y在z= -0.5处的值。小结 曲线的插值和拟合是一个很复杂的工作，但在MATLAB中能由几句轻松的命令来实现，为工程技术人员和科研工作者带来极大的方便，让人不禁感叹它的强大，实为工科学生必备之工具！ 下面的表总结了在MATLAB中所具有的曲线拟合和插值函数。可供同学们参考。曲 线 拟 合 和 插 值 函 数polyfit（x, y, n）对描述n阶多项式y=f（x）的数据进行最小二乘曲线拟合interp1（x, y, xo）1维线性插值interp1（x, y, xo, 1维3次样条插值1维3次插值interp2（x, y, Z, xi, yi）2维线性插值interp2（x, y, Z, xi, yi, 2维3次插值 nearest 2维最近邻插值参考文献精通Matlab综合辅导与指南