数值计算3-插值和曲线拟合

《数值计算3-插值和曲线拟合》由会员分享，可在线阅读，更多相关《数值计算3-插值和曲线拟合（15页珍藏版）》请在装配图网上搜索。

1、数值计算3插值和曲线拟合插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在假设干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。如何根据观测点的值，构造一个比拟简单的函数y=(x)，使函数在观测点的值等于的数值或导数值。用简单函数y=(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数(x)，方法是很多的。(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；(x)可以是任意光滑任意阶导数连续的函数或是分段函数。函数类的不同，自然地有

2、不同的逼近效果。在许多应用中，通常要用一个解析函数一、二元函数来描述观测数据。根据测量数据的类型：1测量值是准确的，没有误差。2测量值与真实值有误差。这时对应地有两种处理观测数据方法：1插值或曲线拟合。2回归分析假定数据测量是精确时，一般用插值法，否那么用曲线拟合。MATLAB中提供了众多的数据处理命令。有插值命令，有拟合命令，有查表命令。一维插值 插值定义为对数据点之间函数的估值方法，这些数据点是由某些集合给定。当人们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况。 interp1(x,y,xi,method)

3、 x和y为既有数据的向量，其长度必须相同。 xi为要插值的数据点向量。 method插值方法,nearest/linear/cubic/spline之一，分别为最近点插值/线性插值/分段三次Hermite插值/三次样条插值。 例x=1.0 2.0 3.0 4.0 5.0; %输入变量数据x y=11.2 16.5 20.4 26.3 30.5; %输入变量数据y x1=2.55; %输入待插值点x y11=interp1(x,y,x1,nearest) %最近点插值方法的插值结果y12=interp1(x,y,x1,linear) %线性插值方法的插值结果y13=interp1(x,y,x1,

4、cubic) %三次Hermite插值方法的插值结果y14=interp1(x,y,x1,spline) %样条插值方法的插值结果 y11 = 20.4000y12 = 18.6450y13 = 18.6028y14 = 18.4874 plot(x,y) 或许最简单插值的例子是MATLAB的作图。按缺省，MATLAB用直线连接所用的数据点以作图。这个线性插值猜想中间值落在数据点之间的直线上。当然，当数据点个数的增加和它们之间距离的减小时，线性插值就更精确。例如： x1=linspace(0, 2*pi, 60);x2=linspace(0, 2*pi, 6);plot(x1, sin(x1)

5、, x2, sin(x2), - )xlabel( x ), ylabel( sin(x) ), title( Linear Interpolation ) 上图是sin函数的两个图，一个在数据点之间用60个点，它比另一个只用6个点更光滑和更精确。 根据所作的假设，有多种插值。而且，可以在一维以上空间中进行插值。即如果有反映两个变量函数的插值，z=f(x, y)，那么就可在x之间和在y之间，找出z的中间值进行插值。MATLAB在一维函数interp1和在二维函数interp2中，提供了许多的插值选择。其中的每个函数将在下面阐述。 为了说明一维插值，考虑以下问题，12小时内，一小时测量一次室外温

6、度。数据存储在两个MATLAB变量中。hours=1:12; % index 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 ) 在线性插值下室外温度曲线正如上图看到的，MATLAB画出了数据点线性插值的直线。为了计算在任意给定时间的温度

7、，人们可试着对可视的图作解释。另外一种方法，可用函数interp1。t=interp1(hours, temps, 9.3) % estimate temperature at hour=9.3 t = 22.9000 t=interp1(hours, temps, 4.7) % estimate temperature at hour=4.7 t = 22 t=interp1(hours, temps, 3.2 6.5 7.1 11.7) % find temp at many points! t = 10.2000 30.0000 30.9000 24.9000 interp1的缺省用法是

8、由interp1(x, y, xo)来描述，这里x是独立变量(横坐标)，y是应变量(纵坐标)，xo是进行插值的一个数值数组。另外，该缺省的使用假定为线性插值。假设不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的头两个导数与该数据点相一致。这种类型的插值被称为3次样条或简称为样条。函数interp1也能执行3次样条插值。 hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24; plot(hours, temps, hours, temps, +

9、 ) title( Temperature )xlabel(时间轴), ylabel(温度变化轴) x1=9.3 2.5;t1=interp1(hours,temps,x1,spline) t1 = 21.8577 8.2364 x2=3.2 6.5 7.1 11.7;t3=interp1(hours,temps,x2,spline) t3 = 9.6734 30.0427 31.1755 25.3820 t3 = 9.6734 30.0427 31.1755 25.3820 注意，样条插值得到的结果，与上面所示的线性插值的结果不同。因为插值是一个估计或猜想的过程，其意义在于，应用不同的估计规

10、那么导致不同的结果。 一个最常用的样条插值是对数据平滑。也就是，给定一组数据，使用样条插值在更细的间隔求值。例如， hours=1:12; temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12; % estimate temperature every 1/10 hourt=interp1(hours,temps,h,spline) ;plot(hours, temps, - , hours, temps, + , h, t) % plot comparative resultstitle( Springfield Temperature )xla

11、bel( Hour ), ylabel( Degrees Celsius ) 在不同插值下室外温度曲线 在上图中，虚线是线性插值，实线是平滑的样条插值，标有 + 的是原始数据。如要求在时间轴上有更细的分辨率，并使用样条插值，我们有一个更平滑、但不一定更精确地对温度的估计。尤其应注意，在数据点，样条解的斜率不突然改变。作为这个平滑插值的回报，3次样条插值要求更大量的计算，因为必须找到3次多项式以描述给定数据之间的特征。 在讨论二维插值之前，了解interp1所强制的二个强约束是很重要的。首先，人们不能要求有独立变量范围以外的结果，例如，interp1(hours, temps, 13.5)导致一

12、个错误，因为hours在1到12之间变化。其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减小的。在我们的例子里，hours是单调的。然而，如果我们已经定义独立变量为一天的实际时间， time_of_day=7:12 1:6 % start at 7AM,end at 6PM time_of_day = 7 8 9 10 11 12 1 2 3 4 5 6那么独立变量将不是单调的，因为time_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours，将会返回一个错误。同样的理由，人们不能对temps插值来找出产生某温度的时间

13、(小时)，因为temps不是单调的。二维插值interp2(x,y,z,xi,yi,method) x,y,z为既有数据的向量，其长度必须相同。 xi,yi为要插值的数据点。 method为插值方法，可为nearest、linear、cubic之一， 分别为最近点插值、双线性插值、双三次插值%例：x=0.0 0.5 1.0 1.5 2.0 2.5 3.0; %输入变量数据x y=0.0 0.5 1.0 1.5 2.0 2.5 3.0; %输入变量数据y %输入变量数据z z=100 99 100 99 100 99 100 100 99 99 99 100 99 100 99 99 98 98

14、 100 99 100 100 98 97 97 99 100 103 101 100 98 98 100 102 104 102 103 101 100 102 106 106 99 102 100 100 103 108 103; x1=2.3;y1=2.8; %输入待插值点x1,y1 z11=interp2(x,y,z,x1,y1,nearest) %最近点插值方法的插值结果z12=interp2(x,y,z,x1,y1,bilinear) %双线性插值插值方法的插值结果z13=interp2(x,y,z,x1,y1,bicubic) %双三次插值插值方法的插值结果 z11 = 108z

15、12 = 105.3600z13 = 105.9744 二维插值是基于与一维插值同样的根本思想。然而，正如名字所隐含的，二维插值是对两变量的函数z=f(x, y) 进行插值。为了说明这个附加的维数，考虑一个问题。设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板外表均匀分布的格栅。 采集了以下的数据：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 80 82 84; 79 63

16、 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-维相联系(见图)。为了估计在中间点的温度，我们必须对它们进行辨识。wi=1:0.2:5; % estimate across width of plated=2; % at a depth of 2zlinear=interp2(width, depth, temps,

17、wi, d) ; %linear interpolationzcubic=interp2(width, depth, temps, wi,d, cubic) ; % cubic interpolationplot(wi, zlinear, - , wi, zcubic,r) % plot resultsxlabel( Width of Plate ), ylabel( Degrees Celsius )title( Temperature at Depth = num2str(d) ) 在深度d=2处的平板温度 另一种方法，我们可以在两个方向插值。先在三维坐标画出原始数据，看一下该数据的粗糙程

18、度(见图11.7)。mesh(width, depth, temps) % use mesh plotxlabel( Width of Plate ), ylabel( Depth of Plate )zlabel( Degrees Celsius ), axis(ij), grid 平板温度 然后在两个方向上插值，以平滑数据。width=1:5; depth=1:3; temps=82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86;di=1:0.2:3; % choose higher resolution for depthwi=1:0.2:5;

19、% choose higher resolution for widthzcubic=interp2(width,depth,temps,wi,di,cubic) ; % cubicmesh(wi, di, zcubic)xlabel( Width of Plate ), ylabel( Depth of Plate )zlabel( Degrees Celsius ),axis(ij),grid 上面的例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持跟踪。interp2的根本形式是interp2(x, y, z, xi, yi, method)。这里x和y是两个独立变量，z是一个

20、应变量矩阵。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-轴插值的一个数值数组。图11.8 二维插值后的平板温度 可选的参数method可以是 linear，cubic或nearest。在这种情况下，cubic不意味着3次样条，而是使用3次多项式的另一种算法。linear方法是线性插值，仅用作连接图上数据点。nearest方法只选择最接近各估计点的粗略数据点。在所有的情况下，假

21、定独立变量x和y是线性间隔和单调的。关于这些方法的更多的信息，可请求在线帮助。二曲线拟合使用多项式曲线拟合 p=polyfit(x,y,n) p为最小二乘意义上拟合多项式的相关系数。 x,y为数据点向量。n为多项式阶次。 y=polyval(p,x) y为自变量为x、系数为p的多项式值. x为矩阵或向量时对其每一点进行操作. p为一向量，各元素为多项式系数(降幂排列). 曲线拟合涉及答复两个根本问题：最正确拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义最正确拟合，并存在无穷数目的曲线。所以，从这里开始，我们走向何方？正如它证实的那样，当最正确拟合被解释为在数据点的最小误差平方和，且

22、所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。数学上，称为多项式的最小二乘曲线拟合。如果这种描述使你混淆，再研究图。虚线和标志的数据点之间的垂直距离是在该点的误差。对各数据点距离求平方，并把平方距离全加起来，就是误差平方和。这条虚线是使误差平方和尽可能小的曲线，即是最正确拟合。最小二乘这个术语仅仅是使误差平方和最小的省略说法。在MATLAB中，函数polyfit求解最小二乘曲线拟合问题。为了阐述这个函数的用法，让我们以上图中的数据开始。x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1; y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56

23、 9.48 9.30 11.2; 为了用polyfit，我们必须给函数赋予上面的数据和我们希望最正确拟合数据的多项式的阶次或度。如果我们选择n=1作为阶次，得到最简单的线性近似。通常称为线性回归。相反，如果我们选择n=2作为阶次，得到一个2阶多项式。现在，我们选择一个2阶多项式。n=2; % polynomial orderp=polyfit(x, y, n) p = -9.8108 20.1293 -0.0317 polyfit 的输出是一个多项式系数的行向量。其解是y = 9.8108x2 20.1293x0.0317。为了将曲线拟合解与数据点比拟，让我们把二者都绘成图。xi=linspa

24、ce(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 ), ylabel( y=f(x) ), title( Second Order Curve Fitting ) 将图作标志。这些步骤的结果表示于前面的图中。例设0.00.30.81.11.62.3，y=0

25、.500.821.141.251.351.40，试求二次多项式拟合系数，并据此计算x1=0.9 1.2时对应的y1。x=0.0 0.3 0.8 1.1 1.6 2.3; %输入变量数据xy=0.50 0.82 1.14 1.25 1.35 1.40; %输入变量数据yp=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数px1=0.9 1.2; %输入点x1y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318y1 = 1.1656 1.2909 使用指定函数进行曲线拟合例2使用试根据土的压缩实验数据数据，用双曲线模型确定模

26、型参数，压力0时的孔隙比，可用天然孔隙比e0代替，p=0,0.05,0.1,0.2,0.4,0.6,0.8,1.2MPa对应的孔隙比分别为e=1.335,1.253,1.180,1.058,0.887,0.803,0.752,0.685。Data =0.0000 1.335;0.0500 1.253;0.1000 1.180;0.2000 1.058;0.4000 0.887;0.6000 0.803;0.8000 0.752;1.2000 0.685;p=Data(:,1);e=Data(:,2);plot(p,e,g-);hold on;x0=1 0;x=lsqnonlin(f0405,x

27、0)e1=e(1);e=e1-p./(x(1)+x(2)*p);plot(p,e,bo-); Optimization terminated: relative function value changing by less than OPTIONS.TolFun.x = 0.4833 1.1104function f = f0405(x,Data)Data =0.0000 1.335;0.0500 1.253;0.1000 1.180;0.2000 1.058;0.4000 0.887;0.6000 0.803;0.8000 0.752;1.2000 0.685;p=Data(:,1);e=Data(:,2);e1=e(1);z=e1-p./(x(1)+x(2)*p);f = z-e; 作业： 0.00.30.81.11.62.3，0.400.561.141.311.562.10，试用MatLab对数据点对进行二次多项式拟合，并使用不同插值方法计算处对应的。 实测数据如下：0.00.51.01.5；0.00.51.01.5；50495049；49494949；49484847；48484746试用MatLab求11.2、11.4处的插值结果。求解方程组