提高化學(xué)反應(yīng)動力學(xué)剛性求解器的準(zhǔn)確度

葛銀海 王雁文 高聚

【摘要】隨著科技的發(fā)展，研究人員對數(shù)值計算的精度和可靠性提出更高的要求。對于特定性問題，特別是在燃燒化學(xué)領(lǐng)域方面，通常需要自己編寫程序計算，然而現(xiàn)在數(shù)值計算結(jié)果分析存在著一些困難：在燃燒化學(xué)反應(yīng)動力學(xué)方程組，一般是非線性常微分方程組，由于沒有符號解，所以很難判斷求解的準(zhǔn)確性；其次通常用不同的求解器求解常微分方程時，會出現(xiàn)不同的計算結(jié)果，就更難以判斷求解的準(zhǔn)確性。列舉幾種常用的求解常微分方程剛性求解器，進行數(shù)值分析，比較它們在求解剛性方程中的優(yōu)劣，以及在計算化學(xué)反應(yīng)動力學(xué)方程組中，如何提高計算的準(zhǔn)確性。

【關(guān)鍵詞】動力學(xué)；數(shù)值分析；化學(xué)反應(yīng)；反應(yīng)動力學(xué)；剛性常微分求解器

引言

許多物理和化學(xué)定律都是以微分形式存在的。在生物學(xué)和經(jīng)濟學(xué)中，用微分方程來建立復(fù)雜的數(shù)學(xué)模型，微分方程在各個領(lǐng)域中得到廣泛的應(yīng)用。隨著計算理論的發(fā)展，許多復(fù)雜的科學(xué)和工程問題得以分析解決，并用來估計未知的參數(shù)和狀態(tài)。在化學(xué)反應(yīng)動力學(xué)方面，對于化學(xué)工作者來說，如何應(yīng)用和利用這些求解器，懂得求解器的優(yōu)劣，并利用這些求解器求得更加準(zhǔn)確的解是非常重要的。本文將介紹我們常用的求解常微分組剛性求解器如MATLAB中的ODE15S，ODE23S，以及求解化學(xué)動力學(xué)方程組中常用的剛性求解器VODE（a variable-coefficient ode solver）。

1、有符號解的常微分方程組

對于常微分方程組是剛性還是非剛性，并沒有明確的界限。我們所定義的剛性方程是指，其數(shù)值解只有在時間間隔很小時才會穩(wěn)定，而當(dāng)時間間隔大時，方程的解就容易不穩(wěn)定，也就是在計算過程中出現(xiàn)方程的解快速變化的程度高。

下面列舉一個有符號解的二級反應(yīng)系統(tǒng)

其中Cx、Cy、Cz表示x、y、z的摩爾溶度；初始條件為Cx=1，Cy=0，Cz=0。在比較這些求解器時，都設(shè)置絕對誤差和相對誤差為1E-6，其他的默認。

當(dāng)求解這樣的化學(xué)反應(yīng)動力學(xué)方程組時，首先要知道它有沒有符號解，因為有符號解，就能得到方程組的精確解。為了求解方便和便于比較，假設(shè)k1=1，k2=1，然后通過MAPLE軟件編寫程序就可以得到方程得精確解，在求解方程組的符號解時，MAPLE是非常優(yōu)秀的軟件。再得到常微分方程組精確解，就能比較幾個剛性求解器的求解準(zhǔn)確度，如下所示。實線表示精確解，符號表示常微分求解器的結(jié)果。

從求解這幾個常用的求解器的求解結(jié)果可知，一般情況下，這幾種流行的求解器的求解的準(zhǔn)確度是很高的，而且容易進行誤差比較。

2、無符號解的常微分方程組

在實際應(yīng)用中，我們所遇到的化學(xué)動力學(xué)問題一般都是沒有符號解的，并且是剛性很強的系統(tǒng)或者說常微分方程的解出現(xiàn)快速變化的程度很高的情況，比如著名的化學(xué)動力學(xué)剛性模型Robertson問題。如果用非剛性常微分求解器如ODE45去求解Robertson問題是得不到結(jié)果的，但許多剛性求解卻能夠求解Robertson問題，而這些剛性求解器也以求解Robertson問題來評價自身的求解器，包括ODE15S，ODE23S，VODE，LSODE等。

Robertson問題描述如下

盡管沒有符號解，不知道方程組的準(zhǔn)確解，但仍可以從其他方面來判斷方程組解的準(zhǔn)確性，在計算化學(xué)反應(yīng)動力學(xué)方程組中有隱藏的等式，比如摩爾溶度守恒等于1、能量守恒等。因為這個方程組不涉及能量方程，用CX+CY+CZ=1去判斷求解的非線性剛性常微分方程組求解器誤差，為了便于分析和比較，我們把時間t的范圍設(shè)置為1E-6到1E+6，用常微分求解器求解和調(diào)用求解器到計算完成所需時間的結(jié)果如下：

從上面的求解結(jié)果可知：這三種求解器求解Robertson問題的值基本一樣，沒有太大偏差，MATLAB的ODE15S和ODE23S求解器在求解的摩爾溶度守恒誤差更小，ODE15S求解器因為采用變階次多步解法在精度上優(yōu)于ODE23S，但計算速度上要多于ODE15S；VODE在計算所用的時間領(lǐng)先于ODE15S和ODE23S，計算這樣的方程組基本不花費時間。用上面的3個求解器求解的結(jié)果與文獻12比較可知，這三個求解器求解溶度守恒誤差的結(jié)果小于C語言編寫的deSolve求解器，而這四個求解器都是雙精度求解器，但求解的準(zhǔn)確度卻不一樣。結(jié)合這兩個例子，可知如果方程組少、計算量小的情況下，可用MATLAB剛性求解器求解，但如果計算量大，建議用FORTRAN或C語言求解器來計算，如VODE、deSolve求解器。

3、提高數(shù)值計算的準(zhǔn)確度

在化學(xué)動力方程組求解中，提高求解的準(zhǔn)確度是非常重要的，因為數(shù)值計算本身與實際運行的系統(tǒng)就存在誤差，如果本身求解就存在很大誤差的話，那么求解的結(jié)果就不可靠了，不能指導(dǎo)實際工作，所以下面將討論如何提高計算的準(zhǔn)確度。

3.1 雅可比矩陣提高準(zhǔn)確度

在求解剛性非線性常微分方程組時，這些剛性求解器都鼓勵給方程組提供雅可比矩陣，因為雅可比矩陣的重要特性是它體現(xiàn)了一個可微方程與給出點的最優(yōu)線性逼近，有助于求得準(zhǔn)確的方程組解。下面以VODE求解器為例子，分析方程組的優(yōu)解。在VODE求解器中，相對誤差和絕對誤差都設(shè)置為1E-6，X1、Y1、Z1表示未加入雅可比矩陣所求得的組分值，X2、Y2、Z2表示未加入雅可比矩陣所求得的組分值，加入雅可比矩陣與未加入雅可比矩陣計算結(jié)果如下：

由圖5可知，在未加入雅可比矩陣時，隨著時間的跨度變大，求解的誤差也隨之變大，特別是組分值比較小的，時間在4.0E+05后，誤差就變得很大，求解剛性問題就很不準(zhǔn)確了，但加入雅可比矩陣后，非常有利于求解器求得準(zhǔn)確值，準(zhǔn)確度更高；從表格中可知，相比于Z組分，當(dāng)小組分X，Y越小時，不加入雅可比矩陣，其誤差越大。

3.2 絕對誤差和相對誤差

影響常微分方程組求解準(zhǔn)確的因素很多，就像上面討論的那樣，如果給求解器提供雅可比矩陣，那么求解器的精度就會有很大程度的提高；而另一方面就是設(shè)置好求解的相對誤差和絕對誤差。準(zhǔn)確度有別于絕對誤差和相對誤差，不是所設(shè)置的相對誤差和絕對誤差高，計算結(jié)果的準(zhǔn)確度就高。絕對誤差是指，真值P與其所求值P*的絕對值差，公式表示為ATOL=|P-P*|；相對誤差是指，求解的值所造成的絕對誤差與真值之比RTOL=|P-P*|/P。下面我們舉個例子，一個是設(shè)置相對誤差和絕對誤差都為1.0E-12，所求組分摩爾溶度值表示為X1、Y1、Z1；一個是設(shè)置相對誤差為1E-4，X、Y、Z組分絕對值誤差分別設(shè)置為，ATOL（X）= 1.0E-8，ATOL（Y）= 1.0E-12，ATOL（Z）= 1.0E-6，所求組分摩爾溶度值表示為X2、Y2、Z2，求解結(jié)果如下。

在求解常微分方程組時，相對誤差和絕對誤差的設(shè)置對結(jié)果的影響是很大，但很多時候無從抉擇，可以從摩爾溶度誤差和能量守恒誤差中提供參考。從圖6和表格三可知，在設(shè)置相對誤差和絕對誤差時，并不是設(shè)置精度越高求解的準(zhǔn)確度越高，誤差越小，這跟所求得組分溶度有關(guān)，也跟你所需要保留的有效位數(shù)有關(guān)。

4、結(jié)論

a）在用不同的常微分求解器進行數(shù)值準(zhǔn)確度時，對于化學(xué)反應(yīng)動力學(xué)常微分方程組來說，如果有符號解，就可以給出每一個方程組的誤差，來判斷求解結(jié)果是否可信；如果沒有符號解，可以通過摩爾溶度和能量守恒中給出誤差值，供參考。

b）在求解常微分方程組時，MATLAB中的剛性求解器相比于VODE和deSolve更加準(zhǔn)確，特別是基于變階次多步解法的ODE15S，但計算時間上多于C語言編寫的deSolve和FOREAN語言編寫的VODE。

c）在提高求解器的準(zhǔn)確度方面，最好要向求解器提供雅可比矩陣，使計算誤差降至最??；求解器求解的準(zhǔn)確度，受所設(shè)置的相對誤差和絕對誤差的影響很大，并不是相對誤差和絕對誤差設(shè)置的精度越高越好，要設(shè)置好相對誤差和絕對誤差，通常需要多次運算才能取得。