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

一种对刚性常微分方程的有效算法 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)式的求解都是用最简单的线性方程逼近法稍加改动而成的