提高化学反应动力学刚性求解器的准确度

葛银海 王雁文 高聚

【摘要】随着科技的发展，研究人员对数值计算的精度和可靠性提出更高的要求。对于特定性问题，特别是在燃烧化学领域方面，通常需要自己编写程序计算，然而现在数值计算结果分析存在着一些困难：在燃烧化学反應动力学方程组，一般是非线性常微分方程组，由于没有符号解，所以很难判断求解的准确性；其次通常用不同的求解器求解常微分方程时，会出现不同的计算结果，就更难以判断求解的准确性。列举几种常用的求解常微分方程刚性求解器，进行数值分析，比较它们在求解刚性方程中的优劣，以及在计算化学反应动力学方程组中，如何提高计算的准确性。

【关键词】动力学；数值分析；化学反应；反应动力学；刚性常微分求解器

引言

许多物理和化学定律都是以微分形式存在的。在生物学和经济学中，用微分方程来建立复杂的数学模型，微分方程在各个领域中得到广泛的应用。随着计算理论的发展，许多复杂的科学和工程问题得以分析解决，并用来估计未知的参数和状态。在化学反应动力学方面，对于化学工作者来说，如何应用和利用这些求解器，懂得求解器的优劣，并利用这些求解器求得更加准确的解是非常重要的。本文将介绍我们常用的求解常微分组刚性求解器如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）在提高求解器的准确度方面，最好要向求解器提供雅可比矩阵，使计算误差降至最小；求解器求解的准确度，受所设置的相对误差和绝对误差的影响很大，并不是相对误差和绝对误差设置的精度越高越好，要设置好相对误差和绝对误差，通常需要多次运算才能取得。