一种对刚性常微分方程的有效算法.doc

一种对刚性常微分方程的有效算法
Marek M. Stabrowski
华沙科技大学，IETiME, Koszykowa 7.5, PLO0661 华沙, 波兰d
收于 1995年7月10日; 1996年4月16日修订

摘要
Brayton-Gustavson-Hatchel法（BGH法）是向后微分公式法(BDF法)的一种，用于求解刚性常微分方程。本文介绍了BGH法的一些基础的细节，且描述了作者原创的一些对该算法的新改动。本文着重关注了减小计算量和改进误差控制两个方面。本文所举的两个例子，包括一个刚性和适度刚性的方程，体现出了BGH法相对于经典的吉尔法巨大的优越性。文中介绍的软件是可以升级的，它在PC/DOS平台和两个UNIX系统下运行良好，没有出现任何问题。1997 Elsevier Science B.V.

关键词：刚性微分方程，刚性微分系统

一、介绍
吉尔法是现在用于求解刚性常微分方程(ODES)的一种经典工具。它是一种相当稳定有效的办法(就求解速度而言)。它也是BDF法的一种。当然，其广泛使用不仅归功于在数学上的优势，也得益于一套完整fortran程序的发表。Byrne-Hindmarsh的吉尔法程序这样就成为了公共软件。
但是，即使在同类的BDF方法里面，也有其他一些相当有竞争力的方法。似乎其中一些算法的潜在优势并没有被深入的研究和探索过，其中BGH法就理应得到更多的关注。其作者声称，BGH法比吉尔法更稳定，尤其是在步长需要频繁变换的时候。并且，BGH法的计算量更小，误差控制也更灵活。在数学上,BGH法和吉尔法(甚至吉尔—诺德萨克法)地位并不相当。但是，吉尔—诺德萨克法按这样的方式进行一些改动之后，就和BGH法旗鼓相当了。
所有这些表明，BGH法是吉尔—诺德萨克法的有力竞争者。
可惜，目前发表的Fortran和C语言下的BGH法程序的特点是非常的刻板和特殊化。他们有些“特别”快速算法的特点，并没有真正库的风格。甚至，其多项式次数k被固定为两次，误差控制也被简化了，以至于对操作量的减少并没有完全的体现出来。

本文将介绍作者自己改进的BGHstiff程序的一些细节。作者本着库风格(就大小而言，完全可以升级)对BGH法进行了一些改进，增加了一些新特点。包括积分进程的有效流动，更加彻底和准确的误差控制，最后还有对减小计算量的完全探索。
第一个版本的BGHstiff程序是在流行的PC/DOS平台下开发的。对于UNIX平台的兼容也正在考虑当中。

二、Brayton-Gustavson-Hatchel法（BGH法）的基本原理
似乎对BGH方法进行一个简要的介绍是明智的。我们设带求解的常微分方程如下：

其中f是一个向量(方程组).t的取值空间[0,T]被分布不均的点t0,t1,t2……tn-k,tn-k+1…,tN分为长度为h= Δtn=tn+1-tn的小区间.在(1)式中的导数可以用下面的k次向后微分公式(修正)进行取代:

这样,常微分方程系统(1)就被转化为了下一个时间点tn+1的非线性方程系统.开始时间t0,初始条件x0下k=1,(2)式可解,于是时间t1下的向量x1也可以求出来.在求解向量x2的时候,多项式的次数k可以取到2.在6步之后,微分方程的k可以在1到6之间变化,其值取决于精确性和计算速度的标准.
方程式(1)的预告向量xPn+1也可以用BDF法计算,公式如下:

其中,预告向量xPn+1最好用做在修正期间迭代求解下式的初始值:

预告公式3)同样适用于求解局部舍位误差.可以证明,局部舍位误差在第k步的方法下可用下式表示:

这里是导数的真实值.同时,xn-1=x(tn-1)+0(hk+1).预告向量xPn+1由(3)式定义,而(4)式可以求出xn+1.
综上所述,BGH法可以归纳为以下五步:
(BGH-1) 建立新的时间点tn+1=tn+h
(BGH-2) 由预告公式(3)求出xPn+1的值.
(BGH-3 ) 用xPn+1作初值,再加上由(2)式求得的,解方程组(4),得到xn+1的值.
(BGH-4) 根据(5)式求得局部舍位误差,如果它不在规定的范围之内,那么需要改变积分步长或者所定k的值,回到步骤BGH-1)
(BGH-5) 把时间n增加到n+1,进入下一步的(BGH-1)

方程式(3)完全由k次多项式的预告系数γi来决定.这种情况下,需要用到如下的k+1次的代数方程:

其中

并且我们规定(tj-tj)0=1
同样的,校正系数αi也由时间点t0,t1,t2……,tn-k,tn-k+1…,tN来决定,用于求解(2)式.它也是由一个k+1词的代数方程来求解的:

其中
在Zimmerman and Baker的程序里面,(6)式和(11)式的求解都是用最简单的线性方程逼近法稍加改动而成的.

三、实际算法对计算量的减少
在Brayton, Gustavson and Hatchel的原文中，他们提到了一种更为有效的对预告参数γi和校正参数αi的计算方法。在这种方法里面，直接计算方程(6)和(11)完全是多余的.这种更为明智的方法缘于他们观察到γi(n+1,k)和αi(n+1,k)的值完全取决于γi(n,k)和αi(n,k)的值.其中n就是时间点tn中的n,而k就是多项式的次数k.在其他地方可以找出下式的具体来历.计算预告参数γi(n+1,k)的最终公式如下:

其中

且
同样的,计算校正参数αi(n+1,k)的公式如下:

初始计算由下式开始:

类似地,

上述(14)到(19)式就是由γi(n,k)和αi(n,k)的值来计算γi(n+1,k)和αi(n+1,k)的值的方法.这样我们就定义了从tn到tn+1时间点逐级递推求解α和γ的公式。
在本文介绍的BGH法程序里面，方程式（14）到(19)被完全的利用了起来。程序BGHstiff定义了一个表示时间点的ti的，有k+1个取值的1维向量t，和有k+2种不同时间点下取值(根据（14）到（17式）)的二维向量Aij。为了使方法的次数k变更的灵活性达到最大，这两个向量的存储空间都按照最高次数k=6来设置。这样，多项式的次数就能在1步之内从最少的1次上升到最多的6次。对这些向量内容的大多数计算都是在基础上简单推导出来的。这是显而易见的。
BGHstiff程序里也有类似直接地对于xi的计算。Brayton等人提到，xi，甚至Δxi的存储方式是一个周期性的缓冲区，里面有不断更新的点。这是一种非常优美，但是其效率却令人存疑的方式。首先，当多项式的次数k改变（上升或者下降）的时候，程序的复杂程度大大增加。其次，它引入了间接指针（双重指针），从而产生了BGHstiff程序所用直接指针所没有的巨大计算量开支。
其他一些BGHstiff程序提高效率的度都属于线性方程求解甚至软件开发里面一些几乎标准的改进方式。首先，二维向量被转化为了一维向量。一维向量的指针得到了“手动”的优化，也就是几乎完全的增减。其次，在校正计算中需要用到的一个Jacobian矩阵，每步积分只需要计算一次。Byrne和Hindmarsh有点过于极端了。校正系数α是由(16)和(17)式计算的，也就是说，取决于预告系数γ的值。最后，矩阵Aij的非对称性得到了充分的利用。我们只计算和更新其上三角部分（包括主对角线）的值。

四、程序BGHstiff中的补充和积分流动
在经典库风格的常微分方程解法里面，积分过程使用的是内部的时间步长，而输出的结果多种多样，完全是由用户指定的步长决定的，我们称为信息步长。一般来说信息步长是个常数，而内部步长是个变量，其范围有的时候还比较大。内部时间步长的变化源于误差控制标准的应用。一些通用的算法试图恰到好处的找到信息步长的大小。这个其实非常简单，可以分为两个步骤实现。第一个步骤，使用大于现有信息步长的时间点来进行积分。第二个步骤，通过对时间上限的简单缩短，减少最后一步积分的步长，就可以恰到好处地找到信息步长的大小。然后再用新的，已经减少的内部积分步长重新积分计算。在第二个步骤里面，没有必要对误差值进行控制。通过步长的减小，误差值自然而然就降低了。接下来的通信步长就正好按照通信时间点来进行了。