左 路， 胡 瑋， 鐘欣欣

(湖北大學(xué)化學(xué)化工學(xué)院，武漢 430062)

1 引 言

化學(xué)動力學(xué)的重要目的就是探索反應(yīng)過程中不同反應(yīng)物質(zhì)量的關(guān)聯(lián)，描述其濃度變化的演變規(guī)律，并預(yù)測反應(yīng)系統(tǒng)宏觀行為.反應(yīng)物質(zhì)量的關(guān)聯(lián)及演變規(guī)律由反應(yīng)速率微分方程描述[1-2]，微分方程組的建立則依賴于化學(xué)動力學(xué)參數(shù)——反應(yīng)速率常數(shù).根據(jù)是否已知速率常數(shù)，化學(xué)反應(yīng)動力學(xué)研究通常呈現(xiàn)兩種互逆的研究路徑：已知速率常數(shù)時，考慮給定的反應(yīng)初始態(tài)，由解軌跡描述并預(yù)測反應(yīng)系統(tǒng)的宏觀行為；未知速率常數(shù)時，研究側(cè)重于如何由實(shí)驗(yàn)路徑估計(jì)動力學(xué)信息.實(shí)際中研究者常面對后一種情形，因此估計(jì)反應(yīng)速率常數(shù)成為反應(yīng)化學(xué)動力學(xué)研究的首要任務(wù).但是估計(jì)方法需要面對若干困難，諸如如何進(jìn)行有效的實(shí)驗(yàn)并獲得精確的實(shí)驗(yàn)路徑，以及當(dāng)反應(yīng)系統(tǒng)復(fù)雜度增加時，如何利用有限的實(shí)驗(yàn)路徑估計(jì)所有的反應(yīng)速率常數(shù).前述情形在催化與酶化學(xué)中很常見，此時的反應(yīng)體系通常為高耦合復(fù)雜系統(tǒng)，因此造成高維度背景下化學(xué)動力學(xué)研究異常困難.

1962年，文獻(xiàn)[3]建立了復(fù)雜體系結(jié)構(gòu)性分析的基本方法，利用速率方程組系數(shù)矩陣的代數(shù)特性構(gòu)建估計(jì)速率常數(shù)的計(jì)算方法，從而開辟了利用矩陣方法研究化學(xué)動力系統(tǒng)的新途徑.基于Wei和Prater結(jié)構(gòu)方法，文獻(xiàn)[2]建立了反應(yīng)速率常數(shù)的正交估計(jì)法.本文將深入分析一級可逆反應(yīng)體系速率方程組，不僅證明了速率常數(shù)矩陣特征系統(tǒng)的固有性質(zhì)，也揭示出這些性質(zhì)如何構(gòu)建得到反應(yīng)體系的結(jié)構(gòu)系統(tǒng).

2 速率常數(shù)矩陣特征體系性質(zhì)

dx(t)/dt=Kx(t)，

其中系數(shù)矩陣K由速率常數(shù)構(gòu)成

dx(t)/dt=Kx(t)中每個微分方程形如

上式右端函數(shù)受控于t時刻所有物質(zhì)的量，即時刻t處物質(zhì)Ai的反應(yīng)速率dxi(t)/dt受到所有與Ai相關(guān)的反應(yīng)通道影響.xi項(xiàng)的系數(shù)則對應(yīng)于所有由物質(zhì)Ai反應(yīng)的通道.由于反應(yīng)物質(zhì)的量必為非負(fù)值，時刻t處各物質(zhì)的摩爾分率必須滿足

(1)

深入分析速率常數(shù)矩陣K的特征體系，不僅揭示出K的結(jié)構(gòu)性質(zhì)有助于建立文獻(xiàn)[2]中的正交估計(jì)法，同時揭示了一級可逆反應(yīng)體系的結(jié)構(gòu)化特征.

2.1 速率常數(shù)矩陣K有唯一零特征值

觀察速率常數(shù)矩陣可知K的每列元素和為零，則K為奇異方陣，必有零特征值.但是僅由K的行列式值無法直接判斷零特征值的代數(shù)重?cái)?shù).記K的特征多項(xiàng)式為

p(λ)=λn+c1λn-1+c2λn-2+…+cn-2λ2+cn-1λ+cn，

其中系數(shù)ck=(-1)kSk,Sk為K的所有k階主子式的和.由于K為奇異方陣，則系數(shù)cn=0，特征多項(xiàng)式可以變化為

p(λ)=λ(λn-1+c1λn-2+c2λn-3+…+cn-2λ+cn-1).

若系數(shù)cn-1≠0，則特征方程p(λ)=0僅有唯一零根.

考慮K的取自i1,i2,…,in-1行、列的n-1階主子式

對于一級可逆反應(yīng)網(wǎng)絡(luò)而言，速率常數(shù)矩陣K的零元素必在對稱位置出現(xiàn)，即若某個kij=0，則必有逆向反應(yīng)的速率常數(shù)kji=0.并且K的每一行、列元素必不全為零，否則零行、列對應(yīng)的物質(zhì)將不在反應(yīng)系統(tǒng)中存在.因此K的任一n-1階主子式對應(yīng)的矩陣K(n-1)至少有一列滿足嚴(yán)格對角占優(yōu)且不可約，則K的所有n-1階主子式值非零[4].但是僅由此點(diǎn)性質(zhì)無法判斷所有n-1階主子式值的和為非零值.

令k=n-1，主子式|K(k)|的展開式中必有一項(xiàng)為主對角線上k個元素的乘積，即為

(3)

化簡得到

此項(xiàng)展開式中至少有一個乘積項(xiàng)kj1i1kj2i2…kjn-1in-1不為零，否則|K(k)|的所有元素必為0.并且當(dāng)k為奇數(shù)時，(3)式值小于零，否則(3)式值將大于零.

|K(k)|展開式中的其他項(xiàng)最多含有k-2個主對角元素，考慮展開式中包含k-2個主對角元素的項(xiàng)，則此項(xiàng)乘積因子中某兩個位置的元素取自非對角元素，形如

(4)

可見(4)式與(3)式符號一致.考慮展開式中包含k-3個主對角元素的項(xiàng)，形如

(5)

此項(xiàng)與(4)式符號相反.形如(5)式的每一項(xiàng)均包含在(4)式形式中，僅僅符號相反.因此符合(5)式結(jié)構(gòu)的項(xiàng)均在形如(4)式的項(xiàng)中抵消掉.依此類推，從k-3開始每一個包含k-(2j+1)個主對角元的項(xiàng)必包含于含k-2j個主對角元的展開項(xiàng)中，且由于符號相反得以抵消.|K(k)|最終值為k個主對角元的乘積項(xiàng)、k-2j個主對角元乘積項(xiàng)的抵消剩余項(xiàng)，以及不含主對角元的乘積項(xiàng)這三部分之和.對于k階主子式而言，前述每一項(xiàng)均符號相同，則K的每一個k階主子式非零.因此任意n-1階主子式|K(n-1)|非零且符號相同，使得所有n-1階主子式和Sn-1非零，從而系數(shù)cn-1≠0，因此特征方程p(λ)=0僅有1個零根.

2.2 零特征值對應(yīng)的特征向量構(gòu)成系統(tǒng)平衡狀態(tài)向量

由統(tǒng)計(jì)熱力學(xué)中的精細(xì)平衡原則[5-7](the principle of detailed balance)可知，可逆反應(yīng)系統(tǒng)經(jīng)歷一定反應(yīng)時間后將達(dá)到唯一平衡態(tài).實(shí)際上，速率常數(shù)矩陣的構(gòu)成就已經(jīng)充分保證了一級可逆反應(yīng)體系唯一平衡態(tài)的存在，且此平衡態(tài)構(gòu)成對應(yīng)于速率常數(shù)矩陣零特征值的特征向量.

dx(t)/dt=Kx(t)

在t→+∞時滿足dxeq/dt=Kxeq=0，即K的零特征值(記為λ1=0)對應(yīng)于平衡態(tài)向量.由于零特征值的代數(shù)重?cái)?shù)為1且?guī)缀沃財(cái)?shù)亦為1，則K的零特征值對應(yīng)的全部特征向量都在n中的一條直線上.但是分量值均為正且分量和為1的特征向量僅有一個(記為v1)，并滿足dv1/dt=Kv1=0.此時反應(yīng)體系達(dá)到平衡態(tài)xeq(xeq=v1)，即v1為滿足反應(yīng)系統(tǒng)約束的對應(yīng)于零特征值唯一特征向量.

2.3 K的所有特征值均小于或等于零

Γ的所有特征值全為實(shí)數(shù)且小于或等于零[2](記為-λj,j=1,2,…,n).由相似關(guān)系不變性，-λ1,-λ2,…,-λn同時也是K的特征值.

滿足Γ=PΛP-1，其中對角矩陣Λ=diag(-λ1,-λ2,…,-λn)，那么有

若記

則得到K=VΛV-1，因此v1,v2,…,vn為K對應(yīng)于-λ1,-λ2,…,-λn的特征向量.

2.4 非零特征值對應(yīng)的特征向量的分量和必為0

令特征向量vj=(v1j,v2j,…,vnj)T,j=1,2,…,n，由于Kvj=-λjvj，則非零特征值-λj,j=2,3,…,n對應(yīng)的特征向量的分量和應(yīng)滿足

3 速率常數(shù)矩陣特征體系與動力學(xué)參數(shù)估計(jì)

3.1 反應(yīng)軌跡結(jié)構(gòu)

利用K=VΛV-1解耦反應(yīng)體系dx(t)/dt=Kx(t)得到微分方程dy(t)/dt=Λy(t)，解軌跡x(t)則變換為虛擬解軌跡

y(t)=V-1x(t)=(y1(t),y2(t),…,yn(t))T.

解耦后的每個微分方程形如dyj(t)/dt=-λjyj(t)，此時的反應(yīng)組分yj(t)(以下稱為虛擬組分)并非真實(shí)反應(yīng)物摩爾組分?jǐn)?shù)，因此不必受到(1)式約束.求解微分方程dy(t)/dt=Λy(t)將得到虛擬解軌跡y(t)=eΛty(0)，其中矩陣指數(shù)函數(shù)eΛt=diag(e-λ1t,e-λ2t,…,e-λnt),y(0)=V-1x(0)為反應(yīng)初始態(tài)x(0)對應(yīng)的初始虛擬組分態(tài).將虛擬解軌跡y(t)還原得到原始反應(yīng)體系dx(t)/dt=Kx(t)的解軌跡

x(t)=VeΛtV-1x(0)=c1v1+c2e-λ2tv2+…+cne-λntvn，

(6)

其中系數(shù)c1,c2,…,cn即為虛擬初始狀態(tài)y(0)的所有分量值.當(dāng)反應(yīng)體系的速率常數(shù)kij未知時，利用(6)式可以估計(jì)矩陣K的特征值與特征向量.

由于v1,v2,…,vn線性無關(guān)可以構(gòu)成n的基，因此初始態(tài)向量表示為

x(0)=c1v1+c2v2+…+cnvn.

而狀態(tài)向量x(0)與平衡態(tài)向量v1必須滿足約束(1)式，且非零特征值的特征向量分量和為0，則上式右端向量的分量和應(yīng)為1，即

得到系數(shù)c1=1.因此(6)式化簡為

x(t)=v1+c2e-λ2tv2+…+cne-λntvn.

(7)

當(dāng)反應(yīng)時間t→+∞時，較大絕對值的非零特征值-λj使得cje-λjtvj快速趨于零向量，而具有最小絕對值的非零特征值對應(yīng)的特征向量將決定x(t)沿什么路徑趨于平衡態(tài).

3.2 特征直線反應(yīng)路徑

實(shí)驗(yàn)反應(yīng)路徑通常選擇從最簡單的初始組分開始，即從反應(yīng)空間的某個頂點(diǎn)對應(yīng)的純物質(zhì)初始態(tài)進(jìn)行反應(yīng)，并測得反應(yīng)軌跡與平衡態(tài)v1.反應(yīng)軌跡趨近平衡態(tài)時的切線就是|-λj|min(-λj≠0)對應(yīng)的特征向量vj所在直線，直線擬合平衡態(tài)附近的反應(yīng)軌跡，并求得擬合直線與反應(yīng)空間邊界的交點(diǎn)，即為沿特征向量vj所在直線的反應(yīng)進(jìn)程初始點(diǎn)(記為xvj(0)).此初始點(diǎn)對應(yīng)于(7)式的系數(shù)分別為1,0,…,0,cj,0,…,0，使得(7)式簡化為特征直線反應(yīng)路徑xvj(t)=v1+cje-λjtvj.由于xvj(0)=v1+cjvj，為了方便將重新放縮后的特征向量cjvj作為與特征值-λj對應(yīng)的特征向量，仍記為vj.由xvj(0)-v1從而求得特征向量vj.

K=Vdiag(-λ1,-λ2,…,-λn)V-1.

4 結(jié) 論

對于一個n組分的一級可逆復(fù)雜反應(yīng)網(wǎng)絡(luò)而言，速率常數(shù)矩陣至少包含2(n-1)個動力學(xué)參數(shù)，至多包含n(n-1)個.隨著組分?jǐn)?shù)n的增加，系統(tǒng)復(fù)雜度與待估計(jì)參數(shù)個數(shù)急速增加.如何利用有限的實(shí)驗(yàn)反應(yīng)路徑估計(jì)動力學(xué)參數(shù)并描述反應(yīng)體系的宏觀行為便成為動力學(xué)研究的困難之處.深入分析速率常數(shù)矩陣的結(jié)構(gòu)特征，不僅建立了基于矩陣特征體系的反應(yīng)模型及反應(yīng)軌跡結(jié)構(gòu)性描述，還提供了利用少量純物質(zhì)反應(yīng)初始態(tài)的實(shí)驗(yàn)路徑估計(jì)動力學(xué)參數(shù)的方法.與此同時，基于矩陣特征體系的反應(yīng)軌跡幾何視角描述，為進(jìn)一步優(yōu)化實(shí)驗(yàn)路徑、尋找最短實(shí)驗(yàn)路徑指明了研究方向.

致謝作者衷心感謝審稿專家提出的寶貴意見，以及相關(guān)文獻(xiàn)對本文的啟發(fā)與激勵.