提高化學反應動力學剛性求解器的準確度

葛銀海 王雁文 高聚

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

【關鍵詞】動力學；數值分析；化學反應；反應動力學；剛性常微分求解器

引言

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

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

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

下面列舉一個有符號解的二級反應系統

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

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

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

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

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

Robertson問題描述如下

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

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

3、提高數值計算的準確度

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

3.1 雅可比矩陣提高準確度

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

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

3.2 絕對誤差和相對誤差

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

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

4、結論

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

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

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