《第X章 MATLAB在拟合与插值中的应用》由会员分享，可在线阅读，更多相关《第X章 MATLAB在拟合与插值中的应用（14页珍藏版）》请在金锄头文库上搜索。
1、MATLAB在拟合与插值中的应用在大量的应用领域中，人们经常面临用一个解析函数描述数撼通常是测量值）的任务。（比 如在土木工程中对实验梁的应力应变（。-）曲线的数据进行拟合，从而得出钢筋混凝土的 弹性模量的计算式。）在这里讨论的方法是曲线拟合与插值。其中包括曲线拟合，一维插值， 二维插值以及如何解决插值中求值时的单调性问题。曲线拟合曲线拟合涉及回答两个基本问题:最佳拟合意味着什么？应该用什么样的曲线？可用许 多不同的方法定义最佳拟合，并存在无穷数目的曲线。我们将最佳拟合解释为在数据点的最 小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。先看看图1在MATLAB中，函数pol
2、yfit求解最小二乘曲线拟合问题。简单阐述这个函数的用法， 让我们以上面图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;% polynomial order p=polyfit（x, y, n）P =-9.8
3、10820.1293-0.0317polyfit的输出是一个多项式系数的行向量。其解是y = 9.8108x2 + 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和z。 xlabel( x )
4、, ylabel( y=f(x) ), title( Second Order Curve Fitting )将图作标志。这些步骤的结果表示于前面的图1中。多项式阶次的选择是有点任意的。两点决定一直线或一阶多项式。三点决定一个平方或 2阶多项式。按此进行，n+1数据点唯一地确定n阶多项式。于是，在上面的情况下，有11 个数据点，我们可选一个高达10阶的多项式。然而，高阶多项式给出很差的数值特性，我 们不应选择比所需的阶次高的多项式。此外，随着多项式阶次的提高，近似变得不够光滑， 因为较高阶次多项式在变零前，可多次求导。不妨选一个10阶多项式 pp=polyfit(x, y, 10); form
5、at 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
6、, : , xi, zz) xlabel( x ), ylabel( y=f(x) ), title( 2nd and 10th Order curve Fitting )在下面的图11.2中，原始数据标以o, 2阶曲线拟合是虚线，10阶拟合是实线。注意， 在10阶拟合中，在左边和右边的极值处，数据点之间出现大的纹波。当企图进行高阶曲线 拟合时，这种纹波现象经常发生。根据图2,显然，越多就越好在这里不适用。维插值正如曲线拟合所描述的那样，插值定义为对数据点之间函数的估值方法，这些数据点是 由某些集合给定。当我们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。 例如，当数据点是某些实
7、验测量的结果或是过长的计算过程时，就有这种情况。举例一维插值，考虑下列问题，(由于手头没有我们专业相关的实验数据，故采用老师 经常提到的测温的例子，只是数据名称不一样)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(hours, temps, hours, temps, + )% view temperatures title( Temper
8、ature ) xlabel( Hour ), ylabel( Degrees Celsius )Temperature0 5 0 53 2 2 1seerosD10 -+-+51111024681012Hour图3在线性插值下室外温度曲线正如图3看到的，MATLAB画出了数据点线性插值的直线。为了计算在任意给定时间 的温度，人们可试着对可视的图作解释。另外一种方法，可用函数interpl o t=interp1(hours, temps, 9.3)% estimate temperature at hour=9.3t =22.9000 t=interp1(hours, temps, 4.7)
9、% estimate temperature at hour=4.7t =22 t=interp1(hours, temps, 3.2 6.5 7.1 11.7)% find temp at many points!t =10.200030.000030.900024.9000若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法 是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的 头两个导数与该数据点相一致。这种类型的插值被称为3次样条或简称为样条。函数interpl 也能执行3次样条插值。（这在工程中经常用到） t=interp1(hou
10、rs, temps, 9.3, spline )% estimate temperature at hour=9.3t =21.8577 t=interp1(hours, temps, 4.7, spline )% estimate temperature at hour=4.7t =22.3143 t=interp1(hours, temps, 3.2 6.5 7.1 11.7, spline )t =9.673430.042731.175525.3820样条插值得到的结果，与上面所示的线性插值的结果不同。因为插值是一个估计或猜测 的过程，其意义在于，应用不同的估计规则导致不同的结果。一个最
11、常用的样条插值是对数据平滑。也就是，给定一组数据，使用样条插值在更细的 间隔求值。例如， 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中，虚线是线性插值，实线是平滑的样
12、条插值，标有+的是原始数据。如要求 在时间轴上有更细的分辨率，并使用样条插值，我们有一个更平滑、但不一定更精确地对温 度的估计。尤其应注意，在数据点，样条解的斜率不突然改变。作为这个平滑插值的回报， 3次样条插值要求更大量的计算，因为必须找到3次多项式以描述给定数据之间的特征。Springfield Temperature 353025C 20D 15105024681012Hour图4在不同插值下室外温度曲线二维插值二维插值是基于与一维插值同样的基本思想。然而，正如名字所隐含的，二维插值是对 两变量的函数z=f（x, y）进行插值（比如钢筋混凝土实验中的正应力和剪应力都对挠度产生 影响）。这
13、里依然考虑温度问题。（数据由课件中改动而成）设人们对平板上的温度分布估计 感兴趣，给定的温度值取自平板表面均匀分布的格栅。采集了下列的数据： width=1:5;% index for width of plate (i.e.,the x-dimension) depth=1:3;% index for depth of plate (i,e,the y-dimension) temps=82 81 8082 84; 79 63 61 65 81; 84 84 82 85 86 %temperature datatemps =828180828479636165818484828586如同在标
14、引点上测量一样，矩阵temps表示整个平板的温度分布。temps的列与下标 depth或y-维相联系，行与下标width或x-维相联系（见图5）。为了估计在中间点的温度， 我们必须对它们进行辨识。 wi=1:0.2:5;% estimate across 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 interpol
15、ation plot(wi, zlinear, - , wi, zcubic) % plot results xlabel( Width of Plate ), ylabel( Degrees Celsius ) title( Temperature at Depth = num2str(d)另一种方法，我们可以在两个方向插值。先在三维坐标画出原始数据，看一下该数据的 粗糙程度(见图6)。 mesh(width, depth, temps)% use mesh plot xlabel( Width of Plate ), ylabel( Depth of Plate ) zlabel( Degrees Celsius ), axis
