實驗數(shù)據(jù)的處理及模型參數(shù)的確定

實驗數(shù)據(jù)的處理及模型參數(shù)的確定第一頁，共五十一頁，2022年，8月28日●數(shù)學(xué)模型中各參數(shù)的確定例：鎳硅藻土上苯加氫合成環(huán)己烷是表面反應(yīng)控制的固體催化劑上的氣相反應(yīng)。在160oC，微分反應(yīng)器中的初始反應(yīng)速率方程為模型參數(shù)

ka

──表觀速率常數(shù)

bH

──H2的吸附系數(shù)

bB──

C6H6的吸附系數(shù)利用實驗得到的全部信息，確定數(shù)學(xué)模型中的待定參數(shù)第二頁，共五十一頁，2022年，8月28日

●線性插值

●

Lagrange插值

●埃米爾特插值

●一元線性回歸

●線性模型的推廣

●

多元回歸

可化為多元線性回歸的問題

●

多項式擬合簡介

●逐次回歸分析函數(shù)關(guān)系插值法回歸分析相關(guān)關(guān)系數(shù)值微分★★★★★引言：2.常用的數(shù)學(xué)方法第三頁，共五十一頁，2022年，8月28日例：72型分光光度計測得某試樣的吸收值如下：2-1-1–1線性插值——問題的提出λ/nm430440450460470480A0.4100.3750.3250.2800.2400.205希望：根據(jù)給定的函數(shù)表作一個既能反應(yīng)f(x)的特性，又便于計算的簡單函數(shù)p(x)，用p(x)近似f(x)，計算出任意

x對應(yīng)的y值求在435，445，455，465，475nm處的吸收值。第四頁，共五十一頁，2022年，8月28日

定義：設(shè)y=f(x)在區(qū)間[a,b]上有意義，且已知在點a<x0<x1<…<xn<b上的值y0,y1,…,yn,若存在一簡單函數(shù)pn(x),使

pn(xi)=yi(i=0,1,…,n)成立，則稱pn(x)為f(x)的插值函數(shù)，x0，x1，…，xn為插值節(jié)點區(qū)間[a,b]為插值區(qū)間，求pn(x)的方法稱為插值法

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

yn幾何意義：2-1-1–2線性插值——方法原理ab第五頁，共五十一頁，2022年，8月28日線性插值原理：兩點間直線方程：xyy=f(x)y=p(x)xi-1yi-1xi

yi2-1-1–2線性插值——方法原理第六頁，共五十一頁，2022年，8月28日分段線性插值：實驗點個數(shù)為n時，求插值結(jié)點x的函數(shù)值。首先確定x在哪兩點間

2-1-1–2線性插值——方法原理第七頁，共五十一頁，2022年，8月28日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線性插值——程序框圖第八頁，共五十一頁，2022年，8月28日開始輸入：數(shù)據(jù)點X(I),Y(I)，未知點X0調(diào)用線性插值子程序求未知點X0對應(yīng)的函數(shù)值Y0輸出：X0,Y0值結(jié)束2-1-1–4線性插值——應(yīng)用示例顯示程序顯示輸入顯示輸出第九頁，共五十一頁，2022年，8月28日2-1-2-1一元三點Lagrange插值——問題的提出例：計算乙醇的平均摩爾體積實驗測得25℃時乙醇溶液的平均摩爾體積（cm2mol-1）與乙醇的物質(zhì)的量分數(shù)的關(guān)系如下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時的。第十頁，共五十一頁，2022年，8月28日線性插值公式：二點(xi-1,yi-1),(xi,yi)(兩點式)Lagrange插值（三點插值,拋物線插值）：xi-1xixi+1即2-1-2-2一元三點Lagrange插值——方法原理第十一頁，共五十一頁，2022年，8月28日xyy=f(x)y=p(x)xi-1yi-1xi+1

yi+1xiyi編程難點：如何確定使用哪三個結(jié)點進行插值xj-1xj+1xjxj-2xj+22-1-2-2一元三點Lagrange插值——方法原理第十二頁，共五十一頁，2022年，8月28日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插值——程序框圖第十三頁，共五十一頁，2022年，8月28日開始輸入：數(shù)據(jù)點X(I),Y(I)，未知點X0調(diào)用lagrange插值子程序求未知點X0對應(yīng)的函數(shù)值Y0輸出：X0,Y0值結(jié)束2-1-2-4一元三點Lagrange插值——應(yīng)用示例顯示程序顯示輸入顯示輸出第十四頁，共五十一頁，2022年，8月28日2-2-1-1一元線性回歸——問題的提出例:銅鉬礦中鉬對銅含量的線性依賴關(guān)系No123456789x(Cu)285290300303310318325335338y(Mo)4.64.74.74.94.95.15.05.35.4一元線性回歸的數(shù)學(xué)模型：y=ax+b+εyi=axi+b+εin個實驗點回歸直線：y=ax+b殘差：εi=yi-（axi+b）第十五頁，共五十一頁，2022年，8月28日最小二乘法：第i點殘差：εi=yi-（axi+b）當(dāng)殘差的平方和為最小時，對應(yīng)的a、b值是最佳值。（正規(guī)方程組）2-2-1-2一元線性回歸——方法原理第十六頁，共五十一頁，2022年，8月28日令平均值離差平方和2-2-1-2一元線性回歸——方法原理第十七頁，共五十一頁，2022年，8月28日1.線性相關(guān)系數(shù)R—衡量回歸方程式與數(shù)據(jù)相符合的程度。若R1，則數(shù)據(jù)點落在直線上。注意：2.加權(quán)最小二乘法3.剔除可疑數(shù)據(jù)2-2-1-2一元線性回歸——方法原理第十八頁，共五十一頁，2022年，8月28日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一元線性回歸——程序框圖第十九頁，共五十一頁，2022年，8月28日開始輸入：數(shù)據(jù)點數(shù)N

銅與鉬的實驗數(shù)據(jù)X(I),Y(I)（I=1，N）調(diào)用一元線性回歸子程序計算A，B，R輸出：A，B，R結(jié)束2-2-1-4一元線性回歸——應(yīng)用示例顯示程序顯示輸入顯示輸出第二十頁，共五十一頁，2022年，8月28日2-2-2-1線性模型的推廣——方法原理變量x與y之間存在某種非線性關(guān)系確定曲線類型（非線性關(guān)系）實際經(jīng)驗散點圖形狀線性關(guān)系最小二乘法確定系數(shù)非線性關(guān)系第二十一頁，共五十一頁，2022年，8月28日曲線類型及變換公式雙曲線型冪指數(shù)型指數(shù)型S型對數(shù)型平方根曲線2-2-2-1線性模型的推廣——方法原理第二十二頁，共五十一頁，2022年，8月28日例1：Arrhenius公式的應(yīng)用根據(jù)k和T數(shù)據(jù)，可確定指前因子A和活化能Ea。2-2-2-2線性模型的推廣——應(yīng)用示例第二十三頁，共五十一頁，2022年，8月28日例2：Clausius-Clapryron方程式的應(yīng)用純組分氣-液（氣-固）兩相平衡的方程式：上式中：p:T/K時液（固）飽和蒸氣壓；ΔH：相變熱不定積分：測定不同溫度下的飽和蒸氣壓，將lnp~1/T進行線性回歸，可算出ΔH，并計算其它溫度下的蒸氣壓2-2-2-2線性模型的推廣——應(yīng)用示例第二十四頁，共五十一頁，2022年，8月28日開始輸入：數(shù)據(jù)點數(shù)N

溫度T與蒸氣壓p的實驗數(shù)據(jù)T(I),P(I)（I=1，N）輸出：A,B,ΔH,Ti,pi結(jié)束調(diào)用線性回歸子程序計算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線性模型的推廣——應(yīng)用示例第二十五頁，共五十一頁，2022年，8月28日例3：幾種常用的吸附等溫式的計算（氣固吸附）2-2-2-2線性模型的推廣——應(yīng)用示例Freundlich經(jīng)驗式：Langmuir方程：B.E.T.方程：第二十六頁，共五十一頁，2022年，8月28日2-2-3-1多元線性回歸——方法原理數(shù)學(xué)模型：函數(shù)y與多個自變量間x1，x2，…，xm的線性相關(guān)關(guān)系設(shè)共進行了n次測定（i=1，2，…，n）的自變量取值分別為xi1，xi2，…，xim，函數(shù)值的測定值為yi，所得的值如表：第二十七頁，共五十一頁，2022年，8月28日選擇回歸系數(shù)，以使殘差的平方和最小，即殘差殘差的平方和Q殘差的平方和最小，相當(dāng)于求解以下方程組（i=1，2，…，n）（k=，0，1，…，m）最小二乘法2-2-3-1多元線性回歸——方法原理第二十八頁，共五十一頁，2022年，8月28日正規(guī)方程組：全相關(guān)系數(shù)R式中（m+1）元線性方程組，Guass消去法求解b0，b1，…，bm2-2-3-1多元線性回歸——方法原理第二十九頁，共五十一頁，2022年，8月28日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多元線性回歸第三十頁，共五十一頁，2022年，8月28日變量xi(i=1,2,…,m)與y之間存在某種非線性關(guān)系確定曲線類型（非線性關(guān)系）實際經(jīng)驗散點圖形狀線性關(guān)系最小二乘法確定系數(shù)非線性關(guān)系2-2-3-3多元線性回歸——可化為多元線性回歸的問題變量代換變量代換第三十一頁，共五十一頁，2022年，8月28日函數(shù)式通式：y’=b0+b1x1+b2x2變換變量變換常數(shù)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多元線性回歸——可化為多元線性回歸的問題第三十二頁，共五十一頁，2022年，8月28日例1：Antoine方程式的應(yīng)用（p：蒸氣壓，T：溫度）令Tlgp=y,T=x1,lgp=x2,b0=ac+b,b1=a,b2=c2-2-3-4多元線性回歸——應(yīng)用示例第三十三頁，共五十一頁，2022年，8月28日例2：用鎳硅藻土作催化劑，苯加氫合成環(huán)己烷。用微分反應(yīng)器測定和分析得到160oC的初始反應(yīng)速率以及相應(yīng)的氫和苯的分壓值pH/kPa，pB/Kpa和r0/(mol/Kg·h)的數(shù)據(jù)初始反應(yīng)速率方程為：請利用上述實驗數(shù)據(jù)擬合出參數(shù)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多元線性回歸——應(yīng)用示例第三十四頁，共五十一頁，2022年，8月28日實驗測得pH,pB和r0數(shù)據(jù)取倒數(shù)移項，開1/4方（線性化）2-2-3-4多元線性回歸——應(yīng)用示例第三十五頁，共五十一頁，2022年，8月28日開始輸入：自變量個數(shù)M，數(shù)據(jù)點個數(shù)N,反應(yīng)溫度T，氣體分壓pH，pB，R0的實驗數(shù)據(jù)X(I,j)輸出：bH、bB和k結(jié)束調(diào)用多元線性回歸子程序計算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多元線性回歸——應(yīng)用示例第三十六頁，共五十一頁，2022年，8月28日例3：化學(xué)反應(yīng)動力學(xué)方程的總級數(shù)n設(shè)：某氣體反應(yīng)A+BC化學(xué)反應(yīng)速率方程式中：ka—反應(yīng)表觀速率常數(shù)；

pA，pB

和pC—參加反應(yīng)各氣體A，B和C的氣相分壓；

，和—參加反應(yīng)各物質(zhì)在化學(xué)反應(yīng)速率方程中的分級數(shù)。此反應(yīng)的總級數(shù)n=++取對數(shù)n=++（線性化）2-2-3-4多元線性回歸——應(yīng)用示例第三十七頁，共五十一頁，2022年，8月28日開始輸入：數(shù)據(jù)點數(shù)N,反應(yīng)溫度T，氣體分壓p的實驗數(shù)據(jù),pa(I),pb(I),pc(I)（I=1，N）輸出：n,ka，，和,結(jié)束調(diào)用多元線性回歸子程序計算ka，，和X1(I)=ln(pa(I)),x2(I)=ln(pb(I)),x3(I)=ln(pb(I))Y(I)=lnv(I)(I=1,N)計算反應(yīng)的總級數(shù)n=++2-2-3-4多元線性回歸——應(yīng)用示例第三十八頁，共五十一頁，2022年，8月28日例4分子結(jié)構(gòu)-性能的多元線性回歸Hammet方程：苯環(huán)上間位或?qū)ξ蝗〈鶎Ψ磻?yīng)速率的影響式中：k——反應(yīng)速率常數(shù)

k0——未取代時母體的反應(yīng)速率常數(shù)ρ——反映取代基電子效應(yīng)的結(jié)構(gòu)常數(shù)σ——與反應(yīng)類型有關(guān)的結(jié)構(gòu)常數(shù)2-2-3-4多元線性回歸——應(yīng)用示例第三十九頁，共五十一頁，2022年，8月28日推廣：考慮其他多種因素對分子性能的影響分子具有的性質(zhì)：“應(yīng)答”具有特定功能的結(jié)構(gòu)參數(shù)π——反映取代基親油/親水能力E——反映取代基的空間效應(yīng)腎上腺阻斷劑：N-N-二甲基-2-溴苯乙胺衍生物結(jié)構(gòu)—性能之間的關(guān)系數(shù)學(xué)模型：2-2-3-4多元線性回歸——應(yīng)用示例第四十頁，共五十一頁，2022年，8月28日

2-2-4多項式擬合簡介——方法原理數(shù)學(xué)模型：殘差的平方和Q殘差的平方和最小，相當(dāng)于求解以下方程組最小二乘法第四十一頁，共五十一頁，2022年，8月28日2-3-1數(shù)值微分——問題的提出實驗測定的一批離散點t/[t]t1t2……tncA/[cA]cA1cA2……cAn有化學(xué)反應(yīng)

aA產(chǎn)物化學(xué)反應(yīng)速率方程：第四十二頁，共五十一頁，2022年，8月28日在恒容反應(yīng)中，化學(xué)反應(yīng)速率r可用反應(yīng)物A濃度隨時間變化率來表示式中，dcA/dt即為反應(yīng)物A的濃度cA~時間t曲線上某點的斜率，數(shù)學(xué)上若cA~t是連續(xù)單值函數(shù)，則dcA/dt是在t點的一階導(dǎo)數(shù)。t/[t]cA/[cA]t2cA22-3-1數(shù)值微分——問題的提出cA1t1第四十三頁，共五十一頁，2022年，8月28日2-3-2數(shù)值微分——方法原理在微積分中，函數(shù)的導(dǎo)數(shù)是通過極限來定義的，有完整的計算方法。而在實際工作中，常常需要求列表函數(shù)在節(jié)點和非節(jié)點處的導(dǎo)數(shù)值，這就是數(shù)值微分要解決的問題?；瘜W(xué)化工中有不少實際問題都需用數(shù)值微分求導(dǎo)來解決。

t/[t]t1t2……tncA/[cA]cA1cA2……cAn……第四十四頁，共五十一頁，2022年，8月28日2-3-2數(shù)值微分——方法原理設(shè)有若干等距離的節(jié)點，欲求節(jié)點處的一階導(dǎo)數(shù)，可用差商代替微商進行計算考慮函數(shù)f(x)在點（x0+h）處的泰勒展開式（向前差商式）CB（向后差商式）AC（中心差商式）ABxyx0-hx0x0+hABCy＝f(x)第四十五頁，共五十一頁，2022年，8月28日中心差商定義：（1）若有函數(shù)y=f(t)，按等間隔在s個t點(t1,t2,…,ts-1,ts)測出相應(yīng)的y值（y1

，y2

，…，ys

），除兩端t1與ts以外，任意點的導(dǎo)數(shù)yi’=f’(ti)可按下式計算：（2）（i=2，3，…，s-1）若y＝f(x)為連續(xù)有界單值函數(shù)，客觀存在在x點的一階導(dǎo)數(shù),可用較為精確的“中心差商”來計算：用插值法求2-3-2數(shù)值微分——方法原理第四十六頁，共五十一頁，2022年，8月28日兩端點的導(dǎo)數(shù)值用以下兩式計算：（3）（4）式中t是獨立變量t的間隔。：用插值法求2-3-2數(shù)值微分——方法原理第四十七頁，共五十一頁，2022年，8月28日CF(N,X,Y,H,R)T=X(1)+H調(diào)用插值子程序計算FYA=FT=X(1)+2H調(diào)用插值子程序計算FYB=FW=2T=X(W)-HR(1)=[4YA-3Y(1)-YB]/2HAB