多约束二阶非线性常微分方程拐点的数值求法.docVIP

多约束二阶非线性 常微分方程拐点的数值求法 杜志明 冯长根 (爆炸灾害预防、控制国家重点实验室，北京理工大学 100081) 摘要：本文提出了一种求解多约束二阶非线性常微分方程拐点的数值解法，并以具有平行反应的化学放热系统为例，给出了一个具体算例。 关键词：非线性; 常微分方程; 拐点; 数值计算; 热爆炸 二阶非线性常微分方程的求解是较为复杂的问题之一，而其极值点和拐点的确定又比求解方程的普通解要复杂得多，特别是一些强非线性的微分方程，即使数值求解也十分困难。在对科学研究和工程问题的描述过程中，包含有许多约束参数的非线性微分方程常常是不可避免的，尤其是方程的一些特殊点（如：极值点、拐点等）上的解往往具有十分重要的物理意义。拐点往往对应于物理上的转向点、转变点、临界分歧点等[1]，它是实际系统性质发生根本变化的关键点，因此，确定多约束微分方程的拐点不仅具有科学意义，也具有实际意义。 拐点条件 拐点是数学上的一类特殊点。对于带有约束参数的微分方程而言，微分方程的解与约束参数密切相关。每个约束参数都有自己的取值范围，当约束参数超出其取值范围时，常常导致方程无解。而当约束参数取一些特殊值（如在其取值范围边界上取值）时，经常对应于微分方程的特殊解，如极值解、拐点处的解等。显然，方程的解是约束参数的函数。本文中微分方程的拐点，指的就是微分方程的一组特殊解，这组解刚好是解随约束参数变化曲线上拐点处的值。下面我们寻求拐点条件。 带有多个约束参数的二阶非线性常微分方程具有以下形式： 边界条件为： 式中，是自变量，。是约束参数。是常数系数。 为了将边值问题转化为初值问题，设 , 若, 则有： 这样，边值问题就变成初值问题了。从积分到，若选得合适，式(3)将得到满足，即有 方程(1)~(3)的极值条件为[2]： 拐点在物理上相应于某一约束参数相应于某一状态变量的二阶导数等于零，即 ，该条件相当于其中一个约束参数相应于另一约束参数取得极值，即 =0。几何图形上拐点代表非线性曲线上极大值与极小值合二而一的点，也就是极值消失的点，根据定态方程组取得极值的条件[3,4]，式(4)和式(5)构成的方程组取得拐点处解的条件为[3,5]： (6) 由于 而 式(6)可简化为 上式就是方程 (1)~(3)的拐点条件。式(4) 、(5)与式(7)联立，即可求得原方程(1)~(3)的拐点。 数值算法 数值计算时，首先须进行数值积分，然后用迭代法求解式(4)、(5)和(7)构成的方程组的解。数值积分可采用Runge-Kutta-Merson法，迭代法选用Newton-Raphson法，其迭代法式为 其中 为了计算 ，需定义9个辅助函数： 各辅助函数及其初值条件可由方程(1)及其初值条件分别对和求导得出。 算例 本算例考虑的实际问题是内部具有放热化学反应的系统，用微分方程描述则是带化学反应项的热传导方程。由于化学反应速率与温度呈强非线性性质, 微分方程具有非线性性质。若考虑系统内同时存在多个平行的化学放热反应，且反应遵循Arrhenius定律，当系统达到热平衡（化学反应放热速率与系统向环境散热速率相等）时可用以下热平衡方程来描述[5] 边界条件可采用Frank-Kamenetskii条件[6] 式中，表示无量纲温度，为无量纲坐标，j是与形状有关的参数，是系统的约束参数，各参数的定义及物理意义可参阅文献[1,4]。 是系统的一个特殊参数，称为表观活化能。当参数超过某一数值时，系统温度将随参数连续变化，系统温度不再有突变性质, 由于系统内部热积累而导致的燃烧爆炸现象将不复存在。该数值记为, 在热爆炸理论中称为转变值，该点即是我们所要求解的拐点[1,5]。则表示系统内第i个化学反应的参数。 利用上面介绍的拐点条件和数值方法，可求方程(9)~(11)的拐点，现将数值积分时所用到的微分方程组罗列如下： 以上各式中 上述各微分方程联立起来数值积分后，再由式(8)进行迭代求解。此时 。 若假定系统内的放热反应只有两个，即n=2，求得的结果见表1。 表1 方程（9）~（11）的拐点数值（取） 0.1 0.253797 1.300368 4.833267 7.715381×10-2 0.2 0.262223 1.294365 4.775442 7.624096?10-2 0.3 0.271015 1.289437 4.725372 7.545112?10-2 0.