实验数据的处理及模型参数的确定

实验数据的处理及模型参数的确定第一页，共五十四页，编辑于2023年，星期二实验数据的处理及模型参数的确定第二页，共五十四页，编辑于2023年，星期二●数学模型中各参数的确定例：镍硅藻土上苯加氢合成环己烷是表面反应控制的固体催化剂上的气相反应。在160oC，微分反应器中的初始反应速率方程为模型参数

ka

──表观速率常数

bH

──H2的吸附系数

bB──

C6H6的吸附系数利用实验得到的全部信息，确定数学模型中的待定参数第三页，共五十四页，编辑于2023年，星期二

●线性插值

●

Lagrange插值

●埃米尔特插值

●一元线性回归

●

线性模型的推广

●

多元回归

可化为多元线性回归的问题

●

多项式拟合简介

●逐次回归分析函数关系插值法回归分析相关关系数值微分★★★★★引言：2.常用的数学方法第四页，共五十四页，编辑于2023年，星期二例：72型分光光度计测得某试样的吸收值如下：2-1-1–1线性插值——问题的提出λ/nm430440450460470480A0.4100.3750.3250.2800.2400.205希望：根据给定的函数表作一个既能反应f(x)的特性，又便于计算的简单函数p(x)，用p(x)近似f(x)，计算出任意

x对应的y值求在435，445，455，465，475nm处的吸收值。第五页，共五十四页，编辑于2023年，星期二

定义：设y=f(x)在区间[a,b]上有意义，且已知在点a<x0<x1<…<xn<b上的值y0,y1,…,yn,若存在一简单函数pn(x),使

pn(xi)=yi(i=0,1,…,n)成立，则称pn(x)为f(x)的插值函数，x0，x1，…，xn为插值节点区间[a,b]为插值区间，求pn(x)的方法称为插值法

xyy=f(x)y=p(x)x1y1xn

yn几何意义：2-1-1–2线性插值——方法原理ab第六页，共五十四页，编辑于2023年，星期二线性插值原理：两点间直线方程：xyy=f(x)y=p(x)xi-1yi-1xi

yi2-1-1–2线性插值——方法原理第七页，共五十四页，编辑于2023年，星期二分段线性插值：实验点个数为n时，求插值结点x的函数值。首先确定x在哪两点间

2-1-1–2线性插值——方法原理第八页，共五十四页，编辑于2023年，星期二LINEPLOT(N,X,Y,X0,Y0)DOJ=1,N-1J1=J+1X0<=X(J1)CONTINUEJ=J-1T=(X0-X(J))/(x(J1)-x(J))Y0=Y(J)+T*(Y(J1)-Y(J))RETURNnoyes2-1-1–3线性插值——程序框图第九页，共五十四页，编辑于2023年，星期二开始输入：数据点X(I),Y(I)，未知点X0调用线性插值子程序求未知点X0对应的函数值Y0输出：X0,Y0值结束2-1-1–4线性插值——应用示例显示程序显示输入显示输出第十页，共五十四页，编辑于2023年，星期二2-1-2-1一元三点Lagrange插值——问题的提出例：计算乙醇的平均摩尔体积实验测得25℃时乙醇溶液的平均摩尔体积（cm2mol-1）与乙醇的物质的量分数的关系如下x/cm2mol-10.089121.220.115322.160.143523.180.173924.320.206825.570.242426.950.281128.470.323430.150.369732.010.420734.070.477136.37计算x=0.1,0.2,0.3,0.4时的。第十一页，共五十四页，编辑于2023年，星期二线性插值公式：二点(xi-1,yi-1),(xi,yi)(两点式)Lagrange插值（三点插值,抛物线插值）：xi-1xixi+1即2-1-2-2一元三点Lagrange插值——方法原理第十二页，共五十四页，编辑于2023年，星期二xyy=f(x)y=p(x)xi-1yi-1xi+1

yi+1xiyi编程难点：如何确定使用哪三个结点进行插值xj-1xj+1xjxj-2xj+22-1-2-2一元三点Lagrange插值——方法原理第十三页，共五十四页，编辑于2023年，星期二LGRG2(X,Y,N,T,Z)DoJ=3,N-1I=JT>X(I)CONTINUEP=(T-X(I))*(T-X(I+1))/(X(I-1)-X(I))/(X(I-1)-X(I+1))Q=(T-X(I-1))*(T-X(I+1))/(X(I)-X(I-1))/(X(I)-X(I+1))R=(T-X(I-1))*(T-X(I))/(X(I+1)-X(I-1))/(X(I+1)-X(I))Z=P*Y(I-1)+Q*Y(I)+R*Y(I+1)RETURNnoyesI=I-1|T-X(I-1)|<=|T-X(I)|yesno2-1-2-3一元三点Lagrange插值——程序框图第十四页，共五十四页，编辑于2023年，星期二开始输入：数据点X(I),Y(I)，未知点X0调用lagrange插值子程序求未知点X0对应的函数值Y0输出：X0,Y0值结束2-1-2-4一元三点Lagrange插值——应用示例显示程序显示输入显示输出第十五页，共五十四页，编辑于2023年，星期二2-2-1-1一元线性回归——问题的提出例:铜钼矿中钼对铜含量的线性依赖关系No123456789x(Cu)285290300303310318325335338y(Mo)4.64.74.74.94.95.15.05.35.4一元线性回归的数学模型：y=ax+b+εyi=axi+b+εin个实验点回归直线：y=ax+b残差：εi=yi-（axi+b）第十六页，共五十四页，编辑于2023年，星期二最小二乘法：第i点残差：εi=yi-（axi+b）当残差的平方和为最小时，对应的a、b值是最佳值。（正规方程组）2-2-1-2一元线性回归——方法原理第十七页，共五十四页，编辑于2023年，星期二令平均值离差平方和2-2-1-2一元线性回归——方法原理第十八页，共五十四页，编辑于2023年，星期二1.线性相关系数R—衡量回归方程式与数据相符合的程度。若R1，则数据点落在直线上。注意：2.加权最小二乘法3.剔除可疑数据2-2-1-2一元线性回归——方法原理第十九页，共五十四页，编辑于2023年，星期二PK(N,X,Y,A,B,R)SX,SY,SXX,SXY,SYY=0DoI=1,NX1=X(I),Y1=Y(I)SX=SX+X1,SY=SY+Y1SXX=SXX+X1*X1SYY=SYY+Y1*Y1,SXY=SXY+X1

LXX=SXX-SX*SX/n,LYY=SYY-SY*SY/n,LXY=SXY-SY*SX/NB=LXY/LXX,A=(SY-B*SX)/N,R=LXY/SQRT(LXX*LYY)

RETURN2-2-1-3一元线性回归——程序框图第二十页，共五十四页，编辑于2023年，星期二开始输入：数据点数N

铜与钼的实验数据X(I),Y(I)（I=1，N）调用一元线性回归子程序计算A，B，R输出：A，B，R结束2-2-1-4一元线性回归——应用示例显示程序显示输入显示输出第二十一页，共五十四页，编辑于2023年，星期二2-2-2-1线性模型的推广——方法原理变量x与y之间存在某种非线性关系确定曲线类型（非线性关系）实际经验散点图形状线性关系最小二乘法确定系数非线性关系第二十二页，共五十四页，编辑于2023年，星期二曲线类型及变换公式双曲线型幂指数型指数型S型对数型平方根曲线2-2-2-1线性模型的推广——方法原理第二十三页，共五十四页，编辑于2023年，星期二例1：Arrhenius公式的应用根据k和T数据，可确定指前因子A和活化能Ea。2-2-2-2线性模型的推广——应用示例第二十四页，共五十四页，编辑于2023年，星期二例2：Clausius-Clapryron方程式的应用纯组分气-液（气-固）两相平衡的方程式：上式中：p:T/K时液（固）饱和蒸气压；ΔH：相变热不定积分：测定不同温度下的饱和蒸气压，将lnp~1/T进行线性回归，可算出ΔH，并计算其它温度下的蒸气压2-2-2-2线性模型的推广——应用示例第二十五页，共五十四页，编辑于2023年，星期二开始输入：数据点数N

温度T与蒸气压p的实验数据T(I),P(I)（I=1，N）输出：A,B,ΔH,Ti,pi结束调用线性回归子程序计算A，B（相变热ΔH=-8.314E-3*B）T(I)=T(I)+273.15X(I)=1/T(I)Y(I)=ln[P(I)](I=1,N)计算其它温度Ti下的蒸气压pi2-2-2-2线性模型的推广——应用示例第二十六页，共五十四页，编辑于2023年，星期二例3：几种常用的吸附等温式的计算（气固吸附）2-2-2-2线性模型的推广——应用示例Freundlich经验式：Langmuir方程：B.E.T.方程：第二十七页，共五十四页，编辑于2023年，星期二2-2-3-1多元线性回归——方法原理数学模型：函数y与多个自变量间x1，x2，…，xm的线性相关关系设共进行了n次测定（i=1，2，…，n）的自变量取值分别为xi1，xi2，…，xim，函数值的测定值为yi，所得的值如表：第二十八页，共五十四页，编辑于2023年，星期二选择回归系数，以使残差的平方和最小，即残差残差的平方和Q残差的平方和最小，相当于求解以下方程组（i=1，2，…，n）（k=，0，1，…，m）最小二乘法2-2-3-1多元线性回归——方法原理第二十九页，共五十四页，编辑于2023年，星期二正规方程组：全相关系数R式中（m+1）元线性方程组，Guass消去法求解b0，b1，…，bm2-2-3-1多元线性回归——方法原理第三十页，共五十四页，编辑于2023年，星期二MLR(M,N,XV,S,B0,R)输入：X（I，J）I=1,2,…,NJ=1,2,…,M+1CALLGS(A,M,M1,EPS)RETURN程序框图DOI=1,MA(I,I)=1.0DOJ=I+1,M1A(I,J)=0DOK=1,NA(I,J)=A(I,J)+X(K,I)*X(K,J)DOI=1,M-1J=I+1,MA(J,I)=A(I,J)2-2-3-2多元线性回归第三十一页，共五十四页，编辑于2023年，星期二变量xi(i=1,2,…,m)与y之间存在某种非线性关系确定曲线类型（非线性关系）实际经验散点图形状线性关系最小二乘法确定系数非线性关系2-2-3-3多元线性回归——可化为多元线性回归的问题变量代换变量代换第三十二页，共五十四页，编辑于2023年，星期二函数式通式：y’=b0+b1x1+b2x2变换变量变换常数y’=x1=x2=b0=b1=b2=y=a+bx+cx2yxx2abcy=alnx+bcosx+cylnxcosxcaby2

=a+bx+cx2y2xx2abcy＝x/

(a+bx+cx2)x/yxx2abcy=a+b/x+c/x2y1/x1/x2abclnyxx2lnabc2-2-3-3多元线性回归——可化为多元线性回归的问题第三十三页，共五十四页，编辑于2023年，星期二例1：Antoine方程式的应用（p：蒸气压，T：温度）令Tlgp=y,T=x1,lgp=x2,b0=ac+b,b1=a,b2=c2-2-3-4多元线性回归——应用示例第三十四页，共五十四页，编辑于2023年，星期二例2：用镍硅藻土作催化剂，苯加氢合成环己烷。用微分反应器测定和分析得到160oC的初始反应速率以及相应的氢和苯的分压值pH/kPa，pB/Kpa和r0/(mol/Kg·h)的数据初始反应速率方程为：请利用上述实验数据拟合出参数ka、bH及bB。74.9467.2157.7650.7592.5626.7034.2443.4250.4310.20218.2220.8223.5189.2117.692.6687.6675.6456.1752.419.9714.7126.0745.0148.77115.1147.2217.8212.2202.4pHpH

pB

pBr0r02-2-3-4多元线性回归——应用示例第三十五页，共五十四页，编辑于2023年，星期二实验测得pH,pB和r0数据取倒数移项，开1/4方（线性化）2-2-3-4多元线性回归——应用示例第三十六页，共五十四页，编辑于2023年，星期二开始输入：自变量个数M，数据点个数N,反应温度T，气体分压pH，pB，R0的实验数据X(I,j)输出：bH、bB和k结束调用多元线性回归子程序计算B0，B1和B2M1=M+1，X1(I,j)＝X(I,j)X1(I,M1)＝DSQRT（DSQRT（X(I,1)**3*X(I,2)/X(I,M1)））（I=1,N）

X(I,J)=X1(I,J)（I=1,N，J=1,M1）计算：bH=A(1,M1)/B0；bB=A(2,M1)/B0；k=1/((B0**4)*(KH**3)*KB)显示程序

显示输入

显示输出2-2-3-4多元线性回归——应用示例第三十七页，共五十四页，编辑于2023年，星期二例3：化学反应动力学方程的总级数n设：某气体反应A+BC化学反应速率方程式中：ka—反应表观速率常数；

pA，pB

和pC—参加反应各气体A，B和C的气相分压；

，和—参加反应各物质在化学反应速率方程中的分级数。此反应的总级数n=++取对数n=++（线性化）2-2-3-4多元线性回归——应用示例第三十八页，共五十四页，编辑于2023年，星期二开始输入：数据点数N,反应温度T，气体分压p的实验数据,pa(I),pb(I),pc(I)（I=1，N）输出：n,ka，，和,结束调用多元线性回归子程序计算ka，，和X1(I)=ln(pa(I)),x2(I)=ln(pb(I)),x3(I)=ln(pb(I))Y(I)=lnv(I)(I=1,N)计算反应的总级数n=++2-2-3-4多元线性回归——应用示例第三十九页，共五十四页，编辑于2023年，星期二例4分子结构-性能的多元线性回归Hammet方程：苯环上间位或对位取代基对反应速率的影响式中：k——反应速率常数

k0——未取代时母体的反应速率常数ρ——反映取代基电子效应的结构常数σ——与反应类型有关的结构常数2-2-3-4多元线性回归——应用示例第四十页，共五十四页，编辑于2023年，星期二推广：考虑其他多种因素对分子性能的影响分子具有的性质：“应答”具有特定功能的结构参数π——反映取代基亲油/亲水能力E——反映取代基的空间效应肾上腺阻断剂：N-N-二甲基-2-溴苯乙胺衍生物结构—性能之间的关系数学模型：2-2-3-4多元线性回归——应用示例第四十一页，共五十四页，编辑于2023年，星期二

2-2-4多项式拟合简介——方法原理数学模型：残差的平方和Q残差的平方和最小，相当于求解以下方程组最小二乘法第四十二页，共五十四页，编辑于2023年，星期二2-3-1数值微分——问题的提出实验测定的一批离散点t/[t]t1t2……tncA/[cA]cA1cA2……cAn有化学反应

aA产物化学反应速率方程：第四十三页，共五十四页，编辑于2023年，星期二在恒容反应中，化学反应速率r可用反应物A浓度随时间变化率来表示式中，dcA/dt即为反应物A的浓度cA~时间t曲线上某点的斜率，数学上若cA~t是连续单值函数，则dcA/dt是在t点的一阶导数。t/[t]cA/[cA]t2cA22-3-1数值微分——问题的提出cA1t1第四十四页，共五十四页，编辑于2023年，星期二2-3-2数值微分——方法原理在微积分中，函数的导数是通过极限来定义的，有完整的计算方法。而在实际工作中，常常需要求列表函数在节点和非节点处的导数值，这就是数值微分要解决的问题。化学化工中有不少实际问题都需用数值微分求导来解决。

t/[t]t1t2……tncA/[cA]cA1cA2……cAn……第四十五页，共五十四页，编辑于2023年，星期二2-3-2数值微分——方法原理设有若干等距离的节点，欲求节点处的一阶导数，可用差商代替微商进行计算考虑函数f(x)在点（x0+h）处的泰勒展开式（向前差商式）CB（向后差商式）AC（中心差商式）ABxyx0-hx0x0+hABCy＝f(x)第四十六页，共五十四页，编辑于2023年，星期二中心差商定义：（1）若有函数y=f(t)，按等间隔在s个t点(t1,t2,…,ts-1,ts)测出相应的y值（y1

，y2

，…，ys

），除两端t1与ts以外，任意点的导数yi’=f’(ti)可按下式计算：（2）（i=2，3，…，s-1）若y＝f(x)为连续有界单值函数，客观存在在x点的一阶导数,可用较为精确的“中心差商”来计算：用插值法求2-3-2数值微分——方法原理第四十七页，共五十四页，编辑于2023年，星期二两端点的导数值用以下两式计算：（3）（4）式中t是独立变量t的间隔。：用插值法求2-3-2数值微分——方法原理第四十八页，共五十四页，编辑于2023年，星期二CF(N,X,Y,H,R)T=X(1)+H调用插值子程序计算FYA=FT=X(1)+2H调用插值子程序计算FYB=FW=2T=X(W)-HR(1)=[4YA-3Y(1)-YB]/2HABAB调用插值子程序计算FYA=FT=X(W)+h调用插值子程序计算FYB=FR(W)=[YA-YB]/2HW=W+1W>N-1noyesT=X(N)-Hcc调用插值子程序计算FYA=FT=X(N)-2H调用插值子程序计算FYB=FR(N)=[3Y(N)-4YA+YB]/2HRETURN2-3-3数值微分——程序框图第四十九页，共五十四页，编辑于2023年，星期二2-3-4数值微分——应用示例某抗菌素在人体血液中分解呈现简单级数反应，如果给病人在上午8点注射一针抗菌素（A），然后在不同时刻t测定抗菌素（A）在血液中的浓度cA（以mg·L-1表示），得到如下数据：

t/h24681012141618cA/mg·L-1

5.764.803.973.262.702.221.821.511.28求抗菌素（A）的消耗速率方程的反应级数和速率常数。（设抗菌素（A）的消耗速率方程为