曲线拟合 教学设计

任务3曲线拟合知识目标：熟悉关联函数的选择和线性化方法，了解线性最小二乘法。能力目标：能用Excel将实验数据拟合成曲线。☆思考：何谓曲线拟合？数据处理的一项十分重要的工作是寻求相关量之间的内在规律，即由已知数据群确立经验或半经验的数学模型。常用的方法是将观测得到的离散数据标记在平面图上，这只是对一个变量情况而言，描成一条光滑曲线（也包括直线，或对数坐标下的直线等）。为了便于进一步分析运算，希望将曲线用一简单的数学表达式加以描述，这种将离散的数据描述成数学表达式的方法称曲线拟合，或者说经验建模。一、关联函数的选择和线性化实测数据关联成数学模型的方法一般有以下几种：①具有一定的理论依据，可直接根据机理选择关联函数的形式。如反应动力学方程通常表示为r=kCAn，其中反应速度常数k与温度T的关系符合阿累尼乌斯（S.A.Arrhenius）方程，k=k0exp(－E/RT)的形式。此法的关键在于确定上述公式中k、n、k0、E等未知系数，以使模型密切逼近实测数据。这种模型称为半经验模型，工作要点在于参数估计。②尚无任何理论依据，但已有一些经验公式可选择。很多物性数据如热容、密度、饱和蒸气压等与温度的关系常表示为：f(T)=b0+b1T+b2T2+b3T3+b4lnT+b5/T当然不一定上述公式中六个系数都很重要，有的物性也只取前三、四项即可满足精度要求，这样可使模型简单化。③没有任何经验可循的情况对于此类情况，通常只能将实验数据画出图形与已知函数图形进行比较，选择图形接近的函数形式作拟合模型。不论上述哪种情况，在选定关联函数的形式之后，就是如何根据实验数据去确定所选关联函数中的待定系数，最常用的方法是线性最小二乘法。此法可用于处理一元或多元的线性模型。☆思考：何谓最小二乘法？一元线性模型：Y=A+BX （1-3-1）多元线性模型：Y=B0+B1X1+B2X2+…… （1-3-2）对于一些非线性模型，应事先将其变换成线性形式，即线性化处理，然后再用线性最小二乘法进行关联。表1-3-1列出了化工中常用的几种函数类型及线性化的方法。此表中所列均为单变量问题，经线性化处理后的线性模型均可统一用式（1-3-1）表示。对于多变量函数关系y=f(x1,x2,…)若采用幂函数的乘积作为关联函数，即将上面的函数关系写成如下形式y=ax1bx2c… （1-3-3）可作如下线性化处理，令Y=lny，X1=lnx1，X2=lnx2，…，B0=lna，B1=b，B2=c，…经线性化处理后的模型即式（1-3-2）。对于一元非线性化方程，如：y=a+bx+cx2+dx3… （1-3-4）令Y=y，X1=x，X2=x2，…，B0=a，B1=b，B2=c，…经线性化处理后的模型也为式（1-3-2）。表1-3-1常用函数线性化方法图形函数及线性化方法图形函数及线性化方法幂函数y=axb令Y=lgyX=lgxA=lgaB=b则Y=A+BX对数函数y=a+blgx令Y=yX=lgxA=aB=b则Y=A+BX指数函数y=aebx令Y=lnyX=xA=lnaB=b则Y=A+BX双曲函数1/y=a+b/x令Y=1/yX=1/xA=aB=b则Y=A+BX负指数函数y=aeb/x令Y=lnyX=1/xA=lnaB=b则Y=A+BXS型曲线函数y=a+be-x令Y=yX=e-xA=aB=b则Y=A+BX二、线性最小二乘法关联函数的形式确定之后，如何由实验数据比较精确地去确定关联函数中的待定系数仍是一个重要问题，最常用的方法就是线性最小二乘法。今有一弹簧秤，用它称重，记下荷重与弹簧伸长的关系如表1-3-2所示。表1-3-2弹簧荷重与弹簧伸长的关系荷重xi(kg)0246810121416长度yi(cm)30.0031.2532.5833.7135.0136.2037.3138.7940.04将数据点画在图1-3-1上，可以看出荷重与伸长两者大致呈直线关系。但并不严格在一条直线上，说明由于读数或其他影响因素造成数据包含有随机误差。根据力学上的虎克定律．弹簧伸长y应该与荷重x成正比，即y是x的线性函数，通过实验确定比例系数（弹簧的弹性系数）。一般地直线方程模型表示为y¢＝a+bx （1-3-5）如果用直尺将图1-3-1上的点连成直线，由于9个点不在一直线上，所以可以画出多条直线。也即式（1-3-5）线性模型中参数a和b可以有多种取值，于是产生这样一个问题，图1-3-1众多的连线中哪一条直线最能体现物理现象的本质呢?换句话说线性模型式（1-3-5）中截距a和斜率b取什么值为最佳选择？为说明这个问题，这里引入“残差”的概念。图1-3-1荷重与弹簧伸长长度关系设有n对实验数据（xi，yi）（i=1，2，…，n），需要寻找一个近似函数模型y¢=f（x）来拟合这一组数据。令第i点实测函数值yi与模型计算值y¢i之差为残差，即di=yi－y¢i=yi－f（x） （1-3-6）显然，di刻划了yi与回归模型计算值y¢i的偏离程度。如果每一个点的残差di=0，说明实验数据（xi，yi）完全可用直线拟合，但出于存在实验误差，di=0是不可能的。也就是说最佳的a和b应使的和最小。但用的和最小原则估计参数a和b，在应用上不很方便，所以，一般采用最小二乘法，其原理可以这样描述：所谓最小二乘原理就是使残差的平方和最小，即 （1-3-7）用最小二乘原理选择最佳拟合模型的物理意义是显见的，即在上例中找一条直线，使它与各实测点的距离（即di）平方加和最小。将式（1-3-5）代入式（1-3-7），可得显然，Q是a和b的函数。由数学分析多元函数求极值的必要条件，使Q最小的a、b必须满足以下方程组 （1） （2）由式（1）得 （3）其中分别表示yi、xi的平均值。由式（2）可推得 （4）由式（3）、式（4）经整理可得回归系数b的计算公式 （1-3-8）或 （1-3-9）具体计算时，先由式（1-3-8）或（1-3-9）求得b，再代入式（3）得a （1-3-10）由于计算a、b的公式中所有的量都可以从观测数据得出，因此回归直线方程y¢＝a+bx便可确定。为了简化公式，将上述公式中的（下同）均用å代替。令任一数据点xi与其平均值之差称为离差（xi－），xi的离差的平方和记为lxx，即 （1-3-11）同样，将yi的离差的平方和记为lyy，即 （1-3-12）将xi的离差与yi的离差的乘积之和记为lxy，即 （1-3-13）用式（1-3-11）、（1-3-12），则式（1-3-9）可表示为 （1-3-14）利用式（1-3-14）、式（1-3-10）计算表1-3-2中的数据，列于表（1-3-3）中。lxx=816－(72)2/9=240lxy=2668.58－72×314.89/9=149.46b=lxy/lxx=149.46/240=0.62275a=314.89/9－0.62275×72/9=30.006表1-3-3弹簧荷重与弹簧伸长的关系及计算数据荷重xi(kg)0246810121416åxi=72长度yi(cm)30.0031.2532.5833.7135.0136.2037.3138.7940.04åyi=314.89xi204163664100144196256åxi2=816xiyi062.50130.32202.26280.08362.00447.72543.06640.64åxiyi=2668.52yi2900.00976.561061.461136.361225.701310.441392.041504.661603.2011110.42故回归方程为y¢＝30.006+0.62275x对于非线性模型作了线性化处理后的回归方程Y¢=A+BX需注意，此模型并非最终模型，最后需恢复线性化处理前的模型原样。例1-3-1某化学反应的速度常数k与绝对温度T的实验数据如表1-3-4第一、二列所示，试用线性最小二乘法进行关联。表1-3-4例1-3-1实验数据表Ti（K）ki（1/min）Xi=1/TiYi=lnkiXi2XiYi363.00.6680×10-22.7548×10-3－5.0097.5890×10-6－1.3798×10-2373.01.3760×10-22.6810×10-3－4.2867.1876×10-6－1.1491×10-2383.02.7170×10-22.6110×10-3－3.6066.8171×10-6－9.4142×10-3393.05.2210×10-22.5445×10-3－2.9526.4746×10-6－7.5127×10-3403.09.6630×10-22.4814×10-3－2.3376.1573×10-6－5.7987×10-3å19.645×10-213.073×10-3－18.19034.226×10-6－4.8014×10-2解：由反应动力学知，反应速度常数与绝对温度的关系一般服从阿累尼乌斯（S.A.Arrhenius）方程，即k=k0exp(－E/RT)的形式。由于此式为非线性函数，需进行线性化处理，两边取对数后得令Y=lnk，X=1/T，A=lnk0，B=－E/R则Y=A+BX将原k~T数据组换算成X~Y组列于表1-3-4第三、四列，并计算出Xi2、XiYi列于第五、六列，将相关数据代入式（1-3-8）得代入（1-3-10）得A=－18.190/5－(－9771.68)×13.073×10-3/5=21.910∵A=lnk0∴k0=3.2779×109∴原关联方程为k=3.2796×109exp(－9771.68/T)☆思考：如何采用Excel进行曲线拟合？三、采用Excel曲线拟合用Excel进行曲线拟合的详细处理方法可参阅有关Excel方面的专著，以下通过两个具体例子加以说明。例1-3-2将表1-3-2所示的弹簧荷重与弹簧伸长之间的关系实验数据，用Excel将其拟合成函数关系。步骤1：打开Excel，将表1-3-2中的数据按列依次输入，如图1-3-2所示左边数据表；图1-3-2例1-3-2数据表及拟合数据散点图步骤2：选中表中数据区，即单元格A31~B39范围，然后选择“插入\图表”，再选择没有线的“XY散点图”，按照图表向导的步骤添加图的标题“例1-3-2直线拟合”、x轴、y轴等，再将所生成的散点图放在同一电子表格内，如图1-3-2右边所示。步骤3：右击图中数据点，选择“添加趋势曲线”，在“类型”中选择“线性”，在“选项”中选择“显示公式”、“显示R平方项”，按“确定”，其结果如图1-3-3所示。步骤4：由图1-3-3拟合直线的上方可得拟合的直线方程和相关系数R2之值，即：y=30.006+0.62275x、R2=0.9995可见，拟合结果与手工计算的结果完全相符。图1-3-3例1-3-2数据表及最终拟合直线图例1-3-3如表1-3-5所示为20℃下SO2在水中的溶解平衡数据，试拟合其函数关系。表1-3-520℃下SO2在水中的溶解平衡数据103x0.1410.2810.5620.8431.401.962.804.20102y0.1580.4211.121.863.425.137.7612.1解：因为SO2在水中属于中等溶解度气体，不能简单地应用亨利定律表示其平衡关系。用表中实验数据作出草图可知该曲线呈幂函数y=axb。处理的具体步骤与例1-3-2相同，所不同的是在步骤3的“类型”中选择“乘幂”，其最终的结果如图1-3-4所示。图1-3-4例1-3-3数据表及最终拟合曲线图从图中即可得回归结果为：y=144.45x1.2767、R2=0.9976四、知识拓展——曲线拟合效果分析（一）线性相关系数与显著性检验需要指出，曲线拟合处理的是随机变量问题，观测值x与y不存在确定性函数关系，而只是—种相关关系。线性最小二乘法只适宜处理变量x与y具有相关的问题，但在线性最小二乘法应用过程中，并不需要限制两个变量之间一定具有线性相关关系，就是说即使平面图上一堆完全杂乱无章的散点，也可用此方法给它们配一条直线方程模型。显然这样做是毫无意义的。只有当两个变量大致呈线性关系时才适宜用直线模型去拟合数据，因此必须给出一个数量性指标描述两个变量线性关系的密切程度，该指标称相关系数，通常记作r，其表达式为 （1-3-15）图（1-3-5）说明了r取各种不同数值时散点的分布情况。图1-3-5r不同时散点分布情况①r=0，此时lxy=0，因此b=0，即根据最小二乘法确定的回归直线平行于x轴，这说明y的变化与x无关，此时x与y毫无线性关系，通常这时散点分布是完全不规则的。如图（1-3-5）中的（A）。②0＜çr÷＜1，这是绝大多数情形，x与y之间存在一定的线性关系。当çr÷越接近于1，说明线性相关越大，也就是散点与回归直线越靠近。当r＞0时，b＞0，y随x增加而增加，称为正相关。当r＜0时，b＜0，y随x增加而减小，称为负相关。如图1-3-5中的（B）、（C）。③çr÷=1，所有数据都在回归直线上，此时，x与y完全相关，实际上此时x，y间存在确定的线性函数关系。如图1-3-5中的（D）。利用式（1-3-15）对表（1-3-2）中弹簧荷重与弹簧伸长长度之间关系的回归方程进行线性相关系数的计算。由式（1-3-13）得lyy=11110.42－(314.89)2/9=93.124由此可见，此例变量间线性相关程度很好。必须指出，相关系数只表示x与y的线性关系的密切程度，当r很小或为零时，并不表示x与y不存在其他关系。如图1-3-5中的（E），x，y呈某种曲线关系。可以对它进行线性化处理，变换成为Y=A+BX直线方程模型，此时可以用相关系数来讨论新变量X与Y之间线性相关程度，但是，新变量X与Y线性相关程度并不能直接说明原始数据x，y与非线性模型拟合效果的优劣。因此，对于非线性模型拟合的效果常用另一指标——相关指数来衡量，记作R2 （1-3-16）式中yi——为未经线性变换的原始数据；yi¢——是非线性模型的计算值；——是原始数据的平均值。显然R2＜1，R2值越接近于1，拟合曲线效果越好，当R2=l时，说明yi与yi¢趋于—致，实测点完全落在拟合曲线上。（二）相关系数r与显著性水平a对例1-3-1关联式进行相关指数R2的计算，以对关联式的拟合效果进行评价。算出各个Ti值所对应的ki¢，再算出å(ki－ki¢)2及å(ki－)2，列于表1-3-6中。表1-3-6例1-3-1计算相关指数R2的数据表Ti（K）363373383393403åki×102（min-1）0.66801.37602.71705.22109.663019.645ki¢×102（min-1）0.66791.37452.72405.21389.6629(ki－ki¢)2×109（min-2）6.6841040.21874.94715.15721.15510310.325(ki－)2×104（min-2）10.6346.51781.46891.669332.87953.169由式（1-3-16）得R2=1－10.325×10-9/(53