《化工数值计算》PPT课件.ppt

《《化工数值计算》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《化工数值计算》PPT课件.ppt（108页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、回归分析和曲线拟合生产过程和科学实验中,常用的变量大体可分两类。一类为确定性变量,另一类为随机变量。确定性变量是指两个或多个变量之间有确定的关系.即其中某个变量的每个值,都与一变量的一个或几个完全确定的值相对应,即它们之间存在着函数关系:例如,理想气体的压力P与摩尔体积V间,存在着确定的函数关系:但在实际问题中,由于变量之间的关系比较复杂,或由于生产或实验过程中不可避免地存在着误差,使变量之间的关系具有不确定性,也就是说,某个变量对应的,不是一个或几个确定的值,而是整个集合的值,这时,变量x和y间的关系,就称为相关关系。例如,流体在圆形直管中做湍流时的情形,通过量纲分析可知,努塞尔特数Nu、普

2、兰特数Pr和雷诺数Re之间存在着如下相关关系:这种关系的不确定性,表现为式中a和b的数值,在每次测量中不尽相同。不确定的原因,首先是影响该过程的因素甚多,有些因素至今尚未弄清;其次是受到实验过程中的偶然因素影响。这种不确定性关系并不说明上述三个量纲为1的数群之间无规律可循。相反,通过大量试验,人们发现,a和b的数值总是围绕着某一定值波动,而且随着试验次数的增多,a、b的数值趋于稳定。a、b的稳定值,可作为a和b的最佳估计值。在一定条件下,a=0.023,b=0.8。由此可见,通过大量试验,是可以找到隐藏在随机性后面的统计规律性的。回归分析和曲线拟合是一种处理变量相关关系的数理统计方法。用它可以

3、寻找隐藏在随机性后面的统计规律性。函数与相关是两种不同类型的变量关系,它们之间并无严格界限。一方面,相关的变量之间,并无确定的关系,但在一定的条件下,从一定的统计意义上看,它们之间又可能存在着某种确定的函数关系。另一方面,由于实际测定的数据中,总存在着误差,即使是确定性变量,也会出现某些非确定性结果。6.1一元线性回归一元线性回归处理的是两个变量之间的线性关系。所用的数学模型为一元线性代数模型,其模型方程式是对这种模型参数的估计,就是根据原始数据点(x1,y1)、(x2,y2)、(xi,yi)、(xn,yn),确定式(6-1)中a、b的估计值。在实际体系中,自变量x与因变量y之间服从线性关系的

4、情况虽然不多,但在不少情况下,x、y之间存在着某种函数组合关系。例如f1(x,y),f2(x,y),设两个函数之间服从线性关系f1与f2是不含待定系数的已知函数。若把f1(x,y)与f2(x,y)分别视为自变量与因变量,则仍可以借用线性模型去估计其参数值。这种方法称为化直法。它在化学化工的实际问题中是常见的。例如单分子基元反应AB的动力学方程式为对上式积分得式中,cA-t是不呈线性关系的函数。若对方程两边取对数,上式可化为lncA-t的线性函数:又例如,按照阿仑尼乌斯定律,反应速率常数k与温度T之间不呈线性关系:但lnk与1/T则呈线性关系:这些都是属于可化为线性关系的例子。一元线性代数模型中

5、的待定参数a和b,称为“估计值”。之所以称为“估计”值,是因为a,b的值是从实验值中通过数理统计方法确定的。图6-1一元线性回归6.1.1方法概述设有一组实验数据(x1,y1)、(x2,y2)、(xn,yn),自变量x与因变量y存在着式(6-1)的关系。当x取值为xi时,y的测定值为yi,计算值为yi*,并有由于参数a,b为未知值,故yi*也是未知值。若将全部实验数据标绘在x-y图中(见图6-1),由于各种因素的影响,它们不会全部落在一条直线上,即n个yi不会与n个yi*完全重合,它们将随机地分布在与xi呈线性关系的yi*的周围。以i表示它们之间的差值,则有这里i就是误差。它反映了xi使yi偏

6、离直线的各种影响因素的总和。现在,要寻找一条最靠近各个数据点的直线,这条直线称为回归直线。由于回归直线是一切直线中最接近各数据点(xi,yi)的,用它代表x与y之间的线性关系,比任何其他直线更为可靠。究竟如何确定回归曲线中的参数a和b呢?目前最常用的方法就是最小二乘法,即残差平方和最小法。式(6-3)中的误差i又称为残差,表示第i个数据与回归直线的偏离程度,则残差平方和Q表示全部数据与回归直线的总偏离程度。显然Q是a和b的函数:不用残差和i的原因是i有正有负,相加时可能彼此抵消,从而不能反映总的偏离程度,而用残差的平方和不会发生这种现象。由多元函数的极值理论可知,要使Q值最小,a、b必须满足下

7、列条件:即得式(6-6)称为一元线性回归的正规方程组,通过求解该方程组,可得:式(6-7)中等号右侧的量全部取自原始数据。因此,就可以确定回归系数a和b,完成参数估计。为了简化a和b的表达式,定义:式中,、分别为xi和yi的平均值。xi与之差(xi-),称为xi的离差;全部xi的离差平方和,称为x的离差平方和,记为Lxx:yi与之差(yi-),称为yi的离差;全部yi的离差平方和,称为y的离差平方和,记为Lyy,同理再令Lxy为全部xi的离差与yi的离差乘积的总和:将以上关系式代入式(6-7),得由式(6-12)第二式可以看出,回归直线是通过点(,)的。从力学观点看,(,)相当于n个实验点(x

8、i,yi)的重心,回归直线是通过重心的。应当指出:残差i只用yi-y*i表示时,表明yi有测量误差,而xi无测量误差;或表示与yi相比,xi的误差很小。因此,测量误差使实验点偏离回归直线,都表现为yi偏离y*i。如果xi的误差与yi的误差相比,不可忽略,则两者都必须考虑。这种情况比较复杂,此处不予介绍。求回归方程的计算过程中,不需要事先假定两个变量之间必须有相关关系。即使是一组杂乱无章的数据,也可以用最小二乘法绘制一条直线,以表示x与y的关系。显然,这种情况下,绘制的直线并无实际意义。为了判断两个变量间线性关系的优劣程度,引入一个新的指标R,称为简单相关系数,它的定义为R值不同时,数据点的分布

9、情况如下。(1)R=0图6-2R=0的数据点分布此时Lxy=0,b=0。即回归直线平行于x轴,y的变化与x无关,表示数据点的分布是无规则的,如图6-2所示。但亦有当R=0时,x与y确实存在明显相关性的情况。这种情形,不能应用线性回归方法,只能用化直线法或曲线拟合法处理。(2)0|R|0时,b0。数据点的y值随着x增加而增加,这种情况称为x与y正相关。R0时,b0。数据点的y值随着x增加而减小,这种情况称为x与y负相关。R的绝对值越小,数据点沿回归直线越分散。图6-301的数据点分布1的数据点分布(3)|R|=1x与y完全相关。全部数据点均落在回归直线上。若x与y为非线性相关,但经变量变换后,用

10、回归直线的方法处理,所求得的回归系数仅对变换后的变量是最佳的,而对原变量来说则并非最佳,但通常还能令人满意,此时应注意原变量的残差平方和并非最小。由以上讨论可知,相关系数R的绝对值在0与1之间,而且越接近于1,其线性关系越密切,那么|R|与1接近到什么程度,才能说明x与y之间存在线性相关关系呢?要回答这个问题,就要对相关系数进行显著性检验。由于篇幅所限,有关相关系数的显著性检验和回归方程的方差分析等问题将不在此讨论。如有需要,可参考有关数理统计方面的书籍。6.1.2程序框图图6-4是一元线性回归的通用计算程序框图。程序框图中的主要变量:N数据点数X、Y一维数组,用于存放原始数据中的x和y值XX

11、Lx离差平方和LxxYYLy离差平方和LyyXYLx离差与y离差乘积总和LxyA回归直线截距aB回归直线斜率bR简单相关系数6.1.3计算实例例6-1已知某反应的速率常数k与热力学温度T的实验数据如下:试求k-T的关系式。解通常反应速率常数与热力学温度的关系,服从阿仑尼乌斯定律:k=Ae-E/RT式中,E为反应活化能;T为热力学温度;R为气体通用常数。上式取对数,且令y=lnk,x=-1/T,可得y=lnA+x。按图6-4,用一元线性回归求得A=1.966109/min,E=79.571kJ/mol。将实验数据点和利用关系式获得的计算点一起绘制在图6-5中。源程序:*Example6-1-Eg

12、6-1.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(50),Y(50)DimXYAAsVariantClsN=5XYA=Array(363,0.00718,373,0.01376,383,0.02701,_393,0.05221,403,0.09718)K=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+2X(I)=-1/X(I)Y(I)=Log(Y(I)NextICallLINEAR1(N,X(),Y(),A,B,R)A=Exp(A):E=B*8.314PrintA=;A:PrintE=;EPrintR=;R

13、EndSub*SubLINEAR1(N,X(),Y(),A,B,R)*XT=0:YT=0:XX=0:YY=0:XY=0ForI=1ToNXT=XT+X(I):YT=YT+Y(I)XX=XX+X(I)*X(I):YY=YY+Y(I)*Y(I)XY=XY+X(I)*Y(I)NextIXXL=XX-XT*XT/NYYL=YY-YT*YT/NXYL=XY-XT*YT/NB=XYL/XXLA=(YT-B*XT)/NR=XYL/Sqr(XXL*YYL)EndSub执行结果:A=1966349283.054212(指前因子)E=79570.97618674007(活化能)R=.999718315533107

14、(相关系数)源程序中将一元线性回归计算安排在子程序LINEAR1中。例6-2某水样BOD测定数据如下:试确定该水样中有机物生物氧化降解反应的经验速率方程表达式。解通常水体中有机物生物氧化降解反应的经验速率方程服从下列方程:式中,BOD、L0为分别为t时和初始时刻的生化需氧量(mg/L);k为BOD的降解系数,即耗氧系数/d-1托马斯(Thomas)提出将1-e-kt按幂级数展开如下:与此展开相似的表达式有:两展开式仅第四项出现微小差别,故可以近似地取即BOD=L0整理得=+t取y=、x=t按一元线性回归计算a=0.25778、b=0.01371,从而解得k=6b/a=0.31907d-1,L0

15、=182.963mg/L。源程序(仅列出主程序,回归子程序LINEAR1同例6-1):*Example6-2-Eg6-2.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(50),Y(50)DimXYAAsVariantClsN=10XYA=Array(1,58,2,85,3,107,4,125,5,138,_6,147,7,155,8,161,9,167,10,170)K=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+2Y(I)=(X(I)/Y(I)(1/3)NextICallLINEAR1(N,X(),Y(),

16、A,B,R)PrintA=;A:PrintB=;BBK=6*B/ABL0=1/BK/A/A/APrintL0=;BL0:PrintK=;BKPrintR=;REndSub执行结果:A=.2577799110470602B=1.370838615288584D-02L0=182.9634248033556K=.3190718647672257R=.9900105997750359此外,塞里奥特(Theriaut)提出了BOD公式的另一种解法:式中,k为待估的k的近似值;h为k的允许偏差量。从而有因h甚小,e-ht1-ht,故上式变为:式中,a=L0,d=L0h,x1=1-e-kt,x2=te-k

17、t。这样可以首先假设k的初始值为k0,利用实验数据通过二元线性回归(见下节例6-5)确定出a和d,并求出L0=a、h=d/L0;若|h|(误差允许值),则令k=k0+h,重新进行线性回归计算,直至|h|m+1才能求出上式中的m+1个回归系数。同样由多元函数的极值理论可知,要使Q值最小,a0和aj必须满足下列条件:式(6-15)经整理可得:式(6-16)称为多元线性回归模型的正规方程组。它是一个m+1元的线性代数方程组。由于xij和yi已知,故可求得m+1个待定系数a0,a1,am。实际计算时,一般作如下处理:先将式(6-16)的第一式写成然后将式(6-17)代入方程组(6-16)的第2至第m+

18、1式,重新组成一个m元线性方程组,其中有a1,a2,am等m个待定系数。通过求解此m元线性方程组,获得系数a1,a2,am,再代回式(6-17),求得a0。为简化计算,用表示第j个x的平均值,表示y的平均值,则用Ljk表示第j个x离差与第k个x离差乘积之和,则用Lyy表示y离差的平方和,则用Ljy表示第j个x离差与y离差乘积之和,则将式(6-17)分别代入式(6-16)的第2至m+1式,经简化整理可得如下m元线性方程组:可用主元素消去法求解此式,然后将求得的a1,a2,am代入式(6-17),求出a0,从而完成对多元线性回归模型的参数估计。多元线性回归的计算中,常用复相关系数衡量数据点之间的线

19、性优劣。复相关系数定义如下:式中,U称为回归平方和:应当指出,并非所有曲线都可以按这种方法处理。例如抛物线就不能通过变量变换把它化为直线。但是如果令x1=x,x2=x2,则上式就化成一个包含两个自变量的线性方程从而将抛物线按二元线性回归计算。对于含多变量的任意多项式也可以通过类似的变换,把它们转化成多元线性回归计算。6.2.2程序框图图6-6是多元线性回归的通用计算程序框图。图6-6(a)多元线性回归的通用计算程序框图(1)图6-6(b)多元线性回归的通用计算程序框图(2)程序框图中的主要变量:N数据点数M多元线性模型元数X二维数组,用于存放原始数据的x值Y一维数组,用于存放原始数据的y值YP

20、值YYLLyy值XP一维数组,用于存放值A二维数组,用于存放m元线性方程组的系数LjkB一维数组,用于存放m元线性方程组的常数项LjyC一维数组,用于存放多元线性模型的系数aj(j=0,1,M)R复相关系数R0U回归平方和Q残差平方和子程序XYF为列主元消去法求解线性方程组的程序,可参见图5-2和图5-3。6.2.3计算实例例6-4已知某溶液由两种物质组成,cA为物质A的浓度(g/L),cB为物质B的浓度(g/L),为溶液的黏度(mPas)。设数学模型为=a0+a1cA+a2cB试根据下列实验数据,确定a0、a1、a2的值。解按图6-6编写计算源程序。源程序:*Example6-4-Eg6-4

21、.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(100,20),Y(100),C(20)DimXYAAsVariantClsN=15:M=2XYA=Array(25.8,98,14.5,15.8,116,9.7,18.1,104,11.3,_13.3,99,26,20.1,153,44.7,10.1,98,21,_17.1,103,25.2,21,112,13.7,23.7,113,38.5,_11.2,80,5.8,10.2,87,17.7,16.4,138,40,_15.9,98,17.1,8,102,3,26,155,37.3)K=0Fo

22、rI=1ToNForJ=1ToM:X(I,J)=XYA(K):K=K+1:NextJY(I)=XYA(K):K=K+1NextICallLINEAR2(N,M,X(),Y(),C(),R)PrintTab(4);*Results*ForJ=0ToMPrintA(;J;)=;Format$(C(J),#.#)NextJPrintR=;Format$(R,#.#)EndSub*SubLINEAR2(N,M,X(),Y(),C(),R)*DimA(20,20),B(20),XP(20)YP=0=yAverageForI=1ToN:YP=YP+Y(I):NextIYP=YP/NYYL=0=LyyAve

23、rageForI=1ToN:YYL=YYL+(Y(I)-YP)*(Y(I)-YP):NextIForJ=1ToM=XjAverageXP(J)=0ForI=1ToN:XP(J)=XP(J)+X(I,J):NextIXP(J)=XP(J)/NNextJForJ=1ToM=LjyXYL=0ForI=1ToNXYL=XYL+(X(I,J)-XP(J)*(Y(I)-YP)NextIB(J)=XYLNextJForJ=1ToM=LjkForK=1ToMXXL=0ForI=1ToNXXL=XXL+(X(I,J)-XP(J)*(X(I,K)-XP(K)NextIA(J,K)=XXLNextKNextJCal

24、lXYF(A(),B(),M,C()C(0)=YP=a0ForJ=1ToMC(0)=C(0)-C(J)*XP(J)NextJU=0:Q=0=RForI=1ToNYI=C(0)ForJ=1ToMYI=YI+C(J)*X(I,J)NextJU=U+(YI-YP)*(YI-YP)Q=Q+(YI-Y(I)*(YI-Y(I)NextIR=Sqr(U/(U+Q)-orR=SQR(U/YYL)EndSub执行结果:*Results*A(0)=-27.4324958A(1)=0.2327103A(2)=0.4095299R=0.7472224源程序中将多元线性回归安排在子程序LINEAR2中,子程序XYF和X

25、LZY见例5-3。例6-5按塞里奥特方法确定例6-2的BOD公式中的参数。解由前叙述可知,塞里奥特方法计算机求解时要应用二元线性回归,且无常数项,即a00。对a00时的多元线性回归,求解回归方程系数的m元线性方程组同样具有方程组(6-18)的形式,不过系数矩阵和常数项必须用下列式子进行计算:将例6-4的源程序加以修改即可得到本例的计算源程序。源程序(子程序XYF和XLZY同例6-4):*Example6-5-Eg6-5.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(100,20),Y(100),C(20),T(100)DimXYAAsVaria

26、ntClsN=10:M=2XYA=Array(1,58,2,85,3,107,4,125,5,138,_6,147,7,155,8,161,9,167,10,170)K=0ForI=1ToNT(I)=XYA(K):Y(I)=XYA(K+1):K=K+2NextIE=0.000001:X0=0.5:H=100*EDoWhileAbs(H)EForI=1ToNX(I,1)=1-Exp(-X0*T(I)X(I,2)=T(I)*Exp(-X0*T(I)NextICallLINEAR20(N,M,X(),Y(),C(),R)CL0=C(1):H=C(2)/C(1)X0=X0+HLoopPrintTab(

27、4);*Results*PrintL0=;CL0:PrintK=;X0:PrintR=;REndSub*SubLINEAR20(N,M,X(),Y(),C(),R)*DimA(20,20),B(20),XP(20)YP=0=yAverageForI=1ToN:YP=YP+Y(I):NextIYP=YP/NForJ=1ToM=SjyXYL=0ForI=1ToNXYL=XYL+X(I,J)*Y(I)NextIB(J)=XYLNextJForJ=1ToM=SjkForK=1ToMXXL=0ForI=1ToNXXL=XXL+X(I,J)*X(I,K)NextIA(J,K)=XXLNextKNextJC

28、allXYF(A(),B(),M,C()U=0:Q=0=RForI=1ToNYI=0ForJ=1ToMYI=YI+C(J)*X(I,J)NextJU=U+(YI-YP)*(YI-YP)Q=Q+(YI-Y(I)*(YI-Y(I)NextIR=Sqr(U/(U+Q)EndSub执行结果:*Results*L0=173.428617756483(mg/L)K=0.3306173539389562(d-1)R=0.9956168注意按这种方法进行计算时,k的初始值必须在真实值附近选取,否则将得出错误的结果。注意本例源程序中的多元线性回归子程序LINEAR20,只适用于无常数项的多元线性方程的回归。6.

29、3剔除可疑数据及其计算程序6.3.1剔除可疑数据的方法在线性回归计算中,假定每个测定数据与回归结果之间的误差均在随机误差允许的范围之内。然而,由于测量误差或过失误差等多种原因,在一组实验值中,误差往往会超出随机误差的允许范围。这些数据,称为可疑数据。为保证回归结果的可靠性,必须剔除这些可疑的数据。剔除可疑数据,应当有一个科学的标准。这个标准就是统计判据,属于统计判据的剔除准则有多种。以一元线性回归为例,其代数模型为y=a+bx。若自变量x无测量误差,则y的标准偏差为式中,n为原始数据点数;m为回归模型中自变量的个数,对一元线性回归m=1;i为残差,即i=yi-a、b是按最小二乘法求出的最佳估计

30、值。根据数理统计分析,合理的数据,其残差不应超出的k倍。若取k=3,便是常用的3准则。据此,可以把残差绝对值超过3的个别数据(xi,yi),判为可疑数据而加以剔除。必须指出,3准则是以数据点数n 为前提的,当n为有限值时,3判据并不十分可靠。下面介绍一种广泛采用的判据,即所谓肖维奈特准则。按肖维奈特准则,若n次等精度测量中,有某个测量值yi,其残差的绝对值超出k,就可以认为是可疑数据而予以剔除。表6-1列出了肖维奈特准则中与n相对应的k值。表6-1肖维奈特准则的n和k值nknknknk51.65152.13252.33802.7461.73162.16262.341002.8171.79172

31、.18272.351502.9381.86182.20282.371853.0091.92192.22292.382003.02101.96202.24302.392503.11112.00212.26352.455003.29122.04222.28402.5010003.48132.07232.30502.5820003.66142.10242.32602.6450003.89使用这个准则时,可根据回归结果,对全部实验值进行逐级检查,把属于可疑数据的实验值选出。若发现不止一个可疑数据,则应把其中残差绝对值最大者剔除,然后重新计算值。根据新的值,再次用肖维奈特准则进行检查。每次只剔除一个可疑

32、数据,其余数据重新进行回归,直至回归所用的数据中不再含有可疑数据为止。6.3.2剔除可疑数据的计算程序框图图6-7是具有剔除可疑数据功能的一元线性回归通用计算程序框图,整个计算过程分为输入原始数据、一元线性回归计算、确定肖维奈特准则的k值、确定残差绝对值最大的数据点、剔除最可疑数据点(即残差绝对值最大的数据点)。图6-7具有剔除可疑数据功能的一元线性回归通用计算程序框图程序框图中的主要变量:N原始数据点数或剔除可疑数据后的合格数据点数N1可疑数据点数X一维数组,用于存放原始数据及合格数据中的x值Y一维数组,用于存放原始数据或合格数据中的y值X1一维数组,用于存放可疑数据点的x值Y1一维数组,用

33、于存放可疑数据点的y值A回归直线截距B回归直线斜率R简单相关系数SD标准偏差ER平均相对误差DALTA绝对值最大的残差ID残差绝对值最大的数据点序号U肖维奈特准则的k值子程序LINEAR1A为一元线性回归计算子程序,比例6-1中的子程序LINEAR1增加了标准偏差和平均相对误差的计算;子程序RULES为肖维奈特准则中k值的计算程序。采用类似的方法,可以编写能剔除可疑数据的多元线性回归计算程序框图。6.3.3计算实例例6-6在定温下某反应的活化能(E)与压力(P)呈直线关系:E=a+bP由实验测得如下数据:试用肖维奈特准则剔除其中的可疑数据,并确定a和b的值,计算相关系数、标准偏差和平均相对偏差

34、。解按图6-7的程序框图编写源程序。源程序:*Example6-6-Eg6-6.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(500),Y(500),X1(500),Y1(500)DimXYAAsVariantClsN=10XYA=Array(1,40.2,2,80,3,40.9,4,41.6,5,41.8,_6,42.6,7,42.6,8,70,9,43.7,10,43.8)K=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+2NextI-IC=1:N1=0DoWhileIC=1CallLINEAR1A(N,X

35、(),Y(),A,B,R,SD,ER)CallRULES(N,U)DALTA=0:ID=0ForI=1ToNDT=Abs(Y(I)-A-B*X(I)IfDTDALTAThenDALTA=DTID=IEndIfNextIIfDALTAU*SDThenINN=0ForI=1ToNIfI=IDThenN1=N1+1X1(N1)=X(I):Y1(N1)=Y(I)ElseINN=INN+1X(INN)=X(I):Y(INN)=Y(I)EndIfNextIN=N-1:IC=1ElseIC=0EndIfLoop-PrintN=;N:PrintN1=;N1PrintI,X1,Y1ForI=1ToN1Print

36、I,X1(I),Y1(I)NextIPrintA=;A:PrintB=;B:PrintR=;RPrintSD=;SD:PrintER=;ER;%EndSub*SubLINEAR1A(N,X(),Y(),A,B,R,SD,ER)*XT=0:YT=0:XX=0:YY=0:XY=0ForI=1ToNXT=XT+X(I):YT=YT+Y(I)XX=XX+X(I)*X(I):YY=YY+Y(I)*Y(I)XY=XY+X(I)*Y(I)NextIXXL=XX-XT*XT/NYYL=YY-YT*YT/NXYL=XY-XT*YT/NB=XYL/XXLA=(YT-B*XT)/NR=XYL/Sqr(XXL*YYL

37、)SD=0:ER=0ForI=1ToNSD=SD+(Y(I)-A-B*X(I)*(Y(I)-A-B*X(I)ER=ER+Abs(Y(I)-A-B*X(I)/(A+B*X(I)NextISD=Sqr(SD/(N-2):ER=(ER/N)*100EndSubSubRULES(N,U)SelectCaseN执行结果:N=8N1=2IX1Y112802870A=39.80313111545988B=0.4172211350293531R=0.9910796148007959SD=0.1830558892195286ER=0.3321330871296535%从源程序的执行结果看,按肖维奈特准则,10个

38、数据点中有两个属于可疑数据,可疑数据的活化能测定值分别为80和70。6.4多项式拟合在化学化工的实验或科研中,经常需要从一组测定数据,例如从n对(xi,yi)数据,去求自变量x和因变量y的近似函数关系式y=p(x)。从图形上看,这是由给定的n个点(xi,yi)(i=1,2,n)作曲线拟合。在曲线拟合中,多项式拟合问题占特殊的地位。任何函数在一个比较小的范围内,可以用多项式任意逼近。因此,在比较复杂的实际问题中,可以不问y与各因素的确切关系,而用多项式拟合进行分析和计算。下面以多项式拟合为例,说明曲线拟合的方法和计算程序。6.4.1方法概述设用下列m次多项式:拟合一组数据(xi,yi)(i=1,

39、2,n),即曲线y=f(x)上已给定n个点,用多项式求作该曲线的近似图形。这一问题与前述的插值问题有类似之处。但插值问题要求近似曲线y=p(x)严格地通过所给的n个点,这一要求将会使近似曲线y=p(x)保留数据的全部测试点的测量误差。如果个别数据的误差很大,那么插值的效果显然是不够理想的。鉴于这种情况,考虑放弃严格通过所有结点(xi,yi)这一要求,而采用别的方法去构造近似曲线,以尽可能反映所给数据的总趋势。曲线拟合的常用方法仍然是最小二乘法,即残差平方和最小法。若以i代表结点处的残差,则残差的平方和为由于xi与yi为已知值,故Q是aj(j=0,1,2,m)的函数。由多元函数的极值理论可知,要

40、使Q最小,则系数ak必须满足下式:即(aj)=yi上式变换后得:令Sl=将上两式代入式(6-24)可得m+1元线性方程组:式(6-25)称为正规方程组。这是m+1元线性方程组,其系数矩阵是对称矩阵,可以证明,上述正规方程组有唯一解。所得的m次多项式Pm(x)确能使残差平方和Q最小,故Pm(x)即为所求的拟合多项式。它是函数y=f(x)的近似表达式,也可用它代替f(x)作微分、积分计算。此外,必须注意的是上式计算中1。6.4.2程序框图图6-8是多项式拟合的通用计算程序框图。图6-8多项式拟合的通用计算程序框图程序框图中的主要变量:M多项式次数N原始数据组数X、Y一维数组,用于存放原始数据的x、

41、y值S一维数组,用于存放正规方程组的Sl数值T一维数组,用于存放正规方程组各常数项Tk的值B二维数组,用于存放正规方程组的系数矩阵A一维数组,用于存放正规方程组的变量值,即多项式系数的值对图6-8需作如下几点说明。由于计算机在做乘幂运算时,是将它转化为对数而进行运算的,所以要考虑当自变量xi取为零时,对Sl及Tk的影响,即应考虑到否则计算机将指出发生运算逻辑错误。子程序XYF0为列主元消去法求解多元线性方程组的计算程序,注意与例6-4中子程序XYF稍有不同,本例中为m+1元线性方程组,而例6-4中为m元线性方程组。在实际计算工作中,当多项式的次数m7时,方程组容易出现病态,所拟合的曲线出现较大

42、的起伏。为了避免此种现象,m的数值不宜取得太大。6.4.3计算实例例6-7已知某化工厂产品的废品率y与该产品中的某种物质百分含量x有关。实际测定的数据如下:试用三次多项式拟合这组数据,并绘制图形。解按图6-8的程序框图,编写拟合的计算程序。源程序:(子程序XYF0见例4-6)*Example6-7-Eg6-7.frm*DefDblA-H,O-ZPrivateSubCommand1_Click()DimX(100),Y(100),A(10)DimXYAAsVariantClsN=16:M=3XYA=Array(0.34,1.3,0.36,1.0,0.37,0.73,0.38,0.9,_0.39,

43、0.81,0.39,0.7,0.39,0.6,0.4,0.5,_0.4,0.44,0.41,0.56,0.42,0.3,0.43,0.42,_0.43,0.35,0.45,0.4,0.47,0.41,0.48,0.6)K=0ForI=1ToNX(I)=XYA(K):Y(I)=XYA(K+1):K=K+2NextICallPOLYNOMIAL(N,M,X(),Y(),A()ForJ=0ToMPrintA(;J;)=;A(J)NextJEndSub*SubPOLYNOMIAL(N,M,X(),Y(),A()*DimS(20),T(10),B(10,10)ForL=0To2*M=S(L)S(L)=0

44、ForI=1ToNIfX(I)0ThenS(L)=S(L)+X(I)LElseIfL=0ThenS(L)=S(L)+1EndIfNextINextLForK=0ToM=T(K)T(K)=0ForI=1ToNIfX(I)0ThenT(K)=T(K)+Y(I)*(X(I)K)ElseIfK=0ThenT(K)=T(K)+Y(I)EndIfNextINextKForI=0ToM=A(I,J)ForJ=0ToMB(I,J)=S(I+J)NextJNextICallXYF0(B(),T(),M,A()EndSub执行结果:A(0)=-13.64297777274281图6-9例6-7的计算结果A(1)=155.8814079905283A(2)=-490.5083888981387A(3)=473.8607827598829用EXCEL办公软件将实验数据点和拟合的三次多项式曲线一起绘制于图6-9中。本例源程序中多项式拟合安排在子程序POLYNOMIAL中;若拟合无常数多项式,则用例5-3中子程序XYF代替子程序XYF0求解方程组即可。