第11章曲线拟合与插值 在大量的应用领域中,人们经常面临用一个解析函数描述数据(通常是测量值)的任务。 对这个问题有两种方法。在插值法里,数据假定是正确的,要求以某种方法描述数据点之间 所发生的情况。这种方法在下一节讨论。这里讨论的方法是曲线拟合或回归。人们设法找出 某条光滑曲线,它最佳地拟合数据,但不必要经过任何数据点。图11.1说明了这两种方法 标有'o’的是数据点;连接数据点的实线描绘了线性内插,虚线是数据的最佳拟合 11.1曲线拟合 曲线拟合涉及回答两个基本问题:最佳拟合意味着什么?应该用什么样的曲线?可用许 多不同的方法定义最佳拟合,并存在无穷数目的曲线。所以,从这里开始,我们走向何方? 正如它证实的那样,当最佳拟合被解释为在数据点的最小误差平方和,且所用的曲线限定为 多项式时,那么曲线拟合是相当简捷的。数学上,称为多项式的最小二乘曲线拟合。如果这 种描述使你混淆,再研宄图11.1。虚线和标志的数据点之间的垂直距离是在该点的误差。对 各数据点距离求平方,并把平方距离全加起来,就是误差平方和。这条虚线是使误差平方和 尽可能小的曲线,即是最佳拟合。最小二乘这个术语仅仅是使误差平方和最小的省略说法 Second Order Curve Fitting 10 86420 0.2 04 0.6 0.8 图1112阶曲线拟合 在 MATLAB中,函数 polyfit求解最小二乘曲线拟合问题。为了阐述这个函数的用法 让我们以上面图11.1中的数据开始
第 11 章 曲线拟合与插值 在大量的应用领域中，人们经常面临用一个解析函数描述数据(通常是测量值)的任务。 对这个问题有两种方法。在插值法里，数据假定是正确的，要求以某种方法描述数据点之间 所发生的情况。这种方法在下一节讨论。这里讨论的方法是曲线拟合或回归。人们设法找出 某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。图 11.1 说明了这两种方法。 标有'o'的是数据点；连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合。 11.1 曲线拟合 曲线拟合涉及回答两个基本问题：最佳拟合意味着什么？应该用什么样的曲线？可用许 多不同的方法定义最佳拟合，并存在无穷数目的曲线。所以，从这里开始，我们走向何方？ 正如它证实的那样，当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为 多项式时，那么曲线拟合是相当简捷的。数学上，称为多项式的最小二乘曲线拟合。如果这 种描述使你混淆，再研究图 11.1。虚线和标志的数据点之间的垂直距离是在该点的误差。对 各数据点距离求平方，并把平方距离全加起来，就是误差平方和。这条虚线是使误差平方和 尽可能小的曲线，即是最佳拟合。最小二乘这个术语仅仅是使误差平方和最小的省略说法。 0 0.2 0.4 0.6 0.8 1 -2 0 2 4 6 8 10 12 x y=f(x) Second Order Curve Fitting 图 11.1 2 阶曲线拟合 在 MATLAB 中，函数 polyfit 求解最小二乘曲线拟合问题。为了阐述这个函数的用法， 让我们以上面图 11.1 中的数据开始
》X 0.1234567891} y=4471.9783.286.167087347.669.569489.30112] 为了用 polypi,我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶 次或度。如果我们选择n=1作为阶次,得到最简单的线性近似。通常称为线性回归。相反 如果我们选择n=2作为阶次,得到一个2阶多项式。现在,我们选择一个2阶多项式。 n=2 polynomial order >p=polyfit(x, y, n) 9.810820.1293-0.0317 polyfit的输出是一个多项式系数的行向量。其解是y=-98108x2+20.293x 0.0317。为了将曲线拟合解与数据点比较,让我们把二者都绘成图。 >)Xilinspace(0, 1, 100); X-axis data for plotting zpolyval(, XI) 为了计算在xi数据点的多项式值,调用 MATLAB的函数 polygala 画出了原始数据x和y,用'o'标出该数据点,在数据点之间,再用直线重画原始数据, 并用点:线,画出多项式数据x和z >xlabel(x), ylabel( y=f(x)), title( Second Order Curve Fitting) 将图作标志。这些步骤的结果表示于前面的图11中 多项式阶次的选择是有点任意的。两点决定一直线或一阶多项式。三点决定一个平方或 2阶多项式。按此进行,n+1数据点唯一地确定n阶多项式。于是,在上面的情况下,有11 个数据点,我们可选一个高达10阶的多项式。然而,高阶多项式给出很差的数值特性,人 们不应选择比所需的阶次高的多项式。此外,随着多项式阶次的提高,近似变得不够光滑, 因为较高阶次多项式在变零前,可多次求导。例如,选一个10阶多项式 lyfit(x, y, 10) display polynomial coefficients as a column
» 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 9.48 9.30 11.2]; 为了用 polyfit，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶 次或度。如果我们选择 n=1 作为阶次，得到最简单的线性近似。通常称为线性回归。相反， 如果我们选择 n=2 作为阶次，得到一个 2 阶多项式。现在，我们选择一个 2 阶多项式。 » n=2; % polynomial order » p=polyfit(x, y, n) p = -9.8108 20.1293 -0.0317 polyfit 的输出是一个多项式系数的行向量。其解是 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 '), ylabel(' y=f(x) '), title(' Second Order Curve Fitting ') 将图作标志。这些步骤的结果表示于前面的图 11.1 中。 多项式阶次的选择是有点任意的。两点决定一直线或一阶多项式。三点决定一个平方或 2 阶多项式。按此进行，n+1 数据点唯一地确定 n 阶多项式。于是，在上面的情况下，有 11 个数据点，我们可选一个高达 10 阶的多项式。然而，高阶多项式给出很差的数值特性，人 们不应选择比所需的阶次高的多项式。此外，随着多项式阶次的提高，近似变得不够光滑， 因为较高阶次多项式在变零前，可多次求导。例如，选一个 10 阶多项式 » pp=polyfit(x, y, 10) ; » format short e % change display format » pp.' % display polynomial coefficients as a column ans =
-4.6436e+005 2.2965e+006 -4.8773e+006 5.8233e+006 -42948e+006 2.021le+006 6.0322e+005 1.0896e+005 -1.0626e+004 4.3599c+002 -4.4700e-00l 要注意在现在情况下,多项式系数的规模与前面的2阶拟合的比较。还要注意在最小 (-44700c-001)和最大(5.8233+006)系数之间有7个数量级的幅度差。将这个解作图,并把此 图与原始数据及2阶曲线拟合相比较,结果如何呢? -polyval(pp, xi); evaluate 1Oth order polynomial > plot(x, y,o, xi, z,:, xi, zz) plot data >)xlabel(x), ylabel(y=f(x)), title( 2nd and 10th Order curve Fitting) 在下面的图112中,原始数据标以'o,2阶曲线拟合是虚线,10阶拟合是实线。注意, 在10阶拟合中,在左边和右边的极值处,数据点之间出现大的纹波。当企图进行高阶曲线 拟合时,这种纹波现象经常发生。根据图11.2,显然,‘越多就越好’的观念在这里不适 用 2nd and 10th Order curve Fitting 10 04 0.6
-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.4700e-001)和最大(5.8233e+006)系数之间有 7 个数量级的幅度差。将这个解作图，并把此 图与原始数据及 2 阶曲线拟合相比较，结果如何呢？ » zz=polyval(pp, xi); % evaluate 10th order polynomial » plot(x, y, ' o ' , xi, z, ' : ' , xi, zz) % plot data » xlabel(' x '), ylabel(' y=f(x) '), title(' 2nd and 10th Order curve Fitting ') 在下面的图 11.2 中，原始数据标以'o'，2 阶曲线拟合是虚线，10 阶拟合是实线。注意， 在 10 阶拟合中，在左边和右边的极值处，数据点之间出现大的纹波。当企图进行高阶曲线 拟合时，这种纹波现象经常发生。根据图 11.2，显然，‘ 越多就越好 ’的观念在这里不适 用。 0 0.2 0.4 0.6 0.8 1 -2 0 2 4 6 8 10 12 14 16 x y=f(x) 2nd and 10th Order curve Fitting
图1122阶和10阶曲线拟合 11.2一维插值 正如在前一节对曲线拟合所描述的那样,插值定义为对数据点之间函数的估值方法,这 些数据点是由某些集合给定。当人们不能很快地求出所需中间点的函数值时,插值是一个有 价值的工具。例如,当数据点是某些实验测量的结果或是过长的计算过程时,就有这种情况。 或许最简单插值的例子是 MATLAB的作图。按缺省, MATLAB用直线连接所用的数 据点以作图。这个线性插值猜测中间值落在数据点之间的直线上。当然,当数据点个数的增 加和它们之间距离的减小时,线性插值就更精确。例如, >)xI=linspace(0, 2*pi, 60); >)x2=linspace(0, 2* pi, 6): >plot(xl, sin(x1), X2, sin(x2),-") >xlabel(x), ylabel( sin(x)), title(' Linear Interpolation) Linear Interpolation 0.5 -0.5 X 图11.3线性插值 图113是sine函数的两个图,一个在数据点之间用60个点,它比另一个只用6个点更 光滑和更精确 如曲线拟合一样,插值要作决策。根据所作的假设,有多种插值。而且,可以在一维以 上空间中进行插值。即如果有反映两个变量函数的插值,z=f(x,y),那么就可在x之间和在
图 11.2 2 阶和 10 阶曲线拟合 11.2 一维插值 正如在前一节对曲线拟合所描述的那样，插值定义为对数据点之间函数的估值方法，这 些数据点是由某些集合给定。当人们不能很快地求出所需中间点的函数值时，插值是一个有 价值的工具。例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况。 或许最简单插值的例子是 MATLAB 的作图。按缺省，MATLAB 用直线连接所用的数 据点以作图。这个线性插值猜测中间值落在数据点之间的直线上。当然，当数据点个数的增 加和它们之间距离的减小时，线性插值就更精确。例如， » x1=linspace(0, 2*pi, 60); » x2=linspace(0, 2*pi, 6); » plot(x1, sin(x1), x2, sin(x2), ' - ') » xlabel(' x '), ylabel(' sin(x) '), title(' Linear Interpolation ') 0 1 2 3 4 5 6 7 -1 -0.5 0 0.5 1 x sin(x) Linear Interpolation 图 11.3 线性插值 图 11.3 是 sine 函数的两个图，一个在数据点之间用 60 个点，它比另一个只用 6 个点更 光滑和更精确。 如曲线拟合一样，插值要作决策。根据所作的假设，有多种插值。而且，可以在一维以 上空间中进行插值。即如果有反映两个变量函数的插值，z=f(x, y)，那么就可在 x 之间和在
y之间,找出z的中间值进行插值。 MATLAB在一维函数 interp和在二维函数 interp.2中, 提供了许多的插值选择。其中的每个函数将在下面阐述 为了说明一维插值,考虑下列问题,12小时内,一小时测量一次室外温度。数据存储 在两个 MATLAB变量中。 >> hours=1: 12: index for hour data was recorded temps=589152529313022252724l% recorded >plot(hours, temps, hours, temps, + view temperatures title( Temperature Temperature 图114在线性插值下室外温度曲线 如图114看到的, MATLAB画出了数据点线性插值的直线。为了计算在任意给定时 间的温度,人们可试着对可视的图作解释。另外一种方法,可用函数 interp >)tinterpl(hours, temps, 9.3) estimate temperature at hour=9.3 22.9000
y 之间，找出 z 的中间值进行插值。MATLAB 在一维函数 interp1 和在二维函数 interp2 中， 提供了许多的插值选择。其中的每个函数将在下面阐述。 为了说明一维插值，考虑下列问题，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(' Temperature ') » xlabel(' Hour '), ylabel(' Degrees Celsius ') 0 2 4 6 8 10 12 5 10 15 20 25 30 35 Hour Degrees Celsius Temperature 图 11.4 在线性插值下室外温度曲线 正如图 11.4 看到的，MATLAB 画出了数据点线性插值的直线。为了计算在任意给定时 间的温度，人们可试着对可视的图作解释。另外一种方法，可用函数 interp1。 » t=interp1(hours, temps, 9.3) % estimate temperature at hour=9.3 t = 22.9000